diff --git a/.gitignore b/.gitignore index e069cf40..e6e759e1 100644 --- a/.gitignore +++ b/.gitignore @@ -31,8 +31,8 @@ Manifest.toml *.code-workspace .vscode benchmarks/matlab*/src -*.txt -!benchmarks/matlab*/iqm/*/*/*/*.txt +#*.txt +#!benchmarks/matlab*/iqm/*/*/*/*.txt timecommand objf*.m sharable diff --git a/ParamterEstimation/Project.toml b/ParamterEstimation/Project.toml new file mode 100644 index 00000000..e158c71e --- /dev/null +++ b/ParamterEstimation/Project.toml @@ -0,0 +1,4 @@ +[deps] +BaryRational = "91aaffc3-5777-4842-85b7-5d3d5d6a3494" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" diff --git a/Project.toml b/Project.toml index 5628cad2..8a847bce 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,13 @@ name = "ParameterEstimation" uuid = "b4cd1eb8-1e24-11e8-3319-93036a3eb9f3" -authors = ["Soo Go ", "Hoon Hong ", "Ilia Ilmer ", "Alexey Ovchinnikov ", "Pedro Soto ", "Chee Yap "] -version = "0.2.1" +authors = ["Oren Bassik ", "Soo Go ", "Hoon Hong ", "Ilia Ilmer ", "Alexey Ovchinnikov ", "Pedro Soto ", "Chee Yap "] +version = "0.3.0" [deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +ArbNumerics = "7e558dbc-694d-5a72-987c-6f4ebed21442" +BaryRational = "91aaffc3-5777-4842-85b7-5d3d5d6a3494" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -14,28 +16,29 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" SIAN = "cf7bdac0-b945-4905-b5ad-bc3f1a757483" StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c" TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04" [compat] -BenchmarkTools = "1" DifferentialEquations = "7" -Groebner = "0.2, 0.3, 0.4, 0.5" +Groebner = "0.2, 0.3, 0.4 - 0.5" HomotopyContinuation = "2" -LinearSolve = "1, 2" +LinearSolve = "2" ModelingToolkit = "8" OrderedCollections = "1" -Oscar = "0.10, 0.11, 0.12, 0.13" +Oscar = "0.10 - 0.20" ProgressMeter = "1" -SIAN = "1.4.1" -StructuralIdentifiability = "0.4.9, 0.5" +SIAN = "1 - 1.5.10" +StructuralIdentifiability = "0.4.9 - 0.5.30" Symbolics = "4, 5" -TaylorSeries = "0.12, 0.13, 0.14, 0.15" +TaylorSeries = "0.12 - 0.30" TestSetExtensions = "2" -julia = "1" diff --git a/README.md b/README.md index 85ed4c68..906e17fc 100644 --- a/README.md +++ b/README.md @@ -17,17 +17,19 @@ Currently is installable via ```julia using Pkg -Pkg.add("ParameterEstimation.jl") +Pkg.add(url="https://github.com/orebas/ParameterEstimation.jl") ``` -or +The production version of this fork is installable via ```julia using Pkg -Pkg.add(url="https://github.com/iliailmer/ParameterEstimation.jl") +Pkg.add("ParameterEstimation.jl") ``` + + ## Toy Example ```julia diff --git a/benchmarks/all_models.jl b/benchmarks/all_models.jl new file mode 100644 index 00000000..5f723de2 --- /dev/null +++ b/benchmarks/all_models.jl @@ -0,0 +1,493 @@ +using ParameterEstimation +using ModelingToolkit, DifferentialEquations#, Plots +using BenchmarkTools + + +struct ParameterEstimationProblem + Name::Any + model::Any + measured_quantities::Any + data_sample::Any + solver::Any + p_true::Any + ic::Any +end + + +function biohydrogenation(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters k5 k6 k7 k8 k9 k10 + @variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t) + D = Differential(t) + states = [x4, x5, x6, x7] + parameters = [k5, k6, k7, k8, k9, k10] + + @named model = ODESystem([ + D(x4) ~ -k5 * x4 / (k6 + x4), + D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6), + D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10, + D(x7) ~ k9 * x6 * (k10 - x6) / k10, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x4, + y2 ~ x5, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("BioHydrogenation", model, measured_quantities, data_sample, solver, p_true, ic) +end + +function crauste(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters mu_N mu_EE mu_LE mu_LL mu_M mu_P mu_PE mu_PL delta_NE delta_EL delta_LM rho_E rho_P + @variables t N(t) E(t) S(t) M(t) P(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [N, E, S, M, P] + parameters = [ + mu_N, + mu_EE, + mu_LE, + mu_LL, + mu_M, + mu_P, + mu_PE, + mu_PL, + delta_NE, + delta_EL, + delta_LM, + rho_E, + rho_P, + ] + @named model = ODESystem( + [ + D(N) ~ -N * mu_N - N * P * delta_NE, + D(E) ~ N * P * delta_NE - E^2 * mu_EE - + E * delta_EL + E * P * rho_E, + D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL - + E * S * mu_LE, + D(M) ~ S * delta_LM - mu_M * M, + D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE - + S * P * mu_PL, + ], t, states, parameters) + measured_quantities = [y1 ~ N, y2 ~ E, y3 ~ S + M, y4 ~ P] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.071, 0.143, 0.214, 0.286, 0.357, 0.429, 0.5, 0.571, 0.643, 0.714, 0.786, 0.857, 0.929] # True Parameters + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Crauste", model, measured_quantities, data_sample, solver, p_true, ic) +end + + +function daisy_ex3(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters p1 p3 p4 p6 p7 + @variables t x1(t) x2(t) x3(t) u0(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, u0] + parameters = [p1, p3, p4, p6, p7] + @named model = ODESystem([ + D(x1) ~ -1 * p1 * x1 + x2 + u0, + D(x2) ~ p3 * x1 - p4 * x2 + x3, + D(x3) ~ p6 * x1 - p7 * x3, + D(u0) ~ 1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ u0, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_ex3", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + + +function daisy_mamil3(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters a12 a13 a21 a31 a01 + @variables t x1(t) x2(t) x3(t) y1(t) y2(t) + D = Differential(t) + + ic = [0.25, 0.5, 0.75] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + states = [x1, x2, x3] + parameters = [a12, a13, a21, a31, a01] + @named model = ODESystem([D(x1) ~ -(a21 + a31 + a01) * x1 + a12 * x2 + a13 * x3, + D(x2) ~ a21 * x1 - a12 * x2, + D(x3) ~ a31 * x1 - a13 * x3], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil3", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + +function daisy_mamil4(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters k01, k12, k13, k14, k21, k31, k41 + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) + D = Differential(t) + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] # True Parameters + + states = [x1, x2, x3, x4] + parameters = [k01, k12, k13, k14, k21, k31, k41] + @named model = ODESystem([D(x1) ~ -k01 * x1 + k12 * x2 + k13 * x3 + k14 * x4 - k21 * x1 - k31 * x1 - k41 * x1, + D(x2) ~ -k12 * x2 + k21 * x1, + D(x3) ~ -k13 * x3 + k31 * x1, + D(x4) ~ -k14 * x4 + k41 * x1], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2, y3 ~ x3 + x4] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil4", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + +function fitzhugh_nagumo(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters g a b + @variables t V(t) R(t) y1(t) y2(t) + D = Differential(t) + states = [V, R] + parameters = [g, a, b] + + ic = [0.333, 0.67] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ V] + + @named model = ODESystem([ + D(V) ~ g * (V - V^3 / 3 + R), + D(R) ~ 1 / g * (V - a + b * R), + ], t, states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("fitzhugh-nagumo", model, measured_quantities, data_sample, solver, p_true, ic) +end + + + +function hiv_local(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters b c d k1 k2 mu1 mu2 q1 q2 s + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, x4] + parameters = [b, c, d, k1, k2, mu1, mu2, q1, q2, s] + + @named model = ODESystem([ + D(x1) ~ -b * x1 * x4 - d * x1 + s, + D(x2) ~ b * q1 * x1 * x4 - k1 * x2 - mu1 * x2, + D(x3) ~ b * q2 * x1 * x4 + k1 * x2 - mu2 * x3, + D(x4) ~ -c * x4 + k2 * x3, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x4, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.091, 0.182, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + time_interval = [-0.5, 0.5] + datasize = 20 + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("hiv_local", model, measured_quantities, data_sample, solver, p_true, ic) +end + + +function hiv(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters lm d beta a k u c q b h + @variables t x(t) y(t) v(t) w(t) z(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [x, y, v, w, z] + parameters = [lm, d, beta, a, k, u, c, q, b, h] + + @named model = ODESystem([ + D(x) ~ lm - d * x - beta * x * v, + D(y) ~ beta * x * v - a * y, + D(v) ~ k * y - u * v, + D(w) ~ c * x * y * w - c * q * y * w - b * w, + D(z) ~ c * q * y * w - h * z, + ], t, states, parameters) + measured_quantities = [y1 ~ w, y2 ~ z, y3 ~ x, y4 ~ y + v] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.091, 0.181, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + + res = ParameterEstimation.estimate(model, measured_quantities, data_sample; + solver = solver) + return ParameterEstimationProblem("hiv", model, measured_quantities, data_sample, solver, p_true, ic) +end + + + +function lotka_volterra(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters k1 k2 k3 + @variables t r(t) w(t) y1(t) + D = Differential(t) + ic = [0.333, 0.667] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ r] + states = [r, w] + parameters = [k1, k2, k3] + + @named model = ODESystem([D(r) ~ k1 * r - k2 * r * w, D(w) ~ k2 * r * w - k3 * w], t, + states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Lotka_Volterra", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + +function seir(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters a b nu + @variables t S(t) E(t) In(t) N(t) y1(t) y2(t) + D = Differential(t) + states = [S, E, In, N] + parameters = [a, b, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N, + D(E) ~ b * S * In / N - nu * E, + D(In) ~ nu * E - a * In, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ In, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.25, 0.5, 0.75] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("SEIR", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + + + +function simple(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ -a * x2, + D(x2) ~ 1 / b * (x1), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("simple", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + + +function sirsforced(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) + + @parameters b0 b1 g M mu nu + @variables t i(t) r(t) s(t) x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [i, r, s, x1, x2] + parameters = [b0, b1, g, M, mu, nu] + + @named model = ODESystem([ + D(s) ~ mu - mu * s - b0 * (1 + b1 * x1) * i * s + g * r, + D(i) ~ b0 * (1 + b1 * x1) * i * s - (nu + mu) * i, + D(r) ~ nu * i - (mu + g) * r, + D(x1) ~ -M * x2, + D(x2) ~ M * x1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ i, + y2 ~ r, + ] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("sirsforced", model, measured_quantities, data_sample, solver, p_true, ic) + +end + +function slowfast(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) # TODO(orebas):in the old code it was CVODE_BDF. should we go back to that? + #solver = CVODE_BDF() + @parameters k1 k2 eB + @variables t xA(t) xB(t) xC(t) eA(t) eC(t) y1(t) y2(t) y3(t) y4(t) #eA(t) eC(t) + D = Differential(t) + states = [xA, xB, xC, eA, eC] + parameters = [k1, k2, eB] + @named model = ODESystem([ + D(xA) ~ -k1 * xA, + D(xB) ~ k1 * xA - k2 * xB, + D(xC) ~ k2 * xB, + D(eA) ~ 0, + D(eC) ~ 0, + ], t, states, parameters) + + measured_quantities = [y1 ~ xC, y2 ~ eA * xA + eB * xB + eC * xC, y3 ~ eA, y4 ~ eC] + ic = [0.166, 0.333, 0.5, 0.666, 0.833] + p_true = [0.25, 0.5, 0.75] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("slowfast", model, measured_quantities, data_sample, solver, p_true, ic) + +end + +function treatment(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) #note the solver. Tsit5 apparently can't handle mass matrices + @parameters a b d g nu + @variables t In(t) N(t) S(t) Tr(t) y1(t) y2(t) + D = Differential(t) + states = [In, N, S, Tr] + parameters = [a, b, d, g, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N - d * b * S * Tr / N, + D(In) ~ b * S * In / N + d * b * S * Tr / N - (a + g) * In, + D(Tr) ~ g * In - nu * Tr, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ Tr, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize, solver = solver) + return ParameterEstimationProblem("treatment", model, measured_quantities, data_sample, solver, p_true, ic) +end + + +function vanderpol(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ a * x2, + D(x2) ~ -(x1) - b * (x1^2 - 1) * (x2), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("vanderpol", model, measured_quantities, data_sample, solver, p_true, ic) +end + + + + + + + + + +function analyze_parameter_estimation_problem(PEP::ParameterEstimationProblem) + @time res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, PEP.data_sample; + solver = PEP.solver) + all_params = vcat(PEP.p_true) + for each in res + estimates = vcat(collect(values(each.parameters))) + println("For model ", PEP.Name, ": Max abs rel. err: ", + maximum(abs.((estimates .- all_params) ./ (all_params)))) + end +end + + +function analyze_parameter_estimation_problem_old(PEP::ParameterEstimationProblem) + res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, PEP.data_sample; + solver = PEP.solver) + all_params = vcat(PEP.ic, PEP.p_true) + for each in res + estimates = vcat(collect(values(each.states)), collect(values(each.parameters))) + println("Max abs rel. err: ", + maximum(abs.((estimates .- all_params) ./ (all_params)))) + end +end + +function main() + for PEP in [ + simple(), + lotka_volterra(), + vanderpol(), + biohydrogenation(), + daisy_ex3(), + daisy_mamil3(), + daisy_mamil4(), + fitzhugh_nagumo(), + hiv_local(), + hiv(), + seir(), + sirsforced(), + slowfast(), + treatment(), + crauste(), + ] + + analyze_parameter_estimation_problem(PEP) + + end +end + +main() diff --git a/benchmarks/all_models_2.jl b/benchmarks/all_models_2.jl new file mode 100644 index 00000000..886365f7 --- /dev/null +++ b/benchmarks/all_models_2.jl @@ -0,0 +1,590 @@ +using ParameterEstimation +using ModelingToolkit, DifferentialEquations#, Plots +using BenchmarkTools + +struct ParameterEstimationProblem + Name::Any + model::Any + measured_quantities::Any + data_sample::Any + solver::Any + p_true::Any + ic::Any +end + +function biohydrogenation(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters k5 k6 k7 k8 k9 k10 + @variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t) + D = Differential(t) + states = [x4, x5, x6, x7] + parameters = [k5, k6, k7, k8, k9, k10] + + @named model = ODESystem([ + D(x4) ~ -k5 * x4 / (k6 + x4), + D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6), + D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10, + D(x7) ~ k9 * x6 * (k10 - x6) / k10, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x4, + y2 ~ x5, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("BioHydrogenation", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function crauste(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters mu_N mu_EE mu_LE mu_LL mu_M mu_P mu_PE mu_PL delta_NE delta_EL delta_LM rho_E rho_P + @variables t N(t) E(t) S(t) M(t) P(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [N, E, S, M, P] + parameters = [ + mu_N, + mu_EE, + mu_LE, + mu_LL, + mu_M, + mu_P, + mu_PE, + mu_PL, + delta_NE, + delta_EL, + delta_LM, + rho_E, + rho_P, + ] + @named model = ODESystem([ + D(N) ~ -N * mu_N - N * P * delta_NE, + D(E) ~ N * P * delta_NE - E^2 * mu_EE - + E * delta_EL + E * P * rho_E, + D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL - + E * S * mu_LE, + D(M) ~ S * delta_LM - mu_M * M, + D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE - + S * P * mu_PL, + ], t, states, parameters) + measured_quantities = [y1 ~ N, y2 ~ E, y3 ~ S + M, y4 ~ P] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [ + 0.071, + 0.143, + 0.214, + 0.286, + 0.357, + 0.429, + 0.5, + 0.571, + 0.643, + 0.714, + 0.786, + 0.857, + 0.929, + ] # True Parameters + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Crauste", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function daisy_ex3(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters p1 p3 p4 p6 p7 + @variables t x1(t) x2(t) x3(t) u0(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, u0] + parameters = [p1, p3, p4, p6, p7] + @named model = ODESystem([ + D(x1) ~ -1 * p1 * x1 + x2 + u0, + D(x2) ~ p3 * x1 - p4 * x2 + x3, + D(x3) ~ p6 * x1 - p7 * x3, + D(u0) ~ 1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ u0, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_ex3", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function daisy_mamil3(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a12 a13 a21 a31 a01 + @variables t x1(t) x2(t) x3(t) y1(t) y2(t) + D = Differential(t) + + ic = [0.25, 0.5, 0.75] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + states = [x1, x2, x3] + parameters = [a12, a13, a21, a31, a01] + @named model = ODESystem([D(x1) ~ -(a21 + a31 + a01) * x1 + a12 * x2 + a13 * x3, + D(x2) ~ a21 * x1 - a12 * x2, + D(x3) ~ a31 * x1 - a13 * x3], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil3", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function daisy_mamil4(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters k01, k12, k13, k14, k21, k31, k41 + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) + D = Differential(t) + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] # True Parameters + + states = [x1, x2, x3, x4] + parameters = [k01, k12, k13, k14, k21, k31, k41] + @named model = ODESystem([ + D(x1) ~ -k01 * x1 + k12 * x2 + k13 * x3 + k14 * x4 - k21 * x1 - k31 * x1 - + k41 * x1, + D(x2) ~ -k12 * x2 + k21 * x1, + D(x3) ~ -k13 * x3 + k31 * x1, + D(x4) ~ -k14 * x4 + k41 * x1], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2, y3 ~ x3 + x4] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil4", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function fitzhugh_nagumo(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters g a b + @variables t V(t) R(t) y1(t) y2(t) + D = Differential(t) + states = [V, R] + parameters = [g, a, b] + + ic = [0.333, 0.67] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ V] + + @named model = ODESystem([ + D(V) ~ g * (V - V^3 / 3 + R), + D(R) ~ 1 / g * (V - a + b * R), + ], t, states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("fitzhugh-nagumo", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function hiv_local(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters b c d k1 k2 mu1 mu2 q1 q2 s + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, x4] + parameters = [b, c, d, k1, k2, mu1, mu2, q1, q2, s] + + @named model = ODESystem([ + D(x1) ~ -b * x1 * x4 - d * x1 + s, + D(x2) ~ b * q1 * x1 * x4 - k1 * x2 - mu1 * x2, + D(x3) ~ b * q2 * x1 * x4 + k1 * x2 - mu2 * x3, + D(x4) ~ -c * x4 + k2 * x3, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x4, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.091, 0.182, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + time_interval = [-0.5, 0.5] + datasize = 20 + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("hiv_local", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function hiv(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters lm d beta a k u c q b h + @variables t x(t) y(t) v(t) w(t) z(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [x, y, v, w, z] + parameters = [lm, d, beta, a, k, u, c, q, b, h] + + @named model = ODESystem([ + D(x) ~ lm - d * x - beta * x * v, + D(y) ~ beta * x * v - a * y, + D(v) ~ k * y - u * v, + D(w) ~ c * x * y * w - c * q * y * w - b * w, + D(z) ~ c * q * y * w - h * z, + ], t, states, parameters) + measured_quantities = [y1 ~ w, y2 ~ z, y3 ~ x, y4 ~ y + v] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.091, 0.181, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + + res = ParameterEstimation.estimate(model, measured_quantities, data_sample; + solver = solver) + return ParameterEstimationProblem("hiv", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function lotka_volterra(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters k1 k2 k3 + @variables t r(t) w(t) y1(t) + D = Differential(t) + ic = [0.333, 0.667] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ r] + states = [r, w] + parameters = [k1, k2, k3] + + @named model = ODESystem([D(r) ~ k1 * r - k2 * r * w, D(w) ~ k2 * r * w - k3 * w], t, + states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Lotka_Volterra", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function seir(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a b nu + @variables t S(t) E(t) In(t) N(t) y1(t) y2(t) + D = Differential(t) + states = [S, E, In, N] + parameters = [a, b, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N, + D(E) ~ b * S * In / N - nu * E, + D(In) ~ nu * E - a * In, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ In, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.25, 0.5, 0.75] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("SEIR", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function simple(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ -a * x2, + D(x2) ~ 1 / b * (x1), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("simple", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function sirsforced(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) + @parameters b0 b1 g M mu nu + @variables t i(t) r(t) s(t) x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [i, r, s, x1, x2] + parameters = [b0, b1, g, M, mu, nu] + + @named model = ODESystem([ + D(s) ~ mu - mu * s - b0 * (1 + b1 * x1) * i * s + g * r, + D(i) ~ b0 * (1 + b1 * x1) * i * s - (nu + mu) * i, + D(r) ~ nu * i - (mu + g) * r, + D(x1) ~ -M * x2, + D(x2) ~ M * x1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ i, + y2 ~ r, + ] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("sirsforced", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function slowfast(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) # TODO(orebas):in the old code it was CVODE_BDF. should we go back to that? + #solver = CVODE_BDF() + @parameters k1 k2 eB + @variables t xA(t) xB(t) xC(t) eA(t) eC(t) y1(t) y2(t) y3(t) y4(t) #eA(t) eC(t) + D = Differential(t) + states = [xA, xB, xC, eA, eC] + parameters = [k1, k2, eB] + @named model = ODESystem([ + D(xA) ~ -k1 * xA, + D(xB) ~ k1 * xA - k2 * xB, + D(xC) ~ k2 * xB, + D(eA) ~ 0, + D(eC) ~ 0, + ], t, states, parameters) + + measured_quantities = [y1 ~ xC, y2 ~ eA * xA + eB * xB + eC * xC, y3 ~ eA, y4 ~ eC] + ic = [0.166, 0.333, 0.5, 0.666, 0.833] + p_true = [0.25, 0.5, 0.75] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, + measured_quantities, + time_interval, + p_true, + ic, + datasize; + solver = solver) + + return ParameterEstimationProblem("slowfast", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function treatment(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) #note the solver. Vern9 apparently can't handle mass matrices + @parameters a b d g nu + @variables t In(t) N(t) S(t) Tr(t) y1(t) y2(t) + D = Differential(t) + states = [In, N, S, Tr] + parameters = [a, b, d, g, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N - d * b * S * Tr / N, + D(In) ~ b * S * In / N + d * b * S * Tr / N - (a + g) * In, + D(Tr) ~ g * In - nu * Tr, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ Tr, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize, solver = solver) + return ParameterEstimationProblem("treatment", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function vanderpol(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ a * x2, + D(x2) ~ -(x1) - b * (x1^2 - 1) * (x2), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("vanderpol", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic, + test_mode = false) +end + +function analyze_parameter_estimation_problem(PEP::ParameterEstimationProblem) + + #interpolators = Dict( + #"AAA" => ParameterEstimation.aaad, + #"FHD3" => ParameterEstimation.fhdn(3), + #"FHD6" => ParameterEstimation.fhdn(6), + #"FHD8" => ParameterEstimation.fhdn(8), + #"Fourier" => ParameterEstimation.FourierInterp, + #) + datasize = 21 #TODO(Orebas) magic number + + #stepsize = max(1, datasize ÷ 8) + #for i in range(1, (datasize - 2), step = stepsize) + # interpolators["RatOld($i)"] = ParameterEstimation.SimpleRationalInterpOld(i) + #end + + @time res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, + PEP.data_sample, + solver = PEP.solver) # , interpolators) + all_params = vcat(PEP.ic, PEP.p_true) + #println("TYPERES: ", typeof(res)) + #println(res) + + #println(res) + besterror = 1e30 + for each in res + #println("TYPE: ", typeof(each)) + #println(each) + estimates = vcat(collect(values(each.states)), collect(values(each.parameters))) + besterror = min(besterror, maximum(abs.((estimates .- all_params) ./ (all_params)))) + end + if (test_mode) + @test besterror < 1e-2 + else + println("For model ", PEP.Name, ": The best max abs rel. err: ", besterror) + end +end + +#function analyze_parameter_estimation_problem_old(PEP::ParameterEstimationProblem) +# res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, PEP.data_sample; +# solver = PEP.solver, interpolator = ("AAA" => aaad)) +# all_params = vcat(PEP.ic, PEP.p_true) +# for each in res +# estimates = vcat(collect(values(each.states)), collect(values(each.parameters))) +# println("Max abs rel. err: ", +# maximum(abs.((estimates .- all_params) ./ (all_params)))) +# end +#end + +function main() + datasize = 21 + solver = Vern9() + #solver = Rodas4P() + time_interval = [-0.5, 0.5] + for PEP in [ + simple(datasize, time_interval, solver), #works + lotka_volterra(datasize, time_interval, solver), #works + vanderpol(datasize, time_interval, solver), #works + biohydrogenation(datasize, time_interval, solver), #works, but one param unidentifiable + #daisy_ex3(datasize, time_interval, solver), + daisy_mamil3(datasize, time_interval, solver), + daisy_mamil4(datasize, time_interval, solver), + #fitzhugh_nagumo(datasize, time_interval, solver), + hiv_local(datasize, time_interval, solver), + hiv(datasize, time_interval, solver), + seir(datasize, time_interval, solver), + sirsforced(datasize, time_interval, Rodas5P()), + slowfast(datasize, time_interval, solver), + treatment(datasize, time_interval, Rodas5P()), + crauste(datasize, time_interval, solver), + ] + analyze_parameter_estimation_problem(PEP) + end +end + +main() diff --git a/benchmarks/backwards_res.txt b/benchmarks/backwards_res.txt new file mode 100644 index 00000000..3b51da4f --- /dev/null +++ b/benchmarks/backwards_res.txt @@ -0,0 +1,488 @@ +WARNING: Method definition isapprox(IntervalSets.AbstractInterval{T} where T, IntervalSets.AbstractInterval{T} where T) in module IntervalSets at /home/orebas/.julia/packages/IntervalSets/6RTOk/src/IntervalSets.jl:144 overwritten in module DomainSets at /home/orebas/.julia/packages/DomainSets/aafhp/src/domains/interval.jl:52. + ** incremental compilation may be fatally broken for this module ** + + DEBUG +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [lm, d, beta, a, k, u, c, q, b, h, x, y, v, w, z] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [u, beta, k, lm, c, w, q, d, h, v, z, a, x, y, b] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] + Computing mixed cells... 2 Time: 0:00:00 + mixed_volume: 2  +  Computing mixed cells... 4 Time: 0:00:00 + mixed_volume: 4 +4 Froissart doublets. Number of residues = 10 +4 Froissart doublets. Number of residues = 10 +6 Froissart doublets. Number of residues = 10 +4 Froissart doublets. Number of residues = 10 +Parameter(s) : lm = 0.091, d = 0.181, beta = 0.273, a = 0.364, k = 0.454, u = 0.545, c = 0.635, q = 0.727, b = 0.818, h = 0.909 +Initial Condition(s): x(t) = 0.167, y(t) = 0.333, v(t) = 0.500, w(t) = 0.667, z(t) = 0.833 +Error: 3.4330e-03 +Time: -5.0000e-01 + +DEBUG +1 Froissart doublets. Number of residues = 10 +5 Froissart doublets. Number of residues = 10 +Parameter(s) : a = 0.333, b = 0.667 +Initial Condition(s): x1(t) = 0.333, x2(t) = 0.667 +Error: 4.3680e-03 +Time: -5.0000e-01 + + 12.353256 seconds (7.31 M allocations: 483.280 MiB, 1.68% gc time, 96.27% compilation time) +For model simple: Max abs rel. err: 8.216821193133987e-9 +DEBUG +3 Froissart doublets. Number of residues = 8 +Parameter(s) : k1 = 0.250, k2 = 0.500, k3 = 0.750 +Initial Condition(s): r(t) = 0.333, w(t) = 0.667 +Error: 5.4874e-04 +Time: -5.0000e-01 + + 6.883697 seconds (5.34 M allocations: 345.916 MiB, 2.97% gc time, 97.90% compilation time: <1% of which was recompilation) +For model Lotka_Volterra: Max abs rel. err: 2.4765554273997736e-6 +DEBUG +3 Froissart doublets. Number of residues = 10 +4 Froissart doublets. Number of residues = 10 +Parameter(s) : a = 0.333, b = 0.667 +Initial Condition(s): x1(t) = 0.333, x2(t) = 0.667 +Error: 3.4165e-03 +Time: -5.0000e-01 + + 5.944709 seconds (3.40 M allocations: 226.845 MiB, 98.44% compilation time) +For model vanderpol: Max abs rel. err: 6.098011933822437e-8 +DEBUG +6 Froissart doublets. Number of residues = 10 +6 Froissart doublets. Number of residues = 10 +Parameter(s) : k5 = 0.143, k6 = 0.286, k7 = 0.111, k8 = -4.472, k9 = 0.224, k10 = 4.850 +Initial Condition(s): x4(t) = 0.200, x5(t) = 0.400, x6(t) = 4.478, x7(t) = 3551765860466.962 +Error: 9.8602e-05 +Time: -5.0000e-01 + +Parameter(s) : k5 = 0.143, k6 = 0.286, k7 = 0.429, k8 = 1.426, k9 = -0.715, k10 = -0.854 +Initial Condition(s): x4(t) = 0.200, x5(t) = 0.400, x6(t) = -0.256, x7(t) = 3551765860466.935 +Error: 9.8602e-05 +Time: -5.0000e-01 + +Parameter(s) : k5 = 0.143, k6 = 0.286, k7 = 0.111, k8 = 0.378, k9 = -0.224, k10 = -4.850 +Initial Condition(s): x4(t) = 0.200, x5(t) = 0.400, x6(t) = -0.372, x7(t) = 3551765860466.962 +Error: 9.8602e-05 +Time: -5.0000e-01 + +Parameter(s) : k5 = 0.143, k6 = 0.286, k7 = 0.429, k8 = 0.572, k9 = 0.715, k10 = 0.854 +Initial Condition(s): x4(t) = 0.200, x5(t) = 0.400, x6(t) = 0.598, x7(t) = 3551765860466.935 +Error: 9.8602e-05 +Time: -5.0000e-01 + + 29.773833 seconds (21.96 M allocations: 1.400 GiB, 2.68% gc time, 83.04% compilation time) +For model BioHydrogenation: Max abs rel. err: 4.439707325582703e12 +For model BioHydrogenation: Max abs rel. err: 4.439707325582669e12 +For model BioHydrogenation: Max abs rel. err: 4.439707325582702e12 +For model BioHydrogenation: Max abs rel. err: 4.439707325582669e12 +DEBUG +1 Froissart doublets. Number of residues = 10 +3 Froissart doublets. Number of residues = 10 + 8.592710 seconds (8.65 M allocations: 574.370 MiB, 3.32% gc time, 87.48% compilation time) +DEBUG +3 Froissart doublets. Number of residues = 10 +5 Froissart doublets. Number of residues = 10 +Parameter(s) : a12 = 0.167, a13 = 0.346, a21 = 0.500, a31 = 0.687, a01 = 0.805 +Initial Condition(s): x1(t) = 0.250, x2(t) = 0.500, x3(t) = 0.717 +Error: 9.5935e-04 +Time: -5.0000e-01 + + 11.823468 seconds (8.61 M allocations: 570.164 MiB, 2.71% gc time, 96.89% compilation time) +For model DAISY_mamil3: Max abs rel. err: 0.04446665857000701 +DEBUG +3 Froissart doublets. Number of residues = 10 +3 Froissart doublets. Number of residues = 10 +3 Froissart doublets. Number of residues = 10 +Parameter(s) : k01 = 0.125, k12 = 0.250, k13 = 0.500, k14 = 0.375, k21 = 0.625, k31 = 0.874, k41 = 0.751 +Initial Condition(s): x1(t) = 0.200, x2(t) = 0.400, x3(t) = 0.800, x4(t) = 0.600 +Error: 1.6178e-03 +Time: -5.0000e-01 + +Parameter(s) : k01 = 0.125, k12 = 0.250, k13 = 0.375, k14 = 0.500, k21 = 0.625, k31 = 0.751, k41 = 0.874 +Initial Condition(s): x1(t) = 0.200, x2(t) = 0.400, x3(t) = 0.600, x4(t) = 0.800 +Error: 1.6178e-03 +Time: -5.0000e-01 + + 14.358593 seconds (10.84 M allocations: 722.351 MiB, 2.43% gc time, 62.31% compilation time) +For model DAISY_mamil4: Max abs rel. err: 0.3332119224172925 +For model DAISY_mamil4: Max abs rel. err: 0.0009308470261812177 +DEBUG +1 Froissart doublets. Number of residues = 10 +Parameter(s) : g = 0.250, a = 0.500, b = 0.750 +Initial Condition(s): V(t) = 0.333, R(t) = 0.670 +Error: 3.9937e-02 +Time: -5.0000e-01 + + 8.068035 seconds (5.87 M allocations: 389.285 MiB, 3.57% gc time, 97.82% compilation time) +For model fitzhugh-nagumo: Max abs rel. err: 1.4541186840371288e-5 +DEBUG +1 Froissart doublets. Number of residues = 10 +4 Froissart doublets. Number of residues = 10 + 36.762913 seconds (17.75 M allocations: 1.142 GiB, 1.47% gc time, 34.73% compilation time) +DEBUG +4 Froissart doublets. Number of residues = 10 +4 Froissart doublets. Number of residues = 10 +6 Froissart doublets. Number of residues = 10 +4 Froissart doublets. Number of residues = 10 +Parameter(s) : lm = 0.091, d = 0.181, beta = 0.273, a = 0.364, k = 0.454, u = 0.545, c = 0.635, q = 0.727, b = 0.818, h = 0.909 +Initial Condition(s): x(t) = 0.167, y(t) = 0.333, v(t) = 0.500, w(t) = 0.667, z(t) = 0.833 +Error: 3.4330e-03 +Time: -5.0000e-01 + + 3.797847 seconds (3.14 M allocations: 472.273 MiB, 5.87% gc time) +For model hiv: Max abs rel. err: 0.0023299982577594204 +DEBUG +4 Froissart doublets. Number of residues = 9 +Parameter(s) : a = 0.750, b = 0.500, nu = 0.250 +Initial Condition(s): S(t) = 0.600, E(t) = 2.400, In(t) = 0.600, N(t) = 0.800 +Error: 1.0245e-03 +Time: -5.0000e-01 + +Parameter(s) : a = 0.250, b = 0.500, nu = 0.750 +Initial Condition(s): S(t) = 0.200, E(t) = 0.400, In(t) = 0.600, N(t) = 0.800 +Error: 1.0245e-03 +Time: -5.0000e-01 + + 10.025146 seconds (9.89 M allocations: 615.467 MiB, 3.63% gc time, 91.50% compilation time) +For model SEIR: Max abs rel. err: 4.999790001393897 +For model SEIR: Max abs rel. err: 7.445569848635314e-5 +DEBUG +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [a, b, x2, x1] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [b, x2, x1, a] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [k1, k2, k3, w, r] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [k3, w, r, k1, k2] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [a, b, x2, x1] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [b, x2, x1, a] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [x7_0] +[ Info: Locally identifiable parameters: [k5, k6, k7, k8, k9, k10, x6, x5, x4] +[ Info: Not identifiable parameters: [x7] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [k6, k7, x5, k5, x4] +[ Info: Locally but not globally identifiable parameters: [k8, x6, k9, k10] +[ Info: Not identifiable parameters: [x7] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [p1, p3, p4, p6, p7, u0, x3, x2, x1] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [p1, p3, u0, x2, x1] +[ Info: Locally but not globally identifiable parameters: [p4, x3, p6, p7] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +┌ Warning: No solution found +└ @ ParameterEstimation ~/.julia/dev/ParameterEstimation/src/estimation/utils.jl:28 +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [a12, a13, a21, a31, a01, x2, x3, x1] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [a01, a13, x2, x3, a31, a12, x1, a21] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [k01, k12, k13, k14, k21, k31, k41, x3, x1, x2, x4] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [k12, k21, x1, k01, x2] +[ Info: Locally but not globally identifiable parameters: [x4, k13, x3, k41, k14, k31] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [g, a, b, R, V] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [b, R, V, g, a] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [k2_0, k1_0] +[ Info: Locally identifiable parameters: [b, c, d, mu2, s, x1, x4] +[ Info: Not identifiable parameters: [k1, k2, mu1, q1, q2, x3, x2] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [b, x1, d, s, x4] +[ Info: Locally but not globally identifiable parameters: [c, mu2] +[ Info: Not identifiable parameters: [q1, mu1, k2, x3, k1, q2, x2] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +┌ Warning: No solution found +└ @ ParameterEstimation ~/.julia/dev/ParameterEstimation/src/estimation/utils.jl:28 +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [lm, d, beta, a, k, u, c, q, b, h, x, y, v, w, z] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [u, beta, k, lm, c, w, q, d, h, v, z, a, x, y, b] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [a, b, nu, E, In, S, N] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [In, b, N] +[ Info: Locally but not globally identifiable parameters: [E, a, nu, S] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [b1_0] +[ Info: Locally identifiable parameters: [b0, g, M, mu, nu, i, s, r] +[ Info: Not identifiable parameters: [b1, x1, x2] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [i, nu, mu, b0, g, s, r] +[ Info: Locally but not globally identifiable parameters: [M] +[ Info: Not identifiable parameters: [x1, x2, b1] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] + Computing mixed cells... 21 Time: 0:00:00 + mixed_volume: 60  +  Computing mixed cells... 52 Time: 0:00:00 + mixed_volume: 118 +2 Froissart doublets. Number of residues = 10 +3 Froissart doublets. Number of residues = 10 +294.187442 seconds (19.28 M allocations: 1.046 GiB, 0.22% gc time, 3.63% compilation time) +DEBUG +1 Froissart doublets. Number of residues = 9 +4 Froissart doublets. Number of residues = 7 +Parameter(s) : k1 = 0.500, k2 = 0.250, eB = 0.708 +Initial Condition(s): xA(t) = -0.167, xB(t) = 0.666, xC(t) = 0.500, eA(t) = 0.666, eC(t) = 0.833 +Error: 2.6110e-04 +Time: -5.0000e-01 + +Parameter(s) : k1 = 0.250, k2 = 0.500, eB = 0.750 +Initial Condition(s): xA(t) = 0.166, xB(t) = 0.333, xC(t) = 0.500, eA(t) = 0.666, eC(t) = 0.833 +Error: 2.6110e-04 +Time: -5.0000e-01 + + 10.135238 seconds (7.58 M allocations: 508.869 MiB, 3.08% gc time, 85.32% compilation time) +For model slowfast: Max abs rel. err: 2.006024059928685 +For model slowfast: Max abs rel. err: 9.966318433180567e-9 +DEBUG +┌ Warning: No solution found +└ @ ParameterEstimation ~/.julia/dev/ParameterEstimation/src/estimation/utils.jl:28 +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [k1, k2, eB, eA, xC, xB, xA, eC] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [eA, xC, eC] +[ Info: Locally but not globally identifiable parameters: [k1, k2, xB, eB, xA] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] +[ Info: Best estimate found +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [a_0] +[ Info: Locally identifiable parameters: [d, Tr, N, In] +[ Info: Not identifiable parameters: [a, b, g, nu, S] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [d, Tr, N, In] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [a, b, S, g, nu] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] + Computing mixed cells... 52 Time: 0:00:00 + mixed_volume: 54  +  Computing mixed cells... 133 Time: 0:00:00 + mixed_volume: 143  +  Computing mixed cells... 180 Time: 0:00:01 + mixed_volume: 221 +3 Froissart doublets. Number of residues = 10 +436.213103 seconds (30.90 M allocations: 2.446 GiB, 0.22% gc time, 3.03% compilation time) +DEBUG +┌ Warning: No solution found +└ @ ParameterEstimation ~/.julia/dev/ParameterEstimation/src/estimation/utils.jl:28 +[ Info: Preproccessing `ModelingToolkit.ODESystem` object +[ Info: Solving the problem +[ Info: Constructing the maximal system +[ Info: Truncating +[ Info: Assessing local identifiability +[ Info: Found Pivots: [] +[ Info: Locally identifiable parameters: [mu_N, mu_EE, mu_LE, mu_LL, mu_M, mu_P, mu_PE, mu_PL, delta_NE, delta_EL, delta_LM, rho_E, rho_P, N, M, P, S, E] +[ Info: Not identifiable parameters: [] +[ Info: Randomizing +[ Info: Gröbner basis computation +[ Info: Remainder computation +[ Info: === Summary === +[ Info: Globally identifiable parameters: [N, delta_EL, mu_PL, M, mu_M, rho_E, mu_LL, mu_PE, mu_LE, mu_EE, E, mu_N, rho_P, P, S, delta_NE, mu_P, delta_LM] +[ Info: Locally but not globally identifiable parameters: [] +[ Info: Not identifiable parameters: [] +[ Info: =============== +[ Info: Estimating via the interpolators: ["AAA"] + Computing mixed cells... 15 Time: 0:00:00 + mixed_volume: 24  +  Computing mixed cells... 29 Time: 0:00:00 + mixed_volume: 48  +  Computing mixed cells... 33 Time: 0:00:00 + mixed_volume: 53 +[ Info: Best estimate found +4 Froissart doublets. Number of residues = 10 +2 Froissart doublets. Number of residues = 8 +3 Froissart doublets. Number of residues = 9 +4 Froissart doublets. Number of residues = 9 +Parameter(s) : mu_N = 0.071, mu_EE = 0.143, mu_LE = -0.082, mu_LL = -0.027, mu_M = -0.398, mu_P = 0.546, mu_PE = 0.036, mu_PL = 0.013, delta_NE = 0.643, delta_EL = 0.714, delta_LM = 0.455, rho_E = 0.857, rho_P = 0.368 +Initial Condition(s): N(t) = 0.167, E(t) = 0.333, S(t) = -11.359, M(t) = 12.526, P(t) = 0.833 +Error: 1.0260e-03 +Time: -5.0000e-01 + + 36.450014 seconds (26.74 M allocations: 2.377 GiB, 2.50% gc time, 42.42% compilation time) +For model Crauste: Max abs rel. err: 23.71880485343265 diff --git a/benchmarks/just_daisy_mamil4.jl b/benchmarks/just_daisy_mamil4.jl new file mode 100644 index 00000000..b93487af --- /dev/null +++ b/benchmarks/just_daisy_mamil4.jl @@ -0,0 +1,479 @@ +using ParameterEstimation +using ModelingToolkit, DifferentialEquations#, Plots +using BenchmarkTools + + +struct ParameterEstimationProblem + Name::Any + model::Any + measured_quantities::Any + data_sample::Any + solver::Any + p_true::Any + ic::Any +end + + +function biohydrogenation(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters k5 k6 k7 k8 k9 k10 + @variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t) + D = Differential(t) + states = [x4, x5, x6, x7] + parameters = [k5, k6, k7, k8, k9, k10] + + @named model = ODESystem([ + D(x4) ~ -k5 * x4 / (k6 + x4), + D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6), + D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10, + D(x7) ~ k9 * x6 * (k10 - x6) / k10, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x4, + y2 ~ x5, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("BioHydrogenation", model, measured_quantities, data_sample, solver, p_true, ic) +end + +function crauste(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters mu_N mu_EE mu_LE mu_LL mu_M mu_P mu_PE mu_PL delta_NE delta_EL delta_LM rho_E rho_P + @variables t N(t) E(t) S(t) M(t) P(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [N, E, S, M, P] + parameters = [ + mu_N, + mu_EE, + mu_LE, + mu_LL, + mu_M, + mu_P, + mu_PE, + mu_PL, + delta_NE, + delta_EL, + delta_LM, + rho_E, + rho_P, + ] + @named model = ODESystem( + [ + D(N) ~ -N * mu_N - N * P * delta_NE, + D(E) ~ N * P * delta_NE - E^2 * mu_EE - + E * delta_EL + E * P * rho_E, + D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL - + E * S * mu_LE, + D(M) ~ S * delta_LM - mu_M * M, + D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE - + S * P * mu_PL, + ], t, states, parameters) + measured_quantities = [y1 ~ N, y2 ~ E, y3 ~ S + M, y4 ~ P] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.071, 0.143, 0.214, 0.286, 0.357, 0.429, 0.5, 0.571, 0.643, 0.714, 0.786, 0.857, 0.929] # True Parameters + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Crauste", model, measured_quantities, data_sample, solver, p_true, ic) +end + + +function daisy_ex3(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters p1 p3 p4 p6 p7 + @variables t x1(t) x2(t) x3(t) u0(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, u0] + parameters = [p1, p3, p4, p6, p7] + @named model = ODESystem([ + D(x1) ~ -1 * p1 * x1 + x2 + u0, + D(x2) ~ p3 * x1 - p4 * x2 + x3, + D(x3) ~ p6 * x1 - p7 * x3, + D(u0) ~ 1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ u0, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_ex3", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + + +function daisy_mamil3(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters a12 a13 a21 a31 a01 + @variables t x1(t) x2(t) x3(t) y1(t) y2(t) + D = Differential(t) + + ic = [0.25, 0.5, 0.75] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + states = [x1, x2, x3] + parameters = [a12, a13, a21, a31, a01] + @named model = ODESystem([D(x1) ~ -(a21 + a31 + a01) * x1 + a12 * x2 + a13 * x3, + D(x2) ~ a21 * x1 - a12 * x2, + D(x3) ~ a31 * x1 - a13 * x3], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil3", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + +function daisy_mamil4(datasize = 20, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters k01, k12, k13, k14, k21, k31, k41 + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) + D = Differential(t) + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] # True Parameters + + states = [x1, x2, x3, x4] + parameters = [k01, k12, k13, k14, k21, k31, k41] + @named model = ODESystem([D(x1) ~ -k01 * x1 + k12 * x2 + k13 * x3 + k14 * x4 - k21 * x1 - k31 * x1 - k41 * x1, + D(x2) ~ -k12 * x2 + k21 * x1, + D(x3) ~ -k13 * x3 + k31 * x1, + D(x4) ~ -k14 * x4 + k41 * x1], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2, y3 ~ x3 + x4] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil4", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + +function fitzhugh_nagumo(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters g a b + @variables t V(t) R(t) y1(t) y2(t) + D = Differential(t) + states = [V, R] + parameters = [g, a, b] + + ic = [0.333, 0.67] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ V] + + @named model = ODESystem([ + D(V) ~ g * (V - V^3 / 3 + R), + D(R) ~ 1 / g * (V - a + b * R), + ], t, states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("fitzhugh-nagumo", model, measured_quantities, data_sample, solver, p_true, ic) +end + + + +function hiv_local(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters b c d k1 k2 mu1 mu2 q1 q2 s + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, x4] + parameters = [b, c, d, k1, k2, mu1, mu2, q1, q2, s] + + @named model = ODESystem([ + D(x1) ~ -b * x1 * x4 - d * x1 + s, + D(x2) ~ b * q1 * x1 * x4 - k1 * x2 - mu1 * x2, + D(x3) ~ b * q2 * x1 * x4 + k1 * x2 - mu2 * x3, + D(x4) ~ -c * x4 + k2 * x3, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x4, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.091, 0.182, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + time_interval = [-0.5, 0.5] + datasize = 20 + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("hiv_local", model, measured_quantities, data_sample, solver, p_true, ic) +end + + +function hiv(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters lm d beta a k u c q b h + @variables t x(t) y(t) v(t) w(t) z(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [x, y, v, w, z] + parameters = [lm, d, beta, a, k, u, c, q, b, h] + + @named model = ODESystem([ + D(x) ~ lm - d * x - beta * x * v, + D(y) ~ beta * x * v - a * y, + D(v) ~ k * y - u * v, + D(w) ~ c * x * y * w - c * q * y * w - b * w, + D(z) ~ c * q * y * w - h * z, + ], t, states, parameters) + measured_quantities = [y1 ~ w, y2 ~ z, y3 ~ x, y4 ~ y + v] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.091, 0.181, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + + res = ParameterEstimation.estimate(model, measured_quantities, data_sample; + solver = solver) + return ParameterEstimationProblem("hiv", model, measured_quantities, data_sample, solver, p_true, ic) +end + + + +function lotka_volterra(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + @parameters k1 k2 k3 + @variables t r(t) w(t) y1(t) + D = Differential(t) + ic = [0.333, 0.667] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ r] + states = [r, w] + parameters = [k1, k2, k3] + + @named model = ODESystem([D(r) ~ k1 * r - k2 * r * w, D(w) ~ k2 * r * w - k3 * w], t, + states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Lotka_Volterra", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + +function seir(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters a b nu + @variables t S(t) E(t) In(t) N(t) y1(t) y2(t) + D = Differential(t) + states = [S, E, In, N] + parameters = [a, b, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N, + D(E) ~ b * S * In / N - nu * E, + D(In) ~ nu * E - a * In, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ In, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.25, 0.5, 0.75] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("SEIR", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + + + +function simple(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ -a * x2, + D(x2) ~ 1 / b * (x1), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("simple", model, measured_quantities, data_sample, solver, p_true, ic) + +end + + + +function sirsforced(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) + + @parameters b0 b1 g M mu nu + @variables t i(t) r(t) s(t) x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [i, r, s, x1, x2] + parameters = [b0, b1, g, M, mu, nu] + + @named model = ODESystem([ + D(s) ~ mu - mu * s - b0 * (1 + b1 * x1) * i * s + g * r, + D(i) ~ b0 * (1 + b1 * x1) * i * s - (nu + mu) * i, + D(r) ~ nu * i - (mu + g) * r, + D(x1) ~ -M * x2, + D(x2) ~ M * x1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ i, + y2 ~ r, + ] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("sirsforced", model, measured_quantities, data_sample, solver, p_true, ic) + +end + +function slowfast(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) # TODO(orebas):in the old code it was CVODE_BDF. should we go back to that? + #solver = CVODE_BDF() + @parameters k1 k2 eB + @variables t xA(t) xB(t) xC(t) eA(t) eC(t) y1(t) y2(t) y3(t) y4(t) #eA(t) eC(t) + D = Differential(t) + states = [xA, xB, xC, eA, eC] + parameters = [k1, k2, eB] + @named model = ODESystem([ + D(xA) ~ -k1 * xA, + D(xB) ~ k1 * xA - k2 * xB, + D(xC) ~ k2 * xB, + D(eA) ~ 0, + D(eC) ~ 0, + ], t, states, parameters) + + measured_quantities = [y1 ~ xC, y2 ~ eA * xA + eB * xB + eC * xC, y3 ~ eA, y4 ~ eC] + ic = [0.166, 0.333, 0.5, 0.666, 0.833] + p_true = [0.25, 0.5, 0.75] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("slowfast", model, measured_quantities, data_sample, solver, p_true, ic) + +end + +function treatment(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) #note the solver. Tsit5 apparently can't handle mass matrices + @parameters a b d g nu + @variables t In(t) N(t) S(t) Tr(t) y1(t) y2(t) + D = Differential(t) + states = [In, N, S, Tr] + parameters = [a, b, d, g, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N - d * b * S * Tr / N, + D(In) ~ b * S * In / N + d * b * S * Tr / N - (a + g) * In, + D(Tr) ~ g * In - nu * Tr, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ Tr, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize, solver = solver) + return ParameterEstimationProblem("treatment", model, measured_quantities, data_sample, solver, p_true, ic) +end + + +function vanderpol(datasize = 21, time_interval = [-0.5, 0.5], solver = Tsit5()) + + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ a * x2, + D(x2) ~ -(x1) - b * (x1^2 - 1) * (x2), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("vanderpol", model, measured_quantities, data_sample, solver, p_true, ic) +end + + + + + + + + + +function analyze_parameter_estimation_problem(PEP::ParameterEstimationProblem) + @time res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, PEP.data_sample; + solver = PEP.solver , interpolators = Dict("AAA" => ParameterEstimation.aaad)) + all_params = vcat(PEP.p_true) + for each in res + estimates = vcat(collect(values(each.parameters))) + println("For model ", PEP.Name, ": Max abs rel. err: ", + maximum(abs.((estimates .- all_params) ./ (all_params)))) + end +end + + +function analyze_parameter_estimation_problem_old(PEP::ParameterEstimationProblem) + res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, PEP.data_sample; + solver = PEP.solver, interpolator = ("AAA" => ParameterEstimation.aaad)) + all_params = vcat(PEP.ic, PEP.p_true) + for each in res + estimates = vcat(collect(values(each.states)), collect(values(each.parameters))) + println("Max abs rel. err: ", + maximum(abs.((estimates .- all_params) ./ (all_params)))) + end +end + +function main() + for PEP in [ + daisy_mamil4(), + ] + + analyze_parameter_estimation_problem(PEP) + + end +end + +main() diff --git a/benchmarks/run.txt b/benchmarks/run.txt new file mode 100644 index 00000000..5ed6c9d1 --- /dev/null +++ b/benchmarks/run.txt @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/benchmarks/supersimple.jl b/benchmarks/supersimple.jl new file mode 100644 index 00000000..a420a793 --- /dev/null +++ b/benchmarks/supersimple.jl @@ -0,0 +1,30 @@ +using ParameterEstimation +using ModelingToolkit, DifferentialEquations +solver = Tsit5() + +@parameters a b c d +@variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) y4(t) +D = Differential(t) +states = [x1, x2, x3, x4] +parameters = [a, b, c, d] + +@named model = ODESystem([ + D(x1) ~ a + x2, + D(x2) ~ b + x3, + D(x3) ~ c + x4, + D(x4) ~ d, + ], t, states, parameters) +measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + y3 ~ x3, + y4 ~ x4, +] + +ic = [0.0, 0.0, 0.0, 0.0] +p_true = [2.0, 3.0, 4.0, 5.0] +time_interval = [-4.0, 4.0] +datasize = 9 +data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) +res = ParameterEstimation.estimate(model, measured_quantities, data_sample) diff --git a/examples/all-global/coupled_pendulum.jl b/examples/all-global/coupled_pendulum.jl deleted file mode 100644 index f31e11e7..00000000 --- a/examples/all-global/coupled_pendulum.jl +++ /dev/null @@ -1,36 +0,0 @@ -using ParameterEstimation -using ModelingToolkit, DifferentialEquations -solver = Tsit5() - -@parameters alpha -@variables t theta_1(t) theta_1_dot(t) theta_2(t) theta_2_dot(t) y1(t) y2(t) -D = Differential(t) - -#coupled pendulum equations -states = [theta_1, theta_1_dot, theta_2, theta_2_dot] -parameters = [alpha] -@named model = ODESystem([ - D(theta_1) ~ theta_1_dot, - D(theta_1_dot) ~ -theta_1 * (alpha + 1) + alpha * theta_2, - D(theta_2) ~ theta_2_dot, - D(theta_2_dot) ~ alpha * theta_1 - theta_2 * (alpha + 1), - ], t, states, parameters) - -#initial conditions -ic = [1.0, 0.0, 0.0, -1.0] -p_true = [1.0] -time_interval = [0.0, 10.0] -datasize = 32 - -v = randn(datasize) -v = sort((v .- minimum(v)) / (maximum(v) - minimum(v))) * time_interval[2] - -measured_quantities = [y1 ~ theta_1, y2 ~ theta_2 + theta_1] -data_sample = Dict{Any, Any}("t" => v) -data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, - p_true, ic, datasize; solver = solver, - uneven_sampling = true, - uneven_sampling_times = data_sample["t"]) - -res = ParameterEstimation.estimate(model, measured_quantities, data_sample; - solver = solver) diff --git a/examples/all-global/crauste.jl b/examples/all-global/crauste.jl index 2f40f231..c2642ae4 100644 --- a/examples/all-global/crauste.jl +++ b/examples/all-global/crauste.jl @@ -1,43 +1,44 @@ using ParameterEstimation using ModelingToolkit, DifferentialEquations -solver = Tsit5() +solver = Vern9() @parameters mu_N mu_EE mu_LE mu_LL mu_M mu_P mu_PE mu_PL delta_NE delta_EL delta_LM rho_E rho_P @variables t N(t) E(t) S(t) M(t) P(t) y1(t) y2(t) y3(t) y4(t) D = Differential(t) states = [N, E, S, M, P] parameters = [ - mu_N, - mu_EE, - mu_LE, - mu_LL, - mu_M, - mu_P, - mu_PE, - mu_PL, - delta_NE, - delta_EL, - delta_LM, - rho_E, - rho_P, + mu_N, + mu_EE, + mu_LE, + mu_LL, + mu_M, + mu_P, + mu_PE, + mu_PL, + delta_NE, + delta_EL, + delta_LM, + rho_E, + rho_P, ] -@named model = ODESystem([ - D(N) ~ -N * mu_N - N * P * delta_NE, - D(E) ~ N * P * delta_NE - E^2 * mu_EE - - E * delta_EL + E * P * rho_E, - D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL - - E * S * mu_LE, - D(M) ~ S * delta_LM - mu_M * M, - D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE - - S * P * mu_PL, - ], t, states, parameters) +@named model = ODESystem( + [ + D(N) ~ -N * mu_N - N * P * delta_NE, + D(E) ~ N * P * delta_NE - E^2 * mu_EE - + E * delta_EL + E * P * rho_E, + D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL - + E * S * mu_LE, + D(M) ~ S * delta_LM - mu_M * M, + D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE - + S * P * mu_PL, + ], t, states, parameters) measured_quantities = [y1 ~ N, y2 ~ E, y3 ~ S + M, y4 ~ P] ic = [1.0, 1.0, 1.0, 1.0, 1.0] -time_interval = [0.0, 1.0] -datasize = 20 +time_interval = [-0.5, 0.5] +datasize = 21 p_true = [1, 1.3, 1.1, 1.2, 1.1, 1, 0.5, 1.0, 1.0, 1.0, 1.0, 0.9, 1.2] # True Parameters data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, - p_true, ic, datasize; solver = solver) + p_true, ic, datasize; solver = solver, abstol = 1.0e-14, reltol = 1.0e-14) res = ParameterEstimation.estimate(model, measured_quantities, data_sample; - solver = solver) + solver = solver) diff --git a/examples/all-global/hiv.jl b/examples/all-global/hiv.jl index 1bfc69cb..308911dd 100644 --- a/examples/all-global/hiv.jl +++ b/examples/all-global/hiv.jl @@ -9,12 +9,12 @@ states = [x, y, v, w, z] parameters = [lm, d, beta, a, k, u, c, q, b, h] @named model = ODESystem([ - D(x) ~ lm - d * x - beta * x * v, - D(y) ~ beta * x * v - a * y, - D(v) ~ k * y - u * v, - D(w) ~ c * x * y * w - c * q * y * w - b * w, - D(z) ~ c * q * y * w - h * z, - ], t, states, parameters) + D(x) ~ lm - d * x - beta * x * v, + D(y) ~ beta * x * v - a * y, + D(v) ~ k * y - u * v, + D(w) ~ c * x * y * w - c * q * y * w - b * w, + D(z) ~ c * q * y * w - h * z, + ], t, states, parameters) measured_quantities = [y1 ~ w, y2 ~ z, y3 ~ x, y4 ~ y + v] ic = [1.0, 1.0, 1.0, 1.0, 1.0] @@ -22,9 +22,9 @@ time_interval = [0.0, 20.0] datasize = 10 p_true = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, - p_true, ic, datasize; solver = solver) + p_true, ic, datasize; solver = solver) -identifiability_result = ParameterEstimation.check_identifiability(model; - measured_quantities = measured_quantities) +#identifiability_result = ParameterEstimation.check_identifiability(model; +# measured_quantities = measured_quantities) res = ParameterEstimation.estimate(model, measured_quantities, data_sample; - solver = solver) + solver = solver) diff --git a/examples/all-global/lotka-volterra.jl b/examples/all-global/lotka-volterra.jl index 7f8e2de8..a5435dc3 100644 --- a/examples/all-global/lotka-volterra.jl +++ b/examples/all-global/lotka-volterra.jl @@ -16,11 +16,12 @@ states = [r, w] parameters = [k1, k2, k3] @named model = ODESystem([D(r) ~ k1 * r - k2 * r * w, - D(w) ~ k2 * r * w - k3 * w], - t, states, parameters) + D(w) ~ k2 * r * w - k3 * w], + t, states, parameters) data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, - p_true, ic, datasize; solver = solver) + p_true, ic, datasize; solver = solver) + res = ParameterEstimation.estimate(model, measured_quantities, data_sample; - solver = solver) + solver = solver) diff --git a/examples/all-global/supersimple.jl b/examples/all-global/supersimple.jl new file mode 100644 index 00000000..a420a793 --- /dev/null +++ b/examples/all-global/supersimple.jl @@ -0,0 +1,30 @@ +using ParameterEstimation +using ModelingToolkit, DifferentialEquations +solver = Tsit5() + +@parameters a b c d +@variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) y4(t) +D = Differential(t) +states = [x1, x2, x3, x4] +parameters = [a, b, c, d] + +@named model = ODESystem([ + D(x1) ~ a + x2, + D(x2) ~ b + x3, + D(x3) ~ c + x4, + D(x4) ~ d, + ], t, states, parameters) +measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + y3 ~ x3, + y4 ~ x4, +] + +ic = [0.0, 0.0, 0.0, 0.0] +p_true = [2.0, 3.0, 4.0, 5.0] +time_interval = [-4.0, 4.0] +datasize = 9 +data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) +res = ParameterEstimation.estimate(model, measured_quantities, data_sample) diff --git a/examples/unidentifiable/biohydrogenation.jl b/examples/unidentifiable/biohydrogenation.jl index f82585a3..960484e8 100644 --- a/examples/unidentifiable/biohydrogenation.jl +++ b/examples/unidentifiable/biohydrogenation.jl @@ -29,3 +29,6 @@ data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_i p_true, ic, datasize; solver = solver) res = ParameterEstimation.estimate(model, measured_quantities, data_sample) + + +alg= Tsit5() \ No newline at end of file diff --git a/src/ParameterEstimation.jl b/src/ParameterEstimation.jl index 995615db..54e93a72 100644 --- a/src/ParameterEstimation.jl +++ b/src/ParameterEstimation.jl @@ -1,17 +1,62 @@ module ParameterEstimation -import DifferentialEquations: Tsit5 -import TaylorSeries: Taylor1 -import OrderedCollections: OrderedDict -using ProgressMeter, Logging, Printf -using ModelingToolkit, LinearSolve, LinearAlgebra -using SIAN, HomotopyContinuation, Groebner, Oscar -using .ReturnCode -import StructuralIdentifiability: eval_at_nemo, ODE +using PrecompileTools +@recompile_invalidations begin -Float = Union{Float64, Float32, Float16} + import DifferentialEquations: Tsit5 + import DifferentialEquations: Vern9 + import TaylorSeries: Taylor1 + import OrderedCollections: OrderedDict + using BaryRational: BaryRational + using ForwardDiff: ForwardDiff + using LinearAlgebra: LinearAlgebra + using TaylorDiff: TaylorDiff + using Suppressor + + using ProgressMeter, Logging, Printf + using ModelingToolkit, LinearSolve, LinearAlgebra + using SIAN, HomotopyContinuation, Groebner, Oscar + using .ReturnCode + import StructuralIdentifiability: eval_at_nemo, ODE + using BaryRational + using ForwardDiff + using ArbNumerics + +end +Float = Union{Float64, Float32, Float16, BigFloat} include("includes.jl") export check_identifiability, estimate, filter_solutions +@recompile_invalidations begin + @compile_workload begin + + @parameters a b c d + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [x1, x2, x3, x4] + parameters = [a, b, c, d] + @named model = ODESystem([ + D(x1) ~ a + x2, + D(x2) ~ b + x3, + D(x3) ~ c + x4, + D(x4) ~ d, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + y3 ~ x3, + y4 ~ x4, + ] + + ic = [0.0, 0.0, 0.0, 0.0] + p_true = [2.0, 3.0, 4.0, 5.0] + time_interval = [-4.0, 4.0] + datasize = 9 + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize;) + res = ParameterEstimation.estimate(model, measured_quantities, data_sample) + end +end end + diff --git a/src/estimation/estimate.jl b/src/estimation/estimate.jl index b0ced3c0..fcf3a6e0 100644 --- a/src/estimation/estimate.jl +++ b/src/estimation/estimate.jl @@ -1,54 +1,61 @@ """ - estimate(model::ModelingToolkit.ODESystem, - measured_quantities::Vector{ModelingToolkit.Equation}, - data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}(); - at_time::T = 0.0, method = :homotopy, solver = Tsit5(), - degree_range = nothing, real_tol = 1e-10, - threaded = Threads.nthreads() > 1) where {T <: Float} + estimate(model::ModelingToolkit.ODESystem, + measured_quantities::Vector{ModelingToolkit.Equation}, + data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + at_time::T = 0.0, method = :homotopy, solver = Tsit5(), + degree_range = nothing, real_tol = 1e-12, + threaded = Threads.nthreads() > 1) where {T <: Float} Run estimation over a range of interpolation degrees. Return the best estimate according to a heuristic: - - the best estimate is the one with the smallest error between sample data and ODE solution with current parameters (estimates); + - the best estimate is the one with the smallest error between sample data and ODE solution with current parameters (estimates); # Arguments - `model::ModelingToolkit.ODESystem`: the ODE model; - `measured_quantities::Vector{ModelingToolkit.Equation}`: the measured quantities (output functions that were sampled experimentally); - `data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}()`: the data sample, a dictionary with keys being the measured quantities and - values being the corresponding data. Must include the time vector; + values being the corresponding data. Must include the time vector; - `at_time::T = 0.0`: the time used for derivative computation; +- `report_time = nothing`: specify a time T, at which the initial conditions (state variables) will be estimated. If "nothing", use the leftmost time. - `method = :homotopy`: the method used for polynomial system solving. Can be one of :homotopy (recommended) or :msolve; -- `solver`: the ODE solver used for ODE solution computation (default: Tsit5()); -- `degree_range = nothing`: the range of interpolation degrees to be used. If `nothing`, the range is computed automatically; -- `real_tol` = 1e-10: the tolerance used for real root finding; +- `solver`: the ODE solver used for ODE solution computation (default: Vern9()); +- `interpolators = nothing`: the set of interpolators to be used. See examples. If `nothing`, a default is used which includes AAA, FLoater-Hormann, and Fourier interpolations; +- `real_tol` = 1e-14: the tolerance used for real root finding; - `threaded = Threads.nthreads() > 1`: whether to use multiple threads for computation (determined automatically). # Returns - `result::Vector{EstimationResult}`: the result of the estimation, a vector of `EstimationResult` objects. """ function estimate(model::ModelingToolkit.ODESystem, - measured_quantities::Vector{ModelingToolkit.Equation}, - data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); - inputs::Vector{ModelingToolkit.Equation} = Vector{ModelingToolkit.Equation - }(), - at_time::T = 0.0, method = :homotopy, solver = Tsit5(), - degree_range = nothing, real_tol = 1e-10, - threaded = Threads.nthreads() > 1) where {T <: Float} - if !(method in [:homotopy, :msolve]) - throw(ArgumentError("Method $method is not supported, must be one of :homotopy or :msolve.")) - end - if threaded - result = estimate_threaded(model, measured_quantities, inputs, data_sample; - at_time = at_time, solver = solver, - degree_range = degree_range, - method = method, - real_tol = real_tol) - else - result = estimate_serial(model, measured_quantities, inputs, data_sample; - solver = solver, at_time = at_time, - degree_range = degree_range, method = method, - real_tol = real_tol) - end - for each in result - display(each) - end - return result + measured_quantities::Vector{ModelingToolkit.Equation}, + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + inputs::Vector{ModelingToolkit.Equation} = Vector{ModelingToolkit.Equation}(), + at_time::T = data_sample["t"][fld(length((data_sample["t"])), 2)], #uses something akin to a midpoint by default + method = :homotopy, solver = Vern9(), + report_time = minimum(data_sample["t"]), + interpolators = nothing, real_tol = 1e-14, + threaded = Threads.nthreads() > 1, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing) where {T <: Float} + + #println("DEBUG") + if !(method in [:homotopy, :msolve]) + throw(ArgumentError("Method $method is not supported, must be one of :homotopy or :msolve.")) + end + if threaded + result = estimate_threaded(model, measured_quantities, inputs, data_sample; + at_time = at_time, solver = solver, + interpolators = interpolators, + method = method, + real_tol = real_tol) + else + result = estimate_serial(model, measured_quantities, + inputs, + data_sample; + solver = solver, at_time = at_time, report_time, + interpolators = interpolators, method = method, + real_tol = real_tol, filtermode, parameter_constraints = parameter_constraints, ic_constraints = ic_constraints) + end + println("Final Results:") + for each in result + display(each) + end + return result end diff --git a/src/estimation/fixed_degree.jl b/src/estimation/fixed_degree.jl index 826eaa48..d1e1e9ee 100644 --- a/src/estimation/fixed_degree.jl +++ b/src/estimation/fixed_degree.jl @@ -1,13 +1,62 @@ """ - estimate_fixed_degree(model::ModelingToolkit.ODESystem, - measured_quantities::Vector{ModelingToolkit.Equation}, - inputs::Vector{ModelingToolkit.Equation}, - data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}(); - identifiability_result = Dict{String, Any}(), - interpolation_degree::Int = 1, - at_time::T = 0.0, - method = :homotopy, - real_tol = 1e-10) where {T <: Float} + backsolve_initial_conditions(model, + E, report_time, inputs::Vector{Equation}, data_sample; + solver = Vern9(), abstol = 1e-14, reltol = 1e-14) + initial_conditions = [E[s] for s in ModelingToolkit.states(model)] + parameter_values = [E[p] for p in ModelingToolkit.parameters(model)] + tspan = (E.at_time, report_time) + + Given a set of estimated state variables at E.at_time, solves ODE backwards to estimate state variables at report_time. In most cases tspan will be backwards. + """ +function backsolve_initial_conditions(model, E, report_time, inputs::Vector{Equation}, data_sample; + solver = Vern9(), abstol = 1e-14, reltol = 1e-14) + initial_conditions = [E[s] for s in ModelingToolkit.states(model)] + parameter_values = [E[p] for p in ModelingToolkit.parameters(model)] + tspan = (E.at_time, report_time) #this is almost always backwards! + + ode_equations = ModelingToolkit.equations(model) + ode_equations = substitute(ode_equations, + Dict(each.lhs => Num(each.rhs) for each in inputs)) + t = ModelingToolkit.get_iv(model) + @named new_model = ODESystem(ode_equations, t, ModelingToolkit.states(model), + ModelingToolkit.parameters(model)) + prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values) + saveat = range(tspan[1], tspan[2], length = length(data_sample["t"])) + + ode_solution = ModelingToolkit.solve(prob, solver, p = parameter_values, saveat = saveat, abstol = abstol, reltol = reltol) + + state_param_map = (Dict(x => replace(string(x), "(t)" => "") + for x in ModelingToolkit.states(model))) + + + newstates = copy(E.states) + + for s in ModelingToolkit.states(model) + temp = ode_solution[Symbol(state_param_map[s])][end] + newstates[s] = temp + end + ER = EstimationResult(E.parameters, newstates, E.degree, report_time, + E.err, E.interpolants, E.return_code, E.datasize, report_time) + return ER + +end + + + + + + +""" + estimate_single_interpolator(model::ModelingToolkit.ODESystem, + measured_quantities::Vector{ModelingToolkit.Equation}, + inputs::Vector{ModelingToolkit.Equation}, + data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + identifiability_result = Dict{String, Any}(), + interpolator = ("AAA" => aaad), + at_time::T, + report_time = minimum(data_sample["t"]), + method = :homotopy, + real_tol = 1e-14) where {T <: Float} time_interval = [minimum(data_sample["t"]), maximum(data_sample["t"]) where {T <: Float} Estimate the parameters of a model using the data sample `data_sample` and the measured quantities `measured_quantities`. @@ -17,56 +66,67 @@ measured quantities `measured_quantities`. - `measured_quantities::Vector{ModelingToolkit.Equation}`: the measured quantities of the model. Used for identifiability assessment. - `inputs::Vector{ModelingToolkit.Equation}`: the input equations of the model. - `data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}()`: the data sample used for estimation (same functions as `measured_quantities`). - The keys of the dictionary are the measured quantities - and the values are the corresponding data samples. + The keys of the dictionary are the measured quantities + and the values are the corresponding data samples. - `identifiability_result = Dict{String, Any}()`: the result of the identifiability assessment. -- `interpolation_degree::Int = 1`: the degree of the polynomial interpolation. +- `interpolator = ("AAA" => aaad)`: the interpolator to use, see examples. - `at_time::T = 0.0`: the time at which the parameters are estimated. +- `report_time::T `: the time at which the state variables are to be reported (by default, the earliest time). These are backsolved from at_time. - `method = :homotopy`: the method used for solving the polynomial system. Can be one of :homotopy (recommended) and :msolve. -- `real_tol = 1e-10`: the tolerance for the real solutions of the polynomial system. +- `real_tol = 1e-14`: the tolerance for the real solutions of the polynomial system. # Returns - `EstimationResult`: the estimated parameters and initial conditions of the model. """ -function estimate_fixed_degree(model::ModelingToolkit.ODESystem, - measured_quantities::Vector{ModelingToolkit.Equation}, - inputs::Vector{ModelingToolkit.Equation}, - data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, - Vector{T}}(); - identifiability_result = Dict{String, Any}(), - interpolation_degree::Int = 1, - at_time::T = 0.0, - method = :homotopy, - real_tol = 1e-10) where {T <: Float} - time_interval = [minimum(data_sample["t"]), maximum(data_sample["t"])] - check_inputs(measured_quantities, data_sample, interpolation_degree) - datasize = length(first(values(data_sample))) - parameters = ModelingToolkit.parameters(model) - states = ModelingToolkit.states(model) - num_parameters = length(parameters) + length(states) - @debug "Interpolating sample data via rational interpolation" - if !haskey(data_sample, "t") - @warn "No sampling time points found in data sample. Assuming uniform sampling t ∈ [$(time_interval[1]), $(time_interval[2])]." - data_sample["t"] = range(time_interval[1], time_interval[2], length = datasize) - end - interpolants = ParameterEstimation.interpolate(identifiability_result, - data_sample, measured_quantities, inputs; - interpolation_degree = interpolation_degree, - diff_order = num_parameters + 1, - at_t = at_time, - method = method) - - if method == :homotopy - all_solutions = solve_via_homotopy(identifiability_result, model; - real_tol = real_tol) - elseif method == :msolve - all_solutions = solve_via_msolve(identifiability_result, model; - real_tol = real_tol) - else - throw(ArgumentError("Method $method not supported")) - end - all_solutions = [EstimationResult(model, each, interpolation_degree, at_time, - interpolants, ReturnCode.Success, datasize) - for each in all_solutions] - return all_solutions +function estimate_single_interpolator(model::ModelingToolkit.ODESystem, + measured_quantities::Vector{ModelingToolkit.Equation}, + inputs::Vector{ModelingToolkit.Equation}, + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, + Vector{T}}(); + identifiability_result = Dict{String, Any}(), + interpolator = ("AAA" => aaad), + at_time::T, report_time = minimum(data_sample["t"]), + method = :homotopy, + real_tol = 1e-14) where {T <: Float} + time_interval = [minimum(data_sample["t"]), maximum(data_sample["t"])] #TODO(orebas) will this break if key is missing? + + + check_inputs(measured_quantities, data_sample) #TODO(orebas): I took out checking the degree. Do we want to check the interpolator otherwise? + datasize = length(first(values(data_sample))) + parameters = ModelingToolkit.parameters(model) + states = ModelingToolkit.states(model) + num_parameters = length(parameters) + length(states) + @debug "Interpolating sample data via interpolation method $(interpolator.first)" + if !haskey(data_sample, "t") + @warn "No sampling time points found in data sample. Assuming uniform sampling t ∈ [$(time_interval[1]), $(time_interval[2])]." + data_sample["t"] = range(time_interval[1], time_interval[2], length = datasize) + end + interpolants = ParameterEstimation.interpolate(identifiability_result, + data_sample, measured_quantities, inputs; + interpolator = interpolator, + diff_order = num_parameters + 1, #todo(OREBAS): is this always forcing num_parameters + 1 derivatives? + at_t = at_time, + method = method) + + if method == :homotopy + all_solutions = solve_via_homotopy(identifiability_result, model; + real_tol = real_tol) + elseif method == :msolve + all_solutions = solve_via_msolve(identifiability_result, model; + real_tol = real_tol) + else + throw(ArgumentError("Method $method not supported")) + end + #println(all_solutions) + #println("HERE") + all_solutions = [EstimationResult(model, each, interpolator.first, at_time, + interpolants, ReturnCode.Success, datasize, report_time) + for each in all_solutions] + + all_solutions_R = [backsolve_initial_conditions(model, each, report_time, inputs, data_sample) + for each in all_solutions] + + + #println(all_solutions_R) + return all_solutions_R end diff --git a/src/estimation/serial.jl b/src/estimation/serial.jl index 2f815618..e7070378 100644 --- a/src/estimation/serial.jl +++ b/src/estimation/serial.jl @@ -1,38 +1,42 @@ function estimate_serial(model::ModelingToolkit.ODESystem, - measured_quantities::Vector{ModelingToolkit.Equation}, - inputs::Vector{ModelingToolkit.Equation}, - data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); - at_time::T = 0.0, solver = Tsit5(), degree_range = nothing, - method = :homotopy, - real_tol::Float64 = 1e-10) where {T <: Float} + measured_quantities::Vector{ModelingToolkit.Equation}, + inputs::Vector{ModelingToolkit.Equation}, + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + at_time::T, solver = Vern9(), interpolators = nothing, + report_time = minimum(data_sample["t"]), + method = :homotopy, + real_tol::Float64 = 1e-14, filtermode = :new, parameter_constraints = nothing, + ic_constraints = nothing) where {T <: Float} check_inputs(measured_quantities, data_sample) datasize = length(first(values(data_sample))) - if degree_range === nothing - degree_range = 1:(length(data_sample[first(keys(data_sample))]) - 1) + + if isnothing(interpolators) + interpolators = default_interpolator(datasize) end id = ParameterEstimation.check_identifiability(model; - measured_quantities = measured_quantities, - inputs = [Num(each.lhs) - for each in inputs]) + measured_quantities = measured_quantities, + inputs = [Num(each.lhs) + for each in inputs]) estimates = Vector{Vector{ParameterEstimation.EstimationResult}}() - @info "Estimating via rational interpolation with degrees between $(degree_range[1]) and $(degree_range[end])" - @showprogress for deg in degree_range - unfiltered = estimate_fixed_degree(model, measured_quantities, inputs, data_sample; - identifiability_result = id, - interpolation_degree = deg, at_time = at_time, - method = method, real_tol = real_tol) + @info "Estimating via the interpolators: $(keys(interpolators))" + @showprogress for interpolator in interpolators + unfiltered = estimate_single_interpolator(model, measured_quantities, inputs, + data_sample; #TODO(orebas) we will rename estimated_fixed_degree to estimate_single_interpolator + identifiability_result = id, + interpolator = interpolator, at_time = at_time, report_time, + method = method, real_tol = real_tol) if length(unfiltered) > 0 filtered = filter_solutions(unfiltered, id, model, inputs, data_sample; - solver = solver) + solver = solver, filtermode) push!(estimates, filtered) else push!(estimates, - [ - EstimationResult(model, Dict(), deg, at_time, - Dict{Any, Interpolant}(), - ReturnCode.Failure, datasize), - ]) + [ + EstimationResult(model, Dict(), interpolator, at_time, + Dict{Any, Interpolant}(), + ReturnCode.Failure, datasize, report_time), + ]) end end - return post_process(estimates) + return post_process(estimates, filtermode, parameter_constraints, ic_constraints) end diff --git a/src/estimation/solve.jl b/src/estimation/solve.jl index 0c8b017b..ab40a79a 100644 --- a/src/estimation/solve.jl +++ b/src/estimation/solve.jl @@ -1,69 +1,71 @@ -function solve_via_homotopy(identifiability_result, model; real_tol = 1e-10) - @debug "Solving $(length(identifiability_result["polynomial_system_to_solve"])) polynomial equations in $(length(identifiability_result["polynomial_system_to_solve"].variables)) variables" +function solve_via_homotopy(identifiability_result, model; real_tol = 1e-12) + @debug "Solving $(length(identifiability_result["polynomial_system_to_solve"])) polynomial equations in $(length(identifiability_result["polynomial_system_to_solve"].variables)) variables" - polynomial_system = identifiability_result["polynomial_system_to_solve"] - state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x - for x in ModelingToolkit.states(model)), - Dict(string(x) => x for x in ModelingToolkit.parameters(model))) - results = HomotopyContinuation.solve(polynomial_system; show_progress = false) - all_solutions = HomotopyContinuation.real_solutions(results) - if length(all_solutions) == 0 - all_solutions = HomotopyContinuation.solutions(results) - if length(all_solutions) == 0 - @debug "Interpolation numerator degree $(interpolation_degree): No solutions found" - return Vector{EstimationResult}() - end - end - all_solutions_ = Vector{Dict}([]) - for sol in all_solutions - tmp = Dict() - sol = map(each -> to_exact(each; tol = real_tol), sol) - for (idx, v) in enumerate(polynomial_system.variables) - if endswith(string(v), "_0") - tmp[state_param_map[string(v)[1:(end - 2)]]] = sol[idx] - end - end - for (key, val) in identifiability_result.transcendence_basis_subs - if endswith(string(key), "_0") - tmp[state_param_map[string(key)[1:(end - 2)]]] = Int(Meta.parse("$val")) - end - end - push!(all_solutions_, tmp) - end - return all_solutions_ + polynomial_system = identifiability_result["polynomial_system_to_solve"] + state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x + for x in ModelingToolkit.states(model)), + Dict(string(x) => x for x in ModelingToolkit.parameters(model))) + results = HomotopyContinuation.solve(polynomial_system; show_progress = false) + all_solutions = HomotopyContinuation.real_solutions(results) + if length(all_solutions) == 0 + all_solutions = HomotopyContinuation.solutions(results) + if length(all_solutions) == 0 + @debug "Interpolation numerator degree $(interpolation_degree): No solutions found" + return Vector{EstimationResult}() + end + end + all_solutions_ = Vector{Dict}([]) + for sol in all_solutions + tmp = Dict() + sol = map(each -> to_exact(each; tol = real_tol), sol) + for (idx, v) in enumerate(polynomial_system.variables) + if endswith(string(v), "_0") + tmp[state_param_map[string(v)[1:(end-2)]]] = sol[idx] + end + end + for (key, val) in identifiability_result.transcendence_basis_subs + if endswith(string(key), "_0") + tmp[state_param_map[string(key)[1:(end-2)]]] = Int(Meta.parse("$val")) + end + end + push!(all_solutions_, tmp) + end + return all_solutions_ end -function solve_via_msolve(identifiability_result, model; real_tol = 1e-10) - polynomial_system = identifiability_result["polynomial_system_to_solve"] - state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x - for x in ModelingToolkit.states(model)), - Dict(string(x) => x for x in ModelingToolkit.parameters(model))) +function solve_via_msolve(identifiability_result, model; real_tol = 1e-12) + polynomial_system = identifiability_result["polynomial_system_to_solve"] + state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x + for x in ModelingToolkit.states(model)), + Dict(string(x) => x for x in ModelingToolkit.parameters(model))) - all_vars = reduce((x, y) -> union(x, y), vars(poly) for poly in polynomial_system) - R, _ = Oscar.PolynomialRing(QQ, string.(all_vars); ordering = :degrevlex) - ps = [SIAN.parent_ring_change(poly, R) for poly in polynomial_system] - i = Oscar.ideal(R, ps) - all_solutions_ = Vector{Dict}([]) + all_vars = reduce((x, y) -> union(x, y), vars(poly) for poly in polynomial_system) + R, _ = Oscar.PolynomialRing(QQ, string.(all_vars); ordering = :degrevlex) + ps = [SIAN.parent_ring_change(poly, R) for poly in polynomial_system] + i = Oscar.ideal(R, ps) + all_solutions_ = Vector{Dict}([]) - if Oscar.VERSION_NUMBER == v"0.11.3" - solutions, rat_param = Oscar.real_solutions(i) - elseif Oscar.VERSION_NUMBER == v"0.10.0" - rat_param, solutions = Oscar.msolve(i) - end - for sol in solutions - tmp = Dict() - sol = map(each -> Float64(each), sol) - for (idx, v) in enumerate(all_vars) - if endswith(string(v), "_0") - tmp[state_param_map[string(v)[1:(end - 2)]]] = sol[idx] - end - end - for (key, val) in identifiability_result.transcendence_basis_subs - if endswith(string(key), "_0") - tmp[state_param_map[string(key)[1:(end - 2)]]] = Int(Meta.parse("$val")) - end - end - push!(all_solutions_, tmp) - end - return all_solutions_ + if Oscar.VERSION_NUMBER == v"0.11.3" + solutions, rat_param = Oscar.real_solutions(i) + elseif Oscar.VERSION_NUMBER == v"0.10.0" + rat_param, solutions = Oscar.msolve(i) + else + solutions, rat_param = Oscar.real_solutions(i) + end + for sol in solutions + tmp = Dict() + sol = map(each -> Float64(each), sol) + for (idx, v) in enumerate(all_vars) + if endswith(string(v), "_0") + tmp[state_param_map[string(v)[1:(end-2)]]] = sol[idx] + end + end + for (key, val) in identifiability_result.transcendence_basis_subs + if endswith(string(key), "_0") + tmp[state_param_map[string(key)[1:(end-2)]]] = Int(Meta.parse("$val")) + end + end + push!(all_solutions_, tmp) + end + return all_solutions_ end diff --git a/src/estimation/threaded.jl b/src/estimation/threaded.jl index 2da16563..0cb4ec52 100644 --- a/src/estimation/threaded.jl +++ b/src/estimation/threaded.jl @@ -1,46 +1,50 @@ -function estimate_threaded(model, measured_quantities, inputs, data_sample; - at_time::Float = 0.0, solver = solver, - degree_range = degree_range, - method = :homotopy, real_tol::Float64 = 1e-10) - @warn "Using threaded estimation." - check_inputs(measured_quantities, data_sample) - datasize = length(first(values(data_sample))) - if isnothing(degree_range) - degree_range = 1:(length(data_sample[first(keys(data_sample))]) - 1) +function estimate_threaded(model::ModelingToolkit.ODESystem, + measured_quantities::Vector{ModelingToolkit.Equation}, + inputs::Vector{ModelingToolkit.Equation}, + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + at_time::T, solver = Vern9(), interpolators = nothing, report_time = minimum(data_sample["t"]), + method = :homotopy, real_tol::Float64 = 1e-14,filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing) where {T <: Float} + @warn "Using threaded estimation." + check_inputs(measured_quantities, data_sample) + datasize = length(first(values(data_sample))) + + if isnothing(interpolators) === nothing + interpolators = default_interpolator(datasize) end - identifiability_result = ParameterEstimation.check_identifiability(model; - measured_quantities = measured_quantities, - inputs = [Num(each.lhs) - for each in inputs]) - n_threads = Threads.nthreads() - N = length(degree_range) - estimates = Vector{Any}(nothing, n_threads) - @info "Estimating via rational interpolation with degrees between $(degree_range[1]) and $(degree_range[end]) using $n_threads threads" - Threads.@threads :static for t in 1:N - deg = degree_range[t] - id = Threads.threadid() - if isnothing(estimates[id]) - estimates[id] = [] - end - unfiltered = estimate_fixed_degree(model, measured_quantities, inputs, data_sample; - identifiability_result = identifiability_result, - interpolation_degree = deg, at_time = at_time, - method = method, - real_tol = real_tol) - if length(unfiltered) > 0 - filtered = filter_solutions(unfiltered, identifiability_result, model, inputs, - data_sample; solver = solver) - push!(estimates[id], filtered) + id = ParameterEstimation.check_identifiability(model; + measured_quantities = measured_quantities, + inputs = [Num(each.lhs) + for each in inputs]) +estimates = Vector{Vector{ParameterEstimation.EstimationResult}}() +@info "Estimating via the interpolators: $(keys(interpolators))" + n_threads = Threads.nthreads() + N = length(interpolators) + estimates = Vector{Any}(nothing, n_threads) + Threads.@threads :static for t in 1:N + interp = interp[t] + id = Threads.threadid() + if isnothing(estimates[id]) + estimates[id] = [] + end + unfiltered = estimate_single_interpolator(model, measured_quantities, inputs, data_sample; + identifiability_result = id, + interpolator=interp, at_time = at_time, report_time, + method = method, + real_tol = real_tol) + if length(unfiltered) > 0 + filtered = filter_solutions(unfiltered, id, model, inputs, + data_sample; solver = solver, filtermode) + push!(estimates[id], filtered) - else - push!(estimates[id], - [ - EstimationResult(model, Dict(), deg, at_time, - Dict{Any, Interpolant}(), ReturnCode.Failure, - datasize), - ]) - end - end - return post_process(estimates) + else + push!(estimates[id], + [ + EstimationResult(model, Dict(), interp, at_time, + Dict{Any, Interpolant}(), ReturnCode.Failure, + datasize, report_time), + ]) + end + end + return post_process(estimates,filtermode,parameter_constraints,ic_constraints) end diff --git a/src/estimation/utils.jl b/src/estimation/utils.jl index 2e5a3189..ce9b0c67 100644 --- a/src/estimation/utils.jl +++ b/src/estimation/utils.jl @@ -1,31 +1,111 @@ -function post_process(estimates) - if Threads.nthreads() > 1 - estimates = filter(x -> !isnothing(x), estimates) - estimates = vcat(estimates...) - end - #filter out the empty vectors - estimates = filter(x -> length(x) > 0, estimates) - estimates = filter(x -> x[1].return_code == ReturnCode.Success, estimates) - - best_solution = nothing - for each in estimates - if all(!isnothing(x.err) for x in each) - if best_solution === nothing - best_solution = each - else - if sum(x.err for x in each) < sum(x.err for x in best_solution) - best_solution = each - end - end - else - best_solution = nothing - end - end - if best_solution !== nothing - @info "Best estimate found" - return best_solution - else - @warn "No solution found" - return estimates - end +function check_constraints(estimate, parameter_constraints, ic_constraints) + sat = true + if (!isnothing(parameter_constraints)) + # println("Here: ", typeof(estimate)) + + # println("Here: ", estimate) + + # println("Here: ", estimate.parameters) + for (k, v) in pairs(estimate.parameters) + if (haskey(parameter_constraints, k)) + if !(parameter_constraints[k][1] <= v && v <= parameter_constraints[k][2]) + sat = false + end + end + end + end + + if (!isnothing(ic_constraints)) + for (k, v) in pairs(estimate.states) + if (haskey(ic_constraints, k)) + if !(ic_constraints[k][1] <= v && v <= ic_constraints[k][2]) + sat = false + end + end + end + end + return sat +end + + +function compare_estimation_results(a, b) + err = 0.0 + akeys = keys(a.states) + bkeys = keys(b.states) + for i in akeys + if i in bkeys + err = max(err, abs(a.states[i] - b.states[i])) + else + return 1.0e10 + end + end + akeysp = keys(a.parameters) + bkeysp = keys(b.parameters) + for j in akeysp + if j in bkeysp + err = max(err, abs(a.parameters[j] - b.parameters[j])) + else + return 1.0e10 + end + end + return err +end + +function new_clustering(estimates, tol = 1e-3) + new_estimates = [] + for e in estimates + keep = true + for f in new_estimates + if (((abs(e.err - f.err)) < tol) && (compare_estimation_results(e, f) < tol)) + keep = false + end + end + if (keep) + push!(new_estimates, e) + end + + end + return new_estimates +end + +function post_process(estimates, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing) + if Threads.nthreads() > 1 + estimates = filter(x -> !isnothing(x), estimates) + estimates = vcat(estimates...) + end + #filter out the empty vectors + estimates = filter(x -> length(x) > 0, estimates) + estimates = filter(x -> x[1].return_code == ReturnCode.Success, estimates) + + if (filtermode == :new) + estimates = collect(Iterators.flatten(estimates)) + estimates = filter(x -> check_constraints(x, parameter_constraints, ic_constraints), estimates) + estimates = sort(estimates, by = x -> x.err) + + estimates_clustered_by_params = new_clustering(estimates, 1e-2) #TODO(orebas): this is a magic number. + estimates_clustered_by_params = filter(x -> (x.err < 1.0), estimates_clustered_by_params) #also a magic number + return estimates_clustered_by_params + else + best_solution = nothing + for each in estimates + if all(!isnothing(x.err) for x in each) + if best_solution === nothing + best_solution = each + else + if sum(x.err for x in each) < sum(x.err for x in best_solution) + best_solution = each + end + end + else + best_solution = nothing + end + end + if best_solution !== nothing + @info "Best estimate found" + return best_solution + else + @warn "No solution found" + return estimates + end + end end diff --git a/src/filtering.jl b/src/filtering.jl index 4ab05b4a..2fdbd58d 100644 --- a/src/filtering.jl +++ b/src/filtering.jl @@ -1,6 +1,6 @@ """ - solve_ode(model, estimate::EstimationResult, data_sample; solver = Tsit5(), - return_ode = false) + solve_ode(model, estimate::EstimationResult, data_sample; solver = Tsit5(), + return_ode = false) Solves the ODE system `model` with the parameters and initial conditions given by `estimate`. Compute the error between the solution and the data sample. The error is recorded in the `EstimationResult`. @@ -9,8 +9,8 @@ Compute the error between the solution and the data sample. The error is recorde - `model`: the ODE system to be solved. - `estimate::EstimationResult`: the parameters and initial conditions of the ODE system. - `data_sample`: the data sample used for estimation (same functions as `measured_quantities`). - The keys of the dictionary are the measured quantities - and the values are the corresponding data samples. + The keys of the dictionary are the measured quantities + and the values are the corresponding data samples. - `solver = Tsit5()`: (optional) the solver used to solve the ODE system, see `DifferentialEquations` for available solvers. - `return_ode = false`: (optional) whether to return the ODE solution. @@ -19,67 +19,67 @@ Compute the error between the solution and the data sample. The error is recorde - `EstimationResult`: the estimated parameters and initial conditions of the model. """ function solve_ode(model, estimate::EstimationResult, inputs::Vector{Equation}, data_sample; - solver = Tsit5(), return_ode = false) - initial_conditions = [estimate[s] for s in ModelingToolkit.states(model)] - parameter_values = [estimate[p] for p in ModelingToolkit.parameters(model)] - tspan = (estimate.at_time, data_sample["t"][end] + estimate.at_time) - ode_equations = ModelingToolkit.equations(model) - ode_equations = substitute(ode_equations, - Dict(each.lhs => Num(each.rhs) for each in inputs)) - t = ModelingToolkit.get_iv(model) - @named new_model = ODESystem(ode_equations, t, ModelingToolkit.states(model), - ModelingToolkit.parameters(model)) - prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values) - ode_solution = ModelingToolkit.solve(prob, solver, - saveat = range(tspan[1], tspan[2], - length = length(data_sample["t"]))) - if ode_solution.retcode == ReturnCode.Success - err = 0 - for (key, sample) in data_sample - if isequal(key, "t") - continue - end - err += ParameterEstimation.mean_abs_err(ode_solution(data_sample["t"])[key], - sample) - end - err /= length(data_sample) - else - err = 1e+10 - end - if return_ode - return ode_solution, - EstimationResult(estimate.parameters, estimate.states, - estimate.degree, estimate.at_time, err, - estimate.interpolants, estimate.return_code, - estimate.datasize) - else - return EstimationResult(estimate.parameters, estimate.states, - estimate.degree, estimate.at_time, err, - estimate.interpolants, estimate.return_code, - estimate.datasize) - end + solver = Tsit5(), return_ode = false, abstol = 1e-12, reltol = 1e-12) + initial_conditions = [estimate[s] for s in ModelingToolkit.states(model)] + parameter_values = [estimate[p] for p in ModelingToolkit.parameters(model)] + tspan = (estimate.at_time, data_sample["t"][end] + estimate.at_time) + ode_equations = ModelingToolkit.equations(model) + ode_equations = substitute(ode_equations, + Dict(each.lhs => Num(each.rhs) for each in inputs)) + t = ModelingToolkit.get_iv(model) + @named new_model = ODESystem(ode_equations, t, ModelingToolkit.states(model), + ModelingToolkit.parameters(model)) + prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values) + ode_solution = ModelingToolkit.solve(prob, solver, + saveat = range(tspan[1], tspan[2], + length = length(data_sample["t"])), abstol = 1e-12, reltol = 1e-12) + if ode_solution.retcode == ReturnCode.Success + err = 0 + for (key, sample) in data_sample + if isequal(key, "t") + continue + end + err += ParameterEstimation.mean_abs_err(ode_solution(data_sample["t"])[key], + sample) + end + err /= length(data_sample) + else + err = 1e+10 + end + if return_ode + return ode_solution, + EstimationResult(estimate.parameters, estimate.states, + estimate.degree, estimate.at_time, err, + estimate.interpolants, estimate.return_code, + estimate.datasize, estimate.report_time) + else + return EstimationResult(estimate.parameters, estimate.states, + estimate.degree, estimate.at_time, err, + estimate.interpolants, estimate.return_code, + estimate.datasize, estimate.report_time) + end end """ - solve_ode!(model, estimates::Vector{EstimationResult}, inputs::Vector{Equation},data_sample; solver = Tsit5()) + solve_ode!(model, estimates::Vector{EstimationResult}, inputs::Vector{Equation},data_sample; solver = Tsit5()) Run solve_ode for multiple estimates and store the results (error between solution and sample) in each estimate. This is done in-place. """ function solve_ode!(model, estimates::Vector{EstimationResult}, inputs::Vector{Equation}, - data_sample; solver = Tsit5()) - estimates[:] = map(each -> solve_ode(model, each, inputs, data_sample, solver = solver), - estimates) + data_sample; solver = Tsit5(), abstol = 1e-12, reltol = 1e-12) + estimates[:] = map(each -> solve_ode(model, each, inputs, data_sample, solver = solver, abstol = abstol, reltol = reltol), + estimates) end """ - filter_solutions(results::Vector{EstimationResult}, - identifiability_result::IdentifiabilityData, - model::ModelingToolkit.ODESystem, - inputs::Vector{ModelingToolkit.Equation}, - data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); - solver = Tsit5(), - topk = 1) where {T <: Float} + filter_solutions(results::Vector{EstimationResult}, + identifiability_result::IdentifiabilityData, + model::ModelingToolkit.ODESystem, + inputs::Vector{ModelingToolkit.Equation}, + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + solver = Tsit5(), + topk = 1) where {T <: Float} Filter estimation results stored in `results` vector based on ODE solving and checking against the sample. In addition, takes into account global and local identifiability of parameters when filtering. @@ -90,8 +90,8 @@ In addition, takes into account global and local identifiability of parameters w - `model::ModelingToolkit.ODESystem`: the ODE system. - `inputs::Vector{ModelingToolkit.Equation}`: the inputs of the ODE system. - `data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}()`: the data sample used for estimation (same functions as `measured_quantities`). - The keys of the dictionary are the measured quantities - and the values are the corresponding data samples. + The keys of the dictionary are the measured quantities + and the values are the corresponding data samples. - `time_interval::Vector{T} = Vector{T}()`: the time interval of the ODE system. - `topk = 1`: (optional) the number of best estimates to return. @@ -99,96 +99,104 @@ In addition, takes into account global and local identifiability of parameters w - `EstimationResult`: the best estimate (if `topk = 1`) or the vector of best estimates (if `topk > 1`). """ function filter_solutions(results::Vector{EstimationResult}, - identifiability_result::IdentifiabilityData, - model::ModelingToolkit.ODESystem, - inputs::Vector{ModelingToolkit.Equation}, - data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); - solver = Tsit5(), - topk = 1) where {T <: Float} - @debug "Filtering" - if length(results) == 0 - @warn "No results to filter." - return results - end - if all(each -> each.return_code == ReturnCode.Failure, results) - return results - end - try - solve_ode!(model, results, inputs, data_sample; solver = solver) # this solves ODE with new parameters and computes err. between sample and solution - catch InexactError - @debug "InexactError when solving the ODE, no filtering was done." - return results - end - if length(identifiability_result["identifiability"]["nonidentifiable"]) > 0 - @debug "The model contains non-identifiable parameters" - filtered_results = Vector{ParameterEstimation.EstimationResult}() - clustered = ParameterEstimation.cluster_estimates(model, results, inputs, - data_sample, solver = solver) - if length(clustered) == 0 - @debug "No clustered data, cannot find results to filter." - return results - end - min_cluster_err, min_cluster_idx = findmin(sum(each.err for each in group) / - length(group) - for (id, group) in pairs(clustered)) - min_cluster = clustered[min_cluster_idx] - @debug "Best estimate yelds ODE solution error $(min_cluster_err)" - filtered_results = min_cluster - return results - end - if length(identifiability_result["identifiability"]["locally_not_globally"]) > 0 - filtered_results = Vector{ParameterEstimation.EstimationResult}() - clustered = ParameterEstimation.cluster_estimates(model, results, inputs, - data_sample, solver = solver) - if length(clustered) == 0 - @debug "No clustered data, cannot find results to filter." - return results - end - # find cluster with smallest error - min_cluster_err, min_cluster_idx = findmin(sum(each.err for each in group) / - length(group) - for (id, group) in pairs(clustered)) - min_cluster = clustered[min_cluster_idx] - @debug "Best estimate yelds ODE solution error $(min_cluster_err)" - filtered_results = min_cluster - else - sorted = sort(results, by = x -> x.err) - if topk == 1 - filtered_results = Vector{EstimationResult}() - @debug "Best estimate yelds ODE solution error $(sorted[1].err)" - push!(filtered_results, sorted[1]) - else - filtered_results = Vector{Vector{ParameterEstimation.EstimationResult}}() - @debug "Best $(topk) estimates yeld ODE solution errors $([s.err for s in sorted[1:topk]])" - push!(filtered_results, sorted[1:topk]) - end - end - return filtered_results -end + identifiability_result::IdentifiabilityData, + model::ModelingToolkit.ODESystem, + inputs::Vector{ModelingToolkit.Equation}, + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(); + solver = Tsit5(), + topk = 1, filtermode = :new, abstol = 1e-12, reltol = 1e-12) where {T <: Float} + @debug "Filtering" + #if (filtermode == :new) + # println("new mode filtering") + #end + if length(results) == 0 + @warn "No results to filter." + return results + end + if all(each -> each.return_code == ReturnCode.Failure, results) + return results + end + try + solve_ode!(model, results, inputs, data_sample; solver = solver, abstol = abstol, reltol = reltol) # this solves ODE with new parameters and computes err. between sample and solution + catch InexactError + println("Inexact error") + @debug "InexactError when solving the ODE, no filtering was done." + return results + end + if (filtermode == :new) + return results + else + if length(identifiability_result["identifiability"]["nonidentifiable"]) > 0 + @debug "The model contains non-identifiable parameters" + filtered_results = Vector{ParameterEstimation.EstimationResult}() + clustered = ParameterEstimation.cluster_estimates(model, results, inputs, + data_sample, solver = solver) + if length(clustered) == 0 + @debug "No clustered data, cannot find results to filter." + return results + end + min_cluster_err, min_cluster_idx = findmin(sum(each.err for each in group) / + length(group) + for (id, group) in pairs(clustered)) + min_cluster = clustered[min_cluster_idx] + @debug "Best estimate yelds ODE solution error $(min_cluster_err)" + filtered_results = min_cluster + return results + end + if length(identifiability_result["identifiability"]["locally_not_globally"]) > 0 + filtered_results = Vector{ParameterEstimation.EstimationResult}() + clustered = ParameterEstimation.cluster_estimates(model, results, inputs, + data_sample, solver = solver) + if length(clustered) == 0 + @debug "No clustered data, cannot find results to filter." + return results + end + # find cluster with smallest error + min_cluster_err, min_cluster_idx = findmin(sum(each.err for each in group) / + length(group) + for (id, group) in pairs(clustered)) + min_cluster = clustered[min_cluster_idx] + @debug "Best estimate yelds ODE solution error $(min_cluster_err)" + filtered_results = min_cluster + else + sorted = sort(results, by = x -> x.err) + if topk == 1 + filtered_results = Vector{EstimationResult}() + @debug "Best estimate yelds ODE solution error $(sorted[1].err)" + push!(filtered_results, sorted[1]) + else + filtered_results = Vector{Vector{ParameterEstimation.EstimationResult}}() + @debug "Best $(topk) estimates yeld ODE solution errors $([s.err for s in sorted[1:topk]])" + push!(filtered_results, sorted[1:topk]) + end + end + return filtered_results + end +end """ - function cluster_estimates(model, res, data_sample; ε = 1e-6, solver = Tsit5()) + function cluster_estimates(model, res, data_sample; ε = 1e-6, solver = Tsit5()) - Clusters the estimates by their error. Used in case of non-identifiability of parameters. + Clusters the estimates by their error. Used in case of non-identifiability of parameters. """ function cluster_estimates(model, res, inputs::Vector{Equation}, data_sample; ε = 1e-6, - solver = Tsit5()) - # clusrers the estimates by their error - ParameterEstimation.solve_ode!(model, res, inputs, data_sample, solver = solver) - clustered = Dict() - #nearest neighbor search by err - for i in eachindex(res) - for j in (i + 1):length(res) - if abs(res[i].err - res[j].err) < ε - if !haskey(clustered, i) - clustered[i] = Vector{ParameterEstimation.EstimationResult}([]) - push!(clustered[i], res[i]) - end - push!(clustered[i], res[j]) - end - end - end - # reset cluster keys - clustered = Dict(i => clustered[key] for (i, key) in enumerate(keys(clustered))) - return clustered + solver = Tsit5()) + # clusrers the estimates by their error + ParameterEstimation.solve_ode!(model, res, inputs, data_sample, solver = solver) + clustered = Dict() + #nearest neighbor search by err + for i in eachindex(res) + for j in (i+1):length(res) + if abs(res[i].err - res[j].err) < ε + if !haskey(clustered, i) + clustered[i] = Vector{ParameterEstimation.EstimationResult}([]) + push!(clustered[i], res[i]) + end + push!(clustered[i], res[j]) + end + end + end + # reset cluster keys + clustered = Dict(i => clustered[key] for (i, key) in enumerate(keys(clustered))) + return clustered end diff --git a/src/identifiability/check_identifiability.jl b/src/identifiability/check_identifiability.jl index d2506dc1..f67b86bc 100644 --- a/src/identifiability/check_identifiability.jl +++ b/src/identifiability/check_identifiability.jl @@ -24,7 +24,7 @@ object that contains the results of the identifiability analysis. [3] - https://github.com/alexeyovchinnikov/SIAN-Julia """ -function check_identifiability(ode::ModelingToolkit.ODESystem; + function check_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities::Vector{Equation} = Array{Equation}[], inputs::Vector{Num} = Array{Num}[], infolevel = 0) diff --git a/src/includes.jl b/src/includes.jl index 94e0d938..7e8fe2c2 100644 --- a/src/includes.jl +++ b/src/includes.jl @@ -1,6 +1,7 @@ include("rational_interpolation/rational_interpolation.jl") include("rational_interpolation/interpolant.jl") include("rational_interpolation/utils.jl") +include("rational_interpolation/bary_derivs.jl") include("identifiability/check_identifiability.jl") include("identifiability/preprocessor.jl") diff --git a/src/rational_interpolation/bary_derivs.jl b/src/rational_interpolation/bary_derivs.jl new file mode 100644 index 00000000..0bf4a1e3 --- /dev/null +++ b/src/rational_interpolation/bary_derivs.jl @@ -0,0 +1,301 @@ +""" + rational_interpolation_coefficients(x, y, n) +CODE COPIED FROM previous version of ParameterEstimation.jl +Perform a rational interpolation of the data `y` at the points `x` with numerator degree `n`. +This function only returns the coefficients of the numerator and denominator polynomials. + +# Arguments +- `x`: the points where the data is sampled (e.g. time points). +- `y`: the data sample. +- `n`: the degree of the numerator. + +# Returns +- `c`: the coefficients of the numerator polynomial. +- `d`: the coefficients of the denominator polynomial. +""" +function rational_interpolation_coefficients(x, y, n) + N = length(x) + m = N - n - 1 + A = zeros(N, N) + if m > 0 + A_left_submatrix = reduce(hcat, [x .^ (i) for i in 0:(n)]) + A_right_submatrix = reduce(hcat, [x .^ (i) for i in 0:(m-1)]) + A = hcat(A_left_submatrix, -y .* A_right_submatrix) + b = y .* (x .^ m) + try + prob = LinearSolve.LinearProblem(A, b) + c = LinearSolve.solve(prob) + return c[1:(n+1)], [c[(n+2):end]; 1] + catch SingularException + lu_res = lu(A) + y = lu_res.L \ lu_res.P * b + c = lu_res.U \ y + return c[1:(n+1)], [c[(n+2):end]; 1] + end + + else + A = reduce(hcat, [x .^ i for i in 0:n]) + b = y + prob = LinearSolve.LinearProblem(A, b) + c = LinearSolve.solve(prob) + return c, [1] + end +end + + + + + +######### TODO(orebas) REFACTOR into bary_derivs or something similiar +#to use the below, you can just pass vectors of xvalues and yvalues like so: +#F = aaad(xdata, ydata) +#and then F will be a callable, i.e. F(0.5) should work. Let's please restrict to real xdata, and if so it should be sorted. F is only defined in the range of the xvalues. + +#To construct a derivative, you can do +#derivf(z) = ForwardDiff.derivative(F, z) +#I hope and suspect that other AD frameworks should work as well. + +function baryEval(z, f::Vector{T}, x::Vector{T}, w::Vector{T}, tol = 1e-13) where {T} + @assert(length(f) == length(x) == length(w)) + num = zero(T) + den = zero(T) + breakflag = false + breakindex = -1 + for j in eachindex(f) + t = w[j] / (z - x[j]) + num += t * f[j] + den += t + if ((z - x[j])^2 < sqrt(tol)) + breakflag = true + breakindex = j + end + end + fz = num / den + if (breakflag) + num = zero(T) + den = zero(T) + for j in eachindex(f) + if (j != breakindex) + t = w[j] / (z - x[j]) + num += t * f[j] + den += t + end + end + m = z - x[breakindex] + fz = (w[breakindex] * f[breakindex] + m * num) / (w[breakindex] + m * den) + end + return fz +end + +struct AAADapprox{T} + internalAAA::T +end + +struct FHDapprox{T} + internalFHD::T +end + +(y::FHDapprox)(z) = baryEval(z, y.internalFHD.f, y.internalFHD.x, y.internalFHD.w) +(y::AAADapprox)(z) = baryEval(z, y.internalAAA.f, y.internalAAA.x, y.internalAAA.w) + +function nth_deriv_at_deprecated(f, n::Int, t) #todo(orebas) make this more efficient. + if (n == 0) + return f(t) + elseif (n == 1) + return ForwardDiff.derivative(f, t) + else + g(t) = nth_deriv_at(f, n - 1, t) + return ForwardDiff.derivative(g, t) + end +end + + +function nth_deriv_at(f, n::Int, t) #todo(orebas) make this more efficient. + if (n == 0) + return f(t) + else + return TaylorDiff.derivative(f, t, n) + end +end + + +function aaad(xs::AbstractArray{T}, ys::AbstractArray{T}) where {T} + @suppress begin + @assert length(xs) == length(ys) + internalApprox = BaryRational.aaa(xs, ys) + return AAADapprox(internalApprox) + end +end + + + +function fhd(xs::AbstractArray{T}, ys::AbstractArray{T}, N::Int) where {T} + @suppress begin + @assert length(xs) == length(ys) + internalApprox = BaryRational.FHInterp(xs, ys, order = N) + return FHDapprox(internalApprox) + end +end + +function fhdn(n) + fh(xs, ys) = fhd(xs, ys, n) + return fh +end + + +#################### + +struct FourierSeries + m::Any + b::Any + K::Any + cosines::Any + sines::Any +end + + +function fourierEval(x, FS) + z = FS.m * x + FS.b + sum = 0.0 + for k in eachindex(FS.cosines) + sum += FS.cosines[k] * cos((k) * z) + end + for k in eachindex(FS.sines) + sum += FS.sines[k] * sin((k) * z) + end + sum += FS.K + return sum +end + +(y::FourierSeries)(z) = fourierEval(z, y) + + + +function FourierInterp(xs, ys) + @assert length(xs) == length(ys) + N = length(xs) + width = xs[end] - xs[begin] + m = pi / width + b = -pi * (xs[begin] / width + 0.5) + f(t) = m * t + b + scaledxs = f.(xs) + sinescount = (N - 1) ÷ 2 + cosinescount = N - 1 - sinescount + A = zeros(Float64, N, N) + for i ∈ 1:N, j ∈ 1:N + if (j == 1) + A[i, 1] = 1 + elseif (j <= (cosinescount + 1)) + A[i, j] = cos((j - 1) * scaledxs[i]) + else + temp = (j - cosinescount - 1) + A[i, j] = sin(temp * scaledxs[i]) + end + end + + prob = LinearProblem(A, ys, LinearSolve.OperatorCondition.VeryIllConditioned) + sol = LinearSolve.solve(prob) + X = sol.u + temp = FourierSeries(m, b, X[begin], X[2:(cosinescount+1)], X[(cosinescount+2):end]) + return temp +end + + + +struct BaryLagrange{T <: AbstractArray} + x::T + f::T + w::T +end + +function BarycentricLagrange(xs, ys) + @assert length(xs) == length(ys) + N = length(xs) + w = ones(Float64, N) + for k in eachindex(xs) + for j in eachindex(xs) + if (k != j) + w[k] *= (xs[k] - xs[j]) + end + end + w[k] = 1 / w[k] + end + return BaryLagrange(xs, ys, w) +end + +(y::BaryLagrange)(z) = baryEval(z, y.f, y.x, y.w) + +struct RationalFunction{T <: AbstractArray} + a::T + b::T +end + +(y::RationalFunction)(z) = rationaleval(z, y.a, y.b) +rationaleval(z, a, b) = evalpoly(z, a) / evalpoly(z, b) + + + +function simpleratinterp(xs, ys, d1::Int) + @assert length(xs) == length(ys) + N = length(xs) + d2 = N - d1 - 1 + A = zeros(Float64, N, N) + for j in 1:N + A[j, 1] = 1 + for k in 1:d1 + A[j, k+1] = xs[j]^k + end + for k in 1:d2 + A[j, d1+1+k] = -1.0 * ys[j] * xs[j]^k + end + end + prob = LinearProblem(A, ys, LinearSolve.OperatorCondition.VeryIllConditioned) + sol = LinearSolve.solve(prob) + X = sol.u + return RationalFunction( + X[1:(d1+1)], + [1; X[(d1+2):end]]) +end + + + + + + +function betterratinterp(xs, ys, d1::Int) + @assert length(xs) == length(ys) + N = length(xs) + (c, d) = rational_interpolation_coefficients(xs, ys, d1) + return RationalFunction(c, d) +end + + + +function SimpleRationalInterp(numerator_degree::Int) + f(xs, ys) = simpleratinterp(xs, ys, numerator_degree) + return f +end + +function SimpleRationalInterpOld(numerator_degree::Int) + f(xs, ys) = betterratinterp(xs, ys, numerator_degree) + return f +end + + +function default_interpolator(datasize) + interpolators = Dict( + "AAA" => aaad, + "FHD3" => fhdn(3), + # "Fourier" => FourierInterp, + # "BaryLagrange" => BarycentricLagrange, + ) + if (datasize > 10) + interpolators["FHD8"] = fhdn(8) + # interpolators["FHD6"] = fhdn(6) + #stepsize = max(1, datasize ÷ 4) + #for i in range(1, (datasize - 2), step = stepsize) + # interpolators["RatOld($i)"] = SimpleRationalInterpOld(i) + #end + end +return interpolators +end diff --git a/src/rational_interpolation/interpolant.jl b/src/rational_interpolation/interpolant.jl index 770d4dbd..3fdeb01c 100644 --- a/src/rational_interpolation/interpolant.jl +++ b/src/rational_interpolation/interpolant.jl @@ -1,5 +1,5 @@ """ - Interpolant + Interpolant A structure that stores information about the interpolation result. @@ -8,6 +8,5 @@ A structure that stores information about the interpolation result. - `dIdt`: the `TaylorSeries` derivative of `I`. """ struct Interpolant - I::Any - dIdt::Any + f::Any end diff --git a/src/rational_interpolation/rational_interpolation.jl b/src/rational_interpolation/rational_interpolation.jl index 63c9a6c5..499386b0 100644 --- a/src/rational_interpolation/rational_interpolation.jl +++ b/src/rational_interpolation/rational_interpolation.jl @@ -1,57 +1,15 @@ -""" - rational_interpolation_coefficients(x, y, n) - -Perform a rational interpolation of the data `y` at the points `x` with numerator degree `n`. -This function only returns the coefficients of the numerator and denominator polynomials. -# Arguments -- `x`: the points where the data is sampled (e.g. time points). -- `y`: the data sample. -- `n`: the degree of the numerator. -# Returns -- `c`: the coefficients of the numerator polynomial. -- `d`: the coefficients of the denominator polynomial. """ -function rational_interpolation_coefficients(x, y, n) - N = length(x) - m = N - n - 1 - A = zeros(N, N) - if m > 0 - A_left_submatrix = reduce(hcat, [x .^ (i) for i in 0:(n)]) - A_right_submatrix = reduce(hcat, [x .^ (i) for i in 0:(m - 1)]) - A = hcat(A_left_submatrix, -y .* A_right_submatrix) - b = y .* (x .^ m) - try - prob = LinearSolve.LinearProblem(A, b) - c = LinearSolve.solve(prob) - return c[1:(n + 1)], [c[(n + 2):end]; 1] - catch SingularException - lu_res = lu(A) - y = lu_res.L \ lu_res.P * b - c = lu_res.U \ y - return c[1:(n + 1)], [c[(n + 2):end]; 1] - end - - else - A = reduce(hcat, [x .^ i for i in 0:n]) - b = y - prob = LinearSolve.LinearProblem(A, b) - c = LinearSolve.solve(prob) - return c, [1] - end -end - -""" - interpolate(identifiability_result, data_sample, - measured_quantities; interpolation_degree::Int = 1, - diff_order::Int = 1, at_t::Float = 0.0, - method::Symbol = :homotopy) + interpolate(identifiability_result, data_sample, + measured_quantities; interpolation_degree::Int = 1, + diff_order::Int = 1, at_t::Float = 0.0, + method::Symbol = :homotopy) This function performs the key step in parameter estimation. - It interpolates the data in `data_sample` and computes the `TaylorSeries` expansion. - These results are stored in the `Interpolant` object and are applied to the polynomial system in `identifiability_result`. + It interpolates the data in `data_sample` and computes the `TaylorSeries` expansion. + These results are stored in the `Interpolant` object and are applied to the polynomial system in `identifiability_result`. # Arguments - `identifiability_result`: the result of the identifiability check. @@ -66,62 +24,50 @@ This function performs the key step in parameter estimation. - `System`: the polynomial system with the interpolated data applied. This system is compatible with `HomotopyContinuation` solving. """ function interpolate(identifiability_result, data_sample, - measured_quantities, inputs; interpolation_degree::Int = 1, - diff_order::Int = 1, at_t::Float = 0.0, - method::Symbol = :homotopy) - polynomial_system = identifiability_result["polynomial_system"] - interpolants = Dict{Any, Interpolant}() - sampling_times = data_sample["t"] - for (key, sample) in pairs(data_sample) - if key == "t" - continue - end - y_function_name = map(x -> replace(string(x.lhs), "(t)" => ""), - filter(x -> string(x.rhs) == string(key), - measured_quantities))[1] - interpolant = ParameterEstimation.interpolate(sampling_times, sample, - interpolation_degree, - diff_order, at_t) - interpolants[key] = interpolant - err = sum(abs.(sample - interpolant.I.(sampling_times))) / length(sampling_times) - @debug "Mean Absolute error in interpolation: $err interpolating $key" - polynomial_system = eval_derivs(polynomial_system, interpolant, y_function_name, - inputs, identifiability_result, method = method) - end - if isequal(method, :homotopy) - try - identifiability_result["polynomial_system_to_solve"] = HomotopyContinuation.System(polynomial_system) - catch KeyError - throw(ArgumentError("HomotopyContinuation threw a KeyError, it is possible that " * - "you are using Unicode characters in your input. Consider " * - "using ASCII characters instead.")) - end - else - identifiability_result["polynomial_system_to_solve"] = polynomial_system - end - return interpolants + measured_quantities, inputs; interpolator, + diff_order::Int = 1, at_t::Float = 0.0, #TODO(orebas)should we remove diff_order? + method::Symbol = :homotopy) + polynomial_system = identifiability_result["polynomial_system"] + interpolants = Dict{Any, Interpolant}() + sampling_times = data_sample["t"] + for (key, sample) in pairs(data_sample) + if key == "t" + continue + end + y_function_name = map(x -> replace(string(x.lhs), "(t)" => ""), + filter(x -> string(x.rhs) == string(key), + measured_quantities))[1] + interpolant = ParameterEstimation.interpolate(sampling_times, sample, + interpolator, + diff_order) + interpolants[key] = interpolant + err = sum(abs.(sample - interpolant.f.(sampling_times))) / length(sampling_times) + @debug "Mean Absolute error in interpolation: $err interpolating $key" + polynomial_system = eval_derivs(polynomial_system, interpolant, y_function_name, + inputs, identifiability_result, at_time = at_t, method = method) + end + if isequal(method, :homotopy) + try + identifiability_result["polynomial_system_to_solve"] = HomotopyContinuation.System(polynomial_system) + catch KeyError + throw(ArgumentError("HomotopyContinuation threw a KeyError, it is possible that " * + "you are using Unicode characters in your input. Consider " * + "using ASCII characters instead.")) + end + else + identifiability_result["polynomial_system_to_solve"] = polynomial_system + end + return interpolants end """ - interpolate(time, sample, numer_degree::Int, diff_order::Int = 1, at_t::Float = 0.0) + interpolate(time, sample, numer_degree::Int, diff_order::Int = 1, at_t::Float = 0.0) This function performs a rational interpolation of the data `sample` at the points `time` with numerator degree `numer_degree`. It returns an `Interpolant` object that contains the interpolated function and its derivatives. """ -function interpolate(time, sample, numer_degree::Int, diff_order::Int = 1, - at_t::Float = 0.0) - numer_coef, denom_coef = rational_interpolation_coefficients(time, sample, - numer_degree) - numer_function(t) = sum(numer_coef[i] * t^(i - 1) for i in eachindex(numer_coef)) - denom_function(t) = sum(denom_coef[i] * t^(i - 1) for i in eachindex(denom_coef)) - interpolated_function(t) = numer_function(t) / denom_function(t) - return Interpolant(interpolated_function, - differentiate_interpolated(interpolated_function, diff_order, at_t)) +function interpolate(time, sample, interpolator, diff_order::Int = 1) + interpolated_function = ((interpolator.second))(time, sample) + return Interpolant(interpolated_function) end -function differentiate_interpolated(interpolated_function, diff_order::Int, - at_t::Float = 0.0) - τ = Taylor1(diff_order + 1) - taylor_expantion = interpolated_function(τ - at_t) - return taylor_expantion -end diff --git a/src/rational_interpolation/utils.jl b/src/rational_interpolation/utils.jl index 9cca134b..fbf42378 100644 --- a/src/rational_interpolation/utils.jl +++ b/src/rational_interpolation/utils.jl @@ -1,69 +1,69 @@ """ - function eval_derivs(polynomial_system, interpolant::Interpolant, - y_function_name, - inputs::Vector{ModelingToolkit.Equation}, - identifiability_result; - method = :homotopy) + function eval_derivs(polynomial_system, interpolant::Interpolant, + y_function_name, + inputs::Vector{ModelingToolkit.Equation}, + identifiability_result; + method = :homotopy) This function evaluates the derivatives of the interpolated function `y_function_name` using the `interpolant` object. - Derivatives are substituted into the polynomial system. + Derivatives are substituted into the polynomial system. """ function eval_derivs(polynomial_system, interpolant::Interpolant, - y_function_name, - inputs::Vector{ModelingToolkit.Equation}, - identifiability_result; - method = :homotopy) - if isequal(method, :homotopy) - for (y_func, y_deriv_order) in pairs(identifiability_result["Y_eq"]) - if occursin(y_function_name, string(y_func)) - y_derivs_vals = Dict(ParameterEstimation.nemo2hc(y_func) => interpolant.dIdt[y_deriv_order] * - factorial(y_deriv_order)) - polynomial_system = HomotopyContinuation.evaluate(ParameterEstimation.nemo2hc.(polynomial_system), - y_derivs_vals) - end - end - for (u_funct, u_deriv_order) in pairs(identifiability_result["u_variables"]) - t = arguments(inputs[1].lhs)[1] - tau = Taylor1(u_deriv_order) - for input_eq in inputs - u_function_name = replace(string(input_eq.lhs), "(t)" => "") - if occursin(u_function_name, string(u_funct)) - u_derivs_vals = Dict(ParameterEstimation.nemo2hc(u_funct) => substitute(input_eq.rhs, - Dict(t => tau))[u_deriv_order] * - factorial(u_deriv_order)) - polynomial_system = HomotopyContinuation.evaluate(ParameterEstimation.nemo2hc.(polynomial_system), - u_derivs_vals) - end - end - end - elseif isequal(method, :msolve) - y_derivs = Vector{SIAN.Nemo.fmpq_mpoly}() - y_vals = Vector{SIAN.Nemo.fmpq}() - u_derivs = Vector{SIAN.Nemo.fmpq_mpoly}() - u_vals = Vector{SIAN.Nemo.fmpq}() - for (y_func, y_deriv_order) in pairs(identifiability_result["Y_eq"]) - if occursin(y_function_name, string(y_func)) - push!(y_derivs, y_func) - push!(y_vals, - rationalize(Float64(interpolant.dIdt[y_deriv_order] * - factorial(y_deriv_order)))) - end - end - for (u_funct, u_deriv_order) in pairs(identifiability_result["u_variables"]) - tau = Taylor1(u_deriv_order) - for input_eq in inputs - u_function_name = replace(string(input_eq.lhs), "(t)" => "") - if occursin(u_function_name, string(u_funct)) - push!(u_derivs, u_funct) - push!(u_vals, - rationalize(Float64(substitute(input_eq.rhs, - Dict(t => tau))[u_deriv_order] * - factorial(u_deriv_order)))) - end - end - end - polynomial_system = [SIAN.Nemo.evaluate(poly, y_derivs, y_vals) - for poly in polynomial_system] - end - return polynomial_system + y_function_name, + inputs::Vector{ModelingToolkit.Equation}, + identifiability_result; + at_time = 0.0, + method = :homotopy) + if isequal(method, :homotopy) + for (y_func, y_deriv_order) in pairs(identifiability_result["Y_eq"]) + if occursin(y_function_name, string(y_func)) + y_derivs_vals = Dict(ParameterEstimation.nemo2hc(y_func) => (nth_deriv_at(interpolant.f, y_deriv_order, at_time))) #check for off-by-one in derivb o + #println(y_derivs_vals) #DEBUG + polynomial_system = HomotopyContinuation.evaluate(ParameterEstimation.nemo2hc.(polynomial_system), + y_derivs_vals) + end + end + for (u_funct, u_deriv_order) in pairs(identifiability_result["u_variables"]) + t = arguments(inputs[1].lhs)[1] + tau = Taylor1(u_deriv_order) #TODO(orebas) replace taylor series with forwardDiff, perhaps? + for input_eq in inputs + u_function_name = replace(string(input_eq.lhs), "(t)" => "") + if occursin(u_function_name, string(u_funct)) + u_derivs_vals = Dict(ParameterEstimation.nemo2hc(u_funct) => substitute(input_eq.rhs, + Dict(t => tau))[u_deriv_order] * + factorial(u_deriv_order)) + polynomial_system = HomotopyContinuation.evaluate(ParameterEstimation.nemo2hc.(polynomial_system), + u_derivs_vals) + end + end + end + elseif isequal(method, :msolve) + y_derivs = Vector{SIAN.Nemo.fmpq_mpoly}() + y_vals = Vector{SIAN.Nemo.fmpq}() + u_derivs = Vector{SIAN.Nemo.fmpq_mpoly}() + u_vals = Vector{SIAN.Nemo.fmpq}() + for (y_func, y_deriv_order) in pairs(identifiability_result["Y_eq"]) + if occursin(y_function_name, string(y_func)) + push!(y_derivs, y_func) + push!(y_vals, + rationalize(Float64(nth_deriv_at(interpolant.f, y_deriv_order, at_time)))) + end + end + for (u_funct, u_deriv_order) in pairs(identifiability_result["u_variables"]) + tau = Taylor1(u_deriv_order) + for input_eq in inputs + u_function_name = replace(string(input_eq.lhs), "(t)" => "") + if occursin(u_function_name, string(u_funct)) + push!(u_derivs, u_funct) + push!(u_vals, + rationalize(Float64(substitute(input_eq.rhs, + Dict(t => tau))[u_deriv_order] * + factorial(u_deriv_order)))) + end + end + end + polynomial_system = [SIAN.Nemo.evaluate(poly, y_derivs, y_vals) + for poly in polynomial_system] + end + return polynomial_system end diff --git a/src/result.jl b/src/result.jl index a214d696..9f3f3c8f 100644 --- a/src/result.jl +++ b/src/result.jl @@ -1,5 +1,5 @@ """ - EstimationResult + EstimationResult A container for the results of an estimation. Contains the estimated parameters and initial conditions (state values at a given time), the degree of the rational interpolation used, @@ -14,78 +14,93 @@ the error between the estimated ODE solution and the sample data, and the return - `interpolants::Union{Nothing, Dict{Any, Interpolant}}`: The rational interpolants used to estimate the parameters and initial conditions. - `return_code::Any`: The return code of the estimation. - `datasize::Int64`: The number of data points used in the estimation. +- `report_time::Any`: The time at which the initial conditions are reported (usually the leftmost point in the time span). + """ struct EstimationResult - parameters::AbstractDict - states::AbstractDict - degree::Int64 - at_time::Float64 - err::Union{Nothing, Float64} - interpolants::Union{Nothing, AbstractDict{Any, Interpolant}} - return_code::Any - datasize::Int64 - function EstimationResult(model::ModelingToolkit.ODESystem, - poly_sol::AbstractDict, degree::Int64, - at_time::Float64, - interpolants::AbstractDict{Any, Interpolant}, - return_code, datasize) - parameters = OrderedDict{Any, Any}() - states = OrderedDict{Any, Any}() - for p in ModelingToolkit.parameters(model) - parameters[ModelingToolkit.Num(p)] = get(poly_sol, p, nothing) - end - for s in ModelingToolkit.states(model) - states[ModelingToolkit.Num(s)] = get(poly_sol, s, nothing) - end - new(parameters, states, degree, at_time, nothing, interpolants, return_code, - datasize) - end - function EstimationResult(parameters::AbstractDict, states::AbstractDict, degree::Int64, - at_time::Float64, err, interpolants, return_code, datasize) - new(parameters, states, degree, at_time, err, interpolants, return_code, datasize) - end + parameters::AbstractDict + states::AbstractDict + degree::Any + at_time::Float64 + err::Union{Nothing, Float64} + interpolants::Union{Nothing, AbstractDict{Any, Interpolant}} + return_code::Any + datasize::Int64 + report_time::Any + function EstimationResult(model::ModelingToolkit.ODESystem, + poly_sol::AbstractDict, degree, + at_time::Float64, + interpolants::AbstractDict{Any, Interpolant}, + return_code, datasize, report_time) + parameters = OrderedDict{Any, Any}() + states = OrderedDict{Any, Any}() + for p in ModelingToolkit.parameters(model) + parameters[ModelingToolkit.Num(p)] = get(poly_sol, p, nothing) + end + for s in ModelingToolkit.states(model) + states[ModelingToolkit.Num(s)] = get(poly_sol, s, nothing) + end + new(parameters, states, degree, at_time, nothing, interpolants, return_code, + datasize, report_time) + end + function EstimationResult(parameters::AbstractDict, states::AbstractDict, degree, + at_time::Float64, err, interpolants, return_code, datasize, report_time::Float64) + new(parameters, states, degree, at_time, err, interpolants, return_code, datasize, report_time) + end end Base.get(sol::EstimationResult, k) = get(sol.parameters, k, get(sol.states, k, nothing)) function Base.getindex(sol::EstimationResult, k) - if haskey(sol.parameters, k) - return sol.parameters[k] - elseif haskey(sol.states, k) - return sol.states[k] - else - throw(KeyError(k)) - end + if haskey(sol.parameters, k) + return sol.parameters[k] + elseif haskey(sol.states, k) + return sol.states[k] + else + throw(KeyError(k)) + end end function Base.show(io::IO, e::EstimationResult) - if any(isnothing.(values(e.parameters))) - println(io, "Parameter(s) :\t", - join([@sprintf("%3s = %3s", k, v) for (k, v) in pairs(e.parameters)], - ", ")) - println(io, "Initial Condition(s):\t", - join([@sprintf("%3s = %3s", k, v) for (k, v) in pairs(e.states)], ", ")) - elseif !all(isreal.(values(e.parameters))) - println(io, "Parameter(s) :\t", - join([@sprintf("%3s = %.3f+%.3fim", k, real(v), imag(v)) - for (k, v) in pairs(e.parameters)], - ", ")) - println(io, "Initial Condition(s):\t", - join([@sprintf("%3s = %.3f+%.3fim", k, real(v), imag(v)) - for (k, v) in pairs(e.states)], ", ")) - else - println(io, "Parameter(s) :\t", - join([@sprintf("%3s = %.3f", k, v) for (k, v) in pairs(e.parameters)], - ", ")) - println(io, "Initial Condition(s):\t", - join([@sprintf("%3s = %.3f", k, v) for (k, v) in pairs(e.states)], ", ")) - end - # println(io, "Interpolation Degree (numerator): ", e.degree) - # println(io, "Interpolation Degree (denominator): ", e.datasize - e.degree - 1) - # if isnothing(e.err) - # println(io, "Error: Not yet calculated") - # else - # println(io, "Error: ", @sprintf("%.4e", e.err)) - # end - # println(io, "Return Code: ", e.return_code) + if (!isnothing(e.report_time)) + report_time_string = @sprintf(", where t = %.3f", e.report_time) + else + report_time_string = "" + end + + if any(isnothing.(values(e.parameters))) + println(io, "Parameter(s) :\t", + join([@sprintf("%3s = %3s", k, v) for (k, v) in pairs(e.parameters)], + ", ")) + println(io, "Initial Condition(s):\t", + join([@sprintf("%3s = %3s", k, v) for (k, v) in pairs(e.states)], ", "), report_time_string) + elseif !all(isreal.(values(e.parameters))) + println(io, "Parameter(s) :\t", + join([@sprintf("%3s = %.3f+%.3fim", k, real(v), imag(v)) + for (k, v) in pairs(e.parameters)], + ", ")) + println(io, "Initial Condition(s):\t", + join([@sprintf("%3s = %.3f+%.3fim", k, real(v), imag(v)) + for (k, v) in pairs(e.states)], ", "), report_time_string) + else + println(io, "Parameter(s) :\t", + join([@sprintf("%3s = %.3f", k, v) for (k, v) in pairs(e.parameters)], + ", ")) + println(io, "Initial Condition(s):\t", + join([@sprintf("%3s = %.3f", k, v) for (k, v) in pairs(e.states)], ", "), report_time_string) + end + # println(io, "Interpolation Degree (numerator): ", e.degree) + # println(io, "Interpolation Degree (denominator): ", e.datasize - e.degree - 1) + if isnothing(e.err) + println(io, "Error: Not yet calculated") + else + println(io, "Error: ", @sprintf("%.4e", e.err)) + end + # if isnothing(e.at_time) + # println(io, "Time: Not specified") + # else + # println(io, "Time: ", @sprintf("%.4e", e.at_time)) + # end + # end + # println(io, "Return Code: ", e.return_code) end diff --git a/src/utils.jl b/src/utils.jl index 084d95f8..8e548ed3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,148 +1,147 @@ """ - nemo2hc(expr_tree::Union{Expr, Symbol}) + nemo2hc(expr_tree::Union{Expr, Symbol}) Converts a symbolic expression from Nemo to HomotopyContinuation format. """ function nemo2hc(expr_tree::Union{Expr, Symbol}) - #traverse expr_tree - if typeof(expr_tree) == Symbol - return HomotopyContinuation.Expression(HomotopyContinuation.variables(expr_tree)[1]) - end - if typeof(expr_tree) == Expr - if expr_tree.head == :call - if expr_tree.args[1] in [:+, :-, :*, :/, :^, ://] - if length(expr_tree.args) == 2 - return eval(expr_tree.args[1])(nemo2hc(expr_tree.args[2])) - else - # println(expr_tree.args[2:end]) - return reduce(eval(expr_tree.args[1]), - map(nemo2hc, expr_tree.args[2:end])) - end - end - end - end + #traverse expr_tree + if typeof(expr_tree) == Symbol + return HomotopyContinuation.Expression(HomotopyContinuation.variables(expr_tree)[1]) + end + if typeof(expr_tree) == Expr + if expr_tree.head == :call + if expr_tree.args[1] in [:+, :-, :*, :/, :^, ://] + if length(expr_tree.args) == 2 + return eval(expr_tree.args[1])(nemo2hc(expr_tree.args[2])) + else + # println(expr_tree.args[2:end]) + return reduce(eval(expr_tree.args[1]), + map(nemo2hc, expr_tree.args[2:end])) + end + end + end + end end function nemo2hc(expr_tree::fmpq_mpoly) - # println(expr_tree) - return nemo2hc(Meta.parse(string(expr_tree))) + # println(expr_tree) + return nemo2hc(Meta.parse(string(expr_tree))) end function nemo2hc(expr_tree::Number) - return expr_tree + return expr_tree end function nemo2hc(expr_tree::Oscar.Generic.Frac) - numer, denom = Oscar.numerator(expr_tree), Oscar.denominator(expr_tree) - return nemo2hc(numer) / nemo2hc(denom) + numer, denom = Oscar.numerator(expr_tree), Oscar.denominator(expr_tree) + return nemo2hc(numer) / nemo2hc(denom) end """ - squarify_system(poly_system::Vector{Expression}) + squarify_system(poly_system::Vector{Expression}) Given a non-square polynomial system in `n` variables, takes first `n - 1` equations and adds a random linear combination of the remaining equations to the system. """ function squarify_system(poly_system::Vector{Expression}) - indets = HomotopyContinuation.variables(poly_system) - M = randn(1, length(poly_system) - length(indets) + 1) - return vcat(poly_system[1:(length(indets) - 1)], M * poly_system[length(indets):end]) + indets = HomotopyContinuation.variables(poly_system) + M = randn(1, length(poly_system) - length(indets) + 1) + return vcat(poly_system[1:(length(indets)-1)], M * poly_system[length(indets):end]) end """ - check_inputs(measured_quantities::Vector{ModelingToolkit.Equation} = Vector{ModelingToolkit.Equation}([]), - data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}(), - time_interval::Vector{T} = Vector{T}(), - interpolation_degree::Int = 1) where {T <: Float} + check_inputs(measured_quantities::Vector{ModelingToolkit.Equation} = Vector{ModelingToolkit.Equation}([]), + data_sample::Dict{Any, Vector{T}} = Dict{Any, Vector{T}}(), + time_interval::Vector{T} = Vector{T}(), + interpolation_degree::Int = 1) where {T <: Float} Checks that the inputs to `estimate` are valid. """ function check_inputs(measured_quantities::Vector{ModelingToolkit.Equation} = Vector{ - ModelingToolkit.Equation - }([]), - data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}(), - interpolation_degree::Int = 1) where {T <: Float} - if length(measured_quantities) == 0 - error("No measured states provided") - end - if length(data_sample) == 0 - error("No data sample provided") - end - if length(first(values(data_sample))) <= 2 - error("Data sample must have at least 3 points") - end - if interpolation_degree < 1 - error("Interpolation degree must be ≥ 1") - end + ModelingToolkit.Equation, + }([]), + data_sample::AbstractDict{Any, Vector{T}} = Dict{Any, Vector{T}}()) where {T <: Float} + if length(measured_quantities) == 0 + error("No measured states provided") + end + if length(data_sample) == 0 + error("No data sample provided") + end + if length(first(values(data_sample))) <= 2 + error("Data sample must have at least 3 points") + end + #if interpolation_degree < 1 + # error("Interpolation degree must be ≥ 1") + #end end -to_real(x::Number; tol = 1e-10) = abs(imag(x)) < tol ? real(x) : x -function to_exact(x::Number; tol = 1e-10) - r = abs(real(x)) < tol ? 0 : real(x) - i = abs(imag(x)) < tol ? 0 : imag(x) - if i == 0 - return r - else - return r + i * im - end +to_real(x::Number; tol = 1e-12) = abs(imag(x)) < tol ? real(x) : x +function to_exact(x::Number; tol = 1e-12) + r = abs(real(x)) < tol ? 0 : real(x) + i = abs(imag(x)) < tol ? 0 : imag(x) + if i == 0 + return r + else + return r + i * im + end end function sample_data(model::ModelingToolkit.ODESystem, - measured_data::Vector{ModelingToolkit.Equation}, - time_interval::Vector{T}, - p_true::Vector{T}, - u0::Vector{T}, - num_points::Int; - uneven_sampling = false, - uneven_sampling_times = Vector{T}(), - solver = Tsit5(), inject_noise = false, mean_noise = 0, - stddev_noise = 1) where {T <: Float} - if uneven_sampling - if length(uneven_sampling_times) == 0 - error("No uneven sampling times provided") - end - if length(uneven_sampling_times) != num_points - error("Uneven sampling times must be of length num_points") - end - sampling_times = uneven_sampling_times - else - sampling_times = range(time_interval[1], time_interval[2], length = num_points) - end - problem = ODEProblem(model, u0, time_interval, p_true) - solution_true = ModelingToolkit.solve(problem, solver, p = p_true, - saveat = sampling_times; - abstol = 1e-10, reltol = 1e-10) - data_sample = OrderedDict{Any, Vector{T}}(Num(v.rhs) => solution_true[Num(v.rhs)] - for v in measured_data) - if inject_noise - for (key, sample) in data_sample - data_sample[key] = sample + randn(num_points) .* stddev_noise .+ mean_noise - end - end - data_sample["t"] = sampling_times - return data_sample + measured_data::Vector{ModelingToolkit.Equation}, + time_interval::Vector{T}, + p_true::Vector{T}, + u0::Vector{T}, + num_points::Int; + uneven_sampling = false, + uneven_sampling_times = Vector{T}(), + solver = Vern9(), inject_noise = false, mean_noise = 0, + stddev_noise = 1, abstol = 1e-14, reltol = 1e-14) where {T <: Number} + if uneven_sampling + if length(uneven_sampling_times) == 0 + error("No uneven sampling times provided") + end + if length(uneven_sampling_times) != num_points + error("Uneven sampling times must be of length num_points") + end + sampling_times = uneven_sampling_times + else + sampling_times = range(time_interval[1], time_interval[2], length = num_points) + end + problem = ODEProblem(model, u0, time_interval, p_true) + solution_true = ModelingToolkit.solve(problem, solver, p = p_true, + saveat = sampling_times; + abstol, reltol) + data_sample = OrderedDict{Any, Vector{T}}(Num(v.rhs) => solution_true[Num(v.rhs)] + for v in measured_data) + if inject_noise + for (key, sample) in data_sample + data_sample[key] = sample + randn(num_points) .* stddev_noise .+ mean_noise + end + end + data_sample["t"] = sampling_times + return data_sample end function write_sample(data_sample; filename = "sample_data.txt") - open(filename, "w") do io - idx_iter = eachindex(first(values(data_sample))) - # print keys - for (key, sample) in data_sample - print(io, key, " ") - end - println(io) - for i in idx_iter - for (key, sample) in data_sample - print(io, sample[i], " ") - end - println(io) - end - end + open(filename, "w") do io + idx_iter = eachindex(first(values(data_sample))) + # print keys + for (key, sample) in data_sample + print(io, key, " ") + end + println(io) + for i in idx_iter + for (key, sample) in data_sample + print(io, sample[i], " ") + end + println(io) + end + end end function max_abs_rel(all_params, estimate::EstimationResult) - estimates = vcat(collect(values(estimate.states)), collect(values(estimate.parameters))) - err = maximum(abs.((result_ode.u .- all_params) ./ (all_params))) - println("Max abs rel. err: ", err) - return err + estimates = vcat(collect(values(estimate.states)), collect(values(estimate.parameters))) + err = maximum(abs.((result_ode.u .- all_params) ./ (all_params))) + println("Max abs rel. err: ", err) + return err end diff --git a/test/homotopy_estimation_tests.jl b/test/homotopy_estimation_tests.jl index ae90a36f..7d0ce89a 100644 --- a/test/homotopy_estimation_tests.jl +++ b/test/homotopy_estimation_tests.jl @@ -14,7 +14,7 @@ x1 + x1^2 => [2.0, 1.56301, 1.22995, 0.97441]) res = ParameterEstimation.estimate(model, outs, data) - @test length(res) == 1 + #@test length(res) == 1 @test Symbol(res[1].return_code) == :Success @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) @@ -23,7 +23,7 @@ data = Dict{Any, Vector{Float64}}(x1^2 + x1 => [2.0, 1.98506, 1.17611, 0.97441], "t" => [0.0, 0.01, 0.73, 1.0]) res = ParameterEstimation.estimate(model, outs, data) - @test length(res) == 1 + #@test length(res) == 1 @test Symbol(res[1].return_code) == :Success @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) diff --git a/test/msolve_estimation_tests.jl b/test/msolve_estimation_tests.jl index 6e9ccb27..98bc85df 100644 --- a/test/msolve_estimation_tests.jl +++ b/test/msolve_estimation_tests.jl @@ -14,7 +14,7 @@ x1 + x1^2 => [2.0, 1.56301, 1.22995, 0.97441]) res = ParameterEstimation.estimate(model, outs, data; method = :msolve) - @test length(res) == 1 + # @test length(res) == 1 @test Symbol(res[1].return_code) == :Success @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) @@ -23,7 +23,7 @@ data = Dict{Any, Vector{Float64}}(x1^2 + x1 => [2.0, 1.98506, 1.17611, 0.97441], "t" => [0.0, 0.01, 0.73, 1.0]) res = ParameterEstimation.estimate(model, outs, data; method = :msolve) - @test length(res) == 1 + #@test length(res) == 1 @test Symbol(res[1].return_code) == :Success @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) diff --git a/test/varied_estimation.jl b/test/varied_estimation.jl new file mode 100644 index 00000000..a760dd1d --- /dev/null +++ b/test/varied_estimation.jl @@ -0,0 +1,590 @@ +@testset "Run larger, slower parameter recovery tests on known ODEs" begin +using ParameterEstimation +using ModelingToolkit, DifferentialEquations#, Plots + +struct ParameterEstimationProblem + Name::Any + model::Any + measured_quantities::Any + data_sample::Any + solver::Any + p_true::Any + ic::Any +end + +function biohydrogenation(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters k5 k6 k7 k8 k9 k10 + @variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t) + D = Differential(t) + states = [x4, x5, x6, x7] + parameters = [k5, k6, k7, k8, k9, k10] + + @named model = ODESystem([ + D(x4) ~ -k5 * x4 / (k6 + x4), + D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6), + D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10, + D(x7) ~ k9 * x6 * (k10 - x6) / k10, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x4, + y2 ~ x5, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("BioHydrogenation", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function crauste(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters mu_N mu_EE mu_LE mu_LL mu_M mu_P mu_PE mu_PL delta_NE delta_EL delta_LM rho_E rho_P + @variables t N(t) E(t) S(t) M(t) P(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [N, E, S, M, P] + parameters = [ + mu_N, + mu_EE, + mu_LE, + mu_LL, + mu_M, + mu_P, + mu_PE, + mu_PL, + delta_NE, + delta_EL, + delta_LM, + rho_E, + rho_P, + ] + @named model = ODESystem([ + D(N) ~ -N * mu_N - N * P * delta_NE, + D(E) ~ N * P * delta_NE - E^2 * mu_EE - + E * delta_EL + E * P * rho_E, + D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL - + E * S * mu_LE, + D(M) ~ S * delta_LM - mu_M * M, + D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE - + S * P * mu_PL, + ], t, states, parameters) + measured_quantities = [y1 ~ N, y2 ~ E, y3 ~ S + M, y4 ~ P] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [ + 0.071, + 0.143, + 0.214, + 0.286, + 0.357, + 0.429, + 0.5, + 0.571, + 0.643, + 0.714, + 0.786, + 0.857, + 0.929, + ] # True Parameters + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Crauste", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function daisy_ex3(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters p1 p3 p4 p6 p7 + @variables t x1(t) x2(t) x3(t) u0(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, u0] + parameters = [p1, p3, p4, p6, p7] + @named model = ODESystem([ + D(x1) ~ -1 * p1 * x1 + x2 + u0, + D(x2) ~ p3 * x1 - p4 * x2 + x3, + D(x3) ~ p6 * x1 - p7 * x3, + D(u0) ~ 1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ u0, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_ex3", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function daisy_mamil3(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a12 a13 a21 a31 a01 + @variables t x1(t) x2(t) x3(t) y1(t) y2(t) + D = Differential(t) + + ic = [0.25, 0.5, 0.75] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] # True Parameters + + states = [x1, x2, x3] + parameters = [a12, a13, a21, a31, a01] + @named model = ODESystem([D(x1) ~ -(a21 + a31 + a01) * x1 + a12 * x2 + a13 * x3, + D(x2) ~ a21 * x1 - a12 * x2, + D(x3) ~ a31 * x1 - a13 * x3], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil3", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function daisy_mamil4(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters k01, k12, k13, k14, k21, k31, k41 + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) y3(t) + D = Differential(t) + + ic = [0.2, 0.4, 0.6, 0.8] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] # True Parameters + + states = [x1, x2, x3, x4] + parameters = [k01, k12, k13, k14, k21, k31, k41] + @named model = ODESystem([ + D(x1) ~ -k01 * x1 + k12 * x2 + k13 * x3 + k14 * x4 - k21 * x1 - k31 * x1 - + k41 * x1, + D(x2) ~ -k12 * x2 + k21 * x1, + D(x3) ~ -k13 * x3 + k31 * x1, + D(x4) ~ -k14 * x4 + k41 * x1], + t, states, parameters) + measured_quantities = [y1 ~ x1, y2 ~ x2, y3 ~ x3 + x4] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("DAISY_mamil4", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function fitzhugh_nagumo(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters g a b + @variables t V(t) R(t) y1(t) y2(t) + D = Differential(t) + states = [V, R] + parameters = [g, a, b] + + ic = [0.333, 0.67] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ V] + + @named model = ODESystem([ + D(V) ~ g * (V - V^3 / 3 + R), + D(R) ~ 1 / g * (V - a + b * R), + ], t, states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("fitzhugh-nagumo", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function hiv_local(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters b c d k1 k2 mu1 mu2 q1 q2 s + @variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2, x3, x4] + parameters = [b, c, d, k1, k2, mu1, mu2, q1, q2, s] + + @named model = ODESystem([ + D(x1) ~ -b * x1 * x4 - d * x1 + s, + D(x2) ~ b * q1 * x1 * x4 - k1 * x2 - mu1 * x2, + D(x3) ~ b * q2 * x1 * x4 + k1 * x2 - mu2 * x3, + D(x4) ~ -c * x4 + k2 * x3, + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x4, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.091, 0.182, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + time_interval = [-0.5, 0.5] + datasize = 20 + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("hiv_local", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function hiv(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters lm d beta a k u c q b h + @variables t x(t) y(t) v(t) w(t) z(t) y1(t) y2(t) y3(t) y4(t) + D = Differential(t) + states = [x, y, v, w, z] + parameters = [lm, d, beta, a, k, u, c, q, b, h] + + @named model = ODESystem([ + D(x) ~ lm - d * x - beta * x * v, + D(y) ~ beta * x * v - a * y, + D(v) ~ k * y - u * v, + D(w) ~ c * x * y * w - c * q * y * w - b * w, + D(z) ~ c * q * y * w - h * z, + ], t, states, parameters) + measured_quantities = [y1 ~ w, y2 ~ z, y3 ~ x, y4 ~ y + v] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.091, 0.181, 0.273, 0.364, 0.455, 0.545, 0.636, 0.727, 0.818, 0.909] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + + res = ParameterEstimation.estimate(model, measured_quantities, data_sample; + solver = solver) + return ParameterEstimationProblem("hiv", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function lotka_volterra(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters k1 k2 k3 + @variables t r(t) w(t) y1(t) + D = Differential(t) + ic = [0.333, 0.667] + sampling_times = range(time_interval[1], time_interval[2], length = datasize) + p_true = [0.25, 0.5, 0.75] # True Parameters + measured_quantities = [y1 ~ r] + states = [r, w] + parameters = [k1, k2, k3] + + @named model = ODESystem([D(r) ~ k1 * r - k2 * r * w, D(w) ~ k2 * r * w - k3 * w], t, + states, parameters) + + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, datasize; solver = solver) + + return ParameterEstimationProblem("Lotka_Volterra", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function seir(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a b nu + @variables t S(t) E(t) In(t) N(t) y1(t) y2(t) + D = Differential(t) + states = [S, E, In, N] + parameters = [a, b, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N, + D(E) ~ b * S * In / N - nu * E, + D(In) ~ nu * E - a * In, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ In, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.25, 0.5, 0.75] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("SEIR", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function simple(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ -a * x2, + D(x2) ~ 1 / b * (x1), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("simple", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function sirsforced(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) + @parameters b0 b1 g M mu nu + @variables t i(t) r(t) s(t) x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [i, r, s, x1, x2] + parameters = [b0, b1, g, M, mu, nu] + + @named model = ODESystem([ + D(s) ~ mu - mu * s - b0 * (1 + b1 * x1) * i * s + g * r, + D(i) ~ b0 * (1 + b1 * x1) * i * s - (nu + mu) * i, + D(r) ~ nu * i - (mu + g) * r, + D(x1) ~ -M * x2, + D(x2) ~ M * x1, + ], t, states, parameters) + measured_quantities = [ + y1 ~ i, + y2 ~ r, + ] + + ic = [0.167, 0.333, 0.5, 0.667, 0.833] + p_true = [0.143, 0.286, 0.429, 0.571, 0.714, 0.857] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("sirsforced", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function slowfast(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) # TODO(orebas):in the old code it was CVODE_BDF. should we go back to that? + #solver = CVODE_BDF() + @parameters k1 k2 eB + @variables t xA(t) xB(t) xC(t) eA(t) eC(t) y1(t) y2(t) y3(t) y4(t) #eA(t) eC(t) + D = Differential(t) + states = [xA, xB, xC, eA, eC] + parameters = [k1, k2, eB] + @named model = ODESystem([ + D(xA) ~ -k1 * xA, + D(xB) ~ k1 * xA - k2 * xB, + D(xC) ~ k2 * xB, + D(eA) ~ 0, + D(eC) ~ 0, + ], t, states, parameters) + + measured_quantities = [y1 ~ xC, y2 ~ eA * xA + eB * xB + eC * xC, y3 ~ eA, y4 ~ eC] + ic = [0.166, 0.333, 0.5, 0.666, 0.833] + p_true = [0.25, 0.5, 0.75] # True Parameters + + data_sample = ParameterEstimation.sample_data(model, + measured_quantities, + time_interval, + p_true, + ic, + datasize; + solver = solver) + + return ParameterEstimationProblem("slowfast", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function treatment(datasize = 21, time_interval = [-0.5, 0.5], solver = Rodas5P()) #note the solver. Vern9 apparently can't handle mass matrices + @parameters a b d g nu + @variables t In(t) N(t) S(t) Tr(t) y1(t) y2(t) + D = Differential(t) + states = [In, N, S, Tr] + parameters = [a, b, d, g, nu] + + @named model = ODESystem([ + D(S) ~ -b * S * In / N - d * b * S * Tr / N, + D(In) ~ b * S * In / N + d * b * S * Tr / N - (a + g) * In, + D(Tr) ~ g * In - nu * Tr, + D(N) ~ 0, + ], t, states, parameters) + measured_quantities = [ + y1 ~ Tr, + y2 ~ N, + ] + + ic = [0.2, 0.4, 0.6, 0.8] + p_true = [0.167, 0.333, 0.5, 0.667, 0.833] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize, solver = solver) + return ParameterEstimationProblem("treatment", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function vanderpol(datasize = 21, time_interval = [-0.5, 0.5], solver = Vern9()) + @parameters a b + @variables t x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + states = [x1, x2] + parameters = [a, b] + + @named model = ODESystem([ + D(x1) ~ a * x2, + D(x2) ~ -(x1) - b * (x1^2 - 1) * (x2), + ], t, states, parameters) + measured_quantities = [ + y1 ~ x1, + y2 ~ x2, + ] + + ic = [0.333, 0.667] + p_true = [0.333, 0.667] + data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, + p_true, ic, + datasize; solver = solver) + return ParameterEstimationProblem("vanderpol", + model, + measured_quantities, + data_sample, + solver, + p_true, + ic) +end + +function analyze_parameter_estimation_problem(PEP::ParameterEstimationProblem; test_mode=false) + + #interpolators = Dict( + #"AAA" => ParameterEstimation.aaad, + #"FHD3" => ParameterEstimation.fhdn(3), + #"FHD6" => ParameterEstimation.fhdn(6), + #"FHD8" => ParameterEstimation.fhdn(8), + #"Fourier" => ParameterEstimation.FourierInterp, + #) + datasize = 21 #TODO(Orebas) magic number + + #stepsize = max(1, datasize ÷ 8) + #for i in range(1, (datasize - 2), step = stepsize) + # interpolators["RatOld($i)"] = ParameterEstimation.SimpleRationalInterpOld(i) + #end + + @time res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, + PEP.data_sample, + solver = PEP.solver) # , interpolators) + all_params = vcat(PEP.ic, PEP.p_true) + #println("TYPERES: ", typeof(res)) + #println(res) + + #println(res) + besterror = 1e30 + for each in res + #println("TYPE: ", typeof(each)) + #println(each) + estimates = vcat(collect(values(each.states)), collect(values(each.parameters))) + besterror = min(besterror, maximum(abs.((estimates .- all_params) ./ (all_params)))) + end + if (test_mode) + @test besterror < 1e-2 + end + println("For model ", PEP.Name, ": The best max abs rel. err: ", besterror) +end + + +#function analyze_parameter_estimation_problem_old(PEP::ParameterEstimationProblem) +# res = ParameterEstimation.estimate(PEP.model, PEP.measured_quantities, PEP.data_sample; +# solver = PEP.solver, interpolator = ("AAA" => aaad)) +# all_params = vcat(PEP.ic, PEP.p_true) +# for each in res +# estimates = vcat(collect(values(each.states)), collect(values(each.parameters))) +# println("Max abs rel. err: ", +# maximum(abs.((estimates .- all_params) ./ (all_params)))) +# end +#end + +function varied_estimation_main() + datasize = 21 + solver = Vern9() + #solver = Rodas4P() + time_interval = [-0.5, 0.5] + for PEP in [ + simple(datasize, time_interval, solver), #works + lotka_volterra(datasize, time_interval, solver), #works + vanderpol(datasize, time_interval, solver), #works + #biohydrogenation(datasize, time_interval, solver), #works, but one param unidentifiable + #daisy_ex3(datasize, time_interval, solver), + daisy_mamil3(datasize, time_interval, solver), + daisy_mamil4(datasize, time_interval, solver), + #fitzhugh_nagumo(datasize, time_interval, solver), + #hiv_local(datasize, time_interval, solver), + hiv(datasize, time_interval, solver), + seir(datasize, time_interval, solver), + #sirsforced(datasize, time_interval, Rodas5P()), + slowfast(datasize, time_interval, solver), + #treatment(datasize, time_interval, Rodas5P()), + crauste(datasize, time_interval, solver), + ] + analyze_parameter_estimation_problem(PEP, test_mode=true) + end +end + +varied_estimation_main() +end \ No newline at end of file