Track variables within DifferentialEquations.jl function

New to Julia
1

Hi,

I have a simple function to solve:

using DifferejntialEquations, Plots, BenchmarkTools

var = []
function coupledd!(du, u, p, t)
    α, β, D1, D2 = p
    c1, c2 = u
    dc1  = α*c1 + D1 * (c2-c1)
    dc2 = β*c2  + D2 * (c1-c2)
    push!(var, c2)
    du .= [dc1, dc2]  # note the broadcast operator!
end

u0 = rand(2)
tspan = (0.,10.)
pars = [-.2, -.1, 1., 1.]
prob = ODEProblem(coupledd!, u0, tspan, pars)
sol = solve(prob, Tsit5())
println("size(var)= ", size(var))
println("nb time steps= ", size(sol.t))

size(var) returns 696
size(sol.t) return 12

However, when I print var, I get the output:
87-element Array{Any,1}

which seems to imply that Tsit5() calls the right-hand side 8 times every time step.

What I am interested in is the value of a variable at every time step (once per time step). So I would like to monitor variables. Is there a facility for that? Thanks.

Gordon

2

You can either specify the keyword saveat, or just call the solution using the interpolator, e.g.

sol(t) # Gives you the solution at time `t` through interpolation of the same order as the solution method

instead of sol.t

3

I was not clear. My mistake for choosing to monitor the actual solution.

Say I defined the variable z=c1+c2 in the coupledd!!() function and wish to store z in a vector at the same time levels as sol.t . Of course, I could recompute z outside the function, but why should I, especially if z is expensive to compute. I just asked a minimal example.

Thanks.

4

You can use a DiscreteCallback for this, and have it fire after every step by making it true. There’s actually already a premade callback in the library for this exact case, the SavingCallback.

https://docs.juliadiffeq.org/latest/features/callback_library/#saving_callback-1

Or you can use the integrator interface for stepwise control of the solution:

https://docs.juliadiffeq.org/latest/basics/integrator/

5

Thank you, Chris. My function is much more complicated. But it does appear that I must recompute the variables of interest, and redefine them in a new function in order to do the monitoring. That violates the DRY principle, and can be error-prone. It would be nice to tag the variables of interest in the function defining the right hand side and have the callbacks generated automatically. Am I making sense?

Thanks.

6

Not necessarily. You’ll only be re-evaluating at the final u if the method is FSAL, which many methods are not (though B-stable methods have it). So the callback, i.e. the value at t + dt, will be a brand new u. You could then Ref these over to p if you really wanted to. I can throw an example together if you think it’s worth it.

But I am not sure it is. You would decrease your calculations by at most 1/8 if this is Tsit5 for example. If this is something for stuff equations, this would be even less. For example, with Rodas5 this would be saving 1/(N+9) where N is the number of states in the system (or if sparse, maximum(colorvec)). I think the case where this factor is largest is BS3 where it is an FSAL algorithm but 4 stages, so you’d reduce the calculation of the explicit algebraic variables by 1/5. So in a worst case scenario, where the calculation of the explicit algebraic variables is 100% of your code’s time and you use BS3, then you at most have a 20% slowdown from this. In total, there are probably a million other things to optimize first.

Is there an interface that you’re used to that allows/tracks explicit algebraic variables in a nice way?

7

Brian2 is such an interface, using Python. I may have to provide an example. Your idea of using p is what I need I think, if by p you mean the parameter array. I am using these variables for debugging purposes with visualization. It is not a out thecomputer time, but about not wishing to introduce new errors.i do not wish to defi e complex e pressions undergoing structural change in two different locations.

Thanks.

Gordon

8

Chris, here are some more details.

function myCallback(u, t, integrator)
	C  = @view(u[1:1:2])  # Calcium
    Ce = @view(u[3:1:4])  # Calcium in ER
    I  = @view(u[5:1:6])  # IP3
    h  = @view(u[7:1:8])
	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, p)
end

I am returning several variables. There are 8 equations, two equations for C, two for Ce, etc.

I setup the callbacks:

saved_values = SavedValues(Float64, Vector::Any)
cb = SavingCallback(callback, saved_values)

prob = ODEProblem(Cresswell!, u0, tspan, pars)  # d is a parameter dictionary
sol = solve(prob, Tsit5(), callback=cb)

I thought that Vector::Any was completely general in the sense that Vector elements could be any type. But that is not the case according to the first line of the error message. Here is a zoomed view:

image
image876×30 4.82 KB

Each of the return elements is an array of two elements (could be more in general). I want the capability of adding more variables without having to change the return type. So I am concerned with usability, not efficiency, since this is a debugging tool for me at this stage.

I did check out https://docs.juliadiffeq.org/release-5.0/features/callback_library.html

image
image2099×278 71.3 KB

(I cannot copy/paste this error. The part in read is not copiable. Why would that be?)

I would like to also know how bring in my parameter array “p” into the myCallback() function. I could access it as a global variable, but that is poor programming.

Thanks!

Here is the full code.

using DifferentialEquations, Plots, BenchmarkTools

function globalParametersEvan()
    # Try to duplicate Evan's model as much as possible
    d = Dict()
    d[:Rgas] = 8.31   # joule / kelvin / mole
    d[:D_I] = 5.3e-10 # Hz * meter^2
    d[:D_C] = 5.3e-10 # Hz * meter^2
    d[:D_CE] = 5.3e-10 # Hz * meter^2
    d[:Cm] = 1.e-6   # Farad * meter^-2
    d[:e] = 1.602176634e-19  # Amp * second
    d[:N_A] = 6.02214076e23 # mole^-1
    d[:k_B] = 1.38064852e-23  # meter^2 * Kg * Hz^2 / kelvin
    d[:F] = 96485.3329  # Amp * second / mole = Cb / mole
    d[:V_T] = 0.026  # voltage_cubic
    d[:P_Ca] = 4.46e-15  # meter * Hz
    d[:Crest] = 0.073e-6 # Molar
    d[:Vrest] = -80. # Volt
    d[:Ce] = 1.800  #  mole/meter**3]
    d[:Ce] = 1.e-4  #  mole/meter**3]  # WHICH ONE TO CHOOSE?

    d[:O_delta] = 0.6e-3  # Hz * mole/meter^3
    d[:k_delta] = 1.5e-3  # mole/meter^3
    d[:K_delta] = 0.1e-3  # mole/meter^3
    d[:K_D] = 0.7e-3  # mole/meter^3
    d[:K_3] = 1.0e-3  # mole/meter^3
    d[:O_3K] = 4.5e-3  # Hz*mole/meter^3
    d[:o_5P]   = 0.05   # Hz
    d[:K_5P] = 10.0e-3  # mole/meter^3
    d[:o_delta] = 0.6e-3  # mole/meter^3
    d[:k_delta] = 1.5e-3  # Hz * mole/meter^3
    d[:K_delta] = 0.1e-3   # mole/meter^3
    d[:Lambda] = 2100e-18  # meter^3
    d[:rho_A] = 0.18   # ER-to-cytoplasm volume ratio ?
    d[:rho] = d[:rho_A] / (1. + d[:rho_A])  # # ER-to-cytoplasm volume ratio ?

    # Multiply value of amplitudes by cytosolic volume (changed rho_A to rho)
    # It does not look as if I did the multiplication
    #self.o_delta = 0.6  # * umolar * Lambda * (1-rho) / second
    d[:o_delta] = 0.6e-3  # * mole/meter**3 * Lambda * (1-rho) / second
    d[:o_3K]    = 4.5e-3  # * mole/meter**3 * Lambda * (1-rho) / second

    d[:Omega_L] = 0.1   #/second         # Maximal rate of Ca^2+ leak from the ER
    d[:Omega_2] = 0.2   #/second      # IP_3R binding rate for Ca^2+ inhibition
    d[:Omega_u] = 0.2-3 #/second      # uptake/dissassociation constant

    # --- IP_3R kinectics
    d[:d_1] = 0.13e-3 #*mole/meter**3            # IP_3 binding affinity
    d[:d_2] = 1.05e-3 #*mole/meter**3            # Ca^2+ inactivation dissociation constant
    d[:O_2] = 0.2e3 #/(mole/meter**3)/second      # IP_3R binding rate for Ca^2+ inhibition
    d[:d_3] = 0.9434e-3 #*mole/meter**3          # IP_3 dissociation constant
    d[:d_5] = 0.08 #*mole/meter**3            # Ca^2+ activation dissociation constant

    d[:p_open] = 1
    d[:P_r]    = 1e-3 #*mole/meter**3/second
    d[:P_CE]    = 1e-3 #*mole/meter**3/second

    d[:eta_p] = 1
    d[:d_ER]  = 1e-6  #* meter
    d[:dv_ER] = 10e-3
    d[:K_P] = 0.05e-3

    d[:L] = 8e-6  # umeter
    d[:R] = 3e-6  # meter
    d[:r] = 2e-6  # meter
    d[:O_delta] = 0.6e-21  # mole/meter**3 * Lambda * (1-rho) / second
    d[:O_5P] = 0.05e-21   # mole/meter**3 * Lambda * (1-rho) / second
    d[:F] = 0.1  # ALREADY DEFINED?
    d[:D] = 0.05

	# Evan parameters
	d[:vIPR] = 6  # sec^-1
	d[:vSERCA] = 4.4e-3  # mole m^-3 s^-1 (1 Molar = 1 Mole/liter = 10^3 mole/meter^3)
	d[:vLeak] = 0.11 # sec^-1
	d[:d1] = 0.13e-3 # mole m^-3
	d[:d5] = 0.08234e-3 # mole m^-3
	# open channels
	d[:d2] = 1.049e-3 # mole m^-3
	d[:d1] = 0.13e-3  # mole m^-3
	d[:d3] = 0.9434e-3 # mole m^-3
	d[:a2] = 0.20e-3 # mole m^-3 s^-1
	# IP3 (J_δ)
	d[:o_δ] = .025e-3 # Mole m^-3 s^-1
	d[:K_δ] = 0.5e-3   # Mole m^-3, Ca2+ affinity of PLCδ
	d[:k_δ] = 1.e-3  # Mole m^-3, IP3 inhibiting affinity   // Oschmann2009: 1.5 | COMPGLIO pg130: 1.
	# End Evan Parameters
    return d
end

lJp = []
lJr = []
lJ1 = []
lJleak = []
lJSERCA = []
lJIPR = []
lJβ = []
lCount = []
lTime = []

function getCurrents(C, Ce, I, h, t, p)
	minf    = @. Hill(I,p[:d1],1) * Hill(C,p[:d5],1)
	Jleak   = @. p[:vLeak] * (Ce-C)
	JIPR    = @. p[:vIPR] * (minf*h)^3 * (C-Ce)
	JSERCA  = @. p[:vSERCA] * Hill(C,p[:d5],2)
	J_β     = sin(2*π*t)^10  #*mole/meter**3/second  : mole/meter**3/second
	J_δ     = @. p[:o_δ]  * Hill(p[:k_δ], I, 1) * Hill(C, p[:K_δ], 2)
	J5P     = @. p[:o_5P] * Hill(I, p[:K_5P], 1)  #: mole/second # o_5P, K_5P
	J3K     = @. p[:o_3K] * Hill(C, p[:K_D], 4) * Hill(I, p[:K_3], 1) #: mole/second # o_3K, K_D, K_3
	Q2      = @. p[:d2] * (I + p[:d1]) / (I + p[:d3])
	OmegaH  = @. p[:a2] * (Q2 + C)
	hinf    = @. p[:d_2] * (I + p[:d_1]) / (p[:d_2]*(I + p[:d_1]) + (I + p[:d_3])*C) #: 1 # d_2, d_1, d_3
	return J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2
end

function Cresswell!(du, u, p, t)
    C  = @view(u[1:1:2])  # Calcium
    Ce = @view(u[3:1:4])  # Calcium in ER
    I  = @view(u[5:1:6])  # IP3
    h  = @view(u[7:1:8])  # fraction of open channels
    push!(lCount, 1)

    # rhs_C, etc. Are arrays, created via implicit vectorization involving numpy arrays
    dR2 = p[:R]^2 - p[:r]^2  #: meter**2
    s = p[:R] + p[:r]
    V = @. p[:Vrest] + (p[:F]*s/p[:Cm]) * (C - p[:Crest]) #: volt
    VVT = -2. * V/p[:V_T] #: 1  # value of 0, which leads to a zero denominator

	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, p)
 	p[:J_δ] = J_δ
    push!(lJleak, Jleak)
    push!(lJIPR, JIPR)
    push!(lJSERCA, JSERCA)
    push!(lJβ, J_β)
	push!(lTime, t)

    dh  = @. OmegaH * (hinf - h) #: 1
    dI  = @. (J_β + J_δ - J3K - J5P)

    # GE: divde Jp by Lambda (2020-03-09)
    dC  = @.                     (JIPR - JSERCA + Jleak)  #: mole / meter**3
    dCE = @. -(1. / p[:rho_A]) * (JIPR - JSERCA + Jleak)  #: mole / meter**3
    dI  = @. (1. *J_β + 1. * J_δ - 1. *J3K - 1. *J5P) # / (p[:Lambda]*(1-p[:rho])) #: mole/meter**3
    dh  = @. OmegaH * (hinf - h) #: 1

    du[1:2] = dC
    du[3:4] = dCE #[0,0]
    du[5:6] = dI # [0,0]  # blows up!
    du[7:8] = dh
end

nb_neurons = 2
eqs_per_neuron = 4
# Initial Conditions
function initialConditions()
    c0 = rand(nb_neurons)
    ce0 = rand(nb_neurons)
    h0 = rand(nb_neurons)
    I0 = rand(nb_neurons)
    # Initial condition vector
    u0 = vcat(c0,ce0,h0,I0)
end

function callback(u, t, integrator)
	C  = @view(u[1:1:2])  # Calcium
    Ce = @view(u[3:1:4])  # Calcium in ER
    I  = @view(u[5:1:6])  # IP3
    h  = @view(u[7:1:8])
	# How to I get my parameter array as a local variable? 
	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, pars)
end

getParams = globalParameters_mod.globalParametersEvan

saved_values = SavedValues(Float64, Vector::Any)
cb = SavingCallback(callback, saved_values)

u0 = initialConditions()
pars = getParams()
tspan = (0., 300.)  # make these real
#u0 = [.1, .2, .3, .4]``
# We have 8 equations. How to collect them?
prob = ODEProblem(Cresswell!, u0, tspan, pars)  # d is a parameter dictionary
sol = solve(prob, Tsit5(), callback=cb)

Here is a modified version of the callback unction:

function callback(u, t, integrator)
	C  = @view(u[1:1:2])  # Calcium
    Ce = @view(u[3:1:4])  # Calcium in ER
    I  = @view(u[5:1:6])  # IP3
    h  = @view(u[7:1:8])
	# How to I get my parameter array as a local variable?
	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, pars)
	cur = Dict()
	cur[:J5P] = J5P
	cur[:J3K] = J3K
	cur[:J_β] = J_β
	cur[:J_δ] = J_δ
	cur[:Jleak] = Jleak
	cur[:JIPR] = JIPR
	cur[:JSERCA] = JSERCA
	cur[:hinf] = hinf
	cur[:OmegaH] = OmegaH
	cur[:minf] = minf
	cur[:Q2] = Q2
	return cur
end

saved_values = SavedValues(Float64, Dict{Any,Any})
cb = SavingCallback(callback, saved_values)

and this “seems” to work. I do not yet know how to use this callback once the solver has completed.
But my immediate comment is that my callback is now substantially longer than it was and more unwieldy. I could have populated a dictionary in the getCurrents() method, but then I would have to access these dictionary values in the ODE solver. I assume that accessing dictionary values is very cheap?

I have another question. If I had an array:

vars = [:J5P, :J3K]

for example, would it be possible to create a macro to execute:

cur[:J5P] = J5P
cur[:J3K] = J3K

?

Thanks,

9

It’s Vector{Any}

10

I was close. Yes, it works. I still do not grok Julia types. But so far so good. Now, how to use the results from the callback.

11

T{X} is the family of T “of” X. I.e. Vector of Float64 is Vector{Float64}. :: is used for declarations, like x::T i.e. “x is a T”.

12

So in your simple example case, you have:

\frac{dc_1}{dt} = \alpha c_1 + D_1(c_2-c_1) \\ \frac{dc_2}{dt} = \beta c_2 + D_2(c_1-c_2)

And you want to compute, say z = c_1 + c_2, in addition, with the same time resolution as for c_1 and c_2? Without having to recompute z in post processing?

Why not pose the problem as a DAE? So you define an error vector e as:

e_1 = -\frac{dc_1}{dt} + \alpha c_1 + D_1(c_2-c_1) \\ e_2 = - \frac{dc_2}{dt} + \beta c_2 + D_2(c_1-c_2) \\ e_3 = -z + c_1 + c_2

and aim to make the errors e_1, e_2, e_3 zero.

The following code does this (I use err instead of out — out is used in the documentation of DifferentialEquations.jl):

using DifferentialEquations, Plots, BenchmarkTools
using Sundials # package of DAE solver(s)

#var = []
function coupledd!(err,du, u, p, t)
    α, β, D1, D2 = p
    #c1, c2 = u
    err[1] = -du[1] + α*u[1] + D1 * (u[2]-u[1])
    err[2] = -du[2] + β*u[2]  + D2 * (u[1]-u[2])
    err[3] = -u[3] + u[1] + u[2]
#    push!(var, c2)
    #du .= [dc1, dc2]  # note the broadcast operator!
end

u0 = rand(3)
u0[3] = u0[1] + u0[2] # making sure that the initial values u0 satisfies the equations
du0 = zeros(3) # Not bothering about checking that du0 is correct or not
d_vars = [true, true, false] # Specifying which of the u variables are differentiable
tspan = (0.,10.)
pars = [-.2, -.1, 1., 1.]
prob = DAEProblem(coupledd!, du0, u0, tspan, pars,differential_vars=d_vars)
sol = solve(prob,IDA());

You can finally plot the result:

plot(sol,vars=[(0,1),(0,2),(0,3)])

which looks as follows:

image
image1442×969 132 KB

Is this what you wanted to achieve?

13

For this simple example, you are correct. But I posted the actual code a few messages back. The DAE approach is not really what I want. My expression z is in reality auxiliary variables used to define du. I simply would like to easuly specify a list of variables to monitor, with an expression such as:

Mon = monitor([t, var1, var2, var3])

Then, solve the equation with solve(), and then I should be to plot using

Plot(Mon.t, Mon.var1)

There seems to be nothing around that can do this. To add a variable to monitor, simply add a variable to the list. Brian2 has this capability in Python.

14

The ODE solvers that Brian2 uses (GSL) don’t have this capability either, so whatever workaround you’re doing is similar to what Brian2 is doing. Going to a DAE might be overkill unless the problem is already stiff (since then you have to use an implicit method), but putting the algebraic part inside of a function and just calling that function in two spots is DRY without performance issues (and most likely exactly what they are doing).

There is an issue to investigate adding this to the solvers, but there is are good reasons for not supporting it directly in the solvers too, and leaving that as a user-level detail.

1 Like
15

I could imagine adding a line inside the rhs routine stating the names of the variables to be monitored along with an indication ad to whether they should be tracked on entry or on exit. In Brian, the equations are written inside a string that is parsed, and so there is more flexibility. The code is translated to Python or C.

The callback mechanism is good for the examples shown, but not for my use case. They are unwieldy when monitoring many variables and it is more difficult than it could be. The notation Mon.var returning a time series vector is very useful.

I enjoyed the discussion and will think some more.

Gordon

16

That seems to me like a modeling layer on top of the differential equation solver (like how Brian2 has it), rather than something inside of a differential equation solver itself. It sounds like something that could have a million variations for every discipline to do something differently. Also anything with “tracking” requires careful compile time work or is slow, so that’s very likely not going to be in the straight solvers which in theory should be something that can compile down to Float64 code to deploy to other languages.

1 Like
17

You make a good point, Chris . However, a good the beauty of ajulia should combin flexibility and ease of use. When my kind of callbacks or monitoring are used, for debugging, ease of use is Paramount. In other cases, efficiency is paramount. Thanks for the insight!

18

I think I can get closer to what I want by using integrate seperately, and integrate step by step, with fixed time step, for debugging, and after every step, place the monitored variables in some type of dataframe.

I just came across the following blog that provides a high level overview of tabular data, which might help me: DataFrames Archives | juliabloggers.com

19

Maybe a few years too late, but I found the same problem. My solution was this

using DifferentialEquations
tspan = (0.0,10.0)
dt = 0.1
sampletimes = tspan[1]:dt:tspan[2]

varz = [ ]
function f(du,u,p,t)
    z = -2*u[1]
    if t in sampletimes
	push!(varz,z)
    end
    du[1] = z
end
u0 = [10.0]
prob = ODEProblem(f,u0,tspan,p=nothing;saveat=dt,tstops=sampletimes)
sol = solve(prob,Tsit5())

#this removes the data at duplicate times, even though I do not know why it appears
k = Int(size(varz)[1]/size(sampletimes)[1])
varz = varz[k:k:size(varz)[1]]
1 Like
20

Hi Luis. Better late than never.
I have got a similar problem, I’m a new Julia user but matlab expert(I am tired of paying them).

I want to call back certain values on my ode as follows from matlab

function [eqn,Qgap,Force,Fr]=navier(t,y)

%parameters
rho=810;
beta=5e8;
r1=0.01;
r2=0.004;
Ap=pi*(0.01^2-0.004^2);
L=0.06;
a=0.01;
f=10;
l=0.01;
Aeff=pi*2*r1*l;
delta=0.007;
N=20;
dx=delta/N;
rad=linspace(r1,r1+delta,N);

vel=y(1:N);
p1=y(N+1);
p2=y(N+2);

x=a*sin(2*pi*f*t);
u=a*2*pi*f*cos(2*pi*f*t);

eqn=zeros(N+2,1);
vel(1)=u;
vel(end)=0;
dudx=[(vel(2)-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-vel(end-1))/dx];
dudx2=(vel(3:end)-2*vel(2:end-1)+vel(1:end-2))/(dx*dx);
dudx2 = [(dudx(1)-dudx(2))/dx;dudx2;(dudx(end-1)-dudx(end))/dx];
mu=[630 * (1+(0.05*(dudx)).^2).^(-0.23)];
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23);

Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta);
% Qgap=trapz(rad,2*pi*transpose(rad).*vel);

eqn(1)=(p1-p2+pf)/l + ((mu(1))/rho)*dudx2(1);
eqn(2:N-1)=(p1-p2+pf)/l + ((mu(2:N-1))./rho).*dudx2(2:N-1);
eqn(N)=(p1-p2+pf)/l + ((mu(N))/rho)*dudx2(N);
eqn(N+1)=(beta/(Ap*(L-l-x)))*(-Qgap+Ap*u);
eqn(N+2)=(beta/(Ap*(L-l+x)))*(Qgap-Ap*u);

Force=(p1-p2+pf)*Ap;
Fr=mu(1)*Aeff*dudx(1);

    
end

In matlab, we just have to put the variables as outputs after the eqn, then write a callback function (a bit complicated version of my own)

[~,Qgap,Force,Fr] = cellfun(@(tsol,ysol)  navier(tsol,ysol.'), num2cell(tsol), num2cell(ysol,2),'uni',0);
Qgap=cell2mat(Qgap);
Force=cell2mat(Force);
Fr=cell2mat(Fr);

where as it gives me back a cell of the values with the time points, and then convert to a matlab simple array.

Can you help me achieve this in Julia?