How to use sympy to locally switch to use of "set of equations" from standard sequence of lines of assignments?

Question

The following Python function has a flaw in the sequence of steps of calculating f(x).

def f(x):
   f = np.zeros((3))
   f[0] = x[0]
   f[1] = f[2] + x[1]
   f[2] = x[2]
   return f

The input x= [1,1,1] gives [1,1,1] but should give [1,2,1] provided we interpret the code mathematically as a system of equations. How can I use sympy to sort out the evaluation order and get the "right" answer?

Note, I am not at all looking for solving the equations system symbolically, but just want the evaluation order correct and avoid work with mapping out that myself.

Answer 1

You can let sympy calculate the formulas for f0, f1 and f2:

import sympy as sym
f0, f1, f2 = sym.symbols('f0,f1,f2')
x0, x1, x2 = sym.symbols('x0,x1,x2')

eq1 = sym.Eq(f0, x0)
eq2 = sym.Eq(f1, f2 + x1)
eq3 = sym.Eq(f2, x2)

result = sym.solve([eq1,eq2, eq3],(f0, f1, f2))
print(result)
>>> {f0: x0, f1: x1 + x2, f2: x2}

Answer 2

Using sympy, you could write the function as:

import numpy as np

def f(x):
    from sympy import symbols, Eq, solve
    f = symbols('f0:3')
    sol = solve([Eq(f[0], x[0]), Eq(f[1], f[2] + x[1]), Eq(f[2], x[2])], f)
    return np.array([sol[fi].evalf() for fi in f], dtype=float)

f(np.array([1, 1, 1]))

This supposes that you are only using a simple example, and that in reality the equations are much more complicated. Note that the approach can get more involved if your equations have more than one valid solution.

lambdify() can convert an expression to numpy format. For this to work here, the function f should get some symbolic input parameters and convert them to output expressions.

Here is how it could look like.

import numpy as np
from sympy import symbols, Eq, solve, lambdify

def f(x):
    y = symbols('y0:3')
    sol = solve([Eq(y[0], x[0]), Eq(y[1], y[2] + x[1]), Eq(y[2], x[2])], y)
    return [sol[yi].simplify() for yi in y]

x0, x1, x2 = symbols('x0:3')
f_np = lambdify((x0, x1, x2), f((x0, x1, x2)))

# calling the function on a 1D numpy array
f_np(*np.array([1, 1, 1]))

# calling the function on a 2D numpy array and convert the result back to a numpy array
a = np.array([[1, 1, 1], [2, 2, 2]])
np.array(f_np(a[:, 0], a[:, 1], a[:, 2])).T  # array([[1, 2, 1], [2, 4, 2]])

Here, f_np is a numpy function that gets three inputs (the * unpacks the array) and has a list of 3 as output. Directly working with arrays as input and output doesn't seem possible at the moment.

In an interactive environment, you can enter f_np? (or f_np??) to see that functions source:

>> f_np?
Signature: f_np(x0, x1, x2)
Docstring:
Created with lambdify. Signature:
func(x0, x1, x2)
Expression:
(x0, x1 + x2, x2)
Source code:
def _lambdifygenerated(x0, x1, x2):
    return ((x0, x1 + x2, x2))

f_np can then be used in scipy optimizers and even be compiled (Cython or numba). This approach assumes the equations stay the same.

Answer 3

You fundamentally misunderstand what = does in Python. It's not mathematical equality (which works both ways), it's assignment. When you say a = b you say "forget whatever a was, it's b from now on".

So going step-by-step we get:

   f[0] = x[0]         # f[0] becomes x[0] = 1
   f[1] = f[2] + x[1]  # f[1] becomes f[2] + x[1] = 0 + 1 = 1
   f[2] = x[2]         # f[2] becomes x[2] = 1

---

Note, I am not at all looking for solving the equations system symbolically, but just want the evaluation order correct and avoid work with mapping out that myself.

You're either going to need to solve the equations symbolically yourself, or do it with a computer algebra system. Regular Python code is imperative and does neither of these things (out of the box).

See e.g. kabooya's answer for setting up such a system of equations using Sympy.