# Numeric value substitution into SymPy Symbolic array or matrix

I am having trouble substituting numerical values into symbolic arrays and matrices.
I have read about `lambdify`, but I still fail to provide numerical values for the arguments in a consistent way.

As an example, consider the following `Julia` code

``````julia> using Sympy
julia> n = 4;
julia> η = [symbols("η\$i", positive=true) for i in 1:n];
julia> α = [(i<3) ? 1 : symbols("α", positive=true) for i in 1:n];
julia> η .*= α
4-element Vector{Sym}:
η₁
η₂
α⋅η₃
α⋅η₄

julia> syms = Symbol.(free_symbols(tensors[:η]))
5-element Vector{Symbol}:
:α
:η1
:η2
:η3
:η4

julia> f = lambdify(η, syms)
#122 (generic function with 1 method)

julia> f.(0.5,1.0,1.0,1.0,1.0)
4-element Vector{Float64}:
1.0
1.0
0.5
0.5
``````

As expected, the outcome is indeed the substitution for `α=0.5` and all `η=1.0`. However, my problem arises from the fact that I would like to feed the expressions from arrays or tuples, i.e. what I would like is:

``````julia> values = (0.5, 1.0, 1.0, 1.0, 1.0);
julia> f.(values)
ERROR: MethodError: no method matching var"##339"(::Float64)

Closest candidates are:
var"##339"(::Any, ::Any, ::Any, ::Any, ::Any)
@ SymPy none:0

Stacktrace:
[1] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Base ./essentials.jl:816
[2] invokelatest(::Any, ::Any, ::Vararg{Any})
@ Base ./essentials.jl:813
[3] (::SymPy.var"#118#119"{SymPy.var"###339"})(::Float64, ::Vararg{Float64})
@ SymPy ~/.julia/packages/SymPy/gkr7m/src/lambdify.jl:223
[4] map(::Function, ::Float64)
@ Base ./number.jl:284
[7] getindex
[8] copy
[9] materialize
[10] (::SymPy.var"#122#124"{Vector{SymPy.var"#118#119"}})(args::Float64)
@ SymPy ~/.julia/packages/SymPy/gkr7m/src/lambdify.jl:253
[13] #31
[14] macro expansion
@ ./ntuple.jl:72 [inlined]
[15] ntuple
@ ./ntuple.jl:69 [inlined]
[16] copy
[18] top-level scope
@ REPL[199]:1
``````

The reason that I want this is so that I can make all of my symbolic arrays into “lambdified” expressions and evaluate them, for example at multiple points in parallel. Additionally, when the array is instead a matrix, I also do not wish to manually input all the matrix elements – especially when the sizes of the matrices grow.
Consider for example,

``````julia> M = [symbols("M\$i\$j") for i in 1:3 for j in 1:3].reshape(3,3)
3×3 Matrix{Sym}:
M₁₁  M₁₂  M₁₃
M₂₁  M₂₂  M₂₃
M₃₁  M₃₂  M₃₃
julia> syms_M = free_symbols.(M)
3×3 Matrix{Vector{Sym}}:
[M11]  [M12]  [M13]
[M21]  [M22]  [M23]
[M31]  [M32]  [M33]
julia> fM = lambdify(M, syms_M)
#122 (generic function with 1 method)
``````

How do I use `fM` to return a numeric valued matrix, e.g. when `Mvals = rand(3,3)`? Any help would be greatly appreciated!

You would write

``````julia> values = (0.5, 1.0, 1.0, 1.0, 1.0);
julia> f.(values)
``````

if `f` is a function of one parameter and you want to call it 5 times and gather the 5 results in a vector.

In your case however `f` is a function that takes 5 parameters. So to call it with a tuple of values you can write

``````julia> f(values...)
4-element Vector{Float64}:
1.0
1.0
0.5
0.5
``````

For the matrix example: here you broadcast the `free_symbols` call. I don’t think that will work: it returns a matrix where each element is the result of a `free_symbols` call, i.e. each element is a vector of symbols. When you call `lambdify` with this matrix of vectors of symbols, for some reason it returns something instead of raising an error, but the returned function doesn’t seem usable. It’s not surprising since you should call `lambdify` with a vector of symbols instead.

To solve this you could convert the matrix of vectors into a simple vector. But the easiest is to call `free_symbols` directly with a vector:

``````julia> M = [symbols("M\$i\$j") for i in 1:3, j in 1:3]
3×3 Matrix{Sym}:
M₁₁  M₁₂  M₁₃
M₂₁  M₂₂  M₂₃
M₃₁  M₃₂  M₃₃

julia> syms_M = free_symbols(vec(M))
9-element Vector{Sym}:
M₁₁
M₁₂
M₁₃
M₂₁
M₂₂
M₂₃
M₃₁
M₃₂
M₃₃

julia> fM = lambdify(M, syms_M)
#122 (generic function with 1 method)

julia> Mvals = rand(3,3)
3×3 Matrix{Float64}:
0.0253515  0.261792  0.640684
0.668682   0.216239  0.433714
0.780887   0.808188  0.315145

julia> fM(Mvals...)
3×3 Matrix{Float64}:
0.0253515  0.668682  0.780887
0.261792   0.216239  0.808188
0.640684   0.433714  0.315145
``````

Note: I think the splatting of a matrix with `...` is not very efficient. If that’s a problem, maybe you can solve it by using a matrix type from StaticArrays.jl (so that the matrix size is known by the compiler).

1 Like

Thanks for the quick response. I was just starting to type up this, as I just figured it out. It makes indeed sense to use splatting. The efficiency of the splatting is not really a deal breaker, but I might need to convert to a `StaticArray` if it appears to be a problem.

Another question I have is that I do not know how to assign values to elements that are a function of multiple variables. I can do it in this way, most likely, but then I need to know the order of symbols that `free_symbols(expr)` returns – which is odd as I do not necessarily know this beforehand.

Consider for example a simple Lotka-Volterra system (which are similar to the systems I need this substitution for):

``````julia> M = reshape([-symbols("a\$i\$j", positive=true) for i in 1:2 for j in 1:2], (2,2))
2×2 Matrix{Sym}:
-a₁₁  -a₂₁
-a₁₂  -a₂₂

julia> r = diagm([symbols("η", positive=true), -symbols("μ", positive=true)])
2×2 Matrix{Sym}:
η   0
0  -μ

julia> A = r+M
2×2 Matrix{Sym}:
-a₁₁ + η      -a₂₁
-a₁₂  -a₂₂ - μ

julia> syms = free_symbols(vec(A))
4-element Vector{Sym}:
a₁₂
a₂₁
η
μ

julia> Avals = vec([0 1.0; 1.0 0]);
julia> ηval = 1.0;
julia> μval = 0.5;
julia> values = vcat(Avals, ηval, μval);
julia> f = lambdify(A, syms)
julia> f(values...)
2×2 Matrix{Float64}:
1.0  -1.0
-1.0  -0.5
``````

as required. However, this forces me to know the order of variables that `free_symbols(expr)` returns, as if I instead have `values = vcat(ηval, μval, Avals)`, then `f(values...)` will understandably be wrong. Is there a way to have `lambdify()` return a function with named arguments, or something like that? Or is there another way that I can safely insert numerical values of each symbol subsequently, so that I know for sure that substituting the value of a symbol results in the evaluation of only that symbol within the expression? I hope my question is clear.

So diving into the source code of `SymPy` we can find the `free_symbols()` function definition. In there, one can see that the sorting is executed by calling `sortperm`:

``````fs[sortperm(string.(fs))]   # some sorting order to rely on
``````

This gives us a way to reliably generate the order, without having to call it beforehand. Using this we are free to specify a dictionary that holds all the values corresponding to the symbols present in the (matrix) expression, order the accordingly, and call the `lambdify`’d function with the appropriate splatted argument list.
For example:

``````julia> using SymPy, LinearAlgebra

julia> A = reshape([-symbols("a\$i\$j", positive=true) for i in 1:2 for j in 1:2], (2,2))
2×2 Matrix{Sym}:
-a₁₁  -a₂₁
-a₁₂  -a₂₂

julia> r = diagm([symbols("η", positive=true), -symbols("μ", positive=true)])
2×2 Matrix{Sym}:
η   0
0  -μ

julia> Ω = r + A
2×2 Matrix{Sym}:
-a₁₁ + η      -a₂₁
-a₁₂  -a₂₂ - μ

julia> Avals = [0.0 1.0; 1.0 0.0]
2×2 Matrix{Float64}:
0.0  1.0
1.0  0.0

julia> adict = Dict([Pair(symbols("a\$i\$j"), Avals[i,j]) for i in 1:2 for j in 1:2])
Dict{Sym, Float64} with 4 entries:
a22 => 0.0
a11 => 0.0
a12 => 1.0
a21 => 1.0

julia> ηdict = Dict(Pair(symbols("η"), 1.0))
Dict{Sym, Float64} with 1 entry:
η => 1.0

julia> μdict = Dict(Pair(symbols("μ"), 0.5))
Dict{Sym, Float64} with 1 entry:
μ => 0.5

julia> var2val = merge(adict, ηdict, μdict)
Dict{Sym, Float64} with 6 entries:
a22 => 0.0
η   => 1.0
μ   => 0.5
a11 => 0.0
a12 => 1.0
a21 => 1.0

julia> order = sortperm(string.(keys(var2val)))
6-element Vector{Int64}:
4
5
6
1
2
3

julia> vals = collect(values(var2val))[order]
6-element Vector{Float64}:
0.0
1.0
1.0
0.0
1.0
0.5

julia> fΩ = lambdify(Ω, free_symbols(vec(Ω)))
#122 (generic function with 1 method)

julia> fΩ(vals...)
2×2 Matrix{Float64}:
1.0  -1.0
-1.0  -0.5
``````

This can probably be done more cleanly considering the available parameters, but the one that defines the problem will probably know the parameters and can therefore construct the `var2val` dictionary without any issue.

The `free_symbols` function is just a convenient way to get a list of symbols, but you don’t have to use it. If you know the symbols you can build the list yourself explicitly, and this way you control the order:

``````julia> A0 = [symbols("a\$i\$j", positive=true) for i in 1:2, j in 1:2]
2×2 Matrix{Sym}:
a₁₁  a₁₂
a₂₁  a₂₂

julia> @vars η μ positive=true
(η, μ)

julia> r = [η 0; 0 μ];

julia> A = r-A0
2×2 Matrix{Sym}:
-a₁₁ + η      -a₁₂
-a₂₁  -a₂₂ + μ

julia> Avals = vec([0 1.0; 1.0 0]);

julia> ηval = 1.0;

julia> μval = 0.5;

julia> syms = [vec(A0); η; μ] # here I choose an order for the symbols
6-element Vector{Sym}:
a₁₁
a₂₁
a₁₂
a₂₂
η
μ

julia> f = lambdify(A, syms);

julia> f(values...)
2×2 Matrix{Float64}:
1.0  -1.0
-1.0   0.5
``````
1 Like

That makes sense as well… and is probably a bit easier. Thanks!

1 Like