# Differentiating through a Jump Problem

**URL:** <https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496>\
**Category:** Numerics\
**Created:** [November 25, 2019, 2:32pm UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496 "2019-11-25T14:32:20Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [November 26, 2019, 8:02pm UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/21 "2019-11-26T20:02:01Z")

</div>

Lots of interesting stuff has been done with regards to sensitivity analysis, I think also involving coupling through multilevel methods in recent years. David Anderson has a lot on this too.

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [November 27, 2019, 11:18am UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/22 "2019-11-27T11:18:17Z")

</div>

Sorry, I posted sloppy code. I’ve cleaned it up now

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [November 27, 2019, 11:21am UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/24 "2019-11-27T11:21:44Z")

</div>

EDIT: I’ve made a less sloppy version of the above code that’s a bit more minimal and actually runs!

> using DifferentialEquations  
> using DiffEqBiological  
> using ForwardDiff  
> using LinearAlgebra
> 
> sir\_model = @reaction\_network SIR begin  
> c1, s + i → 2i  
> c2, i → r  
> end c1 c2
> 
> p = [0.1/1000,0.01]  
> prob = DiscreteProblem([999,1,0],(0.0,250.0),p)
> 
> function forward\_pass(p\_)  
> prob = DiscreteProblem([999,1,0],(0.0,250.0),p\_)  
> prob = remake(prob; u0=convert.(eltype(p\_),prob.u0),p=p\_)  
> jump\_prob = JumpProblem(prob, Direct(), sir\_model)  
> sol = solve(jump\_prob,FunctionMap())  
> loss = norm(sol[end] .- [1,200,700])  
> return loss  
> end
> 
> loss\_gradient = p\_ → ForwardDiff.gradient(forward\_pass,p\_)  
> println(loss\_gradient(p))

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [November 27, 2019, 12:06pm UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/25 "2019-11-27T12:06:43Z")

</div>

It didn’t! Though that’s probably my coding inexperience. Code and error posted below.

However, maybe this whole idea is ill-conceived. Even with tau leaping, we have an everywhere non-differentiable (in the mathematical sense) function from parameters to to trajectory (which is a step function, being a Discrete problem). I’m not seeing how ForwardDiff can propagate a dual through these discrete jumps.

Code:

> using DifferentialEquations  
> using DiffEqBiological  
> using ForwardDiff  
> using LinearAlgebra  
> using Plots  
> pyplot()
> 
> “”" rates of each reaction “”"  
> function rate(out,u,p,t)  
> out[1] = p[1]\*u[1]\*u[2] #infection  
> out[2] = p[2]\*u[2] # recovery  
> end
> 
> u0 = [999.0,1.0,0.0] #susceptible infected recovered  
> p0 = [0.1/1000,0.01]
> 
> “”" change matrix upon reaction “”"  
> function c(dc,u,p,t,mark)  
> “”" first reaction loses susceptible, gains infected “”"  
> dc[1,1] = -1  
> dc[2,1] = 1  
> “”" second reaction loses infected, gains recovered “”"  
> dc[2,2] = -1  
> dc[3,2] = 1  
> end
> 
> “”" naive problem solution “”"  
> dc = zeros(3,2)  
> rj = RegularJump(rate,c,dc;constant\_c=true)

> prob = DiscreteProblem(u0,(0.0,250.0),p0)  
> jump\_prob = JumpProblem(prob,Direct(),rj)  
> sol = solve(jump\_prob,SimpleTauLeaping();dt=1.0)
> 
> function forward\_pass(p\_)  
> prob = DiscreteProblem(u0,(0.0,250.0),p\_)  
> prob = remake(prob; u0=convert.(eltype(p\_),prob.u0),p=p\_)  
> jump\_prob = JumpProblem(prob,Direct(),rj)  
> sol = solve(jump\_prob,SimpleTauLeaping();dt=1.0)  
> return loss(sol)  
> end
> 
> function loss(sol)  
> return norm(sol[end] .- [1,200,700])  
> end
> 
> loss\_gradient = p\_ → ForwardDiff.gradient(forward\_pass,p\_)

Error:

> ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(forward\_pass),Float64},Float64,2})  
> Closest candidates are:  
> Float64(::Real, ::RoundingMode) where T\<:AbstractFloat at rounding.jl:194  
> Float64(::T\<:Number) where T\<:Number at boot.jl:741  
> Float64(::Int8) at float.jl:60  
> …

---

<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:** [November 27, 2019, 4:51pm UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/26 "2019-11-27T16:51:58Z")

</div>

It should be as differentiable as the SDE, in that there is a strong derivative. I’ll look at this a bit later.

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [November 27, 2019, 5:52pm UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/27 "2019-11-27T17:52:08Z")

</div>

Sure, thanks again!

I think there’s a difference: The output trajectory for a reactant is a step function that can only take integer values. You can’t have infinitesimal changes in the output trajectory over a fixed timestep. For an SDE, the output can be changed infinitesimally, as it takes real values.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 28, 2019, 12:31am UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/28 "2019-11-28T00:31:06Z")

</div>

I’m having a stab at this and I can’t quite get the DifEq part to accept my types. I manage to get the `prob` to be

```julia
julia> prob
DiscreteProblem with uType Array{Particles{Float64,500},1} and tType Float64. In-place: true
timespan: (0.0, 250.0)
u0: Particles{Float64,500}[999.0, 1.0, 0.0]

```

but the function `rate` is still called with the signature

```julia
(Array{Float64,1}, Array{Particles{Float64,500},1}, Array{Particles{Float64,500},1}, Float64)

```

i.e., the array to store the derivative in is still `Array{Float64,1}`.

Edit: PR [https://github.com/JuliaDiffEq/DiffEqJump.jl/pull/88](https://github.com/JuliaDiffEq/DiffEqJump.jl/pull/88)

I tried making the problem `inplace=false` but I couldn’t figure out how to do that.

---

<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:** [November 28, 2019, 12:41am UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/29 "2019-11-28T00:41:19Z")

</div>

> [@baggepinnen](#):
>
> i.e., the array to store the derivative in is still `Array{Float64,1}`

Interesting. I’ll take a look at that. Got a lot going on this week, but this is always fun procrastination work.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 28, 2019, 12:41am UTC](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496/30 "2019-11-28T00:41:57Z")

</div>

PR [https://github.com/JuliaDiffEq/DiffEqJump.jl/pull/88](https://github.com/JuliaDiffEq/DiffEqJump.jl/pull/88)

[Previous page](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496.md?page=1)