About

edit

This page assists a talk at the Julia@BGI community event at 11 of May 2023 at the Max Planck Institute for Biogeochemistry.

Logistics

The talk's materials are accessible on-wiki via https://w.wiki/6geW. You can also explore the page in slide form

  • specifying and solving Ordinary Differential Equations (ODE) and variants
    • i.e. Box-models
    • part of Sciml organization
  • specify symbolically
    • MTK simplifies and generates fast code for DifferentialEqations.jl

function sesam3C(;name, k_N=60.0)

    @parameters t 

    D = Differential(t)

    @parameters ϵ ϵ_tvr κ_E a_E  m  τ  

    ...

    @variables (begin

        B(t),  L(t),   R(t),  cumresp(t),

        ...

    end)

    eqs = [

        D(B) ~ dB, dB ~ syn_B - tvr_B,

        D(L) ~ dL, dL ~ -dec_L + i_L,

        D(R) ~ dR, dR ~ -dec_R + ϵ_tvr*tvr_B + (1-κ_E)*syn_Enz,

        syn_Enz ~ a_E*B, tvr_Enz ~ syn_Enz,

        ...

        ]

    ODESystem(eqs; name)

end

function sesam3N(;name, sC = sesam3C(name=:sC))

   @unpack L, R, B, dec_L, dec_R, i_L, ϵ_tvr, tvr_B, tvr_B0, syn_B, syn_Enz, tvr_Enz, r_tvr, κ_E = sC

eqs = [

        β_NL ~ L/L_N, β_NR ~ R/R_N,

        D(L_N) ~ dL_N, dL_N ~ -dec_L/β_NL + i_L/β_Ni,

    ]

    extend(ODESystem(eqs, t, sts, ps; name), sC)

end

SIR:Implementierung in Octave:
edit
% Modifiziertes SIR Modell, Implementierung in Octave
%-------------------
%Parameter 
B=55000 %Kapazität: 2/3 der deutschen Bevölkerung (in Tausend)
%Raten
c=0.3238
w=1/15;
t=0.003;
%Perzentualanteil der erfassten Infizierten
r=0.1

times=(0:0.1:180);
%Anfangswerte
yo=[B-16/1000;16/1000;0];

% Funktion f der rechten Seite des DGL-Systems y'=f(y,t)  
 f=@(y,x) [-c*y(1)*y(2)/(r*B);c*y(1)*y(2)/B-w*y(2);w*y(2) -d*y(2)];
 y = lsode (f, yo, times);
 
 plot (times, y (:,1),'.',times, y (:,2), '.', times, y (:,3),'.',times, y (:,3)+y (:,2),'.')
 legend ("Succeptible","Infected", "Removed", "cummulative  Infected ","location", "east")
 xlabel ('Tage')
 ylabel ('Befoelkerung (in Tausend)')
 axis([0,180, 0 ,55600])

Experiences

edit
  • nice to work with observables
    • syn_Enz unsed for simplifying the system but computed only when quering the solution object
  • Iterative computations not possible any more within derivative function
  • Hard to debug
    • cannot look/step at generated derivative code
  • Position of state variables and parameters only known after simplifying the system
    • Hassle on optimizing a subset of parameters
  • still some problems with Automatic differentiation
    • cannot use all the solvers
    • fast response on issue, but hard to understand the aid
  • Bayesian analysis with flexible statistical model
  • used for inferring parameters and their uncertainty given observations
  • util_sample.jl
    • @model function sppfit(obs, ::Type{T} = Float64) where {T}
a
b
c