# Creating a loss function for sciml\_train with unequal time point distribution of the ODE components

**URL:** <https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204>\
**Category:** General Usage\
**Tags:** question, diffeqflux\
**Created:** [January 25, 2022, 9:42pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204 "2022-01-25T21:42:50Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [January 25, 2022, 9:42pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/1 "2022-01-25T21:42:50Z")

</div>

Hello,

i am wondering how i can build an apropriate loss function for parameter estimation with DiffEqFlux.sciml\_train if i have data points from the ODE at different time points in each component. Can maybe someone help me out? Also i get a ArgumentError whenever i call sciml\_train, my current version of DiffEqFlux is 1.39.0.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 26, 2022, 9:49am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/2 "2022-01-26T09:49:12Z")

</div>

> [@Marius\_123](#):
>
> i am wondering how i can build an apropriate loss function for parameter estimation with DiffEqFlux.sciml\_train if i have data points from the ODE at different time points in each component. Can maybe someone help me out? Also i get a MethodError whenever i call sciml\_train, my current version of DiffEqFlux is 1.30.0.

`saveat` the union and then parse it out after the solve.

---

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [January 26, 2022, 9:59am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/3 "2022-01-26T09:59:34Z")

</div>

Thank you for you response! How would you parse out the soultion vector exactly in order to not come in any conflict with AD?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 26, 2022, 10:02am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/4 "2022-01-26T10:02:18Z")

</div>

`sol[:,[1,4,5]]` etc.

---

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [January 26, 2022, 10:17am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/5 "2022-01-26T10:17:22Z")

</div>

> [@ChrisRackauckas](#):
>
> `sol[:,[1,4,5]]`

So what i am doing at the moment is kind of the same. I iterate through all components of sol, pick the indices in each component where i have my measurements (tmp\_sol = sol[i,[1,2,3,5,6,7,8,9]]) and stack all so obtained tmp\_sol together with vstack to one huge 1-D vector. This is then subtracted to the stacked vector with real mesurements and the mse is calculated. However, this didn’t really work in a previous version of DiffEqFlux. Now with an upgraded version (1.39.0) i seem to not be able to run the copy paste code in the documentation, because sciml\_train produces the error ArgumentError: tuple must be non-empty, so i can’t really test it with sciml\_train. The “old method” that is running for me is with Flux.train! (This is shown in the video in the Optimization of ODEs article in the DiffEqFlux documentation.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 26, 2022, 10:51am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/6 "2022-01-26T10:51:15Z")

</div>

> [@Marius\_123](#):
>
> ArgumentError: tuple must be non-empty, so i can’t really test it with sciml\_train

Can you share code?

---

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [January 26, 2022, 11:22am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/7 "2022-01-26T11:22:59Z")

</div>

```julia
using DifferentialEquations, DiffEqFlux, Plots

function lotka_volterra!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

# Initial condition
u0 = [1.0, 1.0]

# Simulation interval and intermediary points
tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

# LV equation parameter. p = [α, β, δ, γ]
p = [1.5, 1.0, 3.0, 1.0]

# Setup the ODE problem, then solve
prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob, Tsit5())

println(sol[:,[1,4,5]])

# Plot the solution
using Plots
plot(sol)
savefig("LV_ode.png")

function loss(p)
  sol = solve(prob, Tsit5(), p=p, saveat = tsteps)
  loss = sum(abs2, sol.-1)
  return loss, sol
end

callback = function (p, l, pred)
  display(l)
  plt = plot(pred, ylim = (0, 6))
  display(plt)
  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  return false
end

result_ode = DiffEqFlux.sciml_train(loss, p)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 26, 2022, 11:27am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/8 "2022-01-26T11:27:54Z")

</div>

That loss function looks fine but that’s not doing what you say in the title. What are you actually trying to use/do?

---

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [January 26, 2022, 1:49pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/9 "2022-01-26T13:49:50Z")

</div>

This is not the actual loss function i have in my problem, just a copy of the code that is in the documentation and not running. My actual loss function has the following shape

```julia
function data_creator(sol, noise_data)
	combined_sol_vector = []
	GLIs_noise = Array{Float64}[]
	for i in 1:15
		if i<=3
			tmp_sol = sol[i,[1,2,3,5,6,7,8,9]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		elseif i >=7
			tmp_sol = sol[i,[1,2,3,5,6,7,8,9]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		else
			if i == 4
				tmp_sol = sol[i,[4,5,6]]
				for j in 1:size(tmp_sol)[1]
					tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
				end
				push!(GLIs_noise, tmp_sol)
			elseif i == 5
				tmp_sol = sol[i,[4,5,6]]
				for j in 1:size(tmp_sol)[1]
					tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
				end
				push!(GLIs_noise, tmp_sol)
			elseif i == 6
				tmp_sol = sol[i,[4,5,6]]
				for j in 1:size(tmp_sol)[1]
					tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
				end
				push!(GLIs_noise, tmp_sol)
				combined_sol_vector = vcat(combined_sol_vector, (GLIs_noise[1]+GLIs_noise[2])./GLIs_noise[3])
			end
		end
	end

	for i in 1:15
		if i<=3
			tmp_sol = sol[i,[4,6]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		elseif i >=7
			tmp_sol = sol[i,[4,6]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		end
	end

	for i in 16:30
		if i <=18
			tmp_sol = sol[i,[1,2,3,5,6,7,8,9]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		elseif i >=22
			tmp_sol = sol[i,[1,2,3,5,6,7,8,9]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		end
	end

	for i in 31:45
		if i <=33
			tmp_sol = sol[i,[4,6]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		elseif i >=37
			tmp_sol = sol[i,[4,6]]
			for j in 1:size(tmp_sol)[1]
				tmp_sol[j] = rand(Normal(tmp_sol[j], abs(tmp_sol[j]*noise_data)), 1)[1]
			end
			combined_sol_vector = vcat(combined_sol_vector, tmp_sol)
		end
	end
	return combined_sol_vector
end

function loss(p)
  sol = solve(prob, Tsit5(), p=p, saveat = tsteps)
	sol_flatten = data_creator(sol, noise_data)
  loss = sum(abs2, sol.-real_measurements)
  return loss, sol
end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 26, 2022, 1:53pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/10 "2022-01-26T13:53:24Z")

</div>

I’m confused. That will obviously hit a mutation error. What is the code you are actually running that reproduces your error and point to the line that has the error with the full stack trace.

---

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [January 26, 2022, 2:24pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/11 "2022-01-26T14:24:33Z")

</div>

So i am sorry if i have not been clear in my explanations. In the end i have two problems.

The first problem is that the example code from the documentation is not running, the sciml\_train gives me the error ArgumentError: tuple must be non-empty  
first(::Tuple{}) at tuple.jl:95  
\_unapply(::Nothing, ::Tuple{}) at lib.jl:163  
\_unapply(::Tuple{Nothing}, ::Tuple{}) at lib.jl:167  
\_unapply(::Tuple{Tuple{Nothing}}, ::Tuple{}) at lib.jl:167  
\_unapply(::Tuple{NTuple{6,Nothing},Tuple{Nothing}}, ::Tuple{Nothing,Nothing,Nothing,Array{Float64,1},Array{Float64,1},Nothing}) at lib.jl:168  
unapply(::Tuple{NTuple{6,Nothing},Tuple{Nothing}}, ::Tuple{Nothing,Nothing,Nothing,Array{Float64,1},Array{Float64,1},Nothing}) at lib.jl:177  
(::Zygote.var"#188#189"{Zygote.var"#kw\_zpullback#40"{DiffEqSensitivity.var"#adjoint\_sensitivity\_backpass#177"{Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Tsit5,InterpolatingAdjoint{0,true,Val{:central},Bool,Bool},Array{Float64,1},Array{Float64,1},Tuple{},NamedTuple{(),Tuple{}},Colon}},Tuple{NTuple{6,Nothing},Tuple{Nothing}}})(::VectorOfArray{Float64,2,Array{Array{Float64,1},1}}) at lib.jl:195  
(::Zygote.var"#1709#back#190"{Zygote.var"#188#189"{Zygote.var"#kw\_zpullback#40"{DiffEqSensitivity.var"#adjoint\_sensitivity\_backpass#177"{Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Tsit5,InterpolatingAdjoint{0,true,Val{:central},Bool,Bool},Array{Float64,1},Array{Float64,1},Tuple{},NamedTuple{(),Tuple{}},Colon}},Tuple{NTuple{6,Nothing},Tuple{Nothing}}}})(::VectorOfArray{Float64,2,Array{Array{Float64,1},1}}) at adjoint.jl:67  
#solve#57 at solve.jl:70 [inlined]  
(::typeof(∂(#solve#57)))(::VectorOfArray{Float64,2,Array{Array{Float64,1},1}}) at interface2.jl:0  
(::Zygote.var"#188#189"{typeof(∂(#solve#57)),Tuple{NTuple{6,Nothing},Tuple{Nothing}}})(::VectorOfArray{Float64,2,Array{Array{Float64,1},1}}) at lib.jl:194  
(::Zygote.var"#1709#back#190"{Zygote.var"#188#189"{typeof(∂(#solve#57)),Tuple{NTuple{6,Nothing},Tuple{Nothing}}}})(::VectorOfArray{Float64,2,Array{Array{Float64,1},1}}) at adjoint.jl:67  
solve at solve.jl:68 [inlined]  
(::typeof(∂(solve##kw)))(::VectorOfArray{Float64,2,Array{Array{Float64,1},1}}) at interface2.jl:0  
loss at Test\_AD.jl:32 [inlined]  
(::typeof(∂(loss)))(::Tuple{Float64,Nothing}) at interface2.jl:0  
#73 at train.jl:86 [inlined]  
#188 at lib.jl:194 [inlined]  
#1709#back at adjoint.jl:67 [inlined]  
OptimizationFunction at basic\_problems.jl:107 [inlined]  
#188 at lib.jl:194 [inlined]  
(::Zygote.var"#1709#back#190"{Zygote.var"#188#189"{typeof(∂(λ)),Tuple{Tuple{Nothing,Nothing},Int64}}})(::Tuple{Float64,Nothing}) at adjoint.jl:67  
#8 at solve.jl:91 [inlined]  
(::typeof(∂(λ)))(::Float64) at interface2.jl:0  
(::Zygote.var"#69#70"{Params,Zygote.Context,typeof(∂(λ))})(::Float64) at interface.jl:255  
gradient(::Function, ::Params) at interface.jl:59  
\_\_solve(::OptimizationProblem{true,OptimizationFunction{true,GalacticOptim.AutoForwardDiff,DiffEqFlux.var"#73#78"{typeof(loss)},GalacticOptim.var"#118#134"{GalacticOptim.var"#115#131"{OptimizationFunction{true,GalacticOptim.AutoForwardDiff{nothing},DiffEqFlux.var"#73#78"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing},GalacticOptim.var"#116#132"{Array{Float64,1},Int64,GalacticOptim.var"#115#131"{OptimizationFunction{true,GalacticOptim.AutoForwardDiff{nothing},DiffEqFlux.var"#73#78"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}}},GalacticOptim.var"#122#138"{GalacticOptim.var"#115#131"{OptimizationFunction{true,GalacticOptim.AutoForwardDiff{nothing},DiffEqFlux.var"#73#78"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing},GalacticOptim.var"#120#136"{Array{Float64,1},Int64,GalacticOptim.var"#115#131"{OptimizationFunction{true,GalacticOptim.AutoForwardDiff{nothing},DiffEqFlux.var"#73#78"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}}},GalacticOptim.var"#124#140",Nothing,Nothing,Nothing},Array{Float64,1},SciMLBase.NullParameters,Nothing,Nothing,Nothing,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::ADAM, ::Base.Iterators.Cycle{Tuple{GalacticOptim.NullData}}; maxiters::Int64, cb::Function, progress::Bool, save\_best::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at solve.jl:90  
\_\_solve at solve.jl:66 [inlined]  
\_\_solve at solve.jl:66 [inlined]  
#solve#468 at solve.jl:3 [inlined]  
(::CommonSolve.var"#solve##kw")(::NamedTuple{(:maxiters,),Tuple{Int64}}, ::typeof(solve), ::OptimizationProblem{true,OptimizationFunction{true,GalacticOptim.AutoForwardDiff,DiffEqFlux.var"#73#78"{typeof(loss)},GalacticOptim.var"#118#134"{GalacticOptim.var"#115#131"{OptimizationFunction{true,GalacticOptim.AutoForwardDiff{nothing},DiffEqFlux.var"#73#78"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing},GalacticOptim.var"#116#132"{Array{Float64,1},Int64,GalacticOptim.var"#115#131"{OptimizationFunction{true,GalacticOptim.AutoForwardDiff{nothing},DiffEqFlux.var"#73#  
(First code i shared produces this error for me)

The second problem is, that if i will get this code somehow running, how can i define an appropriate loss function for my parameter estimation problem if i have measurements of the different ODE componenets at different time points. At the moment my function data\_creator tries to flatten my ODE solution in order to compare it to flattened measurements. You mentioned a mutation error, how do i need to rewrite the loss function in order to avoid it?

I am sorry for asking these probably basic questions to you, but i am relatively new to all of this and currently learning a lot.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 5, 2022, 12:22pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/12 "2022-02-05T12:22:40Z")

</div>

What code are you running that is erroring?

---

<div class="post-metadata">

**Author:** ![Marius\_123](https://avatars.discourse-cdn.com/v4/letter/m/a5b964/32.png) [@Marius\_123](https://discourse.julialang.org/u/Marius_123)\
**Post date:** [February 25, 2022, 10:50am UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/13 "2022-02-25T10:50:58Z")

</div>

I managed to find the solution, it indeed worked as described by you initially! So i have an additional question regarding sciml\_train, just to be sure that i understand the function correctly: In the end this is an optimizer that computes the gradient wrt. to the parameters or initial conditions of the ODE with the sensitivity method of choice. This is used in an user specified gradient based optimizer to find the optimal parameters/initial conditions of a loss function which compares the ODE output to some data. Is this correct?  
Is there a way to not only return the minimizer of the loss function, but also the current gradient/sensitivity?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 28, 2022, 1:52pm UTC](https://discourse.julialang.org/t/creating-a-loss-function-for-sciml-train-with-unequal-time-point-distribution-of-the-ode-components/75204/14 "2022-02-28T13:52:27Z")

</div>

> [@Marius\_123](#):
>
> Is there a way to not only return the minimizer of the loss function, but also the current gradient/sensitivity?

Not right now. Easiest thing to do is just re-evaluate, but this is something that could be added.