Merge pull request #16 from MatthewStuber/InitCommit

Init commit

Author: MatthewStuber

README.md

@@ -1,2 +1,50 @@
-# EAGODomainReduction
-A Package Containing Multiple Domain Reduction Procedures Used in Global Optimization
+# EAGODomainReduction.jl
+Domain Reduction Procedures in Global Optimization
+
+## Installation
+
+```julia
+julia> Pkg.add("EAGODomainReduction.jl")
+```
+
+## Capabilities
+
+**EAGODomainReduction.jl** provides a series of subroutines for tightening the domains of subproblems
+solved in global optimization (potentially to in-feasibility). Currently, it supports routines for
+nonconvex nonlinear programs:
+- **Interval Contractor Propagation:** A forward-backward interval contractor using **IntervalArthimetic.jl**
+and **IntervalContractors.jl** for the operator library.
+- **Duality-Based Bound Tightening:** Provides algorithms for tightening domains based on duality of solutions
+found for subproblems.
+- **Standard Range Reduction:** Contracts subproblem domain via linear-relaxations generated using McCormick
+relaxations.
+- **Implicit Subroutine support:** Supports domain reduction of reduced space lower-bound problems defined through
+relaxation of implicit functions by fixed-point methods.
+
+The routine are used extensively in the `EAGO.jl`[https://github.com/MatthewStuber/EAGO.jl] package solver.
+Please see the example files for usage cases.
+
+## Future Work
+
+- Update the interval-constraint propagation algorithm to incorporate propagation heuristics in Vu2008.
+- Incorporate control-flow syntax support into constraint propagation algorithm.
+- Incorporate improvements to probing and optimality based-bound tightening.
+- Add support for Mixed-Integer NLP.
+
+## Related Packages
+- **EAGO.jl**[https://github.com/MatthewStuber/EAGO.jl]: A package containing global and robust solvers based mainly on McCormick relaxations.
+This package supports a JuMP and MathProgBase interface.
+- **IntervalConstraintProgramming.jl**[https://github.com/JuliaIntervals/IntervalConstraintProgramming.jl]: Provides algorithms that furnish bounds
+on constraints defined by expressions. The constraint propagation routine in **EAGODomainReduction.jl** can generate tape objects that are
+reusable for generically-defined functions. In addition, we use a `Vector{Interval}` storage object that allows for in-place mutation of intervals.
+- **IntervalContractors.jl**[https://github.com/JuliaIntervals/IntervalContractors.jl]: Provides a library of reverse interval contractors.
+
+## References
+- Benhamou, F., & Older, W.J. (1997). Applying interval arithmetic to real, integer, and boolean constraints. The Journal of Logic Programming, 32, 1–24.
+- Caprara, A., & Locatelli, M. (2010). Global optimization problems and domain reduction strategies. Mathematical Programming, 125, 123–137.
+- Gleixner, A.M., Berthold, T., Müller, B., & Weltge, S. (2016). Three enhancements for optimization-based bound tightening. ZIB Report, 15–16.
+- Ryoo, H.S., & Sahinidis, N.V. (1996). A branch-and-reduce approach to global optimization. Journal of Global Optimization, 8, 107–139.
+- Schichl, H., & Neumaier, A. (2005). Interval analysis on directed acyclic graphs for global optimization. Journal of Global Optimization, 33, 541–562.
+- Tawarmalani, M., & Sahinidis, N.V. (2005). A polyhedral branch-and-cut approach to global optimization. Mathematical Programming, 103, 225–249.
+- Vu, X., Schichl, H., & Sam-Haroud, D. (2009). Interval propagation and search on directed acyclic
+graphs for numerical constraint solving. Journal of Global Optimization, 45, 499–531.

src/src/FBBT/DAGprop.jl

@@ -272,6 +272,18 @@ function ForwardPass!(x::Tape)
   end
 end
 
+"""
+    ForwardPassOne!(x::Tape)
+
+Performs a forward contractor pass using tape `x::Tape` with some uninitialized nodes.
+"""
+function ForwardPassOne!(x::Tape)
+  for i=1:length(x.Edge_List)
+    x.FW_Expr[i].args = vcat([x.Head_List[i]],x.Intv_Storage[x.FW_Arg[i]])
+    x.Intv_Storage[x.Edge_List[i][1][2]] = eval(x.FW_Expr[i])
+  end
+end
+
 """
     ReversePass!(x::Tape)
 
@@ -293,12 +305,15 @@ and the initial interval bounds `X::Vector{Interval{T}}`.
 function DAGContractor!(X::Vector{Interval{T}},x::TapeList,r) where {T}
   Xprev::Vector{Interval} = copy(X) # sets variable bounds on first Array to Box Bounds
   SetConstraintNode!(x)
-  #SetConstantNode!(x)
   for i=1:r
     for j=1:length(x.sto)
       SetVarBounds!(x.sto[j],Xprev) # take refined bounds from previous graph
       SetConstantNode!(x)
-      ForwardPass!(x.sto[j]) # run forward and reverse pass
+      if (i == 1)
+        ForwardPassOne!(x.sto[j])
+      else
+        ForwardPass!(x.sto[j]) # run forward and reverse pass
+      end
       ReversePass!(x.sto[j])
       Xprev = x.sto[j].Intv_Storage[1:length(X)]
     end

src/src/FBBT/NodeFinder.jl

@@ -8,6 +8,8 @@ immutable NodeFinder
   ind::Int64
 end
 
+step(x::T) where T = x>0 ? one(T) : zero(T)
+
 # Defines conversion for NodeFinder Type
 for flt in (Float16,Float32,Float64,Int8,Int16,Int32,Int64,Int128,
             UInt8,UInt16,UInt32,UInt64,UInt128)

src/src/OBBT/STD_LP_Probe.jl

@@ -12,6 +12,9 @@ that indicates the upper bound is finite. The upper bound is `UBD::Float64`.
 """
 function STD_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)
 
+    mu_temp::Int64 = copy(EAGOSmoothMcCormickGrad.MC_param.mu)
+    set_diff_relax(0)
+
     # constructs LP relaxation
     Xlo::Vector{Float64} = [X[i].lo for i=1:opt[1].numVar]
     Xhi::Vector{Float64} = [X[i].hi for i=1:opt[1].numVar]
@@ -27,39 +30,52 @@ function STD_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)
     # probes upper bound
     f_mc::SMCg{opt[1].numVar,Float64} = opt[1].f(x_SMC)
     f_cv::Float64 = f_mc.cv
-    c::Vector{SMCg{opt[1].numVar,Float64}} = opt[1].g(x_SMC)
-    # forms coefficients of linear relaxations
-    dcdx::SparseMatrixCSC{Float64,Int64} = spzeros(length(opt[1].gL_loc)+length(opt[1].gU_loc)+1,opt[1].numVar)
+    if opt[1].numConstr>0
+        c::Vector{SMCg{opt[1].numVar,Float64}} = opt[1].g(x_SMC)
+        dcdx::SparseMatrixCSC{Float64,Int64} = spzeros(length(opt[1].gL_loc)+length(opt[1].gU_loc)+1,opt[1].numVar)
+    else
+        dcdx = spzeros(1,opt[1].numVar)
+    end
+
     dcdx[1,:] = f_mc.cv_grad
-    cx_ind1::Int64 = 2
-    for i in opt[1].gL_loc
-        for j=1:opt[1].numVar
-            if (c[i].cv_grad[j] != 0.0)
-                dcdx[cx_ind1,j] = c[i].cv_grad[j]
+    if opt[1].numConstr>0
+        cx_ind1::Int64 = 2
+        for i in opt[1].gL_loc
+            for j=1:opt[1].numVar
+                if (c[i].cv_grad[j] != 0.0)
+                    dcdx[cx_ind1,j] = c[i].cv_grad[j]
+                end
             end
+            cx_ind1 += 1
         end
-        cx_ind1 += 1
-    end
-    for i in opt[1].gU_loc
-        for j=1:opt[1].numVar
-            if (c[i].cc_grad[j] != 0.0)
-                dcdx[cx_ind1,j] = -c[i].cc_grad[j]
+        for i in opt[1].gU_loc
+            for j=1:opt[1].numVar
+                if (c[i].cc_grad[j] != 0.0)
+                    dcdx[cx_ind1,j] = -c[i].cc_grad[j]
+                end
             end
+            cx_ind1 += 1
         end
-        cx_ind1 += 1
     end
 
     # forms rhs of linear relaxations
-    rhs::Vector{Float64} = zeros(Float64,length(opt[1].gL_loc)+length(opt[1].gU_loc)+1)
-    rhs[1] = UBD + f_mc.cv - sum(x0.*f_mc.cv_grad)
-    cx_ind2::Int64 = 2
-    for i in opt[1].gU_loc
-        rhs[cx_ind2] = sum(x0[:].*c[i].cv_grad[:])+opt[1].gU[i]-c[i].cv
-        cx_ind2 += 1
+    if opt[1].numConstr>0
+        rhs::Vector{Float64} = zeros(Float64,length(opt[1].gL_loc)+length(opt[1].gU_loc)+1)
+    else
+        rhs = zeros(Float64,1)
     end
-    for i in opt[1].gL_loc
-        rhs[cx_ind2] = sum(-x0[:].*c[i].cc_grad[:])-opt[1].gL[i]+c[i].cc
-        cx_ind2 += 1
+    rhs[1] = UBD + f_mc.cv - sum(x0.*f_mc.cv_grad)
+
+    if opt[1].numConstr>0
+        cx_ind2::Int64 = 2
+        for i in opt[1].gU_loc
+            rhs[cx_ind2] = sum(x0[:].*c[i].cv_grad[:])+opt[1].gU[i]-c[i].cv
+            cx_ind2 += 1
+        end
+        for i in opt[1].gL_loc
+            rhs[cx_ind2] = sum(-x0[:].*c[i].cc_grad[:])-opt[1].gL[i]+c[i].cc
+            cx_ind2 += 1
+        end
     end
 
     if (opt[1].numConstr>0)
@@ -78,8 +94,8 @@ function STD_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)
             val::Float64 = result.objval + f_cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])
             pnt::Vector{Float64} = result.sol
             mult::Vector{Float64} = result.attrs[:redcost]
-            mult_lo::Vector{Float64} = [tol_eq(X[i].lo,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
-            mult_hi::Vector{Float64} = [tol_eq(X[i].hi,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
+            mult_lo::Vector{Float64} = [tol_eq(X[i].lo,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
+            mult_hi::Vector{Float64} = [tol_eq(X[i].hi,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
             Variable_DR!(X,mult_lo,mult_hi,val,UBD)
         end
         # probes lower bound
@@ -90,11 +106,18 @@ function STD_LP_Probe!(X::Vector{Interval{Float64}},opt,UBD::Float64)
             val = result.objval + f_cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])
             pnt = result.sol
             mult = result.attrs[:redcost]
-            mult_lo = [tol_eq(X[i].lo,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
-            mult_hi = [tol_eq(X[i].hi,pnt[i],1E-12) ? mult[i] : 0.0 for i=1:opt[1].numVar]
+            mult_lo = [tol_eq(X[i].lo,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
+            mult_hi = [tol_eq(X[i].hi,pnt[i],opt[1].solver.dual_tol) ? mult[i] : 0.0 for i=1:opt[1].numVar]
             Variable_DR!(X,mult_lo,mult_hi,val,UBD)
         end
     end
+
+    # reset McCormick relaxation to original value
+    set_diff_relax(mu_temp)
+
+    # stores outputs
+    X[:] = X
+    return true
 end
 
 """