search

SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.87k stars 230 forks source link

Parameter estimation with callbacks? #425

Closed ayadlin closed 5 years ago

ayadlin commented 5 years ago

I've been looking for a while, but i am not sure if there are any example or anyway to do parameter estimation in a problem with callbacks - The parameter estimation takes the problem, but the callbacks are added in the solve stage. I've been playing with your conditional dose example - Is there anyway around this issue to recover the correct parameters? please let me know if this is a legitimate question, I am new to Julia and ODE programing in general. thanks Here is the code I use:

using DifferentialEquations
using DiffEqParamEstim
using ParameterizedFunctions
using Plots
using Measurements
using RecursiveArrayTools

# Equation and solution with callback
test = @ode_def begin
    dx=-k*x
    end k
x0 = [10.0 ± 1]
const V = 1
k=1

prob = ODEProblem(test,x0,(0.0,10.0),k)

condition(u,t,integrator) = t in 2:2:10# && u[1]/V<1
affect!(integrator) = integrator.u[1] += 10 ± 1

cb = DiscreteCallback(condition,affect!)
sol = solve(jump_prob,Tsit5(),tstops=2:2:10,callback=cb)

mean = [x[1].val for x in sol.u]
band=[x[1].err for x in sol.u]
plot(sol.t,mean, ribbon=band, alpha=0.5, width=3)

# simulate Data
t = collect(range(0,stop=10,length=200))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i])) for i in 1:length(t)])
data = convert(Array,randomized)
plot(sol.t,mean, ribbon=band, alpha=0.5, width=3)
scatter!(t,data')

# CostFunction
#--> this is the gist of the issue, how can the callback be included in the cost function?

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false)

# Checking parameters for cost function:

vals =collect(range(0.0±0,stop=1.1,length=100))
using Plots; plotly()
vals
[cost_function(i) for i in vals]
plot(vals,[cost_function(i) for i in vals],yscale=:log10,
      xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
      lw = 3)

The minimal is found around 0.15, while the real parameter is 1 - this is happening because the cost function doesn't know about the callbacks if I run same code without callback in the solution the minimal is found correctly at 1 thanks for your help, please let me know if this is inappropriate and I will not bother you again! Ale

ayadlin commented 5 years ago

Sorry, i left the # and know this looks awful

ChrisRackauckas commented 5 years ago

You didn't pass the callbacks or the tstops to the cost function?

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                                     callback = cb, tstops = 2:2:10,
                                     maxiters=10000,verbose=false)

all it's doing is feeding those args to the internal solve.

ayadlin commented 5 years ago

Thanks for getting back to me, and doing so qucikly - I passed the cost function as you suggested , but this result in an error when I run the next step of the code (attached below), the error doesn't happen if I leave the callback out - thats one of the reasons I am confused about how to deal with the parameter estimation :

Thanks

# Code that fails when call back is included in cost_function
# Checking parameters for cost function:

vals =collect(range(0.0±0,stop=1.1,length=100))
using Plots; plotly()
vals
[cost_function(i) for i in vals]
plot(vals,[cost_function(i) for i in vals],yscale=:log10,
      xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
      lw = 3)

UndefRefError: access to undefined reference (StackTrace belolw) - The alternative result is that my kernel (in jupyter or VSCode crashes, although I am never close to 100% usage of CPU or RAM (<= 50%)

Stacktrace: [1] getindex at .\array.jl:730 [inlined] [2] getindex at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\adjtrans.jl:130 [inlined] [3] (::L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing})(::ODESolution{Measurement{Float64},2,Array{Array{Measurement{Float64},1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Measurement{Float64},1},1},1},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Measurement{Float64},getfield(Main, Symbol("##399")){getfield(Main, Symbol("##223#227")),getfield(Main, Symbol("##224#228")),getfield(Main, Symbol("##225#229")),Nothing,Nothing,getfield(Main, Symbol("##226#230")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{getfield(Main, Symbol("##399")){getfield(Main, Symbol("##223#227")),getfield(Main, Symbol("##224#228")),getfield(Main, Symbol("##225#229")),Nothing,Nothing,getfield(Main, Symbol("##226#230")),Expr,Expr},Array{Array{Measurement{Float64},1},1},Array{Float64,1},Array{Array{Array{Measurement{Float64},1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Measurement{Float64},1},Array{Measurement{Float64},1},Array{Measurement{Float64},1},OrdinaryDiffEq.Tsit5ConstantCache{Measurement{Float64},Float64}}},DiffEqBase.DEStats{Nothing}}) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\cost_functions.jl:49 [4] (::getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:callback, :tstops, :maxiters, :verbose),Tuple{DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},StepRange{Int64,Int64},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##399")){getfield(Main, Symbol("##223#227")),getfield(Main, Symbol("##224#228")),getfield(Main, Symbol("##225#229")),Nothing,Nothing,getfield(Main, Symbol("##226#230")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing})(::Measurement{Float64}) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\build_loss_objective.jl:29 [5] DiffEqObjective at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\build_loss_objective.jl:8 [inlined] [6] iterate at .\generator.jl:47 [inlined] [7] collect(::Base.Generator{Array{Measurement{Float64},1},DiffEqObjective{getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:callback, :tstops, :maxiters, :verbose),Tuple{DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},StepRange{Int64,Int64},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##399")){getfield(Main, Symbol("##223#227")),getfield(Main, Symbol("##224#228")),getfield(Main, Symbol("##225#229")),Nothing,Nothing,getfield(Main, Symbol("##226#230")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing},getfield(DiffEqParamEstim, Symbol("##17#21"))}}) at .\array.jl:606

ayadlin commented 5 years ago

Doing some extra research I rewrote the code without the using Measurements: and it works well with the callbacks - I can't figure out why it would work in one case and not the other, but maybe you have some comments or ideas what need to be fixed to work with Measurements?

Thanks

using Pkg 
using DifferentialEquations
using DiffEqParamEstim
using ParameterizedFunctions
using Plots
using DiffEqBiological
using Latexify
using BenchmarkTools
using Measurements
using StaticArrays

# using DiffEqCallbacks
# using ParameterizedFunctions
test = @ode_def begin
    dx=-k*x
    end k
x0 = [10.0]# ± 1]
const V = 1
k=1
prob = ODEProblem(test,x0,(0.0,10.0),k)

condition(u,t,integrator) = t in [2,4,6,8] # && u[1]/V<1
affect!(integrator) = integrator.u[1] += 10 #±1

cb = DiscreteCallback(condition,affect!)
sol = solve(prob,Tsit5(),callback=cb,tstops=[2,4,6,8])

# mean = [x[1].val for x in sol.u]
# band=[x[1].err for x in sol.u]

#plot(sol.t,mean, ribbon=band, alpha=0.5, width=3)
plot(sol, width=3)

t = collect(range(0,stop=10,length=1000))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]).+rand(-1:0.01:1)[1]) for i in 1:length(t)])
data = convert(Array,randomized)
#plot(sol.t,mean, ribbon=band, alpha=0.5, width=3)
plot(sol)
scatter!(t,data')

#t= collect(range(0,stop=10,length=20))
cost_function = build_loss_objective(prob,Tsit5(),
     tstops = [2,4,6,8], callback=cb,L2Loss(t,data'), maxiters=100000,verbose=false)

vals =collect(range(0,stop=2,length=1000))
using Plots; plotly()
vals
#[cost_function(i) for i in vals]
plot(vals,[cost_function(i) for i in vals],yscale=:log10,
       xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
       lw = 3)
ChrisRackauckas commented 5 years ago

What code is being ran and what is the full error message?

ayadlin commented 5 years ago

Here is the code - it the one I sent you before with your corrected cost function:


using Pkg 
using DifferentialEquations
using DiffEqParamEstim
using ParameterizedFunctions
using Plots
using DiffEqBiological
using Latexify
using BenchmarkTools
using Measurements
using StaticArrays

# Equation
test = @ode_def begin
    dx=-k*x
    end k
x0 = [10.0 ± 1]
const V = 1
k=1

# Solve Equation
prob = ODEProblem(test,x0,(0.0,10.0),k)
condition(u,t,integrator) = t in [2,4,6,8] # && u[1]/V<1
affect!(integrator) = integrator.u[1] += 10±1
cb = DiscreteCallback(condition,affect!)
sol = solve(prob,Tsit5(),callback=cb,tstops=[2,4,6,8])
plot(sol, width=3)

# Simulate Data
t = collect(range(0,stop=10,length=500))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]).+rand(-2:0.01:2)[1]) for i in 1:length(t)])
data = convert(Array,randomized)
plot(sol)
scatter!(t,data')

# Objective Function
cost_function = build_loss_objective(prob,Tsit5(),
     tstops = [2,4,6,8], callback=cb,L2Loss(t,data'), maxiters=100000,verbose=false)

# Plot cost function Float64
vals =collect(range(0,stop=2,length=1000))
plot(vals,[cost_function(i) for i in vals],yscale=:log10,
       xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
       lw = 3)

# Plot cost function Measurement
vals =collect(range(0±0,stop=2±0,length=1000))
plot(vals,[cost_function(i) for i in vals],yscale=:log10,
       xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
       lw = 3)

Everything runs well until # Plot cost function section I run this section using Float64 or Measurement type for vals.

Full Error Message if run Plot cost function Float64:

MethodError: no method matching Float64(::Measurement{Float64}) Closest candidates are: Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194 Float64(::T<:Number) where T<:Number at boot.jl:741 Float64(!Matched::Int8) at float.jl:60 ...

Stacktrace: [1] _broadcast_getindex_evalf at .\broadcast.jl:578 [inlined] [2] _broadcast_getindex at .\broadcast.jl:551 [inlined] [3] getindex at .\broadcast.jl:511 [inlined] [4] copy at .\broadcast.jl:787 [inlined] [5] materialize at .\broadcast.jl:753 [inlined] [6] STANDARD_PROB_GENERATOR(::ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Float64) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\DiffEqParamEstim.jl:5 [7] (::getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:tstops, :callback, :maxiters, :verbose),Tuple{Array{Int64,1},DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing})(::Float64) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\build_loss_objective.jl:22 [8] (::DiffEqObjective{getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:tstops, :callback, :maxiters, :verbose),Tuple{Array{Int64,1},DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing},getfield(DiffEqParamEstim, Symbol("##17#21"))})(::Float64) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\build_loss_objective.jl:8 [9] iterate at .\generator.jl:47 [inlined] [10] collect(::Base.Generator{Array{Float64,1},DiffEqObjective{getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:tstops, :callback, :maxiters, :verbose),Tuple{Array{Int64,1},DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing},getfield(DiffEqParamEstim, Symbol("##17#21"))}}) at .\array.jl:606

Full Error Message if run Plot cost function Measurement instead:

UndefRefError: access to undefined reference

Stacktrace: [1] getindex at .\array.jl:730 [inlined] [2] getindex at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\adjtrans.jl:130 [inlined] [3] (::L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing})(::ODESolution{Measurement{Float64},2,Array{Array{Measurement{Float64},1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Measurement{Float64},1},1},1},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Measurement{Float64},getfield(Main, Symbol("##365")){getfield(Main, Symbol("##15#19")),getfield(Main, Symbol("##16#20")),getfield(Main, Symbol("##17#21")),Nothing,Nothing,getfield(Main, Symbol("##18#22")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{getfield(Main, Symbol("##365")){getfield(Main, Symbol("##15#19")),getfield(Main, Symbol("##16#20")),getfield(Main, Symbol("##17#21")),Nothing,Nothing,getfield(Main, Symbol("##18#22")),Expr,Expr},Array{Array{Measurement{Float64},1},1},Array{Float64,1},Array{Array{Array{Measurement{Float64},1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Measurement{Float64},1},Array{Measurement{Float64},1},Array{Measurement{Float64},1},OrdinaryDiffEq.Tsit5ConstantCache{Measurement{Float64},Float64}}},DiffEqBase.DEStats{Nothing}}) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\cost_functions.jl:49 [4] (::getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:tstops, :callback, :maxiters, :verbose),Tuple{Array{Int64,1},DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##365")){getfield(Main, Symbol("##15#19")),getfield(Main, Symbol("##16#20")),getfield(Main, Symbol("##17#21")),Nothing,Nothing,getfield(Main, Symbol("##18#22")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing})(::Measurement{Float64}) at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\build_loss_objective.jl:29 [5] DiffEqObjective at C:\Users\awolf-yadlin.julia\packages\DiffEqParamEstim\Rdg1H\src\build_loss_objective.jl:8 [inlined] [6] iterate at .\generator.jl:47 [inlined] [7] collect(::Base.Generator{Array{Measurement{Float64},1},DiffEqObjective{getfield(DiffEqParamEstim, Symbol("##14#18")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:tstops, :callback, :maxiters, :verbose),Tuple{Array{Int64,1},DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)},Int64,Bool}}},ODEProblem{Array{Measurement{Float64},1},Tuple{Float64,Float64},true,Int64,getfield(Main, Symbol("##365")){getfield(Main, Symbol("##15#19")),getfield(Main, Symbol("##16#20")),getfield(Main, Symbol("##17#21")),Nothing,Nothing,getfield(Main, Symbol("##18#22")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},LinearAlgebra.Adjoint{Measurement{Float64},Array{Measurement{Float64},2}},Nothing,Nothing},Nothing},getfield(DiffEqParamEstim, Symbol("##17#21"))}}) at .\array.jl:606

thanks-

A

ayadlin commented 5 years ago

a second question related to callbacks is: Do discrete callbacks work with DAE? take this code for example:


function f(out,du,u,p,t)
    k₁₂, k₂₁, kdeg =p 
    X1=u[1]
    X2=u[2]
    X3=u[3]
    out[1] = -(k₁₂+kdeg)*X1+k₂₁*X2 - du[1]
    out[2] = k₁₂*X1-k₂₁*X2 - du[2]
    out[3] = u[1] - u[2] - u[3]
end

u₀ = [10.0, 0, 0]
du₀ = [0, 0, 1.0]
tspan = (0.0,10.0)
p=[1,0,0]
differential_vars = [true,true, false]
prob = DAEProblem(f,du₀,u₀,tspan,p,differential_vars=differential_vars)

using Sundials
#sol = solve(prob,IDA(linear_solver=:GMRES))

condition(u,t,integrator) = t in [2,4,6,8] # && u[1]/V<1
affect!(integrator) = integrator.u[1] +=10# 10#±1
# affect!(integrator) = integrator.u[1] *=1.01

cb = DiscreteCallback(condition,affect!)
sol = solve(prob,IDA(linear_solver=:GMRES),callback=cb,tstops=[2,4,6,8],alg_hints=[:stiff])

if i solve without callbacks using commented code there is no issue whatsoever:

sol = solve(prob,IDA(linear_solver=:GMRES))

Also if I use callbacks that don't change much (commented code too) there is also not problems

affect!(integrator) = integrator.u[1] *=1.01

if i try to use callbacks that actually modify u[1] (code exactly as above): affect!(integrator) = integrator.u[1] +=10 I get the following: "[IDAS ERROR] IDASolve At t = 2 and h = 4.88281e-007, the error test failed repeatedly or with |h| = hmin."

Is this because by adding the discrete callback I am jumping to far from 0 in out[1]? is there a way around to make the callbacks work with DAEProbelms when you add big jumps?

Thanks,

A

PS- sorry for all the questions and issues - I hope some of this are helpful things and not me just being :Dense

ChrisRackauckas commented 5 years ago

Do discrete callbacks work with DAE?

Yes, the issue here was actually quite subtle. Since you changed u, the residual of f(res,du,u,p,t) was no longer zero, so the new initial condition was not consistent. However, requiring you to change du seems not nice, so what I did was added a consistency change and IC finding step after callback modifications, and released this as Sundials v3.1.0. That fixes this problem and so when you update it should be good.

ayadlin commented 5 years ago

Thanks so much - I am very appreciative of how quickly you respond to questions and fix issues! Going to give it a try tonight :) take care,

A

ChrisRackauckas commented 5 years ago

Here's a version of the cost function that works with measurement types:

using OrdinaryDiffEq
using DiffEqParamEstim
using ParameterizedFunctions
using Measurements

# Equation
test = @ode_def begin
    dx=-k*x
    end k
x0 = [10.0 ± 1]
const V = 1
k=1

# Solve Equation
prob = ODEProblem(test,x0,(0.0,10.0),k)
condition(u,t,integrator) = t in [2,4,6,8] # && u[1]/V<1
affect!(integrator) = integrator.u[1] += 10±1
cb = DiscreteCallback(condition,affect!;save_positions=(false,false))
sol = solve(prob,Tsit5(),callback=cb,tstops=[2,4,6,8])

# Simulate Data
t = collect(range(0,stop=10,length=500))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]).+rand(-2:0.01:2)[1]) for i in 1:length(t)])
data = convert(Array,randomized)

# Objective Function
cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data'),
     tstops = [2,4,6,8], callback=cb, maxiters=100000,verbose=false)

# Plot cost function Measurement
vals =collect(range(0±0,stop=2±0,length=1000))
cost_function(vals[1])

Note that p has to be a measurement if u0 is a measurement because it will auto-uptype u0 to match p in order to handle dual numbers for automatic differentiation. But it will downtype u0 back to a float if p is a float, or at least it tries to and recognize it can't which gives the error. This is then fixed by matching types on u0 and p. The output of the cost function is then an uncertainty type, so I don't know how you want to plot it or handle it, but it's interesting (maybe weigh by the amount of uncertainty?).

I think I answered all of the questions so I'll close this. Generally it's better to ask help questions on https://discourse.julialang.org/ or our chat channel. The issues are mainly for development issues but I'm not overly strict about that.

ayadlin commented 5 years ago

I'll make sure to address questions through the correct channel from now on- Again thanks for helping me,

Ale