integrator-using functions should return NaN on non-successful integrator steps

[Funkerfish]

Hello again!

Now I try to calculate the Lyapunov spectrum in a two-parameter plane. But a warning comes when I choose some of the parameter regions:

Warning: Instability detected. Aborting

From what I googled, the warning comes because it returns a NAN in the trajectory. I am thinking if I can use the 'isnan' function to skip those parameter sets that cause the NAN and continue my loop. Since the lyapunovspectrum function is built-in, I am looking for your help.

Here is my code:

using DynamicalSystems, CairoMakie, GLMakie
using OrdinaryDiffEq
using DelimitedFiles
using ProgressMeter
function chaotictest(u, p, t)
    k1, k2 = p
    x, y, z = u
    dx = y
    dy = z
    dz = k1-7*y+x^2+k2*z
    return SVector(dx, dy, dz)
end
u0 = [0.01, 0.01, 0.01]
para = [-1.0, -1.0] 
ds = ContinuousDynamicalSystem(chaotictest, u0, para)
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.001)
as = -1:0.1:8;  
bs = -1.0:0.05:0.01;
nloop = length(as)*length(bs)
cordiλ = zeros(nloop, 2+3)
λs = zeros(nloop, 3)
paras = zeros(nloop, 2)
cirnum = 0
prog = Progress(nloop; showspeed=true)
for (i, a) in enumerate(as)
    for (j, b) in enumerate(bs)
         global cirnum = cirnum + 1
    set_parameter!(ds, 1, a)
    set_parameter!(ds, 2, b)  
    λs[cirnum, :] .= lyapunovspectrum(ds, 2000; Ttr=3000, diffeq)
    paras[cirnum, :] .= [a, b]
    global cordiλ = [paras λs]
    next!(prog)
end
end

[Funkerfish]

when the diction of parameters is set as:

as = 2:0.1:8;  
bs = -0.9:0.05:0.01;

the warning comes immediately.

[Datseris]

what does lyapunovspectrum return for these cases that there is an instability...?

[Funkerfish]

it seems that the loop would never end when the instability warning comes, I don't know but I guess it would be the NaN, so as the function trajectory() returns.

please check the following codes with the warning：

using DynamicalSystems, CairoMakie, GLMakie
using OrdinaryDiffEq
function chaotictest(u, p, t)
    k1, k2 = p
    x, y, z = u
    dx = y
    dy = z
    dz = k1-7*y+x^2+k2*z
    return SVector(dx, dy, dz)
end
u0 = [0.01, 0.01, 0.01]
para = [2, -0.9] #k2=k6=0
ds = ContinuousDynamicalSystem(chaotictest, u0, para)
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.001)
test111 = trajectory(ds, 20; Ttr=0, diffeq)

I think the warning may be solved by changing the solver. But since it works in some of the parameter sets, we can't do that during a loop, and that is why I need your help.

[Datseris]

it seems that the loop would never end when the instability warning comes, I don't know but I guess it would be the NaN

The answer to my question is not clear though. What is lyapunovspectrum doing if you give it an initial condition you know will give you "instability detected, aborting"? Is the function never returning? is it returning NaNs? Is it returning something else?

[Funkerfish]

lyapunovspectrum never returns once the warning comes. I check the matrix which saves the values and find that it still be zero (already pre-assign 0 values to the matrix)

[Datseris]

Okay. So we need to make lyapunovspectrum terminate when the integrator detects instability. What do you think the lyapunovspectrum function should return when this happens?

[Datseris]

If you want to do a PR that solves this, check for the success field in the integrator sturcture, which becomes false if the step wasn't successful

[Funkerfish]

Because this problem comes from the ODE solver, the same warning would take place in other functions like Interactive Trajectory Evolution when we slide the bar of parameters.

I think these functions should return 'NaN' value so that we can create a filter by applying the isnan function. In old times (when using Fortran), we would let ODE solver do calculations for few steps and check if there are 'Nan' values. If so, we can skip to the next parameter set.

Thank you for your advice, I would try success field. Your work is awesome!

[Datseris]

Because this problem comes from the ODE solver, the same warning would take place in other functions like Interactive Trajectory Evolution when we slide the bar of parameters.

Yes, that is true, but the only way to handle this, as far as I can tell, is to have such a special handling of the integrator for each individual function that uses an integrator. I don't see any other way.

I think these functions should return 'NaN' value so that we can create a filter by applying the isnan function

What's the benefit to that? If the integrator structure already has a success field, then doesn't this already give you the info you need? What is the additional benefit of making the state NaN as well as having a false return code for success?

Your work is awesome!

You are very welcome!

[Funkerfish]

Of course, the success field would work. But my point is, many people like me would use more than one coding language to do the calculation, for example, obtaining the timeseries in Julia and analyzing them in Matlab, so it is important that the structure of data matrix keeps consistent. As far as I know, Functions like ODE45 in Matlab would return 'NaN' values and it would be great if DynamicalSystems.jl do the same.

Besides, the instability thing may not happen if we choose another solver, of course It is our duty to choose the proper one but it would be great to leave a ‘NaN’ mark in the data matrix to help us to evaluate them.

[Datseris]

Okay, I didn't understand you in the first place. You were suggesting that functions like lyapunovspectrum should return NaN in the case of instability, which I completely agree with. What I originally thought is that you were suggesting that the integrator should make the state NaN, which I don't agree with. Of course, now someone has to actually do the pull requests in functions like lyapunovspectrum and make then return NaN in case of instability.

But my point is, many people like me would use more than one coding language to do the calculation, for example, obtaining the timeseries in Julia and analyzing them in Matlab, so it is important that the structure of data matrix keeps consistent.

I need to tell you that for me this is not an argument. Thankfully I agree with you for different reason, that NaN is a sensible design just by itself. But saying "MatLab does it this way" is not convincing for me, not only because I consider Julia vastly superior, but also because by itself "someone else does it different" isn't necessarily a logically convincing argument.

But I am curious. Why are you doing your timeseries analysis in Matlab instead of Julia...? What tools are missing in Julia that you have in Matlab?

[Funkerfish]

What is the additional benefit of making the state NaN as well as having a false return code for success?

To be more specific, please allow me to show you an example:
Now we have a 5x2 Matrix to store the parameters and the corresponding Lyapunov exponent (pre-assigned to 0),

parameters      Lyapunov exponents
  0.0                 0.0
  0.0                 0.0
  0.0                 0.0
  0.0                 0.0
  0.0                 0.0

The parameter set is 1 to 5 and we assume the instability comes when taking the parameter at 3 and 4. Let's see what will happen with NaN compared to false code:

parameters      Lyapunov exponents 
  1.0                 -1.2   
  2.0                 -1.1           
  3.0                  NaN               
  4.0                  NaN             
  5.0                 -0.9  

parameters      Lyapunov exponents 
  1.0                 -1.2   
  2.0                 -1.1           
  5.0                 -0.9              
  0.0                  0.0           
  0.0                  0.0  

Obviously, false will destroy the original data structure. What if the function is not lyapunovspectrum but one defined by ourselves? for example, we define one function to calculate the distance between two trajectories:

dd = x(t) - y(t)

How can we let them return NaN values? So, the final solution is to let the integrator make the state NaN.

[Funkerfish]

But I am curious. Why are you doing your timeseries analysis in Matlab instead of Julia...?

I really agree with you that Julia is superfast when doing calculation and that is the reason I turn to this language and DynamicalSystems.jl.

Just imagine that we got a 300x3 data matrix in Julia and you want to extract a few lines of values as your wish or just sort values by one row. I think it would be more convenient to open the data file in Matlab or even in Microsoft Excel. But we can't say Microsoft Excel is better than Julia, can we?

I hold the opinion that in some ways DynamicalSystems.jl should act like other classical languages so that it could be more friendly and attractive to those starters.

[Datseris]

Let's see what will happen with NaN compared to false code:

That's not what I am suggesting though... Integrator returns false on instability, but lyapunovspectrum returns NaN by detecting the instability of the integrator. In any case, since we agree on the outcome there's no reason to discuss this further. lyapunovspectrum should return NaN on trajectory instability.

Now someone has to actually do the PR that implements this.

Just imagine that we got a 300x3 data matrix in Julia and you want to extract a few lines of values as your wish or just sort values by one row. I think it would be more convenient to open the data file in Matlab or even in Microsoft Excel.

This doesn't make sense to me, why not just use the built in sort! function of Julia...? you'd literally call sort(matrix; dims = 2). From my perspective this is much more convenient than opening Excel :D

[Funkerfish]

Okey, I will look forward to this update! Would we still see the warning message in REPL before returning NaN? because its quite annoying and would block the progress bar.

Thanks for all your warm replies!

[Datseris]

Okey, I will look forward to this update!

This will take a while. I am very busy with other projects, so I won't be doing this. Perhaps you can consider doing a Pull Request? It looks like a simple change.

Would we still see the warning message in REPL before returning NaN?

Yes, this comes from the integrator instability, and I think it is very useful.

[rusandris]

I think lyapunovspectrum does return the values after the warnings, it's just that if your T is large, it takes a long time to print the same warning so many times.
Using your example, and small T:

diffeq = (alg = Tsit5(), adaptive = false, dt = 0.001)
lyapunovspectrum(ds,20;diffeq)

The output:

 Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/byoKQ/src/integrator_interface.jl:372
3-element Vector{Float64}:
 21.606034085778848
 15.725910603390146
 14.312587917169772

But in case of trajectory it never seems to return once the warning comes. I would suggest that in both cases, the function should terminate immediately and return NANs if the integrator step wasn't successful.

[Datseris]

Yes, we have all agreed that in functions that use integrator, termination should stop if integration step was not successful and NaNs are returned, so there is no need to suggest it more :P Now the point is to actually implement it. Pull Requests are very much welcomed!

[rusandris]

There is a check_error(integrator) in SciMLBase which looks at the success fields, but I couldn't get it to work with integrator made from ds.

function chaotictest!(du,u, p, t)
           k1, k2 = p
           x, y, z = u
           du[1] = y
           du[2] = z
           du[3] = k1-7*y+x^2+k2*z
       end

u0 = [0.01, 0.01, 0.01]
para = [-1.0, -1.0] 
ds = ContinuousDynamicalSystem(chaotictest!, u0, para)
integ_from_ds = integrator(ds,ds.u0)
check_error(integ_from_ds)

The output is:
ERROR: type SimpleATsit5Integrator has no field iter

Normally, using an integrator made from a prob it should work like this:

prob = ODEProblem(chaotictest!,rand(3),(0.0,1.0),para)
integ_from_prob = init(prob,Tsit5())
check_error(integ_from_prob)

The output is:
:Success

The fields for an integrator made from ds:

fieldnames(typeof(integ_from_ds))
(:f, :uprev, :u, :tmp, :tprev, :t, :t0, :tf, :dt, :dtnew, :tdir, :p, :u_modified, :ks, 
:cs, :as, :btildes, :rs, :qold, :abstol, :reltol, :internalnorm)

The fields for an integrator made from prob:

fieldnames(typeof(integ_from_prob))
(:sol, :u, :du, :k, :t, :dt, :f, :p, :uprev, :uprev2, :duprev, :tprev, :alg, :dtcache, :dtchangeable, :dtpropose, :tdir, :eigen_est, :EEst,
 :qold, :q11, :erracc, :dtacc, :success_iter, :iter, :saveiter, :saveiter_dense, :cache, :callback_cache, :kshortsize, :force_stepfail, 
:last_stepfail, :just_hit_tstop, :do_error_check, :event_last_time, :vector_event_last_time, :last_event_error, :accept_step, :isout,
 :reeval_fsal, :u_modified, :reinitialize, :isdae, :opts, :destats, :initializealg, :fsalfirst, :fsallast)

What do you suggest? Should integ_from_ds have a success field, or should we just look at the NANs in integ_from_ds.u?

[Datseris]

Haha I see the confusion. In fact, the same integrator would be made from prob or from ds, it is just that the solver algorithm default is different in each case. In ds it is SimpleATsit5 that actually does not have this field. So our default integrator doesn't have a success field, which is fine. If users want more advanced stuff they can use other integrators. In the prob case the default integrator is Tsit5, which does have such a field. The solution is to simply extend the function check_error either in DynamicalSystemsBase.jl or in SimpleDiffEq.jl. The we can use check_error in our code.

[rusandris]

Can we extend check_error by simply omitting the checks that use fields that the default SimpleATsit5 doesn't have?
I mean:

function check_error(integrator::SimpleDiffEq.SimpleATsit5Integrator)
	# This implementation is intended to be used for SimpleDiffEq.SimpleATsit5Integrator
	   if isnan(integrator.dt)
	       @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.")
	       return :DtNaN
	   end

	   if DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK(integrator.dt, integrator.u, integrator.p,
		                             integrator.t)
	       @warn("Instability detected. Aborting")
	       return :Unstable
	   end

	   return :Success
end

And for the other solvers we can leave it like it is.

[Datseris]

That seems totally fine to me. You could probably PR it directly at SimpleDiffEq.jl, and then start another PR here that actually would return NaN at the lyapunovspectrum function

[Funkerfish]

Hello again! I have another question.
Inspired by @rusandris，I use check_error to detect whether the parameter could cause the instability and then pass them to the trajectory.
My question is, should we reinit! the integrator after the check_error process? Here is my code:

using DynamicalSystems, CairoMakie, GLMakie
using OrdinaryDiffEq
function chaotictest!(u, p, t)
    k1, k2 = p
    x, y, z = u
    dx = y
    dy = z
    dz = k1-7*y+x^2+k2*z
    return SVector(dx, dy, dz)
end
u0 = [0.01, 0.01, 0.01]
para = [2, -0.9] #k2=k6=0
ds = ContinuousDynamicalSystem(chaotictest!, u0, para)
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.001)
T, Δt, Ttr = 2200.0, 0.001, 7000
k1s = 1.0:1:2.0;
  for (ik1, vk1) in enumerate(k1s)
      set_parameter!(ds, 1, vk1)
      integNaNtest = integrator(ds, u0; diffeq)
      step!(integNaNtest, 800)
      reinit!(ds)
         if check_error(integNaNtest) != :Success
            continue
         else
            set_parameter!(ds, 1, vk1)
           tr111 = trajectory(ds, T; Δt, Ttr, diffeq)
           reinit!(ds)
         end
  end

[awage]

Hi, in your example it is not necessary since you define a new integrator just for the test. This should work:

using DynamicalSystems, CairoMakie, GLMakie
using OrdinaryDiffEq
function chaotictest!(u, p, t)
    k1, k2 = p
    x, y, z = u
    dx = y
    dy = z
    dz = k1-7*y+x^2+k2*z
    return SVector(dx, dy, dz)
end
u0 = [0.01, 0.01, 0.01]
para = [2, -0.9] #k2=k6=0
ds = ContinuousDynamicalSystem(chaotictest!, u0, para)
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.001)
T, Δt, Ttr = 2200.0, 0.001, 7000
k1s = 1.0:1:2.0;
  for (ik1, vk1) in enumerate(k1s)
      set_parameter!(ds, 1, vk1)
      integNaNtest = integrator(ds, u0; diffeq)
      step!(integNaNtest, 800)
         if check_error(integNaNtest) != :Success
            continue
         else           
           tr111 = trajectory(ds, T; Δt, Ttr, diffeq)
         end
  end