Consistency and QOL Fixes

--Generalized extension name into the abstract type ExtendGPU to allow for alternative GPU-based methods.
--Added manual garbage collection to reduce time spent on GPU GC.
--Corrected variable name generation when variables with indices are used (e.g. x[1:5] instead of x).
--Corrected recognition of variables with indices in some edge cases.
--Corrected mid_expr, which sometimes returned results inconsistent with McCormick.jl.
--Added a ::Nothing relaxation rule.
--Corrected some typos in multiplication relaxation rules.
--Added "cut" operator to multiplication and division operations.
--Added min/max rules (may need further development/verification).
--Changed "factor" to no longer modify input expressions.
--Added some support for "constants", which are inputs that should be treated symbolically but should not be expanded into McCormick tuples.
--Added "sort_vars" function, which groups McCormick tuples together, and places constants at the front, before alphabetically sorting McCormick tuple groups.
--Reordered user-facing inputs and outputs to [cv, cc, lo, hi] instead of [lo, hi, cv, cc].

Author: RXGottlieb

.gitignore

@@ -0,0 +1,2 @@
+.vscode/
+Manifest.toml

Project.toml

@@ -1,18 +1,22 @@
 name = "SourceCodeMcCormick"
 uuid = "a7283dc5-4ecf-47fb-a95b-1412723fc960"
 authors = ["Robert Gottlieb <Robert.x.gottlieb@uconn.edu>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+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"
 
 [compat]
 DocStringExtensions = "0.8 - 0.9"
-ModelingToolkit = "8"
 IfElse = "0.1.0 - 0.1.1"
+ModelingToolkit = "8"
 SymbolicUtils = "1"
 julia = "1.7"
 

examples/ParBB/extension.jl

@@ -5,19 +5,22 @@
 # See also: subroutines.jl, in this folder
 
 using DocStringExtensions, EAGO
+
+abstract type ExtendGPU <: EAGO.ExtensionType end
+
 """
 $(TYPEDEF)
 
-The ExtendGPU integrator is meant to be paired with the SourceCodeMcCormick
-package. A required component of ExtendGPU is the function `convex_func`, 
+The PointwiseGPU integrator is meant to be paired with the SourceCodeMcCormick
+package. A required component of PointwiseGPU is the function `convex_func`, 
 which should take arguments corresponding to the McCormick tuple [cc, cv, hi, lo]
 for each branch variable in the problem and return a vector of convex
 relaxation evaluations of the objective function, of length equal to the
 length of the inputs.
 
 $(TYPEDFIELDS)
 """
-Base.@kwdef mutable struct ExtendGPU <: EAGO.ExtensionType
+Base.@kwdef mutable struct PointwiseGPU <: ExtendGPU
     "A user-defined function taking argument `p` and returning a vector
     of convex evaluations of the objective function"
     convex_func
@@ -49,9 +52,9 @@ Base.@kwdef mutable struct ExtendGPU <: EAGO.ExtensionType
     multistart_points::Int = 1
 end
 
-function ExtendGPU(convex_func, var_count::Int; alpha::Float64 = 0.01, node_limit::Int = 50000, 
+function PointwiseGPU(convex_func, var_count::Int; alpha::Float64 = 0.01, node_limit::Int = 50000, 
                     prepopulate::Bool = true, multistart_points::Int = 1)
-    return ExtendGPU(convex_func, var_count, node_limit, alpha, 
+    return PointwiseGPU(convex_func, var_count, node_limit, alpha, 
                     Vector{Float64}(undef, node_limit), Vector{Float64}(undef, node_limit), Vector{NodeBB}(undef, node_limit), 0,
                     Matrix{Float64}(undef, node_limit, var_count),
                     Matrix{Float64}(undef, node_limit, var_count), prepopulate, multistart_points)

examples/ParBB/subroutines.jl

@@ -50,6 +50,14 @@ function solve_gpu!(m::EAGO.GlobalOptimizer)
     # Run branch and bound; terminate when the stack is empty or when some
     # tolerance or limit is hit
     while !EAGO.termination_check(m)
+        
+        # Garbage collect every 1000 iterations
+        if mod(m._iteration_count, 1000)==0
+            GC.enable(true)
+            GC.gc(false)
+            GC.enable(false)
+        end
+
         # Fathoming step
         EAGO.fathom!(m)
 
@@ -197,7 +205,7 @@ add_to_substack!(m::EAGO.GlobalOptimizer{R,S,Q}) where {R,S,Q<:EAGO.ExtensionTyp
 
 # Solve the lower and upper problem for all nodes simultaneously, using the convex_func
 # function from the ExtendGPU extension
-function lower_and_upper_problem!(t::ExtendGPU, m::EAGO.GlobalOptimizer)
+function lower_and_upper_problem!(t::PointwiseGPU, m::EAGO.GlobalOptimizer)
     # Step 1) Bring the bounds into the GPU
     lvbs_d = CuArray(t.all_lvbs) 
     uvbs_d = CuArray(t.all_uvbs) # [points x num_vars]
@@ -208,8 +216,8 @@ function lower_and_upper_problem!(t::ExtendGPU, m::EAGO.GlobalOptimizer)
     eval_points = Vector{CuArray{Float64}}(undef, 3*w)
     for i = 1:w
         eval_points[3i-2] = CuArray{Float64}(undef, l*np)
-        eval_points[3i-1] = repeat(uvbs_d[:,i], inner=np)
-        eval_points[3i] = repeat(lvbs_d[:,i], inner=np)
+        eval_points[3i-1] = repeat(lvbs_d[:,i], inner=np)
+        eval_points[3i] = repeat(uvbs_d[:,i], inner=np)
     end
     evals_d = CuArray{Float64}(undef, l*np)
     results_d = CuArray{Float64}(undef, l)

examples/explicit_kinetic_problem.jl

@@ -7,7 +7,7 @@
 #############################################################################
 
 # Import the necessary packages
-using EAGO, JuMP, CSV, DataFrames, Symbolics, SourceCodeMcCormick, CUDA, BenchmarkTools
+using EAGO, JuMP, CSV, DataFrames, Symbolics, SourceCodeMcCormick, CUDA
 
 # Import the ParBB algorithm
 include(joinpath(@__DIR__, "ParBB", "extension.jl"))
@@ -255,15 +255,3 @@ JuMP.register(m, :fobj, 3, fobj, autodiff=true)
 @NLobjective(m, Min, fobj(p...))
 
 optimize!(m)
-
-# Save the results
-open("C:/Users/rxg20001/Documents/Research Files/2023-03-25 Kinetic Problem SCMC GPU.txt","a") do io
-    println(io, "SCMC; GPU; Alpha=2e-5; 50k nodes per batch")
-    println(io, "Semi-parallelized lvb/uvb loading")
-    println(io, "Node count, Global Lower Bound, Global Upper Bound, Time (s)")
-    log = m.moi_backend.optimizer.model._global_optimizer._log
-    loglen = length(log.global_lower_bound)
-    for i = 1:loglen
-        println(io, string(log.node_count[i]), ",", string(log.global_lower_bound[i]), ",", string(log.global_upper_bound[i]), ",", string(log.run_time[i]))
-    end
-end

src/SourceCodeMcCormick.jl

@@ -1,10 +1,17 @@
 
 module SourceCodeMcCormick
 
-using ModelingToolkit
+# using ModelingToolkit
+import ModelingToolkit: Equation, ODESystem, @named, toparam, iv_from_nested_derivative, value, collect_vars!
+using Symbolics
 using SymbolicUtils.Code
 using IfElse
 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
 

src/interval/interval.jl

@@ -12,7 +12,7 @@ function var_names(::IntervalTransform, a::BasicSymbolic)
         aU = genvar(Symbol(string(get_name(a))*"_hi"))
         return aL.val, aU.val
     elseif exprtype(a)==TERM
-        if varterm(a)
+        if varterm(a) && typeof(a.f)<:BasicSymbolic
             arg_list = Symbol[]
             for i in a.arguments
                 push!(arg_list, get_name(i))
@@ -56,25 +56,29 @@ end
 """
     get_name
 
-Take a `BasicSymbolic` object such as `x[1,1]` and return a symbol like `:x_1_1`.
+Take a `BasicSymbolic` object such as `x[1,1]` and return a symbol like `:xv1v1`.
 """
 function get_name(a::BasicSymbolic)
     if exprtype(a)==SYM
         return a.name
     elseif exprtype(a)==TERM
         if varterm(a)
-            return a.f.name
-        elseif (a.f==getindex)
-            args = a.arguments
-            new_var = string(args[1])
-            for i in 2:lastindex(args)
-                new_var = new_var * "_" * string(args[i])
+            if a.f==getindex
+                args = a.arguments
+                new_var = string(args[1])
+                for i in 2:lastindex(args)
+                    new_var = new_var * "v" * string(args[i])
+                end
+                return Symbol(new_var)
+            else
+                return a.f.name
             end
-            return Symbol(new_var)
         else
             error("Problem generating variable name. This may happen if the variable is non-standard. Please post an issue if you get this error.")
         end
     end
 end
+get_name(a::Num) = get_name(a.val)
+get_name(a::Real) = a
 
 include(joinpath(@__DIR__, "rules.jl"))

src/relaxation/relaxation.jl

@@ -9,7 +9,7 @@ function var_names(::McCormickTransform, a::BasicSymbolic)
         acc = genvar(Symbol(string(get_name(a))*"_cc"))
         return acv.val, acc.val
     elseif exprtype(a)==TERM
-        if varterm(a)
+        if varterm(a) && typeof(a.f)<:BasicSymbolic
             arg_list = Symbol[]
             for i in a.arguments
                 push!(arg_list, get_name(i))
@@ -84,7 +84,7 @@ end
 line_expr(x, xL, xU, zL, zU) = IfElse.ifelse(zU > zL, (zL*(xU - x) + zU*(x - xL))/(xU - xL), zU)
 
 # A symbolic way of computing the mid of three numbers (returns IfElse block)
-mid_expr(a, b, c) = IfElse.ifelse(a < b, IfElse.ifelse(b < c, b, IfElse.ifelse(c < a, a, c)),
-                        IfElse.ifelse(c < b, b, IfElse.ifelse(a < c, a, c)))
+mid_expr(x, y, z) = IfElse.ifelse(x >= y, IfElse.ifelse(y >= z, y, IfElse.ifelse(y == x, y, IfElse.ifelse(z >= x, x, z))),
+        IfElse.ifelse(z >= y, y, IfElse.ifelse(x >= z, x, z)))
 
 include(joinpath(@__DIR__, "rules.jl"))