Numeric value substitution into SymPy Symbolic array or matrix

Specific Domains Numerics
1

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
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
  [7] getindex
    @ ./broadcast.jl:610 [inlined]
  [8] copy
    @ ./broadcast.jl:912 [inlined]
  [9] materialize
    @ ./broadcast.jl:873 [inlined]
 [10] (::SymPy.var"#122#124"{Vector{SymPy.var"#118#119"}})(args::Float64)
    @ SymPy ~/.julia/packages/SymPy/gkr7m/src/lambdify.jl:253
 [11] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
 [12] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
 [13] #31
    @ ./broadcast.jl:1088 [inlined]
 [14] macro expansion
    @ ./ntuple.jl:72 [inlined]
 [15] ntuple
    @ ./ntuple.jl:69 [inlined]
 [16] copy
    @ ./broadcast.jl:1088 [inlined]
 [17] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, SymPy.var"#122#124"{Vector{SymPy.var"#118#119"}}, Tuple{NTuple{5, Float64}}})
    @ Base.Broadcast ./broadcast.jl:873
 [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!

2

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
3

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.

4

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.

5

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
6

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

1 Like