New noise: IntrinsicTimescales.jl

[behinger]

We could add a new noisetype based on IntrinsicTimescales.jl

struct IntrinsicNoise <: AbstractNoise
    noiselevel
	sfreq
	tau
end
function UnfoldSim.simulate_noise(rng, t::IntrinsicNoise, n::Int)
	noise = generate_ou_process(t.tau, 1, 1/t.sfreq, n/t.sfreq, 1)[1,:]
    return t.noiselevel .* noise

end

This is similar to ExponentialNoise, but I think a better parameterization. Side-effect is, it requires many dependencies, so we could add it as an extension; Problem then is only, that the type has to be defined within UnfoldSim proper.

Alternatively, we use it as a demo to include other noisetypes, with the sideeffect of not being very discoverable.

Let's first include #129 though, the new exponential implementation.

---

This issue depends on duodenum96/IntrinsicTimescales.jl#25 being resolved

[duodenum96]

Hey!

I'm honored that my little package is useful to somebody. But this can be much more elegantly handled without adding a lot of dependencies coming from IntrinsicTimescales.jl, as you mention above. Essentially an OU process is just a linear stochastic differential equation with constant coefficients and has a definite timescale (as opposed to 1/f type of noise with no timescale, i.e. scale-free). So it is a colored noise but not pink or brown noise. My generate_ou_process is just a wrapper around StochasticDiffEq.jl. Anyway, one can take the relevant parts of the code and make a much more concise function:

import StochasticDiffEq as deq

f = (du, u, p, t) -> du .= -u ./ p[1] # p[1]: timescale
g = (du, u, p, t) -> du .= p[2] # p[2]: noise variance

n_trials = 10
dt = 0.01
tspan = (0.0, 10.0)
u0 = randn(n_trials)
saveat = dt:dt:tspan[2]

p = [1.0, 1.0]

prob = deq.SDEProblem(f, g, u0, tspan, p)
sol = deq.solve(prob, deq.SOSRA(), saveat=saveat)

sol_array = Array(sol) # n_trials x n_timesteps

One can make a function out of this as:

f = (du, u, p, t) -> du .= -u ./ p[1] # p[1]: timescale
g = (du, u, p, t) -> du .= p[2] # p[2]: noise variance
function generate_ou_process(timescale, noise_variance, n_trials, dt, t_end)
    p = [timescale, noise_variance]
    tspan = (0.0, t_end)
    u0 = randn(n_trials)
    saveat = dt:dt:t_end
    prob = deq.SDEProblem(f, g, u0, tspan, p)
    sol = deq.solve(prob, deq.SOSRA(), saveat=saveat)
    return Array(sol)
end

Note about deq.SOSRA(): I picked this as the best performing one (at least on my computer) after trying out various solvers in https://docs.sciml.ai/DiffEqDocs/dev/solvers/sde_solve/. In principle one can also do sol = deq.solve(prob) which picks an algorithm automatically but this requires the full DifferentialEquations.jl package which is quite a heavy dependency.

The rest of my code in generate_ou_process just hides the details here and makes certain performance improvements by using static arrays where appropriate, remaking the problem (remake) to avoid creating the problem every time the simulation is run etc. Normally StochasticDiffEq.jl is fast, as in super super-duper fast (IIRC the benchmarks for it was your example in Unfold.jl workshop, and in fact the parent package DifferentialEquations.jl was the reason I got into Julia back in the day, it had the fastest differential equation solvers especially for more complicated models). Most of the time one doesn't need these performance improvements. I put them there because I was solving this thousands of times in Bayesian INT calculations.

Anyway, let me know if I can be of help in this.

[behinger]

Hey!
Thanks a ton for the feedback. I agree that the most parsimonious solution would be to just include the code - but then, I much more prefer someone else maintain the code I use (with the caveat that breaking changes are partially outside of my doing). Also, credit where credit is due :)

The dependencies are not so bad, as I propose to add it as an extension.

Let' see, we first have to merge 20 other PRs 😅