InexactError with TauLeaping()

[oliviaAB]

Hi,

I'm trying to simulate a rather complex model with the Tau Leaping algorithm (110 species, 216 reactions). I'm running into the following error:

ERROR: InexactError: Int64(NaN)
Stacktrace:
 [1] Type at .\float.jl:703 [inlined]
 [2] convert at .\number.jl:7 [inlined]
 [3] setindex! at .\array.jl:767 [inlined]
 [4] step!(::BioSimulator.OTL, ::Array{Int64,1}, ::BioSimulator.SparseReactionSystem) at C:\Users\oangelin\.julia\packages\BioSimulator\cW5tq\src\algorithm\otl.jl:98
 [5] simulate!(::BioSimulator.RegularPath{Float64,Int64}, ::Array{Int64,1}, ::BioSimulator.OTL, ::BioSimulator.SparseReactionSystem) at C:\Users\oangelin\.julia\packages\BioSimulator\cW5tq\src\interface\simulate.jl:118
 [6] simulate_chunk!(::Array{BioSimulator.RegularPath{Float64,Int64},1}, ::Array{Int64,1}, ::Array{Int64,1}, ::BioSimulator.OTL, ::BioSimulator.SparseReactionSystem, ::UnitRange{Int64}) at C:\Users\oangelin\.julia\packages\BioSimulator\cW5tq\src\interface\simulate.jl:108
 [7] simulate_wrapper!(::Array{BioSimulator.RegularPath{Float64,Int64},1}, ::Array{Int64,1}, ::Array{Int64,1}, ::BioSimulator.OTL, ::BioSimulator.SparseReactionSystem) at C:\Users\oangelin\.julia\packages\BioSimulator\cW5tq\src\interface\simulate.jl:90
 [8] #simulate#54(::Float64, ::Int64, ::Int64, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Network, ::TauLeaping, ::Val{:fixed}) at C:\Users\oangelin\.julia\packages\BioSimulator\cW5tq\src\interface\simulate.jl:51
 [9] #simulate at .\none:0 [inlined] (repeats 2 times)
 [10] top-level scope at none:0

Which points to line 98 in the file otl.jl:

        algorithm.events[j] = poisrand(τ * a[j])

I'm not sure how to pinpoint more precisely the source of the error. It does not happens all the time, some simulations are successful. Also simulating the same system with the Direct algorithm works fine. This means that it must be some extreme case in the simulation. Is there any easy/obvious way to correct that?

Thanks!

[alanderos91]

I can reproduce this issue with the following:

using Distributions

rand(Poisson(Inf))

Note that the example above calls poisrand internally. The propensities a[j] ought to be well-behaved, so my guess is that the leap size τ is probably the culprit.

Try starting Julia as: julia --math-mode=ieee to disable the @fastmath annotations. I may have been too eager to throw in those optimizations in otl.jl. See if that prevents the issue from happening. Maybe check that none of the propensities blow up to Inf just to make sure?

[oliviaAB]

Thank you for your answer!

Starting Julia with julia --math-mode=ieee didn't fix the error.

By looking at the propensities I figured out that a[j] becomes negative. Again, not sure if it's to be expected given that this version of TauLeaping is providing a safeguard against negative populations. I tried smaller epsilon but I'm still encountering the error.

[alanderos91]

I was afraid that might be the case. This is more difficult to diagnose. Does this also happen with StepAnticipation()?

It would be great if you can provide a MWE. If that's not possible, a description of the network you're simulating might suffice.

[oliviaAB]

Hi, I tried StepAnticipation, sometimes it works fine, other times it doesn't finish (runs for a long time and I have to terminate the run).

Here is the code I used; I saved the information about the system in the file tauleaping.jld (tauleaping.zip) using the package JLD. The function createBioSimModel comes from my R/Julia package sismonr:

using JLD
using BioSimulator

## --------------------- Function to create the model --------------------- ##

function createBioSimModel(stochmodel, QTLeffects, InitVar, modelname)

  model = BioSimulator.Network(modelname)

  ## Add the species in the model, with their initial abundance
  for i in 1:length(stochmodel["species"])
    i0 = replace(stochmodel["initialconditions"][i], "QTLeffects" => "$QTLeffects")
    i0 = replace(i0, "InitVar" => "$InitVar")
    i0 = eval(Meta.parse(i0))
    #println(stochmodel["species"][i]* "\t"*string(i0))

    if !isa(i0, Number)
      println(stochmodel["initialconditions"][i])
      error("Pb during the construction of the Julia BioSimulator model: initial abundance is not numeric.")
    end
    if typeof(stochmodel["species"][i]) != String
      println(stochmodel["species"][i])
      error("Pb during the construction of the Julia BioSimulator model: species name is not a String.")
    end

    if round(Int, i0)<0
      error("Initial condition of species $(stochmodel["species"][i]) is negative ($(i0)).")
    end

    model <= BioSimulator.Species(stochmodel["species"][i], round(Int, i0))

  end

  ## Add the reactions in the model, with their name and rate
  for i in eachindex(stochmodel["reactions"])
    prop = replace(stochmodel["propensities"][i], "QTLeffects" => "$QTLeffects")
    #println(prop)
    prop = eval(Meta.parse(prop))

    if !isa(prop, Number)
      println(stochmodel["propensities"][i])
      error("Pb during the construction of the Julia BioSimulator model: propensity is not numeric.")
    end
    if typeof(stochmodel["reactionsnames"][i]) != String
      println(stochmodel["reactionsnames"][i])
      error("Pb during the construction of the Julia BioSimulator model: reaction name is not a String.")
    end
    if typeof(stochmodel["reactions"][i]) != String
      println(stochmodel["reactions"][i])
      error("Pb during the construction of the Julia BioSimulator model: reaction is not a String.")
    end

    if prop<0
      error("Reaction rate of reaction $(stochmodel["reactionsnames"][i]) is negative ($(prop)).")
    end

    model <= BioSimulator.Reaction(stochmodel["reactionsnames"][i], prop, stochmodel["reactions"][i])

  end

  return model
end

## --------------------- Simulate the system --------------------- ##

d = load("tauleaping.jld")

QTLeffects = d["QTLeffects"]
stochmodel = d["stochmodel"]
InitVar = d["InitVar"]

model = createBioSimModel(stochmodel, QTLeffects, InitVar, "test_tau_leaping")

simalgorithm = BioSimulator.TauLeaping()

result = simulate(model, simalgorithm, time = convert(Float64, 1000), epochs = round(Int64, 1), trials = convert(Int64, 100))

Thank you for your help!

[alanderos91]

Just a quick update:

The is_badleap check fails to catch a negative population event. I am not sure why this happens with this model.

[alanderos91]

Apologies for the delay. The issue appears to be resolved in the latest release. However, that version does not have all tau-leaping methods working with Gillespie updates yet.

[oliviaAB]

Sorry for the delay in testing it out, both TauLeaping methods seem to work fine! Thanks a lot :)