Search code examples
booleansympyderivative

sympy derivative with boolean

I am trying to take the derivative of a function including a boolean variable with sympy.

My expected result:

Two different derivatives, depending on the boolean being either True or False (i.e. 1 or 0).

Example:

import sympy as sy
c, x = sy.symbols("c x", positive=True, real=True)
bo = sy.Function("bo")
fct1 = sy.Function("fct1")
fct2 = sy.Function("fct2")
FOC2 = sy.Function("FOC2")
y = 5
a = 2
b = 4


def fct1(x):
    return -0.004*x**2 + 0.25*x + 4
# the following gives the smaller positive intercept with the x-axis)
# this intercept is the threshold value for the boolean function, bo
min(sy.solve(fct1(x)-y, x))


def bo(x):
    if fct1(x) <= y:
        return 1
    else:
        return 0


def fct2(c, x):
    return a + b*c + bo(x)*c


def FOC2(c, x):
    return sy.diff(fct2(c, x), c)
print(FOC2(c, x))

The min-function after the comments shows me the threshold of x for bo being True or False would be 4.29..., thus positive and real.

Output:

TypeError: cannot determine truth value of Relation

I understand that the truth value depends on x, which is a symbol. Thus, without knowing x one cannot determine bo.

But how would I get my expected result, where bo is symbolic?

Solution

  • First off, I would advise you to carefully consider what is going on in your code the way it is pasted above. You first define a few sympy functions, e.g.

    fct1 = sy.Function("fct1")

    So after this, fct1 is an undefined sympy.Function - undefined in the sense that it is neither specified what its arguments are, nor what the function looks like.

    However, then you define same-named functions explicitly, as in

    def fct1(x):
    return -0.004*x**2 + 0.25*x + 4

    Note however, that at this point, fct1 ceases to be a sympy.Function, or any sympy object for that matter: you overwrite the old definition, and it is now just a regular python function!

    This is also the reason that you get the error: when you call bo(x), python tries to evaluate

    -0.004*x**2 + 0.25*x + 4 <= 5

    and return a value according to your definition of bo(). But python does not know whether the above is true (or how to make that comparison), so it complains.

    I would suggest 2 changes:

    1. Instead of python functions, as in the code, you could simply use sympy expressions, e.g.

      fct1 = -0.004*x**2 + 0.25*x + 4

    2. To get the truth value of your condition, I would suggest to use the Heaviside function (wiki), which evaluates to 0 for a negative argument, and to 1 for positive. Its implementation in sympy is sympy.Heaviside. Your code could then look as follows:

    import sympy as sy
c, x = sy.symbols("c x", positive=True, real=True)
y = 5
a = 2
b = 4


fct1 = -0.004*x**2 + 0.25*x + 4
bo = sy.Heaviside(y - fct1)
fct2 = a + b*c + bo * c

FOC2 = sy.diff(fct2, c)

print(FOC2)

    Two comments on the line

    bo = sy.Heaviside(y - fct1)

    (1) The current implementation does not evaluate sympy.Heaviside(0)by default; this is beacause there's differing definitions around (some define it to be 1, others 1/2). You'd want it to be 1, to be in accordance with the (weak) inequality in the OP. In sympy 1.1, this can be achieved by passing an additional argument to Heaviside, namely whatever you want Heaviside(0) to evaluate to:

    bo = sy.Heaviside(y - fct1, 1)

    This is not supported in older versions of sympy.

    (2) You will get your FOC2, again involving a Heaviside term. What I like about this, is that you could keep working with this expression, say if you wanted to take a second derivative and so on. If, for the sake of readability, you would prefer a piecewise expression - no problem. Just replace the according line with 

    bo = sy.Heaviside(y - fct1)._eval_rewrite_as_Piecewise(y-fct1)

    Which will translate to a piecewise function automatically. (note that under older versions, this automatically implicitly uses Heaviside(0) = 0.5 - best to use (1) and (2) together:

    bo = sy.Heaviside(y - fct1, 1)._eval_rewrite_as_Piecewise(y-fct1)

    Unfortunately, I don't have a working sympy 1.1 at my hands right now and can only test the old code.

    One more noteconcerning sympy's piecewise functions: they are much more readable if using sympy's latex printing, by inserting

    sy.init_printing()

    early in the code.

    (Disclaimer: I am by no means an expert in sympy, and there might be other, preferable solutions out there. Just trying to make a suggestion!)