Merge pull request #97 from sumiya11/spring-clean

Spring cleaning: update versions, remove Oscar and msolve

Project.toml

@@ -14,8 +14,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 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"
@@ -33,19 +33,26 @@ ArbNumerics = "1"
 BaryRational = "1"
 DifferentialEquations = "7"
 ForwardDiff = "0.10"
-Groebner = "0.2, 0.3, 0.4 - 0.5"
+Groebner = "0.6, 0.7"
 HomotopyContinuation = "2"
 LinearSolve = "2"
-ModelingToolkit = "8"
+ModelingToolkit = "9"
 OrderedCollections = "1"
-Oscar = "0.10 - 0.20"
 PrecompileTools = "1"
 ProgressMeter = "1"
-SIAN = "1 - 1.5.10"
-StructuralIdentifiability = "0.4.9 - 0.5.30"
+SIAN = "1.5"
+StructuralIdentifiability = "0.5.7"
 Suppressor = "0.2"
 Symbolics = "4, 5"
 TaylorDiff = "0.2"
 TaylorSeries = "0.12 - 0.30"
+Test = "1.10"
 TestSetExtensions = "2"
-julia = "1.9"
+julia = "1.10"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
+
+[targets]
+test = ["Test", "TestSetExtensions"]

src/ParameterEstimation.jl

@@ -15,9 +15,10 @@ using Suppressor
 
 using ProgressMeter, Logging, Printf
 using ModelingToolkit, LinearSolve, LinearAlgebra
-using SIAN, HomotopyContinuation, Groebner, Oscar
+using SIAN, HomotopyContinuation, Groebner, Nemo
 using .ReturnCode
-import StructuralIdentifiability: eval_at_nemo, ODE
+import StructuralIdentifiability
+import StructuralIdentifiability: ODE
 using BaryRational
 using ForwardDiff
 using ArbNumerics
@@ -53,10 +54,10 @@ export check_identifiability, estimate, filter_solutions
 		p_true = [2.0, 3.0, 4.0, 5.0]
 		time_interval = [-4.0, 4.0]
 		datasize = 9
+		model = complete(model)
 		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
-

src/estimation/fixed_degree.jl

@@ -2,36 +2,40 @@
 	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)]
+		initial_conditions = [E[s] for s in ModelingToolkit.unknowns(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)]
+	initial_conditions = [E[s] for s in ModelingToolkit.unknowns(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),
+	@named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model),
 		ModelingToolkit.parameters(model))
-	prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values)
+	prob = ODEProblem(
+		ModelingToolkit.complete(new_model), 
+		initial_conditions, 
+		tspan, 
+		Dict(ModelingToolkit.parameters(model) .=> 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)
+	ode_solution = ModelingToolkit.solve(prob, solver, saveat = saveat, abstol = abstol, reltol = reltol)
 
 	state_param_map = (Dict(x => replace(string(x), "(t)" => "")
-							for x in ModelingToolkit.states(model)))
+							for x in ModelingToolkit.unknowns(model)))
 
 
 	newstates = copy(E.states)
 
-	for s in ModelingToolkit.states(model)
+	for s in ModelingToolkit.unknowns(model)
 		temp = ode_solution[Symbol(state_param_map[s])][end]
 		newstates[s] = temp
 	end
@@ -94,7 +98,7 @@ function estimate_single_interpolator(model::ModelingToolkit.ODESystem,
 	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)
+	states = ModelingToolkit.unknowns(model)
 	num_parameters = length(parameters) + length(states)
 	@debug "Interpolating sample data via interpolation method $(interpolator.first)"
 	if !haskey(data_sample, "t")
@@ -107,7 +111,7 @@ function estimate_single_interpolator(model::ModelingToolkit.ODESystem,
 		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)
@@ -117,8 +121,7 @@ function estimate_single_interpolator(model::ModelingToolkit.ODESystem,
 	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]

src/estimation/serial.jl

@@ -38,5 +38,5 @@ function estimate_serial(model::ModelingToolkit.ODESystem,
                 ])
         end
     end
-    return post_process(estimates, filtermode, parameter_constraints, ic_constraints)
+    return post_process(estimates, filtermode, parameter_constraints, ic_constraints; threaded=false)
 end

src/estimation/solve.jl

@@ -3,7 +3,7 @@ function solve_via_homotopy(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)),
+								 for x in ModelingToolkit.unknowns(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)
@@ -34,38 +34,39 @@ function solve_via_homotopy(identifiability_result, model; real_tol = 1e-12)
 end
 
 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)))
+    throw("Not implemented, sorry.")
+    # 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, _ = Nemo.polynomial_ring(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)
-	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_
+    # 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

src/estimation/threaded.jl

@@ -8,7 +8,7 @@ function estimate_threaded(model::ModelingToolkit.ODESystem,
 	check_inputs(measured_quantities, data_sample)
 	datasize = length(first(values(data_sample)))
 
-    if isnothing(interpolators) === nothing
+    if isnothing(interpolators)
         interpolators = default_interpolator(datasize)
     end
 
@@ -20,12 +20,13 @@ estimates = Vector{Vector{ParameterEstimation.EstimationResult}}()
 @info "Estimating via the interpolators: $(keys(interpolators))"
 	n_threads = Threads.nthreads()
 	N = length(interpolators)
+	interpolators = collect(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] = []
+		interp = interpolators[t]
+		tid = Threads.threadid()
+		if isnothing(estimates[tid])
+			estimates[tid] = []
 		end
 		unfiltered = estimate_single_interpolator(model, measured_quantities, inputs, data_sample;
 			identifiability_result = id,
@@ -35,16 +36,16 @@ estimates = Vector{Vector{ParameterEstimation.EstimationResult}}()
 		if length(unfiltered) > 0
 			filtered = filter_solutions(unfiltered, id, model, inputs,
 				data_sample; solver = solver, filtermode)
-			push!(estimates[id], filtered)
+			push!(estimates[tid], filtered)
 
 		else
-			push!(estimates[id],
+			push!(estimates[tid],
 				[
 					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)
+	return post_process(estimates,filtermode,parameter_constraints,ic_constraints,threaded=true)
 end

src/estimation/utils.jl

@@ -68,8 +68,8 @@ function new_clustering(estimates, tol = 1e-3)
 	return new_estimates
 end
 
-function post_process(estimates, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing)
-	if Threads.nthreads() > 1
+function post_process(estimates, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing; threaded=false)
+	if threaded
 		estimates = filter(x -> !isnothing(x), estimates)
 		estimates = vcat(estimates...)
 	end

src/filtering.jl

@@ -20,16 +20,20 @@ Compute the error between the solution and the data sample. The error is recorde
 """
 function solve_ode(model, estimate::EstimationResult, inputs::Vector{Equation}, data_sample;
 	solver = Tsit5(), return_ode = false, abstol = 1e-12, reltol = 1e-12)
-	initial_conditions = [estimate[s] for s in ModelingToolkit.states(model)]
+	initial_conditions = [estimate[s] for s in ModelingToolkit.unknowns(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),
+	@named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model),
 		ModelingToolkit.parameters(model))
-	prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values)
+	prob = ODEProblem(
+		ModelingToolkit.complete(new_model), 
+		initial_conditions, tspan, 
+		Dict(ModelingToolkit.parameters(new_model) .=> 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)