- Join over
**1.5M+ people** - Join over
**100K+ communities** - Free
**without limits** - Create
**your own community**

- 08:46robertfeldt opened #143
- 08:15robertfeldt closed #142
- 05:45amadeusine starred SciML/RecursiveArrayTools.jl
- 00:32ChrisRackauckas commented #345
- Apr 01 21:48Trel725 starred SciML/DifferentialEquations.jl
- Apr 01 20:47tkf commented #241
- Apr 01 19:25ViralBShah closed #238
- Apr 01 19:25ViralBShah commented #238
- Apr 01 19:24ViralBShah closed #227
- Apr 01 19:24ViralBShah commented #227
- Apr 01 17:55aarmey opened #233
- Apr 01 17:46FHell commented #583
- Apr 01 17:45FHell commented #583
- Apr 01 16:38ChrisRackauckas commented #241
- Apr 01 16:33
ChrisRackauckas on master

Fix reduction func in ensembles… Merge pull request #345 from jl… (compare)

- Apr 01 16:33ChrisRackauckas closed #345
- Apr 01 16:32ChrisRackauckas commented #583
- Apr 01 16:12omalleyc edited #584
- Apr 01 16:11omalleyc opened #584
- Apr 01 15:53
ChrisRackauckas on master

Small fixes to the DAE example … Merge pull request #346 from as… (compare)

[slack] <rveltz> :thumbsup:

[slack] <rveltz> I am always worried my question is dumb

[slack] <Thomas Propson> Hi all, I am new to both Julia and the DiffEq package, but I have experience with differential equations and autodiff. I am attempting to do an adjoint sensitivity analysis on a cost function on the final state of my system. Following the documentation [0] it looks like I want to use the Zygote package and the `concrete_solve`

interface. It looks like I could do something similar with the `adjoint_sensitivities`

interface, but I can’t express the derivative of my cost function easily and want to backprop rather than forwardprop or finite difference. I am not sure what I am doing wrong. This is not the problem that I actually want to solve but this is the minimum example I put together to show what I am trying to do.

```# foo.jl

using DiffEqSensitivity

using DifferentialEquations

using Zygote

module Foo

function rhs(u, p, t)

return u * p

end

function cost(u, p, prob)

sol = concrete_solve(prob, Tsit5(), u, p,

sensealg=QuadratureAdjoint())

final_u = sol.u[length(sol.u)]

return final_u * 2

end

function main()

u0 = 1.0

p = 1.0

end_time = 1.0

tspan = (0.0, end_time)

prob = ODEProblem(rhs, u0, tspan, p=p)

res = Zygote.gradient((u, p) -> cost(u, p, prob), u0, p)

println(res)

end

end`The full stack trace is:`

julia> include("foo.jl")

WARNING: replacing module Foo.

Main.Foo

julia> Main.Foo.main()

ERROR: BoundsError: attempt to access 0-element Array{Any,1} at index []

Stacktrace:

[1] throw_boundserror(::Array{Any,1}, ::Tuple{}) at ./abstractarray.jl:537

[2] checkbounds at ./abstractarray.jl:502 [inlined]

[3] _getindex at ./abstractarray.jl:1002 [inlined]

[4] getindex(::Array{Any,1}) at ./abstractarray.jl:980

[5] (::Zygote.var"#back#183"{:u,Zygote.Context,RecursiveArrayTools.DiffEqArray{Float64,1,Array{Float64,1},Array{Float64,1}},Array{Float64,1}})(::Array{Float64,1}) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/lib/lib.jl:199

[6] #355#back at /home/tcpropson/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49 [inlined]

[7] cost at /home/tcpropson/repos/qoc/_test/foo.jl:89 [inlined]

[8] (::typeof(∂(cost)))(::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0

[9] #1 at /home/tcpropson/repos/qoc/_test/foo.jl:99 [inlined]

[10] (::typeof(∂(λ)))(::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0

[11] (::Zygote.var"#36#37"{typeof(∂(λ))})(::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:36

[12] gradient(::Function, ::Float64, ::Vararg{Float64,N} where N) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:45

[13] main() at /home/tcpropson/repos/qoc/_test/foo.jl:99

[14] top-level scope at REPL[94]:1```

[0] https://docs.sciml.ai/dev/analysis/sensitivity/#concrete_solve-Examples-1

[slack] <chrisrackauckas> See if `

fixes it

```
Array(concrete_solve(prob, Tsit5(), u, p,
sensealg=QuadratureAdjoint()))
```

fixes it

[slack] <Thomas Propson> No, unfortunately. Thank you for your quick response. The new error is:

```
julia> Main.Foo.main()
ERROR: MethodError: no method matching vec(::Float64)
Closest candidates are:
vec(::LinearAlgebra.Transpose{T,#s662} where #s662<:(AbstractArray{T,1} where T) where T) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/adjtrans.jl:201
vec(::SparseArrays.AbstractSparseArray{Tv,Ti,1} where Ti where Tv) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsevector.jl:912
vec(::StaticArrays.StaticArray) at /home/tcpropson/.julia/packages/StaticArrays/1g9bq/src/abstractarray.jl:163
...
Stacktrace:
[1] _broadcast_getindex_evalf at ./broadcast.jl:631 [inlined]
[2] _broadcast_getindex at ./broadcast.jl:604 [inlined]
[3] getindex at ./broadcast.jl:564 [inlined]
[4] copy at ./broadcast.jl:854 [inlined]
[5] materialize at ./broadcast.jl:820 [inlined]
[6] Array(::RecursiveArrayTools.DiffEqArray{Float64,1,Array{Float64,1},Array{Float64,1}}) at /home/tcpropson/.julia/packages/RecursiveArrayTools/euiv5/src/vector_of_array.jl:13
[7] adjoint at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/lib/array.jl:7 [inlined]
[8] _pullback(::Zygote.Context, ::Type{Array}, ::RecursiveArrayTools.DiffEqArray{Float64,1,Array{Float64,1},Array{Float64,1}}) at /home/tcpropson/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:47
[9] cost at /home/tcpropson/repos/qoc/_test/foo.jl:87 [inlined]
[10] _pullback(::Zygote.Context, ::typeof(Main.Foo.cost), ::Float64, ::Float64, ::DiffEqBase.ODEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,DiffEqBase.ODEFunction{false,typeof(Main.Foo.rhs),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:p,),Tuple{Float64}}},DiffEqBase.StandardODEProblem}) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0
[11] #1 at /home/tcpropson/repos/qoc/_test/foo.jl:99 [inlined]
[12] _pullback(::Zygote.Context, ::Main.Foo.var"#1#2"{DiffEqBase.ODEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,DiffEqBase.ODEFunction{false,typeof(Main.Foo.rhs),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:p,),Tuple{Float64}}},DiffEqBase.StandardODEProblem}}, ::Float64, ::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0
[13] _pullback(::Function, ::Float64, ::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:29
[14] pullback(::Function, ::Float64, ::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:35
[15] gradient(::Function, ::Float64, ::Vararg{Float64,N} where N) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:44
[16] main() at /home/tcpropson/repos/qoc/_test/foo.jl:99
[17] top-level scope at REPL[96]:1
```

[slack] <chrisrackauckas> oh I see what it is

[slack] <chrisrackauckas> the issue is that the adjoint methods don't work with non-mutable values right now, so you need to use an array instead of a scalar.

[slack] <chrisrackauckas> yes

[slack] <Thomas Propson> I think my problem is that I can’t grab the final state the way I’m currently trying. For instance, this works:

But this doesn’t:

My question is: how do I get the gradients with respect to my cost function evaluated only at the final state.

```
function cost(u0, p, prob)
sol = concrete_solve(prob, Tsit5(), u0, p, saveat=0.1, sensealg=QuadratureAdjoint())
return sum(sol)
end
```

But this doesn’t:

```
function cost(u0, p, prob)
sol = concrete_solve(prob, Tsit5(), u0, p, saveat=0.1, sensealg=QuadratureAdjoint())
final_u = sol.u[length(sol.u)]
return final_u[1]
end
```

My question is: how do I get the gradients with respect to my cost function evaluated only at the final state.

function τ(t::Float64, p::parameters)::Float64

return t<p.t_drug ? 0 : p.ε_D* p.τ_0*

end

function ψ(t::Float64,p::parameters)::Float64

return t < p.t_vac ? 0 : p.ε_V p.ψ_0

end

return t<p.t_drug ? 0 : p.ε_D

end

function ψ(t::Float64,p::parameters)::Float64

return t < p.t_vac ? 0 : p.ε_V

end

results of @btime:

1.960 μs (28 allocations: 1.56 KiB)

1.960 μs (28 allocations: 1.56 KiB)

I know it's not allocating much but I feel like this model should be much faster.

@ChrisRackauckas your diffeq library is awesome and it makes my research with SDEs muuuuch easier, thank you!

use that macro rather than the broadcasting?

[slack] <chrisrackauckas> for long expressions

[slack] <chrisrackauckas> I don't know if that's long enough, but worth a try

[slack] <chrisrackauckas> or maybe you might need to @views

[slack] <chrisrackauckas> I think it may have too tight of a type on it

[slack] <chrisrackauckas> ``` function rober(du,u,p,t)

y₁,y₂,y₃ = u

k₁,k₂,k₃ = p

du[1] = -k₁*y₁+k₃*y₂*y₃*

du[2] = k₁y₁-k₂*y₂^2-k₃*y₂*y₃

du[3] = y₁ + y₂ + y₃ - 1

nothing

end

M = [1. 0 0

0 1. 0

0 0 0]

f = ODEFunction(rober,mass_matrix=M)

p = [0.04,3e7,1e4]

y₁,y₂,y₃ = u

k₁,k₂,k₃ = p

du[1] = -k₁

du[2] = k₁

du[3] = y₁ + y₂ + y₃ - 1

nothing

end

M = [1. 0 0

0 1. 0

0 0 0]

f = ODEFunction(rober,mass_matrix=M)

p = [0.04,3e7,1e4]

```
prob_singular_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,100),p)
sol_singular_mm = solve(prob_singular_mm,Rodas5(),reltol=1e-8,abstol=1e-8)
ts = [50, sol_singular_mm.t[end]]
dg_singular(out,u,p,t,i) = (fill!(out, 0); out[end] = -1)
_, res = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=QuadratureAdjoint())
reference_sol = ForwardDiff.gradient(p->G(p, prob_singular_mm, ts, sol->sum(last, sol.u)), vec(p))
@test res' ≈ reference_sol rtol = 1e-2
_, res_interp = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint())
@test res_interp ≈ res rtol = 1e-2
_, res_interp2 = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint(checkpointing=true),checkpoints=sol_singular_mm.t[1:10:end])
@test res_interp2 ≈ res rtol = 1e-2```
```

[slack] <chrisrackauckas> The other forward mode choice, open an issue on it

[slack] <chrisrackauckas> it's easy to fix but I haven't done it yet.

@ChrisRackauckas I get another error using ODELocalSensitivityProblem with ForwardDiffSensitivity that appears to be caused by the mass matrix: "ERROR: DimensionMismatch("W: (Base.OneTo(12), Base.OneTo(12)), mass matrix: (Base.OneTo(6), Base.OneTo(6))")". I will raise this as an issue. Thanks for your help