How to avoid math domain error when we have sqrt

Question

I wrote the following code, after printing about 10 outputs it gets error

return 1/sqrt(Om*(1+z)**3+omg0*(1+z)**6+(1-Om-omg0))
ValueError: math domain error

the three first functions are correct and work well in other programs please do not be sensitive about their structures, I think the problem is in for-loop and choosing some values which cannot satisfy the condition of sqrt. Could I add some lines to tell Python avoid numbers which lead to math domain error? If yes, How should I do that? I mean passing the step which leads to negative value under sqrt?

Cov is a 31*31 mtrix, xx[n] is a list of 31 numbers.

def ant(z,Om,w):
    return 1/sqrt(Om*(1+z)**3+w*(1+z)**6+(1-Om-w))

def dl(n,Om,w,M):  
    q=quad(ant,0,xx[n],args=(Om,w))[0]
    h=5*log10((1+xx[n])*q)
    fn=(yy[n]-M-h)                  
    return fn

def c(Om,w,M):
    f_list = []
    for i in range(31):  # the value '2' reflects matrix size
        f_list.append(dl(i,Om,w,M))
    A=[f_list]
    B=[[f] for f in f_list]
    C=np.dot(A,Cov)
    D=np.dot(C,B)
    F=np.linalg.det(D)*0.000001
    return F

N=100
for i in range (1,N):
    R3=np.random.uniform(0,1)

    Omn[i]=Omo[i-1]+0.05*np.random.normal()
    wn[i]=wo[i-1]+0.05*np.random.normal()
    Mn[i]=Mo[i-1]+0.1*np.random.normal()

    L=exp(-0.5*(c(Omn[i],wn[i],Mn[i])-c(Omo[i-1],wo[i-1],Mo[i-1])))

    if L>R3:
        wo[i]=wn[i]

    else:
        wo[i]=wo[i-1]

    print(wo[i])

The output is:

0.12059556415714912
0.16292726528216397
0.16644447885609648
0.1067588804671105
0.0321446951572128
0.0321446951572128
0.013169965429457382
Traceback (most recent call last):
......
return 1/sqrt(Om*(1+z)**3+omg0*(1+z)**6+(1-Om-omg0))
ValueError: math domain error

`math.sqrt` doesn't accept negative numbers, which is probably the problem. What do you want `ant` to do when the input is negative? Return a complex number? Return 0? — Alex Hall, Apr 18 '18 at 20:20
If complex values are ok with you and your specific problem use module `cmath` instead of `math` — Paul Panzer, Apr 18 '18 at 20:21
@AlexHall we have random number, some of them lead to negative value under sqrt, How should I tell python to pass them. — Bob, Apr 18 '18 at 20:21
What do you mean by 'pass them'? For example you could do nothing when `f_list.append(dl(i,Om,w,M))` fails (because `dl` fails because `ant` fails) and thus end up with `f_list` having fewer elements. Would that work? — Alex Hall, Apr 18 '18 at 20:26
@AlexHall assume, random number `1`, `2`, `3` lead a negative value, I want to jump of them and calculate with another or next value — Bob, Apr 18 '18 at 20:27
So is it fine if `f_list` has fewer than 31 elements? Or do you want to keep increasing `i` until the size is right? — Alex Hall, Apr 18 '18 at 20:36
@Jean-FrançoisFabre I think in this case this one can be accepted, Thank you, but for guarantee I need some thing more general. But thank you. What are the differences between `sqrt` and `**0.5` — Bob, Apr 18 '18 at 20:36
@AlexHall no it contains a veeeery big matrix `31*31`, `760*760` and so on — Bob, Apr 18 '18 at 20:37
@Bob ask google, which redirects to StackOverflow :) https://stackoverflow.com/questions/33684948/difference-between-1-2-math-sqrt-and-cmath-sqrt — Jean-François Fabre, Apr 18 '18 at 20:41

Alex Hall · Accepted Answer · 2018-04-18T21:06:48.317

2

Here are some options:

Replace sqrt(...) with (...)**0.5. This will produce complex numbers, which may or may not be acceptable. For example (-1)**0.5 produces i which in Python will appear as 1j (ignoring floating point errors).
Catch errors and move on. Since you probably want to catch the errors higher up, I recommend translating the ValueError from sqrt into a custom error:

class SqrtError(ValueError):
    pass

def ant(z, Om, w):
    val = Om * (1 + z) ** 3 + w * (1 + z) ** 6 + (1 - Om - w)
    try:
        sq = sqrt(val)
    except ValueError:
        raise SqrtError
    return 1 / sq

Then it sounds like you want to keep trying new random numbers until you get one that works, which could be done like so:

for i in range(1, N):
    R3 = np.random.uniform(0, 1)
    while True:

        Omn[i] = Omo[i - 1] + 0.05 * np.random.normal()
        wn[i] = wo[i - 1] + 0.05 * np.random.normal()
        Mn[i] = Mo[i - 1] + 0.1 * np.random.normal()

        try:
            L = exp(-0.5 * (c(Omn[i], wn[i], Mn[i]) - c(Omo[i - 1], wo[i - 1], Mo[i - 1])))
        except SqrtError:
            continue
        else:
            break

    if L > R3:
        wo[i] = wn[i]

    else:
        wo[i] = wo[i - 1]

edited Apr 18 '18 at 21:06

answered Apr 18 '18 at 20:52

Alex Hall

34,833
5
57
89

3 questions: `1` you mean using `**0.5` lead to wrong answer? `2` the first box jump of the error? `3` the second box, I am not sure that does what I want. my code are saying: calculate `L` if it is bigger than a random number accept new values but if it is smaller than that random number choose the value of previous step. – Bob Apr 18 '18 at 21:01
the first box works well, but speed of calculating goes down. – Bob Apr 18 '18 at 21:04
@Bob `1` we don't know what the right or wrong answer is. It's not clear what you want to do, we're all guessing. I don't know if putting complex numbers in your code will do the wrong thing. `2` the first box doesn't change much, it just means that in the second box I can write `except SqrtError` instead of the more general `except ValueError` which might catch things I don't want it to (like errors from `log10`). `3` I've edited the second box to clarify and only generate `R3` once. – Alex Hall Apr 18 '18 at 21:11
If you don't understand the meaning of `try/except/raise/else/while True/continue/break` please go read carefully about them. – Alex Hall Apr 18 '18 at 21:12
the first box after printing some output gets this error:`File "C:\.......py", line 68, in ant raise SqrtError SqrtError` – Bob Apr 18 '18 at 21:13
The first box is useless without the second box. It doesn't sound like you understand at all the constructs that I'm using. Please go read about them, that is essential to solving your problem and an important part of learning to code. – Alex Hall Apr 18 '18 at 21:14
No man, I got it entirely. But after too many errors my mind is blowing up. I will read your answer in couple of mins later. sorry for this questions. – Bob Apr 18 '18 at 21:16
It sounds correct, the only problem which I run into it, is, my loop is `N=100` but the out put are 37. maybe 63 numbers lead the sqrt to negative value! – Bob Apr 18 '18 at 21:24

touch my body · Answer 2 · 2018-04-18T20:29:51.037

0

Your code is trying to find the square root of a negative number.

This is not allowed with real-valued numbers (for obvious mathematical reasons), so your best bet is to use a try/except block:

def ant(z,Om,w):
    try:
        return 1/sqrt(Om*(1+z)**3+omg0*(1+z)**6+(1-Om-omg0))
    except ValueError:
        return None

edited Apr 18 '18 at 20:29

answered Apr 18 '18 at 20:23

touch my body

1,634
22
36

should I write this block you wrote under `def name():`? – Bob Apr 18 '18 at 20:28
Hard to say... I don't see where this function is even being used. If you posted where `ant` is being called I could advise – touch my body Apr 18 '18 at 20:29
`ant` is used in the second function `dl(n,Om,w,M): ` in the line of `q=quad(ant,0,xx[n],args=(Om,w))[0]` – Bob Apr 18 '18 at 20:30
This code is difficult to read for formatting reasons. What I can see is you're using `np.random.normal()` to generate gaussian noise. This is the source of your negative numbers, and this is why you're getting negative square roots. – touch my body Apr 18 '18 at 20:36
yes my friend, exactly, because of this, I need something to jump of numbers which lead to negative value. If your recommended code is correct how should I use it under a `def` block? – Bob Apr 18 '18 at 20:39
Hard to say. One way is to just take the absolute value of everything, but I don't know if that's appropriate in your case. Not sure what your code is supposed to do. If that's fine, just add `z = abs(z)` (and the others) under the `def ant` you have now. Again, not sure if that's appropriate since I don't know what your code is supposed to do – touch my body Apr 18 '18 at 20:44
Ok thank you, but my code in the third function is calculating the multiplication of 3 matrixes, the items of these 3 matrixes are calculating in second function, and second function which is a integration calls another function which is the first function. I need a code line which I don't know how, to tell Python, ok, you are checking `if` and` loop` but please check and jump of those numbers which leads the sqrt in the first def a negative value. in the comments a user told me to use `**0.5` – Bob Apr 18 '18 at 20:52

How to avoid math domain error when we have sqrt

2 Answers2