1 /*===========================================================================*
2  * This file is part of the Discrete Conic Optimization (DisCO) Solver.      *
3  *                                                                           *
4  * DisCO is distributed under the Eclipse Public License as part of the      *
5  * COIN-OR repository (http://www.coin-or.org).                              *
6  *                                                                           *
7  * Authors:                                                                  *
8  *          Aykut Bulut, Lehigh University                                   *
9  *          Yan Xu, Lehigh University                                        *
10  *          Ted Ralphs, Lehigh University                                    *
11  *                                                                           *
12  * Copyright (C) 2001-2018, Lehigh University, Aykut Bulut, Yan Xu, and      *
13  *                          Ted Ralphs.                                      *
14  * All Rights Reserved.                                                      *
15  *===========================================================================*/
16 
17 
18 // CoinUtils
19 #include <CoinMpsIO.hpp>
20 
21 // Disco headers
22 #include "DcoModel.hpp"
23 #include "DcoMessage.hpp"
24 #include "DcoTreeNode.hpp"
25 #include "DcoNodeDesc.hpp"
26 #include "DcoVariable.hpp"
27 #include "DcoLinearConstraint.hpp"
28 #include "DcoConicConstraint.hpp"
29 #include "DcoBranchStrategyMaxInf.hpp"
30 #include "DcoBranchStrategyPseudo.hpp"
31 #include "DcoBranchStrategyStrong.hpp"
32 #include "DcoConGenerator.hpp"
33 #include "DcoLinearConGenerator.hpp"
34 #include "DcoConicConGenerator.hpp"
35 #include "DcoSolution.hpp"
36 #include "DcoPresolve.hpp"
37 #include "DcoHeuristic.hpp"
38 #include "DcoHeurRounding.hpp"
39 #include "DcoCbfIO.hpp"
40 
41 // MILP cuts
42 #include <CglCutGenerator.hpp>
43 #include <CglProbing.hpp>
44 #include <CglClique.hpp>
45 #include <CglOddHole.hpp>
46 #include <CglFlowCover.hpp>
47 #include <CglKnapsackCover.hpp>
48 #include <CglMixedIntegerRounding2.hpp>
49 #include <CglGomory.hpp>
50 #include <CglGMI.hpp>
51 #include <CglTwomir.hpp>
52 
53 // Conic cuts
54 #include <CglConicCutGenerator.hpp>
55 #include <CglConicIPM.hpp>
56 #include <CglConicIPMint.hpp>
57 #include <CglConicOA.hpp>
58 
59 // STL headers
60 #include <string>
61 #include <iostream>
62 #include <fstream>
63 #include <sstream>
64 #include <numeric>
65 #include <cmath>
66 #include <iomanip>
67 
68 // ordering of conNames should match the ordering of DcoConstraintType enum
69 // type.
70 char const * conNames[] = {
71   "NotSet",
72   "Core",
73   "Clique",
74   "FCover",
75   "Gomory",
76   //"GMI",
77   "Knap",
78   // Linear MIR
79   "MIR",
80   "OddHole",
81   "Probe",
82   "TwoMIR",
83   // Conic cuts
84   "IPM",
85   "IPMint",
86   "OA",
87   // Conic MIR
88   "CMIR",
89   "GD1"
90 };
91 
92 std::vector<char const *> const
93   dcoConstraintTypeName (conNames, conNames + DcoConstraintTypeEnd);
94 
DcoModel()95 DcoModel::DcoModel() {
96   problemName_ = "";
97   solver_ = NULL;
98   colLB_ = NULL;
99   colUB_ = NULL;
100   rowLB_ = NULL;
101   rowUB_ = NULL;
102   numCols_ = 0;
103   numRows_ = 0;
104   numLinearRows_ = 0;
105   numConicRows_ = 0;
106   objSense_ = 0.0;
107   objCoef_ = NULL;
108   numIntegerCols_ = 0;
109   integerCols_ = NULL;
110   isInteger_ = NULL;
111   matrix_ = NULL;
112   coneStart_ = NULL;
113   coneMembers_ = NULL;
114   coneType_ = NULL;
115 
116   dcoPar_ = new DcoParams();
117   numRelaxedCols_ = 0;
118   relaxedCols_ = NULL;
119   numRelaxedRows_ = 0;
120   relaxedRows_ = NULL;
121   dcoMessageHandler_ = new CoinMessageHandler();
122   dcoMessages_ = new DcoMessage();
123   // set branch strategy
124   branchStrategy_ = NULL;
125   rampUpBranchStrategy_ = NULL;
126   // cut and heuristics objects will be set in setupSelf.
127 
128   initOAcuts_ = 0;
129 
130   dcoMessageHandler_->setPrefix(0);
131   dcoMessageHandler_->message(DISCO_WELCOME, *dcoMessages_)
132     << DISCO_VERSION
133     << __DATE__
134     << CoinMessageEol;
135   dcoMessageHandler_->setPrefix(1);
136 }
137 
~DcoModel()138 DcoModel::~DcoModel() {
139   // solver_ is freed in main function.
140   if (colLB_) {
141     delete[] colLB_;
142     colLB_=NULL;
143   }
144   if (colUB_) {
145     delete[] colUB_;
146     colUB_=NULL;
147   }
148   if (rowLB_) {
149     delete[] rowLB_;
150     rowLB_=NULL;
151   }
152   if (rowUB_) {
153     delete[] rowUB_;
154     rowUB_=NULL;
155   }
156   if (objCoef_) {
157     delete[] objCoef_;
158     objCoef_=NULL;
159   }
160   if (integerCols_) {
161     delete[] integerCols_;
162     integerCols_=NULL;
163   }
164   if (isInteger_) {
165     delete[] isInteger_;
166     isInteger_=NULL;
167   }
168   if (matrix_) {
169     delete matrix_;
170     matrix_=NULL;
171   }
172   if (coneStart_) {
173     delete[] coneStart_;
174     coneStart_=NULL;
175   }
176   if (coneMembers_) {
177     delete[] coneMembers_;
178     coneMembers_=NULL;
179   }
180   if (coneType_) {
181     delete[] coneType_;
182     coneType_=NULL;
183   }
184   if (branchStrategy_) {
185     delete branchStrategy_;
186     branchStrategy_=NULL;
187   }
188   if (rampUpBranchStrategy_) {
189     delete rampUpBranchStrategy_;
190     rampUpBranchStrategy_=NULL;
191   }
192   if (dcoPar_) {
193     delete dcoPar_;
194     dcoPar_=NULL;
195   }
196   if (dcoMessageHandler_) {
197     delete dcoMessageHandler_;
198     dcoMessageHandler_=NULL;
199   }
200   if (dcoMessages_) {
201     delete dcoMessages_;
202     dcoMessages_=NULL;
203   }
204   if (relaxedCols_) {
205     delete[] relaxedCols_;
206     relaxedCols_=NULL;
207   }
208   if (relaxedRows_) {
209     delete[] relaxedRows_;
210     relaxedRows_=NULL;
211   }
212   std::map<DcoConstraintType, DcoConGenerator*>::iterator it;
213   for (it=conGenerators_.begin();
214        it!=conGenerators_.end(); ++it) {
215     delete it->second;
216   }
217   conGenerators_.clear();
218   for (std::vector<DcoHeuristic*>::iterator it=heuristics_.begin();
219        it!=heuristics_.end(); ++it) {
220     delete *it;
221   }
222   heuristics_.clear();
223 }
224 
225 #if defined(__OA__)
setSolver(OsiSolverInterface * solver)226 void DcoModel::setSolver(OsiSolverInterface * solver) {
227   solver_ = solver;
228 }
229 #else
setSolver(OsiConicSolverInterface * solver)230 void DcoModel::setSolver(OsiConicSolverInterface * solver) {
231   solver_ = solver;
232 }
233 #endif
234 
235 // reads problem from the given file and sets the fields required by setupself
236 // only.
237 // setupSelf needs dcoPar, objSense_, variables_ and constraints_
238 // dcoPar_ is already set up by the AlpsKnowledgeBroker::initializeSearch().
readInstance(char const * dataFile)239 void DcoModel::readInstance(char const * dataFile) {
240   // get input file name
241   std::string input_file(dataFile);
242   std::string base_name = input_file.substr(0, input_file.rfind('.'));
243   std::string extension = input_file.substr(input_file.rfind('.')+1);
244   if (!extension.compare("mps")) {
245     readInstanceMps(dataFile);
246   }
247   else if (!extension.compare("cbf")) {
248     problemName_ = base_name;
249     readInstanceCbf(dataFile);
250   }
251   else {
252     dcoMessageHandler_->message(DISCO_READ_MPSCBFFILEONLY,
253                                 *dcoMessages_) << CoinMessageEol;
254   }
255 
256   // == log cone information messages
257   if (numConicRows_) {
258     dcoMessageHandler_->message(DISCO_READ_CONESTATS1,
259                                 *dcoMessages_) << numConicRows_
260                                                << CoinMessageEol;
261     for (int i=0; i<numConicRows_; ++i) {
262       dcoMessageHandler_->message(DISCO_READ_CONESTATS2,
263                                   *dcoMessages_)
264         << i
265         << coneStart_[i+1] - coneStart_[i]
266         << coneType_[i]
267         << CoinMessageEol;
268     }
269   }
270   else {
271     dcoMessageHandler_->message(DISCO_READ_NOCONES,
272                                 *dcoMessages_);
273   }
274 
275   // log problem information
276   std::string sense = (dcoPar_->entry(DcoParams::objSense)==1.0) ? std::string("min") : std::string("min");
277   dcoMessageHandler_->message(DISCO_PROBLEM_INFO,
278                               *dcoMessages_)
279     << problemName_
280     << sense.c_str()
281     << numCols_
282     << numLinearRows_
283     << matrix_->getNumElements()
284     << numConicRows_
285     << numIntegerCols_
286     << CoinMessageEol;
287 }
288 
289 // this should go into OsiConicSolverInterface or CoinUtils?
readInstanceCbf(char const * dataFile)290 void DcoModel::readInstanceCbf(char const * dataFile) {
291   // mps file reader
292   DcoCbfIO * reader = new DcoCbfIO();
293   reader->readCbf(dataFile);
294   // set objective sense
295   objSense_ = reader->objSense();
296   // set dcoPar_
297   dcoPar_->setEntry(DcoParams::objSense, objSense_);
298 
299   reader->getProblem(colLB_, colUB_, rowLB_, rowUB_, matrix_,
300                      numConicRows_, coneStart_, coneMembers_, coneType_);
301   numCols_ = matrix_->getNumCols();
302   numLinearRows_ = matrix_->getNumRows();
303   numRows_ = numLinearRows_ + numConicRows_;
304 
305   // add conic row bounds to rowLB and rowUB
306   double * tempRowLB = new double[numRows_];
307   std::copy(rowLB_, rowLB_+numLinearRows_, tempRowLB);
308   std::fill_n(tempRowLB+numLinearRows_, numConicRows_, 0.0);
309   delete[] rowLB_;
310   rowLB_ = tempRowLB;
311   tempRowLB = NULL;
312   double * tempRowUB = new double[numRows_];
313   std::copy(rowUB_, rowUB_+numLinearRows_, tempRowUB);
314   std::fill_n(tempRowUB+numLinearRows_, numConicRows_, reader->getInfinity());
315   delete[] rowUB_;
316   rowUB_ = tempRowUB;
317   tempRowUB = NULL;
318 
319   objCoef_ = new double [numCols_]();
320 
321   std::copy(reader->objCoef(), reader->objCoef()+reader->getNumCols(),
322             objCoef_);
323 
324   // get integrality info
325   numIntegerCols_ = reader->getNumInteger();
326   integerCols_ = new int[numIntegerCols_];
327   std::copy(reader->integerCols(), reader->integerCols()+numIntegerCols_,
328             integerCols_);
329   isInteger_ = new int[numCols_]();
330   for (int i=0; i<numIntegerCols_; ++i) {
331     isInteger_[integerCols_[i]] = 1;
332   }
333 
334   delete reader;
335 }
336 
readInstanceMps(char const * dataFile)337 void DcoModel::readInstanceMps(char const * dataFile) {
338   // mps file reader
339   CoinMpsIO * reader = new CoinMpsIO;
340   // set reader log level
341   //reader->messageHandler()->setLogLevel(dcoPar_->entry(DcoParams::logLevel));
342   reader->messageHandler()->setLogLevel(0);
343   reader->readMps(dataFile, "");
344   numCols_ = reader->getNumCols();
345 
346   // allocate variable bounds
347   colLB_ = new double [numCols_];
348   colUB_ = new double [numCols_];
349   // set col bounds
350   std::copy(reader->getColLower(), reader->getColLower()+numCols_, colLB_);
351   std::copy(reader->getColUpper(), reader->getColUpper()+numCols_, colUB_);
352 
353   // set objective sense
354   // todo(aykut) we should ask reader about the objective sense
355   objSense_ = dcoPar_->entry(DcoParams::objSense);
356 
357   // set objective coefficients
358   objCoef_ = new double [numCols_];
359   double const * reader_obj = reader->getObjCoefficients();
360   std::copy(reader_obj, reader_obj+numCols_, objCoef_);
361 
362   // set integer columns
363   // get variable integrality constraints
364   numIntegerCols_ = 0;
365   integerCols_ = new int[numCols_];
366   isInteger_ = new int[numCols_];
367   // may return NULL
368   char const * is_integer = reader->integerColumns();
369   if (is_integer!=NULL and strcmp(is_integer, "")) {
370     for (int i=0; i<numCols_; ++i) {
371       if (is_integer[i]) {
372         integerCols_[numIntegerCols_++] = i;
373         isInteger_[i] = 1;
374       }
375       else {
376         isInteger_[i] = 0;
377       }
378     }
379   }
380   else {
381     dcoMessageHandler_->message(3000, "Dco",
382                                 "CoinMpsIO::integerColumns() did return "
383                                 "NULL pointer for "
384                                 " column types! This looks like a CoinMpsIO "
385                                 "bug to me. Please report! I will use "
386                                 " CoinMpsIO::isContinuous(int index) function "
387                                 "instead.",
388                                 'W', 0)
389       << CoinMessageEol;
390     for (int i=0; i<numCols_; ++i) {
391       if (!reader->isContinuous(i)) {
392         integerCols_[numIntegerCols_++] = i;
393         isInteger_[i] = 1;
394       }
395       else {
396         isInteger_[i] = 0;
397       }
398     }
399   }
400   // resize integerCols_
401   int * temp = new int[numIntegerCols_];
402   std::copy(integerCols_, integerCols_+numIntegerCols_, temp);
403   delete[] integerCols_;
404   integerCols_ = temp;
405 
406   // read conic part
407   int reader_return = reader->readConicMps(NULL, coneStart_, coneMembers_,
408                                            coneType_, numConicRows_);
409   // when there is no conic section, status is -3.
410   if (reader_return==-3) {
411     // no cones in problem
412   }
413   else if (reader_return!=0) {
414     dcoMessageHandler_->message(DISCO_READ_MPSERROR,
415                                 *dcoMessages_) << reader_return
416                                                << CoinMessageEol;
417   }
418   // store number of constraints in the problem
419   numLinearRows_ = reader->getNumRows();
420   numRows_ = numLinearRows_ + numConicRows_;
421 
422   // allocate row bounds
423   rowLB_ = new double [numRows_];
424   rowUB_ = new double [numRows_];
425   // set row bounds for linear rows
426   std::copy(reader->getRowLower(), reader->getRowLower()+numLinearRows_, rowLB_);
427   std::copy(reader->getRowUpper(), reader->getRowUpper()+numLinearRows_, rowUB_);
428   // set conic row bounds
429   std::fill_n(rowLB_+numLinearRows_, numConicRows_, 0.0);
430   std::fill_n(rowUB_+numLinearRows_, numConicRows_, reader->getInfinity());
431 
432   matrix_ = new CoinPackedMatrix(*reader->getMatrixByRow());
433   problemName_ = reader->getProblemName();
434 
435   // free Coin MPS reader
436   delete reader;
437 }
438 
439 
readParameters(const int argnum,const char * const * arglist)440 void DcoModel::readParameters(const int argnum,
441                               const char * const * arglist) {
442   AlpsPar()->readFromArglist(argnum, arglist);
443   dcoPar_->readFromArglist(argnum, arglist);
444   setMessageLevel();
445 }
446 
447 // Add variables to *this.
setupAddVariables()448 void DcoModel::setupAddVariables() {
449   // add variables
450   BcpsVariable ** variables = new BcpsVariable*[numCols_];
451   for (int i=0; i<numCols_; ++i) {
452     variables[i] = new DcoVariable(i, colLB_[i], colUB_[i],
453                                    colLB_[i], colUB_[i]);
454     if (isInteger_[i]) {
455       variables[i]->setIntType('I');
456     }
457     else {
458       variables[i]->setIntType('C');
459     }
460     variables[i]->setBroker(broker_);
461   }
462   setVariables(variables, numCols_);
463   // variables[i] are now owned by BcpsModel, do not free them.
464   delete[] variables;
465 }
466 
setupAddLinearConstraints()467 void DcoModel::setupAddLinearConstraints() {
468   // == add constraints to *this
469   BcpsConstraint ** constraints = new BcpsConstraint*[numLinearRows_];
470   int const * indices = matrix_->getIndices();
471   double const * values = matrix_->getElements();
472   int const * lengths = matrix_->getVectorLengths();
473   int const * starts = matrix_->getVectorStarts();
474   for (int i=0; i<numLinearRows_; ++i) {
475     constraints[i] = new DcoLinearConstraint(lengths[i], indices+starts[i],
476                                              values+starts[i], rowLB_[i],
477                                              rowUB_[i]);
478     constraints[i]->setBroker(broker_);
479   }
480   setConstraints(constraints, numLinearRows_);
481   // constraints[i] are owned by BcpsModel. Do not free them here.
482   delete[] constraints;
483 }
484 
setupAddConicConstraints()485 void DcoModel::setupAddConicConstraints() {
486   // iterate over cones and add them to the model
487   for (int i=0; i<numConicRows_; ++i) {
488     if (coneType_[i]!=1 and coneType_[i]!=2) {
489       dcoMessageHandler_->message(DISCO_READ_CONEERROR,
490                                   *dcoMessages_) << CoinMessageEol;
491     }
492     int num_members = coneStart_[i+1]-coneStart_[i];
493     if (coneType_[i]==2 and num_members<3) {
494       dcoMessageHandler_->message(DISCO_READ_ROTATEDCONESIZE,
495                                   *dcoMessages_) << CoinMessageEol;
496     }
497     DcoLorentzConeType type = DcoLorentzCone;
498     if (coneType_[i]==2) {
499       type = DcoRotatedLorentzCone;
500     }
501     else if (coneType_[i]!=1) {
502       dcoMessageHandler_->message(DISCO_UNKNOWN_CONETYPE,
503                                   *dcoMessages_)
504         << __FILE__ << __LINE__ << CoinMessageEol;
505     }
506     DcoConicConstraint * cc =
507       new DcoConicConstraint(type, num_members,
508                              coneMembers_+coneStart_[i]);
509     cc->setBroker(broker_);
510     addConstraint(cc);
511   }
512 }
513 
514 /// Write out parameters.
writeParameters(std::ostream & outstream) const515 void DcoModel::writeParameters(std::ostream& outstream) const {
516   outstream << "\n================================================"
517             <<std::endl;
518   outstream << "ALPS Parameters: " << std::endl;
519   AlpsPar_->writeToStream(outstream);
520   outstream << "\n================================================"
521             <<std::endl;
522   outstream << "DISCO Parameters: " << std::endl;
523   dcoPar_->writeToStream(outstream);
524 }
525 
preprocess()526 void DcoModel::preprocess() {
527   // // set message levels
528   // setMessageLevel();
529 
530 
531   // // some bounds are improved if updated is true
532 
533   // // notes(aykut) this can only improve the bounds if the leading variable is
534   // // bounded.  most of the time it is not. Things can get better if we use some
535   // // other bound improvement first and it does improve the upper bound of
536   // // leading variables.
537   // //bool updated = DcoPresolve::improve_bounds(this);
538 
539   // write parameters used
540   //writeParameters(std::cout);
541 
542   // approximation of cones will update numLinearRows_, numRows_, rowLB_,
543   // rowUB_, matrix_.
544   if (numConicRows_ > 0) {
545     approximateCones();
546   }
547 }
548 
approximateCones()549 void DcoModel::approximateCones() {
550 #ifdef __OA__
551   // need to load problem to the solver.
552 
553   // load problem to the solver
554   solver_->loadProblem(*matrix_, colLB_, colUB_, objCoef_,
555                        rowLB_, rowUB_);
556   bool dual_infeasible = false;
557   int iter = 0;
558   int ipm_iter;
559   int oa_iter;
560   int num_ipm_cuts = 0;
561   int num_oa_cuts = 0;
562   // solve problem
563   solver_->resolve();
564   // get cone data in the required form
565   // todo(aykut) think about updating cut library for the input format
566   OsiLorentzConeType * coneTypes = new OsiLorentzConeType[numConicRows_];
567   int * coneSizes = new int[numConicRows_];
568   int const ** coneMembers = new int const *[numConicRows_];
569   for (int i=0; i<numConicRows_; ++i) {
570     if (coneType_[i]==1) {
571       coneTypes[i] = OSI_QUAD;
572     }
573     else if (coneType_[i]==2) {
574       coneTypes[i] = OSI_RQUAD;
575     }
576     else {
577       dcoMessageHandler_->message(DISCO_UNKNOWN_CONETYPE, *dcoMessages_)
578         << __FILE__ << __LINE__ << CoinMessageEol;
579     }
580     coneSizes[i] = coneStart_[i+1]-coneStart_[i];
581     coneMembers[i] = coneMembers_ + coneStart_[i];
582   }
583   // used to decide on number of iterations in outer approximation
584   int largest_cone_size = *std::max_element(coneSizes,
585                                             coneSizes+numConicRows_);
586   do {
587     // generate cuts
588     OsiCuts * ipm_cuts = new OsiCuts();
589     OsiCuts * oa_cuts = new OsiCuts();
590     CglConicCutGenerator * cg_ipm = new CglConicIPM();
591     CglConicCutGenerator * cg_oa =
592       new CglConicOA(dcoPar_->entry(DcoParams::coneTol));
593     // get cone info
594     cg_ipm->generateCuts(*solver_, *ipm_cuts, numConicRows_, coneTypes,
595                          coneSizes, coneMembers, largest_cone_size);
596     // cg_oa->generateCuts(*solver_, *oa_cuts, numCoreCones_, coneTypes_,
597     //                   coneSizes_, coneMembers_, largest_cone_size);
598     // if we do not get any cuts break the loop
599     if (ipm_cuts->sizeRowCuts()==0 && oa_cuts->sizeRowCuts()==0) {
600       break;
601     }
602     // if problem is unbounded do nothing, add cuts to the problem
603     // this will make lp relaxation infeasible
604     solver_->applyCuts(*ipm_cuts);
605     solver_->applyCuts(*oa_cuts);
606     solver_->resolve();
607     num_ipm_cuts += ipm_cuts->sizeRowCuts();
608     delete ipm_cuts;
609     delete oa_cuts;
610     delete cg_ipm;
611     delete cg_oa;
612     dual_infeasible = solver_->isProvenDualInfeasible();
613     iter++;
614   } while(dual_infeasible);
615   // add outer approximating cuts for DcoParams::approxNumPass (default is 50)
616   // many rounds
617   ipm_iter = iter;
618   iter = 0;
619   int oa_iter_limit = dcoPar_->entry(DcoParams::approxNumPass);
620   while(iter<oa_iter_limit) {
621     OsiCuts * oa_cuts = new OsiCuts();
622     CglConicCutGenerator * cg_oa =
623       new CglConicOA(dcoPar_->entry(DcoParams::coneTol));
624     cg_oa->generateCuts(*solver_, *oa_cuts, numConicRows_, coneTypes,
625                         coneSizes, coneMembers, 1);
626     int num_cuts = oa_cuts->sizeRowCuts();
627     num_oa_cuts += num_cuts;
628     if (num_cuts==0) {
629       // ifno cuts are produced break early
630       delete oa_cuts;
631       delete cg_oa;
632       break;
633     }
634     solver_->applyCuts(*oa_cuts);
635     solver_->resolve();
636     delete oa_cuts;
637     delete cg_oa;
638     iter++;
639   }
640   oa_iter = iter;
641   std::cout << "===== Preprocessing Summary =====" << std::endl;
642   std::cout << "IPM iterations " << ipm_iter << std::endl;
643   std::cout << "IPM cuts " << num_ipm_cuts << std::endl;
644   std::cout << "OA iterations " << oa_iter << std::endl;
645   std::cout << "OA cuts " << num_oa_cuts << std::endl;
646   std::cout << "Linear relaxation objective value "
647             << solver_->getObjValue() << std::endl;
648   std::cout << "=================================" << std::endl;
649   delete[] coneTypes;
650   delete[] coneSizes;
651   delete[] coneMembers;
652 
653   double cutOaSlack = dcoPar_->entry(DcoParams::cutOaSlack1);
654   // print cut activity
655   // {
656   //   int numCuts = solver_->getNumRows() - numLinearRows_;
657   //   double const * activity = solver_->getRowActivity();
658   //   double const * lb = solver_->getRowLower();
659   //   double const * ub = solver_->getRowUpper();
660   //   CoinWarmStartBasis const * ws =
661   //     dynamic_cast<CoinWarmStartBasis*> (solver_->getWarmStart());
662     //for (int i=numLinearRows_; i<solver_->getNumRows(); ++i) {
663       // double norm = 0.0;
664       // {
665       //   // compute norm of cut
666       //   CoinPackedMatrix const * mat = solver_->getMatrixByRow();
667       //   int start = mat->getVectorStarts()[i];
668       //   int size = mat->getVectorLengths()[i];
669       //   double const * value = mat->getElements() + start;
670       //   for (int j=0; j<size; ++j) {
671       //     norm += fabs(value[j]);
672       //   }
673       // }
674 
675       // std::cout << "row: " << std::setw(4) << i
676       //           << " lb: " << std::setw(12) << lb[i]
677       //           << " activity: " << std::setw(12) << activity[i]
678       //           << " ub: " << std::setw(12) << ub[i]
679       //           << " status: " << (ws->getArtifStatus(i)==CoinWarmStartBasis::basic)
680       //           << " remove: " << ((ub[i]-activity[i])>cutOaSlack)
681       //           << " norm: " << norm
682       //           << std::endl;
683     //}
684   //}
685 
686 
687   // remove inactive cuts
688   {
689     int numCuts = solver_->getNumRows() - numLinearRows_;
690     // get warm start basis
691     CoinWarmStartBasis * ws =
692       dynamic_cast<CoinWarmStartBasis*> (solver_->getWarmStart());
693     if (ws==NULL) {
694       // nothing to do if there is no warm start information
695       std::cerr << "Disco warning: No warm start object exists in solver. "
696                 << "Unable to clean cuts." << std::endl;
697     }
698     else if (numCuts==0) {
699       // nothing to do if there no any cuts
700     }
701     else {
702       int numDel = 0;
703       int * delInd = new int[numCuts];
704       double const * activity = solver_->getRowActivity();
705       //double const * lb = solver_->getRowLower();
706       double const * ub = solver_->getRowUpper();
707 
708       // iterate over cuts and colect inactive ones
709       for (int i=0; i<numCuts; ++i) {
710         if (ws->getArtifStatus(numLinearRows_+i) == CoinWarmStartBasis::basic) {
711           // cut is inactive
712           // do we really want to remove? check how much it is inactive
713           if (ub[numLinearRows_+i] - activity[numLinearRows_+i] > cutOaSlack) {
714             // remove since it is inactive too much
715             delInd[numDel++] = i+numLinearRows_;
716           }
717         }
718       }
719       // delete cuts from solver
720       if (numDel) {
721         std::cout << "Approx cones: "
722                   << " removed: " << numDel
723                   << " remain: " << numCuts-numDel
724                   << std::endl;
725         solver_->deleteRows(numDel, delInd);
726         // resolve to correct status
727         solver_->resolve();
728       }
729       delete[] delInd;
730     }
731     delete ws;
732   }
733   initOAcuts_ = solver_->getNumRows() - numLinearRows_;
734 
735   // get updated data from solver
736   // delete matrix_;
737   // matrix_ = new CoinPackedMatrix(*solver_->getMatrixByRow());
738   // // update number of rows
739   // numLinearRows_ = solver_->getNumRows();
740   // numRows_ = numLinearRows_+numConicRows_;
741   // // update row bounds
742   // delete[] rowLB_;
743   // rowLB_ = new double[numRows_];
744   // std::copy(solver_->getRowLower(),
745   //           solver_->getRowLower()+numLinearRows_, rowLB_);
746   // std::fill_n(rowLB_+numLinearRows_, numConicRows_, 0.0);
747   // delete[] rowUB_;
748   // rowUB_ = new double[numRows_];
749   // std::copy(solver_->getRowUpper(),
750   //           solver_->getRowUpper()+numLinearRows_, rowUB_);
751   // std::fill_n(rowUB_+numLinearRows_, numConicRows_, solver_->getInfinity());
752 #endif
753 }
754 
755 //todo(aykut) why does this return to bool?
756 // should be fixed in Alps level.
757 
758 // relase redundant memory once the fields are set.
setupSelf()759 bool DcoModel::setupSelf() {
760   // set message level again. We call this in readInstance() too.
761   // this is so due to parallelization.
762   setMessageLevel();
763   // set relaxed array for integer columns
764   numRelaxedCols_ = numIntegerCols_;
765   relaxedCols_ = new int[numRelaxedCols_];
766   std::copy(integerCols_, integerCols_+numIntegerCols_,
767             relaxedCols_);
768   // set iteration count to 0
769   numRelaxIterations_ = 0;
770 #ifdef __OA__
771   solver_->reset();
772   solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL);
773   // clp specific options for getting unboundedness directions
774   dynamic_cast<OsiClpSolverInterface*>(solver_)->getModelPtr()->setMoreSpecialOptions(0);
775   dynamic_cast<OsiClpSolverInterface*>(solver_)->getModelPtr()->setLogLevel(0);
776 #else
777   solver_->setHintParam(OsiDoReducePrint, true, OsiHintTry);
778 #endif
779   // load problem to the solver
780   solver_->loadProblem(*matrix_, colLB_, colUB_, objCoef_,
781                        rowLB_, rowUB_);
782   // set integer variables
783   solver_->setInteger(integerCols_, numIntegerCols_);
784 
785 #if defined(__OA__)
786   // we relax conic constraints when OA is used.
787   // set relaxed array for conic constraints
788   numRelaxedRows_ = numConicRows_;
789   relaxedRows_ = new int[numRelaxedRows_];
790   // notes(aykut) we assume conic rows start after linear rows.
791   for (int i=0; i<numRelaxedRows_; ++i) {
792     relaxedRows_[i] = numLinearRows_+i;
793   }
794 
795   // set leading variable lower bounds to 0
796   for (int i=0; i<numConicRows_; ++i) {
797     if (coneType_[i]==1) {
798       colLB_[coneMembers_[coneStart_[i]]] = 0.0;
799       solver_->setColLower(coneMembers_[coneStart_[i]], 0.0);
800     }
801     else if (coneType_[i]==2) {
802       colLB_[coneMembers_[coneStart_[i]]] = 0.0;
803       colLB_[coneMembers_[coneStart_[i]+1]] = 0.0;
804       solver_->setColLower(coneMembers_[coneStart_[i]], 0.0);
805       solver_->setColLower(coneMembers_[coneStart_[i]+1], 0.0);
806     }
807     else {
808       dcoMessageHandler_->message(DISCO_UNKNOWN_CONETYPE,
809                                   *dcoMessages_)
810         << coneType_[i] << CoinMessageEol;
811     }
812   }
813 #else
814   // add conic constraints to the solver
815   for (int i=0; i<numConicRows_; ++i) {
816     // do not relax conic constraints, add them to the conic solver
817     OsiLorentzConeType osi_type;
818     if (coneType_[i]==1) {
819       osi_type = OSI_QUAD;
820     }
821     else if (coneType_[i]==2) {
822       osi_type = OSI_RQUAD;
823     }
824     solver_->addConicConstraint(osi_type, coneStart_[i+1]-coneStart_[i],
825                               coneMembers_+coneStart_[i]);
826   }
827 #endif
828 
829   // create disco variables
830   setupAddVariables();
831   // create disco constraints, linear
832   setupAddLinearConstraints();
833   // create disco constraints, conic
834   setupAddConicConstraints();
835 
836   // set branch strategy
837   setBranchingStrategy();
838 
839   // add constraint generators
840 #ifdef __OA__
841   addConstraintGenerators();
842 #endif
843 
844   // add heuristics
845   addHeuristics();
846 
847   // set cutoff for solvers if user provides one.
848   double cutoff = dcoPar_->entry(DcoParams::cutoff);
849   if (cutoff!=ALPS_INC_MAX) {
850     solver_->setDblParam(OsiDualObjectiveLimit, objSense_*cutoff);
851   }
852   return true;
853 }
854 
855 // set message level
setMessageLevel()856 void DcoModel::setMessageLevel() {
857   // get Alps log level
858   int alps_log_level = AlpsPar()->entry(AlpsParams::msgLevel);
859   int dco_log_level = dcoPar_->entry(DcoParams::logLevel);
860 
861 #if defined(DISCO_DEBUG) || defined(DISCO_DEBUG_BRANCH) ||defined(DISCO_DEBUG_CUT) || defined(DISCO_DEBUG_PROCESS)
862   // reset Disco log level, we will set it depending on the debug macros are
863   // defined.
864   dco_log_level = 0;
865 #endif
866 
867 #ifdef DISCO_DEBUG_BRANCH
868   dco_log_level += DISCO_DLOG_BRANCH;
869 #endif
870 
871 #ifdef DISCO_DEBUG_CUT
872   dco_log_level += DISCO_DLOG_CUT;
873 #endif
874 
875 #ifdef DISCO_DEBUG_PROCESS
876   dco_log_level += DISCO_DLOG_PROCESS;
877 #endif
878 
879 #ifdef DISCO_DEBUG
880   // debug branch, cut, process
881   dco_log_level = DISCO_DLOG_BRANCH;
882   dco_log_level += DISCO_DLOG_CUT;
883   dco_log_level += DISCO_DLOG_PROCESS;
884 #endif
885 
886   // todo(aykut) create different parameters for Alps and Bcps.
887   broker()->messageHandler()->setLogLevel(alps_log_level);
888   bcpsMessageHandler()->setLogLevel(alps_log_level);
889   dcoMessageHandler_->setLogLevel(dco_log_level);
890 }
891 
addConstraintGenerators()892 void DcoModel::addConstraintGenerators() {
893   // get cut strategy
894   cutStrategy_ = static_cast<DcoCutStrategy>
895     (dcoPar_->entry(DcoParams::cutStrategy));
896   // get cut generation frequency
897   cutGenerationFrequency_ =
898     dcoPar_->entry(DcoParams::cutGenerationFrequency);
899   if (cutGenerationFrequency_ < 1) {
900     // invalid cut fraquency given, change it to 1.
901     dcoMessageHandler_->message(DISCO_INVALID_CUT_FREQUENCY,
902                                 *dcoMessages_)
903       << cutGenerationFrequency_
904       << 1
905       << CoinMessageEol;
906     cutGenerationFrequency_ = 1;
907   }
908 
909   // get cut strategies from parameters
910   DcoCutStrategy cliqueStrategy = static_cast<DcoCutStrategy>
911     (dcoPar_->entry(DcoParams::cutCliqueStrategy));
912   DcoCutStrategy fCoverStrategy = static_cast<DcoCutStrategy>
913     (dcoPar_->entry(DcoParams::cutFlowCoverStrategy));
914   DcoCutStrategy gomoryStrategy = static_cast<DcoCutStrategy>
915     (dcoPar_->entry(DcoParams::cutGomoryStrategy));
916   DcoCutStrategy knapStrategy = static_cast<DcoCutStrategy>
917     (dcoPar_->entry(DcoParams::cutKnapsackStrategy));
918   DcoCutStrategy mirStrategy = static_cast<DcoCutStrategy>
919     (dcoPar_->entry(DcoParams::cutMirStrategy));
920   DcoCutStrategy oddHoleStrategy = static_cast<DcoCutStrategy>
921     (dcoPar_->entry(DcoParams::cutOddHoleStrategy));
922   DcoCutStrategy probeStrategy = static_cast<DcoCutStrategy>
923     (dcoPar_->entry(DcoParams::cutProbingStrategy));
924   DcoCutStrategy twoMirStrategy = static_cast<DcoCutStrategy>
925     (dcoPar_->entry(DcoParams::cutTwoMirStrategy));
926   DcoCutStrategy ipmStrategy = static_cast<DcoCutStrategy>
927     (dcoPar_->entry(DcoParams::cutIpmStrategy));
928   DcoCutStrategy ipmintStrategy = static_cast<DcoCutStrategy>
929     (dcoPar_->entry(DcoParams::cutIpmIntStrategy));
930   DcoCutStrategy oaStrategy = static_cast<DcoCutStrategy>
931     (dcoPar_->entry(DcoParams::cutOaStrategy));
932 
933   // get cut frequencies from parameters
934   int cliqueFreq = dcoPar_->entry(DcoParams::cutCliqueFreq);
935   int fCoverFreq = dcoPar_->entry(DcoParams::cutFlowCoverFreq);
936   int gomoryFreq = dcoPar_->entry(DcoParams::cutGomoryFreq);
937   int knapFreq = dcoPar_->entry(DcoParams::cutKnapsackFreq);
938   int mirFreq = dcoPar_->entry(DcoParams::cutMirFreq);
939   int oddHoleFreq = dcoPar_->entry(DcoParams::cutOddHoleFreq);
940   int probeFreq = dcoPar_->entry(DcoParams::cutProbingFreq);
941   int twoMirFreq = dcoPar_->entry(DcoParams::cutTwoMirFreq);
942   int ipmFreq = dcoPar_->entry(DcoParams::cutIpmFreq);
943   int ipmintFreq = dcoPar_->entry(DcoParams::cutIpmIntFreq);
944   int oaFreq = dcoPar_->entry(DcoParams::cutOaFreq);
945 
946   //----------------------------------
947   // Add cut generators.
948   //----------------------------------
949   // Add probe cut generator
950   if (probeStrategy == DcoCutStrategyNotSet) {
951     // Disable by default
952     if (cutStrategy_ == DcoCutStrategyNotSet) {
953       probeStrategy = DcoCutStrategyNone;
954     }
955     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
956       probeStrategy = cutStrategy_;
957       probeFreq = cutGenerationFrequency_;
958     }
959     else {
960       probeStrategy = cutStrategy_;
961     }
962   }
963   if (probeStrategy != DcoCutStrategyNone) {
964     CglProbing *probing = new CglProbing;
965     probing->setUsingObjective(true);
966     probing->setMaxPass(1);
967     probing->setMaxPassRoot(5);
968     // Number of unsatisfied variables to look at
969     probing->setMaxProbe(10);
970     probing->setMaxProbeRoot(1000);
971     // How far to follow the consequences
972     probing->setMaxLook(50);
973     probing->setMaxLookRoot(500);
974     // Only look at rows with fewer than this number of elements
975     probing->setMaxElements(200);
976     probing->setRowCuts(3);
977     addConGenerator(probing, DcoConstraintTypeProbe, probeStrategy, probeFreq);
978   }
979 
980   // Add clique cut generator.
981   if (cliqueStrategy == DcoCutStrategyNotSet) {
982     // Disable by default
983     if (cutStrategy_ == DcoCutStrategyNotSet) {
984       cliqueStrategy = DcoCutStrategyNone;
985     }
986     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
987       cliqueFreq = cutGenerationFrequency_;
988       cliqueStrategy = DcoCutStrategyPeriodic;
989     }
990     else { // Root or Auto
991       cliqueStrategy = cutStrategy_;
992     }
993   }
994   if (cliqueStrategy != DcoCutStrategyNone) {
995     CglClique *cliqueCut = new CglClique ;
996     cliqueCut->setStarCliqueReport(false);
997     cliqueCut->setRowCliqueReport(false);
998     addConGenerator(cliqueCut, DcoConstraintTypeClique, cliqueStrategy,
999                     cliqueFreq);
1000   }
1001 
1002   // Add odd hole cut generator.
1003   if (oddHoleStrategy == DcoCutStrategyNotSet) {
1004     if (cutStrategy_ == DcoCutStrategyNotSet) {
1005       // Disable by default
1006       oddHoleStrategy = DcoCutStrategyNone;
1007     }
1008     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1009       oddHoleStrategy = DcoCutStrategyPeriodic;
1010       oddHoleFreq = cutGenerationFrequency_;
1011     }
1012     else {
1013       oddHoleStrategy = cutStrategy_;
1014     }
1015   }
1016   if (oddHoleStrategy != DcoCutStrategyNone) {
1017     CglOddHole *oldHoleCut = new CglOddHole;
1018     oldHoleCut->setMinimumViolation(0.005);
1019     oldHoleCut->setMinimumViolationPer(0.00002);
1020     // try larger limit
1021     oldHoleCut->setMaximumEntries(200);
1022     addConGenerator(oldHoleCut, DcoConstraintTypeOddHole, oddHoleStrategy,
1023                     oddHoleFreq);
1024   }
1025 
1026   // Add flow cover cut generator.
1027   if (fCoverStrategy == DcoCutStrategyNotSet) {
1028     if (cutStrategy_ == DcoCutStrategyNotSet) {
1029       // Disable by default
1030       fCoverStrategy = DcoCutStrategyNone;
1031       fCoverFreq = cutGenerationFrequency_;
1032     }
1033     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1034       fCoverStrategy = cutStrategy_;
1035       fCoverFreq = cutGenerationFrequency_;
1036     }
1037     else {
1038       fCoverStrategy = cutStrategy_;
1039     }
1040   }
1041   if (fCoverStrategy != DcoCutStrategyNone) {
1042     CglFlowCover *flowGen = new CglFlowCover;
1043     addConGenerator(flowGen, DcoConstraintTypeFCover, fCoverStrategy,
1044                     fCoverFreq);
1045   }
1046 
1047   // Add knapsack cut generator.
1048   if (knapStrategy == DcoCutStrategyNotSet) {
1049     if (cutStrategy_ == DcoCutStrategyNotSet) {
1050       // Disable by default
1051       knapStrategy = DcoCutStrategyNone;
1052     }
1053     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1054       knapStrategy = cutStrategy_;
1055       knapFreq = cutGenerationFrequency_;
1056     }
1057     else {
1058       knapStrategy = cutStrategy_;
1059     }
1060   }
1061   if (knapStrategy != DcoCutStrategyNone) {
1062     CglKnapsackCover *knapCut = new CglKnapsackCover;
1063     addConGenerator(knapCut, DcoConstraintTypeKnap, knapStrategy, knapFreq);
1064   }
1065 
1066   // Add MIR cut generator.
1067   if (mirStrategy == DcoCutStrategyNotSet) {
1068     if (cutStrategy_ == DcoCutStrategyNotSet) {
1069       // Disable by default
1070       mirStrategy = DcoCutStrategyNone;
1071     }
1072     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1073       mirStrategy = cutStrategy_;
1074       mirFreq = cutGenerationFrequency_;
1075     }
1076     else {
1077       mirStrategy = cutStrategy_;
1078     }
1079   }
1080   if (mirStrategy != DcoCutStrategyNone) {
1081     CglMixedIntegerRounding2 *mixedGen = new CglMixedIntegerRounding2;
1082     addConGenerator(mixedGen, DcoConstraintTypeMIR, mirStrategy, mirFreq);
1083   }
1084 
1085   // Add Gomory cut generator.
1086   if (gomoryStrategy == DcoCutStrategyNotSet) {
1087     if (cutStrategy_ == DcoCutStrategyNotSet) {
1088       // Disable by default
1089       gomoryStrategy = DcoCutStrategyNone;
1090     }
1091     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1092       gomoryStrategy = cutStrategy_;
1093       gomoryFreq = cutGenerationFrequency_;
1094     }
1095     else {
1096       gomoryStrategy = cutStrategy_;
1097     }
1098   }
1099   if (gomoryStrategy != DcoCutStrategyNone) {
1100     CglGomory *gomoryCut = new CglGomory;
1101     //CglGMI *gomoryCut = new CglGMI;
1102     // try larger limit
1103     gomoryCut->setLimit(40);
1104     addConGenerator(gomoryCut, DcoConstraintTypeGomory, gomoryStrategy,
1105                     gomoryFreq);
1106   }
1107 
1108   // Add Tow MIR cut generator.
1109   // Disable forever, not useful.
1110   twoMirStrategy = DcoCutStrategyNone;
1111   if (twoMirStrategy != DcoCutStrategyNone) {
1112     CglTwomir *twoMirCut =  new CglTwomir;
1113     addConGenerator(twoMirCut, DcoConstraintTypeTwoMIR, twoMirStrategy,
1114                     twoMirFreq);
1115   }
1116 
1117   // Add IPM cut generator
1118   if (ipmStrategy == DcoCutStrategyNotSet) {
1119     if (cutStrategy_ == DcoCutStrategyNotSet) {
1120       // Disable by default
1121       ipmStrategy = DcoCutStrategyNone;
1122     }
1123     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1124       ipmStrategy = cutStrategy_;
1125       ipmFreq = cutGenerationFrequency_;
1126     }
1127     else {
1128       ipmStrategy = cutStrategy_;
1129     }
1130   }
1131   if (ipmStrategy != DcoCutStrategyNone) {
1132     CglConicCutGenerator * ipm_gen = new CglConicIPM();
1133     addConGenerator(ipm_gen, DcoConstraintTypeIPM, ipmStrategy, ipmFreq);
1134   }
1135 
1136   // Add IPM integer cut generator
1137   if (ipmintStrategy == DcoCutStrategyNotSet) {
1138     if (cutStrategy_ == DcoCutStrategyNotSet) {
1139       // Disable by default
1140       ipmintStrategy = DcoCutStrategyNone;
1141     }
1142     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1143       ipmintStrategy = cutStrategy_;
1144       ipmintFreq = cutGenerationFrequency_;
1145     }
1146     else {
1147       ipmintStrategy = cutStrategy_;
1148     }
1149   }
1150   if (ipmintStrategy != DcoCutStrategyNone) {
1151     CglConicCutGenerator * ipm_int_gen = new CglConicIPMint();
1152     addConGenerator(ipm_int_gen, DcoConstraintTypeIPMint, ipmintStrategy,
1153                     ipmintFreq);
1154   }
1155 
1156   // Add Outer approximation cut generator
1157   if (oaStrategy == DcoCutStrategyNotSet) {
1158     if (cutStrategy_ == DcoCutStrategyNotSet) {
1159       // periodic by default
1160       oaStrategy = DcoCutStrategyPeriodic;
1161     }
1162     else if (cutStrategy_ == DcoCutStrategyPeriodic) {
1163       oaStrategy = cutStrategy_;
1164       oaFreq = cutGenerationFrequency_;
1165     }
1166     else {
1167       oaStrategy = cutStrategy_;
1168     }
1169   }
1170   if (oaStrategy != DcoCutStrategyNone && numConicRows_) {
1171     CglConicCutGenerator * oa_gen =
1172       new CglConicOA(dcoPar_->entry(DcoParams::coneTol));
1173     addConGenerator(oa_gen, DcoConstraintTypeOA, oaStrategy, oaFreq);
1174   }
1175 
1176   // Adjust cutStrategy_ according to the strategies of each cut generators.
1177   // set it to the most allowing one.
1178   // if there is at least one periodic strategy, set it to periodic.
1179   // if no periodic and there is at least one root, set it to root
1180   // set it to None otherwise.
1181   cutStrategy_ = DcoCutStrategyNone;
1182   cutGenerationFrequency_ = 100;
1183   bool periodic_exists = false;
1184   bool root_exists = false;
1185   std::map<DcoConstraintType, DcoConGenerator*>::iterator it;
1186   for (it=conGenerators_.begin(); it!=conGenerators_.end(); ++it) {
1187     DcoCutStrategy curr = it->second->strategy();
1188     if (curr==DcoCutStrategyPeriodic) {
1189       periodic_exists = true;
1190       break;
1191     }
1192     else if (curr==DcoCutStrategyRoot) {
1193       root_exists = true;
1194     }
1195   }
1196   if (periodic_exists) {
1197     // greatest divisor of all frequencies
1198     cutStrategy_ = DcoCutStrategyPeriodic;
1199     cutGenerationFrequency_ = 1;
1200   }
1201   else if (root_exists) {
1202     cutStrategy_ = DcoCutStrategyRoot;
1203     // this is not relevant, since we will generate only in root.
1204     cutGenerationFrequency_ = 100;
1205   }
1206 }
1207 
1208 
addConGenerator(CglCutGenerator * cgl_gen,DcoConstraintType type,DcoCutStrategy dco_strategy,int frequency)1209 void DcoModel::addConGenerator(CglCutGenerator * cgl_gen,
1210                                DcoConstraintType type,
1211                                DcoCutStrategy dco_strategy,
1212                                int frequency) {
1213   char const * name = dcoConstraintTypeName[type];
1214   DcoConGenerator * con_gen = new DcoLinearConGenerator(this, cgl_gen, type,
1215                                                         name,
1216                                                         dco_strategy,
1217                                                         frequency);
1218   conGenerators_[type] = con_gen;
1219 }
1220 
1221 /// Add constraint generator.
addConGenerator(CglConicCutGenerator * cgl_gen,DcoConstraintType type,DcoCutStrategy dco_strategy,int frequency)1222 void DcoModel::addConGenerator(CglConicCutGenerator * cgl_gen,
1223                                DcoConstraintType type,
1224                                DcoCutStrategy dco_strategy,
1225                                int frequency) {
1226   char const * name = dcoConstraintTypeName[type];
1227   DcoConGenerator * con_gen = new DcoConicConGenerator(this, cgl_gen, type,
1228                                                        name,
1229                                                        dco_strategy,
1230                                                        frequency);
1231   conGenerators_[type] = con_gen;
1232 }
1233 
addHeuristics()1234 void DcoModel::addHeuristics() {
1235   // todo(aykut) this function (heuristic adding process) can be improved
1236   // since global parameters are ignored with this design.
1237 
1238   heuristics_.clear();
1239   // get global heuristic strategy
1240   heurStrategy_ = static_cast<DcoHeurStrategy>
1241     (dcoPar_->entry(DcoParams::heurStrategy));
1242   // get global heuristic call frequency
1243   heurFrequency_ =
1244     dcoPar_->entry(DcoParams::heurCallFrequency);
1245   if (heurFrequency_ < 1) {
1246     // invalid heur fraquency given, change it to 1.
1247     dcoMessageHandler_->message(DISCO_INVALID_HEUR_FREQUENCY,
1248                                 *dcoMessages_)
1249       << heurFrequency_
1250       << 1
1251       << CoinMessageEol;
1252     heurFrequency_ = 1;
1253   }
1254 
1255   // get rounding heuristics strategy from parameters
1256   DcoHeurStrategy roundingStrategy = static_cast<DcoHeurStrategy>
1257     (dcoPar_->entry(DcoParams::heurRoundStrategy));
1258   // get rounding heuristic frequency from parameters
1259   int roundingFreq = dcoPar_->entry(DcoParams::heurRoundFreq);
1260 
1261   // add heuristics
1262   // == add rounding heuristics
1263   if (roundingStrategy != DcoHeurStrategyNone) {
1264     DcoHeuristic * round = new DcoHeurRounding(this, "rounding",
1265                                                roundingStrategy, roundingFreq);
1266     heuristics_.push_back(round);
1267   }
1268 
1269 
1270   // Adjust heurStrategy_ according to the strategies/frequencies of each
1271   // heuristic. Set it to the most allowing one.
1272   // if there is at least one periodic strategy, set it to periodic.
1273   // if no periodic and there is at least one root, set it to root
1274   // set it to None otherwise.
1275   heurStrategy_ = DcoHeurStrategyNone;
1276   heurFrequency_ = -1;
1277   bool periodic_exists = false;
1278   bool root_exists = false;
1279   std::vector<DcoHeuristic*>::iterator it;
1280   for (it=heuristics_.begin(); it!=heuristics_.end(); ++it) {
1281     DcoHeurStrategy curr = (*it)->strategy();
1282     if (curr==DcoHeurStrategyPeriodic) {
1283       periodic_exists = true;
1284       break;
1285     }
1286     else if (curr==DcoHeurStrategyRoot) {
1287       root_exists = true;
1288     }
1289   }
1290   if (periodic_exists) {
1291     heurStrategy_ = DcoHeurStrategyPeriodic;
1292     heurFrequency_ = 1;
1293   }
1294   else if (root_exists) {
1295     heurStrategy_ = DcoHeurStrategyRoot;
1296     // this is not relevant, since we will generate only in root.
1297     heurFrequency_ = -1;
1298   }
1299 }
1300 
setBranchingStrategy()1301 void DcoModel::setBranchingStrategy() {
1302     // set branching startegy
1303   int brStrategy = static_cast<DcoBranchingStrategy>
1304     (dcoPar_->entry(DcoParams::branchStrategy));
1305   switch(brStrategy) {
1306   case DcoBranchingStrategyMaxInfeasibility:
1307     branchStrategy_ = new DcoBranchStrategyMaxInf(this);
1308     break;
1309   case DcoBranchingStrategyPseudoCost:
1310     branchStrategy_ = new DcoBranchStrategyPseudo(this);
1311     break;
1312   // case DcoBranchingStrategyReliability:
1313   //   branchStrategy_ = new DcoBranchStrategyRel(this, reliability);
1314   //   break;
1315   case DcoBranchingStrategyStrong:
1316      branchStrategy_ = new DcoBranchStrategyStrong(this);
1317      break;
1318   // case DcoBranchingStrategyBilevel:
1319   //   branchStrategy_ = new DcoBranchStrategyBilevel(this);
1320   //   break;
1321   default:
1322     dcoMessageHandler_->message(DISCO_UNKNOWN_BRANCHSTRATEGY,
1323                                 *dcoMessages_) << brStrategy << CoinMessageEol;
1324     throw CoinError("Unknown branch strategy.", "setupSelf","DcoModel");
1325   }
1326 
1327   // set ramp up branch strategy
1328   brStrategy = static_cast<DcoBranchingStrategy>
1329     (dcoPar_->entry(DcoParams::branchStrategyRampUp));
1330   switch(brStrategy) {
1331   case DcoBranchingStrategyMaxInfeasibility:
1332     rampUpBranchStrategy_ = new DcoBranchStrategyMaxInf(this);
1333     break;
1334   case DcoBranchingStrategyPseudoCost:
1335     rampUpBranchStrategy_ = new DcoBranchStrategyPseudo(this);
1336     break;
1337   // case DcoBranchingStrategyReliability:
1338   //   rampUpBranchStrategy_ = new DcoBranchStrategyRel(this, reliability);
1339   //   break;
1340   case DcoBranchingStrategyStrong:
1341      rampUpBranchStrategy_ = new DcoBranchStrategyStrong(this);
1342      break;
1343   // case DcoBranchingStrategyBilevel:
1344   //   rampUpBranchStrategy_ = new DcoBranchStrategyBilevel(this);
1345   //   break;
1346   default:
1347     dcoMessageHandler_->message(DISCO_UNKNOWN_BRANCHSTRATEGY,
1348                                 *dcoMessages_) << brStrategy << CoinMessageEol;
1349     throw std::exception();
1350     throw CoinError("Unknown branch strategy.", "setupSelf","DcoModel");
1351   }
1352 }
1353 
postprocess()1354 void DcoModel::postprocess() {
1355 }
1356 
1357 
createRoot()1358 AlpsTreeNode * DcoModel::createRoot() {
1359   DcoTreeNode * root = new DcoTreeNode();
1360   DcoNodeDesc * desc = new DcoNodeDesc(this);
1361   root->setDesc(desc);
1362   std::vector<BcpsVariable *> & cols = getVariables();
1363   std::vector<BcpsConstraint *> & rows = getConstraints();
1364   int numCols = getNumCoreVariables();
1365   int numRows = getNumCoreConstraints();
1366   int * varIndices1 = new int [numCols];
1367   int * varIndices2 = new int [numCols];
1368   int * varIndices3 = NULL;
1369   int * varIndices4 = NULL;
1370   double * vlhe = new double [numCols];
1371   double * vuhe = new double [numCols];
1372   double * vlse = NULL;
1373   double * vuse = NULL;
1374   int * conIndices1 = new int [numRows];
1375   int * conIndices2 = new int [numRows];
1376   int * conIndices3 = NULL;
1377   int * conIndices4 = NULL;
1378   double * clhe = new double [numRows];
1379   double * cuhe = new double [numRows];
1380   double * clse = NULL;
1381   double * cuse = NULL;
1382   //-------------------------------------------------------------
1383   // Get var bounds and indices.
1384   //-------------------------------------------------------------
1385   for (int i=0; i<numCols; ++i) {
1386     vlhe[i] = cols[i]->getLbHard();
1387     vuhe[i] = cols[i]->getUbHard();
1388     varIndices1[i] = i;
1389     varIndices2[i] = i;
1390   }
1391   //-------------------------------------------------------------
1392   // Get con bounds and indices.
1393   //-------------------------------------------------------------
1394   for (int i=0; i<numRows; ++i) {
1395     clhe[i] = rows[i]->getLbHard();
1396     cuhe[i] = rows[i]->getUbHard();
1397     conIndices1[i] = i;
1398     conIndices2[i] = i;
1399   }
1400   int * tempInd = NULL;
1401   BcpsObject ** tempObj = NULL;
1402   desc->assignVars(0, tempInd,
1403                    0, tempObj,
1404                    false, numCols, varIndices1, vlhe, /*Var hard lb*/
1405                    false, numCols, varIndices2, vuhe, /*Var hard ub*/
1406                    false, 0, varIndices3, vlse,       /*Var soft lb*/
1407                    false, 0, varIndices4, vuse);      /*Var soft ub*/
1408   desc->assignCons(0, tempInd,
1409                    0, tempObj,
1410                    false, numRows, conIndices1, clhe, /*Con hard lb*/
1411                    false, numRows, conIndices2, cuhe, /*Con hard ub*/
1412                    false, 0,conIndices3,clse,         /*Con soft lb*/
1413                    false, 0,conIndices4,cuse);        /*Con soft ub*/
1414   //-------------------------------------------------------------
1415   // Mark it as an explicit node.
1416   //-------------------------------------------------------------
1417   root->setExplicit(1);
1418   return root;
1419 }
1420 
1421 
feasibleSolution(int & numInfColumns,double & colInf,int & numInfRows,double & rowInf)1422 DcoSolution * DcoModel::feasibleSolution(int & numInfColumns,
1423                                          double  & colInf,
1424                                          int & numInfRows,
1425                                          double & rowInf) {
1426   // set stats to 0
1427   numInfColumns = 0;
1428   numInfRows = 0;
1429 
1430   // for debug purposes we will keep the largest column and row infeasibility
1431   colInf = 0.0;
1432   rowInf = 0.0;
1433 
1434   // check feasibility of relxed columns, ie. integrality constraints
1435   // get vector of variables
1436   std::vector<BcpsVariable*> & cols = getVariables();
1437   for (int i=0; i<numRelaxedCols_; ++i) {
1438     // get column relaxedCol_[i]
1439     DcoVariable * curr = dynamic_cast<DcoVariable*> (cols[relaxedCols_[i]]);
1440     // check feasibility
1441     int preferredDir;
1442     double infeas = curr->infeasibility(this, preferredDir);
1443     if (infeas>0) {
1444       numInfColumns++;
1445       if (colInf<infeas) {
1446         colInf = infeas;
1447       }
1448     }
1449   }
1450 
1451   // check feasibility of relaxed rows
1452   // get vector of constraints
1453   std::vector<BcpsConstraint*> & rows = getConstraints();
1454   for (int i=0; i<numRelaxedRows_; ++i) {
1455     // get row relaxedRows_[i]
1456     DcoConstraint * curr = dynamic_cast<DcoConstraint*> (rows[relaxedRows_[i]]);
1457     // check feasibility
1458     int preferredDir;
1459     double infeas = curr->infeasibility(this, preferredDir);
1460     if (infeas>0) {
1461       numInfRows++;
1462       if (rowInf<infeas) {
1463         rowInf = infeas;
1464       }
1465     }
1466   }
1467   // report largest column and row infeasibilities
1468   dcoMessageHandler_->message(DISCO_INFEAS_REPORT, *dcoMessages_)
1469     << broker()->getProcRank()
1470     << colInf
1471     << rowInf
1472     << CoinMessageEol;
1473 
1474   // create DcoSolution instance if feasbile
1475   DcoSolution * dco_sol = 0;
1476   if (numInfColumns==0 && numInfRows==0) {
1477     double const * sol = solver()->getColSolution();
1478     double quality = solver()->getObjValue();
1479     dco_sol = new DcoSolution(numCols_, sol, quality);
1480     dco_sol->setBroker(broker_);
1481 
1482     // log debug information
1483     dcoMessageHandler_->message(DISCO_SOL_FOUND, *dcoMessages_)
1484       << broker()->getProcRank()
1485       << quality
1486       << CoinMessageEol;
1487   }
1488   return dco_sol;
1489 }
1490 
1491 //todo(aykut) When all node bounds are worse than incumbent solution
1492 // this function reports negative gap.
1493 // this happens since Alps takes nodes that will be fathomed into account,
1494 // in (getBestNode function).
nodeLog(AlpsTreeNode * node,bool force)1495 void DcoModel::nodeLog(AlpsTreeNode * node, bool force) {
1496   if ((broker_->getProcType() != AlpsProcessTypeMaster) &&
1497       (broker_->getProcType() != AlpsProcessTypeSerial)) {
1498     // do nothing if not serial code nor master processor
1499     return;
1500   }
1501   // todo(aykut) somehow AlpsKnowledgeBrokerMPI does not call this
1502   // function properly. It calls this once and only header is printed.
1503   // Check this. Disable in parallel mode.
1504   if (broker_->getProcType() == AlpsProcessTypeMaster) {
1505     return ;
1506   }
1507   // number of processed nodes
1508   int num_processed = broker()->getNumNodesProcessed();
1509   // number of nodes left
1510   int num_left  = broker()->updateNumNodesLeft();
1511   // log interval
1512   int interval =
1513     broker()->getModel()->AlpsPar()->entry(AlpsParams::nodeLogInterval);
1514   // need to print header if this is the first call of this function.
1515   dcoMessageHandler_->setPrefix(0);
1516   if (broker()->getNumNodeLog()==0 or
1517       broker()->getNumNodeLog()%50==0) {
1518     //todo(aykut) fix this function's documentation in Alps. It does not
1519     // count how many times this function is called, but how many times
1520     // this function is logged.
1521     broker()->setNumNodeLog(broker()->getNumNodeLog()+1);
1522     dcoMessageHandler_->message(DISCO_NODE_LOG_HEADER,
1523                                 *dcoMessages_)
1524       << CoinMessageEol;
1525   }
1526   else if (force || (num_processed % interval == 0)) {
1527     double lb = ALPS_INFINITY;
1528     AlpsTreeNode * best_node = broker()->getBestNode();
1529     if (best_node) {
1530       lb = best_node->getQuality();
1531     }
1532     else {
1533       // no nodes stored in the broker.
1534     }
1535     if (broker()->hasKnowledge(AlpsKnowledgeTypeSolution)) {
1536       double ub = broker()->getIncumbentValue();
1537       // check whether lb is larger than ub, this might happen if there are
1538       // nodes that are not fathomed yet.
1539       if (lb>ub) {
1540         lb = ub;
1541       }
1542       std::stringstream lb_ss;
1543       lb_ss << std::setw(14) << std::left << std::scientific << lb;
1544       broker()->setNumNodeLog(broker()->getNumNodeLog()+1);
1545       double gap = 100.0*((ub-lb)/fabs(ub));
1546       // we need to convert lb, ub, gap into strings, CoinMessageHandler is
1547       // not capable for this yet.
1548       std::stringstream ub_ss;
1549       ub_ss << std::setw(14) << std::left << std::scientific << ub;
1550       std::stringstream gap_ss;
1551       gap_ss.precision(2);
1552       gap_ss << std::setw(6) << std::fixed << std::left << gap;
1553       dcoMessageHandler_->message(DISCO_NODE_LOG, *dcoMessages_)
1554         << num_processed
1555         << num_left
1556         << lb_ss.str().c_str()
1557         << ub_ss.str().c_str()
1558         << gap_ss.str().c_str()
1559         << static_cast<int>(broker()->timer().getCpuTime())
1560         << CoinMessageEol;
1561     }
1562     else {
1563       std::stringstream lb_ss;
1564       lb_ss << std::setw(14) << std::left << std::scientific << lb;
1565       broker()->setNumNodeLog(broker()->getNumNodeLog()+1);
1566       dcoMessageHandler_->message(DISCO_NODE_LOG_NO_SOL, *dcoMessages_)
1567         << num_processed
1568         << num_left
1569         << lb_ss.str().c_str()
1570         << static_cast<int>(broker()->timer().getCpuTime())
1571         << CoinMessageEol;
1572     }
1573   }
1574   dcoMessageHandler_->setPrefix(1);
1575 }
1576 
1577 /// This is called at the end of the AlpsKnowledgeBroker::rootSearch
1578 /// Prints solution statistics
modelLog()1579 void DcoModel::modelLog() {
1580   if (broker_->getProcType() == AlpsProcessTypeSerial) {
1581 #if defined(__OA__) || defined(__COLA__)
1582     // report solver iterations
1583 
1584     // notes(aykut) This will overflow if number of iterations is large, but I
1585     // have no option since CoinMessageHandler does not overload operator <<
1586     // for long long int.
1587     int num_iter = static_cast<int> (numRelaxIterations_);
1588     dcoMessageHandler_->message(DISCO_SOLVER_ITERATIONS, *dcoMessages_)
1589       << num_iter
1590       << CoinMessageEol;
1591 #endif
1592     // report cut generator statistics
1593     std::map<DcoConstraintType, DcoConGenerator*>::iterator it;
1594     for (it=conGenerators_.begin(); it != conGenerators_.end(); ++it) {
1595       DcoConGenerator * curr = it->second;
1596       if (curr->stats().numCalls() > 0) {
1597         dcoMessageHandler_->message(DISCO_CUT_STATS_FINAL,
1598                                         *dcoMessages_)
1599           << curr->name()
1600           << curr->stats().numCalls()
1601           << curr->stats().numConsGenerated()
1602           << curr->stats().numConsUsed()
1603           << curr->stats().time()
1604           << curr->strategy()
1605           << CoinMessageEol;
1606       }
1607     }
1608     // report heuristic statistics
1609     for (unsigned int k=0; k<heuristics_.size(); ++k) {
1610       if (heuristics(k)->stats().numCalls() > 0) {
1611         dcoMessageHandler_->message(DISCO_HEUR_STATS_FINAL,
1612                                     *dcoMessages_)
1613           << heuristics(k)->name()
1614           << heuristics(k)->stats().numCalls()
1615           << heuristics(k)->stats().numSolutions()
1616           << heuristics(k)->stats().time()
1617           << heuristics(k)->strategy()
1618           << CoinMessageEol;
1619       }
1620     }
1621   }
1622   else if (broker_->getProcType()==AlpsProcessTypeMaster) {
1623     dcoMessageHandler_->message(0, "Dco",
1624                                 "Don't know how to log in parallel mode.",
1625                                 'G', DISCO_DLOG_MPI)
1626       << CoinMessageEol;
1627   }
1628 }
1629 
reportFeasibility()1630 void DcoModel::reportFeasibility() {
1631   // return if there is no solution to report
1632   if (broker_->getNumKnowledges(AlpsKnowledgeTypeSolution)==0) {
1633     return;
1634   }
1635   if (broker_->getProcType()!=AlpsProcessTypeSerial and
1636       broker_->getProcType()!=AlpsProcessTypeMaster) {
1637     return;
1638   }
1639   // get solution
1640   DcoSolution * dcosol =
1641     dynamic_cast<DcoSolution*>
1642     (broker_->getBestKnowledge(AlpsKnowledgeTypeSolution).first);
1643   // conic constraints
1644   double const * sol = dcosol->getValues();
1645 
1646   // integrality
1647   dcoMessageHandler_->message(0, "Dco", "Integrality Report",
1648                               'G', DISCO_DLOG_PROCESS)
1649     << CoinMessageEol;
1650   std::stringstream msg;
1651   for (int i=0; i<numRelaxedCols_; ++i) {
1652     msg << "x["
1653         << relaxedCols_[i]
1654         << "] = "
1655         << sol[relaxedCols_[i]];
1656     dcoMessageHandler_->message(0, "Dco", msg.str().c_str(),
1657                                 'G', DISCO_DLOG_PROCESS)
1658       << CoinMessageEol;
1659     msg.str(std::string());
1660   }
1661 
1662   // conic constraints
1663   dcoMessageHandler_->message(0, "Dco", "Conic Feasibility",
1664                               'G', DISCO_DLOG_PROCESS)
1665     << CoinMessageEol;
1666   for (int i=numLinearRows_; i<numLinearRows_+numConicRows_; ++i) {
1667     DcoConicConstraint * con =
1668       dynamic_cast<DcoConicConstraint*> (constraints_[i]);
1669     int const * members = con->coneMembers();
1670     DcoLorentzConeType type = con->coneType();
1671     int size = con->coneSize();
1672 
1673     double * values = new double[size];
1674     for (int j=0; j<size; ++j) {
1675       values[j] = sol[members[j]];
1676     }
1677     double term1 = 0.0;
1678     double term2 = 0.0;
1679     if (type==DcoLorentzCone) {
1680       term1 = values[0];
1681       term2 = std::inner_product(values+1, values+size, values+1, 0.0);
1682       term2 = sqrt(term2);
1683     }
1684     else if (type==DcoRotatedLorentzCone) {
1685       term1 = 2.0*sol[members[0]]*sol[members[1]];
1686       term2 = std::inner_product(values+2, values+size, values+2, 0.0);
1687     }
1688     else {
1689       dcoMessageHandler_->message(DISCO_UNKNOWN_CONETYPE,
1690                                   *dcoMessages_)
1691         << type << CoinMessageEol;
1692     }
1693     msg << "Cone "
1694         << i - numLinearRows_
1695         << " "
1696         << term1-term2;
1697     dcoMessageHandler_->message(0, "Dco", msg.str().c_str(),
1698                                 'G', DISCO_DLOG_PROCESS)
1699       << CoinMessageEol;
1700     msg.str(std::string());
1701     delete[] values;
1702   }
1703 
1704   // print maximimum violation for integrality
1705   {
1706     double max_violation = 0.0;
1707     for (int i=0; i<numRelaxedCols_; ++i) {
1708       double value = sol[relaxedCols_[i]];
1709       double dist_to_floor = value - floor(value);
1710       double dist_to_ceil = ceil(value) - value;
1711       double viol = std::min(dist_to_floor, dist_to_ceil);
1712       if (viol > max_violation) {
1713         max_violation = viol;
1714       }
1715     }
1716     dcoMessageHandler_->message(DISCO_SOL_INT_FEAS_REPORT,
1717                                 *dcoMessages_)
1718       << max_violation << CoinMessageEol;
1719   }
1720 
1721   // print maximimum violation for conic constraints
1722   {
1723     double max_violation = 0.0;
1724     for (int i=numLinearRows_; i<numLinearRows_+numConicRows_; ++i) {
1725       DcoConicConstraint * con =
1726         dynamic_cast<DcoConicConstraint*> (constraints_[i]);
1727       int const * members = con->coneMembers();
1728       DcoLorentzConeType type = con->coneType();
1729       int size = con->coneSize();
1730       double * values = new double[size];
1731       for (int j=0; j<size; ++j) {
1732         values[j] = sol[members[j]];
1733       }
1734       double term1 = 0.0;
1735       double term2 = 0.0;
1736       if (type==DcoLorentzCone) {
1737         term1 = values[0];
1738         term2 = std::inner_product(values+1, values+size, values+1, 0.0);
1739         term2 = sqrt(term2);
1740       }
1741       else if (type==DcoRotatedLorentzCone) {
1742         term1 = 2.0*sol[members[0]]*sol[members[1]];
1743         term2 = std::inner_product(values+2, values+size, values+2, 0.0);
1744       }
1745       else {
1746         dcoMessageHandler_->message(DISCO_UNKNOWN_CONETYPE,
1747                                     *dcoMessages_)
1748           << type << CoinMessageEol;
1749       }
1750       double viol = term2 - term1;
1751       if (viol > max_violation) {
1752         max_violation = viol;
1753       }
1754       delete[] values;
1755     }
1756     dcoMessageHandler_->message(DISCO_SOL_CONE_FEAS_REPORT,
1757                                 *dcoMessages_)
1758       << max_violation << CoinMessageEol;
1759   }
1760 }
1761 
1762 /// The method that encodes this instance of model into the given
1763 /// #AlpsEncoded object.
encode(AlpsEncoded * encoded) const1764 AlpsReturnStatus DcoModel::encode(AlpsEncoded * encoded) const {
1765   AlpsReturnStatus status;
1766   // encode Alps parts
1767   status = AlpsModel::encode(encoded);
1768   if (status!=AlpsReturnStatusOk) {
1769     dcoMessageHandler_->message(DISCO_UNEXPECTED_ENCODE_STATUS, *dcoMessages_)
1770       << __FILE__ << __LINE__ << CoinMessageEol;
1771   }
1772   // encode problem name
1773   int probname_size = static_cast<int>(problemName_.size());
1774   encoded->writeRep(probname_size);
1775   encoded->writeRep(problemName_.c_str(), probname_size);
1776   // encode number of constraints
1777   encoded->writeRep(numCols_);
1778   encoded->writeRep(colLB_, numCols_);
1779   encoded->writeRep(colUB_, numCols_);
1780   encoded->writeRep(numLinearRows_);
1781   encoded->writeRep(numConicRows_);
1782   encoded->writeRep(rowLB_, numRows_);
1783   encoded->writeRep(rowUB_, numRows_);
1784   encoded->writeRep(objSense_);
1785   encoded->writeRep(objCoef_, numCols_);
1786   encoded->writeRep(numIntegerCols_);
1787   encoded->writeRep(integerCols_, numIntegerCols_);
1788   encoded->writeRep(isInteger_, numCols_);
1789   // encode cone info
1790   if (numConicRows_) {
1791     encoded->writeRep(coneStart_, numConicRows_+1);
1792     encoded->writeRep(coneType_, numConicRows_);
1793     encoded->writeRep(coneMembers_, coneStart_[numConicRows_]);
1794   }
1795   // encode matrix
1796   encoded->writeRep(matrix_->getNumElements());
1797   encoded->writeRep(matrix_->getVectorStarts(), numLinearRows_);
1798   encoded->writeRep(matrix_->getVectorLengths(), numLinearRows_);
1799   encoded->writeRep(matrix_->getIndices(), matrix_->getNumElements());
1800   encoded->writeRep(matrix_->getElements(), matrix_->getNumElements());
1801   encoded->writeRep(initOAcuts_);
1802   // encode parameters
1803   dcoPar_->pack(*encoded);
1804 
1805   // debug stuff
1806   std::stringstream debug_msg;
1807   debug_msg << "Proc[" << broker_->getProcRank() << "]"
1808             << " model " << this << " encoded." << std::endl;
1809   dcoMessageHandler_->message(0, "Dco", debug_msg.str().c_str(),
1810                               'G', DISCO_DLOG_MPI)
1811     << CoinMessageEol;
1812   // end of debug stuff
1813 
1814   return status;
1815 }
1816 
1817 /// The method that decodes the given #AlpsEncoded object into a new #DcoModel
1818 /// instance and returns a pointer to it.
decode(AlpsEncoded & encoded) const1819 AlpsKnowledge * DcoModel::decode(AlpsEncoded & encoded) const {
1820   std::cerr << "not implemented yet." << std::endl;
1821   return NULL;
1822 }
1823 
decodeToSelf(AlpsEncoded & encoded)1824 AlpsReturnStatus DcoModel::decodeToSelf(AlpsEncoded & encoded) {
1825   AlpsReturnStatus status;
1826   // decode Alps parts
1827   status = AlpsModel::decodeToSelf(encoded);
1828   if (status!=AlpsReturnStatusOk) {
1829     dcoMessageHandler_->message(DISCO_UNEXPECTED_DECODE_STATUS, *dcoMessages_)
1830       << __FILE__ << __LINE__ << CoinMessageEol;
1831   }
1832   // decode problem name
1833   int probname_size = 0;
1834   char * probname = NULL;
1835   encoded.readRep(probname_size);
1836   encoded.readRep(probname, probname_size);
1837   problemName_ = probname;
1838   delete[] probname;
1839   // decode rest
1840   encoded.readRep(numCols_);
1841   encoded.readRep(colLB_, numCols_);
1842   encoded.readRep(colUB_, numCols_);
1843   encoded.readRep(numLinearRows_);
1844   encoded.readRep(numConicRows_);
1845   numRows_ = numLinearRows_+numConicRows_;
1846   encoded.readRep(rowLB_, numRows_);
1847   encoded.readRep(rowUB_, numRows_);
1848   encoded.readRep(objSense_);
1849   encoded.readRep(objCoef_, numCols_);
1850   encoded.readRep(numIntegerCols_);
1851   encoded.readRep(integerCols_, numIntegerCols_);
1852   encoded.readRep(isInteger_, numCols_);
1853   // decode cone info
1854   if (numConicRows_) {
1855     int cone_start_size;
1856     encoded.readRep(coneStart_, cone_start_size);
1857     assert(cone_start_size==numConicRows_+1);
1858     encoded.readRep(coneType_, numConicRows_);
1859     encoded.readRep(coneMembers_, coneStart_[numConicRows_]);
1860   }
1861   // decode matrix
1862   int num_elem;
1863   int * starts;
1864   int * lengths;
1865   int * indices;
1866   double * elements;
1867   encoded.readRep(num_elem);
1868   encoded.readRep(starts, numLinearRows_);
1869   encoded.readRep(lengths, numLinearRows_);
1870   encoded.readRep(indices, num_elem);
1871   encoded.readRep(elements, num_elem);
1872   matrix_ = new CoinPackedMatrix(false, numCols_, numLinearRows_, num_elem,
1873                                  elements, indices, starts, lengths, 0.0, 0.0);
1874   encoded.readRep(initOAcuts_);
1875   delete[] starts;
1876   delete[] lengths;
1877   delete[] indices;
1878   delete[] elements;
1879   dcoPar_->unpack(encoded);
1880 
1881   // debug stuff
1882   std::stringstream debug_msg;
1883   debug_msg << "Proc[" << broker_->getProcRank() << "]"
1884             << " model decoded into " << this << "." << std::endl;
1885   dcoMessageHandler_->message(0, "Dco", debug_msg.str().c_str(),
1886                               'G', DISCO_DLOG_MPI)
1887     << CoinMessageEol;
1888   // end of debug stuff
1889 
1890   return status;
1891 }
1892 
addNumRelaxIterations()1893 void DcoModel::addNumRelaxIterations() {
1894   numRelaxIterations_ += solver_->getIterationCount();
1895 }
1896