# StackOverflowError trying to fit a difference equation with Turing

**URL:** <https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349>\
**Category:** Statistics\
**Tags:** diffeq, turing\
**Created:** [August 5, 2020, 3:54pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349 "2020-08-05T15:54:13Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 5, 2020, 3:54pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/1 "2020-08-05T15:54:13Z")

</div>

I am trying to sample a difference equation using Turing but encounter errors.

Basically, it’s a SIR model that solves fine, with a quick solver function written below (I wrote my own very simple solver to understand the basics and see how the solver could be altered/extended for modeling purposes).

Then, I add random noise to the solution and use it as data for the Turing model, in order to find the parameters of the SIR model using Turing. I used this strategy successfully with the Lokta Voltera model or the tutorial (see [my modification of the tutorial here](https://github.com/arnaudmgh/TuringTutorials/blob/master/10_diffeq.ipynb)), but I can’t make it work with my own simple solver.

Any help is appreciated.

## Julia version: 1.5.0.

## Error

```julia-auto
chain = sample(model, NUTS(.65),1000)

ERROR: StackOverflowError:
Stacktrace:
 [1] ScalMat at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:5 [inlined] (repeats 2 times)
 [2] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:13 [inlined]
 [3] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:14 [inlined]
 [4] MvNormal(::Array{Real,1}, ::PDMats.ScalMat{Float64}) at /Users/aamzallag/.julia/packages/Distributions/jFoHB/src/multivariate/mvnormal.jl:202 (repeats 79984 times)

```

## Code

```julia-auto
using Turing, Distributions

# Import MCMCChain, Plots, and StatsPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots

# Set a seed for reproducibility.
using Random
Random.seed!(12);

# Defining the dynamics
function sir_model(du,u,p,t)
    S, I, R = u
    β, γ = p
    du[1] = dS = -β*S*I/N
    du[2] = dI = β*S*I/N - γ*I
    du[3] = dR = γ*I
end

### trying a simple discrete time solver 

# this solver takes a function func_ describing
# the dynamics (a la du/dt = f(u,p,t)) and simulate
# the dynamics with a simple loop, and store the values
# in variable u
function diffsolve(func_, u0, p, n=100)
    du = [0., 0., 0.]
    u_t = copy(u0)
    #u = [u0]
    u = Array{Array{Real,1},1}() # if Real not specified, 
    # I get the error: in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,3}
    push!(u, u0)
    for t in 1:n
        func_(du, u_t, p, 1)
        u_t = u_t + du
        push!(u, u_t)
    end
    #hcat(u...)
    u
end

N=10e3 # total population
conf=10 # confirmed cases at the beginning
u0 = [N-conf,conf,0.0] #.* N

## Trying the solver and plotting
x=diffsolve(sir_model, u0, [0.8, 0.1], 50)
x=hcat(x...)
# add noise
#odedata = x + 0.1*maximum(x) .* randn(size(x))
odedata = x .* exp.(0.1 .* randn(size(x)))

plot(x')
scatter!(odedata')

Turing.setadbackend(:forwarddiff)

@model function fitsir(data)
    σ ~ InverseGamma(1, 1)
    #σ ~ truncated(Normal(0., 0.5),0,5)
    β ~ truncated(Normal(0., 0.5),0,5)
    γ ~ truncated(Normal(0., 0.5),0,5)
    
    p = [β,γ]
    predicted = diffsolve(sir_model, u0, p, 10)

    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end

model = fitsir(odedata)
chain = sample(model, NUTS(.65),1000)

# ERROR: StackOverflowError:
# Stacktrace:
# [1] ScalMat at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:5 [inlined] (repeats 2 times)
# [2] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:13 [inlined]
# [3] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:14 [inlined]
# [4] MvNormal(::Array{Real,1}, ::PDMats.ScalMat{Float64}) at /Users/aamzallag/.julia/packages/Distributions/jFoHB/src/multivariate/mvnormal.jl:202 (repeats 79984 times)

```

---

<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:** [August 5, 2020, 4:18pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/2 "2020-08-05T16:18:56Z")

</div>

> [@arnaudmgh](#):
>
> `MvNormal(::Array{Real,1}, ::PDMats.ScalMat{Float64})`

The way that you write the array made it not concretely typed and then it hits a bad dispatch.

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 5, 2020, 4:48pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/3 "2020-08-05T16:48:38Z")

</div>

Thanks Chris,  
when I try to declare u as

```julia
u = []

```

or

```julia
u = Array{Array{Float64,1},1}()

```

I get  
`ERROR: TypeError: in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,3}`  
I don’t know else how to collect the results from the loop

---

<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:** [August 6, 2020, 6:38am UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/4 "2020-08-06T06:38:37Z")

</div>

> [@arnaudmgh](#):
>
> ```julia
> #u = [u0]
> u = Array{Array{Real,1},1}() # if Real not specified, 
> 
> ```

You want to make this be dual valued, and only your parameters are dual valued, so the easiest thing to do is:

```julia
predicted = diffsolve(sir_model, eltype(p).(u0), p, 10)

```

and

```julia
u = [u0]

```

on the inside. That should propagate the duals in a type-stable way. I’ll just quickly note that DifferentialEquations.jl’s AD integration takes care of these cases for you.

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 6, 2020, 9:24pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/5 "2020-08-06T21:24:12Z")

</div>

> [@ChrisRackauckas](#):
>
> `eltype(p).(u0)`

Thanks Chris,  
unfortunately this makes no difference.

Good to know that DifferentialEquations.jl handles it for the user. It would be nice if Turing could be handling these kinds of problems, or return an informative error, for problems that are not differential equation problems.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [August 7, 2020, 12:43pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/6 "2020-08-07T12:43:48Z")

</div>

Turing handles these issues using a specific syntax. You can find details on turing.ml in the documentation on AD.

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 8, 2020, 2:03am UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/7 "2020-08-08T02:03:34Z")

</div>

Thank you @trappmartin - I am not sure there is something relevant to my problem in the AD section - isn’t it rather in the performance tips, the section on type stability? [Performance Tips](https://turing.ml/dev/docs/using-turing/performancetips)

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 21, 2020, 11:21pm UTC](https://discourse.julialang.org/t/stackoverflowerror-trying-to-fit-a-difference-equation-with-turing/44349/8 "2020-08-21T23:21:51Z")

</div>

This was answered by @wupeifan on Julia Slack. Replacing the definition of `diffsolve` by what follows resolves the problem:

```julia
function diffsolve(func_, u0, p, n=100)
    u_t = deepcopy(u0)
    T = promote_type(eltype(u0), eltype(p))
    u = Array{Array{T,1},1}() 
    du = Vector{T}(undef, 3) 

    push!(u, u0)
    for t in 1:n
        func_(du, u_t, p, 1)
        u_t = u_t + du
        push!(u, u_t)
    end
    u
end

```

This solution converts both `du` _and_ `u0` to the same type as `p` (not only `u0` as in a previous answer) and hence, it works.