SciMLSensitivity.jl icon indicating copy to clipboard operation
SciMLSensitivity.jl copied to clipboard

ReverseDiffVJP and ForwardDiff failing on complex ODE

Open amilsted opened this issue 3 years ago • 11 comments

I want to autodiff through a complex ODE (the Schrödinger equation) with solution v(t) = exp(1.0im * t * H), where H is a matrix. The problem is inherently complex. I encountered issues, so I tried a complex-typed version of this example . This also fails.

using DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff, Zygote

function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0; 1.0]
prob = ODEProblem(fiip,complex(u0),(0.0,10.0),complex(p))

function sum_of_solution(x)
    u0 = complex(x[1:2])
    p = complex(x[3:end])  # also fails (differently) if these are left real
    _prob = remake(prob,u0=u0,p=p)
    real(sum(solve(_prob,Tsit5(),rtol=1e-6,atol=1e-6,saveat=0.1)))
end

dx = ForwardDiff.gradient(sum_of_solution, [u0; p])
# Cannot determine ordering of Dual tags Nothing and ForwardDiff.Tag{typeof(sum_of_solution), Float64}

dx = Zygote.gradient(sum_of_solution, [u0; p])
# TypeError: in TrackedReal, in V, expected V<:Real, got Type{ComplexF64}

Will post full stacktraces below.

amilsted avatar Nov 15 '21 18:11 amilsted

dx = ForwardDiff.gradient(sum_of_solution, [u0; p])

Cannot determine ordering of Dual tags Nothing and ForwardDiff.Tag{typeof(sum_of_solution), Float64}

Stacktrace:
  [1] ≺(a::Type, b::Type)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/aBsxl/src/dual.jl:54
  [2] muladd
    @ ~/.julia/packages/ForwardDiff/aBsxl/src/dual.jl:157 [inlined]
  [3] muladd(z::Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}, x::ForwardDiff.Dual{Nothing, Float64, 6}, w::Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}})
    @ Base ./complex.jl:325
  [4] muladd(x::ForwardDiff.Dual{Nothing, Float64, 6}, z::Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}, y::Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}})
    @ Base ./complex.jl:323
  [5] macro expansion
    @ ~/.julia/packages/OrdinaryDiffEq/Zi9Zh/src/initdt.jl:140 [inlined]
  [6] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [7] ode_determine_initdt(u0::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, t::Float64, tdir::Float64, dtmax::Float64, abstol::ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, reltol::ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), prob::ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Tuple{Float64, Float64}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Nothing, Float64, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, Float64, Float64, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}, ODESolution{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}, 2, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}}, ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Tuple{Float64, Float64}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}, Vector{Float64}, Vector{Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.DEOptions{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Zi9Zh/src/initdt.jl:139
  [8] auto_dt_reset!
    @ ~/.julia/packages/OrdinaryDiffEq/Zi9Zh/src/integrators/integrator_interface.jl:329 [inlined]
  [9] handle_dt!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Nothing, Float64, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, Float64, Float64, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}, ODESolution{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}, 2, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}}, ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Tuple{Float64, Float64}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}, Vector{Float64}, Vector{Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.DEOptions{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Zi9Zh/src/solve.jl:504
 [10] __init(prob::ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Tuple{Float64, Float64}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Zi9Zh/src/solve.jl:466
 [11] #__solve#471
    @ ~/.julia/packages/OrdinaryDiffEq/Zi9Zh/src/solve.jl:4 [inlined]
 [12] #solve_call#42
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:61 [inlined]
 [13] solve_up(prob::ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, Tuple{Float64, Float64}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, p::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:rtol, :atol, :saveat), Tuple{Float64, Float64, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:87
 [14] #solve#43
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:73 [inlined]
 [15] sum_of_solution(x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}})
    @ Main ./In[1]:14
 [16] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/aBsxl/src/apiutils.jl:37 [inlined]
 [17] vector_mode_gradient(f::typeof(sum_of_solution), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/aBsxl/src/gradient.jl:106
 [18] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/aBsxl/src/gradient.jl:19
 [19] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(sum_of_solution), Float64}, Float64, 6}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/aBsxl/src/gradient.jl:17
 [20] top-level scope
    @ In[2]:1
 [21] eval
    @ ./boot.jl:360 [inlined]
 [22] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

amilsted avatar Nov 15 '21 18:11 amilsted

dx = Zygote.gradient(sum_of_solution, [u0; p])

TypeError: in TrackedReal, in V, expected V<:Real, got Type{ComplexF64}

Stacktrace:
  [1] ReverseDiff.TrackedArray(value::Vector{ComplexF64}, deriv::Vector{ComplexF64}, tape::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/tracked.jl:86
  [2] track(x::Vector{ComplexF64}, ::Type{ComplexF64}, tp::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/tracked.jl:452
  [3] (::ReverseDiff.var"#661#662"{ComplexF64, Vector{ReverseDiff.AbstractInstruction}})(x::Vector{ComplexF64})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:46
  [4] map
    @ ./tuple.jl:215 [inlined]
  [5] ReverseDiff.GradientConfig(input::Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{Float64}}, ::Type{ComplexF64}, tp::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:46
  [6] ReverseDiff.GradientConfig(input::Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{Float64}}, tp::Vector{ReverseDiff.AbstractInstruction}) (repeats 2 times)
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:37
  [7] ReverseDiff.GradientTape(f::Function, input::Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{Float64}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/tape.jl:204
  [8] adjointdiffcache(g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, dg::Nothing, f::ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}; quad::Bool, noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/adjoint_common.jl:144
  [9] DiffEqSensitivity.ODEInterpolatingAdjointSensitivityFunction(g::Function, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, dg::Nothing, f::Function, checkpoints::Vector{Float64}, tols::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, tstops::Nothing; noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/interpolating_adjoint.jl:72
 [10] DiffEqSensitivity.ODEInterpolatingAdjointSensitivityFunction(g::Function, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, dg::Nothing, f::Function, checkpoints::Vector{Float64}, tols::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, tstops::Nothing)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/interpolating_adjoint.jl:23
 [11] ODEAdjointProblem(sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; checkpoints::Vector{Float64}, callback::Nothing, reltol::Float64, abstol::Float64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/interpolating_adjoint.jl:255
 [12] _adjoint_sensitivities(sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; abstol::Float64, reltol::Float64, checkpoints::Vector{Float64}, corfunc_analytical::Nothing, callback::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/sensitivity_interface.jl:17
 [13] adjoint_sensitivities(::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, ::Vararg{Any, N} where N; sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, kwargs::Base.Iterators.Pairs{Symbol, Union{Nothing, Float64}, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :rtol, :atol), Tuple{Nothing, Float64, Float64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/sensitivity_interface.jl:6
 [14] (::DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/concrete_solve.jl:272
 [15] ZBack
    @ ~/.julia/packages/Zygote/Y1xbK/src/compiler/chainrules.jl:168 [inlined]
 [16] (::Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}})(dy::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/chainrules.jl:194
 [17] #203
    @ ~/.julia/packages/Zygote/Y1xbK/src/lib/lib.jl:203 [inlined]
 [18] (::Zygote.var"#1734#back#205"{Zygote.var"#203#204"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{true}, Bool}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}}}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [19] Pullback
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:73 [inlined]
 [20] (::typeof(∂(#solve#43)))(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface2.jl:0
 [21] (::Zygote.var"#203#204"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/lib/lib.jl:203
 [22] (::Zygote.var"#1734#back#205"{Zygote.var"#203#204"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [23] Pullback
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:68 [inlined]
 [24] (::typeof(∂(solve##kw)))(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface2.jl:0
 [25] Pullback
    @ ./In[1]:14 [inlined]
 [26] (::typeof(∂(sum_of_solution)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface2.jl:0
 [27] (::Zygote.var"#50#51"{typeof(∂(sum_of_solution))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface.jl:41
 [28] gradient(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface.jl:76
 [29] top-level scope
    @ In[3]:1
 [30] eval
    @ ./boot.jl:360 [inlined]
 [31] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

amilsted avatar Nov 15 '21 18:11 amilsted

What happens with EnzymeVJP? We should make it avoid ReverseDiffVJP when it's complex since it seems like that's the reason for the Zygote issue, but if EnzymeVJP doesn't work then there isn't a valid VJP choice 😅

ChrisRackauckas avatar Nov 15 '21 20:11 ChrisRackauckas

Well, just autojacvec=true I guess.

ChrisRackauckas avatar Nov 15 '21 21:11 ChrisRackauckas

Enzyme crashed Julia. I will try again with Julia 1.7-rc2 at some point...

amilsted avatar Nov 15 '21 21:11 amilsted

Try with autojacvec=true?

ChrisRackauckas avatar Nov 15 '21 22:11 ChrisRackauckas

Tried autojacvec=true, but nothing appears to change. Zygote stacktrace looks the same. I also tried setting sensealg=QuadratureAdjoint(), but there was no obvious change.

TypeError: in TrackedReal, in V, expected V<:Real, got Type{ComplexF64}

Stacktrace:
  [1] ReverseDiff.TrackedArray(value::Vector{ComplexF64}, deriv::Vector{ComplexF64}, tape::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/tracked.jl:86
  [2] track(x::Vector{ComplexF64}, ::Type{ComplexF64}, tp::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/tracked.jl:452
  [3] (::ReverseDiff.var"#661#662"{ComplexF64, Vector{ReverseDiff.AbstractInstruction}})(x::Vector{ComplexF64})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:46
  [4] map
    @ ./tuple.jl:215 [inlined]
  [5] ReverseDiff.GradientConfig(input::Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}}, ::Type{ComplexF64}, tp::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:46
  [6] ReverseDiff.GradientConfig(input::Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}}, tp::Vector{ReverseDiff.AbstractInstruction}) (repeats 2 times)
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:37
  [7] ReverseDiff.GradientTape(f::Function, input::Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/tape.jl:204
  [8] adjointdiffcache(g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, sensealg::QuadratureAdjoint{0, true, Val{:central}, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, dg::Nothing, f::ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}; quad::Bool, noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/adjoint_common.jl:144
  [9] DiffEqSensitivity.ODEQuadratureAdjointSensitivityFunction(g::Function, sensealg::QuadratureAdjoint{0, true, Val{:central}, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, dg::Nothing)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/quadrature_adjoint.jl:12
 [10] ODEAdjointProblem(sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, sensealg::QuadratureAdjoint{0, true, Val{:central}, Bool}, g::Function, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing, callback::Nothing)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/quadrature_adjoint.jl:68
 [11] _adjoint_sensitivities(sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, sensealg::QuadratureAdjoint{0, true, Val{:central}, Bool}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; abstol::Float64, reltol::Float64, callback::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/quadrature_adjoint.jl:252
 [12] adjoint_sensitivities(::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ::Vararg{Any, N} where N; sensealg::QuadratureAdjoint{0, true, Val{:central}, Bool}, kwargs::Base.Iterators.Pairs{Symbol, Union{Nothing, Real}, NTuple{4, Symbol}, NamedTuple{(:callback, :rtol, :atol, :autojacvec), Tuple{Nothing, Float64, Float64, Bool}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/sensitivity_interface.jl:6
 [13] (::DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, QuadratureAdjoint{0, true, Val{:central}, Bool}, Vector{ComplexF64}, Vector{Float64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/concrete_solve.jl:272
 [14] ZBack
    @ ~/.julia/packages/Zygote/AlLTp/src/compiler/chainrules.jl:204 [inlined]
 [15] (::Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, QuadratureAdjoint{0, true, Val{:central}, Bool}, Vector{ComplexF64}, Vector{Float64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}}})(dy::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/chainrules.jl:230
 [16] #208
    @ ~/.julia/packages/Zygote/AlLTp/src/lib/lib.jl:203 [inlined]
 [17] (::Zygote.var"#1734#back#210"{Zygote.var"#208#209"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, QuadratureAdjoint{0, true, Val{:central}, Bool}, Vector{ComplexF64}, Vector{Float64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol, :autojacvec), Tuple{Float64, Float64, Bool}}}}}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [18] Pullback
    @ ~/.julia/packages/DiffEqBase/b1nST/src/solve.jl:73 [inlined]
 [19] (::typeof(∂(#solve#43)))(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface2.jl:0
 [20] (::Zygote.var"#208#209"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/lib/lib.jl:203
 [21] (::Zygote.var"#1734#back#210"{Zygote.var"#208#209"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [22] Pullback
    @ ~/.julia/packages/DiffEqBase/b1nST/src/solve.jl:68 [inlined]
 [23] (::typeof(∂(solve##kw)))(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface2.jl:0
 [24] Pullback
    @ ./In[18]:12 [inlined]
 [25] (::typeof(∂(sum_of_solution)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface2.jl:0
 [26] (::Zygote.var"#55#56"{typeof(∂(sum_of_solution))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface.jl:41
 [27] gradient(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface.jl:76
 [28] top-level scope
    @ In[19]:1
 [29] eval
    @ ./boot.jl:360 [inlined]
 [30] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

amilsted avatar Nov 17 '21 17:11 amilsted

That stack trace is impossible with autojacvec=true. How did you do it?

ChrisRackauckas avatar Nov 20 '21 13:11 ChrisRackauckas

I was using autojacvec=true as an argument to solve(). From a search of the repo, that was probably incorrect? Here is what I am doing now:

function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0; 1.0]
prob = ODEProblem(fiip,complex(u0),(0.0,10.0),complex(p))

function sum_of_solution(u0, p)
    _prob = remake(prob,u0=u0,p=p)
    real(sum(solve(_prob,Tsit5(),rtol=1e-6,atol=1e-6,saveat=0.1, sensealg=InterpolatingAdjoint(autojacvec=true))))
end

dx = Zygote.gradient(sum_of_solution_2, complex(u0), complex(p))

The response is:

TypeError: in TrackedReal, in V, expected V<:Real, got Type{ComplexF64}

Stacktrace:
  [1] ReverseDiff.TrackedArray(value::Vector{ComplexF64}, deriv::Vector{ComplexF64}, tape::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/tracked.jl:86
  [2] track(x::Vector{ComplexF64}, ::Type{ComplexF64}, tp::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/tracked.jl:452
  [3] (::ReverseDiff.var"#661#662"{ComplexF64, Vector{ReverseDiff.AbstractInstruction}})(x::Vector{ComplexF64})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:46
  [4] map
    @ ./tuple.jl:215 [inlined]
  [5] ReverseDiff.GradientConfig(input::Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{Float64}}, ::Type{ComplexF64}, tp::Vector{ReverseDiff.AbstractInstruction})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:46
  [6] ReverseDiff.GradientConfig(input::Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{Float64}}, tp::Vector{ReverseDiff.AbstractInstruction}) (repeats 2 times)
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/Config.jl:37
  [7] ReverseDiff.GradientTape(f::Function, input::Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{Float64}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/OciDO/src/api/tape.jl:204
  [8] adjointdiffcache(g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, dg::Nothing, f::ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}; quad::Bool, noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/adjoint_common.jl:144
  [9] DiffEqSensitivity.ODEInterpolatingAdjointSensitivityFunction(g::Function, sensealg::InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, dg::Nothing, f::Function, checkpoints::Vector{Float64}, tols::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, tstops::Nothing; noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/interpolating_adjoint.jl:72
 [10] DiffEqSensitivity.ODEInterpolatingAdjointSensitivityFunction(g::Function, sensealg::InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, discrete::Bool, sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, dg::Nothing, f::Function, checkpoints::Vector{Float64}, tols::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, tstops::Nothing)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/interpolating_adjoint.jl:23
 [11] ODEAdjointProblem(sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; checkpoints::Vector{Float64}, callback::Nothing, reltol::Float64, abstol::Float64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/interpolating_adjoint.jl:255
 [12] _adjoint_sensitivities(sol::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, g::DiffEqSensitivity.var"#df#235"{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; abstol::Float64, reltol::Float64, checkpoints::Vector{Float64}, corfunc_analytical::Nothing, callback::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/sensitivity_interface.jl:17
 [13] adjoint_sensitivities(::ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ::Vararg{Any, N} where N; sensealg::InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, kwargs::Base.Iterators.Pairs{Symbol, Union{Nothing, Float64}, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :rtol, :atol), Tuple{Nothing, Float64, Float64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/sensitivity_interface.jl:6
 [14] (::DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/uakCr/src/concrete_solve.jl:272
 [15] ZBack
    @ ~/.julia/packages/Zygote/AlLTp/src/compiler/chainrules.jl:204 [inlined]
 [16] (::Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}})(dy::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/chainrules.jl:230
 [17] #208
    @ ~/.julia/packages/Zygote/AlLTp/src/lib/lib.jl:203 [inlined]
 [18] (::Zygote.var"#1734#back#210"{Zygote.var"#208#209"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#234"{Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, InterpolatingAdjoint{0, true, Val{:central}, Bool, Bool}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{}, Colon, NamedTuple{(:rtol, :atol), Tuple{Float64, Float64}}}}}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [19] Pullback
    @ ~/.julia/packages/DiffEqBase/b1nST/src/solve.jl:73 [inlined]
 [20] (::typeof(∂(#solve#43)))(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface2.jl:0
 [21] (::Zygote.var"#208#209"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/lib/lib.jl:203
 [22] (::Zygote.var"#1734#back#210"{Zygote.var"#208#209"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))}})(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [23] Pullback
    @ ~/.julia/packages/DiffEqBase/b1nST/src/solve.jl:68 [inlined]
 [24] (::typeof(∂(solve##kw)))(Δ::FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface2.jl:0
 [25] Pullback
    @ ./In[7]:17 [inlined]
 [26] (::typeof(∂(sum_of_solution_2)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface2.jl:0
 [27] (::Zygote.var"#55#56"{typeof(∂(sum_of_solution_2))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface.jl:41
 [28] gradient(::Function, ::Vector{ComplexF64}, ::Vararg{Vector{ComplexF64}, N} where N)
    @ Zygote ~/.julia/packages/Zygote/AlLTp/src/compiler/interface.jl:76
 [29] top-level scope
    @ In[11]:1
 [30] eval
    @ ./boot.jl:360 [inlined]
 [31] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

amilsted avatar Nov 30 '21 23:11 amilsted

Enzyme works:

using OrdinaryDiffEq, Zygote, DiffEqSensitivity
function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0; 1.0]
prob = ODEProblem(fiip,complex(u0),(0.0,10.0),complex(p))

function sum_of_solution(u0, p)
    _prob = remake(prob,u0=u0,p=p)
    real(sum(solve(_prob,Tsit5(),reltol=1e-6,abstol=1e-6,saveat=0.1, sensealg=InterpolatingAdjoint(autojacvec=EnzymeVJP()))))
end

dx = Zygote.gradient(sum_of_solution, complex(u0), complex(p))
(ComplexF64[-39.12610372553748 - 0.0im, -8.787926840582909 - 0.0im], ComplexF64[8.307609744273716 + 0.0im, -159.48459534795904 + 0.0im, 75.20354257564327 + 0.0im, -339.19349921463834 + 0.0im])

Oh oops I meant autojacvec=false

using OrdinaryDiffEq, Zygote, DiffEqSensitivity
function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0; 1.0]
prob = ODEProblem(fiip,complex(u0),(0.0,10.0),complex(p))

function sum_of_solution(u0, p)
    _prob = remake(prob,u0=u0,p=p)
    real(sum(solve(_prob,Tsit5(),reltol=1e-6,abstol=1e-6,saveat=0.1, sensealg=InterpolatingAdjoint(autojacvec=false))))
end

dx = Zygote.gradient(sum_of_solution, complex(u0), complex(p))
ArgumentError: Cannot create a dual over scalar type ComplexF64. If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{ComplexF64}) = true.
throw_cannot_dual(V::Type) at dual.jl:41
ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2}(value::ComplexF64, partials::ForwardDiff.Partials{2, ComplexF64}) at dual.jl:18
_broadcast_getindex_evalf at broadcast.jl:648 [inlined]
_broadcast_getindex at broadcast.jl:621 [inlined]
getindex at broadcast.jl:575 [inlined]
macro expansion at broadcast.jl:984 [inlined]
macro expansion at simdloop.jl:77 [inlined]
copyto! at broadcast.jl:983 [inlined]
copyto! at broadcast.jl:936 [inlined]
materialize! at broadcast.jl:894 [inlined]
materialize! at broadcast.jl:891 [inlined]
seed!(duals::Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2}}, x::Vector{ComplexF64}, seeds::Tuple{ForwardDiff.Partials{2, ComplexF64}, ForwardDiff.Partials{2, ComplexF64}}) at apiutils.jl:65
vector_mode_dual_eval!(f!::SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2}}}}, y::Vector{ComplexF64}, x::Vector{ComplexF64}) at apiutils.jl:42
vector_mode_jacobian!(result::Matrix{ComplexF64}, f!::SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, y::Vector{ComplexF64}, x::Vector{ComplexF64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2}}}}) at jacobian.jl:172
jacobian!(result::Matrix{ComplexF64}, f!::Function, y::Vector{ComplexF64}, x::Vector{ComplexF64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{ComplexF64}}, ComplexF64}, ComplexF64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Not...

Okay, so the first error is because ReverseDiff.jl cannot do complex numbers, and now this one is because ForwardDiff.jl cannot.

And you can always do numerical vjps with adjoints, which is equivalent to SUNDIALS but 😅

using OrdinaryDiffEq, Zygote, DiffEqSensitivity
function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0; 1.0]
prob = ODEProblem(fiip,complex(u0),(0.0,10.0),complex(p))

function sum_of_solution(u0, p)
    _prob = remake(prob,u0=u0,p=p)
    real(sum(solve(_prob,Tsit5(),reltol=1e-6,abstol=1e-6,saveat=0.1, sensealg=InterpolatingAdjoint(autodiff=false,autojacvec=false))))
end

dx = Zygote.gradient(sum_of_solution, complex(u0), complex(p)) # works!
(ComplexF64[-39.12610254292564 - 0.0im, -8.78792677656453 - 0.0im], ComplexF64[8.30761239779402 + 0.0im, -159.48459555162577 + 0.0im, 75.20354307391543 + 0.0im, -339.1934977276008 + 0.0im])

So okay, the solution is easy then: find out why we aren't defaulting to Enzyme here.

ChrisRackauckas avatar Dec 02 '21 02:12 ChrisRackauckas

Thanks. This found a nice solution of just fixing part of the defaulting system in https://github.com/SciML/DiffEqSensitivity.jl/pull/511 . Now if you don't specify any sensealgs it'll just work using a fast autodiff. I'll leave this open because it would be good to get other solutions, but ultimately getting those other ones working is an upstream problem.

ChrisRackauckas avatar Dec 02 '21 02:12 ChrisRackauckas