1/* $Id: ClpNode.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2008, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include "ClpSimplex.hpp"
8#include "ClpNode.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpDualRowSteepest.hpp"
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpNode::ClpNode () :
20 branchingValue_(0.5),
21 objectiveValue_(0.0),
22 sumInfeasibilities_(0.0),
23 estimatedSolution_(0.0),
24 factorization_(NULL),
25 weights_(NULL),
26 status_(NULL),
27 primalSolution_(NULL),
28 dualSolution_(NULL),
29 lower_(NULL),
30 upper_(NULL),
31 pivotVariables_(NULL),
32 fixed_(NULL),
33 sequence_(1),
34 numberInfeasibilities_(0),
35 depth_(0),
36 numberFixed_(0),
37 flags_(0),
38 maximumFixed_(0),
39 maximumRows_(0),
40 maximumColumns_(0),
41 maximumIntegers_(0)
42{
43 branchState_.firstBranch = 0;
44 branchState_.branch = 0;
45}
46//-------------------------------------------------------------------
47// Useful Constructor from model
48//-------------------------------------------------------------------
49ClpNode::ClpNode (ClpSimplex * model, const ClpNodeStuff * stuff, int depth) :
50 branchingValue_(0.5),
51 objectiveValue_(0.0),
52 sumInfeasibilities_(0.0),
53 estimatedSolution_(0.0),
54 factorization_(NULL),
55 weights_(NULL),
56 status_(NULL),
57 primalSolution_(NULL),
58 dualSolution_(NULL),
59 lower_(NULL),
60 upper_(NULL),
61 pivotVariables_(NULL),
62 fixed_(NULL),
63 sequence_(1),
64 numberInfeasibilities_(0),
65 depth_(0),
66 numberFixed_(0),
67 flags_(0),
68 maximumFixed_(0),
69 maximumRows_(0),
70 maximumColumns_(0),
71 maximumIntegers_(0)
72{
73 branchState_.firstBranch = 0;
74 branchState_.branch = 0;
75 gutsOfConstructor(model, stuff, 0, depth);
76}
77
78//-------------------------------------------------------------------
79// Most of work of constructor from model
80//-------------------------------------------------------------------
81void
82ClpNode::gutsOfConstructor (ClpSimplex * model, const ClpNodeStuff * stuff,
83 int arraysExist, int depth)
84{
85 int numberRows = model->numberRows();
86 int numberColumns = model->numberColumns();
87 int numberTotal = numberRows + numberColumns;
88 int maximumTotal = maximumRows_ + maximumColumns_;
89 depth_ = depth;
90 // save stuff
91 objectiveValue_ = model->objectiveValue() * model->optimizationDirection();
92 estimatedSolution_ = objectiveValue_;
93 flags_ = 1; //say scaled
94 if (!arraysExist) {
95 maximumRows_ = CoinMax(maximumRows_, numberRows);
96 maximumColumns_ = CoinMax(maximumColumns_, numberColumns);
97 maximumTotal = maximumRows_ + maximumColumns_;
98 assert (!factorization_);
99 factorization_ = new ClpFactorization(*model->factorization(), numberRows);
100 status_ = CoinCopyOfArrayPartial(model->statusArray(), maximumTotal, numberTotal);
101 primalSolution_ = CoinCopyOfArrayPartial(model->solutionRegion(), maximumTotal, numberTotal);
102 dualSolution_ = CoinCopyOfArrayPartial(model->djRegion(), maximumTotal, numberTotal); //? has duals as well?
103 pivotVariables_ = CoinCopyOfArrayPartial(model->pivotVariable(), maximumRows_, numberRows);
104 ClpDualRowSteepest* pivot =
105 dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
106 if (pivot) {
107 assert (!weights_);
108 weights_ = new ClpDualRowSteepest(*pivot);
109 }
110 } else {
111 if (arraysExist == 2)
112 assert(lower_);
113 if (numberRows <= maximumRows_ && numberColumns <= maximumColumns_) {
114 CoinMemcpyN(model->statusArray(), numberTotal, status_);
115 if (arraysExist == 1) {
116 *factorization_ = *model->factorization();
117 CoinMemcpyN(model->solutionRegion(), numberTotal, primalSolution_);
118 CoinMemcpyN(model->djRegion(), numberTotal, dualSolution_); //? has duals as well?
119 ClpDualRowSteepest* pivot =
120 dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
121 if (pivot) {
122 if (weights_) {
123 //if (weights_->numberRows()==pivot->numberRows()) {
124 weights_->fill(*pivot);
125 //} else {
126 //delete weights_;
127 //weights_ = new ClpDualRowSteepest(*pivot);
128 //}
129 } else {
130 weights_ = new ClpDualRowSteepest(*pivot);
131 }
132 }
133 CoinMemcpyN(model->pivotVariable(), numberRows, pivotVariables_);
134 } else {
135 CoinMemcpyN(model->primalColumnSolution(), numberColumns, primalSolution_);
136 CoinMemcpyN(model->dualColumnSolution(), numberColumns, dualSolution_);
137 flags_ = 0;
138 CoinMemcpyN(model->dualRowSolution(), numberRows, dualSolution_ + numberColumns);
139 }
140 } else {
141 // size has changed
142 maximumRows_ = CoinMax(maximumRows_, numberRows);
143 maximumColumns_ = CoinMax(maximumColumns_, numberColumns);
144 maximumTotal = maximumRows_ + maximumColumns_;
145 delete weights_;
146 weights_ = NULL;
147 delete [] status_;
148 delete [] primalSolution_;
149 delete [] dualSolution_;
150 delete [] pivotVariables_;
151 status_ = CoinCopyOfArrayPartial(model->statusArray(), maximumTotal, numberTotal);
152 primalSolution_ = new double [maximumTotal*sizeof(double)];
153 dualSolution_ = new double [maximumTotal*sizeof(double)];
154 if (arraysExist == 1) {
155 *factorization_ = *model->factorization(); // I think this is OK
156 CoinMemcpyN(model->solutionRegion(), numberTotal, primalSolution_);
157 CoinMemcpyN(model->djRegion(), numberTotal, dualSolution_); //? has duals as well?
158 ClpDualRowSteepest* pivot =
159 dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
160 if (pivot) {
161 assert (!weights_);
162 weights_ = new ClpDualRowSteepest(*pivot);
163 }
164 } else {
165 CoinMemcpyN(model->primalColumnSolution(), numberColumns, primalSolution_);
166 CoinMemcpyN(model->dualColumnSolution(), numberColumns, dualSolution_);
167 flags_ = 0;
168 CoinMemcpyN(model->dualRowSolution(), numberRows, dualSolution_ + numberColumns);
169 }
170 pivotVariables_ = new int [maximumRows_];
171 if (model->pivotVariable() && model->numberRows() == numberRows)
172 CoinMemcpyN(model->pivotVariable(), numberRows, pivotVariables_);
173 else
174 CoinFillN(pivotVariables_, numberRows, -1);
175 }
176 }
177 numberFixed_ = 0;
178 const double * lower = model->columnLower();
179 const double * upper = model->columnUpper();
180 const double * solution = model->primalColumnSolution();
181 const char * integerType = model->integerInformation();
182 const double * columnScale = model->columnScale();
183 if (!flags_)
184 columnScale = NULL; // as duals correct
185 int iColumn;
186 sequence_ = -1;
187 double integerTolerance = stuff->integerTolerance_;
188 double mostAway = 0.0;
189 int bestPriority = COIN_INT_MAX;
190 sumInfeasibilities_ = 0.0;
191 numberInfeasibilities_ = 0;
192 int nFix = 0;
193 double gap = CoinMax(model->dualObjectiveLimit() - objectiveValue_, 1.0e-4);
194#define PSEUDO 3
195#if PSEUDO==1||PSEUDO==2
196 // Column copy of matrix
197 ClpPackedMatrix * matrix = model->clpScaledMatrix();
198 const double *objective = model->costRegion() ;
199 if (!objective) {
200 objective = model->objective();
201 //if (!matrix)
202 matrix = dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
203 } else if (!matrix) {
204 matrix = dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
205 }
206 const double * element = matrix->getElements();
207 const int * row = matrix->getIndices();
208 const CoinBigIndex * columnStart = matrix->getVectorStarts();
209 const int * columnLength = matrix->getVectorLengths();
210 double direction = model->optimizationDirection();
211 const double * dual = dualSolution_ + numberColumns;
212#if PSEUDO==2
213 double * activeWeight = new double [numberRows];
214 const double * rowLower = model->rowLower();
215 const double * rowUpper = model->rowUpper();
216 const double * rowActivity = model->primalRowSolution();
217 double tolerance = 1.0e-6;
218 for (int iRow = 0; iRow < numberRows; iRow++) {
219 // could use pi to see if active or activity
220 if (rowActivity[iRow] > rowUpper[iRow] - tolerance
221 || rowActivity[iRow] < rowLower[iRow] + tolerance) {
222 activeWeight[iRow] = 0.0;
223 } else {
224 activeWeight[iRow] = -1.0;
225 }
226 }
227 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
228 if (integerType[iColumn]) {
229 double value = solution[iColumn];
230 if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
231 CoinBigIndex start = columnStart[iColumn];
232 CoinBigIndex end = start + columnLength[iColumn];
233 for (CoinBigIndex j = start; j < end; j++) {
234 int iRow = row[j];
235 if (activeWeight[iRow] >= 0.0)
236 activeWeight[iRow] += 1.0;
237 }
238 }
239 }
240 }
241 for (int iRow = 0; iRow < numberRows; iRow++) {
242 if (activeWeight[iRow] > 0.0) {
243 // could use pi
244 activeWeight[iRow] = 1.0 / activeWeight[iRow];
245 } else {
246 activeWeight[iRow] = 0.0;
247 }
248 }
249#endif
250#endif
251 const double * downPseudo = stuff->downPseudo_;
252 const int * numberDown = stuff->numberDown_;
253 const int * numberDownInfeasible = stuff->numberDownInfeasible_;
254 const double * upPseudo = stuff->upPseudo_;
255 const int * priority = stuff->priority_;
256 const int * numberUp = stuff->numberUp_;
257 const int * numberUpInfeasible = stuff->numberUpInfeasible_;
258 int numberBeforeTrust = stuff->numberBeforeTrust_;
259 int stateOfSearch = stuff->stateOfSearch_;
260 int iInteger = 0;
261 // weight at 1.0 is max min (CbcBranch was 0.8,0.1) (ClpNode was 0.9,0.9)
262#define WEIGHT_AFTER 0.9
263#define WEIGHT_BEFORE 0.2
264 //Stolen from Constraint Integer Programming book (with epsilon change)
265#define WEIGHT_PRODUCT
266#ifdef WEIGHT_PRODUCT
267 double smallChange = stuff->smallChange_;
268#endif
269#ifndef INFEAS_MULTIPLIER
270#define INFEAS_MULTIPLIER 1.0
271#endif
272 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
273 if (integerType[iColumn]) {
274 double value = solution[iColumn];
275 value = CoinMax(value, static_cast<double> (lower[iColumn]));
276 value = CoinMin(value, static_cast<double> (upper[iColumn]));
277 double nearest = floor(value + 0.5);
278 if (fabs(value - nearest) > integerTolerance) {
279 numberInfeasibilities_++;
280 sumInfeasibilities_ += fabs(value - nearest);
281#if PSEUDO==1 || PSEUDO ==2
282 double upValue = 0.0;
283 double downValue = 0.0;
284 double value2 = direction * objective[iColumn];
285 //double dj2=value2;
286 if (value2) {
287 if (value2 > 0.0)
288 upValue += 1.5 * value2;
289 else
290 downValue -= 1.5 * value2;
291 }
292 CoinBigIndex start = columnStart[iColumn];
293 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
294 for (CoinBigIndex j = start; j < end; j++) {
295 int iRow = row[j];
296 value2 = -dual[iRow];
297 if (value2) {
298 value2 *= element[j];
299 //dj2 += value2;
300#if PSEUDO==2
301 assert (activeWeight[iRow] > 0.0 || fabs(dual[iRow]) < 1.0e-6);
302 value2 *= activeWeight[iRow];
303#endif
304 if (value2 > 0.0)
305 upValue += value2;
306 else
307 downValue -= value2;
308 }
309 }
310 //assert (fabs(dj2)<1.0e-4);
311 int nUp = numberUp[iInteger];
312 double upValue2 = (upPseudo[iInteger] / (1.0 + nUp));
313 // Extra for infeasible branches
314 if (nUp) {
315 double ratio = 1.0 + INFEAS_MULTIPLIER*static_cast<double>(numberUpInfeasible[iInteger]) /
316 static_cast<double>(nUp);
317 upValue2 *= ratio;
318 }
319 int nDown = numberDown[iInteger];
320 double downValue2 = (downPseudo[iInteger] / (1.0 + nDown));
321 if (nDown) {
322 double ratio = 1.0 + INFEAS_MULTIPLIER*static_cast<double>(numberDownInfeasible[iInteger]) /
323 static_cast<double>(nDown);
324 downValue2 *= ratio;
325 }
326 //printf("col %d - downPi %g up %g, downPs %g up %g\n",
327 // iColumn,upValue,downValue,upValue2,downValue2);
328 upValue = CoinMax(0.1 * upValue, upValue2);
329 downValue = CoinMax(0.1 * downValue, downValue2);
330 //upValue = CoinMax(upValue,1.0e-8);
331 //downValue = CoinMax(downValue,1.0e-8);
332 upValue *= ceil(value) - value;
333 downValue *= value - floor(value);
334 double infeasibility;
335 //if (depth>1000)
336 //infeasibility = CoinMax(upValue,downValue)+integerTolerance;
337 //else
338 if (stateOfSearch <= 2) {
339 // no solution
340 infeasibility = (1.0 - WEIGHT_BEFORE) * CoinMax(upValue, downValue) +
341 WEIGHT_BEFORE * CoinMin(upValue, downValue) + integerTolerance;
342 } else {
343#ifndef WEIGHT_PRODUCT
344 infeasibility = (1.0 - WEIGHT_AFTER) * CoinMax(upValue, downValue) +
345 WEIGHT_AFTER * CoinMin(upValue, downValue) + integerTolerance;
346#else
347 infeasibility = CoinMax(CoinMax(upValue, downValue), smallChange) *
348 CoinMax(CoinMin(upValue, downValue), smallChange);
349#endif
350 }
351 estimatedSolution_ += CoinMin(upValue2, downValue2);
352#elif PSEUDO==3
353 int nUp = numberUp[iInteger];
354 int nDown = numberDown[iInteger];
355 // Extra 100% for infeasible branches
356 double upValue = (ceil(value) - value) * (upPseudo[iInteger] /
357 (1.0 + nUp));
358 if (nUp) {
359 double ratio = 1.0 + INFEAS_MULTIPLIER*static_cast<double>(numberUpInfeasible[iInteger]) /
360 static_cast<double>(nUp);
361 upValue *= ratio;
362 }
363 double downValue = (value - floor(value)) * (downPseudo[iInteger] /
364 (1.0 + nDown));
365 if (nDown) {
366 double ratio = 1.0 + INFEAS_MULTIPLIER*static_cast<double>(numberDownInfeasible[iInteger]) /
367 static_cast<double>(nDown);
368 downValue *= ratio;
369 }
370 if (nUp < numberBeforeTrust || nDown < numberBeforeTrust) {
371 upValue *= 10.0;
372 downValue *= 10.0;
373 }
374
375 double infeasibility;
376 //if (depth>1000)
377 //infeasibility = CoinMax(upValue,downValue)+integerTolerance;
378 //else
379 if (stateOfSearch <= 2) {
380 // no solution
381 infeasibility = (1.0 - WEIGHT_BEFORE) * CoinMax(upValue, downValue) +
382 WEIGHT_BEFORE * CoinMin(upValue, downValue) + integerTolerance;
383 } else {
384#ifndef WEIGHT_PRODUCT
385 infeasibility = (1.0 - WEIGHT_AFTER) * CoinMax(upValue, downValue) +
386 WEIGHT_AFTER * CoinMin(upValue, downValue) + integerTolerance;
387#else
388 infeasibility = CoinMax(CoinMax(upValue, downValue), smallChange) *
389 CoinMax(CoinMin(upValue, downValue), smallChange);
390 //infeasibility += CoinMin(upValue,downValue)*smallChange;
391#endif
392 }
393 //infeasibility = 0.1*CoinMax(upValue,downValue)+
394 //0.9*CoinMin(upValue,downValue) + integerTolerance;
395 estimatedSolution_ += CoinMin(upValue, downValue);
396#else
397 double infeasibility = fabs(value - nearest);
398#endif
399 assert (infeasibility > 0.0);
400 if (priority[iInteger] < bestPriority) {
401 mostAway = 0.0;
402 bestPriority = priority[iInteger];
403 } else if (priority[iInteger] > bestPriority) {
404 infeasibility = 0.0;
405 }
406 if (infeasibility > mostAway) {
407 mostAway = infeasibility;
408 sequence_ = iColumn;
409 branchingValue_ = value;
410 branchState_.branch = 0;
411#if PSEUDO>0
412 if (upValue <= downValue)
413 branchState_.firstBranch = 1; // up
414 else
415 branchState_.firstBranch = 0; // down
416#else
417 if (value <= nearest)
418 branchState_.firstBranch = 1; // up
419 else
420 branchState_.firstBranch = 0; // down
421#endif
422 }
423 } else if (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound) {
424 bool fix = false;
425 if (columnScale) {
426 if (dualSolution_[iColumn] > gap * columnScale[iColumn])
427 fix = true;
428 } else {
429 if (dualSolution_[iColumn] > gap)
430 fix = true;
431 }
432 if (fix) {
433 nFix++;
434 //printf("fixed %d to zero gap %g dj %g %g\n",iColumn,
435 // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
436 model->setColumnStatus(iColumn, ClpSimplex::isFixed);
437 }
438 } else if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
439 bool fix = false;
440 if (columnScale) {
441 if (-dualSolution_[iColumn] > gap * columnScale[iColumn])
442 fix = true;
443 } else {
444 if (-dualSolution_[iColumn] > gap)
445 fix = true;
446 }
447 if (fix) {
448 nFix++;
449 //printf("fixed %d to one gap %g dj %g %g\n",iColumn,
450 // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
451 model->setColumnStatus(iColumn, ClpSimplex::isFixed);
452 }
453 }
454 iInteger++;
455 }
456 }
457 //printf("Choosing %d inf %g pri %d\n",
458 // sequence_,mostAway,bestPriority);
459#if PSEUDO == 2
460 delete [] activeWeight;
461#endif
462 if (lower_) {
463 // save bounds
464 if (iInteger > maximumIntegers_) {
465 delete [] lower_;
466 delete [] upper_;
467 maximumIntegers_ = iInteger;
468 lower_ = new int [maximumIntegers_];
469 upper_ = new int [maximumIntegers_];
470 }
471 iInteger = 0;
472 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
473 if (integerType[iColumn]) {
474 lower_[iInteger] = static_cast<int> (lower[iColumn]);
475 upper_[iInteger] = static_cast<int> (upper[iColumn]);
476 iInteger++;
477 }
478 }
479 }
480 // Could omit save of fixed if doing full save of bounds
481 if (sequence_ >= 0 && nFix) {
482 if (nFix > maximumFixed_) {
483 delete [] fixed_;
484 fixed_ = new int [nFix];
485 maximumFixed_ = nFix;
486 }
487 numberFixed_ = 0;
488 unsigned char * status = model->statusArray();
489 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
490 if (status[iColumn] != status_[iColumn]) {
491 if (solution[iColumn] <= lower[iColumn] + 2.0 * integerTolerance) {
492 model->setColumnUpper(iColumn, lower[iColumn]);
493 fixed_[numberFixed_++] = iColumn;
494 } else {
495 assert (solution[iColumn] >= upper[iColumn] - 2.0 * integerTolerance);
496 model->setColumnLower(iColumn, upper[iColumn]);
497 fixed_[numberFixed_++] = iColumn | 0x10000000;
498 }
499 }
500 }
501 //printf("%d fixed\n",numberFixed_);
502 }
503}
504
505//-------------------------------------------------------------------
506// Copy constructor
507//-------------------------------------------------------------------
508ClpNode::ClpNode (const ClpNode & )
509{
510 printf("ClpNode copy not implemented\n");
511 abort();
512}
513
514//-------------------------------------------------------------------
515// Destructor
516//-------------------------------------------------------------------
517ClpNode::~ClpNode ()
518{
519 delete factorization_;
520 delete weights_;
521 delete [] status_;
522 delete [] primalSolution_;
523 delete [] dualSolution_;
524 delete [] lower_;
525 delete [] upper_;
526 delete [] pivotVariables_;
527 delete [] fixed_;
528}
529
530//----------------------------------------------------------------
531// Assignment operator
532//-------------------------------------------------------------------
533ClpNode &
534ClpNode::operator=(const ClpNode& rhs)
535{
536 if (this != &rhs) {
537 printf("ClpNode = not implemented\n");
538 abort();
539 }
540 return *this;
541}
542// Create odd arrays
543void
544ClpNode::createArrays(ClpSimplex * model)
545{
546 int numberColumns = model->numberColumns();
547 const char * integerType = model->integerInformation();
548 int iColumn;
549 int numberIntegers = 0;
550 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
551 if (integerType[iColumn])
552 numberIntegers++;
553 }
554 if (numberIntegers > maximumIntegers_ || !lower_) {
555 delete [] lower_;
556 delete [] upper_;
557 maximumIntegers_ = numberIntegers;
558 lower_ = new int [numberIntegers];
559 upper_ = new int [numberIntegers];
560 }
561}
562// Clean up as crunch is different model
563void
564ClpNode::cleanUpForCrunch()
565{
566 delete weights_;
567 weights_ = NULL;
568}
569/* Applies node to model
570 0 - just tree bounds
571 1 - tree bounds and basis etc
572 2 - saved bounds and basis etc
573*/
574void
575ClpNode::applyNode(ClpSimplex * model, int doBoundsEtc )
576{
577 int numberColumns = model->numberColumns();
578 const double * lower = model->columnLower();
579 const double * upper = model->columnUpper();
580 if (doBoundsEtc < 2) {
581 // current bound
582 int way = branchState_.firstBranch;
583 if (branchState_.branch > 0)
584 way = 1 - way;
585 if (!way) {
586 // This should also do underlying internal bound
587 model->setColumnUpper(sequence_, floor(branchingValue_));
588 } else {
589 // This should also do underlying internal bound
590 model->setColumnLower(sequence_, ceil(branchingValue_));
591 }
592 // apply dj fixings
593 for (int i = 0; i < numberFixed_; i++) {
594 int iColumn = fixed_[i];
595 if ((iColumn & 0x10000000) != 0) {
596 iColumn &= 0xfffffff;
597 model->setColumnLower(iColumn, upper[iColumn]);
598 } else {
599 model->setColumnUpper(iColumn, lower[iColumn]);
600 }
601 }
602 } else {
603 // restore bounds
604 assert (lower_);
605 int iInteger = -1;
606 const char * integerType = model->integerInformation();
607 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
608 if (integerType[iColumn]) {
609 iInteger++;
610 if (lower_[iInteger] != static_cast<int> (lower[iColumn]))
611 model->setColumnLower(iColumn, lower_[iInteger]);
612 if (upper_[iInteger] != static_cast<int> (upper[iColumn]))
613 model->setColumnUpper(iColumn, upper_[iInteger]);
614 }
615 }
616 }
617 if (doBoundsEtc && doBoundsEtc < 3) {
618 //model->copyFactorization(*factorization_);
619 model->copyFactorization(*factorization_);
620 ClpDualRowSteepest* pivot =
621 dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
622 if (pivot && weights_) {
623 pivot->fill(*weights_);
624 }
625 int numberRows = model->numberRows();
626 int numberTotal = numberRows + numberColumns;
627 CoinMemcpyN(status_, numberTotal, model->statusArray());
628 if (doBoundsEtc < 2) {
629 CoinMemcpyN(primalSolution_, numberTotal, model->solutionRegion());
630 CoinMemcpyN(dualSolution_, numberTotal, model->djRegion());
631 CoinMemcpyN(pivotVariables_, numberRows, model->pivotVariable());
632 CoinMemcpyN(dualSolution_ + numberColumns, numberRows, model->dualRowSolution());
633 } else {
634 CoinMemcpyN(primalSolution_, numberColumns, model->primalColumnSolution());
635 CoinMemcpyN(dualSolution_, numberColumns, model->dualColumnSolution());
636 CoinMemcpyN(dualSolution_ + numberColumns, numberRows, model->dualRowSolution());
637 if (model->columnScale()) {
638 // See if just primal will work
639 double * solution = model->primalColumnSolution();
640 const double * columnScale = model->columnScale();
641 int i;
642 for (i = 0; i < numberColumns; i++) {
643 solution[i] *= columnScale[i];
644 }
645 }
646 }
647 model->setObjectiveValue(objectiveValue_);
648 }
649}
650// Choose a new variable
651void
652ClpNode::chooseVariable(ClpSimplex * , ClpNodeStuff * /*info*/)
653{
654#if 0
655 int way = branchState_.firstBranch;
656 if (branchState_.branch > 0)
657 way = 1 - way;
658 assert (!branchState_.branch);
659 // We need to use pseudo costs to choose a variable
660 int numberColumns = model->numberColumns();
661#endif
662}
663// Fix on reduced costs
664int
665ClpNode::fixOnReducedCosts(ClpSimplex * )
666{
667
668 return 0;
669}
670/* Way for integer variable -1 down , +1 up */
671int
672ClpNode::way() const
673{
674 int way = branchState_.firstBranch;
675 if (branchState_.branch > 0)
676 way = 1 - way;
677 return way == 0 ? -1 : +1;
678}
679// Return true if branch exhausted
680bool
681ClpNode::fathomed() const
682{
683 return branchState_.branch >= 1
684 ;
685}
686// Change state of variable i.e. go other way
687void
688ClpNode::changeState()
689{
690 branchState_.branch++;
691 assert (branchState_.branch <= 2);
692}
693//#############################################################################
694// Constructors / Destructor / Assignment
695//#############################################################################
696
697//-------------------------------------------------------------------
698// Default Constructor
699//-------------------------------------------------------------------
700ClpNodeStuff::ClpNodeStuff () :
701 integerTolerance_(1.0e-7),
702 integerIncrement_(1.0e-8),
703 smallChange_(1.0e-8),
704 downPseudo_(NULL),
705 upPseudo_(NULL),
706 priority_(NULL),
707 numberDown_(NULL),
708 numberUp_(NULL),
709 numberDownInfeasible_(NULL),
710 numberUpInfeasible_(NULL),
711 saveCosts_(NULL),
712 nodeInfo_(NULL),
713 large_(NULL),
714 whichRow_(NULL),
715 whichColumn_(NULL),
716#ifndef NO_FATHOM_PRINT
717 handler_(NULL),
718#endif
719 nBound_(0),
720 saveOptions_(0),
721 solverOptions_(0),
722 maximumNodes_(0),
723 numberBeforeTrust_(0),
724 stateOfSearch_(0),
725 nDepth_(-1),
726 nNodes_(0),
727 numberNodesExplored_(0),
728 numberIterations_(0),
729 presolveType_(0)
730#ifndef NO_FATHOM_PRINT
731 ,startingDepth_(-1),
732 nodeCalled_(-1)
733#endif
734{
735
736}
737
738//-------------------------------------------------------------------
739// Copy constructor
740//-------------------------------------------------------------------
741ClpNodeStuff::ClpNodeStuff (const ClpNodeStuff & rhs)
742 : integerTolerance_(rhs.integerTolerance_),
743 integerIncrement_(rhs.integerIncrement_),
744 smallChange_(rhs.smallChange_),
745 downPseudo_(NULL),
746 upPseudo_(NULL),
747 priority_(NULL),
748 numberDown_(NULL),
749 numberUp_(NULL),
750 numberDownInfeasible_(NULL),
751 numberUpInfeasible_(NULL),
752 saveCosts_(NULL),
753 nodeInfo_(NULL),
754 large_(NULL),
755 whichRow_(NULL),
756 whichColumn_(NULL),
757#ifndef NO_FATHOM_PRINT
758 handler_(rhs.handler_),
759#endif
760 nBound_(0),
761 saveOptions_(rhs.saveOptions_),
762 solverOptions_(rhs.solverOptions_),
763 maximumNodes_(rhs.maximumNodes_),
764 numberBeforeTrust_(rhs.numberBeforeTrust_),
765 stateOfSearch_(rhs.stateOfSearch_),
766 nDepth_(rhs.nDepth_),
767 nNodes_(rhs.nNodes_),
768 numberNodesExplored_(rhs.numberNodesExplored_),
769 numberIterations_(rhs.numberIterations_),
770 presolveType_(rhs.presolveType_)
771#ifndef NO_FATHOM_PRINT
772 ,startingDepth_(rhs.startingDepth_),
773 nodeCalled_(rhs.nodeCalled_)
774#endif
775{
776}
777//----------------------------------------------------------------
778// Assignment operator
779//-------------------------------------------------------------------
780ClpNodeStuff &
781ClpNodeStuff::operator=(const ClpNodeStuff& rhs)
782{
783 if (this != &rhs) {
784 integerTolerance_ = rhs.integerTolerance_;
785 integerIncrement_ = rhs.integerIncrement_;
786 smallChange_ = rhs.smallChange_;
787 downPseudo_ = NULL;
788 upPseudo_ = NULL;
789 priority_ = NULL;
790 numberDown_ = NULL;
791 numberUp_ = NULL;
792 numberDownInfeasible_ = NULL;
793 numberUpInfeasible_ = NULL;
794 saveCosts_ = NULL;
795 nodeInfo_ = NULL;
796 large_ = NULL;
797 whichRow_ = NULL;
798 whichColumn_ = NULL;
799 nBound_ = 0;
800 saveOptions_ = rhs.saveOptions_;
801 solverOptions_ = rhs.solverOptions_;
802 maximumNodes_ = rhs.maximumNodes_;
803 numberBeforeTrust_ = rhs.numberBeforeTrust_;
804 stateOfSearch_ = rhs.stateOfSearch_;
805 int n = maximumNodes();
806 if (n) {
807 for (int i = 0; i < n; i++)
808 delete nodeInfo_[i];
809 }
810 delete [] nodeInfo_;
811 nodeInfo_ = NULL;
812 nDepth_ = rhs.nDepth_;
813 nNodes_ = rhs.nNodes_;
814 numberNodesExplored_ = rhs.numberNodesExplored_;
815 numberIterations_ = rhs.numberIterations_;
816 presolveType_ = rhs.presolveType_;
817#ifndef NO_FATHOM_PRINT
818 handler_ = rhs.handler_;
819 startingDepth_ = rhs.startingDepth_;
820 nodeCalled_ = rhs.nodeCalled_;
821#endif
822 }
823 return *this;
824}
825// Zaps stuff 1 - arrays, 2 ints, 3 both
826void
827ClpNodeStuff::zap(int type)
828{
829 if ((type & 1) != 0) {
830 downPseudo_ = NULL;
831 upPseudo_ = NULL;
832 priority_ = NULL;
833 numberDown_ = NULL;
834 numberUp_ = NULL;
835 numberDownInfeasible_ = NULL;
836 numberUpInfeasible_ = NULL;
837 saveCosts_ = NULL;
838 nodeInfo_ = NULL;
839 large_ = NULL;
840 whichRow_ = NULL;
841 whichColumn_ = NULL;
842 }
843 if ((type & 2) != 0) {
844 nBound_ = 0;
845 saveOptions_ = 0;
846 solverOptions_ = 0;
847 maximumNodes_ = 0;
848 numberBeforeTrust_ = 0;
849 stateOfSearch_ = 0;
850 nDepth_ = -1;
851 nNodes_ = 0;
852 presolveType_ = 0;
853 numberNodesExplored_ = 0;
854 numberIterations_ = 0;
855 }
856}
857
858//-------------------------------------------------------------------
859// Destructor
860//-------------------------------------------------------------------
861ClpNodeStuff::~ClpNodeStuff ()
862{
863 delete [] downPseudo_;
864 delete [] upPseudo_;
865 delete [] priority_;
866 delete [] numberDown_;
867 delete [] numberUp_;
868 delete [] numberDownInfeasible_;
869 delete [] numberUpInfeasible_;
870 int n = maximumNodes();
871 if (n) {
872 for (int i = 0; i < n; i++)
873 delete nodeInfo_[i];
874 }
875 delete [] nodeInfo_;
876#ifdef CLP_INVESTIGATE
877 // Should be NULL - find out why not?
878 assert (!saveCosts_);
879#endif
880 delete [] saveCosts_;
881}
882// Return maximum number of nodes
883int
884ClpNodeStuff::maximumNodes() const
885{
886 int n = 0;
887#if 0
888 if (nDepth_ != -1) {
889 if ((solverOptions_ & 32) == 0)
890 n = (1 << nDepth_);
891 else if (nDepth_)
892 n = 1;
893 }
894 assert (n == maximumNodes_ - (1 + nDepth_) || n == 0);
895#else
896 if (nDepth_ != -1) {
897 n = maximumNodes_ - (1 + nDepth_);
898 assert (n > 0);
899 }
900#endif
901 return n;
902}
903// Return maximum space for nodes
904int
905ClpNodeStuff::maximumSpace() const
906{
907 return maximumNodes_;
908}
909/* Fill with pseudocosts */
910void
911ClpNodeStuff::fillPseudoCosts(const double * down, const double * up,
912 const int * priority,
913 const int * numberDown, const int * numberUp,
914 const int * numberDownInfeasible,
915 const int * numberUpInfeasible,
916 int number)
917{
918 delete [] downPseudo_;
919 delete [] upPseudo_;
920 delete [] priority_;
921 delete [] numberDown_;
922 delete [] numberUp_;
923 delete [] numberDownInfeasible_;
924 delete [] numberUpInfeasible_;
925 downPseudo_ = CoinCopyOfArray(down, number);
926 upPseudo_ = CoinCopyOfArray(up, number);
927 priority_ = CoinCopyOfArray(priority, number);
928 numberDown_ = CoinCopyOfArray(numberDown, number);
929 numberUp_ = CoinCopyOfArray(numberUp, number);
930 numberDownInfeasible_ = CoinCopyOfArray(numberDownInfeasible, number);
931 numberUpInfeasible_ = CoinCopyOfArray(numberUpInfeasible, number);
932 // scale
933 for (int i = 0; i < number; i++) {
934 int n;
935 n = numberDown_[i];
936 if (n)
937 downPseudo_[i] *= n;
938 n = numberUp_[i];
939 if (n)
940 upPseudo_[i] *= n;
941 }
942}
943// Update pseudo costs
944void
945ClpNodeStuff::update(int way, int sequence, double change, bool feasible)
946{
947 assert (numberDown_[sequence] >= numberDownInfeasible_[sequence]);
948 assert (numberUp_[sequence] >= numberUpInfeasible_[sequence]);
949 if (way < 0) {
950 numberDown_[sequence]++;
951 if (!feasible)
952 numberDownInfeasible_[sequence]++;
953 downPseudo_[sequence] += CoinMax(change, 1.0e-12);
954 } else {
955 numberUp_[sequence]++;
956 if (!feasible)
957 numberUpInfeasible_[sequence]++;
958 upPseudo_[sequence] += CoinMax(change, 1.0e-12);
959 }
960}
961//#############################################################################
962// Constructors / Destructor / Assignment
963//#############################################################################
964
965//-------------------------------------------------------------------
966// Default Constructor
967//-------------------------------------------------------------------
968ClpHashValue::ClpHashValue () :
969 hash_(NULL),
970 numberHash_(0),
971 maxHash_(0),
972 lastUsed_(-1)
973{
974}
975//-------------------------------------------------------------------
976// Useful Constructor from model
977//-------------------------------------------------------------------
978ClpHashValue::ClpHashValue (ClpSimplex * model) :
979 hash_(NULL),
980 numberHash_(0),
981 maxHash_(0),
982 lastUsed_(-1)
983{
984 maxHash_ = 1000;
985 int numberColumns = model->numberColumns();
986 const double * columnLower = model->columnLower();
987 const double * columnUpper = model->columnUpper();
988 int numberRows = model->numberRows();
989 const double * rowLower = model->rowLower();
990 const double * rowUpper = model->rowUpper();
991 const double * objective = model->objective();
992 CoinPackedMatrix * matrix = model->matrix();
993 const int * columnLength = matrix->getVectorLengths();
994 const CoinBigIndex * columnStart = matrix->getVectorStarts();
995 const double * elementByColumn = matrix->getElements();
996 int i;
997 int ipos;
998
999 hash_ = new CoinHashLink[maxHash_];
1000
1001 for ( i = 0; i < maxHash_; i++ ) {
1002 hash_[i].value = -1.0e-100;
1003 hash_[i].index = -1;
1004 hash_[i].next = -1;
1005 }
1006 // Put in +0
1007 hash_[0].value = 0.0;
1008 hash_[0].index = 0;
1009 numberHash_ = 1;
1010 /*
1011 * Initialize the hash table. Only the index of the first value that
1012 * hashes to a value is entered in the table; subsequent values that
1013 * collide with it are not entered.
1014 */
1015 for ( i = 0; i < numberColumns; i++ ) {
1016 int length = columnLength[i];
1017 CoinBigIndex start = columnStart[i];
1018 for (CoinBigIndex i = start; i < start + length; i++) {
1019 double value = elementByColumn[i];
1020 ipos = hash ( value);
1021 if ( hash_[ipos].index == -1 ) {
1022 hash_[ipos].index = numberHash_;
1023 numberHash_++;
1024 hash_[ipos].value = elementByColumn[i];
1025 }
1026 }
1027 }
1028
1029 /*
1030 * Now take care of the values that collided in the preceding loop,
1031 * Also do other stuff
1032 */
1033 for ( i = 0; i < numberRows; i++ ) {
1034 if (numberHash_ * 2 > maxHash_)
1035 resize(true);
1036 double value;
1037 value = rowLower[i];
1038 ipos = index(value);
1039 if (ipos < 0)
1040 addValue(value);
1041 value = rowUpper[i];
1042 ipos = index(value);
1043 if (ipos < 0)
1044 addValue(value);
1045 }
1046 for ( i = 0; i < numberColumns; i++ ) {
1047 int length = columnLength[i];
1048 CoinBigIndex start = columnStart[i];
1049 if (numberHash_ * 2 > maxHash_)
1050 resize(true);
1051 double value;
1052 value = objective[i];
1053 ipos = index(value);
1054 if (ipos < 0)
1055 addValue(value);
1056 value = columnLower[i];
1057 ipos = index(value);
1058 if (ipos < 0)
1059 addValue(value);
1060 value = columnUpper[i];
1061 ipos = index(value);
1062 if (ipos < 0)
1063 addValue(value);
1064 for (CoinBigIndex j = start; j < start + length; j++) {
1065 if (numberHash_ * 2 > maxHash_)
1066 resize(true);
1067 value = elementByColumn[j];
1068 ipos = index(value);
1069 if (ipos < 0)
1070 addValue(value);
1071 }
1072 }
1073 resize(false);
1074}
1075
1076//-------------------------------------------------------------------
1077// Copy constructor
1078//-------------------------------------------------------------------
1079ClpHashValue::ClpHashValue (const ClpHashValue & rhs) :
1080 hash_(NULL),
1081 numberHash_(rhs.numberHash_),
1082 maxHash_(rhs.maxHash_),
1083 lastUsed_(rhs.lastUsed_)
1084{
1085 if (maxHash_) {
1086 CoinHashLink * newHash = new CoinHashLink[maxHash_];
1087 int i;
1088 for ( i = 0; i < maxHash_; i++ ) {
1089 newHash[i].value = rhs.hash_[i].value;
1090 newHash[i].index = rhs.hash_[i].index;
1091 newHash[i].next = rhs.hash_[i].next;
1092 }
1093 }
1094}
1095
1096//-------------------------------------------------------------------
1097// Destructor
1098//-------------------------------------------------------------------
1099ClpHashValue::~ClpHashValue ()
1100{
1101 delete [] hash_;
1102}
1103
1104//----------------------------------------------------------------
1105// Assignment operator
1106//-------------------------------------------------------------------
1107ClpHashValue &
1108ClpHashValue::operator=(const ClpHashValue& rhs)
1109{
1110 if (this != &rhs) {
1111 numberHash_ = rhs.numberHash_;
1112 maxHash_ = rhs.maxHash_;
1113 lastUsed_ = rhs.lastUsed_;
1114 delete [] hash_;
1115 if (maxHash_) {
1116 CoinHashLink * newHash = new CoinHashLink[maxHash_];
1117 int i;
1118 for ( i = 0; i < maxHash_; i++ ) {
1119 newHash[i].value = rhs.hash_[i].value;
1120 newHash[i].index = rhs.hash_[i].index;
1121 newHash[i].next = rhs.hash_[i].next;
1122 }
1123 } else {
1124 hash_ = NULL;
1125 }
1126 }
1127 return *this;
1128}
1129// Return index or -1 if not found
1130int
1131ClpHashValue::index(double value) const
1132{
1133 if (!value)
1134 return 0;
1135 int ipos = hash ( value);
1136 int returnCode = -1;
1137 while ( hash_[ipos].index >= 0 ) {
1138 if (value == hash_[ipos].value) {
1139 returnCode = hash_[ipos].index;
1140 break;
1141 } else {
1142 int k = hash_[ipos].next;
1143 if ( k == -1 ) {
1144 break;
1145 } else {
1146 ipos = k;
1147 }
1148 }
1149 }
1150 return returnCode;
1151}
1152// Add value to list and return index
1153int
1154ClpHashValue::addValue(double value)
1155{
1156 int ipos = hash ( value);
1157
1158 assert (value != hash_[ipos].value);
1159 if (hash_[ipos].index == -1) {
1160 // can put in here
1161 hash_[ipos].index = numberHash_;
1162 numberHash_++;
1163 hash_[ipos].value = value;
1164 return numberHash_ - 1;
1165 }
1166 int k = hash_[ipos].next;
1167 while (k != -1) {
1168 ipos = k;
1169 k = hash_[ipos].next;
1170 }
1171 while ( true ) {
1172 ++lastUsed_;
1173 assert (lastUsed_ <= maxHash_);
1174 if ( hash_[lastUsed_].index == -1 ) {
1175 break;
1176 }
1177 }
1178 hash_[ipos].next = lastUsed_;
1179 hash_[lastUsed_].index = numberHash_;
1180 numberHash_++;
1181 hash_[lastUsed_].value = value;
1182 return numberHash_ - 1;
1183}
1184
1185namespace {
1186 /*
1187 Originally a local static variable in ClpHashValue::hash.
1188 Local static variables are a problem when building DLLs on Windows, but
1189 file-local constants seem to be ok. -- lh, 101016 --
1190 */
1191 const int mmult_for_hash[] = {
1192 262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247,
1193 241667, 239179, 236609, 233983, 231289, 228859, 226357, 223829,
1194 221281, 218849, 216319, 213721, 211093, 208673, 206263, 203773,
1195 201233, 198637, 196159, 193603, 191161, 188701, 186149, 183761,
1196 181303, 178873, 176389, 173897, 171469, 169049, 166471, 163871,
1197 161387, 158941, 156437, 153949, 151531, 149159, 146749, 144299,
1198 141709, 139369, 136889, 134591, 132169, 129641, 127343, 124853,
1199 122477, 120163, 117757, 115361, 112979, 110567, 108179, 105727,
1200 103387, 101021, 98639, 96179, 93911, 91583, 89317, 86939, 84521,
1201 82183, 79939, 77587, 75307, 72959, 70793, 68447, 66103
1202 };
1203}
1204int
1205ClpHashValue::hash ( double value) const
1206{
1207
1208 union {
1209 double d;
1210 char c[8];
1211 } v1;
1212 assert (sizeof(double) == 8);
1213 v1.d = value;
1214 int n = 0;
1215 int j;
1216
1217 for ( j = 0; j < 8; ++j ) {
1218 int ichar = v1.c[j];
1219 n += mmult_for_hash[j] * ichar;
1220 }
1221 return ( abs ( n ) % maxHash_ ); /* integer abs */
1222}
1223void
1224ClpHashValue::resize(bool increaseMax)
1225{
1226 int newSize = increaseMax ? ((3 * maxHash_) >> 1) + 1000 : maxHash_;
1227 CoinHashLink * newHash = new CoinHashLink[newSize];
1228 int i;
1229 for ( i = 0; i < newSize; i++ ) {
1230 newHash[i].value = -1.0e-100;
1231 newHash[i].index = -1;
1232 newHash[i].next = -1;
1233 }
1234 // swap
1235 CoinHashLink * oldHash = hash_;
1236 hash_ = newHash;
1237 int oldSize = maxHash_;
1238 maxHash_ = newSize;
1239 /*
1240 * Initialize the hash table. Only the index of the first value that
1241 * hashes to a value is entered in the table; subsequent values that
1242 * collide with it are not entered.
1243 */
1244 int ipos;
1245 int n = 0;
1246 for ( i = 0; i < oldSize; i++ ) {
1247 if (oldHash[i].index >= 0) {
1248 ipos = hash ( oldHash[i].value);
1249 if ( hash_[ipos].index == -1 ) {
1250 hash_[ipos].index = n;
1251 n++;
1252 hash_[ipos].value = oldHash[i].value;
1253 // unmark
1254 oldHash[i].index = -1;
1255 }
1256 }
1257 }
1258 /*
1259 * Now take care of the values that collided in the preceding loop,
1260 * by finding some other entry in the table for them.
1261 * Since there are as many entries in the table as there are values,
1262 * there must be room for them.
1263 */
1264 lastUsed_ = -1;
1265 for ( i = 0; i < oldSize; ++i ) {
1266 if (oldHash[i].index >= 0) {
1267 double value = oldHash[i].value;
1268 ipos = hash ( value);
1269 int k;
1270 while ( true ) {
1271 assert (value != hash_[ipos].value);
1272 k = hash_[ipos].next;
1273 if ( k == -1 ) {
1274 while ( true ) {
1275 ++lastUsed_;
1276 assert (lastUsed_ <= maxHash_);
1277 if ( hash_[lastUsed_].index == -1 ) {
1278 break;
1279 }
1280 }
1281 hash_[ipos].next = lastUsed_;
1282 hash_[lastUsed_].index = n;
1283 n++;
1284 hash_[lastUsed_].value = value;
1285 break;
1286 } else {
1287 ipos = k;
1288 }
1289 }
1290 }
1291 }
1292 assert (n == numberHash_);
1293 delete [] oldHash;
1294}
1295