1 /* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
2 
3 /*
4  *  Main authors:
5  *     Gleb Belov <gleb.belov@monash.edu>
6  */
7 
8 /* This Source Code Form is subject to the terms of the Mozilla Public
9  * License, v. 2.0. If a copy of the MPL was not distributed with this
10  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
11 
12 #pragma once
13 
14 #include <utility>
15 #include <vector>
16 //#include <map>
17 #include <minizinc/solver_instance_defs.hh>
18 
19 #include <cassert>
20 #include <chrono>
21 #include <iostream>
22 #include <sstream>
23 #include <string>
24 #include <unordered_map>
25 
26 /// Facilitate lhs computation of a cut
compute_sparse(int n,const int * ind,const double * coef,const double * dense,int nVarsDense)27 inline double compute_sparse(int n, const int* ind, const double* coef, const double* dense,
28                              int nVarsDense) {
29   assert(ind && coef && dense);
30   double val = 0.0;
31   for (int i = 0; i < n; ++i) {
32     assert(ind[i] >= 0);
33     assert(ind[i] < nVarsDense);
34     val += coef[i] * dense[ind[i]];
35   }
36   return val;
37 }
38 
39 class MIPWrapper;
40 
41 /// An abstract MIP wrapper.
42 /// Does not include MZN stuff so can be used independently
43 /// although it's limited to the MZN solving needs.
44 class MIPWrapper {
45 public:
46   typedef int VarId;  // CPLEX uses int
47   enum VarType { REAL, INT, BINARY };
48   enum LinConType { LQ = -1, EQ = 0, GQ = 1 };
49 
50   // CPLEX 12.6.2 advises anti-symmetry constraints to be user+lazy
51   static const int MaskConsType_Normal = 1;
52   /// User cut. Only cuts off fractional points, no integer feasible points
53   static const int MaskConsType_Usercut = 2;
54   /// Lazy cut. Can cut off otherwise feasible integer solutions.
55   /// Callback should be able to produce previously generated cuts again if needed [Gurobi]
56   static const int MaskConsType_Lazy = 4;
57   enum Status { OPT, SAT, UNSAT, UNBND, UNSATorUNBND, UNKNOWN, ERROR_STATUS };
58 
59   /// Search strategy for the solver
60   enum SearchType { FIXED_SEARCH = 0, FREE_SEARCH = 1, UNIFORM_SEARCH = 2 };
61 
62   /// Columns for SCIP upfront and with obj coefs:
63   std::vector<double> colObj, colLB, colUB;
64   std::vector<VarType> colTypes;
65   std::vector<std::string> colNames;
66   //     , rowLB, rowUB, elements;
67   //     veci whichInt
68   //     , starts, column;
69   //     double objUB;
70   //     double qpu;
71 
72   /// Parameter
73   bool fVerbose = false;
74 
75   int nProbType = -2;  // +-1: max/min; 0: sat
76 
77   struct Output {
78     Status status;
79     std::string statusName = "Untouched";
80     double objVal = 1e308;
81     double bestBound = 1e308;
82     int nCols = 0;
83     int nObjVarIndex = -1;
84     const double* x = nullptr;
85     int nNodes = 0;
86     int nOpenNodes = 0;
87     double dWallTime = 0.0;
88     std::chrono::time_point<std::chrono::steady_clock> dWallTime0;
89     double dCPUTime = 0;
90     std::clock_t cCPUTime0 = 0;
91   };
92   Output output;
93 
94   /// General cut definition, could be used for addRow() too
95   class CutDef {
CutDef()96     CutDef() {}
97 
98   public:
CutDef(LinConType s,int m)99     CutDef(LinConType s, int m) : sense(s), mask(m) {}
100     std::vector<int> rmatind;
101     std::vector<double> rmatval;
102     LinConType sense = LQ;
103     double rhs = 0.0;
104     int mask = 0;  // need to know what type of cuts are registered before solve()  TODO
105     std::string rowName;
addVar(int i,double c)106     void addVar(int i, double c) {
107       rmatind.push_back(i);
108       rmatval.push_back(c);
109     }
computeViol(const double * x,int nCols)110     double computeViol(const double* x, int nCols) {
111       double lhs = compute_sparse(static_cast<int>(rmatind.size()), rmatind.data(), rmatval.data(),
112                                   x, nCols);
113       if (LQ == sense) {
114         return lhs - rhs;
115       }
116       if (GQ == sense) {
117         return rhs - lhs;
118       }
119       assert(0);
120 
121       return 0.0;
122     }
123   };
124   /// Cut callback fills one
125   typedef std::vector<CutDef> CutInput;
126 
127   /// solution callback handler, the wrapper might not have these callbacks implemented
128   typedef void (*SolCallbackFn)(const Output&, void*);
129   /// cut callback handler, the wrapper might not have these callbacks implemented
130   typedef void (*CutCallbackFn)(const Output&, CutInput&, void*,
131                                 bool fMIPSol  // if with a MIP feas sol - lazy cuts only
132   );
133   struct CBUserInfo {
134     MIPWrapper* wrapper = nullptr;
135     MIPWrapper::Output* pOutput = nullptr;
136     MIPWrapper::Output* pCutOutput = nullptr;
137     void* psi = nullptr;  // external info. Intended to keep MIPSolverinstance
138     SolCallbackFn solcbfn = nullptr;
139     CutCallbackFn cutcbfn = nullptr;
140     /// Union of all flags used for the registered callback cuts
141     /// See MaskConstrType_..
142     /// Solvers need to know this
143     /// In MIPSolverinstance, class CutGen defines getMask() which should return that
144     int cutMask = 0;             // can be any combination of User/Lazy
145     bool fVerb = false;          // used in Gurobi
146     bool printed = false;        // whether any solution was output
147     double nTimeoutFeas = -1.0;  // >=0 => stop that long after 1st feas
148     double nTime1Feas = -1e100;  // time of the 1st feas
149   };
150   CBUserInfo cbui;
151 
MIPWrapper()152   MIPWrapper() { cbui.wrapper = this; }
~MIPWrapper()153   virtual ~MIPWrapper() { /* cleanup(); */
154   }
155 
156   /// derived should overload and call the ancestor
157   //     virtual void cleanup() {
158   //       colObj.clear(); colLB.clear(); colUB.clear();
159   //       colTypes.clear(); colNames.clear();
160   //     }
161   /// re-create solver object. Called from the base class constructor
162   //     virtual void resetModel() { };
163 
164   //     virtual void printVersion(ostream& os) { os << "Abstract MIP wrapper"; }
165   //     virtual void printHelp(ostream& ) { }
166 
167   bool fPhase1Over = false;
168 
169 private:
170   /// adding a variable just internally (in Phase 1 only that). Not to be used directly.
addVarLocal(double obj,double lb,double ub,VarType vt,const std::string & name="")171   virtual VarId addVarLocal(double obj, double lb, double ub, VarType vt,
172                             const std::string& name = "") {
173     //       cerr << "  addVarLocal: colObj.size() == " << colObj.size()
174     //         << " obj == " <<obj
175     //         << " lb == " << lb
176     //         << " ub == " << ub
177     //         << " vt == " << vt
178     //         << " nm == " << name
179     //         << endl;
180     colObj.push_back(obj);
181     colLB.push_back(lb);
182     colUB.push_back(ub);
183     colTypes.push_back(vt);
184     colNames.push_back(name);
185     return static_cast<VarId>(colObj.size() - 1);
186   }
187   /// add the given var to the solver. Asserts all previous are added. Phase >=2. No direct use
addVar(int j)188   virtual void addVar(int j) {
189     assert(j == getNCols());
190     assert(fPhase1Over);
191     doAddVars(1, &colObj[j], &colLB[j], &colUB[j], &colTypes[j], &colNames[j]);
192   }
193   /// actual adding new variables to the solver. "Updates" the model (e.g., Gurobi). No direct use
194   virtual void doAddVars(size_t n, double* obj, double* lb, double* ub, VarType* vt,
195                          std::string* names) = 0;
196 
197 public:
198   /// debugging stuff
199   //     set<double> sLitValues;
200   std::unordered_map<double, VarId> sLitValues;
201 
setProbType(int t)202   void setProbType(int t) { nProbType = t; }
203 
204   /// adding a variable, at once to the solver, this is for the 2nd phase
addVar(double obj,double lb,double ub,VarType vt,const std::string & name="")205   virtual VarId addVar(double obj, double lb, double ub, VarType vt, const std::string& name = "") {
206     //       cerr << "  AddVar: " << lb << ":   ";
207     VarId res = addVarLocal(obj, lb, ub, vt, name);
208     if (fPhase1Over) {
209       addVar(res);
210     }
211     return res;
212   }
213   int nLitVars = 0;
214   /// adding a literal as a variable. Should not happen in feasible models
addLitVar(double v)215   virtual VarId addLitVar(double v) {
216     // Cannot do this: at least CBC does not support duplicated indexes    TODO??
217     //       ++nLitVars;
218     //       auto itFound = sLitValues.find(v);
219     //       if (sLitValues.end() != itFound)
220     //         return itFound->second;
221     std::ostringstream oss;
222     oss << "lit_" << v << "__" << (nLitVars++);
223     std::string name = oss.str();
224     size_t pos = name.find('.');
225     if (std::string::npos != pos) {
226       name.replace(pos, 1, "p");
227     }
228     VarId res = addVarLocal(0.0, v, v, REAL, name);
229     if (fPhase1Over) {
230       addVar(res);
231     }
232     //       cerr << "  AddLitVar " << v << "   (PROBABLY WRONG)" << endl;
233     sLitValues[v] = res;
234     return res;
235   }
236   /// adding all local variables upfront. Makes sure it's called only once
addPhase1Vars()237   virtual void addPhase1Vars() {
238     assert(0 == getNColsModel());
239     assert(!fPhase1Over);
240     if (fVerbose) {
241       std::cerr << "  MIPWrapper: adding the " << colObj.size() << " Phase-1 variables..."
242                 << std::flush;
243     }
244     if (!colObj.empty()) {
245       doAddVars(colObj.size(), &colObj[0], &colLB[0], &colUB[0], &colTypes[0], &colNames[0]);
246     }
247     if (fVerbose) {
248       std::cerr << " done." << std::endl;
249     }
250     fPhase1Over = true;  // SCIP needs after adding
251   }
252 
253   /// var bounds
setVarBounds(int iVar,double lb,double ub)254   virtual void setVarBounds(int iVar, double lb, double ub) { throw 0; }
setVarLB(int iVar,double lb)255   virtual void setVarLB(int iVar, double lb) { throw 0; }
setVarUB(int iVar,double ub)256   virtual void setVarUB(int iVar, double ub) { throw 0; }
257   /// adding a linear constraint
258   virtual void addRow(int nnz, int* rmatind, double* rmatval, LinConType sense, double rhs,
259                       int mask = MaskConsType_Normal, const std::string& rowName = "") = 0;
260   /// Indicator constraint: x[iBVar]==bVal -> lin constr
addIndicatorConstraint(int iBVar,int bVal,int nnz,int * rmatind,double * rmatval,LinConType sense,double rhs,const std::string & rowName="")261   virtual void addIndicatorConstraint(int iBVar, int bVal, int nnz, int* rmatind, double* rmatval,
262                                       LinConType sense, double rhs,
263                                       const std::string& rowName = "") {
264     throw std::runtime_error("Indicator constraints not supported. ");
265   }
addMinimum(int iResultVar,int nnz,int * ind,const std::string & rowName="")266   virtual void addMinimum(int iResultVar, int nnz, int* ind, const std::string& rowName = "") {
267     throw std::runtime_error("This backend does not support the Minimum constraint");
268   }
269 
270   /// Bounds disj for SCIP
addBoundsDisj(int n,double * fUB,double * bnd,int * vars,int nF,double * fUBF,double * bndF,int * varsF,const std::string & rowName="")271   virtual void addBoundsDisj(int n, double* fUB, double* bnd, int* vars, int nF, double* fUBF,
272                              double* bndF, int* varsF, const std::string& rowName = "") {
273     throw std::runtime_error("Bounds disjunctions not supported. ");
274   }
275   /// Times constraint: var[x]*var[y] == var[z]
addTimes(int x,int y,int z,const std::string & rowName="")276   virtual void addTimes(int x, int y, int z, const std::string& rowName = "") {
277     throw std::runtime_error("Backend: [int/float]_times not supported. ");
278   }
279 
280   /// Cumulative, currently SCIP only
addCumulative(int nnz,int * rmatind,double * d,double * r,double b,const std::string & rowName="")281   virtual void addCumulative(int nnz, int* rmatind, double* d, double* r, double b,
282                              const std::string& rowName = "") {
283     throw std::runtime_error("Cumulative constraints not supported. ");
284   }
285 
286   /// Lex-lesseq binary, currently SCIP only
addLexLesseq(int nnz,int * rmatind1,int * rmatind2,bool isModelCons,const std::string & rowName="")287   virtual void addLexLesseq(int nnz, int* rmatind1, int* rmatind2, bool isModelCons,
288                             const std::string& rowName = "") {
289     throw std::runtime_error("MIP: lex_lesseq built-in not supported. ");
290   }
291 
292   /// Lex-chain-lesseq binary, currently SCIP only
addLexChainLesseq(int m,int n,int * rmatind,int nOrbitopeType,bool resolveprop,bool isModelCons,const std::string & rowName="")293   virtual void addLexChainLesseq(int m, int n, int* rmatind, int nOrbitopeType, bool resolveprop,
294                                  bool isModelCons, const std::string& rowName = "") {
295     throw std::runtime_error("MIP: lex_chain_lesseq built-in not supported. ");
296   }
297 
298   /// 0: model-defined level, 1: free, 2: uniform search
getFreeSearch()299   virtual int getFreeSearch() { return SearchType::FREE_SEARCH; }
300   /// Return 0 if ignoring searches
addSearch(const std::vector<VarId> & vars,const std::vector<int> & pri)301   virtual bool addSearch(const std::vector<VarId>& vars, const std::vector<int>& pri) {
302     return false;
303   }
304   /// Return 0 if ignoring warm starts
addWarmStart(const std::vector<VarId> & vars,const std::vector<double> & vals)305   virtual bool addWarmStart(const std::vector<VarId>& vars, const std::vector<double>& vals) {
306     return false;
307   }
308 
309   using MultipleObjectives = MiniZinc::MultipleObjectivesTemplate<VarId>;
defineMultipleObjectives(const MultipleObjectives & mo)310   virtual bool defineMultipleObjectives(const MultipleObjectives& mo) { return false; }
311 
312   int nAddedRows = 0;  // for name counting
313   int nIndicatorConstr = 0;
314   /// adding an implication
315   //     virtual void addImpl() = 0;
316   virtual void setObjSense(int s) = 0;  // +/-1 for max/min
317 
318   virtual double getInfBound() = 0;
319 
320   virtual int getNCols() = 0;
getNColsModel()321   virtual int getNColsModel() { return getNCols(); }  // from the solver
322   virtual int getNRows() = 0;
323 
324   //     void setObjUB(double ub) { objUB = ub; }
325   //     void addQPUniform(double c) { qpu = c; } // also sets problem type to MIQP unless c=0
326 
327   /// Set solution callback. Thread-safety??
328   /// solution callback handler, the wrapper might not have these callbacks implemented
provideSolutionCallback(SolCallbackFn cbfn,void * info)329   virtual void provideSolutionCallback(SolCallbackFn cbfn, void* info) {
330     assert(cbfn);
331     cbui.pOutput = &output;
332     cbui.psi = info;
333     cbui.solcbfn = cbfn;
334   }
335   /// solution callback handler, the wrapper might not have these callbacks implemented
provideCutCallback(CutCallbackFn cbfn,void * info)336   virtual void provideCutCallback(CutCallbackFn cbfn, void* info) {
337     assert(cbfn);
338     cbui.pCutOutput = nullptr;  // &outpCuts;   thread-safety: caller has to provide this
339     cbui.psi = info;
340     cbui.cutcbfn = cbfn;
341   }
342 
343   virtual void solve() = 0;
344 
345   /// OUTPUT, should also work in a callback
346   virtual const double* getValues() = 0;
347   virtual double getObjValue() = 0;
348   virtual double getBestBound() = 0;
getWallTimeElapsed()349   virtual double getWallTimeElapsed() { return output.dWallTime; }
350   virtual double getCPUTime() = 0;
351 
352   virtual Status getStatus() = 0;
353   virtual std::string getStatusName() = 0;
354 
355   virtual int getNNodes() = 0;
356   virtual int getNOpen() = 0;
357 
358   /// Default MZN library for MIP
359   static std::string getMznLib();
360 };
361