DifferentialEquations.jl package, variable rate jumps with complex variables

General Usage
,
1

Hello, I’m trying to get the following example to work:

using DifferentialEquations

function f(dx, x,p,t)
    dx[1] = x[1]
end

function rate(x,p,t)
    t
end

function affect!(integrator)
    integrator.u[1] = integrator.u[1] * 0.5
end

function main()
    jump = VariableRateJump(rate, affect!)

    x₀ = 1.0 + 0.0im
    Δt = (0.0, 6.0)
    prob = ODEProblem(f, [x₀], Δt)
    jumpProblem = JumpProblem(prob, Direct(), jump)
    sol = solve(jumpProblem, Tsit5())
end

main()

For some reason changing x₀ from real to imaginary, as in here by adding the 0.0im, breaks the code. What could be the problem? If I would like to work with complex numbers and variable rate jumps, what should I do? Curiously changing the type of the jump from variable rate to constant rate works, even with complex numbers.

I understand that I could possibly circumvent this problem by decomposing my complex matrix to a vector with real numbers. But I am curious, is this intended behaviour or is it somekind of a bug?

Here is the error:

ERROR: MethodError: no method matching randexp(::Random.TaskLocalRNG, ::Type{ComplexF64})
Closest candidates are:
  randexp(::Random.AbstractRNG, ::Type{T}, ::Tuple{Vararg{Int64, N}} where N) where T at E:\Julia-1.7.2\share\julia\stdlib\v1.7\Random\src\normal.jl:235
  randexp(::Random.AbstractRNG, ::Type{T}, ::Integer, ::Integer...) where T at E:\Julia-1.7.2\share\julia\stdlib\v1.7\Random\src\normal.jl:238
  randexp(::Random.AbstractRNG) at E:\Julia-1.7.2\share\julia\stdlib\v1.7\Random\src\normal.jl:118
  ...
Stacktrace:
  [1] randexp!
    @ E:\Julia-1.7.2\share\julia\stdlib\v1.7\Random\src\normal.jl:196 [inlined]
  [2] randexp!
    @ E:\Julia-1.7.2\share\julia\stdlib\v1.7\Random\src\normal.jl:232 [inlined]
  [3] reset_jump_problem!(jump_prob::JumpProblem{true, ODEProblem{ExtendedJumpArray{ComplexF64, 1, Vector{ComplexF64}, Vector{ComplexF64}}, Tuple{Float64, Float64}, true, 
SciMLBase.NullParameters, ODEFunction{true, DiffEqJump.var"#jump_f#130"{ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, JumpSet{Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Tuple{}, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Direct, CallbackSet{Tuple{ContinuousCallback{DiffEqJump.var"#151#153"{Int64}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}}, Tuple{}}, Nothing, Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Nothing, Nothing}, seed::Nothing)
    @ DiffEqJump C:\Users\Work\.julia\packages\DiffEqJump\b19T5\src\solve.jl:61
  [4] #__init#156
    @ C:\Users\Work\.julia\packages\DiffEqJump\b19T5\src\solve.jl:19 [inlined]
  [5] __init(_jump_prob::JumpProblem{true, ODEProblem{ExtendedJumpArray{ComplexF64, 1, Vector{ComplexF64}, Vector{ComplexF64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, DiffEqJump.var"#jump_f#130"{ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, JumpSet{Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Tuple{}, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Direct, CallbackSet{Tuple{ContinuousCallback{DiffEqJump.var"#151#153"{Int64}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}}, Tuple{}}, Nothing, Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Nothing, Nothing}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}})
    @ DiffEqJump C:\Users\Work\.julia\packages\DiffEqJump\b19T5\src\solve.jl:17
  [6] #init_call#37
    @ C:\Users\Work\.julia\packages\DiffEqBase\U3LtB\src\solve.jl:121 [inlined]
  [7] init_call
    @ C:\Users\Work\.julia\packages\DiffEqBase\U3LtB\src\solve.jl:108 [inlined]
  [8] #init#38
    @ C:\Users\Work\.julia\packages\DiffEqBase\U3LtB\src\solve.jl:134 [inlined]
  [9] init
    @ C:\Users\Work\.julia\packages\DiffEqBase\U3LtB\src\solve.jl:126 [inlined]
 [10] __solve(jump_prob::JumpProblem{true, ODEProblem{ExtendedJumpArray{ComplexF64, 1, Vector{ComplexF64}, Vector{ComplexF64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, DiffEqJump.var"#jump_f#130"{ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, JumpSet{Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Tuple{}, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Direct, CallbackSet{Tuple{ContinuousCallback{DiffEqJump.var"#151#153"{Int64}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}}, Tuple{}}, Nothing, Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Nothing, Nothing}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqJump C:\Users\Work\.julia\packages\DiffEqJump\b19T5\src\solve.jl:6
 [11] __solve (repeats 5 times)
    @ C:\Users\Work\.julia\packages\DiffEqJump\b19T5\src\solve.jl:6 [inlined]
 [12] #solve#44
    @ C:\Users\Work\.julia\packages\DiffEqBase\U3LtB\src\solve.jl:205 [inlined]
 [13] solve(prob::JumpProblem{true, ODEProblem{ExtendedJumpArray{ComplexF64, 1, Vector{ComplexF64}, Vector{ComplexF64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, DiffEqJump.var"#jump_f#130"{ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, JumpSet{Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Tuple{}, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Direct, CallbackSet{Tuple{ContinuousCallback{DiffEqJump.var"#151#153"{Int64}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, DiffEqJump.var"#152#154"{Int64, VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}}, Tuple{}}, Nothing, Tuple{VariableRateJump{typeof(rate), typeof(affect!), Nothing, Float64, Int64}}, Nothing, Nothing}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ DiffEqBase C:\Users\Work\.julia\packages\DiffEqBase\U3LtB\src\solve.jl:205
 [14] main()
    @ Main e:\KouluHommia\Gradu\Koodeja\JumpProblemExample.jl:22
 [15] top-level scope
    @ e:\KouluHommia\Gradu\Koodeja\JumpProblemExample.jl:26

Addional context:
I am planning on solving differential equations similar to
dp = Xdt + YdN,
where dp, X and Y are complex valued matrices and dN is the increment of a counting process. The rate of the jumps depends on the state of the system, or more precisely on p.

2

This looks like a bug in the code that handles solutions in the variable rate case. Can you open an issue at DiffEqJump.jl (Issues · SciML/JumpProcesses.jl · GitHub)?

@ChrisRackauckas I’m not so familiar with how it works, but it looks like this is related to the extended u0 that is constructed?

1 Like
3

I opened an issue on DiffEqJump about this. This might involve issues involving interfaces of several packages, so could take a bit to fix. Sorry about that!

1 Like
4

Thanks for the help! I hope fixing it won’t cause too much of a headache. Thanks for the quick reply!

5

Turned out to be easier to fix than I originally thought. There should be a new DiffEqJump release, version 8.4, shortly with the fix.

2 Likes
6

@isaacsas hello again! I tried out the update and now the example works. But when I tried to do a calculation with a complex matrix I ran into an issue.

So let’s say I have a rate function with the arguments rate(u, p, t). Now I am getting an error due to u being the type of an ExtendedJumpArray. Is this intended? I haven’t ran into this before and I thought that u was simply the value of the thing I am calculating at the given moment. I tried looking through the DifferentialEquations documentation but I couldn’t find a definition for what u exactly is in the rate function. I guess u has the value of my matrix somewhere, how can I grab it from it?

If you need, I can type out somekind of an minimal reproducible example.

7

A MWE to add to our tests would be great. I think I see the issue though so will try to put a fix together later today.

8

Just out of curiosity, what is the application you are interested in that is using complex matrices for the state? Something QM related?

9

Yes exactly! This is connected to continuous measurement and more precisely on observing the emission of photons from an atom, and how this observation affects the state of the atom. Thanks for asking :slight_smile:

Please shoot me a message tomorrow if you need the MWE to figure this one out. I’ll provide it then.

1 Like
10

Try version 8.4.1 later today when it is released and see if that fixes it. If that still doesn’t work please open an issue on DiffEqJump with a MWE we can use for testing. Thanks!

1 Like
11

I ran into something new, so I opened an issue on it here https://github.com/SciML/DiffEqJump.jl/issues/229. Thanks for your help!

12

Here are some things I had trouble with incase someone else is trying to get variable rate jumps to work and stumbles here.

VariableRateJumps change the type of integrator.u to the type of ExtendedJumpArray, so you need to use integrator.u.u to get your matrices. Also it changes the type of the solution to ExtendedJumpArray so that your sol.u is now a vector of ExtendedJumpArrays. Hope this someone!

Also, saveat seems to be broken for VariableRateJumps.