Skip to content

Commit

Permalink
Merge pull request #143 from JuliaDiffEq/solve
Browse files Browse the repository at this point in the history
Remove solve warning
  • Loading branch information
ChrisRackauckas committed Jul 30, 2019
2 parents 3a726dd + 195c006 commit 923f8d9
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 13 deletions.
4 changes: 3 additions & 1 deletion src/DiffEqBiological.jl
Expand Up @@ -5,7 +5,9 @@ module DiffEqBiological
using Compat, DataStructures, MacroTools, Parameters, Reexport, SparseArrays, SymEngine
using DiffEqBase, DiffEqJump
@reexport using DiffEqBase, DiffEqJump
using HomotopyContinuation, DynamicPolynomials, LinearAlgebra, RecipesBase
using DynamicPolynomials, LinearAlgebra, RecipesBase

import HomotopyContinuation

import Base: (==)

Expand Down
24 changes: 12 additions & 12 deletions src/equilibrate_utils.jl
Expand Up @@ -13,9 +13,9 @@ mutable struct EquilibrateContent
end

function EquilibrateContent(make_polynomial,n_vars,n_params)
pvus = (@polyvar internal___polyvar___u[1:n_vars])[1]
pvps = (@polyvar internal___polyvar___p[1:n_params])[1]
pvt = (@polyvar internal___polyvar___t)[1]
pvus = (HomotopyContinuation.@polyvar internal___polyvar___u[1:n_vars])[1]
pvps = (HomotopyContinuation.@polyvar internal___polyvar___p[1:n_params])[1]
pvt = (HomotopyContinuation.@polyvar internal___polyvar___t)[1]
is_polynomial_system = try
make_polynomial(pvus,pvt,pvps,internal___var___paramvals=fill(1,length(pvps)))
true
Expand Down Expand Up @@ -197,7 +197,7 @@ function make_hc_template!(rn::DiffEqBase.AbstractReactionNetwork)
check_has_polynomial(rn)
p_template = randn(ComplexF64, length(rn.params))
f_template = DynamicPolynomials.subs.(get_equi_poly(rn), Ref(get_polyvars(rn).p => p_template))
solution_template = solutions(HomotopyContinuation.solve(f_template, show_progress=false))
solution_template = HomotopyContinuation.solutions(HomotopyContinuation.solve(f_template, show_progress=false))
reset_hc_templates!(rn)
push_hc_template!(rn,p_template,solution_template)
end
Expand All @@ -213,7 +213,7 @@ end
function add_hc_template!(rn::DiffEqBase.AbstractReactionNetwork)
p_template = randn(ComplexF64, length(rn.params))
f_template = DynamicPolynomials.subs.(get_equi_poly(rn), Ref(get_polyvars(rn).p => p_template))
solution_template = solutions(HomotopyContinuation.solve(f_template, show_progress=false))
solution_template = HomotopyContinuation.solutions(HomotopyContinuation.solve(f_template, show_progress=false))
if (isempty(get_hc_templates(rn))) || (length(solution_template) == length(get_hc_templates(rn)[1]))
push_hc_template!(rn,p_template,solution_template)
elseif length(solution_template) > length(get_hc_templates(rn)[1])
Expand Down Expand Up @@ -272,21 +272,21 @@ end
#Finds steady states of a system using homotopy continuation.
function steady_states(rn::DiffEqBase.AbstractReactionNetwork, solver::HCSteadyStateSolver, p=Vector{Float64}())
using_temp_poly = initialise_solver!(rn,p)
(length(p)==0) && (return positive_real_solutions(solutions(HomotopyContinuation.solve(get_equi_poly(rn),show_progress=false))))
(length(p)==0) && (return positive_real_solutions(HomotopyContinuation.solutions(HomotopyContinuation.solve(get_equi_poly(rn),show_progress=false))))
result = hc_solve_at(rn,p)
using_temp_poly && finalise_solver!(rn)
return positive_real_solutions(result)
end
# Solves the reaction networks at the specific parameter values using Homotopy Continuation.
function hc_solve_at(rn::DiffEqBase.AbstractReactionNetwork,p::Vector{Float64})
for hc_template in get_hc_templates(rn)
result = solutions(HomotopyContinuation.solve(get_equi_poly(rn), hc_template.sol, parameters=get_polyvars(rn).p, p₁=hc_template.p, p₀=p, show_progress=false))
result = HomotopyContinuation.solutions(HomotopyContinuation.solve(get_equi_poly(rn), hc_template.sol, parameters=get_polyvars(rn).p, p₁=hc_template.p, p₀=p, show_progress=false))
(length(result) == length(hc_template.sol)) && (return result)
end
@warn "While solving the system using homotopy continuation some solutions were lost."
best_length = 0; best_solution = [];
for hc_template in get_hc_templates(rn)
result = solutions(HomotopyContinuation.solve(get_equi_poly(rn), hc_template.sol, parameters=get_polyvars(rn).p, p₁=hc_template.p, p₀=p, show_progress=false))
result = HomotopyContinuation.solutions(HomotopyContinuation.solve(get_equi_poly(rn), hc_template.sol, parameters=get_polyvars(rn).p, p₁=hc_template.p, p₀=p, show_progress=false))
(length(result) > best_length) && (best_length = length(result); best_solution = result;)
end
return best_solution
Expand Down Expand Up @@ -560,7 +560,7 @@ function solve_bifurcation(
paths_incomplete = Vector{NamedTuple{(:p, :u),Tuple{Vector{Float64},Vector{Vector{Complex{Float64}}}}}}()
for sol in sol1
path = track_path(sol,tracker1,range...)
if (tracker1.state.status == CoreTrackerStatus.success)
if (tracker1.state.status == HomotopyContinuation.CoreTrackerStatus.success)
remove_sol!(sol2,path.u[end],d_sol)
push!(paths_complete,path)
else
Expand All @@ -569,7 +569,7 @@ function solve_bifurcation(
end
for sol in sol2
path = track_path(sol,tracker2,reverse(range)...)
(tracker2.state.status == CoreTrackerStatus.success) && remove_path!(paths_incomplete,path.u[end],d_sol)
(tracker2.state.status == HomotopyContinuation.CoreTrackerStatus.success) && remove_path!(paths_incomplete,path.u[end],d_sol)
push!(paths_complete,path)
end
append!(paths_complete,paths_incomplete)
Expand Down Expand Up @@ -663,12 +663,12 @@ end
#--- Bifurcation functions used by both solvers ---#
#Makes a coretracker from the given information
function make_coretracker(rn,sol,p1,p2,Δt)
return coretracker(get_equi_poly(rn), sol, parameters=get_polyvars(rn).p, p₁=p1, p₀=p2, max_step_size=Δt)
return HomotopyContinuation.coretracker(get_equi_poly(rn), sol, parameters=get_polyvars(rn).p, p₁=p1, p₀=p2, max_step_size=Δt)
end
#Tracks a path for the given coretracker, returns the path.
function track_path(solution,coretracker,p_start,p_end)
P = Vector{Float64}(); U = Vector{Vector{ComplexF64}}();
for (u,t) in iterator(coretracker, solution)
for (u,t) in HomotopyContinuation.iterator(coretracker, solution)
push!(P,t2p(t,p_start,p_end)); push!(U,u);
end
return (p=P,u=U)
Expand Down

0 comments on commit 923f8d9

Please sign in to comment.