Allow parameters (and initial conditions) to be distributions

[TorkelE]

It would be really cool, and useful, to have a parameter be a distribution:

using ModelingToolkit, Distributions
using Plots, OrdinaryDiffEq

@variables t x(t)
@parameters p d
D = Differential(t) 

@named bd = ODESystem([D(x) ~ p - d * x])

prob = ODEProblem(bd, [x => 2.0], (0.0, 20.0), [p => Uniform(0.8,1.2), d => 0.2])
sol = solve(prob, Tsit5())
plot(sol)

then, each time one run sol = solve(prob, Tsit5()), a p is drawn from Uniform(0.8,1.2) and used for the simulation. In systems biology this could e.g. be used to model extrinsic noise, and I'm sure there are tons of applications for this all over the place.

I've opened a related issue in Catalyst on the same topic: https://github.com/SciML/Catalyst.jl/issues/697

[baggepinnen]

using ModelingToolkit, Distributions
using Plots, OrdinaryDiffEq
using MonteCarloMeasurements

@variables t x(t)
@parameters p d
D = Differential(t) 

@named bd = ODESystem([D(x) ~ p - d * x])

prob = ODEProblem(bd, [x => 2.0], (0.0, 20.0), [p => Particles(Uniform(0.8,1.2)), d => 0.2])
sol = solve(prob, Tsit5())
plot(sol)

This approach has some limitations though, it does not play well with events for one.

[baggepinnen]

I do have some experimental code that would translate the uncertain model representation into an EnsembleProblem which could be located in DiffEqBaseMonteCarloMeasurementsExt.jl. I recal the code below working in at least some situations, but I couldn't quite figure out how to properly construct the ODESolution in postprocess in the general case

const ParticleODEProblem = ODEProblem{<:Any, <:Any, <:Any, <:MonteCarloMeasurements.SomeKindOfParticles}

function SciMLBase.EnsembleProblem(prob::ParticleODEProblem, args...; kwargs...)
    p = copy(prob.p)
    u0 = copy(prob.u0)
    N = max(nparticles(eltype(p)), nparticles(eltype(u0)))
    p0 = MonteCarloMeasurements.vecindex.(p, 1)
    u00 = MonteCarloMeasurements.vecindex.(u0, 1)
    prob0 = remake(prob, p = p0, u0 = u00)

    prob_func = function(probi,i,repeat)
        for j in eachindex(probi.u0)
            probi.u0[j] = MonteCarloMeasurements.vecindex(u0[j], i)
        end
        for j in eachindex(probi.p)
            probi.p[j] = MonteCarloMeasurements.vecindex(p[j], i)
        end
        probi
    end
    # reduction = function(u,data,I)
    #     (append!(u,data),false)
    # end
    
    eprob = EnsembleProblem(prob0; prob_func)
end

function DiffEqBase.solve(prob::ParticleODEProblem, alg, args...; kwargs...)
    N = max(nparticles(eltype(prob.p)), nparticles(eltype(prob.u0)))
    eprob = EnsembleProblem(prob)
    esol = DiffEqBase.solve(eprob, alg, EnsembleThreads(); trajectories = N, kwargs...)
    postprocess(esol)
end

function postprocess(esol)
    soli = esol[1]
    t = soli.t
    nx = length(soli[1])

    # [state_index][mc_index]
    utraj = map(t) do t
        data = reduce(hcat, OrdinaryDiffEq.EnsembleAnalysis.componentwise_vectors_timepoint(esol,t)) # nmc × nx
        xparts = Particles(data)
    end

    u_analytic  = nothing
    errors      = nothing
    k           = nothing
    prob        = soli.prob
    alg         = soli.alg
    interp      = nothing
    dense       = false
    tslocation  = soli.tslocation
    destats     = nothing
    alg_choice  = soli.alg_choice
    retcode     = soli.retcode
    ODESolution{eltype(soli.u[1]), 2}(utraj, u_analytic, errors, t, k, prob, alg, interp, dense,
    tslocation, destats, alg_choice, retcode)
end

[TorkelE]

This is both really cool, I didn't know things were so far progressed!

I think both your first example, and something like in the second (which I understand, if trajectories=1 will correspond to something that I suggested) would be really useful. Is there something inherently preventing the simple notation, and then users can create their own EnsembleProblems?

In my application, I would actually like to create an SDEProblem. These models the intrinsic noise of biological processes, and having parameters being distributions would model extrinsic noise. If I created such an SDEProblem, by modulating the two types of noise I could investigate how they affect biological systems in a way I don't think have been done.

Next, if we can create a nice interface for it, I would like to document and advertise it, since I don't think there are any good tools fo extrinsic noise modulation (although it is something you are able to set up yourself).