source: trunk/Cbc/examples/CbcBranchLink.cpp @ 2496

Last change on this file since 2496 was 2469, checked in by unxusr, 8 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.2 KB
Line 
1// $Id: CbcBranchLink.cpp 2469 2019-01-06 23:17:46Z forrest $
2// Copyright (C) 2005, 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 <cassert>
7#include <cmath>
8#include <cfloat>
9//#define CBC_DEBUG
10
11#include "CoinPragma.hpp"
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcBranchLink.hpp"
16#include "CoinError.hpp"
17#include "CoinPackedMatrix.hpp"
18
19// Default Constructor
20CbcLink::CbcLink()
21  : CbcObject()
22  , weights_(NULL)
23  , numberMembers_(0)
24  , numberLinks_(0)
25  , which_(NULL)
26  , sosType_(1)
27{
28}
29
30// Useful constructor (which are indices)
31CbcLink::CbcLink(CbcModel *model, int numberMembers,
32  int numberLinks, int first, const double *weights, int identifier)
33  : CbcObject(model)
34  , numberMembers_(numberMembers)
35  , numberLinks_(numberLinks)
36  , which_(NULL)
37  , sosType_(1)
38{
39  id_ = identifier;
40  if (numberMembers_) {
41    weights_ = new double[numberMembers_];
42    which_ = new int[numberMembers_ * numberLinks_];
43    if (weights) {
44      memcpy(weights_, weights, numberMembers_ * sizeof(double));
45    } else {
46      for (int i = 0; i < numberMembers_; i++)
47        weights_[i] = i;
48    }
49    // weights must be increasing
50    int i;
51    for (i = 0; i < numberMembers_; i++)
52      assert(i == 0 || weights_[i] > weights_[i - 1] + 1.0e-12);
53    for (i = 0; i < numberMembers_ * numberLinks_; i++) {
54      which_[i] = first + i;
55    }
56  } else {
57    weights_ = NULL;
58  }
59}
60
61// Useful constructor (which are indices)
62CbcLink::CbcLink(CbcModel *model, int numberMembers,
63  int numberLinks, int sosType, const int *which, const double *weights, int identifier)
64  : CbcObject(model)
65  , numberMembers_(numberMembers)
66  , numberLinks_(numberLinks)
67  , which_(NULL)
68  , sosType_(sosType)
69{
70  id_ = identifier;
71  if (numberMembers_) {
72    weights_ = new double[numberMembers_];
73    which_ = new int[numberMembers_ * numberLinks_];
74    if (weights) {
75      memcpy(weights_, weights, numberMembers_ * sizeof(double));
76    } else {
77      for (int i = 0; i < numberMembers_; i++)
78        weights_[i] = i;
79    }
80    // weights must be increasing
81    int i;
82    for (i = 0; i < numberMembers_; i++)
83      assert(i == 0 || weights_[i] > weights_[i - 1] + 1.0e-12);
84    for (i = 0; i < numberMembers_ * numberLinks_; i++) {
85      which_[i] = which[i];
86    }
87  } else {
88    weights_ = NULL;
89  }
90}
91
92// Copy constructor
93CbcLink::CbcLink(const CbcLink &rhs)
94  : CbcObject(rhs)
95{
96  numberMembers_ = rhs.numberMembers_;
97  numberLinks_ = rhs.numberLinks_;
98  sosType_ = rhs.sosType_;
99  if (numberMembers_) {
100    weights_ = CoinCopyOfArray(rhs.weights_, numberMembers_);
101    which_ = CoinCopyOfArray(rhs.which_, numberMembers_ * numberLinks_);
102  } else {
103    weights_ = NULL;
104    which_ = NULL;
105  }
106}
107
108// Clone
109CbcObject *
110CbcLink::clone() const
111{
112  return new CbcLink(*this);
113}
114
115// Assignment operator
116CbcLink &
117CbcLink::operator=(const CbcLink &rhs)
118{
119  if (this != &rhs) {
120    CbcObject::operator=(rhs);
121    delete[] weights_;
122    delete[] which_;
123    numberMembers_ = rhs.numberMembers_;
124    numberLinks_ = rhs.numberLinks_;
125    sosType_ = rhs.sosType_;
126    if (numberMembers_) {
127      weights_ = CoinCopyOfArray(rhs.weights_, numberMembers_);
128      which_ = CoinCopyOfArray(rhs.which_, numberMembers_ * numberLinks_);
129    } else {
130      weights_ = NULL;
131      which_ = NULL;
132    }
133  }
134  return *this;
135}
136
137// Destructor
138CbcLink::~CbcLink()
139{
140  delete[] weights_;
141  delete[] which_;
142}
143
144// Infeasibility - large is 0.5
145double
146CbcLink::infeasibility(int &preferredWay) const
147{
148  int j;
149  int firstNonZero = -1;
150  int lastNonZero = -1;
151  OsiSolverInterface *solver = model_->solver();
152  const double *solution = model_->testSolution();
153  //const double * lower = solver->getColLower();
154  const double *upper = solver->getColUpper();
155  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
156  double weight = 0.0;
157  double sum = 0.0;
158
159  // check bounds etc
160  double lastWeight = -1.0e100;
161  int base = 0;
162  for (j = 0; j < numberMembers_; j++) {
163    for (int k = 0; k < numberLinks_; k++) {
164      int iColumn = which_[base + k];
165      //if (lower[iColumn])
166      //throw CoinError("Non zero lower bound in CBCLink","infeasibility","CbcLink");
167      if (lastWeight >= weights_[j] - 1.0e-7)
168        throw CoinError("Weights too close together in CBCLink", "infeasibility", "CbcLink");
169      double value = CoinMax(0.0, solution[iColumn]);
170      sum += value;
171      if (value > integerTolerance && upper[iColumn]) {
172        // Possibly due to scaling a fixed variable might slip through
173        if (value > upper[iColumn] + 1.0e-8) {
174          // Could change to #ifdef CBC_DEBUG
175#ifndef NDEBUG
176          if (model_->messageHandler()->logLevel() > 1)
177            printf("** Variable %d (%d) has value %g and upper bound of %g\n",
178              iColumn, j, value, upper[iColumn]);
179#endif
180        }
181        value = CoinMin(value, upper[iColumn]);
182        weight += weights_[j] * value;
183        if (firstNonZero < 0)
184          firstNonZero = j;
185        lastNonZero = j;
186      }
187    }
188    base += numberLinks_;
189  }
190  double valueInfeasibility;
191  preferredWay = 1;
192  if (lastNonZero - firstNonZero >= sosType_) {
193    // find where to branch
194    assert(sum > 0.0);
195    weight /= sum;
196    valueInfeasibility = lastNonZero - firstNonZero + 1;
197    valueInfeasibility *= 0.5 / ((double)numberMembers_);
198    //#define DISTANCE
199#ifdef DISTANCE
200    assert(sosType_ == 1); // code up
201    /* may still be satisfied.
202       For LOS type 2 we might wish to move coding around
203       and keep initial info in model_ for speed
204    */
205    int iWhere;
206    bool possible = false;
207    for (iWhere = firstNonZero; iWhere <= lastNonZero; iWhere++) {
208      if (fabs(weight - weights_[iWhere]) < 1.0e-8) {
209        possible = true;
210        break;
211      }
212    }
213    if (possible) {
214      // One could move some of this (+ arrays) into model_
215      const CoinPackedMatrix *matrix = solver->getMatrixByCol();
216      const double *element = matrix->getMutableElements();
217      const int *row = matrix->getIndices();
218      const CoinBigIndex *columnStart = matrix->getVectorStarts();
219      const int *columnLength = matrix->getVectorLengths();
220      const double *rowSolution = solver->getRowActivity();
221      const double *rowLower = solver->getRowLower();
222      const double *rowUpper = solver->getRowUpper();
223      int numberRows = matrix->getNumRows();
224      double *array = new double[numberRows];
225      CoinZeroN(array, numberRows);
226      int *which = new int[numberRows];
227      int n = 0;
228      int base = numberLinks_ * firstNonZero;
229      for (j = firstNonZero; j <= lastNonZero; j++) {
230        for (int k = 0; k < numberLinks_; k++) {
231          int iColumn = which_[base + k];
232          double value = CoinMax(0.0, solution[iColumn]);
233          if (value > integerTolerance && upper[iColumn]) {
234            value = CoinMin(value, upper[iColumn]);
235            for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
236              int iRow = row[j];
237              double a = array[iRow];
238              if (a) {
239                a += value * element[j];
240                if (!a)
241                  a = 1.0e-100;
242              } else {
243                which[n++] = iRow;
244                a = value * element[j];
245                assert(a);
246              }
247              array[iRow] = a;
248            }
249          }
250        }
251        base += numberLinks_;
252      }
253      base = numberLinks_ * iWhere;
254      for (int k = 0; k < numberLinks_; k++) {
255        int iColumn = which_[base + k];
256        const double value = 1.0;
257        for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
258          int iRow = row[j];
259          double a = array[iRow];
260          if (a) {
261            a -= value * element[j];
262            if (!a)
263              a = 1.0e-100;
264          } else {
265            which[n++] = iRow;
266            a = -value * element[j];
267            assert(a);
268          }
269          array[iRow] = a;
270        }
271      }
272      for (j = 0; j < n; j++) {
273        int iRow = which[j];
274        // moving to point will increase row solution by this
275        double distance = array[iRow];
276        if (distance > 1.0e-8) {
277          if (distance + rowSolution[iRow] > rowUpper[iRow] + 1.0e-8) {
278            possible = false;
279            break;
280          }
281        } else if (distance < -1.0e-8) {
282          if (distance + rowSolution[iRow] < rowLower[iRow] - 1.0e-8) {
283            possible = false;
284            break;
285          }
286        }
287      }
288      for (j = 0; j < n; j++)
289        array[which[j]] = 0.0;
290      delete[] array;
291      delete[] which;
292      if (possible) {
293        valueInfeasibility = 0.0;
294        printf("possible %d %d %d\n", firstNonZero, lastNonZero, iWhere);
295      }
296    }
297#endif
298  } else {
299    valueInfeasibility = 0.0; // satisfied
300  }
301  return valueInfeasibility;
302}
303
304// This looks at solution and sets bounds to contain solution
305void CbcLink::feasibleRegion()
306{
307  int j;
308  int firstNonZero = -1;
309  int lastNonZero = -1;
310  OsiSolverInterface *solver = model_->solver();
311  const double *solution = model_->testSolution();
312  const double *upper = solver->getColUpper();
313  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
314  double weight = 0.0;
315  double sum = 0.0;
316
317  int base = 0;
318  for (j = 0; j < numberMembers_; j++) {
319    for (int k = 0; k < numberLinks_; k++) {
320      int iColumn = which_[base + k];
321      double value = CoinMax(0.0, solution[iColumn]);
322      sum += value;
323      if (value > integerTolerance && upper[iColumn]) {
324        weight += weights_[j] * value;
325        if (firstNonZero < 0)
326          firstNonZero = j;
327        lastNonZero = j;
328      }
329    }
330    base += numberLinks_;
331  }
332#ifdef DISTANCE
333  if (lastNonZero - firstNonZero > sosType_ - 1) {
334    /* may still be satisfied.
335       For LOS type 2 we might wish to move coding around
336       and keep initial info in model_ for speed
337    */
338    int iWhere;
339    bool possible = false;
340    for (iWhere = firstNonZero; iWhere <= lastNonZero; iWhere++) {
341      if (fabs(weight - weights_[iWhere]) < 1.0e-8) {
342        possible = true;
343        break;
344      }
345    }
346    if (possible) {
347      // One could move some of this (+ arrays) into model_
348      const CoinPackedMatrix *matrix = solver->getMatrixByCol();
349      const double *element = matrix->getMutableElements();
350      const int *row = matrix->getIndices();
351      const CoinBigIndex *columnStart = matrix->getVectorStarts();
352      const int *columnLength = matrix->getVectorLengths();
353      const double *rowSolution = solver->getRowActivity();
354      const double *rowLower = solver->getRowLower();
355      const double *rowUpper = solver->getRowUpper();
356      int numberRows = matrix->getNumRows();
357      double *array = new double[numberRows];
358      CoinZeroN(array, numberRows);
359      int *which = new int[numberRows];
360      int n = 0;
361      int base = numberLinks_ * firstNonZero;
362      for (j = firstNonZero; j <= lastNonZero; j++) {
363        for (int k = 0; k < numberLinks_; k++) {
364          int iColumn = which_[base + k];
365          double value = CoinMax(0.0, solution[iColumn]);
366          if (value > integerTolerance && upper[iColumn]) {
367            value = CoinMin(value, upper[iColumn]);
368            for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
369              int iRow = row[j];
370              double a = array[iRow];
371              if (a) {
372                a += value * element[j];
373                if (!a)
374                  a = 1.0e-100;
375              } else {
376                which[n++] = iRow;
377                a = value * element[j];
378                assert(a);
379              }
380              array[iRow] = a;
381            }
382          }
383        }
384        base += numberLinks_;
385      }
386      base = numberLinks_ * iWhere;
387      for (int k = 0; k < numberLinks_; k++) {
388        int iColumn = which_[base + k];
389        const double value = 1.0;
390        for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
391          int iRow = row[j];
392          double a = array[iRow];
393          if (a) {
394            a -= value * element[j];
395            if (!a)
396              a = 1.0e-100;
397          } else {
398            which[n++] = iRow;
399            a = -value * element[j];
400            assert(a);
401          }
402          array[iRow] = a;
403        }
404      }
405      for (j = 0; j < n; j++) {
406        int iRow = which[j];
407        // moving to point will increase row solution by this
408        double distance = array[iRow];
409        if (distance > 1.0e-8) {
410          if (distance + rowSolution[iRow] > rowUpper[iRow] + 1.0e-8) {
411            possible = false;
412            break;
413          }
414        } else if (distance < -1.0e-8) {
415          if (distance + rowSolution[iRow] < rowLower[iRow] - 1.0e-8) {
416            possible = false;
417            break;
418          }
419        }
420      }
421      for (j = 0; j < n; j++)
422        array[which[j]] = 0.0;
423      delete[] array;
424      delete[] which;
425      if (possible) {
426        printf("possible feas region %d %d %d\n", firstNonZero, lastNonZero, iWhere);
427        firstNonZero = iWhere;
428        lastNonZero = iWhere;
429      }
430    }
431  }
432#else
433  assert(lastNonZero - firstNonZero < sosType_);
434#endif
435  base = 0;
436  for (j = 0; j < firstNonZero; j++) {
437    for (int k = 0; k < numberLinks_; k++) {
438      int iColumn = which_[base + k];
439      solver->setColUpper(iColumn, 0.0);
440    }
441    base += numberLinks_;
442  }
443  // skip
444  base += numberLinks_;
445  for (j = lastNonZero + 1; j < numberMembers_; j++) {
446    for (int k = 0; k < numberLinks_; k++) {
447      int iColumn = which_[base + k];
448      solver->setColUpper(iColumn, 0.0);
449    }
450    base += numberLinks_;
451  }
452}
453
454// Creates a branching object
455CbcBranchingObject *
456CbcLink::createCbcBranch(OsiSolverInterface * /*solver*/, const OsiBranchingInformation * /*info*/, int way)
457{
458  int j;
459  const double *solution = model_->testSolution();
460  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
461  OsiSolverInterface *solver = model_->solver();
462  const double *upper = solver->getColUpper();
463  int firstNonFixed = -1;
464  int lastNonFixed = -1;
465  int firstNonZero = -1;
466  int lastNonZero = -1;
467  double weight = 0.0;
468  double sum = 0.0;
469  int base = 0;
470  for (j = 0; j < numberMembers_; j++) {
471    for (int k = 0; k < numberLinks_; k++) {
472      int iColumn = which_[base + k];
473      if (upper[iColumn]) {
474        double value = CoinMax(0.0, solution[iColumn]);
475        sum += value;
476        if (firstNonFixed < 0)
477          firstNonFixed = j;
478        lastNonFixed = j;
479        if (value > integerTolerance) {
480          weight += weights_[j] * value;
481          if (firstNonZero < 0)
482            firstNonZero = j;
483          lastNonZero = j;
484        }
485      }
486    }
487    base += numberLinks_;
488  }
489  assert(lastNonZero - firstNonZero >= sosType_);
490  // find where to branch
491  assert(sum > 0.0);
492  weight /= sum;
493  int iWhere;
494  double separator = 0.0;
495  for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
496    if (weight < weights_[iWhere + 1])
497      break;
498  if (sosType_ == 1) {
499    // SOS 1
500    separator = 0.5 * (weights_[iWhere] + weights_[iWhere + 1]);
501  } else {
502    // SOS 2
503    if (iWhere == firstNonFixed)
504      iWhere++;
505    ;
506    if (iWhere == lastNonFixed - 1)
507      iWhere = lastNonFixed - 2;
508    separator = weights_[iWhere + 1];
509  }
510  // create object
511  CbcBranchingObject *branch;
512  branch = new CbcLinkBranchingObject(model_, this, way, separator);
513  return branch;
514}
515// Useful constructor
516CbcLinkBranchingObject::CbcLinkBranchingObject(CbcModel *model,
517  const CbcLink *set,
518  int way,
519  double separator)
520  : CbcBranchingObject(model, set->id(), way, 0.5)
521{
522  set_ = set;
523  separator_ = separator;
524}
525
526// Copy constructor
527CbcLinkBranchingObject::CbcLinkBranchingObject(const CbcLinkBranchingObject &rhs)
528  : CbcBranchingObject(rhs)
529{
530  set_ = rhs.set_;
531  separator_ = rhs.separator_;
532}
533
534// Assignment operator
535CbcLinkBranchingObject &
536CbcLinkBranchingObject::operator=(const CbcLinkBranchingObject &rhs)
537{
538  if (this != &rhs) {
539    CbcBranchingObject::operator=(rhs);
540    set_ = rhs.set_;
541    separator_ = rhs.separator_;
542  }
543  return *this;
544}
545CbcBranchingObject *
546CbcLinkBranchingObject::clone() const
547{
548  return (new CbcLinkBranchingObject(*this));
549}
550
551// Destructor
552CbcLinkBranchingObject::~CbcLinkBranchingObject()
553{
554}
555double
556CbcLinkBranchingObject::branch()
557{
558  decrementNumberBranchesLeft();
559  int numberMembers = set_->numberMembers();
560  int numberLinks = set_->numberLinks();
561  const double *weights = set_->weights();
562  const int *which = set_->which();
563  OsiSolverInterface *solver = model_->solver();
564  //const double * lower = solver->getColLower();
565  //const double * upper = solver->getColUpper();
566  // *** for way - up means fix all those in down section
567  if (way_ < 0) {
568    int i;
569    for (i = 0; i < numberMembers; i++) {
570      if (weights[i] > separator_)
571        break;
572    }
573    assert(i < numberMembers);
574    int base = i * numberLinks;
575    ;
576    for (; i < numberMembers; i++) {
577      for (int k = 0; k < numberLinks; k++) {
578        int iColumn = which[base + k];
579        solver->setColUpper(iColumn, 0.0);
580      }
581      base += numberLinks;
582    }
583    way_ = 1; // Swap direction
584  } else {
585    int i;
586    int base = 0;
587    for (i = 0; i < numberMembers; i++) {
588      if (weights[i] >= separator_) {
589        break;
590      } else {
591        for (int k = 0; k < numberLinks; k++) {
592          int iColumn = which[base + k];
593          solver->setColUpper(iColumn, 0.0);
594        }
595        base += numberLinks;
596      }
597    }
598    assert(i < numberMembers);
599    way_ = -1; // Swap direction
600  }
601  return 0.0;
602}
603// Print what would happen
604void CbcLinkBranchingObject::print()
605{
606  int numberMembers = set_->numberMembers();
607  int numberLinks = set_->numberLinks();
608  const double *weights = set_->weights();
609  const int *which = set_->which();
610  OsiSolverInterface *solver = model_->solver();
611  const double *upper = solver->getColUpper();
612  int first = numberMembers;
613  int last = -1;
614  int numberFixed = 0;
615  int numberOther = 0;
616  int i;
617  int base = 0;
618  for (i = 0; i < numberMembers; i++) {
619    for (int k = 0; k < numberLinks; k++) {
620      int iColumn = which[base + k];
621      double bound = upper[iColumn];
622      if (bound) {
623        first = CoinMin(first, i);
624        last = CoinMax(last, i);
625      }
626    }
627    base += numberLinks;
628  }
629  // *** for way - up means fix all those in down section
630  base = 0;
631  if (way_ < 0) {
632    printf("SOS Down");
633    for (i = 0; i < numberMembers; i++) {
634      if (weights[i] > separator_)
635        break;
636      for (int k = 0; k < numberLinks; k++) {
637        int iColumn = which[base + k];
638        double bound = upper[iColumn];
639        if (bound)
640          numberOther++;
641      }
642      base += numberLinks;
643    }
644    assert(i < numberMembers);
645    for (; i < numberMembers; i++) {
646      for (int k = 0; k < numberLinks; k++) {
647        int iColumn = which[base + k];
648        double bound = upper[iColumn];
649        if (bound)
650          numberFixed++;
651      }
652      base += numberLinks;
653    }
654  } else {
655    printf("SOS Up");
656    for (i = 0; i < numberMembers; i++) {
657      if (weights[i] >= separator_)
658        break;
659      for (int k = 0; k < numberLinks; k++) {
660        int iColumn = which[base + k];
661        double bound = upper[iColumn];
662        if (bound)
663          numberFixed++;
664      }
665      base += numberLinks;
666    }
667    assert(i < numberMembers);
668    for (; i < numberMembers; i++) {
669      for (int k = 0; k < numberLinks; k++) {
670        int iColumn = which[base + k];
671        double bound = upper[iColumn];
672        if (bound)
673          numberOther++;
674      }
675      base += numberLinks;
676    }
677  }
678  assert((numberFixed % numberLinks) == 0);
679  assert((numberOther % numberLinks) == 0);
680  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
681    separator_, first, weights[first], last, weights[last], numberFixed / numberLinks,
682    numberOther / numberLinks);
683}
684/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
685    same type and must have the same original object, but they may have
686    different feasible regions.
687    Return the appropriate CbcRangeCompare value (first argument being the
688    sub/superset if that's the case). In case of overlap (and if \c
689    replaceIfOverlap is true) replace the current branching object with one
690    whose feasible region is the overlap.
691*/
692CbcRangeCompare
693CbcLinkBranchingObject::compareBranchingObject(const CbcBranchingObject *brObj, const bool replaceIfOverlap)
694{
695  throw("must implement");
696}
Note: See TracBrowser for help on using the repository browser.