Add the option to dump polynomial system to files

.JuliaFormatter.toml

@@ -5,4 +5,4 @@ always_for_in = true
 whitespace_typedefs = true
 normalize_line_endings = "unix"
 short_to_long_function_def = false
-format_markdown = true
+format_markdown = true

src/estimation/estimate.jl

@@ -33,7 +33,7 @@ function estimate(model::ModelingToolkit.ODESystem,
 	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}
+	threaded = Threads.nthreads() > 1, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing, dump_systems_to=nothing) where {T <: Float}
 
 	#println("DEBUG")
 	if !(method in [:homotopy, :msolve])
@@ -51,7 +51,7 @@ function estimate(model::ModelingToolkit.ODESystem,
 			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)
+			real_tol = real_tol, filtermode, parameter_constraints = parameter_constraints, ic_constraints = ic_constraints, dump_systems_to=dump_systems_to)
 	end
 	println("Final Results:")
 	for each in result

src/estimation/fixed_degree.jl

@@ -8,48 +8,43 @@
 
 		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.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.unknowns(model),
-		ModelingToolkit.parameters(model))
-	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, saveat = saveat, abstol = abstol, reltol = reltol)
-
-	state_param_map = (Dict(x => replace(string(x), "(t)" => "")
-							for x in ModelingToolkit.unknowns(model)))
-
-
-	newstates = copy(E.states)
-
-	for s in ModelingToolkit.unknowns(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
-
+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.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.unknowns(model),
+        ModelingToolkit.parameters(model))
+    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, saveat = saveat, abstol = abstol, reltol = reltol)
+
+    state_param_map = (Dict(x => replace(string(x), "(t)" => "")
+    for x in ModelingToolkit.unknowns(model)))
+
+    newstates = copy(E.states)
+
+    for s in ModelingToolkit.unknowns(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},
@@ -83,53 +78,120 @@ measured quantities `measured_quantities`.
 - `EstimationResult`: the estimated parameters and initial conditions of the model.
 """
 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.unknowns(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
-
-	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
+        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,
+        dump_systems_to = nothing) 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.unknowns(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 dump_systems_to != nothing
+        if !isdir(dump_systems_to)
+            error("Directory $dump_systems_to does not exist")
+        end
+        polynomial_system = identifiability_result["polynomial_system_to_solve"]
+
+	tag = interpolator[1]
+        # HC.jl
+	io = open("$(dump_systems_to)/$(tag)_hc.txt", "w")
+        println(io, polynomial_system)
+        println(io)
+        close(io)
+
+        # PHCpack
+        io = open("$(dump_systems_to)/$(tag)_phc.txt", "w")
+        input = String[]
+        push!(input,
+            "$(HomotopyContinuation.length(polynomial_system)) $(HomotopyContinuation.nvariables(polynomial_system))")
+        n = length(polynomial_system)
+        for (i, fi) in enumerate(polynomial_system)
+            push!(input, "$(fi);")
+        end
+        # replace exponentiation
+        input = map(line -> replace(line, "^" => "**"), input)
+        # replace complex numbers
+        input = map(line -> replace(line, r"(^|\W)(im)(\W|$)" => s"\1I\3"), input)
+        println(io, join(input, "\n"))
+        close(io)
+
+        # Bertini
+        io = open("$(dump_systems_to)/$(tag)_bertini.txt", "w")
+        n = length(polynomial_system)
+        funcs = ["f$i" for i in 1:n]
+        var_group = map(string, HomotopyContinuation.variables(polynomial_system))
+        input = String[]
+        for (i, fi) in enumerate(polynomial_system)
+            push!(input, "$(funcs[i]) = $(fi);")
+        end
+        println(io, """
+CONFIG
+MPTYPE: 2;
+END;""")
+        println(io, "variable_group $(join(var_group, ", "));")
+        println(io, "function $(join(funcs, ", "));")
+        println(io, join(input, "\n"))
+        println(io, "END")
+        close(io)
+
+        # AA.jl
+        io = open("$(dump_systems_to)/$(tag)_aa.txt", "w")
+        polynomial_system_aa = hc2nemo(polynomial_system)
+        v = gens(parent(polynomial_system_aa[1]))
+        println(io)
+        println(io,
+            "R, ($(join(string.(v), ","))) = polynomial_ring(QQ, [$(join(map(_v -> "\""*_v*"\"", string.(v)), ","))], internal_ordering=:degrevlex)")
+        println(io)
+        println(io, "sys = [")
+        for eq in polynomial_system_aa
+            println(io, eq, ",")
+        end
+        println(io, "]")
+        println(io)
+        close(io)
+    end
+
+    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)
+    elseif method == :rur
+        all_solutions = solve_via_rur(identifiability_result, model; real_tol = real_tol)
+    else
+        throw(ArgumentError("Method $method not supported"))
+    end
+
+    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

src/estimation/serial.jl

@@ -6,7 +6,8 @@ function estimate_serial(model::ModelingToolkit.ODESystem,
         report_time = minimum(data_sample["t"]),
         method = :homotopy,
         real_tol::Float64 = 1e-14, filtermode = :new, parameter_constraints = nothing,
-        ic_constraints = nothing) where {T <: Float}
+        ic_constraints = nothing,
+	dump_systems_to=nothing) where {T <: Float}
     check_inputs(measured_quantities, data_sample)
     datasize = length(first(values(data_sample)))
 
@@ -24,7 +25,7 @@ function estimate_serial(model::ModelingToolkit.ODESystem,
             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)
+            method = method, real_tol = real_tol, dump_systems_to=dump_systems_to)
         if length(unfiltered) > 0
             filtered = filter_solutions(unfiltered, id, model, inputs, data_sample;
                 solver = solver, filtermode)

src/estimation/solve.jl

@@ -1,36 +1,37 @@
 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"
+    @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.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)
-	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.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)
+    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-12)

src/utils.jl

@@ -9,6 +9,9 @@ function nemo2hc(expr_tree::Union{Expr, Symbol})
 		return HomotopyContinuation.Expression(HomotopyContinuation.variables(expr_tree)[1])
 	end
 	if typeof(expr_tree) == Expr
+		if expr_tree.head == :macrocall
+			return eval(expr_tree)
+		end
 		if expr_tree.head == :call
 			if expr_tree.args[1] in [:+, :-, :*, :/, :^, ://]
 				if length(expr_tree.args) == 2
@@ -37,6 +40,18 @@ function nemo2hc(expr_tree::Nemo.Generic.FracFieldElem)
 	return nemo2hc(numer) / nemo2hc(denom)
 end
 
+function hc2nemo(sys)
+        vars = variables(HomotopyContinuation.System(sys))
+        ring, aa_vars = polynomial_ring(QQ, map(string, vars), internal_ordering=:degrevlex)
+        dicts = map(f -> HomotopyContinuation.to_dict(HomotopyContinuation.expand(f), vars), sys)
+        aa_sys = Vector{elem_type(ring)}()
+        for dict in dicts
+                poly = sum(e_c -> rationalize(BigInt, to_number(e_c[2]); tol=0) * prod(aa_vars .^ e_c[1]), collect(dict))
+                push!(aa_sys, poly)
+        end
+        aa_sys
+end
+
 """
 	squarify_system(poly_system::Vector{Expression})