mwe for sbml test suite case 00375

mwe_case_00375.jl

################################################
# Critical dependencies from Manifest.toml:
# [[SBML]]
# deps = ["DocStringExtensions", "IfElse", "Libdl", "SBML_jll", "SparseArrays", "Symbolics", "Unitful"]
# git-tree-sha1 = "79cff18661b3a649ba852e23d70c8e9415f933f3"
# repo-rev = "pl-extensivemath_compID"
# repo-url = "https://github.com/paulflang/SBML.jl.git"
# uuid = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
# version = "0.10.1"

# [[SBMLToolkit]]
# deps = ["Catalyst", "MathML", "SBML", "SymbolicUtils"]
# git-tree-sha1 = "5a99f7d036c414ee51c4b9276359a32e6d8dd2eb"
# repo-rev = "pl-handle_boundarycondition_constant"
# repo-url = "https://github.com/SciML/SBMLToolkit.jl.git"
# uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
# version = "0.1.13"
################################################

using SBML, SBMLToolkit, Catalyst

using Downloads, DifferentialEquations, Sundials
using CSV, DataFrames
using Plots

case = "00369"  # Works nicely
case = "00375" # 12345ERROR: ArgumentError: unable to check bounds for indices of type Nothing
# 00375 fails, as S2 has boundary_condition=true, is assigned to `2`, gets eliminated in `structural_simplify` so that the event cannot assign `1` to it.
logdir = joinpath(@__DIR__, "logs")

function setup_settings_txt(fn)
    ls = split(fn, "\n")
    spls = split.(ls, ": ")
    filter!(x->length(x) == 2, spls)
    Dict(map(x -> x[1] => Meta.parse(x[2]), spls))
end

function to_concentrations(sol, ml)
    volumes = [1.]
    sol_df = DataFrame(sol)
    for sn in names(sol_df)[2:end]
        if haskey(ml.species, sn[1:3-end])
            spec = ml.species[sn[1:end-3]]
            comp = ml.compartments[spec.compartment]
            ic = spec.initial_concentration
            isnothing(ic) ? push!(volumes, 1.) : push!(volumes, comp.size)
        end
    end
    sol_df./Array(volumes)'
end

base_url = "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/$case/$case"
sbml_url = base_url*"-sbml-l3v2.xml"
settings_url = base_url*"-settings.txt"
results_url = base_url*"-results.csv"

sbml = String(take!(Downloads.download(sbml_url, IOBuffer())))
settings = String(take!(Downloads.download(settings_url, IOBuffer())))
results = String(take!(Downloads.download(results_url, IOBuffer())))

# Read results
settings = setup_settings_txt(settings)
results = CSV.read(IOBuffer(results), DataFrame)
ts = results[:, 1]  # LinRange(settings["start"], settings["duration"], settings["steps"])

# Read SBML
# SBMLToolkit.checksupport(sbml)
ml = readSBMLFromString(sbml, doc -> begin
            set_level_and_version(3, 2)(doc)
            convert_simplify_math(doc)
        end)
print(1)

# Convert to ReactionSystem
rs = ReactionSystem(ml)
print(2)

# Convert to ODESystem
sys = convert(ODESystem, rs; include_zero_odes=false, combinatoric_ratelaws=false)
if length(ml.events) > 0
    sys = ODESystem(ml)
end
n_dvs = length(states(sys))
n_ps = length(parameters(sys))
print(3)

# Convert to ODEProblem
ssys = structural_simplify(sys)
print(4)

prob = ODEProblem(ssys, Pair[], (settings["start"], Float64(settings["duration"])); saveat=ts)
print(5)

# Solve ODEProblem
sol = solve(prob, Sundials.CVODE_BDF())
println(6)

# Evaluate results
sol_df = to_concentrations(sol, ml)
idx = [sol.t[i] in results[:, 1] ? true : false for i in 1:length(sol.t)]
sol_df = sol_df[idx, :]
cols = names(sol_df)[2:end]
cols = [c for c in cols if c[1:end-3] in names(results)]
res_df = results[:, [c[1:end-3] for c in cols]]
solm = Matrix(sol_df[:, cols])
resm = Matrix(res_df)
res = isapprox(solm, resm; atol=1e-9, rtol=3e-2)
println("isapprox: $(res).")
atol = maximum(solm .- resm)

plt = plot(ts, solm)
plt = plot!(ts, resm, linestyle=:dot)

rm(logdir,recursive=true)
mkdir(logdir)
savefig(joinpath(logdir, case*".png"))
plt