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