Comlex Valued Parameter Estimation in Julia using DifferentialEquations

Question

I want to estimate my complex parameters of an ODE using Optim and DifferentialEquations. I built an example case by changing the parameters in the example from the docs: https://diffeq.sciml.ai/latest/analysis/parameter_estimation/#Optimization-Based-ODE-Inference-Examples-1

I come to the point where I can built my loss objective, which then returns complex values.

using DifferentialEquations
using Optim
using DiffEqParamEstim
using RecursiveArrayTools # for VectorOfArray

# ODE function
function f(du,u,p,t)
  du[1] = p[1]*u[1] - u[1]*u[2]
  du[2] = dy = -3*u[2] + u[1]*u[2]
end

# A set of comlex Initial conditions and parameters
u0 = [1.0 + 1im ;1.0 + 1im]
tspan = (0.0,10.0)
p = [1.5 + 1.5im]
prob = ODEProblem(f,u0,tspan,p)

sol = solve(prob,Tsit5())
t = collect(range(0,stop=10,length=200))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false)
cost_function(1.0im+0.1)

result = optimize(cost_function, 0.0+ 0im, 10.0+ 10im)

Then I tried to use optimize() to optimize my cost function which returns the following Error:

LoadError: MethodError: no method matching optimize(::DiffEqObjective{DiffEqParamEstim.var"#43#48"{Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Integer,Tuple{Symbol,Symbol},Name
dTuple{(:maxiters, :verbose),Tuple{Int64,Bool}}},ODEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},true,Array{Complex{Float64},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},Array{Complex{Float64},2},Nothing,Nothing,Nothing},Nothing},DiffEqParamEstim.var"#47#53"{DiffEqParamEstim.var"#43#48"{Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Integer,Tuple{Symbol,Symbol},NamedTuple{(:maxiters, :verbose),Tuple{Int64,Bool}}},ODEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},true,Array{Complex{Float64},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},Array{Complex{Float64},2},Nothing,Nothing,Nothing},Nothing}}}, ::Complex{Float64}, ::Complex{Float64})
Closest candidates are:
  optimize(::Any, ::Number, ::Number, ::AbstractArray{T,N} where N; kwargs...) where T at /home/.julia/packages/Optim/L5T76/src/multivariate/solvers/constrained/fminbox.jl:156
  optimize(::Any, ::Number, ::Number, ::AbstractArray{T,N} where N, ::Optim.AbstractConstrainedOptimizer) where T at /home/.julia/packages/Optim/L5T76/src/multivariate/solvers/constrained/fminbox.jl:160
  optimize(::Any, ::Number, ::Number, ::AbstractArray{T,N} where N, ::Optim.AbstractConstrainedOptimizer, ::Optim.Options; kwargs...) where T at /home/.julia/packages/Optim/L5T76/src/multivariate/solvers/constrained/fminbox.jl:160
  ...
Stacktrace:
 [1] top-level scope at /home/Documents/random stuff/minimal_example.jl:27
in expression starting at /home/Documents/random stuff/minimal_example.jl:27

Is it possible to use this method on complex valued functions? And if yes what would be the correct way of doing so?

EDIT:

using DifferentialEquations
using Optim
using DiffEqParamEstim
using RecursiveArrayTools # for VectorOfArray

# ODE function
function f(du::AbstractArray{T}, u::AbstractArray{T},p::AbstractArray{T},t) where T<:Complex
  du[1] = p[1]*u[1] - u[1]*u[2]
  du[2] = dy = -3*u[2] + u[1]*u[2]
end

# A set of complex Initial conditions and parameters
u0 = [1.0 + 1im ;1.0 + 1im]
tspan = (0.0,10.0)
p = [1.5 + 1.5im]
prob = ODEProblem(f,u0,tspan,p)

sol = solve(prob,Tsit5())
sol.u
t = collect(range(0,stop=10,length=200))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data[1,:]),
                                     maxiters=10000,verbose=false)

result = optimize(cost_function, u0,BFGS())


InexactError: Float64(-567409.1782943244 - 312256.811660436im)
in top-level scope at minimal_example.jl:27
in optimize at Optim/L5T76/src/multivariate/optimize/interface.jl:115
in optimize at Optim/L5T76/src/multivariate/optimize/interface.jl:115
in #optimize#93 at Optim/L5T76/src/multivariate/optimize/interface.jl:116 
in optimize at Optim/L5T76/src/multivariate/optimize/optimize.jl:33 
in initial_state at Optim/L5T76/src/multivariate/solvers/first_order/bfgs.jl:66
in value_gradient!! at NLSolversBase/mGaJg/src/interface.jl:82
in setproperty! at base/Base.jl:21
in convert at base/number.jl:7
in Real at base/complex.jl:37

Brent cannot be used on Complex numbers because there is no ordering. Try BFGS on an array of 1 complex number — Chris Rackauckas, May 13 '20 at 10:14
What do you mean by array of 1 complex number? I tried using: ` cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data[1,:]), maxiters=10000,verbose=false) ` result = optimize(cost_function, BFGS()) But it still throws an error — , May 13 '20 at 17:17
I updated the question and wrote it there because it was too long to put in a comment. — , May 14 '20 at 06:45
Thanks again! I added the initial conditions but now I get a InexactError in the NLSolverBase. I thought it would be enougth to add type assertions to my function f?( I edited the question again) — , May 14 '20 at 08:00
Looks like Optim just doesn't support Complex numbers. Make it BFGS optimization where your initial condition is two values, and inside of your function turn those two real values into a complex value. — Chris Rackauckas, May 14 '20 at 08:29

Comlex Valued Parameter Estimation in Julia using DifferentialEquations

0 Answers0