Avoiding "complex division by zero" error in Numpy (related to cmath)

Question

I would like to avoid ZeroDivisionError: complex division by zero in my calculation, getting nan at that exception.

Let me lay out my question with a simple example:

from numpy import *
from cmath import *

def f(x) :
    x = float_(x) 
    return 1./(1.-x)

def g(x) :
    x = float_(x)
    return sqrt(-1.)/(1.-x)

f(1.)   # This gives 'inf'.
g(1.)   # This gives 'ZeroDivisionError: complex division by zero'.

It is my intention to get g(1.) = nan, or at least anything but an error that interrupts the calculation. First question: how can I do it?

Importantly, I would not like to modify the code inside the functions (for example, inserting conditions for exceptions, as is done in this answer), but rather keep it in its current form (or even deleting the x = float_(x) line if possible, as I mention below). The reason is that I am working with a long code with dozens of functions: I would like them all to avoid ZeroDivisionError without the need of making lots of changes.

I was forced to insert x = float_(x) to avoid a ZeroDivisionError in f(1.). Second question: would there be a way of suppressing this line but still get f(1.) = inf without modifying at all the code defining f?

---

EDIT:

I have realized that using cmath (from cmath import *) is responsible for the error. Without it, I get g(1.) = nan, which is what I want. However, I need it in my code. So now the first question turns into the following: how can I avoid "complex division by zero" when using cmath?

---

EDIT 2:

After reading the answers, I have made some changes and I simplify the question, getting closer to the point:

import numpy as np                            
import cmath as cm                            

def g(x) :                                    
    x = np.float_(x)                         
    return cm.sqrt(x+1.)/(x-1.)    # I want 'g' to be defined in R-{1}, 
                                   # so I have to use 'cm.sqrt'. 

print 'g(1.) =', g(1.)             # This gives 'ZeroDivisionError: 
                                   # complex division by zero'.

Question: how can I avoid the ZeroDivisionError not modifying the code of my function g?

Answer 1

I would expect that the function handled the errors by itself, but if not, then you could do it when you call the function (I do beleive this is more work than modifying the fnctions though).

try:
    f(1.)   # This gives 'inf'.
except ZeroDifivisionError:
    None    #Or whatever
try:
    g(1.)   # This gives 'ZeroDivisionError: complex division by zero'.
except ZeroDifivisionError:
    None    #Or whatever

Or, as the the answer you referred to tells, call your function with a numpy number:

f(np.float_(1.))
g(np.float_(1.))

=======
On my machine, this exact code gives me a warning and no error. I don't think it's possible to get what you want without catching errors...

def f(x):
    return 1./(1.-x)

def g(x):
    return np.sqrt(-1.)/(1.-x)

print f(np.float_(1.))
>> __main__:2: RuntimeWarning: divide by zero encountered in double_scalars
>>inf

print g(np.float_(1.))
>> __main__:5: RuntimeWarning: invalid value encountered in sqrt
>>nan

Answer 2

I still don't understand why you need to use cmath. Type-cast x into np.complex_ when you expect complex output, and then use np.sqrt.

import numpy as np

def f(x):
    x = np.float_(x)
    return 1. / (1. - x)

def g(x):
    x = np.complex_(x)
    return np.sqrt(x + 1.) / (x - 1.)

This yields:

>>> f(1.)
/usr/local/bin/ipython3:3: RuntimeWarning: divide by zero encountered in double_scalars
  # -*- coding: utf-8 -*-
Out[131]: inf

>>> g(-3.)
Out[132]: -0.35355339059327379j

>>> g(1.)
/usr/local/bin/ipython3:3: RuntimeWarning: divide by zero encountered in cdouble_scalars
  # -*- coding: utf-8 -*-
/usr/local/bin/ipython3:3: RuntimeWarning: invalid value encountered in cdouble_scalars
  # -*- coding: utf-8 -*-
Out[133]: (inf+nan*j)

The drawback being, of course, that the function g will always end up giving you complex output, which could later cause problems if you feed its result back into f, for instance, because that now type-casts to float, and so on... Maybe you should just type-cast to complex everywhere. But this will depend on what you need to achieve on some bigger scale.

EDIT

It turns out that there's a way to get g to return complex only if it needs to be. Use numpy.emath.

import numpy as np

def f(x):
    x = np.float_(x)
    return 1. / (1. - x)

def g(x):
    x = np.float_(x)
    return np.emath.sqrt(x + 1.) / (x - 1.)

This now gives what you would expect, converting to complex only when necessary.

>>> f(1)
/usr/local/bin/ipython3:8: RuntimeWarning: divide by zero encountered in double_scalars
Out[1]: inf

>>> g(1)
/usr/local/bin/ipython3:12: RuntimeWarning: divide by zero encountered in double_scalars
Out[2]: inf

>>> g(-3)
Out[3]: -0.35355339059327379j

Answer 3

The error arises because your program is using the sqrt function from the cmath library. This is a good example of why you should avoid (or at least be careful with) importing the entire library using from library import *.

For instance, if you reverse the import statements then you will have no problem:

from cmath import *
from numpy import *

def f(x):
    x = float_(x) 
    return 1./(1.-x)


def g(x) :
    x = float_(x)
    return sqrt(-1.)/(1.-x)

Now g(1.) returns nan, because the sqrt function of numpy is used (which can handle negative numbers). But that's still bad practice: if you have a large file, then it is not clear which sqrt is being used.

I recommend always using named imports such as import numpy as np and import cmath as cm. Then the function sqrt is not imported, and to define g(x) you need to write np.sqrt (or cm.sqrt).

However, you noted that you don't want to change your function. In that case you should only import those functions from each library that you need, i.e.

from cmath import sin  # plus whatever functions you are using in that file
from numpy import sqrt, float_

def f(x):
    x = float_(x) 
    return 1./(1.-x)


def g(x) :
    x = float_(x)
    return sqrt(-1.)/(1.-x)

Unfortunately, you cannot easily get rid of the conversion to numpy float_, because python floats are not the same as numpy floats so this conversion needs to take place at some point.