From e402b95c88e86ff5695458d63b3e704394ebf98a Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 12 Aug 2024 11:49:23 +1200 Subject: [PATCH 1/2] Remove unnecessary Alp module reference --- src/Alpine.jl | 1 - src/MOI_wrapper/MOI_wrapper.jl | 12 +- src/bounding_model.jl | 125 +++++++++---------- src/embedding.jl | 21 ++-- src/heuristics.jl | 4 +- src/log.jl | 79 ++++++------ src/main_algorithm.jl | 218 ++++++++++++++++----------------- src/multilinear.jl | 121 ++++++++---------- src/nlexpr.jl | 115 +++++++++-------- src/operators.jl | 144 +++++++++++----------- src/presolve.jl | 130 ++++++++++---------- src/utility.jl | 80 ++++++------ src/variable_bounds.jl | 56 ++++----- 13 files changed, 530 insertions(+), 576 deletions(-) diff --git a/src/Alpine.jl b/src/Alpine.jl index 2f19318c..325c8752 100644 --- a/src/Alpine.jl +++ b/src/Alpine.jl @@ -10,7 +10,6 @@ import Combinatorics import Pkg const ALPINE_DEBUG = false -const Alp = Alpine include("const.jl") diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 8cd54ff6..ed0103af 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -96,7 +96,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer bound_sol_pool::Dict{Any,Any} # A pool of solutions from solving model_mip # Linking constraints info for Multilinear terms - linking_constraints_info::Union{Nothing,Dict{Any,Any}} # Stored multilinear linking constraints info + linking_constraints_info::Union{Nothing,Dict{Any,Any}} # Stored multilinear linking constraints info # Logging information and status logs::Dict{Symbol,Any} # Logging information @@ -108,7 +108,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer # Constructor for Alpine.Optimizer function Optimizer() m = new() - m.options = Alp.get_default_options() + m.options = get_default_options() MOI.empty!(m) return m end @@ -254,10 +254,10 @@ MOI.get(::Optimizer, ::MOI.SolverName) = "Alpine" MOI.get(::Optimizer, ::MOI.SolverVersion) = _ALPINE_VERSION function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value) - return Alp.set_option(model, Symbol(param.name), value) + return set_option(model, Symbol(param.name), value) end function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute) - return Alp.get_option(model, Symbol(param.name)) + return get_option(model, Symbol(param.name)) end function MOI.add_variables(model::Optimizer, n::Int) @@ -395,10 +395,10 @@ is_max_sense(model::Optimizer) = model.sense_orig == MOI.MAX_SENSE function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense) model.sense_orig = sense - if Alp.is_max_sense(model) + if is_max_sense(model) model.best_obj = -Inf model.best_bound = Inf - elseif Alp.is_min_sense(model) + elseif is_min_sense(model) model.best_obj = Inf model.best_bound = -Inf else diff --git a/src/bounding_model.jl b/src/bounding_model.jl index 191e9115..fd041ea1 100644 --- a/src/bounding_model.jl +++ b/src/bounding_model.jl @@ -3,10 +3,10 @@ Set up a MILP bounding model base on variable domain partitioning information stored in `use_disc`. By default, if `use_disc` is not provided, it will use `m.discretizations` store in the Alpine model. -The basic idea of this MILP bounding model is to use piecewise polyhedral/convex relaxations to tighten the -basic relaxations of the original non-convex region. Among all presented partitions, the bounding model -will choose one specific partition as the lower bound solution. The more partitions there are, the -better or finer bounding model relax the original MINLP while the more efforts required to solve +The basic idea of this MILP bounding model is to use piecewise polyhedral/convex relaxations to tighten the +basic relaxations of the original non-convex region. Among all presented partitions, the bounding model +will choose one specific partition as the lower bound solution. The more partitions there are, the +better or finer bounding model relax the original MINLP while the more efforts required to solve this MILP is required. """ function create_bounding_mip(m::Optimizer; use_disc = nothing) @@ -14,27 +14,27 @@ function create_bounding_mip(m::Optimizer; use_disc = nothing) if (m.logs[:n_iter] == 1) && (m.status[:local_solve] in union(STATUS_OPT, STATUS_LIMIT, STATUS_WARM_START)) # Setting up an initial partition - Alp.add_partition(m, use_solution = m.best_sol) + add_partition(m, use_solution = m.best_sol) elseif m.logs[:n_iter] >= 2 # Add subsequent iteration partitions - Alp.add_partition(m) + add_partition(m) end discretization = m.discretization else discretization = use_disc end - m.model_mip = JuMP.Model(Alp.get_option(m, :mip_solver)) # Construct JuMP Model + m.model_mip = JuMP.Model(get_option(m, :mip_solver)) # Construct JuMP Model start_build = time() # ------- Model Construction ------ # - Alp.amp_post_vars(m) # Post original and lifted variables - Alp.amp_post_lifted_constraints(m) # Post original and lifted constraints - Alp.amp_post_lifted_objective(m) # Post original and/or objective - Alp.amp_post_convexification(m, use_disc = discretization) # Post convexification constraints + amp_post_vars(m) # Post original and lifted variables + amp_post_lifted_constraints(m) # Post original and lifted constraints + amp_post_lifted_objective(m) # Post original and/or objective + amp_post_convexification(m, use_disc = discretization) # Post convexification constraints # --------------------------------- # cputime_build = time() - start_build m.logs[:total_time] += cputime_build - m.logs[:time_left] = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + m.logs[:time_left] = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) return end @@ -48,8 +48,8 @@ to finish the last step required during the construction of bounding model. function amp_post_convexification(m::Optimizer; use_disc = nothing) use_disc === nothing ? discretization = m.discretization : discretization = use_disc - Alp.amp_post_convhull(m, use_disc = discretization) # convex hull representation - Alp._is_fully_convexified(m) # Ensure if all the non-linear terms are convexified + amp_post_convhull(m, use_disc = discretization) # convex hull representation + _is_fully_convexified(m) # Ensure if all the non-linear terms are convexified return end @@ -101,16 +101,16 @@ end function amp_post_lifted_constraints(m::Optimizer) for i in 1:m.num_constr_orig if m.constr_structure[i] == :affine - Alp.amp_post_affine_constraint(m.model_mip, m.bounding_constr_mip[i]) + amp_post_affine_constraint(m.model_mip, m.bounding_constr_mip[i]) elseif m.constr_structure[i] == :convex - Alp.amp_post_convex_constraint(m.model_mip, m.bounding_constr_mip[i]) + amp_post_convex_constraint(m.model_mip, m.bounding_constr_mip[i]) else error("Unknown constr_structure type $(m.constr_structure[i])") end end for i in keys(m.linear_terms) - Alp.amp_post_linear_lift_constraints(m.model_mip, m.linear_terms[i]) + amp_post_linear_lift_constraints(m.model_mip, m.linear_terms[i]) end return @@ -234,21 +234,18 @@ function add_partition(m::Optimizer; kwargs...) haskey(options, :use_solution) ? point_vec = options[:use_solution] : point_vec = m.best_bound_sol - if isa(Alp.get_option(m, :disc_add_partition_method), Function) - # m.discretization = eval(Alp.get_option(m, :disc_add_partition_method))(m, use_disc=discretization, use_solution=point_vec) - m.discretization = Alp.get_option(m, :disc_add_partition_method)( + if isa(get_option(m, :disc_add_partition_method), Function) + # m.discretization = eval(get_option(m, :disc_add_partition_method))(m, use_disc=discretization, use_solution=point_vec) + m.discretization = get_option(m, :disc_add_partition_method)( m, use_disc = discretization, use_solution = point_vec, ) - elseif Alp.get_option(m, :disc_add_partition_method) == "adaptive" - m.discretization = Alp.add_adaptive_partition( - m, - use_disc = discretization, - use_solution = point_vec, - ) - elseif Alp.get_option(m, :disc_add_partition_method) == "uniform" - m.discretization = Alp.add_uniform_partition(m, use_disc = discretization) + elseif get_option(m, :disc_add_partition_method) == "adaptive" + m.discretization = + add_adaptive_partition(m, use_disc = discretization, use_solution = point_vec) + elseif get_option(m, :disc_add_partition_method) == "uniform" + m.discretization = add_uniform_partition(m, use_disc = discretization) else error("Unknown input on how to add partitions.") end @@ -263,7 +260,7 @@ end A built-in method used to add a new partition on feasible domains of variables chosen for partitioning. -This can be illustrated by the following example. Let the previous iteration's partition vector on +This can be illustrated by the following example. Let the previous iteration's partition vector on variable "x" be given by [0, 3, 7, 9]. And say, the lower bounding solution has a value of 4 for variable "x". In the case when `partition_scaling_factor = 4`, this function creates the new partition vector as follows: [0, 3, 3.5, 4, 4.5, 7, 9] @@ -283,11 +280,11 @@ function add_adaptive_partition(m::Optimizer; kwargs...) haskey(options, :use_solution) ? point_vec = copy(options[:use_solution]) : point_vec = copy(m.best_bound_sol) haskey(options, :use_ratio) ? ratio = options[:use_ratio] : - ratio = Alp.get_option(m, :partition_scaling_factor) + ratio = get_option(m, :partition_scaling_factor) haskey(options, :branching) ? branching = options[:branching] : branching = false if length(point_vec) < m.num_var_orig + m.num_var_linear_mip + m.num_var_nonlinear_mip - point_vec = Alp.resolve_lifted_var_value(m, point_vec) # Update the solution vector for lifted variable + point_vec = resolve_lifted_var_value(m, point_vec) # Update the solution vector for lifted variable end if branching @@ -304,11 +301,11 @@ function add_adaptive_partition(m::Optimizer; kwargs...) # Built-in method based-on variable type if m.var_type[i] == :Cont i in processed && continue - point = Alp.correct_point(m, discretization[i], point, i) + point = correct_point(m, discretization[i], point, i) for j in 1:λCnt if point >= discretization[i][j] && point <= discretization[i][j+1] # Locating the right location - radius = Alp.calculate_radius(discretization[i], j, ratio) - Alp.insert_partition_helper(m, i, j, point, radius, discretization[i]) + radius = calculate_radius(discretization[i], j, ratio) + insert_partition_helper(m, i, j, point, radius, discretization[i]) push!(processed, i) break end @@ -332,7 +329,7 @@ end This function targets to address unexpected numerical issues when adding partitions in tight regions. """ function correct_point(m::Optimizer, partvec::Vector, point::Float64, var::Int) - tol = Alp.get_option(m, :tol) + tol = get_option(m, :tol) if point < partvec[1] - tol || point > partvec[end] + tol @warn " Warning: VAR$(var) SOL=$(point) out of discretization [$(partvec[1]),$(partvec[end])]. Hence, taking the middle point." @@ -370,7 +367,7 @@ function insert_partition_helper( partvec::Vector, ) abstol, reltol = - Alp.get_option(m, :disc_abs_width_tol), Alp.get_option(m, :disc_rel_width_tol) + get_option(m, :disc_abs_width_tol), get_option(m, :disc_rel_width_tol) lb_local, ub_local = partvec[partidx], partvec[partidx+1] ub_touch, lb_touch = true, true @@ -390,28 +387,28 @@ function insert_partition_helper( lb_touch = false end - if (ub_touch && lb_touch) || Alp.check_solution_history(m, var) + if (ub_touch && lb_touch) || check_solution_history(m, var) distvec = [(j, partvec[j+1] - partvec[j]) for j in 1:length(partvec)-1] sort!(distvec, by = x -> x[2]) point_orig = point pos = distvec[end][1] lb_local = partvec[pos] ub_local = partvec[pos+1] - isapprox(lb_local, ub_local; atol = Alp.get_option(m, :tol)) && return - chunk = (ub_local - lb_local) / Alp.get_option(m, :disc_divert_chunks) - point = lb_local + (ub_local - lb_local) / Alp.get_option(m, :disc_divert_chunks) - for i in 2:Alp.get_option(m, :disc_divert_chunks) + isapprox(lb_local, ub_local; atol = get_option(m, :tol)) && return + chunk = (ub_local - lb_local) / get_option(m, :disc_divert_chunks) + point = lb_local + (ub_local - lb_local) / get_option(m, :disc_divert_chunks) + for i in 2:get_option(m, :disc_divert_chunks) insert!( partvec, pos + 1, - lb_local + chunk * (Alp.get_option(m, :disc_divert_chunks) - (i - 1)), + lb_local + chunk * (get_option(m, :disc_divert_chunks) - (i - 1)), ) end - (Alp.get_option(m, :log_level) > 199) && println( - "[DEBUG] !D! VAR$(var): SOL=$(round(point_orig; digits = 4))=>$(point) |$(round(lb_local; digits = 4)) | $(Alp.get_option(m, :disc_divert_chunks)) SEGMENTS | $(round(ub_local; digits = 4))|", + (get_option(m, :log_level) > 199) && println( + "[DEBUG] !D! VAR$(var): SOL=$(round(point_orig; digits = 4))=>$(point) |$(round(lb_local; digits = 4)) | $(get_option(m, :disc_divert_chunks)) SEGMENTS | $(round(ub_local; digits = 4))|", ) else - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[DEBUG] VAR$(var): SOL=$(round(point; digits = 4)) RADIUS=$(radius), PARTITIONS=$(length(partvec)-1) |$(round(lb_local; digits = 4)) |$(round(lb_new; digits = 6)) <- * -> $(round(ub_new; digits = 6))| $(round(ub_local; digits = 4))|", ) end @@ -428,11 +425,11 @@ function add_uniform_partition(m::Optimizer; kwargs...) lb_local = discretization[i][1] ub_local = discretization[i][end] distance = ub_local - lb_local - # chunk = distance / ((m.logs[:n_iter] + 1) * Alp.get_option(m, :disc_uniform_rate)) - chunk = distance / (m.logs[:n_iter] * Alp.get_option(m, :disc_uniform_rate)) + # chunk = distance / ((m.logs[:n_iter] + 1) * get_option(m, :disc_uniform_rate)) + chunk = distance / (m.logs[:n_iter] * get_option(m, :disc_uniform_rate)) discretization[i] = [ lb_local + chunk * (j - 1) for - j in 1:(m.logs[:n_iter])*Alp.get_option(m, :disc_uniform_rate) + j in 1:(m.logs[:n_iter])*get_option(m, :disc_uniform_rate) ] push!(discretization[i], ub_local) # Safety Scheme end @@ -441,20 +438,20 @@ function add_uniform_partition(m::Optimizer; kwargs...) end function update_partition_scaling_factor(m::Optimizer, presolve = false) - m.logs[:n_iter] > 2 && return Alp.get_option(m, :partition_scaling_factor) # Stop branching after the second iterations + m.logs[:n_iter] > 2 && return get_option(m, :partition_scaling_factor) # Stop branching after the second iterations ratio_pool = [8:2:20;] # Built-in try range convertor = Dict(MOI.MAX_SENSE => :<, MOI.MIN_SENSE => :>) # revconvertor = Dict(MOI.MAX_SENSE => :>, MOI.MIN_SENSE => :<) incumb_ratio = ratio_pool[1] - Alp.is_min_sense(m) ? incumb_res = -Inf : incumb_res = Inf + is_min_sense(m) ? incumb_res = -Inf : incumb_res = Inf res_collector = Float64[] for r in ratio_pool st = time() if !isempty(m.best_sol) - branch_disc = Alp.add_adaptive_partition( + branch_disc = add_adaptive_partition( m, use_disc = m.discretization, branching = true, @@ -462,15 +459,15 @@ function update_partition_scaling_factor(m::Optimizer, presolve = false) use_solution = m.best_sol, ) else - branch_disc = Alp.add_adaptive_partition( + branch_disc = add_adaptive_partition( m, use_disc = m.discretization, branching = true, use_ratio = r, ) end - Alp.create_bounding_mip(m, use_disc = branch_disc) - res = Alp.disc_branch_solve(m) + create_bounding_mip(m, use_disc = branch_disc) + res = disc_branch_solve(m) push!(res_collector, res) if eval(convertor[m.sense_orig])(res, incumb_res) # && abs(abs(collector[end]-res)/collector[end]) > 1e-1 # %1 of difference incumb_res = res @@ -481,18 +478,18 @@ function update_partition_scaling_factor(m::Optimizer, presolve = false) println("Expensive disc branching pass... Fixed at 8") return 8 end - Alp.get_option(m, :log_level) > 0 && + get_option(m, :log_level) > 0 && println("BRANCH RATIO = $(r), METRIC = $(res) || TIME = $(time()-st)") end if Statistics.std(res_collector) >= 1e-2 # Detect if all solutions are similar to each other - Alp.get_option(m, :log_level) > 0 && + get_option(m, :log_level) > 0 && println("RATIO BRANCHING OFF due to solution variance test passed.") - Alp.set_option(m, :partition_scaling_factor_branch, false) # If an incumbent ratio is selected, then stop the branching scheme + set_option(m, :partition_scaling_factor_branch, false) # If an incumbent ratio is selected, then stop the branching scheme end if !isempty(m.best_sol) - m.discretization = Alp.add_adaptive_partition( + m.discretization = add_adaptive_partition( m, use_disc = m.discretization, branching = true, @@ -500,7 +497,7 @@ function update_partition_scaling_factor(m::Optimizer, presolve = false) use_solution = m.best_sol, ) else - m.discretization = Alp.add_adaptive_partition( + m.discretization = add_adaptive_partition( m, use_disc = m.discretization, branching = true, @@ -508,7 +505,7 @@ function update_partition_scaling_factor(m::Optimizer, presolve = false) ) end - Alp.get_option(m, :log_level) > 0 && println("INCUMB_RATIO = $(incumb_ratio)") + get_option(m, :log_level) > 0 && println("INCUMB_RATIO = $(incumb_ratio)") return incumb_ratio end @@ -516,13 +513,13 @@ end function disc_branch_solve(m::Optimizer) # ================= Solve Start ================ # - Alp.set_mip_time_limit(m) + set_mip_time_limit(m) start_bounding_solve = time() JuMP.optimize!(m.model_mip) status = MOI.get(m.model_mip, MOI.TerminationStatus()) cputime_branch_bounding_solve = time() - start_bounding_solve m.logs[:total_time] += cputime_branch_bounding_solve - m.logs[:time_left] = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + m.logs[:time_left] = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) # ================= Solve End ================ # if status in STATUS_OPT || status in STATUS_LIMIT @@ -532,9 +529,9 @@ function disc_branch_solve(m::Optimizer) end # Safety scheme - if Alp.is_min_sense(m) + if is_min_sense(m) return -Inf - elseif Alp.is_max_sense(m) + elseif is_max_sense(m) return Inf end end diff --git a/src/embedding.jl b/src/embedding.jl index 2b26c2e5..5a83a0cb 100644 --- a/src/embedding.jl +++ b/src/embedding.jl @@ -5,14 +5,14 @@ function embedding_map(λCnt::Int, encoding::Any = ebd_gray, ibs::Bool = false) map = Dict() - encoding = Alp.resolve_encoding_key(encoding) + encoding = resolve_encoding_key(encoding) L = Int(ceil(log(2, λCnt - 1))) for i in 1:L*2 map[i] = Set() end H = [encoding(i, L) for i in 0:max(1, (2^L - 1))] map[:H_orig] = H - map[:H] = [Alp.ebd_support_binary_vec(H[i]) for i in 1:length(H)] # Store the map + map[:H] = [ebd_support_binary_vec(H[i]) for i in 1:length(H)] # Store the map map[:L] = L !is_compatible_encoding(H) && error("Encodign method is not SOS-2 compatible...") @@ -57,7 +57,7 @@ end This function is the same σ() function described in Vielma and Nemhauser 2011. """ function ebd_σ(b::String) - sv = Alp.ebd_support_bool_vec(b) + sv = ebd_support_bool_vec(b) return [i for i in 1:length(sv) if sv[i]] end @@ -82,10 +82,7 @@ end function is_compatible_encoding(code_seq::Vector) for i in 1:(length(code_seq)-1) sum( - abs.( - Alp.ebd_support_bool_vec(code_seq[i]) - - Alp.ebd_support_bool_vec(code_seq[i+1]) - ), + abs.(ebd_support_bool_vec(code_seq[i]) - ebd_support_bool_vec(code_seq[i+1])), ) != 1 && return false end return true @@ -132,8 +129,8 @@ function ebd_link_xα( # Expression expansion for i in 1:P - code_vec = Alp.ebd_support_bool_vec(code_seq[i]) - lifters, exprs = Alp.ebd_link_expression(code_vec, lifters, exprs, i) + code_vec = ebd_support_bool_vec(code_seq[i]) + lifters, exprs = ebd_link_expression(code_vec, lifters, exprs, i) end # Construct Variable Vector @@ -145,11 +142,7 @@ function ebd_link_xα( base_name = "αA$(var_idx)" ) for i in keys(lifters) # Build first-level evaluation - Alp.relaxation_multilinear_binary( - m.model_mip, - α_A[lifters[i]-L], - [α[j] for j in i], - ) + relaxation_multilinear_binary(m.model_mip, α_A[lifters[i]-L], [α[j] for j in i]) end α_R = [α; α_A] # Initialize/re-arrgange the variable sequence diff --git a/src/heuristics.jl b/src/heuristics.jl index 1f37a32e..402a788d 100644 --- a/src/heuristics.jl +++ b/src/heuristics.jl @@ -21,9 +21,9 @@ function update_disc_cont_var(m::Optimizer) (var_diffs = _eval_var_diffs!(m.nonconvex_terms, var_diffs)) distance = Dict(zip(var_idxs, var_diffs)) - Alp.weighted_min_vertex_cover(m, distance) + weighted_min_vertex_cover(m, distance) - (Alp.get_option(m, :log_level) > 100) && + (get_option(m, :log_level) > 100) && println("updated partition var selection => $(m.disc_vars)") return end diff --git a/src/log.jl b/src/log.jl index bf020320..4e456441 100644 --- a/src/log.jl +++ b/src/log.jl @@ -5,7 +5,7 @@ function create_logs!(m) # Timers logs[:presolve_time] = 0.0 # Total presolve-time of the algorithm logs[:total_time] = 0.0 # Total run-time of the algorithm - logs[:time_left] = Alp.get_option(m, :time_limit) # Total remaining time of the algorithm if time-out is specified + logs[:time_left] = get_option(m, :time_limit) # Total remaining time of the algorithm if time-out is specified # Values logs[:obj] = [] # Iteration-based objective @@ -23,15 +23,15 @@ end function reset_timer(m::Optimizer) m.logs[:total_time] = 0.0 - m.logs[:time_left] = Alp.get_option(m, :time_limit) + m.logs[:time_left] = get_option(m, :time_limit) return m end function logging_summary(m::Optimizer) - if Alp.get_option(m, :log_level) > 0 + if get_option(m, :log_level) > 0 printstyled("\nPROBLEM STATISTICS\n", color = :cyan, bold = true) - Alp.is_min_sense(m) && (println(" Objective sense = Min")) - Alp.is_max_sense(m) && (println(" Objective sense = Max")) + is_min_sense(m) && (println(" Objective sense = Min")) + is_max_sense(m) && (println(" Objective sense = Max")) println( " # Variables = ", length([i for i in 1:m.num_var_orig if m.var_type[i] == :Cont]) + @@ -45,7 +45,7 @@ function logging_summary(m::Optimizer) println(" # Constraints = ", m.num_constr_orig) println(" # NL Constraints = ", m.num_nlconstr_orig) println(" # Linear Constraints = ", m.num_lconstr_orig) - Alp.get_option(m, :recognize_convex) && println( + get_option(m, :recognize_convex) && println( " # Detected convex constraints = $(length([i for i in m.constr_structure if i == :convex]))", ) println(" # Detected nonlinear terms = ", length(m.nonconvex_terms)) @@ -56,7 +56,7 @@ function logging_summary(m::Optimizer) println(" # Potential variables for partitioning = ", length(m.disc_vars)) printstyled("SUB-SOLVERS USED BY ALPINE\n", color = :cyan, bold = true) - if Alp.get_option(m, :minlp_solver) === nothing + if get_option(m, :minlp_solver) === nothing println(" NLP local solver = ", m.nlp_solver_id) else println(" MINLP local solver = ", m.minlp_solver_id) @@ -69,66 +69,59 @@ function logging_summary(m::Optimizer) Pkg.TOML.parse(read(string(pkgdir(Alpine), "/Project.toml"), String))["version"], ) - if Alp.is_min_sense(m) + if is_min_sense(m) println( " Maximum iterations (lower-bounding MIPs) = ", - Alp.get_option(m, :max_iter), + get_option(m, :max_iter), ) - elseif Alp.is_max_sense(m) + elseif is_max_sense(m) println( " Maximum iterations (upper-bounding MIPs) = ", - Alp.get_option(m, :max_iter), + get_option(m, :max_iter), ) else - println( - " Maximum iterations (bounding MIPs) = ", - Alp.get_option(m, :max_iter), - ) + println(" Maximum iterations (bounding MIPs) = ", get_option(m, :max_iter)) end - println( - " Relative global optimality gap = ", - Alp.get_option(m, :rel_gap) * 100, - "%", - ) + println(" Relative global optimality gap = ", get_option(m, :rel_gap) * 100, "%") - if Alp.get_option(m, :disc_var_pick) == 0 + if get_option(m, :disc_var_pick) == 0 println(" Potential variables chosen for partitioning = All") - elseif Alp.get_option(m, :disc_var_pick) == 1 + elseif get_option(m, :disc_var_pick) == 1 println( " Potential variables chosen for partitioning = Minimum vertex cover", ) end - if Alp.get_option(m, :partition_scaling_factor_branch) + if get_option(m, :partition_scaling_factor_branch) println(" Partition scaling factor branch activated") else println( " Partition scaling factor = ", - Alp.get_option(m, :partition_scaling_factor), + get_option(m, :partition_scaling_factor), ) end - (Alp.get_option(m, :convhull_ebd)) && println(" Using convhull_ebd formulation") - (Alp.get_option(m, :convhull_ebd)) && - println(" Encoding method = $(Alp.get_option(m, :convhull_ebd_encode))") - (Alp.get_option(m, :convhull_ebd)) && println( - " Independent branching scheme = $(Alp.get_option(m, :convhull_ebd_ibs))", + (get_option(m, :convhull_ebd)) && println(" Using convhull_ebd formulation") + (get_option(m, :convhull_ebd)) && + println(" Encoding method = $(get_option(m, :convhull_ebd_encode))") + (get_option(m, :convhull_ebd)) && println( + " Independent branching scheme = $(get_option(m, :convhull_ebd_ibs))", ) - println(" Bound-tightening presolve = ", Alp.get_option(m, :presolve_bt)) - Alp.get_option(m, :presolve_bt) && println( + println(" Bound-tightening presolve = ", get_option(m, :presolve_bt)) + get_option(m, :presolve_bt) && println( " Maximum iterations (OBBT) = ", - Alp.get_option(m, :presolve_bt_max_iter), + get_option(m, :presolve_bt_max_iter), ) end end function logging_head(m::Optimizer) - if Alp.is_min_sense(m) + if is_min_sense(m) printstyled("LOWER-BOUNDING ITERATIONS", color = :cyan, bold = true) UB_iter = "Incumbent" UB = "Best Incumbent" LB = "Lower Bound" - elseif Alp.is_max_sense(m) + elseif is_max_sense(m) printstyled("UPPER-BOUNDING ITERATIONS", color = :cyan, bold = true) UB_iter = "Incumbent" UB = "Best Incumbent" @@ -246,7 +239,7 @@ function summary_status(m::Optimizer) if m.detected_bound && m.detected_incumbent m.alpine_status = - m.best_rel_gap > Alp.get_option(m, :rel_gap) ? MOI.OTHER_LIMIT : MOI.OPTIMAL + m.best_rel_gap > get_option(m, :rel_gap) ? MOI.OTHER_LIMIT : MOI.OPTIMAL elseif m.status[:bounding_solve] == MOI.INFEASIBLE m.alpine_status = MOI.INFEASIBLE elseif m.detected_bound && !m.detected_incumbent @@ -269,11 +262,11 @@ end # cnt = length([1 for j in keys(m.nonconvex_terms) if m.nonconvex_terms[j][:nonlinear_type] == i]) # cnt > 0 && println("\tTerm $(i) Count = $(cnt) ") # end -# println(" Maximum solution time = ", Alp.get_option(m, :time_limit)) -# println(" Basic bound propagation = ", Alp.get_option(m, :presolve_bp)) -# println(" Conseuctive solution rejection = after ", Alp.get_option(m, :disc_consecutive_forbid), " times") -# Alp.get_option(m, :presolve_bt) && println("bound tightening presolve algorithm = ", Alp.get_option(m, :presolve_bt)_algo) -# Alp.get_option(m, :presolve_bt) && println("bound tightening presolve width tolerance = ", Alp.get_option(m, :presolve_bt)_width_tol) -# Alp.get_option(m, :presolve_bt) && println("bound tightening presolve output tolerance = ", Alp.get_option(m, :presolve_bt)_output_tol) -# Alp.get_option(m, :presolve_bt) && println("bound tightening presolve relaxation = ", Alp.get_option(m, :presolve_bt)_relax) -# Alp.get_option(m, :presolve_bt) && println("bound tightening presolve mip regulation time = ", Alp.get_option(m, :presolve_bt)_mip_time_limit) +# println(" Maximum solution time = ", get_option(m, :time_limit)) +# println(" Basic bound propagation = ", get_option(m, :presolve_bp)) +# println(" Conseuctive solution rejection = after ", get_option(m, :disc_consecutive_forbid), " times") +# get_option(m, :presolve_bt) && println("bound tightening presolve algorithm = ", get_option(m, :presolve_bt)_algo) +# get_option(m, :presolve_bt) && println("bound tightening presolve width tolerance = ", get_option(m, :presolve_bt)_width_tol) +# get_option(m, :presolve_bt) && println("bound tightening presolve output tolerance = ", get_option(m, :presolve_bt)_output_tol) +# get_option(m, :presolve_bt) && println("bound tightening presolve relaxation = ", get_option(m, :presolve_bt)_relax) +# get_option(m, :presolve_bt) && println("bound tightening presolve mip regulation time = ", get_option(m, :presolve_bt)_mip_time_limit) diff --git a/src/main_algorithm.jl b/src/main_algorithm.jl index a5da7b67..b6b58a60 100644 --- a/src/main_algorithm.jl +++ b/src/main_algorithm.jl @@ -28,12 +28,12 @@ end function load!(m::Optimizer) # Initialize NLP interface - requested_features = Alp.features_available(m) + requested_features = features_available(m) if m.d_orig !== nothing MOI.initialize(m.d_orig, requested_features::Vector{Symbol}) end for feat in requested_features - if !(feat in Alp.features_available(m)) + if !(feat in features_available(m)) error("Unsupported feature $feat") end end @@ -41,11 +41,11 @@ function load!(m::Optimizer) # Collect objective & constraint expressions if m.has_nl_objective m.obj_expr_orig = - Alp.expr_isolate_const(_variable_index_to_index(MOI.objective_expr(m.d_orig))) # see in nlexpr.jl if this expr isolation has any issue + expr_isolate_const(_variable_index_to_index(MOI.objective_expr(m.d_orig))) # see in nlexpr.jl if this expr isolation has any issue elseif m.objective_function isa Nothing m.obj_expr_orig = Expr(:call, :+) else - m.obj_expr_orig = Alp._moi_function_to_expr(m.objective_function) + m.obj_expr_orig = _moi_function_to_expr(m.objective_function) end # Collect original variable type and build dynamic variable type space @@ -54,7 +54,7 @@ function load!(m::Optimizer) m.bin_vars = [i for i in 1:m.num_var_orig if m.var_type[i] == :Bin] if !isempty(m.int_vars) || !isempty(m.bin_vars) - (Alp.get_option(m, :minlp_solver) === nothing) && (error( + (get_option(m, :minlp_solver) === nothing) && (error( "No MINLP local solver specified; use option 'minlp_solver' to specify a MINLP local solver", )) end @@ -95,51 +95,51 @@ function load!(m::Optimizer) isa(m.obj_expr_orig, Number) && (m.obj_structure = :constant) # Populate data to create the bounding MIP model - Alp.recategorize_var(m) # Initial round of variable re-categorization + recategorize_var(m) # Initial round of variable re-categorization :Int in m.var_type_orig && error( "Alpine does not support MINLPs with generic integer (non-binary) variables yet!", ) # Solver-dependent detection - Alp._fetch_mip_solver_identifier(m) - Alp._fetch_nlp_solver_identifier(m) - Alp._fetch_minlp_solver_identifier(m) + _fetch_mip_solver_identifier(m) + _fetch_nlp_solver_identifier(m) + _fetch_minlp_solver_identifier(m) # Main Algorithmic Initialization - Alp.process_expr(m) # Compact process of every expression - Alp.init_tight_bound(m) # Initialize bounds for algorithmic processes - Alp.resolve_var_bounds(m) # resolve lifted var bounds - Alp.pick_disc_vars(m) # Picking variables to be discretized - Alp.init_disc(m) # Initialize discretization dictionaries + process_expr(m) # Compact process of every expression + init_tight_bound(m) # Initialize bounds for algorithmic processes + resolve_var_bounds(m) # resolve lifted var bounds + pick_disc_vars(m) # Picking variables to be discretized + init_disc(m) # Initialize discretization dictionaries # Turn-on bt presolve if variables are not discrete if isempty(m.int_vars) && length(m.bin_vars) <= 50 && m.num_var_orig <= 10000 && length(m.candidate_disc_vars) <= 300 && - Alp.get_option(m, :presolve_bt) == nothing - Alp.set_option(m, :presolve_bt, true) + get_option(m, :presolve_bt) == nothing + set_option(m, :presolve_bt, true) println("Automatically turning on bound-tightening presolve") - elseif Alp.get_option(m, :presolve_bt) == nothing # If no use indication - Alp.set_option(m, :presolve_bt, false) + elseif get_option(m, :presolve_bt) == nothing # If no use indication + set_option(m, :presolve_bt, false) end if length(m.bin_vars) > 200 || m.num_var_orig > 2000 println( "Automatically turning OFF 'partition_scaling_factor_branch' due to the size of the problem", ) - Alp.set_option(m, :partition_scaling_factor_branch, false) + set_option(m, :partition_scaling_factor_branch, false) end # Initialize the solution pool - m.bound_sol_pool = Alp.initialize_solution_pool(m, 0) # Initialize the solution pool + m.bound_sol_pool = initialize_solution_pool(m, 0) # Initialize the solution pool # Check if any illegal term exist in the warm-solution any(isnan, m.best_sol) && (m.best_sol = zeros(length(m.best_sol))) # Initialize log - Alp.logging_summary(m) + logging_summary(m) return end @@ -148,23 +148,23 @@ end High-level Function """ function MOI.optimize!(m::Optimizer) - Alp.load!(m) + load!(m) if getproperty(m, :presolve_infeasible) - Alp.summary_status(m) + summary_status(m) return end - Alp.presolve(m) + presolve(m) - if !Alp.check_exit(m) && Alp.get_option(m, :apply_partitioning) - Alp.global_solve(m) - Alp.get_option(m, :log_level) > 0 && Alp.logging_row_entry(m, finish_entry = true) + if !check_exit(m) && get_option(m, :apply_partitioning) + global_solve(m) + get_option(m, :log_level) > 0 && logging_row_entry(m, finish_entry = true) println( "====================================================================================================", ) end - Alp.summary_status(m) + summary_status(m) return end @@ -179,19 +179,19 @@ Each [`local_solve`](@ref) provides an incumbent feasible solution. The algorith """ function global_solve(m::Optimizer) - Alp.get_option(m, :log_level) > 0 && Alp.logging_head(m) - Alp.get_option(m, :presolve_track_time) || Alp.reset_timer(m) - while !Alp.check_exit(m) + get_option(m, :log_level) > 0 && logging_head(m) + get_option(m, :presolve_track_time) || reset_timer(m) + while !check_exit(m) m.logs[:n_iter] += 1 - Alp.create_bounding_mip(m) # Build the relaxation model - Alp.bounding_solve(m) # Solve the relaxation model - Alp.update_opt_gap(m) # Update optimality gap - Alp.check_exit(m) && break # Feasibility check - Alp.get_option(m, :log_level) > 0 && Alp.logging_row_entry(m) # Logging - Alp.local_solve(m) # Solve local model for feasible solution - Alp.update_opt_gap(m) # Update optimality gap - Alp.check_exit(m) && break # Detect optimality termination - Alp.algorithm_automation(m) # Automated adjustments + create_bounding_mip(m) # Build the relaxation model + bounding_solve(m) # Solve the relaxation model + update_opt_gap(m) # Update optimality gap + check_exit(m) && break # Feasibility check + get_option(m, :log_level) > 0 && logging_row_entry(m) # Logging + local_solve(m) # Solve local model for feasible solution + update_opt_gap(m) # Update optimality gap + check_exit(m) && break # Detect optimality termination + algorithm_automation(m) # Automated adjustments end return @@ -202,11 +202,11 @@ end """ function presolve(m::Optimizer) start_presolve = time() - Alp.get_option(m, :log_level) > 0 && + get_option(m, :log_level) > 0 && printstyled("PRESOLVE \n", color = :cyan, bold = true) - Alp.get_option(m, :log_level) > 0 && println(" Doing local search") + get_option(m, :log_level) > 0 && println(" Doing local search") - if Alp.get_option(m, :use_start_as_incumbent) + if get_option(m, :use_start_as_incumbent) # Check for bound feasibility of the warm start value if !(m.l_var_orig <= m.warm_start_value <= m.u_var_orig) error( @@ -222,44 +222,44 @@ function presolve(m::Optimizer) m.objective_function, ) end - Alp.update_incumbent(m, obj_warmval, m.warm_start_value) + update_incumbent(m, obj_warmval, m.warm_start_value) m.status[:local_solve] = MOI.OPTIMIZE_NOT_CALLED - Alp.get_option(m, :log_level) > 0 && println( + get_option(m, :log_level) > 0 && println( " Using warm starting point as a local incumbent solution with value $(round(m.best_obj, digits=4))", ) else - Alp.local_solve(m, presolve = true) + local_solve(m, presolve = true) end # Solver status if m.status[:local_solve] in union(STATUS_OPT, STATUS_LIMIT, STATUS_WARM_START) - if !Alp.get_option(m, :use_start_as_incumbent) - Alp.get_option(m, :log_level) > 0 && println( + if !get_option(m, :use_start_as_incumbent) + get_option(m, :log_level) > 0 && println( " Local solver returns a feasible point with value $(round(m.best_obj, digits=4))", ) end - Alp.bound_tightening_wrapper(m, use_bound = true) # performs bound-tightening with the local solve objective value - Alp.get_option(m, :presolve_bt) && Alp.init_disc(m) # Re-initialize discretization dictionary on tight bounds - Alp.get_option(m, :partition_scaling_factor_branch) && (Alp.set_option( + bound_tightening_wrapper(m, use_bound = true) # performs bound-tightening with the local solve objective value + get_option(m, :presolve_bt) && init_disc(m) # Re-initialize discretization dictionary on tight bounds + get_option(m, :partition_scaling_factor_branch) && (set_option( m, :partition_scaling_factor, - Alp.update_partition_scaling_factor(m, true), + update_partition_scaling_factor(m, true), )) elseif m.status[:local_solve] in STATUS_INF - (Alp.get_option(m, :log_level) > 0) && + (get_option(m, :log_level) > 0) && println(" Bound tightening without objective bounds (OBBT)") - Alp.bound_tightening_wrapper(m, use_bound = false) # do bound tightening without objective value - (Alp.get_option(m, :partition_scaling_factor_branch)) && (Alp.set_option( + bound_tightening_wrapper(m, use_bound = false) # do bound tightening without objective value + (get_option(m, :partition_scaling_factor_branch)) && (set_option( m, :partition_scaling_factor, - Alp.update_partition_scaling_factor(m, true), + update_partition_scaling_factor(m, true), )) - Alp.get_option(m, :presolve_bt) && Alp.init_disc(m) + get_option(m, :presolve_bt) && init_disc(m) elseif m.status[:local_solve] == MOI.INVALID_MODEL @warn " Warning: Presolve ends with local solver yielding $(m.status[:local_solve]). \n This may come from Ipopt's `:Not_Enough_Degrees_Of_Freedom`." - elseif !Alp.get_option(m, :use_start_as_incumbent) + elseif !get_option(m, :use_start_as_incumbent) @warn " Warning: Presolve ends with local solver yielding $(m.status[:local_solve])." end @@ -268,12 +268,12 @@ function presolve(m::Optimizer) m.logs[:total_time] = m.logs[:presolve_time] m.logs[:time_left] -= m.logs[:presolve_time] - if Alp.get_option(m, :presolve_bt) - (Alp.get_option(m, :log_level) > 0) && println( + if get_option(m, :presolve_bt) + (get_option(m, :log_level) > 0) && println( " Post-presolve optimality gap: $(round(m.presolve_best_rel_gap; digits = 3))%", ) end - (Alp.get_option(m, :log_level) > 0) && + (get_option(m, :log_level) > 0) && println(" Completed presolve in $(round.(m.logs[:total_time]; digits = 2))s") return @@ -285,16 +285,12 @@ end function algorithm_automation(m::Optimizer) # Only if a feasible solution is detected - if (Alp.get_option(m, :disc_var_pick) == 3) && m.detected_incumbent - Alp.update_disc_cont_var(m) + if (get_option(m, :disc_var_pick) == 3) && m.detected_incumbent + update_disc_cont_var(m) end - if Alp.get_option(m, :partition_scaling_factor_branch) - Alp.set_option( - m, - :partition_scaling_factor, - Alp.update_partition_scaling_factor(m, true), - ) # Only perform for a maximum three times + if get_option(m, :partition_scaling_factor_branch) + set_option(m, :partition_scaling_factor, update_partition_scaling_factor(m, true)) # Only perform for a maximum three times end return @@ -306,7 +302,7 @@ end function check_exit(m::Optimizer) # constant objective with feasible local solve check - if Alp.expr_isconst(m.obj_expr_orig) && + if expr_isconst(m.obj_expr_orig) && (m.status[:local_solve] == MOI.OPTIMAL || m.status == MOI.LOCALLY_SOLVED) # m.best_bound = eval(m.obj_expr_orig) m.best_bound = m.obj_expr_orig @@ -325,15 +321,15 @@ function check_exit(m::Optimizer) m.status[:bounding_solve] == MOI.DUAL_INFEASIBLE && return true # Optimality check - if m.best_rel_gap <= Alp.get_option(m, :rel_gap) + if m.best_rel_gap <= get_option(m, :rel_gap) m.detected_bound = true return true end - m.logs[:n_iter] >= Alp.get_option(m, :max_iter) && return true - m.best_abs_gap <= Alp.get_option(m, :abs_gap) && return true + m.logs[:n_iter] >= get_option(m, :max_iter) && return true + m.best_abs_gap <= get_option(m, :abs_gap) && return true # User-limits check - m.logs[:time_left] < Alp.get_option(m, :tol) && return true + m.logs[:time_left] < get_option(m, :tol) && return true return false end @@ -341,7 +337,7 @@ end function load_nonlinear_model(m::Optimizer, model::MOI.ModelLike, l_var, u_var) x = MOI.add_variables(model, m.num_var_orig) for i in eachindex(x) - set = Alp._bound_set(l_var[i], u_var[i]) + set = _bound_set(l_var[i], u_var[i]) if set !== nothing MOI.add_constraint(model, x[i], set) end @@ -380,7 +376,7 @@ function set_variable_type(model::MOI.ModelLike, xs, variable_types) end """ - Alp.local_solve(m::Optimizer, presolve::Bool=false) + local_solve(m::Optimizer, presolve::Bool=false) Perform a local NLP or MINLP solve to obtain a feasible solution. The `presolve` option is set to `true` when the function is invoked in [`presolve`](@ref). @@ -393,34 +389,28 @@ function local_solve(m::Optimizer; presolve = false) var_type_screener = [i for i in m.var_type_orig if i in [:Bin, :Int]] if presolve - if !isempty(var_type_screener) && Alp.get_option(m, :minlp_solver) !== nothing - local_solve_model = MOI.instantiate( - Alp.get_option(m, :minlp_solver), - with_bridge_type = Float64, - ) + if !isempty(var_type_screener) && get_option(m, :minlp_solver) !== nothing + local_solve_model = + MOI.instantiate(get_option(m, :minlp_solver), with_bridge_type = Float64) elseif !isempty(var_type_screener) - local_solve_model = MOI.instantiate( - Alp.get_option(m, :nlp_solver), - with_bridge_type = Float64, - ) + local_solve_model = + MOI.instantiate(get_option(m, :nlp_solver), with_bridge_type = Float64) else - local_solve_model = MOI.instantiate( - Alp.get_option(m, :nlp_solver), - with_bridge_type = Float64, - ) + local_solve_model = + MOI.instantiate(get_option(m, :nlp_solver), with_bridge_type = Float64) end else local_solve_model = - MOI.instantiate(Alp.get_option(m, :nlp_solver), with_bridge_type = Float64) + MOI.instantiate(get_option(m, :nlp_solver), with_bridge_type = Float64) end if presolve == false - l_var, u_var = Alp.fix_domains(m) + l_var, u_var = fix_domains(m) else l_var, u_var = m.l_var_tight[1:m.num_var_orig], m.u_var_tight[1:m.num_var_orig] end - x = Alp.load_nonlinear_model(m, local_solve_model, l_var, u_var) + x = load_nonlinear_model(m, local_solve_model, l_var, u_var) if m.d_orig !== nothing MOI.initialize(m.d_orig, [:Grad, :Jac, :Hess, :HessVec, :ExprGraph]) # Safety scheme for sub-solvers re-initializing the NLPEvaluator end @@ -436,11 +426,11 @@ function local_solve(m::Optimizer; presolve = false) # The only case when MINLP solver is actually used if presolve && !isempty(var_type_screener) - if Alp.get_option(m, :minlp_solver) === nothing + if get_option(m, :minlp_solver) === nothing error("Provide a valid MINLP solver") # do_heuristic = true else - Alp.set_variable_type(local_solve_model, x, m.var_type_orig) + set_variable_type(local_solve_model, x, m.var_type_orig) end end @@ -454,7 +444,7 @@ function local_solve(m::Optimizer; presolve = false) cputime_local_solve = time() - start_local_solve m.logs[:total_time] += cputime_local_solve - m.logs[:time_left] = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + m.logs[:time_left] = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) # if do_heuristic # m.status[:local_solve] = heu_basic_rounding(m, MOI.get(local_solve_model, MOI.VariablePrimal(), x)) @@ -480,12 +470,12 @@ function local_solve(m::Optimizer; presolve = false) end @assert length(candidate_sol) == length(sol_temp) - Alp.update_incumbent(m, candidate_obj, candidate_sol) + update_incumbent(m, candidate_obj, candidate_sol) m.status[:local_solve] = local_nlp_status return # elseif local_nlp_status in STATUS_INF - # Alp.heu_pool_multistart(m) == MOI.LOCALLY_SOLVED && return + # heu_pool_multistart(m) == MOI.LOCALLY_SOLVED && return # push!(m.logs[:obj], "INF") # m.status[:local_solve] = MOI.LOCALLY_INFEASIBLE # return @@ -515,7 +505,7 @@ function local_solve(m::Optimizer; presolve = false) end """ - Alp.bounding_solve(m::Optimizer; kwargs...) + bounding_solve(m::Optimizer; kwargs...) This step usually solves a convex MILP/MIQCP/MIQCQP problem for lower bounding the given minimization problem. It solves the problem built upon a piecewise convexification based on the discretization sictionary of some variables. @@ -525,7 +515,7 @@ function bounding_solve(m::Optimizer) convertor = Dict(MOI.MAX_SENSE => :<, MOI.MIN_SENSE => :>) # Updates time metric and the termination bounds - Alp.set_mip_time_limit(m) + set_mip_time_limit(m) # update_boundstop_options(m) # ================= Solve Start ================ # @@ -533,7 +523,7 @@ function bounding_solve(m::Optimizer) JuMP.optimize!(m.model_mip) status = JuMP.termination_status(m.model_mip) m.logs[:total_time] += time() - start_bounding_solve - m.logs[:time_left] = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + m.logs[:time_left] = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) # ================= Solve End ================ # if status in STATUS_OPT || status in STATUS_LIMIT @@ -546,11 +536,11 @@ function bounding_solve(m::Optimizer) ] # Experimental code - Alp.measure_relaxed_deviation(m, sol = candidate_bound_sol) - if Alp.get_option(m, :disc_consecutive_forbid) > 0 + measure_relaxed_deviation(m, sol = candidate_bound_sol) + if get_option(m, :disc_consecutive_forbid) > 0 m.bound_sol_history[mod( m.logs[:n_iter] - 1, - Alp.get_option(m, :disc_consecutive_forbid), + get_option(m, :disc_consecutive_forbid), )+1] = copy(candidate_bound_sol) # Requires proper offseting end push!(m.logs[:bound], candidate_bound) @@ -590,17 +580,17 @@ This function helps pick the variables for discretization. The method chosen dep In case when `indices::Int` is provided, the method is chosen as built-in method. Currently, there are two built-in options for users as follows: -* `max_cover (Alp.get_option(m, :disc_var_pick)=0, default)`: pick all variables involved in the non-linear term for discretization -* `min_vertex_cover (Alp.get_option(m, :disc_var_pick)=1)`: pick a minimum vertex cover for variables involved in non-linear terms so that each non-linear term is at least convexified +* `max_cover (get_option(m, :disc_var_pick)=0, default)`: pick all variables involved in the non-linear term for discretization +* `min_vertex_cover (get_option(m, :disc_var_pick)=1)`: pick a minimum vertex cover for variables involved in non-linear terms so that each non-linear term is at least convexified -For advanced usage, `Alp.get_option(m, :disc_var_pick)` allows `::Function` inputs. User can provide his/her own function to choose the variables for discretization. +For advanced usage, `get_option(m, :disc_var_pick)` allows `::Function` inputs. User can provide his/her own function to choose the variables for discretization. """ function pick_disc_vars(m::Optimizer) - disc_var_pick = Alp.get_option(m, :disc_var_pick) + disc_var_pick = get_option(m, :disc_var_pick) if isa(disc_var_pick, Function) - # eval(Alp.get_option(m, :disc_var_pick))(m) + # eval(get_option(m, :disc_var_pick))(m) disc_var_pick(m) length(m.disc_vars) == 0 && length(m.nonconvex_terms) > 0 && @@ -609,15 +599,15 @@ function pick_disc_vars(m::Optimizer) ) elseif isa(disc_var_pick, Int) if disc_var_pick == 0 - Alp.get_candidate_disc_vars(m) + get_candidate_disc_vars(m) elseif disc_var_pick == 1 - Alp.min_vertex_cover(m) + min_vertex_cover(m) elseif disc_var_pick == 2 - (length(m.candidate_disc_vars) > 15) ? Alp.min_vertex_cover(m) : - Alp.get_candidate_disc_vars(m) + (length(m.candidate_disc_vars) > 15) ? min_vertex_cover(m) : + get_candidate_disc_vars(m) elseif disc_var_pick == 3 # Initial - (length(m.candidate_disc_vars) > 15) ? Alp.min_vertex_cover(m) : - Alp.get_candidate_disc_vars(m) + (length(m.candidate_disc_vars) > 15) ? min_vertex_cover(m) : + get_candidate_disc_vars(m) else error( "Unsupported default indicator for picking variables for discretization", diff --git a/src/multilinear.jl b/src/multilinear.jl index d7539383..8217a04e 100644 --- a/src/multilinear.jl +++ b/src/multilinear.jl @@ -13,26 +13,26 @@ function amp_post_convhull(m::Optimizer; kwargs...) nl_type = m.nonconvex_terms[k][:nonlinear_type] if ((nl_type == :MULTILINEAR) || (nl_type == :BILINEAR)) && (m.nonconvex_terms[k][:convexified] == false) - λ, α = Alp.amp_convexify_multilinear(m, k, λ, α, d) + λ, α = amp_convexify_multilinear(m, k, λ, α, d) contains_multilinear = true elseif nl_type == :MONOMIAL && !m.nonconvex_terms[k][:convexified] - λ, α = Alp.amp_convexify_quadratic_univariate(m, k, λ, α, d) + λ, α = amp_convexify_quadratic_univariate(m, k, λ, α, d) elseif nl_type == :BINLIN && !m.nonconvex_terms[k][:convexified] - β = Alp.amp_convexify_binlin(m, k, β) + β = amp_convexify_binlin(m, k, β) elseif nl_type == :BINPROD && !m.nonconvex_terms[k][:convexified] - β = Alp.amp_convexify_multilinear_binary(m, k, β) + β = amp_convexify_multilinear_binary(m, k, β) end end # Add lambda linking constraints if m.options.linking_constraints && contains_multilinear - Alp._add_multilinear_linking_constraints(m, λ) + _add_multilinear_linking_constraints(m, λ) end # Code for warm starting bounding MIP iterations - Alp.get_option(m, :convhull_warmstart) && + get_option(m, :convhull_warmstart) && !isempty(m.best_bound_sol) && - Alp.amp_warmstart_α(m, α) + amp_warmstart_α(m, α) return end @@ -46,9 +46,9 @@ function amp_convexify_multilinear( ) m.nonconvex_terms[k][:convexified] = true # Bookeeping the convexified terms - ml_indices, dim, extreme_point_cnt = Alp.amp_convhull_prepare(m, discretization, k) # convert key to easy read mode - λ = Alp.amp_convhull_λ(m, k, ml_indices, λ, extreme_point_cnt, dim) - λ = Alp.populate_convhull_extreme_values( + ml_indices, dim, extreme_point_cnt = amp_convhull_prepare(m, discretization, k) # convert key to easy read mode + λ = amp_convhull_λ(m, k, ml_indices, λ, extreme_point_cnt, dim) + λ = populate_convhull_extreme_values( m, discretization, ml_indices, @@ -56,16 +56,8 @@ function amp_convexify_multilinear( dim, ones(Int, length(dim)), ) - α = Alp.amp_convhull_α(m, ml_indices, α, dim, discretization) - Alp.amp_post_convhull_constrs( - m, - λ, - α, - ml_indices, - dim, - extreme_point_cnt, - discretization, - ) + α = amp_convhull_α(m, ml_indices, α, dim, discretization) + amp_post_convhull_constrs(m, λ, α, ml_indices, dim, extreme_point_cnt, discretization) return λ, α end @@ -80,11 +72,11 @@ function amp_convexify_quadratic_univariate( m.nonconvex_terms[k][:convexified] = true # Bookeeping the convexified terms monomial_index, dim, extreme_point_cnt = - Alp.amp_convhull_prepare(m, discretization, k, monomial = true) - λ = Alp.amp_convhull_λ(m, k, monomial_index, λ, extreme_point_cnt, dim) - λ = Alp.populate_convhull_extreme_values(m, discretization, monomial_index, λ, 2) - α = Alp.amp_convhull_α(m, [monomial_index], α, dim, discretization) - Alp.amp_post_convhull_constrs(m, λ, α, monomial_index, discretization) + amp_convhull_prepare(m, discretization, k, monomial = true) + λ = amp_convhull_λ(m, k, monomial_index, λ, extreme_point_cnt, dim) + λ = populate_convhull_extreme_values(m, discretization, monomial_index, λ, 2) + α = amp_convhull_α(m, [monomial_index], α, dim, discretization) + amp_post_convhull_constrs(m, λ, α, monomial_index, discretization) return λ, α end @@ -110,7 +102,7 @@ function amp_convexify_binlin(m::Optimizer, k::Any, β::Dict) bin_idx = bin_idx[1] cont_idx = cont_idx[1] - Alp.relaxation_bilinear( + relaxation_bilinear( m.model_mip, _index_to_variable_ref(m.model_mip, lift_idx), _index_to_variable_ref(m.model_mip, bin_idx), @@ -138,7 +130,7 @@ function amp_convexify_multilinear_binary(m::Optimizer, k::Any, β::Dict) x = [_index_to_variable_ref(m.model_mip, i) for i in m.nonconvex_terms[k][:var_idxs]] if length(x) >= 1 - Alp.relaxation_multilinear_binary(m.model_mip, z, x) + relaxation_multilinear_binary(m.model_mip, z, x) end return β @@ -255,7 +247,7 @@ function populate_convhull_extreme_values( else for i in 1:dim[level] locator[level] = i - λ = Alp.populate_convhull_extreme_values( + λ = populate_convhull_extreme_values( m, discretization, indices, @@ -284,7 +276,7 @@ function amp_convhull_α( if !(i in keys(α)) lambda_cnt = length(discretization[i]) partition_cnt = length(discretization[i]) - 1 - if Alp.get_option(m, :convhull_ebd) && partition_cnt > 2 + if get_option(m, :convhull_ebd) && partition_cnt > 2 αCnt = Int(ceil(log(2, partition_cnt))) α[i] = JuMP.@variable( m.model_mip, @@ -326,7 +318,7 @@ function amp_warmstart_α(m::Optimizer, α::Dict) if m.bound_sol_pool[:cnt] >= 2 # can only warm-start the problem when pool is large enough ws_idx = -1 - Alp.is_min_sense(m) ? ws_obj = Inf : ws_obj = -Inf + is_min_sense(m) ? ws_obj = Inf : ws_obj = -Inf comp_opr = Dict(MOI.MIN_SENSE => :<, MOI.MAX_SENSE => :>) # Search for the pool for incumbent warm starter @@ -344,14 +336,14 @@ function amp_warmstart_α(m::Optimizer, α::Dict) for v in m.bound_sol_pool[:vars] partition_cnt = length(d[v]) - 1 active_j = - Alp.get_active_partition_idx(d, m.bound_sol_pool[:sol][ws_idx][v], v) + get_active_partition_idx(d, m.bound_sol_pool[:sol][ws_idx][v], v) for j in 1:partition_cnt j == active_j ? JuMP.set_start_value(α[v][j], 1.0) : JuMP.set_start_value(α[v][j], 0.0) end end m.bound_sol_pool[:stat][ws_idx] = :Warmstarter - Alp.get_option(m, :log_level) > 0 && println( + get_option(m, :log_level) > 0 && println( "!! WARM START bounding MIP using POOL SOL $(ws_idx) OBJ=$(m.bound_sol_pool[:obj][ws_idx])", ) end @@ -385,12 +377,12 @@ function amp_post_convhull_constrs( for (cnt, i) in enumerate(indices) l_cnt = length(d[i]) if m.var_type[i] == :Cont - Alp.amp_post_inequalities_cont(m, d, λ, α, indices, dim, i, cnt) # Add links between λ and α + amp_post_inequalities_cont(m, d, λ, α, indices, dim, i, cnt) # Add links between λ and α else error("EXCEPTION: unexpected variable type during integer related realxation") end sliced_indices = - [Alp.collect_indices(λ[indices][:indices], cnt, [k], dim) for k in 1:l_cnt] # Add x = f(λ) for convex representation of x value + [collect_indices(λ[indices][:indices], cnt, [k], dim) for k in 1:l_cnt] # Add x = f(λ) for convex representation of x value JuMP.@constraint( m.model_mip, _index_to_variable_ref(m.model_mip, i) == sum( @@ -432,11 +424,11 @@ function amp_post_convhull_constrs( ) # Add SOS-2 Constraints with basic encoding - if Alp.get_option(m, :convhull_ebd) && partition_cnt > 2 - ebd_map = Alp.embedding_map( + if get_option(m, :convhull_ebd) && partition_cnt > 2 + ebd_map = embedding_map( lambda_cnt, - Alp.get_option(m, :convhull_ebd_encode), - Alp.get_option(m, :convhull_ebd_ibs), + get_option(m, :convhull_ebd_encode), + get_option(m, :convhull_ebd_ibs), ) YCnt = Int(ebd_map[:L]) @assert YCnt == length(α[monomial_idx]) @@ -515,24 +507,20 @@ function amp_post_inequalities_cont( partition_cnt = lambda_cnt - 1 # Embedding formulation - if Alp.get_option(m, :convhull_formulation) == "sos2" && - Alp.get_option(m, :convhull_ebd) && + if get_option(m, :convhull_formulation) == "sos2" && + get_option(m, :convhull_ebd) && partition_cnt > 2 - ebd_map = Alp.embedding_map( + ebd_map = embedding_map( lambda_cnt, - Alp.get_option(m, :convhull_ebd_encode), - Alp.get_option(m, :convhull_ebd_ibs), + get_option(m, :convhull_ebd_encode), + get_option(m, :convhull_ebd_ibs), ) YCnt = Int(ebd_map[:L]) @assert YCnt == length(α[var_ind]) for i in 1:YCnt - p_sliced_indices = Alp.collect_indices( - λ[ml_indices][:indices], - cnt, - collect(ebd_map[i]), - dim, - ) - n_sliced_indices = Alp.collect_indices( + p_sliced_indices = + collect_indices(λ[ml_indices][:indices], cnt, collect(ebd_map[i]), dim) + n_sliced_indices = collect_indices( λ[ml_indices][:indices], cnt, collect(ebd_map[i+YCnt]), @@ -547,7 +535,7 @@ function amp_post_inequalities_cont( sum(λ[ml_indices][:vars][n_sliced_indices]) <= 1 - α[var_ind][i] ) end - Alp.get_option(m, :convhull_ebd_link) && Alp.ebd_link_xα( + get_option(m, :convhull_ebd_link) && ebd_link_xα( m, α[var_ind], lambda_cnt, @@ -559,9 +547,9 @@ function amp_post_inequalities_cont( end # SOS-2 Formulation - if Alp.get_option(m, :convhull_formulation) == "sos2" + if get_option(m, :convhull_formulation) == "sos2" for j in 1:lambda_cnt - sliced_indices = Alp.collect_indices(λ[ml_indices][:indices], cnt, [j], dim) + sliced_indices = collect_indices(λ[ml_indices][:indices], cnt, [j], dim) if (j == 1) JuMP.@constraint( m.model_mip, @@ -581,10 +569,9 @@ function amp_post_inequalities_cont( end end return - elseif Alp.get_option(m, :convhull_formulation) == "facet" + elseif get_option(m, :convhull_formulation) == "facet" for j in 1:(partition_cnt-1) # Constraint cluster of α >= f(λ) - sliced_indices = - Alp.collect_indices(λ[ml_indices][:indices], cnt, [1:j;], dim) + sliced_indices = collect_indices(λ[ml_indices][:indices], cnt, [1:j;], dim) JuMP.@constraint( m.model_mip, sum(α[var_ind][1:j]) >= sum(λ[ml_indices][:vars][sliced_indices]) @@ -592,7 +579,7 @@ function amp_post_inequalities_cont( end for j in 1:(partition_cnt-1) # Constraint cluster of α <= f(λ) sliced_indices = - Alp.collect_indices(λ[ml_indices][:indices], cnt, [1:(j+1);], dim) + collect_indices(λ[ml_indices][:indices], cnt, [1:(j+1);], dim) JuMP.@constraint( m.model_mip, sum(α[var_ind][1:j]) <= sum(λ[ml_indices][:vars][sliced_indices]) @@ -625,26 +612,26 @@ end _add_multilinear_linking_constraints(m::Optimizer, λ::Dict) This internal function adds linking constraints between λ multipliers corresponding to multilinear terms -that share more than two variables and are partitioned. For example, suppose we have λ[i], λ[j], and -λ[k] where i=(1,2,3), j=(1,2,4), and k=(1,2,5). λ[i] contains all multipliers for the extreme points +that share more than two variables and are partitioned. For example, suppose we have λ[i], λ[j], and +λ[k] where i=(1,2,3), j=(1,2,4), and k=(1,2,5). λ[i] contains all multipliers for the extreme points in the space of (x1,x2,x3). λ[j] contains all multipliers for the extreme points in the space of (x1,x2,x4). λ[k] contains all multipliers for the extreme points in the space of (x1,x2,x5). Using λ[i], λ[j], or λ[k], we can express multilinear function x1*x2. We define a linking variable μ(1,2) that represents the value of x1*x2. Linking constraints are - μ(1,2) == convex combination expr for x1*x2 using λ[i], - μ(1,2) == convex combination expr for x1*x2 using λ[j], and - μ(1,2) == convex combination expr for x1*x2 using λ[k]. + μ(1,2) == convex combination expr for x1*x2 using λ[i], + μ(1,2) == convex combination expr for x1*x2 using λ[j], and + μ(1,2) == convex combination expr for x1*x2 using λ[k]. Thus, these constraints link between λ[i], λ[j], and λ[k] variables. -Reference: J. Kim, J.P. Richard, M. Tawarmalani, Piecewise Polyhedral Relaxations of Multilinear Optimization, +Reference: J. Kim, J.P. Richard, M. Tawarmalani, Piecewise Polyhedral Relaxations of Multilinear Optimization, http://www.optimization-online.org/DB_HTML/2022/07/8974.html """ function _add_multilinear_linking_constraints(m::Optimizer, λ::Dict) if isnothing(m.linking_constraints_info) - m.linking_constraints_info = Alp._get_shared_multilinear_terms_info( + m.linking_constraints_info = _get_shared_multilinear_terms_info( λ, m.options.linking_constraints_degree_limit, ) @@ -684,8 +671,8 @@ end """ _get_shared_multilinear_terms_info(λ, linking_constraints_degree_limit) -This function checks to see if linking constraints are -necessary for a given vector of each multilinear terms and returns the approapriate +This function checks to see if linking constraints are +necessary for a given vector of each multilinear terms and returns the approapriate linking constraints information. """ function _get_shared_multilinear_terms_info( @@ -700,7 +687,7 @@ function _get_shared_multilinear_terms_info( return (linking_constraints_info = nothing) end - # Limit the linking constraints to a prescribed multilinear degree + # Limit the linking constraints to a prescribed multilinear degree if !isnothing(linking_constraints_degree_limit) && (linking_constraints_degree_limit < max_degree) max_degree = linking_constraints_degree_limit diff --git a/src/nlexpr.jl b/src/nlexpr.jl index 302c2bfa..9bff622f 100644 --- a/src/nlexpr.jl +++ b/src/nlexpr.jl @@ -4,11 +4,11 @@ process_expr(expr; kwargs...) High-level wrapper for processing expression with sub-tree operators """ function process_expr(m::Optimizer) - Alp.expr_initialization(m) # S0 : initialize the space for parsing and analyzing - Alp.expr_preprocess(m) # S1 : pre-process the negative sign in expressions - Alp.expr_parsing(m) # S2 : parsing the expressions for nonlinear information - Alp.expr_conversion(m) # S3 : convert lifted(linear) expressions into affine function - Alp.expr_finalized(m) # S4 : finalize process by extracting some measurements + expr_initialization(m) # S0 : initialize the space for parsing and analyzing + expr_preprocess(m) # S1 : pre-process the negative sign in expressions + expr_parsing(m) # S2 : parsing the expressions for nonlinear information + expr_conversion(m) # S3 : convert lifted(linear) expressions into affine function + expr_finalized(m) # S4 : finalize process by extracting some measurements return end @@ -34,13 +34,13 @@ end STEP 2: preprocess expression for trivial sub-trees and nasty pieces for easier later process """ function expr_preprocess(m::Optimizer) - Alp.expr_resolve_const(m.bounding_obj_expr_mip) - Alp.expr_resolve_sign(m.bounding_obj_expr_mip) - Alp.expr_flatten(m.bounding_obj_expr_mip) + expr_resolve_const(m.bounding_obj_expr_mip) + expr_resolve_sign(m.bounding_obj_expr_mip) + expr_flatten(m.bounding_obj_expr_mip) for i in 1:m.num_constr_orig - Alp.expr_resolve_const(m.bounding_constr_expr_mip[i]) - Alp.expr_resolve_sign(m.bounding_constr_expr_mip[i]) - Alp.expr_flatten(m.bounding_constr_expr_mip[i].args[2]) + expr_resolve_const(m.bounding_constr_expr_mip[i]) + expr_resolve_sign(m.bounding_constr_expr_mip[i]) + expr_flatten(m.bounding_constr_expr_mip[i].args[2]) end return @@ -52,30 +52,29 @@ STEP 3: parse expression for patterns on either the generic level or term level function expr_parsing(m::Optimizer) # Throw an error if the objective expression has fractional exponents - Alp.expr_is_fractional_exponent(m.bounding_obj_expr_mip) - is_structural = Alp.expr_constr_parsing(m.bounding_obj_expr_mip, m) + expr_is_fractional_exponent(m.bounding_obj_expr_mip) + is_structural = expr_constr_parsing(m.bounding_obj_expr_mip, m) if !is_structural - m.bounding_obj_expr_mip = Alp.expr_term_parsing(m.bounding_obj_expr_mip, 0, m) + m.bounding_obj_expr_mip = expr_term_parsing(m.bounding_obj_expr_mip, 0, m) m.obj_structure = :generic_linear end - (Alp.get_option(m, :log_level) > 199) && println("[OBJ] $(m.obj_expr_orig)") + (get_option(m, :log_level) > 199) && println("[OBJ] $(m.obj_expr_orig)") for i in 1:m.num_constr_orig expr = m.bounding_constr_expr_mip[i] # Throw an error if the constraint expression has fractional exponents if expr.args[1] in [:(<=), :(>=), :(==)] - Alp.expr_is_fractional_exponent(expr.args[2]) + expr_is_fractional_exponent(expr.args[2]) end - is_structural = Alp.expr_constr_parsing(expr, m, i) + is_structural = expr_constr_parsing(expr, m, i) if !is_structural - m.bounding_constr_expr_mip[i] = Alp.expr_term_parsing(expr, i, m) + m.bounding_constr_expr_mip[i] = expr_term_parsing(expr, i, m) m.constr_structure[i] = :generic_linear end - (Alp.get_option(m, :log_level) > 199) && - println("[CONSTR] $(m.constr_expr_orig[i])") + (get_option(m, :log_level) > 199) && println("[CONSTR] $(m.constr_expr_orig[i])") end return @@ -86,36 +85,35 @@ STEP 4: convert the parsed expressions into affine-based function that can be us """ function expr_conversion(m::Optimizer) if m.obj_structure == :generic_linear - m.bounding_obj_mip = Alp.expr_linear_to_affine(m.bounding_obj_expr_mip) + m.bounding_obj_mip = expr_linear_to_affine(m.bounding_obj_expr_mip) m.obj_structure = :affine end - Alp.get_option(m, :log_level) > 199 && println("type :: ", m.obj_structure) - Alp.get_option(m, :log_level) > 199 && println("lifted ::", m.bounding_obj_expr_mip) - Alp.get_option(m, :log_level) > 199 && - println("coeffs ::", m.bounding_obj_mip[:coefs]) - Alp.get_option(m, :log_level) > 199 && println("vars ::", m.bounding_obj_mip[:vars]) - Alp.get_option(m, :log_level) > 199 && println("sense ::", m.bounding_obj_mip[:sense]) - Alp.get_option(m, :log_level) > 199 && println("rhs ::", m.bounding_obj_mip[:rhs]) - Alp.get_option(m, :log_level) > 199 && println("----------------") + get_option(m, :log_level) > 199 && println("type :: ", m.obj_structure) + get_option(m, :log_level) > 199 && println("lifted ::", m.bounding_obj_expr_mip) + get_option(m, :log_level) > 199 && println("coeffs ::", m.bounding_obj_mip[:coefs]) + get_option(m, :log_level) > 199 && println("vars ::", m.bounding_obj_mip[:vars]) + get_option(m, :log_level) > 199 && println("sense ::", m.bounding_obj_mip[:sense]) + get_option(m, :log_level) > 199 && println("rhs ::", m.bounding_obj_mip[:rhs]) + get_option(m, :log_level) > 199 && println("----------------") for i in 1:m.num_constr_orig if m.constr_structure[i] == :generic_linear m.bounding_constr_mip[i] = - Alp.expr_linear_to_affine(m.bounding_constr_expr_mip[i]) + expr_linear_to_affine(m.bounding_constr_expr_mip[i]) m.constr_structure[i] = :affine end - Alp.get_option(m, :log_level) > 199 && println("type :: ", m.constr_structure[i]) - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("type :: ", m.constr_structure[i]) + get_option(m, :log_level) > 199 && println("lifted ::", m.bounding_constr_expr_mip[i]) - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("coeffs ::", m.bounding_constr_mip[i][:coefs]) - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("vars ::", m.bounding_constr_mip[i][:vars]) - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("sense ::", m.bounding_constr_mip[i][:sense]) - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("rhs ::", m.bounding_constr_mip[i][:rhs]) - Alp.get_option(m, :log_level) > 199 && println("----------------") + get_option(m, :log_level) > 199 && println("----------------") end return @@ -125,7 +123,7 @@ end STEP 5: collect information as needed for handy operations in the algorithm section """ function expr_finalized(m::Optimizer) - Alp.collect_nonconvex_vars(m) + collect_nonconvex_vars(m) m.candidate_disc_vars = sort(m.candidate_disc_vars) m.num_var_linear_mip = length(m.linear_terms) m.num_var_nonlinear_mip = length(m.nonconvex_terms) @@ -179,7 +177,7 @@ function expr_strip_const(expr, subs = [], rhs = 0.0) elseif expr.args[i].head == :ref continue elseif expr.args[i].head == :call - subs, rhs = Alp.expr_strip_const(expr.args[i], subs, rhs) + subs, rhs = expr_strip_const(expr.args[i], subs, rhs) end end @@ -205,22 +203,22 @@ function build_constr_block(y_idx::Int, var_idxs::Vector, operator::Symbol) end """ -Alp.expr_constr_parsing(expr, m::Optimizer) +expr_constr_parsing(expr, m::Optimizer) Recognize structural constraints. """ function expr_constr_parsing(expr::Expr, m::Optimizer, idx::Int = 0) # First process user-defined structures in-cases of over-ride - # for i in 1:length(Alp.get_option(m, :constr_patterns)) - # is_structural = eval(Alp.get_option(m, :constr_patterns)[i])(expr, m, idx) + # for i in 1:length(get_option(m, :constr_patterns)) + # is_structural = eval(get_option(m, :constr_patterns)[i])(expr, m, idx) # return # end isa(expr, Number) && return false # Recognize built-in special structural pattern - if Alp.get_option(m, :recognize_convex) - is_convex = Alp.resolve_convex_constr(expr, m, idx) + if get_option(m, :recognize_convex) + is_convex = resolve_convex_constr(expr, m, idx) is_convex && return true end @@ -241,7 +239,7 @@ function expr_is_axn(expr, scalar = 1.0, var_idxs = [], power = []; N = nothing) push!(power, 1) elseif (expr.args[i].head == :call) scalar, var_idxs, power = - Alp.expr_is_axn(expr.args[i], scalar, var_idxs, power) + expr_is_axn(expr.args[i], scalar, var_idxs, power) scalar === nothing && return nothing, nothing, nothing end end @@ -284,8 +282,7 @@ function expr_linear_to_affine(expr) @assert isa(expr.args[3], Float64) || isa(expr.args[3], Int) @assert isa(expr.args[2], Expr) # non are buffer spaces, not used anywhere - lhscoeff, lhsvars, rhs, non, non = - Alp.traverse_expr_linear_to_affine(expr.args[2]) + lhscoeff, lhsvars, rhs, non, non = traverse_expr_linear_to_affine(expr.args[2]) rhs = -rhs + expr.args[3] affdict[:sense] = expr.args[1] @@ -296,7 +293,7 @@ function expr_linear_to_affine(expr) affdict[:sense] = nothing else # For an objective expression - lhscoeff, lhsvars, rhs, non, non = Alp.traverse_expr_linear_to_affine(expr) + lhscoeff, lhsvars, rhs, non, non = traverse_expr_linear_to_affine(expr) affdict[:sense] = nothing end @@ -446,7 +443,7 @@ function expr_resolve_sign(expr, level = 0; kwargs...) expr.args[i] = expr.args[i].args[2] end elseif expr.args[i].head == :call - Alp.expr_resolve_sign(expr.args[i], level + 1) + expr_resolve_sign(expr.args[i], level + 1) end end end @@ -462,7 +459,7 @@ Most issues can be caused by this function. function expr_flatten(expr, level = 0; kwargs...) isa(expr, Number) && return if level > 0 # No trivial constraint is allowed "3>5" - flat = Alp.expr_arrangeargs(expr.args) + flat = expr_arrangeargs(expr.args) if isa(flat, Number) return flat else @@ -477,7 +474,7 @@ function expr_flatten(expr, level = 0; kwargs...) end if level > 0 #Root level process, no additional processes - expr.args = Alp.expr_arrangeargs(expr.args) + expr.args = expr_arrangeargs(expr.args) end return expr @@ -591,14 +588,14 @@ function expr_resolve_const(expr) if isa(expr.args[i], Number) || isa(expr.args[i], Symbol) continue elseif expr.args[i].head == :call - if Alp.expr_is_emptysum(expr.args[i]) + if expr_is_emptysum(expr.args[i]) expr.args[i] = :(0) end (expr_isconst(expr.args[i])) && (expr.args[i] = eval(expr.args[i])) if isa(expr.args[i], Number) || isa(expr.args[i], Symbol) continue else - Alp.expr_resolve_const(expr.args[i]) + expr_resolve_const(expr.args[i]) end end end @@ -615,7 +612,7 @@ function expr_resolve_const_denominator(args) resolvable = true for i in 3:length(args) - resolvable *= Alp.expr_isconst(args[i]) + resolvable *= expr_isconst(args[i]) end if resolvable @@ -663,7 +660,7 @@ function expr_is_fractional_exponent(expr) end for i in 1:length(expr.args) if typeof(expr.args[i]) == Expr - Alp.expr_is_fractional_exponent(expr.args[i]) # Recursively search for fractional exponents + expr_is_fractional_exponent(expr.args[i]) # Recursively search for fractional exponents end end end @@ -671,7 +668,7 @@ end # Returns true if the expression is a constant, linear or affine function expr_isaffine(expr) - Alp.expr_isconst(expr) && return true + expr_isconst(expr) && return true expr.head == :ref && return true is_affine = false @@ -689,7 +686,7 @@ function expr_isaffine(expr) if expr.args[i].head == :ref k += 1 elseif expr.args[i].head == :call - status = Alp.expr_isaffine(expr.args[i]) + status = expr_isaffine(expr.args[i]) status && (k += 1) end end @@ -704,7 +701,7 @@ Signs in the summation can be +/- Note: This function does not support terms of type (a*(x[1] + x[2]))^2 yet. """ function expr_isolate_const(expr) - Alp.expr_isaffine(expr) && return expr + expr_isaffine(expr) && return expr # Check nonlinear expressions if (expr.head == :call && expr.args[1] in [:+, :-]) @@ -770,7 +767,7 @@ function expr_isolate_const(expr) # Handle negative sign in the remaining terms elseif (expr_i.args[1] in [:+, :-]) && (length(expr.args[i].args) > 2) - expr_rec = Alp.expr_isolate_const(expr_i) #recursion + expr_rec = expr_isolate_const(expr_i) #recursion push!(expr_array, expr_rec) # For any other terms diff --git a/src/operators.jl b/src/operators.jl index daff8d49..b3c93560 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -11,7 +11,7 @@ function expr_term_parsing(expr::Any, constr_id::Int, m::Optimizer, level = 0; o if isa(node, Float64) || isa(node, Int64) || isa(node, Symbol) continue elseif node.head == :call - expr.args[cnt] = Alp.expr_term_parsing(node, constr_id, m, level + 1) + expr.args[cnt] = expr_term_parsing(node, constr_id, m, level + 1) elseif node.head == :ref continue else @@ -19,7 +19,7 @@ function expr_term_parsing(expr::Any, constr_id::Int, m::Optimizer, level = 0; o end end - return Alp.detect_nonconvex_terms(expr, constr_id, m) + return detect_nonconvex_terms(expr, constr_id, m) end """ @@ -35,8 +35,8 @@ Specific structure pattern information will be described formally. function detect_nonconvex_terms(expr::Any, constr_id::Int, m::Optimizer; kwargs...) # First process user-defined structures in-cases of over-ride - # for i in 1:length(Alp.get_option(m, :term_patterns)) - # skip, expr = eval(Alp.get_option(m, :term_patterns)[i])(expr, constr_id, m) + # for i in 1:length(get_option(m, :term_patterns)) + # skip, expr = eval(get_option(m, :term_patterns)[i])(expr, constr_id, m) # skip && return expr # end @@ -46,20 +46,20 @@ function detect_nonconvex_terms(expr::Any, constr_id::Int, m::Optimizer; kwargs. # general multilinear terms to apply better outter approxiamtion. # LEVEL 1 : : Recognize all built-in structural patterns - skip, expr = Alp.detect_discretemulti_term(expr, constr_id, m) #L1 : General case : discrete variables * continous variables + skip, expr = detect_discretemulti_term(expr, constr_id, m) #L1 : General case : discrete variables * continous variables skip && return expr - skip, expr = Alp.detect_binprod_term(expr, constr_id, m) #L1 : binprod = binary products + skip, expr = detect_binprod_term(expr, constr_id, m) #L1 : binprod = binary products skip && return expr # LEVEL 2 : : Recognize all built-in structural patterns - skip, expr = Alp.detect_bilinear_term(expr, constr_id, m) #L2 + skip, expr = detect_bilinear_term(expr, constr_id, m) #L2 skip && return expr - skip, expr = Alp.detect_monomial_term(expr, constr_id, m) #L2 : must go before multilinear + skip, expr = detect_monomial_term(expr, constr_id, m) #L2 : must go before multilinear skip && return expr - skip, expr = Alp.detect_multilinear_term(expr, constr_id, m) #L2 + skip, expr = detect_multilinear_term(expr, constr_id, m) #L2 skip && return expr return expr # if no structure is detected, simply return the original tree @@ -80,14 +80,14 @@ function store_nonconvex_term( y_idx = m.num_var_orig + nl_cnt + l_cnt + 1 # y is always lifted var lifted_var_ref = Expr(:ref, :x, y_idx) - lifted_constr_ref = Alp.build_constr_block(y_idx, var_idxs, operator) + lifted_constr_ref = build_constr_block(y_idx, var_idxs, operator) m.nonconvex_terms[nl_key] = Dict( :lifted_var_ref => lifted_var_ref, :id => nl_cnt + 1, :y_idx => y_idx, :y_type => - Alp.resolve_lifted_var_type([m.var_type[k] for k in var_idxs], operator), + resolve_lifted_var_type([m.var_type[k] for k in var_idxs], operator), :var_idxs => var_idxs, :ref => nl_key, :evaluator => evaluator, @@ -103,7 +103,7 @@ function store_nonconvex_term( # push!(m.var_type, :Cont) # TODO check if this replacement is good since additional constraints should be able to sufficiently constraint the type push!(m.var_type, m.nonconvex_terms[nl_key][:y_type]) # Keep track of the lifted var type - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("found lifted $(term_type) term $(lifted_constr_ref)") return y_idx end @@ -122,7 +122,7 @@ function store_linear_term(m::Optimizer, term_key::Any, expr::Any)#, bound_resol :id => length(keys(m.linear_terms)) + 1, :ref => term_key, :y_idx => y_idx, - :y_type => Alp.resolve_lifted_var_type( + :y_type => resolve_lifted_var_type( [m.var_type[k[2]] for k in term_key[:coef_var]], :+, ), @@ -134,7 +134,7 @@ function store_linear_term(m::Optimizer, term_key::Any, expr::Any)#, bound_resol m.term_seq[l_cnt+nl_cnt+1] = term_key push!(m.var_type, m.linear_terms[term_key][:y_type]) # Keep track of the lifted var type - Alp.get_option(m, :log_level) > 199 && + get_option(m, :log_level) > 199 && println("found lifted linear term $(lifted_var_ref) = $expr") return y_idx @@ -162,9 +162,9 @@ function detect_linear_term(expr::Any, constr_id::Int, m::Optimizer) coef_fetch = Dict(:+ => 1.0, :- => -1.0) # Re-process the expression sub-tree [REQUIRED] - Alp.expr_resolve_const(expr) - Alp.expr_resolve_sign(expr) - Alp.expr_flatten(expr) + expr_resolve_const(expr) + expr_resolve_sign(expr) + expr_flatten(expr) if expr.args[1] in [:+, :-] scalar = 0.0 @@ -214,15 +214,15 @@ function detect_linear_term(expr::Any, constr_id::Int, m::Optimizer) push!(coef_var, (coef_fetch[expr.args[1]] * sub_coef[1], sub_vars[1])) else down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) # Recursive detection + detect_linear_term(expr.args[i], constr_id, m) # Recursive detection down_check ? expr.args[i] = linear_lift_var : return false, expr push!(coef_var, (1.0, expr.args[i].args[2])) end end # By reaching here, it is already certain that we have found the term, always treat with :+ term_key = Dict(:scalar => scalar, :coef_var => coef_var, :sign => :+) - term_key in keys(m.linear_terms) || Alp.store_linear_term(m, term_key, expr) - return true, Alp.lift_linear_term(m, term_key, constr_id) + term_key in keys(m.linear_terms) || store_linear_term(m, term_key, expr) + return true, lift_linear_term(m, term_key, constr_id) elseif expr.args[1] in [:*] && length(expr.args) == 3 # For terms like (3*x)*y scalar = 0.0 @@ -239,8 +239,8 @@ function detect_linear_term(expr::Any, constr_id::Int, m::Optimizer) # By reaching here, it is already certain that we have found the term, always treat with :+ term_key = Dict(:scalar => scalar, :coef_var => coef_var, :sign => :+) - term_key in keys(m.linear_terms) || Alp.store_linear_term(m, term_key, expr) - return true, Alp.lift_linear_term(m, term_key, constr_id) + term_key in keys(m.linear_terms) || store_linear_term(m, term_key, expr) + return true, lift_linear_term(m, term_key, constr_id) end return false, expr @@ -260,10 +260,10 @@ function basic_linear_bounds(m::Optimizer, k::Any, linear_terms = nothing) end lb += linear_terms[k][:ref][:scalar] ub += linear_terms[k][:ref][:scalar] - if lb > m.l_var_tight[lifted_idx] + Alp.get_option(m, :tol) + if lb > m.l_var_tight[lifted_idx] + get_option(m, :tol) m.l_var_tight[lifted_idx] = lb end - if ub < m.u_var_tight[lifted_idx] - Alp.get_option(m, :tol) + if ub < m.u_var_tight[lifted_idx] - get_option(m, :tol) m.u_var_tight[lifted_idx] = ub end @@ -281,10 +281,10 @@ function basic_linear_bounds(m::Optimizer, k::Any, d::Dict) lb += m.linear_terms[k][:ref][:scalar] ub += m.linear_terms[k][:ref][:scalar] - if lb > d[lifted_idx][1] + Alp.get_option(m, :tol) + if lb > d[lifted_idx][1] + get_option(m, :tol) d[lifted_idx][1] = lb end - if ub < d[lifted_idx][end] - Alp.get_option(m, :tol) + if ub < d[lifted_idx][end] - get_option(m, :tol) d[lifted_idx][end] = ub end @@ -335,7 +335,7 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) push!(var_types, m.var_type[expr.args[i].args[2]]) if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr push!(var_idxs, linear_lift_var.args[2]) push!(var_types, m.var_type[linear_lift_var.args[2]]) @@ -360,7 +360,7 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) for idx in cont_var_idxs push!(ml_term_expr.args, Expr(:ref, :x, idx)) end - ml_lift_term = Alp.detect_nonconvex_terms(ml_term_expr, constr_id, m) + ml_lift_term = detect_nonconvex_terms(ml_term_expr, constr_id, m) ml_idx = ml_lift_term.args[2] else ml_idx = cont_var_idxs[1] @@ -373,7 +373,7 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) for idx in bin_var_idxs push!(bp_term_expr.args, Expr(:ref, :x, idx)) end - bp_lift_term = Alp.detect_nonconvex_terms(bp_term_expr, constr_id, m) + bp_lift_term = detect_nonconvex_terms(bp_term_expr, constr_id, m) bp_idx = bp_lift_term.args[2] else isempty(bin_var_idxs) ? bp_idx = -1 : bp_idx = bin_var_idxs[1] @@ -386,7 +386,7 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) for idx in int_var_idxs push!(ip_term_expr.args, Expr(:ref, :x, idx)) end - ip_lift_term = Alp.detect_nonconvex_terms(ip_term_expr, constr_id, m) + ip_lift_term = detect_nonconvex_terms(ip_term_expr, constr_id, m) ip_idx = ip_lift_term.args[2] else isempty(int_var_idxs) ? ip_idx = -1 : ip_idx = int_var_idxs[1] @@ -395,7 +395,7 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) if ip_idx < 0 # binlin term if no integer variable binlin_key = [Expr(:ref, :x, bp_idx), Expr(:ref, :x, ml_idx)] binlin_idxs = [bp_idx; ml_idx] - binlin_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + binlin_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, binlin_key, binlin_idxs, @@ -405,14 +405,14 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) basic_binlin_bounds, collect_binlin_discvar, ) - return true, Alp.lift_nonconvex_term(m, binlin_key, constr_id, scalar) + return true, lift_nonconvex_term(m, binlin_key, constr_id, scalar) end # binlin_key = [Expr(:ref, :x, bp_idx), Expr(:ref, :x, intlin_idx)] binlin_key = [Expr(:ref, :x, bp_idx)] # binlin_idxs = [bp_idx; intlin_idx] binlin_idxs = [bp_idx] - binlin_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + binlin_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, binlin_key, binlin_idxs, @@ -422,7 +422,7 @@ function detect_discretemulti_term(expr::Any, constr_id::Int, m::Optimizer) basic_binlin_bounds, collect_binlin_discvar, ) - return true, Alp.lift_nonconvex_term(m, binlin_key, constr_id, scalar) + return true, lift_nonconvex_term(m, binlin_key, constr_id, scalar) end return false, expr @@ -496,7 +496,7 @@ function detect_binprod_term(expr::Any, constr_id::Int, m::Optimizer) !isempty(var_idxs) && m.var_type[var_idxs[end]] != :Bin && return false, expr if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr m.var_type[linear_lift_var.args[2]] != :Bin && return false, expr push!(var_idxs, linear_lift_var.args[2]) @@ -505,7 +505,7 @@ function detect_binprod_term(expr::Any, constr_id::Int, m::Optimizer) end if length(var_idxs) >= 2 term_key = [Expr(:ref, :x, idx) for idx in var_idxs] - term_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + term_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, term_key, var_idxs, @@ -515,7 +515,7 @@ function detect_binprod_term(expr::Any, constr_id::Int, m::Optimizer) basic_binprod_bounds, collect_binprod_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end elseif (expr.args[1] == :^) && length(expr.args) == 3 # Pattern: (x)^(>2), where x is binary variable @@ -534,7 +534,7 @@ function detect_binprod_term(expr::Any, constr_id::Int, m::Optimizer) !isempty(var_idxs) && m.var_type[var_idxs[end]] != :Bin && return false, expr if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr m.var_type[linear_lift_var.args[2]] != :Bin && return false, expr push!(var_idxs, linear_lift_var.args[2]) @@ -543,7 +543,7 @@ function detect_binprod_term(expr::Any, constr_id::Int, m::Optimizer) end if length(var_idxs) == 1 && power_scalar > 2.0 && mod(power_scalar, 1.0) == 0.0 term_key = [Expr(:ref, :x, var_idxs[1]) for i in 1:power_scalar] - term_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + term_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, term_key, var_idxs, @@ -553,7 +553,7 @@ function detect_binprod_term(expr::Any, constr_id::Int, m::Optimizer) basic_binprod_bounds, collect_binprod_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end end @@ -608,7 +608,7 @@ function detect_bilinear_term(expr::Any, constr_id::Int, m::Optimizer) return false, expr # Don't consider the discrete variable if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr push!(var_idxs, linear_lift_var.args[2]) m.var_type[var_idxs[end]] in [:Bin, :Int] && return false, expr # Don't consider the discrete variable @@ -623,19 +623,19 @@ function detect_bilinear_term(expr::Any, constr_id::Int, m::Optimizer) reverse(term_key) in keys(m.nonconvex_terms) term_key in keys(m.nonconvex_terms) ? term_key = term_key : term_key = reverse(term_key) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) else - Alp.store_nonconvex_term( + store_nonconvex_term( m, term_key, var_idxs, :BILINEAR, :*, bilinear, - Alp.basic_monomial_bounds, + basic_monomial_bounds, collect_monomial_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end end end @@ -662,7 +662,7 @@ function detect_multilinear_term(expr::Any, constr_id::Int, m::Optimizer) return false, expr if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr push!(var_idxs, linear_lift_var.args[2]) m.var_type[var_idxs[end]] in [:Bin, :Int] && return false, expr @@ -671,17 +671,17 @@ function detect_multilinear_term(expr::Any, constr_id::Int, m::Optimizer) end if length(var_idxs) > 2 term_key = [Expr(:ref, :x, idx) for idx in var_idxs] - term_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + term_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, term_key, var_idxs, :MULTILINEAR, :*, multilinear, - Alp.basic_monomial_bounds, + basic_monomial_bounds, collect_monomial_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end elseif (expr.args[1] == :^) && length(expr.args) == 3 # Pattern: (x)^(>2) var_idxs = [] @@ -701,7 +701,7 @@ function detect_multilinear_term(expr::Any, constr_id::Int, m::Optimizer) return false, expr # Avoid fetching terms with discrete variables if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr push!(var_idxs, linear_lift_var.args[2]) m.var_type[var_idxs[end]] in [:Bin, :Int] && return false, expr @@ -710,17 +710,17 @@ function detect_multilinear_term(expr::Any, constr_id::Int, m::Optimizer) end if length(var_idxs) == 1 && power_scalar > 2.0 && mod(power_scalar, 1.0) == 0.0 term_key = [Expr(:ref, :x, var_idxs[1]) for i in 1:power_scalar] - term_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + term_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, term_key, var_idxs, :MULTILINEAR, :*, multilinear, - Alp.basic_monomial_bounds, + basic_monomial_bounds, collect_monomial_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end end @@ -747,7 +747,7 @@ function detect_monomial_term(expr::Any, constr_id::Int, m::Optimizer) return false, expr # Avoid fetching terms with discrete variables if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr push!(var_idxs, linear_lift_var.args[2]) m.var_type[var_idxs[end]] in [:Bin, :Int] && return false, expr # Avoid fetching terms with discrete variables @@ -762,17 +762,17 @@ function detect_monomial_term(expr::Any, constr_id::Int, m::Optimizer) if length(var_idxs) == 1 && power_scalar == 2.0 term_key = [Expr(:ref, :x, var_idxs[1]) for i in 1:2] - term_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + term_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, term_key, var_idxs, :MONOMIAL, :*, monomial, - Alp.basic_monomial_bounds, + basic_monomial_bounds, collect_monomial_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end end @@ -795,7 +795,7 @@ function detect_monomial_term(expr::Any, constr_id::Int, m::Optimizer) return false, expr # Avoid fetching terms with discrete variables if (expr.args[i].head == :call) down_check, linear_lift_var = - Alp.detect_linear_term(expr.args[i], constr_id, m) + detect_linear_term(expr.args[i], constr_id, m) !down_check && return false, expr push!(var_idxs, linear_lift_var.args[2]) m.var_type[var_idxs[end]] in [:Bin, :Int] && return false, expr # Avoid fetching terms with discrete variables @@ -805,17 +805,17 @@ function detect_monomial_term(expr::Any, constr_id::Int, m::Optimizer) # Confirm detection of pattern A and perform store & lifting procedures if (length(var_idxs) == 2) && (length(Set(var_idxs)) == 1) term_key = [Expr(:ref, :x, var_idxs[1]) for i in 1:2] - term_key in keys(m.nonconvex_terms) || Alp.store_nonconvex_term( + term_key in keys(m.nonconvex_terms) || store_nonconvex_term( m, term_key, var_idxs, :MONOMIAL, :*, monomial, - Alp.basic_monomial_bounds, + basic_monomial_bounds, collect_monomial_discvar, ) - return true, Alp.lift_nonconvex_term(m, term_key, constr_id, scalar) + return true, lift_nonconvex_term(m, term_key, constr_id, scalar) end end @@ -838,13 +838,13 @@ function basic_monomial_bounds(m::Optimizer, k::Any) bound = vec(bound) * var_bounds' end end - if minimum(bound) > m.l_var_tight[lifted_idx] + Alp.get_option(m, :tol) + if minimum(bound) > m.l_var_tight[lifted_idx] + get_option(m, :tol) m.l_var_tight[lifted_idx] = minimum(bound) if m.nonconvex_terms[k][:nonlinear_type] == :MONOMIAL m.l_var_tight[lifted_idx] = 0.0 end end - if maximum(bound) < m.u_var_tight[lifted_idx] - Alp.get_option(m, :tol) + if maximum(bound) < m.u_var_tight[lifted_idx] - get_option(m, :tol) m.u_var_tight[lifted_idx] = maximum(bound) end @@ -868,11 +868,11 @@ function basic_monomial_bounds(m::Optimizer, nlk::Any, d::Dict) end end - if minimum(bound) > d[lifted_idx][1] + Alp.get_option(m, :tol) + if minimum(bound) > d[lifted_idx][1] + get_option(m, :tol) d[lifted_idx][1] = minimum(bound) end - if maximum(bound) < d[lifted_idx][end] - Alp.get_option(m, :tol) + if maximum(bound) < d[lifted_idx][end] - get_option(m, :tol) d[lifted_idx][end] = maximum(bound) end @@ -909,13 +909,13 @@ function resolve_convex_constr( expr_orig = :constr sense = expr.args[1] rhs = expr.args[3] - subs, rhs_strip = Alp.expr_strip_const(expr.args[2]) # Focus on the regularized subtree (stripped with constants) + subs, rhs_strip = expr_strip_const(expr.args[2]) # Focus on the regularized subtree (stripped with constants) rhs += rhs_strip # TODO: check if sign is flipped elseif expr.args[1] == :(==) return false elseif idx == 0 expr_orig = :obj - subs, rhs = Alp.expr_strip_const(expr) # Focus on the regularized subtree (stripped with constants) + subs, rhs = expr_strip_const(expr) # Focus on the regularized subtree (stripped with constants) end for sub in subs @@ -932,7 +932,7 @@ function resolve_convex_constr( elseif (sub.args[i].head == :ref) return false elseif (sub.args[i].head == :call) - scalar, idxs, powers = Alp.expr_is_axn(sub.args[i]) + scalar, idxs, powers = expr_is_axn(sub.args[i]) (scalar === nothing) && return false if length(Set(idxs)) == 1 push!(idxs_bin, idxs[1]) @@ -947,7 +947,7 @@ function resolve_convex_constr( end end elseif sub.args[1] in [:*] - scalar, idxs, powers = Alp.expr_is_axn(sub) + scalar, idxs, powers = expr_is_axn(sub) (scalar === nothing) && return false if length(Set(idxs)) == 1 push!(idxs_bin, idxs[1]) @@ -1045,7 +1045,7 @@ function resolve_convex_constr( :powers => power_bin, ) - Alp.get_option(m, :log_level) > 99 && println("CONVEX Constraint $(idx): $(expr)") + get_option(m, :log_level) > 99 && println("CONVEX Constraint $(idx): $(expr)") return true elseif expr_orig == :obj @@ -1054,7 +1054,7 @@ function resolve_convex_constr( (expr_isconst(m.obj_expr_orig)) && return true # Follows the same mapping to convex constraints - # Important: This is a fix to return obj is non-convex if the exponents have value greater than 2. + # Important: This is a fix to return obj is non-convex if the exponents have value greater than 2. # This is needs to be removed once outer-approximation for convex functions is implemented (maximum(power_bin) > 2) && return false @@ -1109,7 +1109,7 @@ function resolve_convex_constr( :powers => power_bin, ) - Alp.get_option(m, :log_level) > 99 && println("CONVEX Objective: $(expr)") + get_option(m, :log_level) > 99 && println("CONVEX Objective: $(expr)") return true end diff --git a/src/presolve.jl b/src/presolve.jl index abc04a6b..f9ca770f 100644 --- a/src/presolve.jl +++ b/src/presolve.jl @@ -11,15 +11,15 @@ Currently, two OBBT methods are implemented in [`optimization_based_bound_tighte If no local feasible solution is obtained, the algorithm defaults to OBBT without partitions """ function bound_tightening_wrapper(m::Optimizer; use_bound = true, kwargs...) - Alp.get_option(m, :presolve_bt) || return - - if Alp.get_option(m, :presolve_bt_algo) == 1 - Alp.optimization_based_bound_tightening(m, use_bound = use_bound) - elseif Alp.get_option(m, :presolve_bt_algo) == 2 - Alp.optimization_based_bound_tightening(m, use_bound = use_bound, use_tmc = true) - elseif isa(Alp.get_option(m, :presolve_bt_algo), Function) - # eval(Alp.get_option(m, :presolve_bt_algo))(m) - Alp.get_option(m, :presolve_bt_algo)(m) + get_option(m, :presolve_bt) || return + + if get_option(m, :presolve_bt_algo) == 1 + optimization_based_bound_tightening(m, use_bound = use_bound) + elseif get_option(m, :presolve_bt_algo) == 2 + optimization_based_bound_tightening(m, use_bound = use_bound, use_tmc = true) + elseif isa(get_option(m, :presolve_bt_algo), Function) + # eval(get_option(m, :presolve_bt_algo))(m) + get_option(m, :presolve_bt_algo)(m) else error("Unrecognized bound tightening algorithm") end @@ -64,7 +64,7 @@ function optimization_based_bound_tightening( st = time() # Track start time if time_limit == Inf - time_limit = Alp.get_option(m, :presolve_bt_time_limit) + time_limit = get_option(m, :presolve_bt_time_limit) end # Regulating special input conditions: default use best feasible solution objective value @@ -72,36 +72,35 @@ function optimization_based_bound_tightening( l_var_orig = copy(m.l_var_tight) u_var_orig = copy(m.u_var_tight) - discretization = Alp._get_discretization_dict(m, m.l_var_tight, m.u_var_tight) + discretization = _get_discretization_dict(m, m.l_var_tight, m.u_var_tight) if use_bound == false && haskey(options, :use_tmc) - (Alp.get_option(m, :log_level) > 0) && + (get_option(m, :log_level) > 0) && @warn " Local solve infeasible; defaulting to doing bound-tightening without partitions." end if use_bound == true && haskey(options, :use_tmc) - discretization = Alp.add_adaptive_partition( + discretization = add_adaptive_partition( m, use_solution = m.best_sol, use_disc = discretization, ) end - discretization = Alp.resolve_var_bounds(m, discretization) # recomputation of bounds for lifted_variables + discretization = resolve_var_bounds(m, discretization) # recomputation of bounds for lifted_variables - (Alp.get_option(m, :log_level) > 0) && println(" Starting bound-tightening") + (get_option(m, :log_level) > 0) && println(" Starting bound-tightening") - width_tol = Alp.get_option(m, :presolve_bt_width_tol) - bound_tol = Alp.get_option(m, :presolve_bt_bound_tol) - improv_tol = Alp.get_option(m, :presolve_bt_improv_tol) + width_tol = get_option(m, :presolve_bt_width_tol) + bound_tol = get_option(m, :presolve_bt_bound_tol) + improv_tol = get_option(m, :presolve_bt_improv_tol) keep_tightening = true # start of the solve while keep_tightening && - (m.logs[:time_left] > Alp.get_option(m, :tol)) && - (m.logs[:bt_iter] < Alp.get_option(m, :presolve_bt_max_iter)) # Stopping criteria + (m.logs[:time_left] > get_option(m, :tol)) && + (m.logs[:bt_iter] < get_option(m, :presolve_bt_max_iter)) # Stopping criteria keep_tightening = false m.logs[:bt_iter] += 1 - Alp.get_option(m, :log_level) > 199 && - println(" Iteration - $(m.logs[:bt_iter])") + get_option(m, :log_level) > 199 && println(" Iteration - $(m.logs[:bt_iter])") temp_bounds = Dict() avg_reduction = 0.0 @@ -115,14 +114,14 @@ function optimization_based_bound_tightening( temp_bounds[var_idx] = [discretization[var_idx][1], discretization[var_idx][end]] if (discretization[var_idx][end] - discretization[var_idx][1]) > width_tol - Alp.create_obbt_model(m, discretization, bound) + create_obbt_model(m, discretization, bound) for sense in both_senses JuMP.@objective( m.model_mip, sense, _index_to_variable_ref(m.model_mip, var_idx) ) - stats = Alp.solve_obbt_model(m) + stats = solve_obbt_model(m) if stats["status"] in STATUS_OPT temp_bounds[var_idx][tell_side[sense]] = tell_round[sense]( @@ -200,12 +199,12 @@ function optimization_based_bound_tightening( max_reduction = max(bound_reduction, max_reduction) max_width = max(new_range, max_width) - (Alp.get_option(m, :log_level) > 99) && print("+") - (Alp.get_option(m, :log_level) > 99) && println( + (get_option(m, :log_level) > 99) && print("+") + (get_option(m, :log_level) > 99) && println( " VAR $(var_idx) LB contracted $(discretization[var_idx][1])=>$(temp_bounds[var_idx][1])", ) - (Alp.get_option(m, :log_level) > 99) && print("+") - (Alp.get_option(m, :log_level) > 99) && println( + (get_option(m, :log_level) > 99) && print("+") + (get_option(m, :log_level) > 99) && println( " VAR $(var_idx) UB contracted $(discretization[var_idx][end])=>$(temp_bounds[var_idx][end])", ) @@ -219,45 +218,44 @@ function optimization_based_bound_tightening( bound_max_reduction = (max_reduction > improv_tol) bound_max_width = (max_width > width_tol) - # Deactivate this termination criterion if it slows down the OBBT convergence - stats = Alp.relaxation_model_obbt(m, discretization, bound) - if Alp.is_min_sense(m) - current_rel_gap = Alp.eval_opt_gap(m, stats["relaxed_obj"], bound) - elseif Alp.is_max_sense(m) - current_rel_gap = Alp.eval_opt_gap(m, bound, stats["relaxed_obj"]) + # Deactivate this termination criterion if it slows down the OBBT convergence + stats = relaxation_model_obbt(m, discretization, bound) + if is_min_sense(m) + current_rel_gap = eval_opt_gap(m, stats["relaxed_obj"], bound) + elseif is_max_sense(m) + current_rel_gap = eval_opt_gap(m, bound, stats["relaxed_obj"]) end keep_tightening = (bound_avg_reduction) && (bound_max_reduction) && (bound_max_width) && - (current_rel_gap > Alp.get_option(m, :rel_gap)) + (current_rel_gap > get_option(m, :rel_gap)) if !isinf(current_rel_gap) m.presolve_best_rel_gap = current_rel_gap * 100 end - discretization = Alp.resolve_var_bounds(m, discretization) + discretization = resolve_var_bounds(m, discretization) if haskey(options, :use_tmc) - discretization = Alp.add_adaptive_partition( + discretization = add_adaptive_partition( m, use_solution = m.best_sol, - use_disc = Alp.flatten_discretization(discretization), + use_disc = flatten_discretization(discretization), ) end time() - st > time_limit && break end - (Alp.get_option(m, :log_level) > 1) && - println(" Variables whose bounds were tightened:") - (Alp.get_option(m, :log_level) > 0) && + (get_option(m, :log_level) > 1) && println(" Variables whose bounds were tightened:") + (get_option(m, :log_level) > 0) && println(" Actual iterations (OBBT): ", (m.logs[:bt_iter])) - m.l_var_tight, m.u_var_tight = Alp.update_var_bounds(discretization) + m.l_var_tight, m.u_var_tight = update_var_bounds(discretization) if haskey(options, :use_tmc) - m.discretization = Alp.add_adaptive_partition(m, use_solution = m.best_sol) + m.discretization = add_adaptive_partition(m, use_solution = m.best_sol) end for i in m.disc_vars @@ -268,8 +266,8 @@ function optimization_based_bound_tightening( abs(l_var_orig[i] - u_var_orig[i]); digits = 2, ) * 100 - if Alp.get_option(m, :log_level) > 0 && contract_ratio > 0.0001 - (Alp.get_option(m, :log_level) > 1) && (println( + if get_option(m, :log_level) > 0 && contract_ratio > 0.0001 + (get_option(m, :log_level) > 1) && (println( " VAR $(i): $(contract_ratio)% contraction |$(round(l_var_orig[i]; digits=4)) --> | $(round(m.l_var_tight[i]; digits=4)) - $(round(m.u_var_tight[i]; digits=4)) | <-- $(round(u_var_orig[i]; digits=4)) |", )) end @@ -287,26 +285,26 @@ function create_obbt_model(m::Optimizer, discretization, bound::Number; kwargs.. # options = Dict(kwargs) start_build = time() - m.model_mip = Model(Alp.get_option(m, :mip_solver)) # Construct JuMP model - Alp.amp_post_vars(m, use_disc = discretization) - Alp.amp_post_lifted_constraints(m) - Alp.amp_post_convexification(m, use_disc = discretization) # Convexify problem + m.model_mip = Model(get_option(m, :mip_solver)) # Construct JuMP model + amp_post_vars(m, use_disc = discretization) + amp_post_lifted_constraints(m) + amp_post_convexification(m, use_disc = discretization) # Convexify problem if !(isinf(bound)) - Alp.post_objective_bound(m, bound) + post_objective_bound(m, bound) else @warn "Dropping the objective bound constraint in presolve bound tightening" end cputime_build = time() - start_build m.logs[:total_time] += cputime_build - m.logs[:time_left] = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + m.logs[:time_left] = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) return end function relaxation_model_obbt(m::Optimizer, discretization, bound::Number) - Alp.create_obbt_model(m, discretization, bound) + create_obbt_model(m, discretization, bound) obj_expr = sum( m.bounding_obj_mip[:coefs][j] * @@ -314,13 +312,13 @@ function relaxation_model_obbt(m::Optimizer, discretization, bound::Number) j in 1:m.bounding_obj_mip[:cnt] ) - if Alp.is_min_sense(m) + if is_min_sense(m) JuMP.@objective(m.model_mip, Min, obj_expr) - elseif Alp.is_max_sense(m) + elseif is_max_sense(m) JuMP.@objective(m.model_mip, Max, obj_expr) end - return Alp.solve_obbt_model(m) + return solve_obbt_model(m) end """ @@ -334,19 +332,19 @@ function solve_obbt_model(m::Optimizer; kwargs...) # ========= MILP Solve ========= # time_limit = max( 0.0, - if Alp.get_option(m, :presolve_bt_mip_time_limit) < Inf + if get_option(m, :presolve_bt_mip_time_limit) < Inf min( - Alp.get_option(m, :presolve_bt_mip_time_limit), - Alp.get_option(m, :time_limit) - m.logs[:total_time], + get_option(m, :presolve_bt_mip_time_limit), + get_option(m, :time_limit) - m.logs[:total_time], ) else - Alp.get_option(m, :time_limit) - m.logs[:total_time] + get_option(m, :time_limit) - m.logs[:total_time] end, ) MOI.set(m.model_mip, MOI.TimeLimitSec(), time_limit) start_solve = time() - if Alp.get_option(m, :presolve_bt_relax_integrality) + if get_option(m, :presolve_bt_relax_integrality) JuMP.relax_integrality(m.model_mip) end @@ -358,7 +356,7 @@ function solve_obbt_model(m::Optimizer; kwargs...) cputime_solve = time() - start_solve m.logs[:total_time] += cputime_solve - m.logs[:time_left] = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + m.logs[:time_left] = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) # ========= MILP Solve ========= # return stats @@ -367,8 +365,8 @@ end """ post_objective_bound(m::Optimizer, bound::Float64; kwargs...) -This function adds the upper/lower bounding constraint on the objective function -for the optimization models solved within the OBBT algorithm. +This function adds the upper/lower bounding constraint on the objective function +for the optimization models solved within the OBBT algorithm. """ function post_objective_bound(m::Optimizer, bound::Number; kwargs...) obj_expr = sum( @@ -377,14 +375,14 @@ function post_objective_bound(m::Optimizer, bound::Number; kwargs...) j in 1:m.bounding_obj_mip[:cnt] ) - obj_bound_tol = Alp.get_option(m, :presolve_bt_obj_bound_tol) + obj_bound_tol = get_option(m, :presolve_bt_obj_bound_tol) - if Alp.is_max_sense(m) + if is_max_sense(m) JuMP.@constraint( m.model_mip, obj_expr >= floor(bound / obj_bound_tol) * obj_bound_tol ) - elseif Alp.is_min_sense(m) + elseif is_min_sense(m) JuMP.@constraint( m.model_mip, obj_expr <= ceil(bound / obj_bound_tol) * obj_bound_tol diff --git a/src/utility.jl b/src/utility.jl index 5d6cb5ad..ac01d22a 100644 --- a/src/utility.jl +++ b/src/utility.jl @@ -15,13 +15,13 @@ The absolute gap calculation is ``` """ function update_opt_gap(m::Optimizer) - tol = Alp.get_option(m, :tol) + tol = get_option(m, :tol) if m.best_obj in [Inf, -Inf] m.best_rel_gap = Inf return else - p = convert(Int, round(abs(log(10, Alp.get_option(m, :rel_gap))))) + p = convert(Int, round(abs(log(10, get_option(m, :rel_gap))))) n = round(abs(m.best_obj - m.best_bound); digits = p) # dn = round(abs(1e-12+abs(m.best_obj)); digits=p) if isapprox(n, 0.0; atol = tol) && isapprox(m.best_obj, 0.0; atol = tol) @@ -29,10 +29,10 @@ function update_opt_gap(m::Optimizer) return end - if Alp.is_min_sense(m) - m.best_rel_gap = Alp.eval_opt_gap(m, m.best_bound, m.best_obj) - elseif Alp.is_max_sense(m) - m.best_rel_gap = Alp.eval_opt_gap(m, m.best_obj, m.best_bound) + if is_min_sense(m) + m.best_rel_gap = eval_opt_gap(m, m.best_bound, m.best_obj) + elseif is_max_sense(m) + m.best_rel_gap = eval_opt_gap(m, m.best_obj, m.best_bound) end end @@ -42,16 +42,16 @@ function update_opt_gap(m::Optimizer) end function eval_opt_gap(m::Optimizer, lower_bound::Number, upper_bound::Number) - tol = Alp.get_option(m, :tol) + tol = get_option(m, :tol) if isinf(lower_bound) || isinf(upper_bound) error("Infinite bounds detected in optimality gap evalutation") else - m.best_rel_gap = Alp._eval_best_rel_gap(lower_bound, upper_bound, tol) + m.best_rel_gap = _eval_best_rel_gap(lower_bound, upper_bound, tol) # If LB > UB - if (m.best_rel_gap > Alp.get_option(m, :rel_gap)) && + if (m.best_rel_gap > get_option(m, :rel_gap)) && ((lower_bound - upper_bound) > 1E-1) error( "Lower bound cannot exceed the upper bound in optimality gap evaluation", @@ -95,7 +95,7 @@ function measure_relaxed_deviation(m::Optimizer; sol = nothing) sort!(dev, by = x -> x[1]) for i in dev - Alp.get_option(m, :log_level) > 199 && println( + get_option(m, :log_level) > 199 && println( "Y-VAR$(i[1]): DIST=$(i[2]) || Y-hat = $(i[3]), Y-val = $(i[4]) || COMP $(i[5])", ) end @@ -108,7 +108,7 @@ end Same as [`update_var_bounds`](@ref) """ -discretization_to_bounds(d::Dict, l::Int) = Alp.update_var_bounds(d, len = l) +discretization_to_bounds(d::Dict, l::Int) = update_var_bounds(d, len = l) """ Update the data structure with feasible solution and its associated objective (if better) @@ -178,24 +178,24 @@ end Check if the solution is always the same within the last disc_consecutive_forbid iterations. Return `true` if solution has stalled. """ function check_solution_history(m::Optimizer, ind::Int) - Alp.get_option(m, :disc_consecutive_forbid) == 0 && return false - ((m.logs[:n_iter] - 1) < Alp.get_option(m, :disc_consecutive_forbid)) && return false + get_option(m, :disc_consecutive_forbid) == 0 && return false + ((m.logs[:n_iter] - 1) < get_option(m, :disc_consecutive_forbid)) && return false sol_val = m.bound_sol_history[mod( m.logs[:n_iter] - 2, - Alp.get_option(m, :disc_consecutive_forbid), + get_option(m, :disc_consecutive_forbid), )+1][ind] - for i in 1:(Alp.get_option(m, :disc_consecutive_forbid)-1) + for i in 1:(get_option(m, :disc_consecutive_forbid)-1) search_pos = - mod(m.logs[:n_iter] - 2 - i, Alp.get_option(m, :disc_consecutive_forbid)) + 1 + mod(m.logs[:n_iter] - 2 - i, get_option(m, :disc_consecutive_forbid)) + 1 !isapprox( sol_val, m.bound_sol_history[search_pos][ind]; - atol = Alp.get_option(m, :disc_rel_width_tol), + atol = get_option(m, :disc_rel_width_tol), ) && return false end - Alp.get_option(m, :log_level) > 99 && + get_option(m, :log_level) > 99 && println("Consecutive bounding solution on VAR$(ind) obtained. Diverting...") return true @@ -223,8 +223,8 @@ function fix_domains(m::Optimizer; discrete_sol = nothing, use_orig = false) end PCnt = length(m.discretization[i]) - 1 for j in 1:PCnt - if point >= (m.discretization[i][j] - Alp.get_option(m, :tol)) && - (point <= m.discretization[i][j+1] + Alp.get_option(m, :tol)) + if point >= (m.discretization[i][j] - get_option(m, :tol)) && + (point <= m.discretization[i][j+1] + get_option(m, :tol)) @assert j < length(m.discretization[i]) use_orig ? l_var[i] = m.discretization[i][1] : l_var[i] = m.discretization[i][j] @@ -376,8 +376,8 @@ end function build_discvar_graph(m::Optimizer) # Collect the information of nonlinear terms in terms of arcs and nodes - nodes = Alp.get_candidate_disc_vars(m, getoutput = true) - arcs = Alp.ncvar_collect_arcs(m, nodes) + nodes = get_candidate_disc_vars(m, getoutput = true) + arcs = ncvar_collect_arcs(m, nodes) # Collect integer variables for i in 1:m.num_var_orig @@ -401,10 +401,10 @@ nonlinear terms which are adaptively partitioned for global optimization. This o setting `disc_var_pick = 1`. """ function min_vertex_cover(m::Optimizer) - nodes, arcs = Alp.build_discvar_graph(m) + nodes, arcs = build_discvar_graph(m) # Set up minimum vertex cover problem - minvertex = Model(Alp.get_option(m, :mip_solver)) + minvertex = Model(get_option(m, :mip_solver)) MOI.set(minvertex, MOI.TimeLimitSec(), 60.0) # Time limit for min vertex cover formulation JuMP.@variable(minvertex, 0 <= x[nodes] <= 1, Bin) JuMP.@constraint(minvertex, [a in arcs], x[a[1]] + x[a[2]] >= 1) @@ -417,8 +417,8 @@ function min_vertex_cover(m::Optimizer) # Collecting required information m.num_var_disc_mip = Int(sum(xVal)) m.disc_vars = [ - i for i in nodes if xVal[i] > Alp.get_option(m, :tol) && - abs(m.u_var_tight[i] - m.l_var_tight[i]) >= Alp.get_option(m, :tol) + i for i in nodes if xVal[i] > get_option(m, :tol) && + abs(m.u_var_tight[i] - m.l_var_tight[i]) >= get_option(m, :tol) ] return @@ -427,7 +427,7 @@ end function weighted_min_vertex_cover(m::Optimizer, distance::Dict) # Collect the graph information - nodes, arcs = Alp.build_discvar_graph(m) + nodes, arcs = build_discvar_graph(m) # A little bit redundancy before disvec = [distance[i] for i in keys(distance) if i in m.candidate_disc_vars] @@ -437,12 +437,12 @@ function weighted_min_vertex_cover(m::Optimizer, distance::Dict) for i in m.candidate_disc_vars isapprox(distance[i], 0.0; atol = 1e-6) ? weights[i] = heavy : (weights[i] = (1 / distance[i])) - (Alp.get_option(m, :log_level) > 100) && + (get_option(m, :log_level) > 100) && println("VAR$(i) WEIGHT -> $(weights[i]) ||| DISTANCE -> $(distance[i])") end # Set up minimum vertex cover problem - minvertex = Model(Alp.get_option(m, :mip_solver)) + minvertex = Model(get_option(m, :mip_solver)) MOI.set(minvertex, MOI.TimeLimitSec(), 60.0) # Set a timer to avoid waste of time in proving optimality JuMP.@variable(minvertex, 0 <= x[nodes] <= 1, Bin) for arc in arcs @@ -456,29 +456,29 @@ function weighted_min_vertex_cover(m::Optimizer, distance::Dict) xVal = JuMP.value.(x) m.num_var_disc_mip = Int(sum(xVal)) m.disc_vars = [ - i for i in nodes if xVal[i] > 0 && - abs(m.u_var_tight[i] - m.l_var_tight[i]) >= Alp.get_option(m, :tol) + i for i in nodes if + xVal[i] > 0 && abs(m.u_var_tight[i] - m.l_var_tight[i]) >= get_option(m, :tol) ] - Alp.get_option(m, :log_level) >= 99 && + get_option(m, :log_level) >= 99 && println("UPDATED DISC-VAR COUNT = $(length(m.disc_vars)) : $(m.disc_vars)") return end function _fetch_mip_solver_identifier(m::Optimizer; override = "") - (Alp.get_option(m, :mip_solver) === nothing) && return + (get_option(m, :mip_solver) === nothing) && return m.mip_solver_id = _get_solver_name(m, :mip_solver, override) return end function _fetch_nlp_solver_identifier(m::Optimizer; override = "") - (Alp.get_option(m, :nlp_solver) === nothing) && return + (get_option(m, :nlp_solver) === nothing) && return m.nlp_solver_id = _get_solver_name(m, :nlp_solver, override) return end function _fetch_minlp_solver_identifier(m::Optimizer; override = "") - (Alp.get_option(m, :minlp_solver) === nothing) && return + (get_option(m, :minlp_solver) === nothing) && return m.minlp_solver_id = _get_solver_name(m, :minlp_solver, override) return end @@ -486,7 +486,7 @@ end function _get_solver_name(m::Optimizer, key::Symbol, override::String) solver_string = override if isempty(override) - solver = MOI.instantiate(Alp.get_option(m, key)) + solver = MOI.instantiate(get_option(m, key)) solver_string = try MOI.get(solver, MOI.SolverName()) catch @@ -507,7 +507,7 @@ end An utility function used to dynamically regulate MILP solver time limits to fit Alpine solver's time limits. """ function set_mip_time_limit(m::Optimizer) - time_limit = max(0.0, Alp.get_option(m, :time_limit) - m.logs[:total_time]) + time_limit = max(0.0, get_option(m, :time_limit) - m.logs[:total_time]) return MOI.set(m.model_mip, MOI.TimeLimitSec(), time_limit) end @@ -535,7 +535,7 @@ function resolve_lifted_var_value(m::Optimizer, sol_vec::Array) return sol_vec end -# Unused functions +# Unused functions # function amp_post_λ_upperbound( # m::Optimizer, # λ::Dict, @@ -564,7 +564,7 @@ end # for i in 1:length(tregions[level+1]) # push!(reg, tregions[level+1][i]) -# Alp.amp_post_λ_upperbound(m, λ, indices, dim, d, tregions, reg, level + 1) +# amp_post_λ_upperbound(m, λ, indices, dim, d, tregions, reg, level + 1) # length(reg) < level && error("Something is wrong") # length(reg) > level && pop!(reg) # end @@ -593,7 +593,7 @@ end # sum(α[v][m.bound_sol_pool[:disc][i][v]] for v in no_good_idxs) <= # no_good_size # ) -# Alp.get_option(m, :log_level) > 0 && println( +# get_option(m, :log_level) > 0 && println( # "!! GLOBAL cuts off POOL_SOL-$(i) POOL_OBJ=$(m.bound_sol_pool[:obj][i])!", # ) # m.bound_sol_pool[:stat][i] = :Cutoff diff --git a/src/variable_bounds.jl b/src/variable_bounds.jl index 20c8bbd4..2766108d 100644 --- a/src/variable_bounds.jl +++ b/src/variable_bounds.jl @@ -104,7 +104,7 @@ x >= 5, x <= 5 or x == 5 and fetch the information to m.l_var_tight and m.u_var_ function bound_propagation(m::Optimizer) exhausted = false infeasible = false - tol = Alp.get_option(m, :tol) + tol = get_option(m, :tol) while !exhausted exhausted = true @@ -146,11 +146,11 @@ function bound_propagation(m::Optimizer) if eval_l_bound > m.l_var_tight[var_idx] + tol exhausted = false m.l_var_tight[var_idx] = eval_l_bound - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] LB $(m.l_var_tight[var_idx]) evaluated from constraint", ) elseif eval_l_bound > m.u_var_tight[var_idx] + tol - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] Infeasibility detection during bound propagation", ) infeasible = true @@ -160,11 +160,11 @@ function bound_propagation(m::Optimizer) if eval_u_bound < m.u_var_tight[var_idx] - tol exhausted = false m.u_var_tight[var_idx] = eval_u_bound - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] UB $(m.u_var_tight[var_idx]) evaluated from constraints", ) elseif eval_u_bound < m.l_var_tight[var_idx] - tol - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] Infeasibility detection during bound propagation", ) infeasible = true @@ -190,11 +190,11 @@ function bound_propagation(m::Optimizer) if eval_bound > m.l_var_tight[var_idx] + tol exhausted = false m.l_var_tight[var_idx] = eval_bound - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] LB $(m.l_var_tight[var_idx]) evaluated from constraints", ) elseif eval_bound > m.u_var_tight[var_idx] + tol - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] Infeasibility detection during bound propagation", ) infeasible = true @@ -219,11 +219,11 @@ function bound_propagation(m::Optimizer) if eval_bound < m.u_var_tight[var_idx] - tol exhausted = false m.u_var_tight[var_idx] = eval_bound - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] UB $(m.u_var_tight[var_idx]) evaluated from constraints", ) elseif eval_bound < m.l_var_tight[var_idx] - tol - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] Infeasibility detection during bound propagation", ) infeasible = true @@ -248,11 +248,11 @@ function bound_propagation(m::Optimizer) if eval_bound < m.u_var_tight[var_idx] - tol exhausted = false m.u_var_tight[var_idx] = eval_bound - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] UB $(m.u_var_tight[var_idx]) evaluated from constraints", ) elseif eval_bound < m.l_var_tight[var_idx] - tol - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] Infeasibility detection during bound propagation", ) infeasible = true @@ -276,11 +276,11 @@ function bound_propagation(m::Optimizer) if eval_bound > m.l_var_tight[var_idx] + tol exhausted = false m.l_var_tight[var_idx] = eval_bound - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] LB $(m.l_var_tight[var_idx]) evaluated from constraints", ) elseif eval_bound > m.u_var_tight[var_idx] + tol - (Alp.get_option(m, :log_level) > 199) && println( + (get_option(m, :log_level) > 199) && println( "[VAR$(var_idx)] Infeasibility detection during bound propagation", ) infeasible = true @@ -289,7 +289,7 @@ function bound_propagation(m::Optimizer) end end end - (exhausted == true && Alp.get_option(m, :log_level) > 99) && + (exhausted == true && get_option(m, :log_level) > 99) && println("Initial constraint-based bound evaluation exhausted...") end @@ -324,12 +324,12 @@ in known or trivial bounds information to deduce lifted variable bounds and to p function resolve_var_bounds(m::Optimizer) # Basic Bound propagation - if Alp.get_option(m, :presolve_bp) - setproperty!(m, :presolve_infeasible, Alp.bound_propagation(m)) # Fetch bounds from constraints + if get_option(m, :presolve_bp) + setproperty!(m, :presolve_infeasible, bound_propagation(m)) # Fetch bounds from constraints end - # Resolve unbounded variables in the original formulation - Alp.resolve_inf_bounds(m) + # Resolve unbounded variables in the original formulation + resolve_inf_bounds(m) # Added sequential bound resolving process base on DFS process, which ensures all bounds are secured. # Increased complexity from linear to square but a reasonable amount @@ -339,7 +339,7 @@ function resolve_var_bounds(m::Optimizer) if haskey(m.nonconvex_terms, k) m.nonconvex_terms[k][:bound_resolver](m, k) elseif haskey(m.linear_terms, k) - Alp.basic_linear_bounds(m, k) + basic_linear_bounds(m, k) else error("Found homeless term key $(k) during bound resolution.") end @@ -360,23 +360,23 @@ function resolve_inf_bounds(m::Optimizer) for i in 1:length(m.l_var_orig) if m.l_var_tight[i] == -Inf warnuser = true - m.l_var_tight[i] = -Alp.get_option(m, :large_bound) + m.l_var_tight[i] = -get_option(m, :large_bound) infcount_l += 1 end if m.u_var_tight[i] == Inf warnuser = true - m.u_var_tight[i] = Alp.get_option(m, :large_bound) + m.u_var_tight[i] = get_option(m, :large_bound) infcount_u += 1 end end infcount = min(infcount_l, infcount_u) if infcount == 1 warnuser && println( - "Warning: -/+Inf bounds detected on at least $infcount variable. Initializing with values -/+$(Alp.get_option(m, :large_bound)). This may affect global optimal values and run times.", + "Warning: -/+Inf bounds detected on at least $infcount variable. Initializing with values -/+$(get_option(m, :large_bound)). This may affect global optimal values and run times.", ) elseif infcount > 1 warnuser && println( - "Warning: -/+Inf bounds detected on at least $infcount variables. Initializing with values -/+$(Alp.get_option(m, :large_bound)). This may affect global optimal values and run times.", + "Warning: -/+Inf bounds detected on at least $infcount variables. Initializing with values -/+$(get_option(m, :large_bound)). This may affect global optimal values and run times.", ) end @@ -400,16 +400,16 @@ function resolve_var_bounds(m::Optimizer, d::Dict; kwargs...) if haskey(m.nonconvex_terms, k) nlk = k if m.nonconvex_terms[nlk][:nonlinear_type] in ALPINE_C_MONOMIAL - d = Alp.basic_monomial_bounds(m, nlk, d) + d = basic_monomial_bounds(m, nlk, d) elseif m.nonconvex_terms[nlk][:nonlinear_type] in [:BINPROD] - d = Alp.basic_binprod_bounds(m, nlk, d) + d = basic_binprod_bounds(m, nlk, d) elseif m.nonconvex_terms[nlk][:nonlinear_type] in [:BINLIN] - d = Alp.basic_binlin_bounds(m, nlk, d) + d = basic_binlin_bounds(m, nlk, d) else error("EXPECTED ERROR : NEED IMPLEMENTATION") end elseif haskey(m.linear_terms, k) - d = Alp.basic_linear_bounds(m, k, d) + d = basic_linear_bounds(m, k, d) else error( "Found a homeless term key $(k) during bound resolution im resolve_var_bounds.", @@ -457,7 +457,7 @@ end # function resolve_closed_var_bounds(m::Optimizer; kwargs...) # for var in m.candidate_disc_vars # if abs(m.l_var_tight[var] - m.u_var_tight[var]) < -# Alp.get_option(m, :presolve_bt_width_tol) # Closed Bound Criteria +# get_option(m, :presolve_bt_width_tol) # Closed Bound Criteria # deleteat!(m.disc_vars, findfirst(m.disc_vars, var)) # Clean nonconvex_terms by deleting the info # m.discretization[var] = [m.l_var_tight[var], m.u_var_tight[var]] # Clean up the discretization for basic McCormick if necessary # end From 18a4adb4def2f12b1057ffc9f79df0a8ad09affe Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 12 Aug 2024 14:47:10 +1200 Subject: [PATCH 2/2] Update src/log.jl --- src/log.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/log.jl b/src/log.jl index c4473e21..6bf1c4ee 100644 --- a/src/log.jl +++ b/src/log.jl @@ -67,7 +67,7 @@ function logging_summary(m::Optimizer) project = read(joinpath(dirname(@__DIR__), "Project.toml"), String) m_version = match(r"version \= \"(.+?)\"", project) println(" Alpine version = ", m_version[1]) - if Alp.is_min_sense(m) + if is_min_sense(m) println( " Maximum iterations (lower-bounding MIPs) = ", get_option(m, :max_iter),