Problems with Parameter optimization of an ODE

[MaAl13]

Hi everybody, i already posted the same issue on the Optimization GitHub page, but Chris asked me to transfer it here. I get the following domain error while trying to fit the parameters of an ODE to artificial created data. The problem is that i always run into a DomainError. I already set the ODE values that are negative to zero and punish negative parameter values. Nevertheless, the error still occurs. Changing the algorithm to PolyOpt() helps in the way, that there are more iterations until the error pops up.

using DifferentialEquations, Optimization, OptimizationPolyalgorithms, OptimizationOptimJL,
      DiffEqSensitivity, Zygote, Plots
using Distributions
using DiffEqCallbacks
using OptimizationFlux


u0 = [0.501, 1.955, 0.198, 0.148, 0.161, 0.161, 0.064]
parameters_real = [2.5, 100.0, 6.0, 16.0, 100.0, 1.28, 12.0, 1.8, 13.0, 4.0, 0.52, 0.1, 1.0, 4.0]

function Glucose_Yeast2(du,u,p,t)
    du[1] = p[1] - p[2] * u[1] * u[6] / (1 + (u[6] / p[11])^p[10])
    du[2] = 2 * p[2] * u[1] * u[6] / (1 + (u[6] / p[11])^p[10]) - p[3]*u[2]*(p[13]-u[5]) - p[7]*u[2]*u[5]
    du[3] = p[3]*u[2]*(p[13]-u[5]) - p[4]*u[3]*(p[14] - u[6])
    du[4] = p[4]*u[3]*(p[14] - u[6]) - p[5]*u[4]*u[5] - p[9]*(u[4]-u[7])
    du[5] = p[3]*u[2]*(p[13]-u[5]) - p[5]*u[4]*u[5] - p[7] * u[2] * u[5]
    du[6] = -2 * p[2] * u[1] * u[6] / (1 + (u[6] / p[11])^p[10]) + 2 * p[4]*u[3]*(p[14] - u[6]) - p[6] * u[6]
    du[7] = p[12]*p[9] * (u[4] - u[7]) - p[9] * u[7]
end

time_points = [[0.25, 0.75, 1.15, 1.5, 1.75, 2.5, 3.25, 3.75, 4.95, 6.0, 7.0, 7.5, 8.5, 8.75, 9.9, 10.0],
[0.25, 0.75, 1.15, 1.5, 1.75, 2.5, 3.25, 3.75, 4.95, 6.0, 7.2, 7.5, 8.5, 8.75, 9.9],
[0.25, 0.75, 1.15, 1.5, 1.75, 2.5, 3.25, 3.75, 4.95, 6.0, 7.0, 7.5, 8.5, 8.75, 9.9, 10.0],
[0.25, 0.75, 1.15, 1.5, 1.75, 2.5, 3.25, 3.75, 4.95, 6.0, 7.0, 7.5, 8.5, 8.75, 9.9, 10.0],
[0.25, 0.75, 1.15, 1.65, 1.75, 2.5, 3.25, 3.75, 4.95, 6.0, 7.25, 7.55, 7.75, 7.9, 8.5, 8.75, 9.9, 10.0],
[0.25, 0.75, 1.15, 0.75, 1.75, 2.5, 3.25, 3.75, 5.0, 6.0, 7.0, 7.5, 8.6, 8.85, 9.9, 10.0],
[0.25, 0.75, 1.15, 1.45, 1.75, 2.5, 3.25, 3.5, 5.1, 6.0, 7.0, 7.5, 8.5, 8.75, 9.75, 10.0]]

tstart = 0
tend = 10


function t_points_creator(time_points)
    no_of_measurements = length(time_points)#length(string_measurements)
    union_dataset_t = []
    dataset_t_indices = []
    total_no_measurements = 0
    indices_in_stacked_form = []

    for i in 1:no_of_measurements
        tmp =  [total_no_measurements+1]
        total_no_measurements += length(time_points[i])
        tmp = push!(tmp, total_no_measurements)
        indices_in_stacked_form = push!(indices_in_stacked_form, tmp)
        union_dataset_t = union(union_dataset_t, time_points[i])
    end


    union_dataset_t = sort!(union_dataset_t)


    for i in 1:no_of_measurements
        tmp = []
        for j in 1:length(time_points[i])
            push!(tmp, findall(x->x==time_points[i][j], union_dataset_t)[1])
        end
        push!(dataset_t_indices, tmp)
    end
    return union_dataset_t, dataset_t_indices, total_no_measurements, indices_in_stacked_form
end

function sol_ODE(model, u0, tspan, model_params, saves_t)
  prob=ODEProblem(model,u0,tspan, model_params)
  return Array(solve(prob,Tsit5(),
              saveat=saves_t,abstol=1e-8,reltol=1e-6))
end





function data_creator(solution, noise_data, time_indices)
  solution_combined = [solution[1,:], solution[2,:], solution[3,:], solution[4,:], solution[5,:], solution[6,:], solution[7,:]]

  sol_noise = zeros((length(solution_combined), length(solution_combined[1])))
  for i in 1:length(solution_combined)
        tmp_sol = solution_combined[i][time_indices[i]]
        dim = length(tmp_sol)
        noise = zeros(dim)
        for j in 1:dim
            d = Normal(tmp_sol[j], maximum([noise_data*tmp_sol[j],1e-8]))
            
            # Lower limit must be >0 in order that the neuron still size_filters
            # and that they can be distinguished by  0 values, which come from padding or missing values
            # THIS MEANS THAT IF THE REAL DATA HAS 0 VALUES THEY MUST BE SHIFTED UPWARDS BY 1E-8
            td = truncated(d, 1e-8, Inf)
            noise[j] = rand(td, 1)[1]
        end
        sol_noise[i, time_indices[i]] = noise
  end
  return sol_noise
end



function loss_ODE(model, u0, tspan, model_params::Vector, indices_components, saves_t, data)

    # Using PositiveDomain callback in DiffEqCallbacks.jl, which is inspired by Shampine's et al. paper
    # about non-negative ODE solutions.

    #-----------------------------------------------
    # Check if parameter is complex

    prediction_raw = solve(ODEProblem(model,u0,tspan,model_params), Tsit5(), saveat=saves_t, isoutofdomain = (u,p,t)->any(x->x<0,u))
#    		   callback=PositiveDomain(ones(7)),verbose=false) 
    prediction = Array(prediction_raw)
    prediction_components = [prediction[1,:],prediction[2,:], prediction[3,:], 
    prediction[4,:], prediction[5,:], prediction[6,:], prediction[7,:]]
    
    if all(<(0), model_params)
        return 1e+30, prediction
    else
        loss = 0.0
        for i in 1:length(indices_components)
            prediction_component=prediction_components[i][indices_components[i]]
            data_component = data[i,:][indices_components[i]]
            prediction_component = prediction_component
            #println(prediction_component)
            loss += sum(abs2,log.(prediction_component.+1e-8) .- log.(data_component.+1e-8))
        end
        return loss, prediction
    end
end

callback = function (p, l, pred)
  println(l)
  return false
end


function train(model, u0, tspan, indices_components, saves_t, data, p_guess)
    dummy = zeros(20)
    adtype = Optimization.AutoZygote()
    optf = Optimization.OptimizationFunction((x,dummy)->loss_ODE(model, u0, tspan, x, indices_components, saves_t, data), adtype)
    optprob = Optimization.OptimizationProblem(optf, p_guess)

    result_ode = Optimization.solve(optprob, BFGS(),
                                  callback = callback,
                                  maxiters = 10000)
    println(result_ode)
end

saves_t, time_indices, total_no_measurements, indices_components = t_points_creator(time_points)
tspan = (tstart, tend)
data = sol_ODE(Glucose_Yeast2, u0, tspan, parameters_real, saves_t)
data_noise = data_creator(data, 0.1, time_indices)
# Simulation interval and intermediary points

p_guess = ones(14)
println(loss_ODE(Glucose_Yeast2, u0, tspan, p_guess, time_indices, saves_t, data_noise)[1])
train(Glucose_Yeast2, u0, tspan, time_indices, saves_t, data_noise, p_guess)

ERROR: LoadError: DomainError with -0.015165612843735379: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar. Stacktrace: [1] throw_exp_domainerror(x::Float64) @ Base.Math ./math.jl:37 [2] ^ @ ./math.jl:909 [inlined] [3] Glucose_Yeast2(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, t::Float64) @ Main ~/Documents/julia_GPU_solver/Test5.jl:12 [4] ODEFunction @ ~/.julia/packages/SciMLBase/UEAKN/src/scimlfunctions.jl:1613 [inlined] [5] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, var"#3#5", 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{}, Vector{Any}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/ZgJ9s/src/perform_step/low_order_rk_perform_step.jl:627 [6] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Vector{Any}, 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::var"#3#5", 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.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/ZgJ9s/src/solve.jl:456 [7] #__solve#502 @ ~/.julia/packages/OrdinaryDiffEq/ZgJ9s/src/solve.jl:4 [inlined] [8] #solve_call#28 @ ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:428 [inlined] [9] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :isoutofdomain), Tuple{Vector{Any}, var"#3#5"}}}) @ DiffEqBase ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:726 [10] #solve#29 @ ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:710 [inlined] [11] _concrete_solve_adjoint(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ::ForwardDiffSensitivity{0, nothing}, ::Vector{Float64}, ::Vector{Float64}; saveat::Vector{Any}, kwargs::Base.Pairs{Symbol, var"#3#5", Tuple{Symbol}, NamedTuple{(:isoutofdomain,), Tuple{var"#3#5"}}}) @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/60SaR/src/concrete_solve.jl:362 [12] _concrete_solve_adjoint(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ::Nothing, ::Vector{Float64}, ::Vector{Float64}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :isoutofdomain), Tuple{Vector{Any}, var"#3#5"}}}) @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/60SaR/src/concrete_solve.jl:63 [13] _solve_adjoint(prob::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; merge_callbacks::Bool, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :isoutofdomain), Tuple{Vector{Any}, var"#3#5"}}}) @ DiffEqBase ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:1069 [14] #rrule#50 @ ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:1032 [inlined] [15] chain_rrule_kw @ ~/.julia/packages/Zygote/DkIUK/src/compiler/chainrules.jl:229 [inlined] [16] macro expansion @ ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [inlined] [17] _pullback @ ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:9 [inlined] [18] _apply @ ./boot.jl:814 [inlined] [19] adjoint @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:204 [inlined] [20] _pullback @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined] [21] _pullback @ ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:710 [inlined] [22] _pullback(::Zygote.Context, ::DiffEqBase.var"##solve#29", ::Nothing, ::Nothing, ::Nothing, ::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :isoutofdomain), Tuple{Vector{Any}, var"#3#5"}}}, ::typeof(solve), ::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [23] _apply(::Function, ::Vararg{Any}) @ Core ./boot.jl:814 [24] adjoint @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:204 [inlined] [25] _pullback @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined] [26] _pullback @ ~/.julia/packages/DiffEqBase/hZncn/src/solve.jl:705 [inlined] [27] _pullback(::Zygote.Context, ::CommonSolve.var"#solve##kw", ::NamedTuple{(:saveat, :isoutofdomain), Tuple{Vector{Any}, var"#3#5"}}, ::typeof(solve), ::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, ODEFunction{true, typeof(Glucose_Yeast2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [28] _pullback @ ~/Documents/julia_GPU_solver/Test5.jl:104 [inlined] [29] _pullback(::Zygote.Context, ::typeof(loss_ODE), ::typeof(Glucose_Yeast2), ::Vector{Float64}, ::Tuple{Int64, Int64}, ::Vector{Float64}, ::Vector{Any}, ::Vector{Any}, ::Matrix{Float64}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [30] _pullback @ ~/Documents/julia_GPU_solver/Test5.jl:134 [inlined] [31] _pullback(::Zygote.Context, ::var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, ::Vector{Float64}, ::SciMLBase.NullParameters) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [32] _apply(::Function, ::Vararg{Any}) @ Core ./boot.jl:814 [33] adjoint @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:204 [inlined] [34] _pullback @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined] [35] _pullback @ ~/.julia/packages/SciMLBase/UEAKN/src/scimlfunctions.jl:2697 [inlined] [36] _pullback(::Zygote.Context, ::OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, ::Vector{Float64}, ::SciMLBase.NullParameters) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [37] _apply(::Function, ::Vararg{Any}) @ Core ./boot.jl:814 [38] adjoint @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:204 [inlined] [39] _pullback @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined] [40] _pullback @ ~/.julia/packages/Optimization/RUgSr/src/function/zygote.jl:30 [inlined] [41] _pullback(ctx::Zygote.Context, f::Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}, args::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [42] _apply(::Function, ::Vararg{Any}) @ Core ./boot.jl:814 [43] adjoint @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:204 [inlined] [44] _pullback @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined] [45] _pullback @ ~/.julia/packages/Optimization/RUgSr/src/function/zygote.jl:32 [inlined] [46] _pullback(ctx::Zygote.Context, f::Optimization.var"#131#141"{Tuple{}, Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}}, args::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0 [47] _pullback(f::Function, args::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:34 [48] pullback(f::Function, args::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:40 [49] gradient(f::Function, args::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:75 [50] (::Optimization.var"#129#139"{Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}})(::Vector{Float64}, ::Vector{Float64}) @ Optimization ~/.julia/packages/Optimization/RUgSr/src/function/zygote.jl:32 [51] (::OptimizationOptimJL.var"#5#13"{OptimizationProblem{true, OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, OptimizationOptimJL.var"#4#12"{OptimizationProblem{true, OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#129#139"{Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}}, Optimization.var"#132#142"{Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}}, Optimization.var"#137#147", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}}, OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#129#139"{Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}}, Optimization.var"#132#142"{Optimization.var"#128#138"{OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, SciMLBase.NullParameters}}, Optimization.var"#137#147", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}})(G::Vector{Float64}, θ::Vector{Float64}) @ OptimizationOptimJL ~/.julia/packages/OptimizationOptimJL/fdrJg/src/OptimizationOptimJL.jl:106 [52] value_gradient!!(obj::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, x::Vector{Float64}) @ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/interface.jl:82 [53] value_gradient!(obj::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, x::Vector{Float64}) @ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/interface.jl:69 [54] value_gradient!(obj::Optim.ManifoldObjective{TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}}, x::Vector{Float64}) @ Optim ~/.julia/packages/Optim/6Lpjy/src/Manifolds.jl:50 [55] (::LineSearches.var"#ϕdϕ#6"{Optim.ManifoldObjective{TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}})(α::Float64) @ LineSearches ~/.julia/packages/LineSearches/Ki4c5/src/LineSearches.jl:84 [56] (::LineSearches.HagerZhang{Float64, Base.RefValue{Bool}})(ϕ::Function, ϕdϕ::LineSearches.var"#ϕdϕ#6"{Optim.ManifoldObjective{TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, c::Float64, phi_0::Float64, dphi_0::Float64) @ LineSearches ~/.julia/packages/LineSearches/Ki4c5/src/hagerzhang.jl:139 [57] HagerZhang @ ~/.julia/packages/LineSearches/Ki4c5/src/hagerzhang.jl:101 [inlined] [58] perform_linesearch!(state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}, method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}, d::Optim.ManifoldObjective{TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}}) @ Optim ~/.julia/packages/Optim/6Lpjy/src/utilities/perform_linesearch.jl:59 [59] update_state!(d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}, method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}) @ Optim ~/.julia/packages/Optim/6Lpjy/src/multivariate/solvers/first_order/bfgs.jl:139 [60] optimize(d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}, options::Optim.Options{Float64, OptimizationOptimJL.var"#_cb#11"{var"#7#8", BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}, Base.Iterators.Cycle{Tuple{Optimization.NullData}}}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}) @ Optim ~/.julia/packages/Optim/6Lpjy/src/multivariate/optimize/optimize.jl:54 [61] optimize(d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}, options::Optim.Options{Float64, OptimizationOptimJL.var"#_cb#11"{var"#7#8", BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}, Base.Iterators.Cycle{Tuple{Optimization.NullData}}}}) @ Optim ~/.julia/packages/Optim/6Lpjy/src/multivariate/optimize/optimize.jl:36 [62] ___solve(prob::OptimizationProblem{true, OptimizationFunction{true, Optimization.AutoZygote, var"#9#10"{typeof(Glucose_Yeast2), Vector{Float64}, Tuple{Int64, Int64}, Vector{Any}, Vector{Any}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, opt::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Flat}, data::Base.Iterators.Cycle{Tuple{Optimization.NullData}}; callback::Function, maxiters::Int64, maxtime::Nothing, abstol::Nothing, reltol::Nothing, progress::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ OptimizationOptimJL ~/.julia/packages/OptimizationOptimJL/fdrJg/src/OptimizationOptimJL.jl:143 [63] #__solve#2 @ ~/.julia/packages/OptimizationOptimJL/fdrJg/src/OptimizationOptimJL.jl:56 [inlined] [64] #solve#492 @ ~/.julia/packages/SciMLBase/UEAKN/src/solve.jl:56 [inlined] [65] train(model::Function, u0::Vector{Float64}, tspan::Tuple{Int64, Int64}, indices_components::Vector{Any}, saves_t::Vector{Any}, data::Matrix{Float64}, p_guess::Vector{Float64}) @ Main ~/Documents/julia_GPU_solver/Test5.jl:137 in expression starting at /Users/malmansto/Documents/julia_GPU_solver/Test5.jl:151

[ChrisRackauckas]

Did you try decreasing the tolerances? Did you try isoutofdomain? I'd try those first.

Are you sure this ODE has a guaranteed positive solution? It doesn't look like it would for arbitrary parameters. Since that's the case, you'd need to change x^y to abs(x)^y for this to be a well-defined model, since otherwise true negative values would just error.

[ChrisRackauckas]

This ended up just being a tolerance thing.