Protected Division Using Sympy

Question

A protected division is a normal division but when you divide by 0 it returns a fixed constant (usually 1).

def protected_div(x, y):
    if y == 0:
        return 1
    return x/y

Is there a way to use this as an operator on sympy (For example replacing the standard division)?

Here is an example of what I want:

>>> import sympy as sym
>>> x = sym.Symbol('x')
>>> expr = 1/x #(protected division goes here?)
>>> expr.subs(x, 0)
1

The division has to be protected at evaluation time.

EDIT 1:

What I've tried:

1. Using sym.lambidify with the modules parameter set:

>>> x = sym.Symbol('x')
>>> expr = 1/x
>>> lamb = sym.lambdify(x, expr, modules={'/':protected_div})
>>> print(lamb(0))
ZeroDivisionError: 0.0 cannot be raised to a negative power

This does not work because sympy converts 1/x to x**(-1) when lambidifying. I tried overriding the power operator but I don't know the function name. I've tried 'Pow', 'pow', '**' and none worked.

However if i declare the expression as expr = 1.0/x it actually does not convert to a negative power, however it does not use my custom division function. I think these types of functions are not overridable using the module parameter.

2. @Zaz suggestion:

class floatsafe(float):
    def __truediv__(self, __x):
        if __x == 0:
            return floatsafe(1)
        return super().__truediv__(__x)

x = sym.Symbol('x')
expr = floatsafe(1)/x
print(expr.subs(x, floatsafe(0)))

Returns

zoo

Which is complex infinity.

I tried combining this approach with sym.lambdify, but the dividend is converted to a float after I lambdify the function.

In the case that the dividend is variable it also does not work:

x = sym.Symbol('x')
expr = x/0.0
a = sym.lambdify(x, expr, modules={'/':floatsafe.__truediv__})
print(inspect.getsource(a))
print(a(floatsafe(0)))

Outputs

def _lambdifygenerated(x):
    return nan*x

nan

EDIT: There seems to some confusion around why I'd want that. It's for a genetic programming algorithm using sympy. A protected division is a common operator in GP so that the created solutions are valid.

You've already shown one way to do this so it's unclear what you want. Maybe you want to use a Piecewise? — Oscar Benjamin, May 14 '22 at 22:56
If i use this function with say sympy.Symbols('x y') it just returns x/y. Then If I try to evaluate for y=0 it does not work. I need it so the division is protected in a sympy expression, so that x/0 = 1. — user, May 15 '22 at 01:55

Zaz · Answer 1 · 2022-05-14T23:33:38.303

The regular mathematics we use on the day-to-day is a ring on the set of real numbers, ℝ: The properties of a ring are that you have two operations (such as multiplication and addition) and one of them (such as addition) will always produce another number within the set.

You can create a more specific notion of a field (such that both operations will always produce another member in the set) by removing 0 or expanding the set to the hyperreals.

My point being, without knowing what problem exactly you're trying to solve, I would guess that instead of redefining division, it makes more sense to redefine the number system that you're using: For whatever reason, you have some system of numbers that should return 1 when divided by zero, so why not create a subclass of float, for example?

class floatD01(float):
    def __truediv__(self, divisor):
        if divisor == 0:
            return 1
        return self/divisor

You may also want to scan help(float) for any other methods related to division that you may want to change such as __divmod__, __floordiv__ (7//3 == 2), etc, and have a hard think about how you want this new mathematical group that you're creating to work and why.

Other options that may potentially be more robust would be to go nuclear and try catching all ZeroDivisionErrors and replace them with one (either by modifying the class) or within whatever code you're running or, if appropriate, implementing something like what the language R extensively uses: NA values. I'm sure there's some way (I believe in numpy) to do something along the lines of: C = [1/3, 2/2, 3/1, 4/0] # == [1/3, 2/2, 3/1, NA] sum(C) = 4.333

score 0 · Accepted Answer · answered May 20 '22 at 22:12

The solution was pretty simple actually, although I was not able to actualy overload the division operator all I had to do was create a sympy function for the protected division and use that instead.

class protected_division(sym.Function):
    @classmethod
    def eval(cls, x, y):
        if y.is_Number:
            if y.is_zero:
                return sym.S.One
            else:
                return x/y

Then just use that in an expression:

>>> expr = protected_division(1, sym.Symbol('x'))
protected_division(1, x)
>>> expr.subs(sym.Symbol('x'), 0)
1
>>> expr.subs(sym.Symbol('x'), 3)
1/3

I did not find out how to make the class tell sym.lambdify what to do in case of a "lambdification", but you can use the modules parameters for that:

>>> def pd(x, y):
...     if y == 0:
...         return 1
...     return x/y
... 
>>> l = sym.lambdify(sym.Symbol('x'), expr, modules={'protected_division': pd})
>>> l(3)
1.6666666666666667
>>> l(0)
1

Protected Division Using Sympy

2 Answers2