From 8bcc45491bb1015f512d8247b084816bca6f831b Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Mon, 3 Jun 2024 12:01:23 -0400 Subject: [PATCH] Remove ModelingToolkit from dependencies --- Project.toml | 2 - src/SourceCodeMcCormick.jl | 8 +- src/interval/interval.jl | 40 +++---- src/relaxation/relaxation.jl | 84 +++++++------- src/relaxation/rules.jl | 2 +- src/transform/transform.jl | 104 +++++++++--------- src/transform/utilities.jl | 207 +++++++++++++++++------------------ 7 files changed, 219 insertions(+), 228 deletions(-) diff --git a/Project.toml b/Project.toml index 72b888a..3d6feda 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" @@ -18,7 +17,6 @@ CUDA = "5" DocStringExtensions = "0.8 - 0.9" Graphs = "1" IfElse = "0.1.0 - 0.1.1" -ModelingToolkit = "8" SymbolicUtils = "1" Symbolics = "5" julia = "1.6" diff --git a/src/SourceCodeMcCormick.jl b/src/SourceCodeMcCormick.jl index 1e0a828..f2e0706 100644 --- a/src/SourceCodeMcCormick.jl +++ b/src/SourceCodeMcCormick.jl @@ -1,8 +1,6 @@ module SourceCodeMcCormick -# using ModelingToolkit -import ModelingToolkit: Equation, ODESystem, @named, toparam, iv_from_nested_derivative, value, collect_vars! using Symbolics using SymbolicUtils.Code using IfElse @@ -10,9 +8,6 @@ using DocStringExtensions using Graphs using CUDA import Dates - -import SymbolicUtils: BasicSymbolic, exprtype, SYM, TERM, ADD, MUL, POW, DIV - import SymbolicUtils: BasicSymbolic, exprtype, SYM, TERM, ADD, MUL, POW, DIV """ @@ -29,10 +24,9 @@ abstract type AbstractTransform end """ transform_rule(::AbstractTransform, ::Equation) transform_rule(::AbstractTransform, ::Vector{Equation}) - transform_rule(::AbstractTransform, ::ModelingToolkit.ODESystem) transform_rule(::AbstractTransform, ::Num) -Apply the desired transformation to a given expression, equation, set of equations, or ODESystem. +Apply the desired transformation to a given expression, equation, or set of equations. If the `AbstractTransform` is an `IntervalTransform`, all symbols are split into lower and upper bounds, and the number of expressions/equations is initially doubled. `SourceCodeMcCormick` uses the auxiliary variable method, so factorable expressions may be separated into multiple expressions diff --git a/src/interval/interval.jl b/src/interval/interval.jl index f4576dd..b1813f0 100644 --- a/src/interval/interval.jl +++ b/src/interval/interval.jl @@ -30,27 +30,27 @@ function var_names(::IntervalTransform, a::BasicSymbolic) end end -function translate_initial_conditions(::IntervalTransform, prob::ODESystem, new_eqs::Vector{Equation}) - vars, params = extract_terms(new_eqs) - var_defaults = Dict{Any, Any}() - param_defaults = Dict{Any, Any}() +# function translate_initial_conditions(::IntervalTransform, prob::ODESystem, new_eqs::Vector{Equation}) +# vars, params = extract_terms(new_eqs) +# var_defaults = Dict{Any, Any}() +# param_defaults = Dict{Any, Any}() - for (key, val) in prob.defaults - name_lo = String(get_name(key))*"_"*"lo" - name_hi = String(get_name(key))*"_"*"hi" - for i in vars - if in(String(i.f.name), (name_lo, name_hi)) - var_defaults[i] = val - end - end - for i in params - if in(String(i.name), (name_lo, name_hi)) - param_defaults[i] = val - end - end - end - return var_defaults, param_defaults -end +# for (key, val) in prob.defaults +# name_lo = String(get_name(key))*"_"*"lo" +# name_hi = String(get_name(key))*"_"*"hi" +# for i in vars +# if in(String(i.f.name), (name_lo, name_hi)) +# var_defaults[i] = val +# end +# end +# for i in params +# if in(String(i.name), (name_lo, name_hi)) +# param_defaults[i] = val +# end +# end +# end +# return var_defaults, param_defaults +# end """ diff --git a/src/relaxation/relaxation.jl b/src/relaxation/relaxation.jl index c460d6d..9f85cb1 100644 --- a/src/relaxation/relaxation.jl +++ b/src/relaxation/relaxation.jl @@ -34,51 +34,51 @@ function var_names(::McCormickIntervalTransform, a::Any) end -function translate_initial_conditions(::McCormickTransform, prob::ODESystem, new_eqs::Vector{Equation}) - vars, params = extract_terms(new_eqs) - var_defaults = Dict{Any, Any}() - param_defaults = Dict{Any, Any}() +# function translate_initial_conditions(::McCormickTransform, prob::ODESystem, new_eqs::Vector{Equation}) +# vars, params = extract_terms(new_eqs) +# var_defaults = Dict{Any, Any}() +# param_defaults = Dict{Any, Any}() - for (key, val) in prob.defaults - name_cv = String(get_name(key))*"_"*"cv" - name_cc = String(get_name(key))*"_"*"cc" - for i in vars - if String(i.f.name)==name_cv || String(i.f.name)==name_cc - var_defaults[i] = val - end - end - for i in params - if String(i.name)==name_cv || String(i.name)==name_cc - param_defaults[i] = val - end - end - end - return var_defaults, param_defaults -end +# for (key, val) in prob.defaults +# name_cv = String(get_name(key))*"_"*"cv" +# name_cc = String(get_name(key))*"_"*"cc" +# for i in vars +# if String(i.f.name)==name_cv || String(i.f.name)==name_cc +# var_defaults[i] = val +# end +# end +# for i in params +# if String(i.name)==name_cv || String(i.name)==name_cc +# param_defaults[i] = val +# end +# end +# end +# return var_defaults, param_defaults +# end -function translate_initial_conditions(::McCormickIntervalTransform, prob::ODESystem, new_eqs::Vector{Equation}) - vars, params = extract_terms(new_eqs) - var_defaults = Dict{Any, Any}() - param_defaults = Dict{Any, Any}() +# function translate_initial_conditions(::McCormickIntervalTransform, prob::ODESystem, new_eqs::Vector{Equation}) +# vars, params = extract_terms(new_eqs) +# var_defaults = Dict{Any, Any}() +# param_defaults = Dict{Any, Any}() - for (key, val) in prob.defaults - name_lo = String(get_name(key))*"_"*"lo" - name_hi = String(get_name(key))*"_"*"hi" - name_cv = String(get_name(key))*"_"*"cv" - name_cc = String(get_name(key))*"_"*"cc" - for i in vars - if in(String(i.f.name), (name_lo, name_hi, name_cv, name_cc)) - var_defaults[i] = val - end - end - for i in params - if in(String(i.name), (name_lo, name_hi, name_cv, name_cc)) - param_defaults[i] = val - end - end - end - return var_defaults, param_defaults -end +# for (key, val) in prob.defaults +# name_lo = String(get_name(key))*"_"*"lo" +# name_hi = String(get_name(key))*"_"*"hi" +# name_cv = String(get_name(key))*"_"*"cv" +# name_cc = String(get_name(key))*"_"*"cc" +# for i in vars +# if in(String(i.f.name), (name_lo, name_hi, name_cv, name_cc)) +# var_defaults[i] = val +# end +# end +# for i in params +# if in(String(i.name), (name_lo, name_hi, name_cv, name_cc)) +# param_defaults[i] = val +# end +# end +# end +# return var_defaults, param_defaults +# end # A symbolic way of evaluating the line segment between (xL, zL) and (xU, zU) at x (returns IfElse block) line_expr(x, xL, xU, zL, zU) = IfElse.ifelse(zU > zL, (zL*(xU - x) + zU*(x - xL))/(xU - xL), zU) diff --git a/src/relaxation/rules.jl b/src/relaxation/rules.jl index 8f40862..be88d03 100644 --- a/src/relaxation/rules.jl +++ b/src/relaxation/rules.jl @@ -132,7 +132,7 @@ function transform_rule(::McCormickTransform, ::typeof(*), zL, zU, zcv, zcc, xL, # # This returns mult_kernel(-x,-y). Swap x with flipped x and y with flipped # # y, and since both are flipped, no need to worry about negative signs # cv = max(yL*xcc + xL*ycc - xL*yL, yU*xcc + xU*ycc - xU*yU) - # cc = min(yU*xcv + xL*ycv - xL*yU, yL*xcv + xL*ycv - xU*yL) + # cc = min(yU*xcv + xL*ycv - xL*yU, yL*xcv + xU*ycv - xU*yL) # # [x-,ym] diff --git a/src/transform/transform.jl b/src/transform/transform.jl index 9fc4850..7060ddd 100644 --- a/src/transform/transform.jl +++ b/src/transform/transform.jl @@ -3,58 +3,58 @@ include(joinpath(@__DIR__, "binarize.jl")) include(joinpath(@__DIR__, "factor.jl")) -function apply_transform(transform::T, prob::ODESystem; constants::Vector{Num}=Num[]) where T<:AbstractTransform - - # Factorize all model equations to generate a new set of equations - - genparam(get_name(prob.iv.val)) - - equations = Equation[] - for eqn in prob.eqs - current = length(equations) - factor(eqn.rhs, eqs=equations) - if length(equations) > current - push!(equations, Equation(eqn.lhs, equations[end].rhs)) - deleteat!(equations, length(equations)-1) - else - index = findall(x -> isequal(x.rhs, eqn.rhs), equations) - push!(equations, Equation(eqn.lhs, equations[index[1]].lhs)) - end - end - - # Apply transform rules to the factored equations to make the final equation set - new_equations = Equation[] - for a in equations - zn = var_names(transform, zstr(a)) - if string(xstr(a)) in string.(constants) - xn = (xstr(a), xstr(a), xstr(a), xstr(a)) - else - xn = var_names(transform, xstr(a)) - end - if isone(arity(a)) - targs = (transform, op(a), zn..., xn...) - else - if string(ystr(a)) in string.(constants) - yn = (ystr(a), ystr(a), ystr(a), ystr(a)) - else - yn = var_names(transform, ystr(a)) - end - targs = (transform, op(a), zn..., xn..., yn...) - end - new = transform_rule(targs...) - for i in new - push!(new_equations, i) - end - end - - # Copy model start points to the newly transformed variables - var_defaults, param_defaults = translate_initial_conditions(transform, prob, new_equations) - - # Use the transformed equations and new start points to generate a new ODE system - @named new_sys = ODESystem(new_equations, defaults=merge(var_defaults, param_defaults)) - - return new_sys -end +# function apply_transform(transform::T, prob::ODESystem; constants::Vector{Num}=Num[]) where T<:AbstractTransform + +# # Factorize all model equations to generate a new set of equations + +# genparam(get_name(prob.iv.val)) + +# equations = Equation[] +# for eqn in prob.eqs +# current = length(equations) +# factor(eqn.rhs, eqs=equations) +# if length(equations) > current +# push!(equations, Equation(eqn.lhs, equations[end].rhs)) +# deleteat!(equations, length(equations)-1) +# else +# index = findall(x -> isequal(x.rhs, eqn.rhs), equations) +# push!(equations, Equation(eqn.lhs, equations[index[1]].lhs)) +# end +# end + +# # Apply transform rules to the factored equations to make the final equation set +# new_equations = Equation[] +# for a in equations +# zn = var_names(transform, zstr(a)) +# if string(xstr(a)) in string.(constants) +# xn = (xstr(a), xstr(a), xstr(a), xstr(a)) +# else +# xn = var_names(transform, xstr(a)) +# end +# if isone(arity(a)) +# targs = (transform, op(a), zn..., xn...) +# else +# if string(ystr(a)) in string.(constants) +# yn = (ystr(a), ystr(a), ystr(a), ystr(a)) +# else +# yn = var_names(transform, ystr(a)) +# end +# targs = (transform, op(a), zn..., xn..., yn...) +# end +# new = transform_rule(targs...) +# for i in new +# push!(new_equations, i) +# end +# end + +# # Copy model start points to the newly transformed variables +# var_defaults, param_defaults = translate_initial_conditions(transform, prob, new_equations) + +# # Use the transformed equations and new start points to generate a new ODE system +# @named new_sys = ODESystem(new_equations, defaults=merge(var_defaults, param_defaults)) + +# return new_sys +# end function apply_transform(transform::T, eqn_vector::Vector{Equation}; constants::Vector{Num}=Num[]) where T<:AbstractTransform diff --git a/src/transform/utilities.jl b/src/transform/utilities.jl index ec7be73..3fadc08 100644 --- a/src/transform/utilities.jl +++ b/src/transform/utilities.jl @@ -100,7 +100,7 @@ genvar(a::Symbol) = genparam(a) function genparam(a::Symbol) params = Symbol[] ex = Expr(:block) - var_name, expr = Symbolics.construct_vars(:parameters, a, Real, nothing, nothing, nothing, toparam, false) + var_name, expr = Symbolics.construct_vars(:parameters, a, Real, nothing, nothing, nothing, identity, false) push!(params, var_name) push!(ex.args, expr) rhs = Symbolics.build_expr(:vect, params) @@ -109,108 +109,108 @@ function genparam(a::Symbol) end -# A function to extract terms from a set of equations, for use in dynamic systems -function extract_terms(eqs::Vector{Equation}) - allstates = SymbolicUtils.OrderedSet() - ps = SymbolicUtils.OrderedSet() - for eq in eqs - if ~(eq.lhs isa Number) - iv = iv_from_nested_derivative(eq.lhs) - break - end - end - iv = value(iv) - for eq in eqs - eq.lhs isa Union{SymbolicUtils.Symbolic,Number} || (push!(compressed_eqs, eq); continue) - collect_vars!(allstates, ps, eq.lhs, iv) - collect_vars!(allstates, ps, eq.rhs, iv) - end - - return allstates, ps -end - -# Interprets existing start points in an ODESystem, then applies the provided bounds to the -# given term and adds to (or overwrites) existing start points. Returns an updated ODESystem -function set_bounds(sys::ODESystem, term::Num, bounds::Tuple{Float64, Float64}) - base_name = get_name(Symbolics.value(term)) - name_lo = String(base_name)*"_"*"lo" - name_hi = String(base_name)*"_"*"hi" - - model_terms = Vector{Union{Term,Sym}}() - for i in sys.states - push!(model_terms, Symbolics.value(i)) - end - for i in sys.ps - push!(model_terms, Symbolics.value(i)) - end - real_lo = nothing - real_hi = nothing - for i in model_terms - if String(get_name(i))==name_lo - real_lo = i - elseif String(get_name(i))==name_hi - real_hi = i - end - end - if real_lo in keys(sys.defaults) - delete!(sys.defaults, real_lo) - sys.defaults[real_lo] = bounds[1] - else - sys.defaults[real_lo] = bounds[1] - end - if real_hi in keys(sys.defaults) - delete!(sys.defaults, real_hi) - sys.defaults[real_hi] = bounds[2] - else - sys.defaults[real_hi] = bounds[2] - end - return sys -end - -function set_bounds(sys::ODESystem, terms::Vector{Num}, bounds::Vector{Tuple{Float64, Float64}}) - for i in 1:length(terms) - sys = set_bounds(sys, terms[i], bounds[i]) - end - return sys -end - -function get_cvcc_start_dict(sys::ODESystem, term::Num, start_point::Float64) - base_name = get_name(Symbolics.value(term)) - name_cv = String(base_name)*"_"*"cv" - name_cc = String(base_name)*"_"*"cc" - - model_terms = Vector{Union{Term,Sym}}() - for i in sys.states - push!(model_terms, Symbolics.value(i)) - end - for i in sys.ps - push!(model_terms, Symbolics.value(i)) - end - real_cv = nothing - real_cc = nothing - for i in model_terms - if String(get_name(i))==name_cv - real_cv = i - elseif String(get_name(i))==name_cc - real_cc = i - end - end - - new_dict = copy(sys.defaults) - if real_cv in keys(new_dict) - delete!(new_dict, real_cv) - new_dict[real_cv] = start_point - else - new_dict[real_cv] = start_point - end - if real_cc in keys(new_dict) - delete!(new_dict, real_cc) - new_dict[real_cc] = start_point - else - new_dict[real_cc] = start_point - end - return new_dict -end +# # A function to extract terms from a set of equations, for use in dynamic systems +# function extract_terms(eqs::Vector{Equation}) +# allstates = SymbolicUtils.OrderedSet() +# ps = SymbolicUtils.OrderedSet() +# for eq in eqs +# if ~(eq.lhs isa Number) +# iv = iv_from_nested_derivative(eq.lhs) +# break +# end +# end +# iv = Symbolics.value(iv) +# for eq in eqs +# eq.lhs isa Union{SymbolicUtils.Symbolic,Number} || (push!(compressed_eqs, eq); continue) +# collect_vars!(allstates, ps, eq.lhs, iv) +# collect_vars!(allstates, ps, eq.rhs, iv) +# end + +# return allstates, ps +# end + +# # Interprets existing start points in an ODESystem, then applies the provided bounds to the +# # given term and adds to (or overwrites) existing start points. Returns an updated ODESystem +# function set_bounds(sys::ODESystem, term::Num, bounds::Tuple{Float64, Float64}) +# base_name = get_name(Symbolics.value(term)) +# name_lo = String(base_name)*"_"*"lo" +# name_hi = String(base_name)*"_"*"hi" + +# model_terms = Vector{Union{Term,Sym}}() +# for i in sys.states +# push!(model_terms, Symbolics.value(i)) +# end +# for i in sys.ps +# push!(model_terms, Symbolics.value(i)) +# end +# real_lo = nothing +# real_hi = nothing +# for i in model_terms +# if String(get_name(i))==name_lo +# real_lo = i +# elseif String(get_name(i))==name_hi +# real_hi = i +# end +# end +# if real_lo in keys(sys.defaults) +# delete!(sys.defaults, real_lo) +# sys.defaults[real_lo] = bounds[1] +# else +# sys.defaults[real_lo] = bounds[1] +# end +# if real_hi in keys(sys.defaults) +# delete!(sys.defaults, real_hi) +# sys.defaults[real_hi] = bounds[2] +# else +# sys.defaults[real_hi] = bounds[2] +# end +# return sys +# end + +# function set_bounds(sys::ODESystem, terms::Vector{Num}, bounds::Vector{Tuple{Float64, Float64}}) +# for i in 1:length(terms) +# sys = set_bounds(sys, terms[i], bounds[i]) +# end +# return sys +# end + +# function get_cvcc_start_dict(sys::ODESystem, term::Num, start_point::Float64) +# base_name = get_name(Symbolics.value(term)) +# name_cv = String(base_name)*"_"*"cv" +# name_cc = String(base_name)*"_"*"cc" + +# model_terms = Vector{Union{Term,Sym}}() +# for i in sys.states +# push!(model_terms, Symbolics.value(i)) +# end +# for i in sys.ps +# push!(model_terms, Symbolics.value(i)) +# end +# real_cv = nothing +# real_cc = nothing +# for i in model_terms +# if String(get_name(i))==name_cv +# real_cv = i +# elseif String(get_name(i))==name_cc +# real_cc = i +# end +# end + +# new_dict = copy(sys.defaults) +# if real_cv in keys(new_dict) +# delete!(new_dict, real_cv) +# new_dict[real_cv] = start_point +# else +# new_dict[real_cv] = start_point +# end +# if real_cc in keys(new_dict) +# delete!(new_dict, real_cc) +# new_dict[real_cc] = start_point +# else +# new_dict[real_cc] = start_point +# end +# return new_dict +# end """ @@ -337,7 +337,6 @@ function sort_vars(strings::Vector{String}) end end end - return vars, strings end # Step 3) Attach simplified variable names to each element, and then