1// Copyright (C) 2006, International Business Machines
2// Corporation and others. All Rights Reserved.
3// This code is licensed under the terms of the Eclipse Public License (EPL).
4
5#if defined(_MSC_VER)
6// Turn off compiler warning about long names
7# pragma warning(disable:4786)
8#endif
9#include <cassert>
10#include <cstdlib>
11#include <cmath>
12#include <cfloat>
13//#define OSI_DEBUG
14#include "OsiSolverInterface.hpp"
15#include "OsiBranchingObject.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "CoinPackedMatrix.hpp"
18#include "CoinSort.hpp"
19#include "CoinError.hpp"
20#include "CoinFinite.hpp"
21
22// Default Constructor
23OsiObject::OsiObject()
24 :infeasibility_(0.0),
25 whichWay_(0),
26 numberWays_(2),
27 priority_(1000)
28{
29}
30
31
32// Destructor
33OsiObject::~OsiObject ()
34{
35}
36
37// Copy constructor
38OsiObject::OsiObject ( const OsiObject & rhs)
39{
40 infeasibility_ = rhs.infeasibility_;
41 whichWay_ = rhs.whichWay_;
42 priority_ = rhs.priority_;
43 numberWays_ = rhs.numberWays_;
44}
45
46// Assignment operator
47OsiObject &
48OsiObject::operator=( const OsiObject& rhs)
49{
50 if (this!=&rhs) {
51 infeasibility_ = rhs.infeasibility_;
52 whichWay_ = rhs.whichWay_;
53 priority_ = rhs.priority_;
54 numberWays_ = rhs.numberWays_;
55 }
56 return *this;
57}
58// Return "up" estimate (default 1.0e-5)
59double
60OsiObject::upEstimate() const
61{
62 return 1.0e-5;
63}
64// Return "down" estimate (default 1.0e-5)
65double
66OsiObject::downEstimate() const
67{
68 return 1.0e-5;
69}
70// Column number if single column object -1 otherwise
71int
72OsiObject::columnNumber() const
73{
74 return -1;
75}
76// Infeasibility - large is 0.5
77double
78OsiObject::infeasibility(const OsiSolverInterface * solver, int & preferredWay) const
79{
80 // Can't guarantee has matrix
81 OsiBranchingInformation info(solver,false,false);
82 return infeasibility(&info,preferredWay);
83}
84// This does NOT set mutable stuff
85double
86OsiObject::checkInfeasibility(const OsiBranchingInformation * info) const
87{
88 int way;
89 double saveInfeasibility = infeasibility_;
90 short int saveWhichWay = whichWay_ ;
91 double value = infeasibility(info,way);
92 infeasibility_ = saveInfeasibility;
93 whichWay_ = saveWhichWay;
94 return value;
95}
96
97/* For the variable(s) referenced by the object,
98 look at the current solution and set bounds to match the solution.
99 Returns measure of how much it had to move solution to make feasible
100*/
101double
102OsiObject::feasibleRegion(OsiSolverInterface * solver) const
103{
104 // Can't guarantee has matrix
105 OsiBranchingInformation info(solver,false,false);
106 return feasibleRegion(solver,&info);
107}
108
109// Default Constructor
110OsiObject2::OsiObject2()
111 : OsiObject(),
112 preferredWay_(-1),
113 otherInfeasibility_(0.0)
114{
115}
116
117
118// Destructor
119OsiObject2::~OsiObject2 ()
120{
121}
122
123// Copy constructor
124OsiObject2::OsiObject2 ( const OsiObject2 & rhs)
125 : OsiObject(rhs),
126 preferredWay_(rhs.preferredWay_),
127 otherInfeasibility_ (rhs.otherInfeasibility_)
128{
129}
130
131// Assignment operator
132OsiObject2 &
133OsiObject2::operator=( const OsiObject2& rhs)
134{
135 if (this!=&rhs) {
136 OsiObject::operator=(rhs);
137 preferredWay_ = rhs.preferredWay_;
138 otherInfeasibility_ = rhs.otherInfeasibility_;
139 }
140 return *this;
141}
142// Default Constructor
143OsiBranchingObject::OsiBranchingObject()
144{
145 originalObject_=NULL;
146 branchIndex_=0;
147 value_=0.0;
148 numberBranches_=2;
149}
150
151// Useful constructor
152OsiBranchingObject::OsiBranchingObject (OsiSolverInterface * ,
153 double value)
154{
155 originalObject_=NULL;
156 branchIndex_=0;
157 value_=value;
158 numberBranches_=2;
159}
160
161// Copy constructor
162OsiBranchingObject::OsiBranchingObject ( const OsiBranchingObject & rhs)
163{
164 originalObject_=rhs.originalObject_;
165 branchIndex_=rhs.branchIndex_;
166 value_=rhs.value_;
167 numberBranches_=rhs.numberBranches_;
168}
169
170// Assignment operator
171OsiBranchingObject &
172OsiBranchingObject::operator=( const OsiBranchingObject& rhs)
173{
174 if (this != &rhs) {
175 originalObject_=rhs.originalObject_;
176 branchIndex_=rhs.branchIndex_;
177 value_=rhs.value_;
178 numberBranches_=rhs.numberBranches_;
179 }
180 return *this;
181}
182
183// Destructor
184OsiBranchingObject::~OsiBranchingObject ()
185{
186}
187// For debug
188int
189OsiBranchingObject::columnNumber() const
190{
191 if (originalObject_)
192 return originalObject_->columnNumber();
193 else
194 return -1;
195}
196/** Default Constructor
197
198*/
199OsiBranchingInformation::OsiBranchingInformation ()
200 : objectiveValue_(COIN_DBL_MAX),
201 cutoff_(COIN_DBL_MAX),
202 direction_(COIN_DBL_MAX),
203 integerTolerance_(1.0e-7),
204 primalTolerance_(1.0e-7),
205 timeRemaining_(COIN_DBL_MAX),
206 defaultDual_(-1.0),
207 solver_(NULL),
208 numberColumns_(0),
209 lower_(NULL),
210 solution_(NULL),
211 upper_(NULL),
212 hotstartSolution_(NULL),
213 pi_(NULL),
214 rowActivity_(NULL),
215 objective_(NULL),
216 rowLower_(NULL),
217 rowUpper_(NULL),
218 elementByColumn_(NULL),
219 columnStart_(NULL),
220 columnLength_(NULL),
221 row_(NULL),
222 usefulRegion_(NULL),
223 indexRegion_(NULL),
224 numberSolutions_(0),
225 numberBranchingSolutions_(0),
226 depth_(0),
227 owningSolution_(false)
228{
229}
230
231/** Useful constructor
232*/
233OsiBranchingInformation::OsiBranchingInformation (const OsiSolverInterface * solver,
234 bool /*normalSolver*/,
235 bool owningSolution)
236 : timeRemaining_(COIN_DBL_MAX),
237 defaultDual_(-1.0),
238 solver_(solver),
239 hotstartSolution_(NULL),
240 usefulRegion_(NULL),
241 indexRegion_(NULL),
242 numberSolutions_(0),
243 numberBranchingSolutions_(0),
244 depth_(0),
245 owningSolution_(owningSolution)
246{
247 direction_ = solver_->getObjSense();
248 objectiveValue_ = solver_->getObjValue();
249 objectiveValue_ *= direction_;
250 solver_->getDblParam(OsiDualObjectiveLimit,cutoff_) ;
251 cutoff_ *= direction_;
252 integerTolerance_ = solver_->getIntegerTolerance();
253 solver_->getDblParam(OsiPrimalTolerance,primalTolerance_) ;
254 numberColumns_ = solver_->getNumCols();
255 lower_ = solver_->getColLower();
256 if (owningSolution_)
257 solution_ = CoinCopyOfArray(solver_->getColSolution(),numberColumns_);
258 else
259 solution_ = solver_->getColSolution();
260 upper_ = solver_->getColUpper();
261 pi_ = solver_->getRowPrice();
262 rowActivity_ = solver_->getRowActivity();
263 objective_ = solver_->getObjCoefficients();
264 rowLower_ = solver_->getRowLower();
265 rowUpper_ = solver_->getRowUpper();
266 const CoinPackedMatrix* matrix = solver_->getMatrixByCol();
267 if (matrix) {
268 // Column copy of matrix if matrix exists
269 elementByColumn_ = matrix->getElements();
270 row_ = matrix->getIndices();
271 columnStart_ = matrix->getVectorStarts();
272 columnLength_ = matrix->getVectorLengths();
273 } else {
274 // Matrix does not exist
275 elementByColumn_ = NULL;
276 row_ = NULL;
277 columnStart_ = NULL;
278 columnLength_ = NULL;
279 }
280}
281// Copy constructor
282OsiBranchingInformation::OsiBranchingInformation ( const OsiBranchingInformation & rhs)
283{
284 objectiveValue_ = rhs.objectiveValue_;
285 cutoff_ = rhs.cutoff_;
286 direction_ = rhs.direction_;
287 integerTolerance_ = rhs.integerTolerance_;
288 primalTolerance_ = rhs.primalTolerance_;
289 timeRemaining_ = rhs.timeRemaining_;
290 defaultDual_ = rhs.defaultDual_;
291 solver_ = rhs.solver_;
292 numberColumns_ = rhs.numberColumns_;
293 lower_ = rhs.lower_;
294 owningSolution_ = rhs.owningSolution_;
295 if (owningSolution_)
296 solution_ = CoinCopyOfArray(rhs.solution_,numberColumns_);
297 else
298 solution_ = rhs.solution_;
299 upper_ = rhs.upper_;
300 hotstartSolution_ = rhs.hotstartSolution_;
301 pi_ = rhs.pi_;
302 rowActivity_ = rhs.rowActivity_;
303 objective_ = rhs.objective_;
304 rowLower_ = rhs.rowLower_;
305 rowUpper_ = rhs.rowUpper_;
306 elementByColumn_ = rhs.elementByColumn_;
307 row_ = rhs.row_;
308 columnStart_ = rhs.columnStart_;
309 columnLength_ = rhs.columnLength_;
310 usefulRegion_ = rhs.usefulRegion_;
311 assert (!usefulRegion_);
312 indexRegion_ = rhs.indexRegion_;
313 numberSolutions_ = rhs.numberSolutions_;
314 numberBranchingSolutions_ = rhs.numberBranchingSolutions_;
315 depth_ = rhs.depth_;
316}
317
318// Clone
319OsiBranchingInformation *
320OsiBranchingInformation::clone() const
321{
322 return new OsiBranchingInformation(*this);
323}
324
325// Assignment operator
326OsiBranchingInformation &
327OsiBranchingInformation::operator=( const OsiBranchingInformation& rhs)
328{
329 if (this!=&rhs) {
330 objectiveValue_ = rhs.objectiveValue_;
331 cutoff_ = rhs.cutoff_;
332 direction_ = rhs.direction_;
333 integerTolerance_ = rhs.integerTolerance_;
334 primalTolerance_ = rhs.primalTolerance_;
335 timeRemaining_ = rhs.timeRemaining_;
336 defaultDual_ = rhs.defaultDual_;
337 numberColumns_ = rhs.numberColumns_;
338 lower_ = rhs.lower_;
339 owningSolution_ = rhs.owningSolution_;
340 if (owningSolution_) {
341 solution_ = CoinCopyOfArray(rhs.solution_,numberColumns_);
342 delete [] solution_;
343 } else {
344 solution_ = rhs.solution_;
345 }
346 upper_ = rhs.upper_;
347 hotstartSolution_ = rhs.hotstartSolution_;
348 pi_ = rhs.pi_;
349 rowActivity_ = rhs.rowActivity_;
350 objective_ = rhs.objective_;
351 rowLower_ = rhs.rowLower_;
352 rowUpper_ = rhs.rowUpper_;
353 elementByColumn_ = rhs.elementByColumn_;
354 row_ = rhs.row_;
355 columnStart_ = rhs.columnStart_;
356 columnLength_ = rhs.columnLength_;
357 usefulRegion_ = rhs.usefulRegion_;
358 assert (!usefulRegion_);
359 indexRegion_ = rhs.indexRegion_;
360 numberSolutions_ = rhs.numberSolutions_;
361 numberBranchingSolutions_ = rhs.numberBranchingSolutions_;
362 depth_ = rhs.depth_;
363 }
364 return *this;
365}
366
367// Destructor
368OsiBranchingInformation::~OsiBranchingInformation ()
369{
370 if (owningSolution_)
371 delete[] solution_;
372}
373// Default Constructor
374OsiTwoWayBranchingObject::OsiTwoWayBranchingObject()
375 :OsiBranchingObject()
376{
377 firstBranch_=0;
378}
379
380// Useful constructor
381OsiTwoWayBranchingObject::OsiTwoWayBranchingObject (OsiSolverInterface * solver,
382 const OsiObject * object,
383 int way , double value)
384 :OsiBranchingObject(solver,value)
385{
386 originalObject_ = object;
387 firstBranch_=way;
388}
389
390
391// Copy constructor
392OsiTwoWayBranchingObject::OsiTwoWayBranchingObject ( const OsiTwoWayBranchingObject & rhs) :OsiBranchingObject(rhs)
393{
394 firstBranch_=rhs.firstBranch_;
395}
396
397// Assignment operator
398OsiTwoWayBranchingObject &
399OsiTwoWayBranchingObject::operator=( const OsiTwoWayBranchingObject& rhs)
400{
401 if (this != &rhs) {
402 OsiBranchingObject::operator=(rhs);
403 firstBranch_=rhs.firstBranch_;
404 }
405 return *this;
406}
407
408// Destructor
409OsiTwoWayBranchingObject::~OsiTwoWayBranchingObject ()
410{
411}
412
413/********* Simple Integers *******************************/
414/** Default Constructor
415
416 Equivalent to an unspecified binary variable.
417*/
418OsiSimpleInteger::OsiSimpleInteger ()
419 : OsiObject2(),
420 originalLower_(0.0),
421 originalUpper_(1.0),
422 columnNumber_(-1)
423{
424}
425
426/** Useful constructor
427
428 Loads actual upper & lower bounds for the specified variable.
429*/
430OsiSimpleInteger::OsiSimpleInteger (const OsiSolverInterface * solver, int iColumn)
431 : OsiObject2()
432{
433 columnNumber_ = iColumn ;
434 originalLower_ = solver->getColLower()[columnNumber_] ;
435 originalUpper_ = solver->getColUpper()[columnNumber_] ;
436}
437
438
439// Useful constructor - passed solver index and original bounds
440OsiSimpleInteger::OsiSimpleInteger ( int iColumn, double lower, double upper)
441 : OsiObject2()
442{
443 columnNumber_ = iColumn ;
444 originalLower_ = lower;
445 originalUpper_ = upper;
446}
447
448// Copy constructor
449OsiSimpleInteger::OsiSimpleInteger ( const OsiSimpleInteger & rhs)
450 :OsiObject2(rhs)
451
452{
453 columnNumber_ = rhs.columnNumber_;
454 originalLower_ = rhs.originalLower_;
455 originalUpper_ = rhs.originalUpper_;
456}
457
458// Clone
459OsiObject *
460OsiSimpleInteger::clone() const
461{
462 return new OsiSimpleInteger(*this);
463}
464
465// Assignment operator
466OsiSimpleInteger &
467OsiSimpleInteger::operator=( const OsiSimpleInteger& rhs)
468{
469 if (this!=&rhs) {
470 OsiObject2::operator=(rhs);
471 columnNumber_ = rhs.columnNumber_;
472 originalLower_ = rhs.originalLower_;
473 originalUpper_ = rhs.originalUpper_;
474 }
475 return *this;
476}
477
478// Destructor
479OsiSimpleInteger::~OsiSimpleInteger ()
480{
481}
482/* Reset variable bounds to their original values.
483
484Bounds may be tightened, so it may be good to be able to reset them to
485their original values.
486*/
487void
488OsiSimpleInteger::resetBounds(const OsiSolverInterface * solver)
489{
490 originalLower_ = solver->getColLower()[columnNumber_] ;
491 originalUpper_ = solver->getColUpper()[columnNumber_] ;
492}
493// Redoes data when sequence numbers change
494void
495OsiSimpleInteger::resetSequenceEtc(int numberColumns, const int * originalColumns)
496{
497 int i;
498 for (i=0;i<numberColumns;i++) {
499 if (originalColumns[i]==columnNumber_)
500 break;
501 }
502 if (i<numberColumns)
503 columnNumber_=i;
504 else
505 abort(); // should never happen
506}
507
508// Infeasibility - large is 0.5
509double
510OsiSimpleInteger::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
511{
512 double value = info->solution_[columnNumber_];
513 value = CoinMax(value, info->lower_[columnNumber_]);
514 value = CoinMin(value, info->upper_[columnNumber_]);
515 double nearest = floor(value+(1.0-0.5));
516 if (nearest>value) {
517 whichWay=1;
518 } else {
519 whichWay=0;
520 }
521 infeasibility_ = fabs(value-nearest);
522 double returnValue = infeasibility_;
523 if (infeasibility_<=info->integerTolerance_) {
524 otherInfeasibility_ = 1.0;
525 returnValue = 0.0;
526 } else if (info->defaultDual_<0.0) {
527 otherInfeasibility_ = 1.0-infeasibility_;
528 } else {
529 const double * pi = info->pi_;
530 const double * activity = info->rowActivity_;
531 const double * lower = info->rowLower_;
532 const double * upper = info->rowUpper_;
533 const double * element = info->elementByColumn_;
534 const int * row = info->row_;
535 const CoinBigIndex * columnStart = info->columnStart_;
536 const int * columnLength = info->columnLength_;
537 double direction = info->direction_;
538 double downMovement = value - floor(value);
539 double upMovement = 1.0-downMovement;
540 double valueP = info->objective_[columnNumber_]*direction;
541 CoinBigIndex start = columnStart[columnNumber_];
542 CoinBigIndex end = start + columnLength[columnNumber_];
543 double upEstimate = 0.0;
544 double downEstimate = 0.0;
545 if (valueP>0.0)
546 upEstimate = valueP*upMovement;
547 else
548 downEstimate -= valueP*downMovement;
549 double tolerance = info->primalTolerance_;
550 for (CoinBigIndex j=start;j<end;j++) {
551 int iRow = row[j];
552 if (lower[iRow]<-1.0e20)
553 assert (pi[iRow]<=1.0e-4);
554 if (upper[iRow]>1.0e20)
555 assert (pi[iRow]>=-1.0e-4);
556 valueP = pi[iRow]*direction;
557 double el2 = element[j];
558 double value2 = valueP*el2;
559 double u=0.0;
560 double d=0.0;
561 if (value2>0.0)
562 u = value2;
563 else
564 d = -value2;
565 // if up makes infeasible then make at least default
566 double newUp = activity[iRow] + upMovement*el2;
567 if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
568 u = CoinMax(u,info->defaultDual_);
569 upEstimate += u*upMovement;
570 // if down makes infeasible then make at least default
571 double newDown = activity[iRow] - downMovement*el2;
572 if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
573 d = CoinMax(d,info->defaultDual_);
574 downEstimate += d*downMovement;
575 }
576 if (downEstimate>=upEstimate) {
577 infeasibility_ = CoinMax(1.0e-12,upEstimate);
578 otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
579 whichWay = 1;
580 } else {
581 infeasibility_ = CoinMax(1.0e-12,downEstimate);
582 otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
583 whichWay = 0;
584 }
585 returnValue = infeasibility_;
586 }
587 if (preferredWay_>=0&&returnValue)
588 whichWay = preferredWay_;
589 whichWay_ = static_cast<short int>(whichWay) ;
590 return returnValue;
591}
592
593// This looks at solution and sets bounds to contain solution
594/** More precisely: it first forces the variable within the existing
595 bounds, and then tightens the bounds to fix the variable at the
596 nearest integer value.
597*/
598double
599OsiSimpleInteger::feasibleRegion(OsiSolverInterface * solver,
600 const OsiBranchingInformation * info) const
601{
602 double value = info->solution_[columnNumber_];
603 double newValue = CoinMax(value, info->lower_[columnNumber_]);
604 newValue = CoinMin(newValue, info->upper_[columnNumber_]);
605 newValue = floor(newValue+0.5);
606 solver->setColLower(columnNumber_,newValue);
607 solver->setColUpper(columnNumber_,newValue);
608 return fabs(value-newValue);
609}
610/* Column number if single column object -1 otherwise,
611 so returns >= 0
612 Used by heuristics
613*/
614int
615OsiSimpleInteger::columnNumber() const
616{
617 return columnNumber_;
618}
619// Creates a branching object
620OsiBranchingObject *
621OsiSimpleInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
622{
623 double value = info->solution_[columnNumber_];
624 value = CoinMax(value, info->lower_[columnNumber_]);
625 value = CoinMin(value, info->upper_[columnNumber_]);
626 assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
627#ifndef NDEBUG
628 double nearest = floor(value+0.5);
629 assert (fabs(value-nearest)>info->integerTolerance_);
630#endif
631 OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
632 value);
633 return branch;
634}
635// Return "down" estimate
636double
637OsiSimpleInteger::downEstimate() const
638{
639 if (whichWay_)
640 return 1.0-infeasibility_;
641 else
642 return infeasibility_;
643}
644// Return "up" estimate
645double
646OsiSimpleInteger::upEstimate() const
647{
648 if (!whichWay_)
649 return 1.0-infeasibility_;
650 else
651 return infeasibility_;
652}
653
654// Default Constructor
655OsiIntegerBranchingObject::OsiIntegerBranchingObject()
656 :OsiTwoWayBranchingObject()
657{
658 down_[0] = 0.0;
659 down_[1] = 0.0;
660 up_[0] = 0.0;
661 up_[1] = 0.0;
662}
663
664// Useful constructor
665OsiIntegerBranchingObject::OsiIntegerBranchingObject (OsiSolverInterface * solver,
666 const OsiSimpleInteger * object,
667 int way , double value)
668 :OsiTwoWayBranchingObject(solver,object, way, value)
669{
670 int iColumn = object->columnNumber();
671 down_[0] = solver->getColLower()[iColumn];
672 down_[1] = floor(value_);
673 up_[0] = ceil(value_);
674 up_[1] = solver->getColUpper()[iColumn];
675}
676/* Create a standard floor/ceiling branch object
677 Specifies a simple two-way branch in a more flexible way. One arm of the
678 branch will be lb <= x <= downUpperBound, the other upLowerBound <= x <= ub.
679 Specify way = -1 to set the object state to perform the down arm first,
680 way = 1 for the up arm.
681*/
682OsiIntegerBranchingObject::OsiIntegerBranchingObject (OsiSolverInterface * solver,
683 const OsiSimpleInteger * object,
684 int way , double value, double downUpperBound,
685 double upLowerBound)
686 :OsiTwoWayBranchingObject(solver,object, way, value)
687{
688 int iColumn = object->columnNumber();
689 down_[0] = solver->getColLower()[iColumn];
690 down_[1] = downUpperBound;
691 up_[0] = upLowerBound;
692 up_[1] = solver->getColUpper()[iColumn];
693}
694
695
696// Copy constructor
697OsiIntegerBranchingObject::OsiIntegerBranchingObject ( const OsiIntegerBranchingObject & rhs) :OsiTwoWayBranchingObject(rhs)
698{
699 down_[0] = rhs.down_[0];
700 down_[1] = rhs.down_[1];
701 up_[0] = rhs.up_[0];
702 up_[1] = rhs.up_[1];
703}
704
705// Assignment operator
706OsiIntegerBranchingObject &
707OsiIntegerBranchingObject::operator=( const OsiIntegerBranchingObject& rhs)
708{
709 if (this != &rhs) {
710 OsiTwoWayBranchingObject::operator=(rhs);
711 down_[0] = rhs.down_[0];
712 down_[1] = rhs.down_[1];
713 up_[0] = rhs.up_[0];
714 up_[1] = rhs.up_[1];
715 }
716 return *this;
717}
718OsiBranchingObject *
719OsiIntegerBranchingObject::clone() const
720{
721 return (new OsiIntegerBranchingObject(*this));
722}
723
724
725// Destructor
726OsiIntegerBranchingObject::~OsiIntegerBranchingObject ()
727{
728}
729
730/*
731 Perform a branch by adjusting the bounds of the specified variable. Note
732 that each arm of the branch advances the object to the next arm by
733 advancing the value of branchIndex_.
734
735 Providing new values for the variable's lower and upper bounds for each
736 branching direction gives a little bit of additional flexibility and will
737 be easily extensible to multi-way branching.
738 Returns change in guessed objective on next branch
739*/
740double
741OsiIntegerBranchingObject::branch(OsiSolverInterface * solver)
742{
743 const OsiSimpleInteger * obj =
744 dynamic_cast <const OsiSimpleInteger *>(originalObject_) ;
745 assert (obj);
746 int iColumn = obj->columnNumber();
747 double olb,oub ;
748 olb = solver->getColLower()[iColumn] ;
749 oub = solver->getColUpper()[iColumn] ;
750 int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
751 if (0) {
752 printf("branching %s on %d bounds %g %g / %g %g\n",
753 (way==-1) ? "down" :"up",iColumn,
754 down_[0],down_[1],up_[0],up_[1]);
755 const double * lower = solver->getColLower();
756 const double * upper = solver->getColUpper();
757 for (int i=0;i<8;i++)
758 printf(" [%d (%g,%g)]",i,lower[i],upper[i]);
759 printf("\n");
760 }
761 if (way<0) {
762#ifdef OSI_DEBUG
763 { double olb,oub ;
764 olb = solver->getColLower()[iColumn] ;
765 oub = solver->getColUpper()[iColumn] ;
766 printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
767 iColumn,olb,oub,down_[0],down_[1]) ; }
768#endif
769 solver->setColLower(iColumn,down_[0]);
770 solver->setColUpper(iColumn,down_[1]);
771 } else {
772#ifdef OSI_DEBUG
773 { double olb,oub ;
774 olb = solver->getColLower()[iColumn] ;
775 oub = solver->getColUpper()[iColumn] ;
776 printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
777 iColumn,olb,oub,up_[0],up_[1]) ; }
778#endif
779 solver->setColLower(iColumn,up_[0]);
780 solver->setColUpper(iColumn,up_[1]);
781 }
782 double nlb = solver->getColLower()[iColumn];
783 if (nlb<olb) {
784#ifndef NDEBUG
785 printf("bad lb change for column %d from %g to %g\n",iColumn,olb,nlb);
786#endif
787 solver->setColLower(iColumn,olb);
788 }
789 double nub = solver->getColUpper()[iColumn];
790 if (nub>oub) {
791#ifndef NDEBUG
792 printf("bad ub change for column %d from %g to %g\n",iColumn,oub,nub);
793#endif
794 solver->setColUpper(iColumn,oub);
795 }
796#ifndef NDEBUG
797 if (nlb<olb+1.0e-8&&nub>oub-1.0e-8)
798 printf("bad null change for column %d - bounds %g,%g\n",iColumn,olb,oub);
799#endif
800 branchIndex_++;
801 return 0.0;
802}
803// Print what would happen
804void
805OsiIntegerBranchingObject::print(const OsiSolverInterface * solver)
806{
807 const OsiSimpleInteger * obj =
808 dynamic_cast <const OsiSimpleInteger *>(originalObject_) ;
809 assert (obj);
810 int iColumn = obj->columnNumber();
811 int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
812 if (way<0) {
813 { double olb,oub ;
814 olb = solver->getColLower()[iColumn] ;
815 oub = solver->getColUpper()[iColumn] ;
816 printf("OsiInteger would branch down on var %d : [%g,%g] => [%g,%g]\n",
817 iColumn,olb,oub,down_[0],down_[1]) ; }
818 } else {
819 { double olb,oub ;
820 olb = solver->getColLower()[iColumn] ;
821 oub = solver->getColUpper()[iColumn] ;
822 printf("OsiInteger would branch up on var %d : [%g,%g] => [%g,%g]\n",
823 iColumn,olb,oub,up_[0],up_[1]) ; }
824 }
825}
826// Default Constructor
827OsiSOS::OsiSOS ()
828 : OsiObject2(),
829 members_(NULL),
830 weights_(NULL),
831 numberMembers_(0),
832 sosType_(-1),
833 integerValued_(false)
834{
835}
836
837// Useful constructor (which are indices)
838OsiSOS::OsiSOS (const OsiSolverInterface * , int numberMembers,
839 const int * which, const double * weights, int type)
840 : OsiObject2(),
841 numberMembers_(numberMembers),
842 sosType_(type)
843{
844 integerValued_ = type==1; // not strictly true - should check problem
845 if (numberMembers_) {
846 members_ = new int[numberMembers_];
847 weights_ = new double[numberMembers_];
848 memcpy(members_,which,numberMembers_*sizeof(int));
849 if (weights) {
850 memcpy(weights_,weights,numberMembers_*sizeof(double));
851 } else {
852 for (int i=0;i<numberMembers_;i++)
853 weights_[i]=i;
854 }
855 // sort so weights increasing
856 CoinSort_2(weights_,weights_+numberMembers_,members_);
857 double last = -COIN_DBL_MAX;
858 int i;
859 for (i=0;i<numberMembers_;i++) {
860 double possible = CoinMax(last+1.0e-10,weights_[i]);
861 weights_[i] = possible;
862 last=possible;
863 }
864 } else {
865 members_ = NULL;
866 weights_ = NULL;
867 }
868 assert (sosType_>0&&sosType_<3);
869}
870
871// Copy constructor
872OsiSOS::OsiSOS ( const OsiSOS & rhs)
873 :OsiObject2(rhs)
874{
875 numberMembers_ = rhs.numberMembers_;
876 sosType_ = rhs.sosType_;
877 integerValued_ = rhs.integerValued_;
878 if (numberMembers_) {
879 members_ = new int[numberMembers_];
880 weights_ = new double[numberMembers_];
881 memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
882 memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
883 } else {
884 members_ = NULL;
885 weights_ = NULL;
886 }
887}
888
889// Clone
890OsiObject *
891OsiSOS::clone() const
892{
893 return new OsiSOS(*this);
894}
895
896// Assignment operator
897OsiSOS &
898OsiSOS::operator=( const OsiSOS& rhs)
899{
900 if (this!=&rhs) {
901 OsiObject2::operator=(rhs);
902 delete [] members_;
903 delete [] weights_;
904 numberMembers_ = rhs.numberMembers_;
905 sosType_ = rhs.sosType_;
906 integerValued_ = rhs.integerValued_;
907 if (numberMembers_) {
908 members_ = new int[numberMembers_];
909 weights_ = new double[numberMembers_];
910 memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
911 memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
912 } else {
913 members_ = NULL;
914 weights_ = NULL;
915 }
916 }
917 return *this;
918}
919
920// Destructor
921OsiSOS::~OsiSOS ()
922{
923 delete [] members_;
924 delete [] weights_;
925}
926
927// Infeasibility - large is 0.5
928double
929OsiSOS::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
930{
931 int j;
932 int firstNonZero=-1;
933 int lastNonZero = -1;
934 int firstNonFixed=-1;
935 int lastNonFixed = -1;
936 const double * solution = info->solution_;
937 //const double * lower = info->lower_;
938 const double * upper = info->upper_;
939 //double largestValue=0.0;
940 double integerTolerance = info->integerTolerance_;
941 double primalTolerance = info->primalTolerance_;
942 double weight = 0.0;
943 double sum =0.0;
944
945 // check bounds etc
946 double lastWeight=-1.0e100;
947 for (j=0;j<numberMembers_;j++) {
948 int iColumn = members_[j];
949 if (lastWeight>=weights_[j]-1.0e-12)
950 throw CoinError("Weights too close together in SOS","infeasibility","OsiSOS");
951 lastWeight = weights_[j];
952 if (upper[iColumn]) {
953 double value = CoinMax(0.0,solution[iColumn]);
954 if (value>integerTolerance) {
955 // Possibly due to scaling a fixed variable might slip through
956#ifdef COIN_DEVELOP
957 if (value>upper[iColumn]+10.0*primalTolerance)
958 printf("** Variable %d (%d) has value %g and upper bound of %g\n",
959 iColumn,j,value,upper[iColumn]);
960#endif
961 if (value>upper[iColumn]) {
962 value=upper[iColumn];
963 }
964 sum += value;
965 weight += weights_[j]*value;
966 if (firstNonZero<0)
967 firstNonZero=j;
968 lastNonZero=j;
969 }
970 if (firstNonFixed<0)
971 firstNonFixed=j;
972 lastNonFixed=j;
973 }
974 }
975 whichWay=1;
976 whichWay_=1;
977 if (lastNonZero-firstNonZero>=sosType_) {
978 // find where to branch
979 assert (sum>0.0);
980 // probably best to use pseudo duals
981 double value = lastNonZero-firstNonZero+1;
982 value *= 0.5/static_cast<double> (numberMembers_);
983 infeasibility_=value;
984 otherInfeasibility_=1.0-value;
985 if (info->defaultDual_>=0.0) {
986 // Using pseudo shadow prices
987 weight /= sum;
988 int iWhere;
989 for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
990 if (weight<weights_[iWhere+1])
991 break;
992 assert (iWhere!=lastNonZero);
993 /* Complicated - infeasibility is being used for branching so we
994 don't want estimate of satisfying set but of each way on branch.
995 So let us suppose that all on side being fixed to 0 goes to closest
996 */
997 int lastDown=iWhere;
998 int firstUp=iWhere+1;
999 if (sosType_==2) {
1000 // SOS 2 - choose nearest
1001 if (weight-weights_[iWhere]>=weights_[iWhere+1]-weight)
1002 lastDown++;
1003 // But make sure OK
1004 if (lastDown==firstNonFixed) {
1005 lastDown ++;
1006 } else if (lastDown==lastNonFixed) {
1007 lastDown --;
1008 }
1009 firstUp=lastDown;
1010 }
1011 // Now get current contribution and compute weight for end points
1012 double weightDown = 0.0;
1013 double weightUp = 0.0;
1014 const double * element = info->elementByColumn_;
1015 const int * row = info->row_;
1016 const CoinBigIndex * columnStart = info->columnStart_;
1017 const int * columnLength = info->columnLength_;
1018 double direction = info->direction_;
1019 const double * objective = info->objective_;
1020 // Compute where we would move to
1021 double objValue=0.0;
1022 double * useful = info->usefulRegion_;
1023 int * index = info->indexRegion_;
1024 int n=0;
1025 for (j=firstNonZero;j<=lastNonZero;j++) {
1026 int iColumn = members_[j];
1027 double multiplier = solution[iColumn];
1028 if (j>=lastDown)
1029 weightDown += multiplier;
1030 if (j<=firstUp)
1031 weightUp += multiplier;
1032 if (multiplier>0.0) {
1033 objValue += objective[iColumn]*multiplier;
1034 CoinBigIndex start = columnStart[iColumn];
1035 CoinBigIndex end = start + columnLength[iColumn];
1036 for (CoinBigIndex j=start;j<end;j++) {
1037 int iRow = row[j];
1038 double value = element[j]*multiplier;
1039 if (useful[iRow]) {
1040 value += useful[iRow];
1041 if (!value)
1042 value = 1.0e-100;
1043 } else {
1044 assert (value);
1045 index[n++]=iRow;
1046 }
1047 useful[iRow] = value;
1048 }
1049 }
1050 }
1051 if (sosType_==2)
1052 assert (fabs(weightUp+weightDown-sum-solution[members_[lastDown]])<1.0e-4);
1053 int startX[2];
1054 int endX[2];
1055 startX[0]=firstNonZero;
1056 startX[1]=firstUp;
1057 endX[0]=lastDown;
1058 endX[1]=lastNonZero;
1059 double fakeSolution[2];
1060 int check[2];
1061 fakeSolution[0]=weightDown;
1062 check[0]=members_[lastDown];
1063 fakeSolution[1]=weightUp;
1064 check[1]=members_[firstUp];
1065 const double * pi = info->pi_;
1066 const double * activity = info->rowActivity_;
1067 const double * lower = info->rowLower_;
1068 const double * upper = info->rowUpper_;
1069 int numberRows = info->solver_->getNumRows();
1070 double * useful2 = useful+numberRows;
1071 int * index2 = index+numberRows;
1072 for (int i=0;i<2;i++) {
1073 double obj=0.0;
1074 int n2=0;
1075 for (j=startX[i];j<=endX[i];j++) {
1076 int iColumn = members_[j];
1077 double multiplier = solution[iColumn];
1078 if (iColumn==check[i])
1079 multiplier=fakeSolution[i];
1080 if (multiplier>0.0) {
1081 obj += objective[iColumn]*multiplier;
1082 CoinBigIndex start = columnStart[iColumn];
1083 CoinBigIndex end = start + columnLength[iColumn];
1084 for (CoinBigIndex j=start;j<end;j++) {
1085 int iRow = row[j];
1086 double value = element[j]*multiplier;
1087 if (useful2[iRow]) {
1088 value += useful2[iRow];
1089 if (!value)
1090 value = 1.0e-100;
1091 } else {
1092 assert (value);
1093 index2[n2++]=iRow;
1094 }
1095 useful2[iRow] = value;
1096 }
1097 }
1098 }
1099 // movement in objective
1100 obj = (obj-objValue) * direction;
1101 double estimate = (obj>0.0) ? obj : 0.0;
1102 for (j=0;j<n;j++) {
1103 int iRow = index[j];
1104 // movement
1105 double movement = useful2[iRow]-useful[iRow];
1106 useful[iRow]=0.0;
1107 useful2[iRow]=0.0;
1108 double valueP = pi[iRow]*direction;
1109 if (lower[iRow]<-1.0e20)
1110 assert (valueP<=1.0e-4);
1111 if (upper[iRow]>1.0e20)
1112 assert (valueP>=-1.0e-4);
1113 double value2 = valueP*movement;
1114 double thisEstimate = (value2>0.0) ? value2 : 0;
1115 // if makes infeasible then make at least default
1116 double newValue = activity[iRow] + movement;
1117 if (newValue>upper[iRow]+primalTolerance||newValue<lower[iRow]-primalTolerance)
1118 thisEstimate = CoinMax(thisEstimate,info->defaultDual_);
1119 estimate += thisEstimate;
1120 }
1121 for (j=0;j<n2;j++) {
1122 int iRow = index2[j];
1123 // movement
1124 double movement = useful2[iRow]-useful[iRow];
1125 useful[iRow]=0.0;
1126 useful2[iRow]=0.0;
1127 if (movement) {
1128 double valueP = pi[iRow]*direction;
1129 if (lower[iRow]<-1.0e20)
1130 assert (valueP<=1.0e-4);
1131 if (upper[iRow]>1.0e20)
1132 assert (valueP>=-1.0e-4);
1133 double value2 = valueP*movement;
1134 double thisEstimate = (value2>0.0) ? value2 : 0;
1135 // if makes infeasible then make at least default
1136 double newValue = activity[iRow] + movement;
1137 if (newValue>upper[iRow]+primalTolerance||newValue<lower[iRow]-primalTolerance)
1138 thisEstimate = CoinMax(thisEstimate,info->defaultDual_);
1139 estimate += thisEstimate;
1140 }
1141 }
1142 // store in fakeSolution
1143 fakeSolution[i]=estimate;
1144 }
1145 double downEstimate = fakeSolution[0];
1146 double upEstimate = fakeSolution[1];
1147 if (downEstimate>=upEstimate) {
1148 infeasibility_ = CoinMax(1.0e-12,upEstimate);
1149 otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
1150 whichWay = 1;
1151 } else {
1152 infeasibility_ = CoinMax(1.0e-12,downEstimate);
1153 otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
1154 whichWay = 0;
1155 }
1156 whichWay_=static_cast<short>(whichWay);
1157 value=infeasibility_;
1158 }
1159 return value;
1160 } else {
1161 infeasibility_=0.0;
1162 otherInfeasibility_=1.0;
1163 return 0.0; // satisfied
1164 }
1165}
1166
1167// This looks at solution and sets bounds to contain solution
1168double
1169OsiSOS::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
1170{
1171 int j;
1172 int firstNonZero=-1;
1173 int lastNonZero = -1;
1174 const double * solution = info->solution_;
1175 //const double * lower = info->lower_;
1176 const double * upper = info->upper_;
1177 double sum =0.0;
1178 // Find largest one or pair
1179 double movement=0.0;
1180 if (sosType_==1) {
1181 for (j=0;j<numberMembers_;j++) {
1182 int iColumn = members_[j];
1183 double value = CoinMax(0.0,solution[iColumn]);
1184 if (value>sum&&upper[iColumn]) {
1185 firstNonZero=j;
1186 sum=value;
1187 }
1188 }
1189 lastNonZero=firstNonZero;
1190 } else {
1191 // type 2
1192 for (j=1;j<numberMembers_;j++) {
1193 int iColumn = members_[j];
1194 int jColumn = members_[j-1];
1195 double value1 = CoinMax(0.0,solution[iColumn]);
1196 double value0 = CoinMax(0.0,solution[jColumn]);
1197 double value = value0+value1;
1198 if (value>sum) {
1199 if (upper[iColumn]||upper[jColumn]) {
1200 firstNonZero=upper[jColumn] ? j-1 : j;
1201 lastNonZero=upper[iColumn] ? j : j-1;
1202 sum=value;
1203 }
1204 }
1205 }
1206 }
1207 for (j=0;j<numberMembers_;j++) {
1208 if (j<firstNonZero||j>lastNonZero) {
1209 int iColumn = members_[j];
1210 double value = CoinMax(0.0,solution[iColumn]);
1211 movement += value;
1212 solver->setColUpper(iColumn,0.0);
1213 }
1214 }
1215 return movement;
1216}
1217// Redoes data when sequence numbers change
1218void
1219OsiSOS::resetSequenceEtc(int numberColumns, const int * originalColumns)
1220{
1221 int n2=0;
1222 for (int j=0;j<numberMembers_;j++) {
1223 int iColumn = members_[j];
1224 int i;
1225 for (i=0;i<numberColumns;i++) {
1226 if (originalColumns[i]==iColumn)
1227 break;
1228 }
1229 if (i<numberColumns) {
1230 members_[n2]=i;
1231 weights_[n2++]=weights_[j];
1232 }
1233 }
1234 if (n2<numberMembers_) {
1235 printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
1236 numberMembers_=n2;
1237 }
1238}
1239// Return "down" estimate
1240double
1241OsiSOS::downEstimate() const
1242{
1243 if (whichWay_)
1244 return otherInfeasibility_;
1245 else
1246 return infeasibility_;
1247}
1248// Return "up" estimate
1249double
1250OsiSOS::upEstimate() const
1251{
1252 if (!whichWay_)
1253 return otherInfeasibility_;
1254 else
1255 return infeasibility_;
1256}
1257
1258// Creates a branching object
1259OsiBranchingObject *
1260OsiSOS::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
1261{
1262 int j;
1263 const double * solution = info->solution_;
1264 double tolerance = info->primalTolerance_;
1265 const double * upper = info->upper_;
1266 int firstNonFixed=-1;
1267 int lastNonFixed=-1;
1268 int firstNonZero=-1;
1269 int lastNonZero = -1;
1270 double weight = 0.0;
1271 double sum =0.0;
1272 for (j=0;j<numberMembers_;j++) {
1273 int iColumn = members_[j];
1274 if (upper[iColumn]) {
1275 double value = CoinMax(0.0,solution[iColumn]);
1276 sum += value;
1277 if (firstNonFixed<0)
1278 firstNonFixed=j;
1279 lastNonFixed=j;
1280 if (value>tolerance) {
1281 weight += weights_[j]*value;
1282 if (firstNonZero<0)
1283 firstNonZero=j;
1284 lastNonZero=j;
1285 }
1286 }
1287 }
1288 assert (lastNonZero-firstNonZero>=sosType_) ;
1289 // find where to branch
1290 assert (sum>0.0);
1291 weight /= sum;
1292 int iWhere;
1293 double separator=0.0;
1294 for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
1295 if (weight<weights_[iWhere+1])
1296 break;
1297 if (sosType_==1) {
1298 // SOS 1
1299 separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
1300 } else {
1301 // SOS 2
1302 if (iWhere==lastNonFixed-1)
1303 iWhere = lastNonFixed-2;
1304 separator = weights_[iWhere+1];
1305 }
1306 // create object
1307 OsiBranchingObject * branch;
1308 branch = new OsiSOSBranchingObject(solver,this,way,separator);
1309 return branch;
1310}
1311// Default Constructor
1312OsiSOSBranchingObject::OsiSOSBranchingObject()
1313 :OsiTwoWayBranchingObject()
1314{
1315}
1316
1317// Useful constructor
1318OsiSOSBranchingObject::OsiSOSBranchingObject (OsiSolverInterface * solver,
1319 const OsiSOS * set,
1320 int way ,
1321 double separator)
1322 :OsiTwoWayBranchingObject(solver, set,way,separator)
1323{
1324}
1325
1326// Copy constructor
1327OsiSOSBranchingObject::OsiSOSBranchingObject ( const OsiSOSBranchingObject & rhs) :OsiTwoWayBranchingObject(rhs)
1328{
1329}
1330
1331// Assignment operator
1332OsiSOSBranchingObject &
1333OsiSOSBranchingObject::operator=( const OsiSOSBranchingObject& rhs)
1334{
1335 if (this != &rhs) {
1336 OsiTwoWayBranchingObject::operator=(rhs);
1337 }
1338 return *this;
1339}
1340OsiBranchingObject *
1341OsiSOSBranchingObject::clone() const
1342{
1343 return (new OsiSOSBranchingObject(*this));
1344}
1345
1346
1347// Destructor
1348OsiSOSBranchingObject::~OsiSOSBranchingObject ()
1349{
1350}
1351double
1352OsiSOSBranchingObject::branch(OsiSolverInterface * solver)
1353{
1354 const OsiSOS * set =
1355 dynamic_cast <const OsiSOS *>(originalObject_) ;
1356 assert (set);
1357 int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
1358 branchIndex_++;
1359 int numberMembers = set->numberMembers();
1360 const int * which = set->members();
1361 const double * weights = set->weights();
1362 //const double * lower = solver->getColLower();
1363 //const double * upper = solver->getColUpper();
1364 // *** for way - up means fix all those in down section
1365 if (way<0) {
1366 int i;
1367 for ( i=0;i<numberMembers;i++) {
1368 if (weights[i] > value_)
1369 break;
1370 }
1371 assert (i<numberMembers);
1372 for (;i<numberMembers;i++)
1373 solver->setColUpper(which[i],0.0);
1374 } else {
1375 int i;
1376 for ( i=0;i<numberMembers;i++) {
1377 if (weights[i] >= value_)
1378 break;
1379 else
1380 solver->setColUpper(which[i],0.0);
1381 }
1382 assert (i<numberMembers);
1383 }
1384 return 0.0;
1385}
1386// Print what would happen
1387void
1388OsiSOSBranchingObject::print(const OsiSolverInterface * solver)
1389{
1390 const OsiSOS * set =
1391 dynamic_cast <const OsiSOS *>(originalObject_) ;
1392 assert (set);
1393 int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
1394 int numberMembers = set->numberMembers();
1395 const int * which = set->members();
1396 const double * weights = set->weights();
1397 //const double * lower = solver->getColLower();
1398 const double * upper = solver->getColUpper();
1399 int first=numberMembers;
1400 int last=-1;
1401 int numberFixed=0;
1402 int numberOther=0;
1403 int i;
1404 for ( i=0;i<numberMembers;i++) {
1405 double bound = upper[which[i]];
1406 if (bound) {
1407 first = CoinMin(first,i);
1408 last = CoinMax(last,i);
1409 }
1410 }
1411 // *** for way - up means fix all those in down section
1412 if (way<0) {
1413 printf("SOS Down");
1414 for ( i=0;i<numberMembers;i++) {
1415 double bound = upper[which[i]];
1416 if (weights[i] > value_)
1417 break;
1418 else if (bound)
1419 numberOther++;
1420 }
1421 assert (i<numberMembers);
1422 for (;i<numberMembers;i++) {
1423 double bound = upper[which[i]];
1424 if (bound)
1425 numberFixed++;
1426 }
1427 } else {
1428 printf("SOS Up");
1429 for ( i=0;i<numberMembers;i++) {
1430 double bound = upper[which[i]];
1431 if (weights[i] >= value_)
1432 break;
1433 else if (bound)
1434 numberFixed++;
1435 }
1436 assert (i<numberMembers);
1437 for (;i<numberMembers;i++) {
1438 double bound = upper[which[i]];
1439 if (bound)
1440 numberOther++;
1441 }
1442 }
1443 printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
1444 value_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
1445}
1446/** Default Constructor
1447
1448*/
1449OsiLotsize::OsiLotsize ()
1450 : OsiObject2(),
1451 columnNumber_(-1),
1452 rangeType_(0),
1453 numberRanges_(0),
1454 largestGap_(0),
1455 bound_(NULL),
1456 range_(0)
1457{
1458}
1459
1460/** Useful constructor
1461
1462 Loads actual upper & lower bounds for the specified variable.
1463*/
1464OsiLotsize::OsiLotsize (const OsiSolverInterface * ,
1465 int iColumn, int numberPoints,
1466 const double * points, bool range)
1467 : OsiObject2()
1468{
1469 assert (numberPoints>0);
1470 columnNumber_ = iColumn ;
1471 // sort ranges
1472 int * sort = new int[numberPoints];
1473 double * weight = new double [numberPoints];
1474 int i;
1475 if (range) {
1476 rangeType_=2;
1477 } else {
1478 rangeType_=1;
1479 }
1480 for (i=0;i<numberPoints;i++) {
1481 sort[i]=i;
1482 weight[i]=points[i*rangeType_];
1483 }
1484 CoinSort_2(weight,weight+numberPoints,sort);
1485 numberRanges_=1;
1486 largestGap_=0;
1487 if (rangeType_==1) {
1488 bound_ = new double[numberPoints+1];
1489 bound_[0]=weight[0];
1490 for (i=1;i<numberPoints;i++) {
1491 if (weight[i]!=weight[i-1])
1492 bound_[numberRanges_++]=weight[i];
1493 }
1494 // and for safety
1495 bound_[numberRanges_]=bound_[numberRanges_-1];
1496 for (i=1;i<numberRanges_;i++) {
1497 largestGap_ = CoinMax(largestGap_,bound_[i]-bound_[i-1]);
1498 }
1499 } else {
1500 bound_ = new double[2*numberPoints+2];
1501 bound_[0]=points[sort[0]*2];
1502 bound_[1]=points[sort[0]*2+1];
1503 double lo=bound_[0];
1504 double hi=bound_[1];
1505 assert (hi>=lo);
1506 for (i=1;i<numberPoints;i++) {
1507 double thisLo =points[sort[i]*2];
1508 double thisHi =points[sort[i]*2+1];
1509 assert (thisHi>=thisLo);
1510 if (thisLo>hi) {
1511 bound_[2*numberRanges_]=thisLo;
1512 bound_[2*numberRanges_+1]=thisHi;
1513 numberRanges_++;
1514 lo=thisLo;
1515 hi=thisHi;
1516 } else {
1517 //overlap
1518 hi=CoinMax(hi,thisHi);
1519 bound_[2*numberRanges_-1]=hi;
1520 }
1521 }
1522 // and for safety
1523 bound_[2*numberRanges_]=bound_[2*numberRanges_-2];
1524 bound_[2*numberRanges_+1]=bound_[2*numberRanges_-1];
1525 for (i=1;i<numberRanges_;i++) {
1526 largestGap_ = CoinMax(largestGap_,bound_[2*i]-bound_[2*i-1]);
1527 }
1528 }
1529 delete [] sort;
1530 delete [] weight;
1531 range_=0;
1532}
1533
1534// Copy constructor
1535OsiLotsize::OsiLotsize ( const OsiLotsize & rhs)
1536 :OsiObject2(rhs)
1537
1538{
1539 columnNumber_ = rhs.columnNumber_;
1540 rangeType_ = rhs.rangeType_;
1541 numberRanges_ = rhs.numberRanges_;
1542 range_ = rhs.range_;
1543 largestGap_ = rhs.largestGap_;
1544 if (numberRanges_) {
1545 assert (rangeType_>0&&rangeType_<3);
1546 bound_= new double [(numberRanges_+1)*rangeType_];
1547 memcpy(bound_,rhs.bound_,(numberRanges_+1)*rangeType_*sizeof(double));
1548 } else {
1549 bound_=NULL;
1550 }
1551}
1552
1553// Clone
1554OsiObject *
1555OsiLotsize::clone() const
1556{
1557 return new OsiLotsize(*this);
1558}
1559
1560// Assignment operator
1561OsiLotsize &
1562OsiLotsize::operator=( const OsiLotsize& rhs)
1563{
1564 if (this!=&rhs) {
1565 OsiObject2::operator=(rhs);
1566 columnNumber_ = rhs.columnNumber_;
1567 rangeType_ = rhs.rangeType_;
1568 numberRanges_ = rhs.numberRanges_;
1569 largestGap_ = rhs.largestGap_;
1570 delete [] bound_;
1571 range_ = rhs.range_;
1572 if (numberRanges_) {
1573 assert (rangeType_>0&&rangeType_<3);
1574 bound_= new double [(numberRanges_+1)*rangeType_];
1575 memcpy(bound_,rhs.bound_,(numberRanges_+1)*rangeType_*sizeof(double));
1576 } else {
1577 bound_=NULL;
1578 }
1579 }
1580 return *this;
1581}
1582
1583// Destructor
1584OsiLotsize::~OsiLotsize ()
1585{
1586 delete [] bound_;
1587}
1588/* Finds range of interest so value is feasible in range range_ or infeasible
1589 between hi[range_] and lo[range_+1]. Returns true if feasible.
1590*/
1591bool
1592OsiLotsize::findRange(double value, double integerTolerance) const
1593{
1594 assert (range_>=0&&range_<numberRanges_+1);
1595 int iLo;
1596 int iHi;
1597 double infeasibility=0.0;
1598 if (rangeType_==1) {
1599 if (value<bound_[range_]-integerTolerance) {
1600 iLo=0;
1601 iHi=range_-1;
1602 } else if (value<bound_[range_]+integerTolerance) {
1603 return true;
1604 } else if (value<bound_[range_+1]-integerTolerance) {
1605 return false;
1606 } else {
1607 iLo=range_+1;
1608 iHi=numberRanges_-1;
1609 }
1610 // check lo and hi
1611 bool found=false;
1612 if (value>bound_[iLo]-integerTolerance&&value<bound_[iLo+1]+integerTolerance) {
1613 range_=iLo;
1614 found=true;
1615 } else if (value>bound_[iHi]-integerTolerance&&value<bound_[iHi+1]+integerTolerance) {
1616 range_=iHi;
1617 found=true;
1618 } else {
1619 range_ = (iLo+iHi)>>1;
1620 }
1621 //points
1622 while (!found) {
1623 if (value<bound_[range_]) {
1624 if (value>=bound_[range_-1]) {
1625 // found
1626 range_--;
1627 break;
1628 } else {
1629 iHi = range_;
1630 }
1631 } else {
1632 if (value<bound_[range_+1]) {
1633 // found
1634 break;
1635 } else {
1636 iLo = range_;
1637 }
1638 }
1639 range_ = (iLo+iHi)>>1;
1640 }
1641 if (value-bound_[range_]<=bound_[range_+1]-value) {
1642 infeasibility = value-bound_[range_];
1643 } else {
1644 infeasibility = bound_[range_+1]-value;
1645 if (infeasibility<integerTolerance)
1646 range_++;
1647 }
1648 return (infeasibility<integerTolerance);
1649 } else {
1650 // ranges
1651 if (value<bound_[2*range_]-integerTolerance) {
1652 iLo=0;
1653 iHi=range_-1;
1654 } else if (value<bound_[2*range_+1]+integerTolerance) {
1655 return true;
1656 } else if (value<bound_[2*range_+2]-integerTolerance) {
1657 return false;
1658 } else {
1659 iLo=range_+1;
1660 iHi=numberRanges_-1;
1661 }
1662 // check lo and hi
1663 bool found=false;
1664 if (value>bound_[2*iLo]-integerTolerance&&value<bound_[2*iLo+2]-integerTolerance) {
1665 range_=iLo;
1666 found=true;
1667 } else if (value>=bound_[2*iHi]-integerTolerance) {
1668 range_=iHi;
1669 found=true;
1670 } else {
1671 range_ = (iLo+iHi)>>1;
1672 }
1673 //points
1674 while (!found) {
1675 if (value<bound_[2*range_]) {
1676 if (value>=bound_[2*range_-2]) {
1677 // found
1678 range_--;
1679 break;
1680 } else {
1681 iHi = range_;
1682 }
1683 } else {
1684 if (value<bound_[2*range_+2]) {
1685 // found
1686 break;
1687 } else {
1688 iLo = range_;
1689 }
1690 }
1691 range_ = (iLo+iHi)>>1;
1692 }
1693 if (value>=bound_[2*range_]-integerTolerance&&value<=bound_[2*range_+1]+integerTolerance)
1694 infeasibility=0.0;
1695 else if (value-bound_[2*range_+1]<bound_[2*range_+2]-value) {
1696 infeasibility = value-bound_[2*range_+1];
1697 } else {
1698 infeasibility = bound_[2*range_+2]-value;
1699 }
1700 return (infeasibility<integerTolerance);
1701 }
1702}
1703/* Returns floor and ceiling
1704 */
1705void
1706OsiLotsize::floorCeiling(double & floorLotsize, double & ceilingLotsize, double value,
1707 double tolerance) const
1708{
1709 bool feasible=findRange(value,tolerance);
1710 if (rangeType_==1) {
1711 floorLotsize=bound_[range_];
1712 ceilingLotsize=bound_[range_+1];
1713 // may be able to adjust
1714 if (feasible&&fabs(value-floorLotsize)>fabs(value-ceilingLotsize)) {
1715 floorLotsize=bound_[range_+1];
1716 ceilingLotsize=bound_[range_+2];
1717 }
1718 } else {
1719 // ranges
1720 assert (value>=bound_[2*range_+1]);
1721 floorLotsize=bound_[2*range_+1];
1722 ceilingLotsize=bound_[2*range_+2];
1723 }
1724}
1725
1726// Infeasibility - large is 0.5
1727double
1728OsiLotsize::infeasibility(const OsiBranchingInformation * info, int & preferredWay) const
1729{
1730 const double * solution = info->solution_;
1731 const double * lower = info->lower_;
1732 const double * upper = info->upper_;
1733 double value = solution[columnNumber_];
1734 value = CoinMax(value, lower[columnNumber_]);
1735 value = CoinMin(value, upper[columnNumber_]);
1736 double integerTolerance = info->integerTolerance_;
1737 /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1738 solution[columnNumber_],upper[columnNumber_]);*/
1739 assert (value>=bound_[0]-integerTolerance
1740 &&value<=bound_[rangeType_*numberRanges_-1]+integerTolerance);
1741 infeasibility_=0.0;
1742 bool feasible = findRange(value,integerTolerance);
1743 if (!feasible) {
1744 if (rangeType_==1) {
1745 if (value-bound_[range_]<bound_[range_+1]-value) {
1746 preferredWay=-1;
1747 infeasibility_ = value-bound_[range_];
1748 otherInfeasibility_ = bound_[range_+1] - value ;
1749 } else {
1750 preferredWay=1;
1751 infeasibility_ = bound_[range_+1]-value;
1752 otherInfeasibility_ = value-bound_[range_];
1753 }
1754 } else {
1755 // ranges
1756 if (value-bound_[2*range_+1]<bound_[2*range_+2]-value) {
1757 preferredWay=-1;
1758 infeasibility_ = value-bound_[2*range_+1];
1759 otherInfeasibility_ = bound_[2*range_+2]-value;
1760 } else {
1761 preferredWay=1;
1762 infeasibility_ = bound_[2*range_+2]-value;
1763 otherInfeasibility_ = value-bound_[2*range_+1];
1764 }
1765 }
1766 } else {
1767 // always satisfied
1768 preferredWay=-1;
1769 otherInfeasibility_ = 1.0;
1770 }
1771 if (infeasibility_<integerTolerance)
1772 infeasibility_=0.0;
1773 else
1774 infeasibility_ /= largestGap_;
1775 return infeasibility_;
1776}
1777/* Column number if single column object -1 otherwise,
1778 so returns >= 0
1779 Used by heuristics
1780*/
1781int
1782OsiLotsize::columnNumber() const
1783{
1784 return columnNumber_;
1785}
1786/* Set bounds to contain the current solution.
1787 More precisely, for the variable associated with this object, take the
1788 value given in the current solution, force it within the current bounds
1789 if required, then set the bounds to fix the variable at the integer
1790 nearest the solution value. Returns amount it had to move variable.
1791*/
1792double
1793OsiLotsize::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
1794{
1795 const double * lower = solver->getColLower();
1796 const double * upper = solver->getColUpper();
1797 const double * solution = info->solution_;
1798 double value = solution[columnNumber_];
1799 value = CoinMax(value, lower[columnNumber_]);
1800 value = CoinMin(value, upper[columnNumber_]);
1801 findRange(value,info->integerTolerance_);
1802 double nearest;
1803 if (rangeType_==1) {
1804 nearest = bound_[range_];
1805 solver->setColLower(columnNumber_,nearest);
1806 solver->setColUpper(columnNumber_,nearest);
1807 } else {
1808 // ranges
1809 solver->setColLower(columnNumber_,bound_[2*range_]);
1810 solver->setColUpper(columnNumber_,bound_[2*range_+1]);
1811 if (value>bound_[2*range_+1])
1812 nearest=bound_[2*range_+1];
1813 else if (value<bound_[2*range_])
1814 nearest = bound_[2*range_];
1815 else
1816 nearest = value;
1817 }
1818 // Scaling may have moved it a bit
1819 // Lotsizing variables could be a lot larger
1820#ifndef NDEBUG
1821 assert (fabs(value-nearest)<=(100.0+10.0*fabs(nearest))*info->integerTolerance_);
1822#endif
1823 return fabs(value-nearest);
1824}
1825
1826// Creates a branching object
1827// Creates a branching object
1828OsiBranchingObject *
1829OsiLotsize::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
1830{
1831 const double * solution = info->solution_;
1832 const double * lower = solver->getColLower();
1833 const double * upper = solver->getColUpper();
1834 double value = solution[columnNumber_];
1835 value = CoinMax(value, lower[columnNumber_]);
1836 value = CoinMin(value, upper[columnNumber_]);
1837 assert (!findRange(value,info->integerTolerance_));
1838 return new OsiLotsizeBranchingObject(solver,this,way,
1839 value);
1840}
1841
1842
1843/*
1844 Bounds may be tightened, so it may be good to be able to refresh the local
1845 copy of the original bounds.
1846 */
1847void
1848OsiLotsize::resetBounds(const OsiSolverInterface * )
1849{
1850}
1851// Return "down" estimate
1852double
1853OsiLotsize::downEstimate() const
1854{
1855 if (whichWay_)
1856 return otherInfeasibility_;
1857 else
1858 return infeasibility_;
1859}
1860// Return "up" estimate
1861double
1862OsiLotsize::upEstimate() const
1863{
1864 if (!whichWay_)
1865 return otherInfeasibility_;
1866 else
1867 return infeasibility_;
1868}
1869// Redoes data when sequence numbers change
1870void
1871OsiLotsize::resetSequenceEtc(int numberColumns, const int * originalColumns)
1872{
1873 int i;
1874 for (i=0;i<numberColumns;i++) {
1875 if (originalColumns[i]==columnNumber_)
1876 break;
1877 }
1878 if (i<numberColumns)
1879 columnNumber_=i;
1880 else
1881 abort(); // should never happen
1882}
1883
1884
1885// Default Constructor
1886OsiLotsizeBranchingObject::OsiLotsizeBranchingObject()
1887 :OsiTwoWayBranchingObject()
1888{
1889 down_[0] = 0.0;
1890 down_[1] = 0.0;
1891 up_[0] = 0.0;
1892 up_[1] = 0.0;
1893}
1894
1895// Useful constructor
1896OsiLotsizeBranchingObject::OsiLotsizeBranchingObject (OsiSolverInterface * solver,
1897 const OsiLotsize * originalObject,
1898 int way , double value)
1899 :OsiTwoWayBranchingObject(solver,originalObject,way,value)
1900{
1901 int iColumn = originalObject->columnNumber();
1902 down_[0] = solver->getColLower()[iColumn];
1903 double integerTolerance = solver->getIntegerTolerance();
1904 originalObject->floorCeiling(down_[1],up_[0],value,integerTolerance);
1905 up_[1] = solver->getColUpper()[iColumn];
1906}
1907
1908// Copy constructor
1909OsiLotsizeBranchingObject::OsiLotsizeBranchingObject ( const OsiLotsizeBranchingObject & rhs) :OsiTwoWayBranchingObject(rhs)
1910{
1911 down_[0] = rhs.down_[0];
1912 down_[1] = rhs.down_[1];
1913 up_[0] = rhs.up_[0];
1914 up_[1] = rhs.up_[1];
1915}
1916
1917// Assignment operator
1918OsiLotsizeBranchingObject &
1919OsiLotsizeBranchingObject::operator=( const OsiLotsizeBranchingObject& rhs)
1920{
1921 if (this != &rhs) {
1922 OsiTwoWayBranchingObject::operator=(rhs);
1923 down_[0] = rhs.down_[0];
1924 down_[1] = rhs.down_[1];
1925 up_[0] = rhs.up_[0];
1926 up_[1] = rhs.up_[1];
1927 }
1928 return *this;
1929}
1930OsiBranchingObject *
1931OsiLotsizeBranchingObject::clone() const
1932{
1933 return (new OsiLotsizeBranchingObject(*this));
1934}
1935
1936
1937// Destructor
1938OsiLotsizeBranchingObject::~OsiLotsizeBranchingObject ()
1939{
1940}
1941
1942/*
1943 Perform a branch by adjusting the bounds of the specified variable. Note
1944 that each arm of the branch advances the object to the next arm by
1945 advancing the value of way_.
1946
1947 Providing new values for the variable's lower and upper bounds for each
1948 branching direction gives a little bit of additional flexibility and will
1949 be easily extensible to multi-way branching.
1950*/
1951double
1952OsiLotsizeBranchingObject::branch(OsiSolverInterface * solver)
1953{
1954 const OsiLotsize * obj =
1955 dynamic_cast <const OsiLotsize *>(originalObject_) ;
1956 assert (obj);
1957 int iColumn = obj->columnNumber();
1958 int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
1959 if (way<0) {
1960#ifdef OSI_DEBUG
1961 { double olb,oub ;
1962 olb = solver->getColLower()[iColumn] ;
1963 oub = solver->getColUpper()[iColumn] ;
1964 printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
1965 iColumn,olb,oub,down_[0],down_[1]) ; }
1966#endif
1967 solver->setColLower(iColumn,down_[0]);
1968 solver->setColUpper(iColumn,down_[1]);
1969 } else {
1970#ifdef OSI_DEBUG
1971 { double olb,oub ;
1972 olb = solver->getColLower()[iColumn] ;
1973 oub = solver->getColUpper()[iColumn] ;
1974 printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
1975 iColumn,olb,oub,up_[0],up_[1]) ; }
1976#endif
1977 solver->setColLower(iColumn,up_[0]);
1978 solver->setColUpper(iColumn,up_[1]);
1979 }
1980 branchIndex_++;
1981 return 0.0;
1982}
1983// Print
1984void
1985OsiLotsizeBranchingObject::print(const OsiSolverInterface * solver)
1986{
1987 const OsiLotsize * obj =
1988 dynamic_cast <const OsiLotsize *>(originalObject_) ;
1989 assert (obj);
1990 int iColumn = obj->columnNumber();
1991 int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
1992 if (way<0) {
1993 { double olb,oub ;
1994 olb = solver->getColLower()[iColumn] ;
1995 oub = solver->getColUpper()[iColumn] ;
1996 printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
1997 iColumn,olb,oub,down_[0],down_[1]) ; }
1998 } else {
1999 { double olb,oub ;
2000 olb = solver->getColLower()[iColumn] ;
2001 oub = solver->getColUpper()[iColumn] ;
2002 printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
2003 iColumn,olb,oub,up_[0],up_[1]) ; }
2004 }
2005}
2006