Symbolic interval matrix power

sympower_imatrix.jl

# =====================================
# Version using MacroTools + SymEngine
# =====================================

using IntervalMatrices
using SymEngine
using MacroTools: postwalk

subidx(i) = join(Char.(0x2080 .+ convert.(UInt16, digits(i)[end:-1:1])))

function symbolic_matrix(n; name="M")
    M = Matrix{SymEngine.Basic}(undef, n, n)
    for i in 1:n, j in 1:n
        M[i, j] = name * subidx(i) * subidx(j)
    end
    return M
end

function subs(ex, imat; name="M")
    n = size(imat, 1) # imat assumed square
    for i in 1:n, j in 1:n
        name_ij = Symbol(name * subidx(i) * subidx(j))
        ex == name_ij && return imat[i, j]
    end
    return ex
end

Msym_eval(Msym, M, k) = [postwalk(ex -> subs(ex, M), convert(Expr, m)) for m in Msym]

# this version gives a more accurate result than just M^k because lots of intermediate
# interval evaluations are saved. however, it is very bad in terms of performance because
# the expression tree is consdered for each term..
function power(M, k)
    @assert size(M, 1) == size(M, 2)
    Msym = symbolic_matrix(size(M, 1))
    X = Msym_eval(Msym^k, M, k)
    return eval.(X)
end

# =====================================
# Version using DynamicPolynomials
# =====================================
using DynamicPolynomials

# this version is as fast as M^k where M is an inteval, but it is not accurate. the reason is that 
# DynamicPolynomials expands the symbolic polynomial products 
function power_poly(M::IntervalMatrix, k::Int)
    n = size(M, 1)
    @polyvar Msym[1:n, 1:n]
    vars = vec(Msym) # symbolic variables, column-major
    Msym_k = Msym^k # symbolic matrix multiplication << can we "hold" this product?
    Mk = map(x -> subs(x, vars => vec(M)), Msym_k) # substitution of the intervals
    return Mk
end


using SymEngine: toString

# this works but it is slow
# it has been reported here: 
# https://discourse.julialang.org/t/creating-julia-function-from-a-symbolic-expression-using-symengine/21476
# using lambify2 didn't help, it gives world age issues (see below)
function power2(M, k)
    n = size(M, 1)
    @assert n == size(M, 2)
    Msym = symbolic_matrix(n)
    Msymᵏ = Msym^k
    vars = Tuple(Msym)
    Msym_callable = map(x -> lambdify(x, vars), Msymᵏ)
    vals = Tuple(M)
    Mᵏ = map(Mij-> Mij(vals...), Msym_callable)
    return Mᵏ
end

#lambdify2(ex, vars=free_symbols(ex)) = eval(Expr(:function, Expr(:call, gensym(), map(Symbol,vars)...),
#                                       convert(Expr, ex)))
#=
MethodError: no method matching ##371(::Interval{Float64}, ::Interval{Float64}, ::Interval{Float64}, ::Interval{Float64})
The applicable method may be too new: running in world age 26092, while current world is 26096.
Closest candidates are:
  ##371(::Any, ::Any, ::Any, ::Any) at none:0 (method too new to be called from this world context.)
=#
lambdify3(ex, vars=free_symbols(ex)) = eval(Meta.parse("ex(x)="*SymEngine.toString(ex...)))

function power3(M, k)
    n = size(M, 1)
    @assert n == size(M, 2)
    Msym = symbolic_matrix(n)
    vars = Tuple(Msym)
    Msymᵏ = Msym^k
    
    # create callable entries
    Msymᵏ_callable = []
    for j in 1:n, i in 1:n
        push!(Msymᵏ_callable, eval(Meta.parse("M$i$j("* vars_str_concat * ") = " * SymEngine.toString(Msymᵏ[i, j]))))
    end

    # evaluate them
    vals = Tuple(M)
    Mᵏ = map(Mij-> Mij(vals...), Msymᵏ_callable)
    Mᵏ = reshape(Mᵏ, n, n)
    return Mᵏ
end

# ===========================
# Version using TaylorSeries
# ===========================

function power4(M, k)
    n = size(M, 1)
    @assert n == size(M, 2)

    # the order should be sufficiently high (>= k ?)
    
    vars = set_variables("M", numvars=n^2, order=k+1) # coeffs are float
    #vars = set_variables(Interval{Float64}, "M", numvars=n^2, order=k+1) # coeffs are intervals
    
    Msym = reshape(vars, n, n)
    Msymᵏ = Msym^k

    vals = vec(M)
    Mᵏ = map(Mij -> evaluate(Mij, vals), Msymᵏ)
    
    return Mᵏ
end