CODEWITHSUNDEEP
pythonnumpysympy
Ian Roberts
Ian Roberts

Reputation: 1

How to substitute 2D IndexedBase variable in Sympy

Can I substitute a 2D IndexedBase that is in an expression for a numpy array using sympy.subs() or sympy.evalf(subs=())?

So something like:

i, j, m, n = sp.symbols('i j m n', integer=True)
x = sp.IndexedBase('x')
a = sp.IndexedBase('a')
b = sp.IndexedBase('b')

f = sp.ln(sp.Sum(sp.exp(sp.Sum(a[i, j]*x[j]+b[i], (j, 1, n))), (i, 1, m)))

which currently outputs: expression output

Note that A is a 2D IndexedBase and B is a 1D IndexedBase. I would like to be able to substitute the IndexedBases defined above for numpy arrays, something like

A = np.random.uniform(MIN, MAX, (M, N))
B = np.random.uniform(MIN, MAX, M)

f.evalf(subs={a:A, b:B, x:(1,5), m:M, n:N}

How would I be able to achieve something like this?

Upvotes: 0

Views: 109

Answers (1)

Oscar Benjamin
Oscar Benjamin

Reputation: 14480

You can do this with lambdify like this:

In [36]: i, j, m, n = sp.symbols('i j m n', integer=True)
    ...: x = sp.IndexedBase('x')
    ...: a = sp.IndexedBase('a')
    ...: b = sp.IndexedBase('b')
    ...: 
    ...: f = sp.ln(sp.Sum(sp.exp(sp.Sum(a[i, j]*x[j]+b[i], (j, 0, n-1))), (i, 0, m-1)))

In [37]: f_lam = lambdify((a, b, x, m, n), f, 'numpy')

In [38]: MIN, MAX, M, N = 0, 1, 3, 4

In [39]: A = np.random.uniform(MIN, MAX, (M, N))
    ...: B = np.random.uniform(MIN, MAX, M)

In [40]: X = np.random.uniform(MIN, MAX, N)

In [41]: f_lam(A, B, X, M, N)
Out[41]: 4.727558334863294

Note that I shifted the limits of the Sum go from e.g. 0 to n-1 to match up with numpy indexing rules.

Also note that this uses the ordinary Python sum and a double-loop generator expression rather than more efficient numpy operations:

In [42]: import inspect

In [44]: print(inspect.getsource(f_lam))
def _lambdifygenerated(Dummy_34, Dummy_33, Dummy_32, m, n):
    return log((builtins.sum(exp((builtins.sum(Dummy_32[j]*Dummy_34[i, j] + Dummy_33[i] for j in range(0, n - 1+1)))) for i in range(0, m - 1+1))))

The generated code is more efficient if you set this up using matrix symbols rather than explicit expressions involving IndexedBase:

In [22]: m, n = symbols('m, n')

In [23]: a = MatrixSymbol('a', m, n)

In [24]: b = MatrixSymbol('b', m, 1)

In [25]: x = MatrixSymbol('x', n, 1)

In [26]: one = OneMatrix(1, m)

In [28]: f = log(one*(a*x + b))

In [29]: f_lam = lambdify((a, b, x, m, n), f, 'numpy')

In [30]: f_lam(A, B, X, M, N)
Out[30]: array([1.52220638])

In [33]: print(inspect.getsource(f_lam))
def _lambdifygenerated(a, b, x, m, n):
    return log((ones((1, m))).dot((a).dot(x) + b))

Upvotes: 1

Related Questions