diff --git a/Project.toml b/Project.toml
index 11fe28a..16c189c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,6 +24,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
OpenAPI = "d5e62ea6-ddf3-4d43-8e4c-ad5e6c8bfd7d"
Oxygen = "df9a0d86-3283-4920-82dc-4555fc0d1d8b"
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
SwaggerMarkdown = "1b6eb727-ad4b-44eb-9669-b9596a6e760f"
SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
diff --git a/examples/request-calibrate-no-integration.json b/examples/request-calibrate-no-integration.json
index 1cb9b4a..06824c2 100644
--- a/examples/request-calibrate-no-integration.json
+++ b/examples/request-calibrate-no-integration.json
@@ -1,8 +1,21 @@
{
- "model": "{\"name\": \"Giordano2020 - SIDARTHE model of COVID-19 spread in Italy\",\"schema\": \"https://raw.githubusercontent.com/DARPA-ASKEM/Model-Representations/petrinet_v0.5/petrinet/petrinet_schema.json\",\"schema_name\": \"petrinet\",\"description\": \"Giordano2020 - SIDARTHE model of COVID-19 spread in Italy\",\"model_version\": \"0.1\",\"properties\": {},\"model\": { \"states\": [ { \"id\": \"Susceptible\", \"name\": \"Susceptible\", \"grounding\": { \"identifiers\": { \"ido\": \"0000514\" }, \"modifiers\": {} }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Diagnosed\", \"name\": \"Diagnosed\", \"grounding\": { \"identifiers\": { \"ido\": \"0000511\" }, \"modifiers\": { \"diagnosis\": \"ncit:C15220\" } }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Infected\", \"name\": \"Infected\", \"grounding\": { \"identifiers\": { \"ido\": \"0000511\" }, \"modifiers\": {} }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Ailing\", \"name\": \"Ailing\", \"grounding\": { \"identifiers\": { \"ido\": \"0000511\" }, \"modifiers\": { \"disease_severity\": \"ncit:C25269\", \"diagnosis\": \"ncit:C113725\" } }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Recognized\", \"name\": \"Recognized\", \"grounding\": { \"identifiers\": { \"ido\": \"0000511\" }, \"modifiers\": { \"diagnosis\": \"ncit:C15220\" } }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Healed\", \"name\": \"Healed\", \"grounding\": { \"identifiers\": { \"ido\": \"0000592\" }, \"modifiers\": {} }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Threatened\", \"name\": \"Threatened\", \"grounding\": { \"identifiers\": { \"ido\": \"0000511\" }, \"modifiers\": { \"disease_severity\": \"ncit:C25467\" } }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } }, { \"id\": \"Extinct\", \"name\": \"Extinct\", \"grounding\": { \"identifiers\": { \"ncit\": \"C28554\" }, \"modifiers\": {} }, \"units\": { \"expression\": \"1\", \"expression_mathml\": \"1\" } } ], \"transitions\": [ { \"id\": \"t1\", \"input\": [ \"Diagnosed\", \"Susceptible\" ], \"output\": [ \"Diagnosed\", \"Infected\" ], \"properties\": { \"name\": \"t1\" } }, { \"id\": \"t2\", \"input\": [ \"Ailing\", \"Susceptible\" ], \"output\": [ \"Ailing\", \"Infected\" ], \"properties\": { \"name\": \"t2\" } }, { \"id\": \"t3\", \"input\": [ \"Recognized\", \"Susceptible\" ], \"output\": [ \"Recognized\", \"Infected\" ], \"properties\": { \"name\": \"t3\" } }, { \"id\": \"t4\", \"input\": [ \"Infected\", \"Susceptible\" ], \"output\": [ \"Infected\", \"Infected\" ], \"properties\": { \"name\": \"t4\" } }, { \"id\": \"t5\", \"input\": [ \"Infected\" ], \"output\": [ \"Diagnosed\" ], \"properties\": { \"name\": \"t5\" } }, { \"id\": \"t6\", \"input\": [ \"Infected\" ], \"output\": [ \"Ailing\" ], \"properties\": { \"name\": \"t6\" } }, { \"id\": \"t7\", \"input\": [ \"Infected\" ], \"output\": [ \"Healed\" ], \"properties\": { \"name\": \"t7\" } }, { \"id\": \"t8\", \"input\": [ \"Diagnosed\" ], \"output\": [ \"Recognized\" ], \"properties\": { \"name\": \"t8\" } }, { \"id\": \"t9\", \"input\": [ \"Diagnosed\" ], \"output\": [ \"Healed\" ], \"properties\": { \"name\": \"t9\" } }, { \"id\": \"t10\", \"input\": [ \"Ailing\" ], \"output\": [ \"Recognized\" ], \"properties\": { \"name\": \"t10\" } }, { \"id\": \"t11\", \"input\": [ \"Ailing\" ], \"output\": [ \"Healed\" ], \"properties\": { \"name\": \"t11\" } }, { \"id\": \"t12\", \"input\": [ \"Ailing\" ], \"output\": [ \"Threatened\" ], \"properties\": { \"name\": \"t12\" } }, { \"id\": \"t13\", \"input\": [ \"Recognized\" ], \"output\": [ \"Threatened\" ], \"properties\": { \"name\": \"t13\" } }, { \"id\": \"t14\", \"input\": [ \"Recognized\" ], \"output\": [ \"Healed\" ], \"properties\": { \"name\": \"t14\" } }, { \"id\": \"t15\", \"input\": [ \"Threatened\" ], \"output\": [ \"Extinct\" ], \"properties\": { \"name\": \"t15\" } }, { \"id\": \"t16\", \"input\": [ \"Threatened\" ], \"output\": [ \"Healed\" ], \"properties\": { \"name\": \"t16\" } } ] }, \"semantics\": { \"ode\": { \"rates\": [ { \"target\": \"t1\", \"expression\": \"Diagnosed*Susceptible*beta\", \"expression_mathml\": \"DiagnosedSusceptiblebeta\" }, { \"target\": \"t2\", \"expression\": \"Ailing*Susceptible*gamma\", \"expression_mathml\": \"AilingSusceptiblegamma\" }, { \"target\": \"t3\", \"expression\": \"Recognized*Susceptible*delta\", \"expression_mathml\": \"RecognizedSusceptibledelta\" }, { \"target\": \"t4\", \"expression\": \"Infected*Susceptible*alpha\", \"expression_mathml\": \"InfectedSusceptiblealpha\" }, { \"target\": \"t5\", \"expression\": \"Infected*epsilon\", \"expression_mathml\": \"Infectedepsilon\" }, { \"target\": \"t6\", \"expression\": \"Infected*zeta\", \"expression_mathml\": \"Infectedzeta\" }, { \"target\": \"t7\", \"expression\": \"Infected*lambda\", \"expression_mathml\": \"Infectedlambda\" }, { \"target\": \"t8\", \"expression\": \"Diagnosed*eta\", \"expression_mathml\": \"Diagnosedeta\" }, { \"target\": \"t9\", \"expression\": \"Diagnosed*rho\", \"expression_mathml\": \"Diagnosedrho\" }, { \"target\": \"t10\", \"expression\": \"Ailing*theta\", \"expression_mathml\": \"Ailingtheta\" }, { \"target\": \"t11\", \"expression\": \"Ailing*kappa\", \"expression_mathml\": \"Ailingkappa\" }, { \"target\": \"t12\", \"expression\": \"Ailing*mu\", \"expression_mathml\": \"Ailingmu\" }, { \"target\": \"t13\", \"expression\": \"Recognized*nu\", \"expression_mathml\": \"Recognizednu\" }, { \"target\": \"t14\", \"expression\": \"Recognized*xi\", \"expression_mathml\": \"Recognizedxi\" }, { \"target\": \"t15\", \"expression\": \"Threatened*tau\", \"expression_mathml\": \"Threatenedtau\" }, { \"target\": \"t16\", \"expression\": \"Threatened*sigma\", \"expression_mathml\": \"Threatenedsigma\" } ], \"initials\": [ { \"target\": \"Susceptible\", \"expression\": \"0.999996300000000\", \"expression_mathml\": \"0.99999629999999995\" }, { \"target\": \"Diagnosed\", \"expression\": \"3.33333333000000e-7\", \"expression_mathml\": \"3.33333333e-7\" }, { \"target\": \"Infected\", \"expression\": \"3.33333333000000e-6\", \"expression_mathml\": \"3.3333333299999999e-6\" }, { \"target\": \"Ailing\", \"expression\": \"1.66666666000000e-8\", \"expression_mathml\": \"1.6666666599999999e-8\" }, { \"target\": \"Recognized\", \"expression\": \"3.33333333000000e-8\", \"expression_mathml\": \"3.33333333e-8\" }, { \"target\": \"Healed\", \"expression\": \"0.0\", \"expression_mathml\": \"0.0\" }, { \"target\": \"Threatened\", \"expression\": \"0.0\", \"expression_mathml\": \"0.0\" }, { \"target\": \"Extinct\", \"expression\": \"0.0\", \"expression_mathml\": \"0.0\" } ], \"parameters\": [ { \"id\": \"beta\", \"value\": 0.011, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.008799999999999999, \"maximum\": 0.0132 } } }, { \"id\": \"gamma\", \"value\": 0.456, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.3648, \"maximum\": 0.5472 } } }, { \"id\": \"delta\", \"value\": 0.011, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.008799999999999999, \"maximum\": 0.0132 } } }, { \"id\": \"alpha\", \"value\": 0.57, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.45599999999999996, \"maximum\": 0.6839999999999999 } } }, { \"id\": \"epsilon\", \"value\": 0.171, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.1368, \"maximum\": 0.20520000000000002 } } }, { \"id\": \"zeta\", \"value\": 0.125, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.1, \"maximum\": 0.15 } } }, { \"id\": \"lambda\", \"value\": 0.034, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.027200000000000002, \"maximum\": 0.0408 } } }, { \"id\": \"eta\", \"value\": 0.125, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.1, \"maximum\": 0.15 } } }, { \"id\": \"rho\", \"value\": 0.034, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.027200000000000002, \"maximum\": 0.0408 } } }, { \"id\": \"theta\", \"value\": 0.371, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.2968, \"maximum\": 0.4452 } } }, { \"id\": \"kappa\", \"value\": 0.017, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.013600000000000001, \"maximum\": 0.0204 } } }, { \"id\": \"mu\", \"value\": 0.017, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.013600000000000001, \"maximum\": 0.0204 } } }, { \"id\": \"nu\", \"value\": 0.027, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.0216, \"maximum\": 0.0324 } } }, { \"id\": \"xi\", \"value\": 0.017, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.013600000000000001, \"maximum\": 0.0204 } } }, { \"id\": \"tau\", \"value\": 0.01, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.008, \"maximum\": 0.012 } } }, { \"id\": \"sigma\", \"value\": 0.017, \"distribution\": { \"type\": \"StandardUniform1\", \"parameters\": { \"minimum\": 0.013600000000000001, \"maximum\": 0.0204 } } } ], \"observables\": [ { \"id\": \"Cases\", \"name\": \"Cases\", \"expression\": \"Diagnosed + Recognized + Threatened\", \"expression_mathml\": \"DiagnosedRecognizedThreatened\" }, { \"id\": \"Hospitalizations\", \"name\": \"Hospitalizations\", \"expression\": \"Recognized + Threatened\", \"expression_mathml\": \"RecognizedThreatened\" }, { \"id\": \"Deaths\", \"name\": \"Deaths\", \"expression\": \"Extinct\", \"expression_mathml\": \"Extinct\" } ], \"time\": { \"id\": \"t\", \"units\": { \"expression\": \"day\", \"expression_mathml\": \"day\" } } } }, \"metadata\": { \"annotations\": { \"license\": \"CC0\", \"authors\": [], \"references\": [ \"pubmed:32322102\" ], \"time_scale\": null, \"time_start\": null, \"time_end\": null, \"locations\": [], \"pathogens\": [ \"ncbitaxon:2697049\" ], \"diseases\": [ \"doid:0080600\" ], \"hosts\": [ \"ncbitaxon:9606\" ], \"model_types\": [ \"mamo:0000028\" ] } }}",
- "dataset": {
- "id": "2ea2d39f-866f-46f6-beec-972ed2136ed5",
- "filename": "dataset.csv"
- },
- "timespan": {"start": 101, "end": 190}
+ "engine": "sciml",
+ "model_config_id": "c1cd941a-047d-11ee-be56",
+ "dataset": {
+ "id": "cd339570-047d-11ee-be55",
+ "filename": "dataset.csv",
+ "mappings": {
+ "postive_tests": "infected"
+ }
+ },
+ "timespan": {
+ "start": 0,
+ "end": 90
+ },
+ "extra": {
+ "num_chains": 4,
+ "num_iterations": 100,
+ "ode_method": "default",
+ "calibrate_method": "bayesian"
}
+}
diff --git a/src/SimulationService.jl b/src/SimulationService.jl
index c5924b0..eb4d30a 100644
--- a/src/SimulationService.jl
+++ b/src/SimulationService.jl
@@ -25,6 +25,7 @@ import SwaggerMarkdown
import SymbolicUtils
import UUIDs
import YAML
+import Statistics
export start!, stop!
diff --git a/src/operations.jl b/src/operations.jl
index 445773f..700e9c9 100644
--- a/src/operations.jl
+++ b/src/operations.jl
@@ -58,41 +58,80 @@ struct Calibrate <: Operation
sys::ODESystem
timespan::Tuple{Float64, Float64}
priors::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Uniform{Float64}}}
- data::Any # ???
+ data::Any
+ num_chains::Int
+ num_iterations::Int
+ calibrate_method::String
+ ode_method::Any
end
function Calibrate(o::OperationRequest)
sys = amr_get(o.model, ODESystem)
priors = amr_get(o.model, sys, Val(:priors))
data = amr_get(o.df, sys, Val(:data))
- Calibrate(sys, o.timespan, priors, data)
+
+ num_chains = 4
+ num_iterations = 100
+ calibrate_method = "bayesian"
+ ode_method = nothing
+
+ if :extra in keys(o)
+ extrakeys = keys(o.extra)
+ :num_chains in extrakeys && (num_chains = o.extra.num_chains)
+ :num_iterations in extrakeys && (num_iterations = o.extra.num_iterations)
+ :calibrate_method in extrakeys && (calibrate_method = o.extra.calibrate_method)
+ end
+ Calibrate(sys, o.timespan, priors, data, num_chains, num_iterations, calibrate_method, ode_method)
end
function solve(o::Calibrate; callback)
prob = ODEProblem(o.sys, [], o.timespan)
statenames = [states(o.sys);getproperty.(observed(o.sys), :lhs)]
- p_posterior = EasyModelAnalysis.bayesian_datafit(prob, o.priors, o.data;
- nchains = 2,
- niter = 100,
- mcmcensemble = SimulationService.EasyModelAnalysis.Turing.MCMCSerial())
-
- pvalues = last.(p_posterior)
-
- probs = [EasyModelAnalysis.remake(prob, p = Pair.(first.(p_posterior), getindex.(pvalues,i))) for i in 1:length(p_posterior[1][2])]
- enprob = EasyModelAnalysis.EnsembleProblem(probs)
- ensol = solve(enprob, saveat = 1)
- outs = map(1:length(probs)) do i
- mats = stack(ensol[i][statenames])'
- headers = string.("ensemble",i,"_", statenames)
- mats, headers
- end
- dfsim = DataFrame(hcat(ensol[1].t, reduce(hcat, first.(outs))), :auto)
- rename!(dfsim, ["timestamp";reduce(vcat, last.(outs))])
-
- dfparam = DataFrame(last.(p_posterior), :auto)
- rename!(dfparam, Symbol.(first.(p_posterior)))
- dfsim, dfparam
+ if o.calibrate_method == "bayesian"
+ p_posterior = EasyModelAnalysis.bayesian_datafit(prob, o.priors, o.data;
+ nchains = 2,
+ niter = 100,
+ mcmcensemble = SimulationService.EasyModelAnalysis.Turing.MCMCSerial())
+
+ pvalues = last.(p_posterior)
+
+ probs = [EasyModelAnalysis.remake(prob, p = Pair.(first.(p_posterior), getindex.(pvalues,i))) for i in 1:length(p_posterior[1][2])]
+ enprob = EasyModelAnalysis.EnsembleProblem(probs)
+ ensol = solve(enprob, saveat = 1)
+ outs = map(1:length(probs)) do i
+ mats = stack(ensol[i][statenames])'
+ headers = string.("ensemble",i,"_", statenames)
+ mats, headers
+ end
+ dfsim = DataFrame(hcat(ensol[1].t, reduce(hcat, first.(outs))), :auto)
+ rename!(dfsim, ["timestamp";reduce(vcat, last.(outs))])
+
+ dfparam = DataFrame(last.(p_posterior), :auto)
+ rename!(dfparam, Symbol.(first.(p_posterior)))
+
+ dfsim, dfparam
+ elseif o.calibrate_method == "local" || o.calibrate_method == "global"
+ if o.calibrate_method == "local"
+ init_params = Pair.(EasyModelAnalysis.ModelingToolkit.Num.(first.(o.priors)), Statistics.mean.(last.(o.priors)))
+ fit = EasyModelAnalysis.datafit(prob, init_params, o.data)
+ else
+ init_params = Pair.(EasyModelAnalysis.ModelingToolkit.Num.(first.(o.priors)), tuple.(minimum.(last.(o.priors)), maximum.(last.(o.priors))))
+ fit = global_datafit(prob, init_params, o.data)
+ end
+
+ newprob = remake(prob, p=fit)
+ sol = solve(newprob)
+ dfsim = DataFrame(hcat(sol.t,stack(sol[statenames])'), :auto)
+ rename!(dfsim, ["timestamp";string.(statenames)])
+
+ dfparam = DataFrame(Matrix(last.(fit)'), :auto)
+ rename!(dfparam, Symbol.(first.(fit)))
+
+ dfsim, dfparam
+ else
+ error("$(o.calibrate_method) is not a valid choice of calibration method")
+ end
end
#-----------------------------------------------------------------------------# Ensemble
diff --git a/test/runtests.jl b/test/runtests.jl
index 374f49f..42f74da 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -84,7 +84,13 @@ end
priors = SimulationService.amr_get(amr, sys, Val(:priors))
df = CSV.read(here("examples", "dataset.csv"), DataFrame)
data = SimulationService.amr_get(df, sys, Val(:data))
- o = SimulationService.Calibrate(sys, (0.0, 89.0), priors, data)
+ num_chains = 4
+ num_iterations = 100
+ calibrate_method = "bayesian"
+ ode_method = nothing
+ o = SimulationService.Calibrate(sys, (0.0, 89.0), priors, data, num_chains, num_iterations, calibrate_method, ode_method)
+
+
dfsim, dfparam = SimulationService.solve(o; callback = nothing)
statenames = [states(o.sys);getproperty.(observed(o.sys), :lhs)]
@@ -105,7 +111,19 @@ end
end
@testset "calibrate" begin
- # TODO
+ json_url = here("examples", "request-calibrate-no-integration.json")
+ obj = JSON3.read(read(json_url), Config)
+
+ amr_url = here("examples", "BIOMD0000000955_askenet.json")
+ amr = JSON3.read(read(amr_url), Config)
+ priors = SimulationService.amr_get(amr, sys, Val(:priors))
+ df = CSV.read(here("examples", "dataset.csv"), DataFrame)
+
+ sys = SimulationService.amr_get(obj, ODESystem)
+ op = Simulate(sys, (0.0, 99.0))
+ df = solve(op)
+ @test df isa DataFrame
+ @test extrema(df.timestamp) == (0.0, 99.0)
end
@testset "ensemble" begin