DifferentialEquations.jl and Measurements.jl

[giordano]

Today I was asked whether it was possible to solve in Julia differential equations involving numbers with uncertainties. Of course the answer is yes. What I find really amazing about Julia is that the two packages don’t know anything about each other, yet they can work together without any effort. Here is an short example based on this tutorial: Jupyter Notebook Viewer

The only think I found weird was that the extrema of the time interval had to be Measurement objects, otherwise solve would complain (because it can’t convert Measurement type to Float64). Why it’s that? I can see why u0 should have the same type as that returned by f, but not why also the time should match the same type.

[ChrisRackauckas]

@giordano

The only think I found weird was that the extrema of the time interval had to be Measurement objects, otherwise solve would complain (because it can’t convert Measurement type to Float64). Why it’s that?

Adaptive timestepping. When using adaptive timestepping (the default), the times that the solver is internally doing have an uncertainty as well since its steps are determined by the ODE system itself, so it needs to keep track of the uncertainty in time.

Cool stuff, I like it when Julia makes new features for us . Would you mind doing a PR to GitHub - SciML/SciMLTutorials.jl: Tutorials for doing scientific machine learning (SciML) and high-performance differential equation solving with open source software. ?

[giordano]

@ChrisRackauckas

Adaptive timestepping. When using adaptive timestepping (the default), the times that the solver is internally doing have an uncertainty as well since its steps are determined by the ODE system itself, so it needs to keep track of the uncertainty in time.

Thank you for the explanation!

@ChrisRackauckas

Cool stuff, I like it when Julia makes new features for us . Would you mind doing a PR to GitHub - SciML/SciMLTutorials.jl: Tutorials for doing scientific machine learning (SciML) and high-performance differential equation solving with open source software. ?

Yes, I was thinking about adding a tutorial, but I’d like to come up with a different example rather than just copying another tutorial

[TS-CUBED]

Stumbled across this after doing an ensemble based uncertainty study.

This could save a lot of time!

Unfortunately I couldn’t get it to work on a ModelingToolkit.jl ODEProblem (acausal model - structural_simplfy).

Am I missing something and this should work, or is there a problem with the Measuremnt-Type going through the more complex model?

[Paulo_Jabardo]

Did you reproduce the example above at least? I’m asking because the the interface of DiiferentialEquations.jl has changed (inverted some of the parameters).That might be it but I haven’t done that with MTK.

If that is the case, checkout the docs for the correct syntax and don’t forget that initial values and tspan should include uncertainty.

[TS-CUBED]

Actually, I hadn’t tried the example. I get the same error, so assume that this does not work anymore.

Pity. It would be neat and impressive to be able to carry the uncertainty and propagate it through. Although I cannot imagine how it is going to work on a non-linear equation …

[giordano]

@TS-CUBED

I get the same error

You didn’t clearly show what you have tried, and exactly what error you got.

Also, have you seen this tutorial (from the DiffEqTutorials mentioned above)?

[Paulo_Jabardo]

The example works if you correct for the new syntax, basically instead of f(t,u), f(u,p,t) where p stores the parameter α

julia> f(u,p,t) = p[1]*u
f (generic function with 1 method)

julia> u0 = 0.5 ± 0.0
0.5 ± 0.0

julia> tspan = (0.0±0.0, 1.0±0.0)
(0.0 ± 0.0, 1.0 ± 0.0)

julia> prob = ODEProblem(f, u0, tspan, [1.01±0.01])
ODEProblem with uType Measurement{Float64} and tType Measurement{Float64}. In-place: false
timespan: (0.0 ± 0.0, 1.0 ± 0.0)
u0: 0.5 ± 0.0

julia> sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 17-element Vector{Measurement{Float64}}:
                  0.0 ± 0.0
 0.012407826196308189 ± 0.0
 0.042501278333560696 ± 0.0
   0.0817804940926822 ± 0.0
  0.12887384498704435 ± 0.0
  0.18409796286152927 ± 0.0
  0.24627457447456758 ± 0.0
  0.31479297816557983 ± 0.0
   0.3885963515160237 ± 0.0
   0.4668617724420117 ± 0.0
   0.5487161305960653 ± 0.0
   0.6334346972152323 ± 0.0
   0.7203630000154827 ± 0.0
    0.808957991167541 ± 0.0
   0.8987655040395068 ± 0.0
   0.9894161889652783 ± 0.0
                  1.0 ± 0.0
u: 17-element Vector{Measurement{Float64}}:
      0.5 ± 0.0
 0.506305 ± 6.3e-5
  0.52193 ± 0.00022
  0.54305 ± 0.00044
  0.56951 ± 0.00073
   0.6022 ± 0.0011
   0.6412 ± 0.0016
   0.6871 ± 0.0022
   0.7403 ± 0.0029
   0.8012 ± 0.0037
   0.8703 ± 0.0048
    0.948 ± 0.006
    1.035 ± 0.0075
   1.1319 ± 0.0092
    1.239 ± 0.011
    1.358 ± 0.013
    1.373 ± 0.014

I couldn’t get it to work with ModelingToolkit but it does work using the lower level functions.

[TS-CUBED]

I hadn’t seen the tutorial you reference here, but the one in the original post works if the function definition is adjusted to reflect the changes in DifferentialEquations.jl - as in the DiffEq tutorials you mentioned.

So it does work with directly defined ODEs.

My original question was, if it is supposed or expected to work with MTK models that have a symbolic component. If not, then we can close the case. If yes, then we can start digging into it.

[ChrisRackauckas]

It would work if you create an ODEProblem from ModelingToolkit with a u0 that is a Measurement. I am not sure about the other potential use case of tracing through Measurements generated code.

[TS-CUBED]

I have recreated my model from scratch (since all the default values were Float64, I needed to replace all of them with defaults of type Measurement{Float64}) and now the error messages are a bit clearer and fit for discussion

When calling structural_simplify(rc_model) I get:

“”"
sys

TypeError: in Measurement, in T, expected T<:AbstractFloat, got Type{Real}

1. promote_rule @ conversions.jl:46 [inlined]
2. promote_type @ promotion.jl:233 [inlined]
3. promote_symtype (::typeof(+), ::Type{Measurements.Measurement{Float64}}, ::Type{Real})@ methods.jl:93
4. add_t @ types.jl:685 [inlined]
5. + (::Measurements.Measurement{Float64}, ::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})@ types.jl:700
6. - (::Measurements.Measurement{Float64}, ::SymbolicUtils.Term{Real, Base.ImmutableDict{DataType, Any}})@ types.jl:724
7. initialize_system_structure (::ModelingToolkit.ODESystem)@ systemstructure.jl:114
8. alias_elimination (::ModelingToolkit.ODESystem)@ alias_elimination.jl:6
9. var"#structural_simplify#70" (::Bool, ::typeof(ModelingToolkit.structural_simplify), ::ModelingToolkit.ODESystem)@ abstractsystem.jl:806
10. structural_simplify @ abstractsystem.jl:806 [inlined]
11. top-level scope @ Local: 1 [inlined]
    “”"

The Julia code for a simple model is here. I’ll need to triple check that I haven’t overlooked a Float64 initialisation in there, but couldn’t find any.

CR_measurements.jl

[TS-CUBED]

I was trying to have all the parameters having uncertainty, rather than u0.

We want to know

1. how important the uncertainties are on parameter estimates that come from clinical measurements
2. how close we need to get with our parameter estimate in order to get a pressure response that is within physiological parameters (i.e. not dangerous to a patient)

At the moment I define the uncertainties on the parameter set and run an ensemble analysis, i.e. a few thousand to ten thousand runs on the whole model** and had hoped to use this as an alternative approach.

** this is now, due to a 200x speed-up over the original MATLAB code, not as much of a problem as it used to be and can be done in about 10-20 seconds, but still, getting that down further could allow real-time predictions