DifferentialEquations.jl - how to calculate an 'assignment' and pass back to an ODE?

bgsmall · October 7, 2021, 7:24pm

Hi.

I have hit a snag in that I wish to calculate a quantity u[3] as an ‘assignment’ (not sure that is the correct terminology) within the ODE function and pass back to one of my ODE’s du[3] My code is below:

# Packages (import Pkg; Pkg.add("name of package")
	using DifferentialEquations
	using Plots 

	function perfadap(du,u,p,t)
	## Reactions
	# Reaction: id = r1
		reaction_r1 = p[1]*p[2]*u[3]; 
	# Reaction: id = r2
		reaction_r2 = p[1]*p[3]*u[1]*u[2];
	# Reaction: id = r3
		reaction_r3 = p[1]*p[4]*u[3]; 
	# Reaction: id = r4
		reaction_r4 = p[1]*p[5]*u[2];

	## Assignment	
	# assignmentRule: variable = S
		u[3] = 1*floor(tspan/p[6]);
	
	## Species	
	# Species:   id = R, name = u[1], affected by kineticLaw
		du[1] = (1/(p[1]))*(( 1.0 * reaction_r1) + (-1.0 * reaction_r2));
	# Species:   id = X = u, name = X = u[2], affected by kineticLaw
		du[2] = (1/(p[1]))*(( 1.0 * reaction_r3) + (-1.0 * reaction_r4));
	# Species:   id = S, name = u[3], involved in a rule 	
	    du[3] = u[3]; 
	end
	
	### PARAMETERS (p)
	p     = [1.0, 2.0, 2.0, 1.0, 1.0, 4.0]  # compartment vol, k1, k2, k3, k4, tau
	### TIME (tspan)
	tspan = (0.0, 20.0) #10.0
	### INITIAL CONDITIONS (u0)
	u0    = [1.0,0.0,0.0] #du[1] = R, du[2] = X, du[3] = S respectively 
    ### PROB SPECIFICATION (op)	
	op    = ODEProblem(perfadap, u0, tspan, p)
	### SYS SOLUTION (sol)
	sol   = solve(op, Tsit5())       # use Tsit5 ODE solver
	
	### PLOT
	plot(sol, title = "PerfAdap.jl") # seriestype = :scatter,

This is the full ‘ERROR’ and ‘Stacktrace’ appended below

ERROR: LoadError: MethodError: no method matching /(::Tuple{Float64, Float64}, ::Float64)
Closest candidates are:
  /(::Any, ::ChainRulesCore.AbstractThunk) at C:\Users\bsmall\.julia\packages\ChainRulesCore\1LqRD\src\differentials\thunks.jl:33
  /(::StridedArray{P, N} where N, ::Real) where P<:Dates.Period at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\Dates\src\deprecated.jl:44
  /(::Union{SparseArrays.SparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}}, false} where var"#s832"<:SparseArrays.AbstractSparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false} where var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}} where {Tv, Ti}, ::Number) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\SparseArrays\src\sparsevector.jl:1450
  ...
Stacktrace:
  [1] perfadap(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, t::Float64)
    @ Main z:\Documents\Julia\tutorial\Tyson2003_PerfAdap_V2.jl:23
  [2] ODEFunction
    @ C:\Users\bsmall\.julia\packages\SciMLBase\UIp7W\src\scimlfunctions.jl:334 [inlined]
  [3] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, 
typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), 
Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}})
    @ OrdinaryDiffEq C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\perform_step\low_order_rk_perform_step.jl:623
  [4] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:456
  [5] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:67
  [6] #__solve#471
    @ C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:4 [inlined]
  [7] __solve
    @ C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:4 [inlined]
  [8] #solve_call#42
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:61 [inlined]
  [9] solve_call
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:48 [inlined]
 [10] #solve_up#44
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:87 [inlined]
 [11] solve_up
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:78 [inlined]
 [12] #solve#43
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:73 [inlined]
 [13] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5)
    @ DiffEqBase C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:68
 [14] top-level scope
    @ z:\Documents\Julia\tutorial\Tyson2003_PerfAdap_V2.jl:43
in expression starting at z:\Documents\Julia\tutorial\Tyson2003_PerfAdap_V2.jl:43

Any and all help / pointers appreciated.

Thanks

Ben

BLI · October 7, 2021, 7:48pm

It is not quite clear to me what you try to do…

u[3] = 1*floor(tspan/p[6])

seems pretty weird – and won’t work. As a minimum, you would have to do tspan./p[6] just to make your code work. But it doesn’t seem to make sense… normally, tspan should not be available inside perfadap – and why would you take the floor of a tuple such as tspan? Also, you cannot apply the floor function to a tuple – you would need to do floor.(tspan./p[6]) – but why would you do that?

Also, assigning u[3] like you do kind of ruins the entire differential equation – at least in a conventional sense.

It would be easier to suggest what to do if you explain what you try to achieve.

ChrisRackauckas · October 7, 2021, 9:44pm

Yeah, this is not an ODE… I’m confused.

bgsmall · October 8, 2021, 8:43am

Dear BLI.

Thanks for your rapid response! The calculation I assign to u[3] is a signal (an input) that should trigger the evolution of species u[1] and u[2] in the context of the model (Both reaction_r1 and reaction_r2 are triggered by the output of u[3].

Is this (assignment) what is also called a discontinuity in respect of the system of ODE’s - i.e. the calculation is explicit with respect to time tspan but is not explicitly calculated WITHIN the system of ODE’s?

So, in summary what I want to do is calculate the quantity u[3] and pass this back to the ODE. I guess if there is a way around it, u[3] does not need to remain in the function perfadap.

Thanks again for your help, much appreciated!

Ben

BLI · October 8, 2021, 11:40am

But what do you mean by the following expression?

What is tspan meant to be here? Is it the same tspan that you use to define the simulation time span, e.g., tspan = (0,20.0)? Or do you mean the current time t? Or do you mean the final time tspan[2] = 20.0. Or what?

And: why should the value of u[3] that triggers the reactions be different from u[3] = 1*floor(tspan/p[6])?

bgsmall · October 8, 2021, 4:29pm

Hi BLI (Bernt?).

Thanks again for your persistence. tspan' in this context is the current time (as you correctly point out - t`).

I have worked on this problem further now and apart from not really understanding how to ‘assign’ the value of a derivative to the expression, i.e. u[3] = 1*floor(tspan/p[6]) I decided to try a different approach.

To that end, I create a function sig(t) that uses nested ifelse. to determine when the input is on or off and then re-use this function inside my perfadap function.

# Function that turns signal on and off
        function sig(t) 
            ifelse.(t<4,0,
                ifelse.(t>=4&&t<5,1,
                    ifelse.(t>=8&&t<9,1,
                        ifelse.(t>=12&&t<13,1,
                            ifelse.(t>=16&&t<17,1,0)
                                            )
                                     )
                              )
                       )
        end

Output

Thanks again (and to Chris R)

Ben

BLI · October 8, 2021, 5:45pm

A couple of additional comments…

It is not clear to me how you used your function sig, so it is not clear to me how you eventuelly produced the plot.
You original model function heading is function perfadap(du,u,p,t). Here, t denotes the current time, and you cannot suddenly use tspan to denote the current time — so the original code was rather odd.
Your reactions are not very clear to me. From your comments, let u[1] be the concentration of R, u[2] be the concentration of X, and u[3] be the concentration (?) of S — what are the “chemical” reactions?

vettert · October 8, 2021, 8:12pm

Maybe bgsmall is actually looking for the callback functionality of differentialequations.jl?

bgsmall · October 11, 2021, 9:04am

Hi BLI.

sig just replaces instances of u[3] in the code.
I think (?) that tspan represents the time over which I want to see the solution (i.e. the boundaries of the solution to perfadap) versus t which is current ‘model’ time.
The reactions are really just a conceptual model I’m working on.

Thanks again for your help.

Ben

bgsmall · October 11, 2021, 9:16am

Hi vettert,

Thank you for your prompt.

When I said ‘assignment’, what I meant in Julia terminology was ‘event’. This example may be of use to others:

https://diffeq.sciml.ai/stable/features/callback_functions/#DiscreteCallback-Examples

I will work on this further.

Thanks

Ben

BLI · October 11, 2021, 9:29am

So – I assume you have removed the above line, and replaced u[3] with a call to your sig function, i.e., sig(t).

function perfadap(du,u,p,t)
    # input u computed from sig(t) -- inputs are often denoted u
    u = sig(t)
	## Reactions
	# Reaction: id = r1
		reaction_r1 = p[1]*p[2]*u 
	# Reaction: id = r2
		reaction_r2 = p[1]*p[3]*u[1]*u[2];
	# Reaction: id = r3
		reaction_r3 = p[1]*p[4]*u; 
	# Reaction: id = r4
		reaction_r4 = p[1]*p[5]*u[2];

	## Assignment	
	# assignmentRule: variable = S
	#	u[3] = 1*floor(tspan/p[6]);
	
	## Species	
	# Species:   id = R, name = u[1], affected by kineticLaw
		du[1] = (1/(p[1]))*(( 1.0 * reaction_r1) + (-1.0 * reaction_r2));
	# Species:   id = X = u, name = X = u[2], affected by kineticLaw
		du[2] = (1/(p[1]))*(( 1.0 * reaction_r3) + (-1.0 * reaction_r4));
	# Species:   id = S, name = u[3], involved in a rule 	
	    du[3] = u; 
	end

In this case, sig(t) is not a state; it’s an input, thus u = sig(t). [Control engineers tend to denote the state by x, and not u, and inputs by u — see above]. Thus, adding a “state” of form du[3] = u just integrates up sig(t) – which is rarely of interest.

If you had computed the input u inside of the function instead of by an external function sig, you would probably want to store u so that you can view it after the simulation. This can not be done easily in ODEProblem formulations, but is straightforward in the slightly more complex form DAEProblem.