source: trunk/Cbc/src/CbcFathomDynamicProgramming.cpp @ 2464

Last change on this file since 2464 was 2464, checked in by unxusr, 9 months ago

.clang-format with proposal for formatting code

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.2 KB
Line 
1/*
2  $Id: CbcFathomDynamicProgramming.cpp 2464 2019-01-03 19:03:23Z unxusr $
3*/
4// Copyright (C) 2004, International Business Machines
5// Corporation and others.  All Rights Reserved.
6// This code is licensed under the terms of the Eclipse Public License (EPL).
7
8#if defined(_MSC_VER)
9// Turn off compiler warning about long names
10#pragma warning(disable : 4786)
11#endif
12#include <cassert>
13#include <cstdlib>
14#include <cmath>
15#include <cfloat>
16
17#include "OsiSolverInterface.hpp"
18#include "CbcModel.hpp"
19#include "CbcMessage.hpp"
20#include "CbcFathomDynamicProgramming.hpp"
21#include "CoinHelperFunctions.hpp"
22#include "CoinPackedMatrix.hpp"
23#include "CoinSort.hpp"
24// Default Constructor
25CbcFathomDynamicProgramming::CbcFathomDynamicProgramming()
26  : CbcFathom()
27  , size_(0)
28  , type_(-1)
29  , cost_(NULL)
30  , back_(NULL)
31  , lookup_(NULL)
32  , indices_(NULL)
33  , numberActive_(0)
34  , maximumSizeAllowed_(1000000)
35  , startBit_(NULL)
36  , numberBits_(NULL)
37  , rhs_(NULL)
38  , coefficients_(NULL)
39  , target_(0)
40  , numberNonOne_(0)
41  , bitPattern_(0)
42  , algorithm_(-1)
43{
44}
45
46// Constructor from model
47CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(CbcModel &model)
48  : CbcFathom(model)
49  , cost_(NULL)
50  , back_(NULL)
51  , lookup_(NULL)
52  , indices_(NULL)
53  , numberActive_(0)
54  , maximumSizeAllowed_(1000000)
55  , startBit_(NULL)
56  , numberBits_(NULL)
57  , rhs_(NULL)
58  , coefficients_(NULL)
59  , target_(0)
60  , numberNonOne_(0)
61  , bitPattern_(0)
62  , algorithm_(-1)
63{
64  type_ = checkPossible();
65}
66
67// Destructor
68CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming()
69{
70  gutsOfDelete();
71}
72// Does deleteions
73void CbcFathomDynamicProgramming::gutsOfDelete()
74{
75  delete[] cost_;
76  delete[] back_;
77  delete[] lookup_;
78  delete[] indices_;
79  delete[] startBit_;
80  delete[] numberBits_;
81  delete[] rhs_;
82  delete[] coefficients_;
83  cost_ = NULL;
84  back_ = NULL;
85  lookup_ = NULL;
86  indices_ = NULL;
87  startBit_ = NULL;
88  numberBits_ = NULL;
89  rhs_ = NULL;
90  coefficients_ = NULL;
91}
92// Clone
93CbcFathom *
94CbcFathomDynamicProgramming::clone() const
95{
96  return new CbcFathomDynamicProgramming(*this);
97}
98
99// Copy constructor
100CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(const CbcFathomDynamicProgramming &rhs)
101  : CbcFathom(rhs)
102  , size_(rhs.size_)
103  , type_(rhs.type_)
104  , cost_(NULL)
105  , back_(NULL)
106  , lookup_(NULL)
107  , indices_(NULL)
108  , numberActive_(rhs.numberActive_)
109  , maximumSizeAllowed_(rhs.maximumSizeAllowed_)
110  , startBit_(NULL)
111  , numberBits_(NULL)
112  , rhs_(NULL)
113  , coefficients_(NULL)
114  , target_(rhs.target_)
115  , numberNonOne_(rhs.numberNonOne_)
116  , bitPattern_(rhs.bitPattern_)
117  , algorithm_(rhs.algorithm_)
118{
119  if (size_) {
120    cost_ = CoinCopyOfArray(rhs.cost_, size_);
121    back_ = CoinCopyOfArray(rhs.back_, size_);
122    int numberRows = model_->getNumRows();
123    lookup_ = CoinCopyOfArray(rhs.lookup_, numberRows);
124    startBit_ = CoinCopyOfArray(rhs.startBit_, numberActive_);
125    indices_ = CoinCopyOfArray(rhs.indices_, numberActive_);
126    numberBits_ = CoinCopyOfArray(rhs.numberBits_, numberActive_);
127    rhs_ = CoinCopyOfArray(rhs.rhs_, numberActive_);
128    coefficients_ = CoinCopyOfArray(rhs.coefficients_, numberActive_);
129  }
130}
131// Returns type
132int CbcFathomDynamicProgramming::checkPossible(int allowableSize)
133{
134  algorithm_ = -1;
135  assert(model_->solver());
136  OsiSolverInterface *solver = model_->solver();
137  const CoinPackedMatrix *matrix = solver->getMatrixByCol();
138
139  int numberIntegers = model_->numberIntegers();
140  int numberColumns = solver->getNumCols();
141  size_ = 0;
142  if (numberIntegers != numberColumns)
143    return -1; // can't do dynamic programming
144
145  const double *lower = solver->getColLower();
146  const double *upper = solver->getColUpper();
147  const double *rowUpper = solver->getRowUpper();
148
149  int numberRows = model_->getNumRows();
150  int i;
151
152  // First check columns to see if possible
153  double *rhs = new double[numberRows];
154  CoinCopyN(rowUpper, numberRows, rhs);
155
156  // Column copy
157  const double *element = matrix->getElements();
158  const int *row = matrix->getIndices();
159  const CoinBigIndex *columnStart = matrix->getVectorStarts();
160  const int *columnLength = matrix->getVectorLengths();
161  bool bad = false;
162  /* It is just possible that we could say okay as
163       variables may get fixed but seems unlikely */
164  for (i = 0; i < numberColumns; i++) {
165    CoinBigIndex j;
166    double lowerValue = lower[i];
167    assert(lowerValue == floor(lowerValue));
168    for (j = columnStart[i];
169         j < columnStart[i] + columnLength[i]; j++) {
170      int iRow = row[j];
171      double value = element[j];
172      if (upper[i] > lowerValue && (value <= 0.0 || value != floor(value)))
173        bad = true;
174      if (lowerValue)
175        rhs[iRow] -= lowerValue * value;
176    }
177  }
178  // check possible (at present do not allow covering)
179  int numberActive = 0;
180  bool infeasible = false;
181  bool saveBad = bad;
182  for (i = 0; i < numberRows; i++) {
183    if (rhs[i] < 0)
184      infeasible = true;
185    else if (rhs[i] > 1.0e5 || fabs(rhs[i] - floor(rhs[i] + 0.5)) > 1.0e-7)
186      bad = true;
187    else if (rhs[i] > 0.0)
188      numberActive++;
189  }
190  if (bad || infeasible) {
191    delete[] rhs;
192    if (!saveBad && infeasible)
193      return -2;
194    else
195      return -1;
196  }
197  // check size of array needed
198  double size = 1.0;
199  double check = COIN_INT_MAX;
200  for (i = 0; i < numberRows; i++) {
201    int n = static_cast<int>(floor(rhs[i] + 0.5));
202    if (n) {
203      n++; // allow for 0,1... n
204      if (numberActive != 1) {
205        // power of 2
206        int iBit = 0;
207        int k = n;
208        k &= ~1;
209        while (k) {
210          iBit++;
211          k &= ~(1 << iBit);
212        }
213        // See if exact power
214        if (n != (1 << iBit)) {
215          // round up to next power of 2
216          n = 1 << (iBit + 1);
217        }
218        size *= n;
219        if (size >= check)
220          break;
221      } else {
222        size = n; // just one constraint
223      }
224    }
225  }
226  // set size needed
227  if (size >= check)
228    size_ = COIN_INT_MAX;
229  else
230    size_ = static_cast<int>(size);
231
232  int n01 = 0;
233  int nbadcoeff = 0;
234  // See if we can tighten bounds
235  for (i = 0; i < numberColumns; i++) {
236    CoinBigIndex j;
237    double lowerValue = lower[i];
238    double gap = upper[i] - lowerValue;
239    for (j = columnStart[i];
240         j < columnStart[i] + columnLength[i]; j++) {
241      int iRow = row[j];
242      double value = element[j];
243      if (value != 1.0)
244        nbadcoeff++;
245      if (gap * value > rhs[iRow] + 1.0e-8)
246        gap = rhs[iRow] / value;
247    }
248    gap = lowerValue + floor(gap + 1.0e-7);
249    if (gap < upper[i])
250      solver->setColUpper(i, gap);
251    if (gap <= 1.0)
252      n01++;
253  }
254  if (allowableSize && size_ <= allowableSize) {
255    if (n01 == numberColumns && !nbadcoeff)
256      algorithm_ = 0; // easiest
257    else
258      algorithm_ = 1;
259  }
260  if (allowableSize && size_ <= allowableSize) {
261    numberActive_ = numberActive;
262    indices_ = new int[numberActive_];
263    cost_ = new double[size_];
264    CoinFillN(cost_, size_, COIN_DBL_MAX);
265    // but do nothing is okay
266    cost_[0] = 0.0;
267    back_ = new int[size_];
268    CoinFillN(back_, size_, -1);
269    startBit_ = new int[numberActive_];
270    numberBits_ = new int[numberActive_];
271    lookup_ = new int[numberRows];
272    rhs_ = new int[numberActive_];
273    numberActive = 0;
274    int kBit = 0;
275    for (i = 0; i < numberRows; i++) {
276      int n = static_cast<int>(floor(rhs[i] + 0.5));
277      if (n) {
278        lookup_[i] = numberActive;
279        rhs_[numberActive] = n;
280        startBit_[numberActive] = kBit;
281        n++; // allow for 0,1... n
282        int iBit = 0;
283        // power of 2
284        int k = n;
285        k &= ~1;
286        while (k) {
287          iBit++;
288          k &= ~(1 << iBit);
289        }
290        // See if exact power
291        if (n != (1 << iBit)) {
292          // round up to next power of 2
293          iBit++;
294        }
295        if (numberActive != 1) {
296          n = 1 << iBit;
297          size *= n;
298          if (size >= check)
299            break;
300        } else {
301          size = n; // just one constraint
302        }
303        numberBits_[numberActive++] = iBit;
304        kBit += iBit;
305      } else {
306        lookup_[i] = -1;
307      }
308    }
309    const double *rowLower = solver->getRowLower();
310    if (algorithm_ == 0) {
311      // rhs 1 and coefficients 1
312      // Get first possible solution for printing
313      target_ = -1;
314      int needed = 0;
315      int numberActive = 0;
316      for (i = 0; i < numberRows; i++) {
317        int newRow = lookup_[i];
318        if (newRow >= 0) {
319          if (rowLower[i] == rowUpper[i]) {
320            needed += 1 << numberActive;
321            numberActive++;
322          }
323        }
324      }
325      for (i = 0; i < size_; i++) {
326        if ((i & needed) == needed) {
327          break;
328        }
329      }
330      target_ = i;
331    } else {
332      coefficients_ = new int[numberActive_];
333      // If not too many general rhs then we can be more efficient
334      numberNonOne_ = 0;
335      for (i = 0; i < numberActive_; i++) {
336        if (rhs_[i] != 1)
337          numberNonOne_++;
338      }
339      if (numberNonOne_ * 2 < numberActive_) {
340        // put rhs >1 every second
341        int *permute = new int[numberActive_];
342        int *temp = new int[numberActive_];
343        // try different ways
344        int k = 0;
345        for (i = 0; i < numberRows; i++) {
346          int newRow = lookup_[i];
347          if (newRow >= 0 && rhs_[newRow] > 1) {
348            permute[newRow] = k;
349            k += 2;
350          }
351        }
352        // adjust so k points to last
353        k -= 2;
354        // and now rest
355        int k1 = 1;
356        for (i = 0; i < numberRows; i++) {
357          int newRow = lookup_[i];
358          if (newRow >= 0 && rhs_[newRow] == 1) {
359            permute[newRow] = k1;
360            k1++;
361            if (k1 <= k)
362              k1++;
363          }
364        }
365        for (i = 0; i < numberActive_; i++) {
366          int put = permute[i];
367          temp[put] = rhs_[i];
368        }
369        memcpy(rhs_, temp, numberActive_ * sizeof(int));
370        for (i = 0; i < numberActive_; i++) {
371          int put = permute[i];
372          temp[put] = numberBits_[i];
373        }
374        memcpy(numberBits_, temp, numberActive_ * sizeof(int));
375        k = 0;
376        for (i = 0; i < numberActive_; i++) {
377          startBit_[i] = k;
378          k += numberBits_[i];
379        }
380        for (i = 0; i < numberRows; i++) {
381          int newRow = lookup_[i];
382          if (newRow >= 0)
383            lookup_[i] = permute[newRow];
384        }
385        delete[] permute;
386        delete[] temp;
387        // mark new method
388        algorithm_ = 2;
389      }
390      // Get first possible solution for printing
391      target_ = -1;
392      int needed = 0;
393      int *lower2 = new int[numberActive_];
394      for (i = 0; i < numberRows; i++) {
395        int newRow = lookup_[i];
396        if (newRow >= 0) {
397          int gap = static_cast<int>(rowUpper[i] - CoinMax(0.0, rowLower[i]));
398          lower2[newRow] = rhs_[newRow] - gap;
399          int numberBits = numberBits_[newRow];
400          int startBit = startBit_[newRow];
401          if (numberBits == 1 && !gap) {
402            needed |= 1 << startBit;
403          }
404        }
405      }
406      for (i = 0; i < size_; i++) {
407        if ((i & needed) == needed) {
408          // this one may do
409          bool good = true;
410          for (int kk = 0; kk < numberActive_; kk++) {
411            int numberBits = numberBits_[kk];
412            int startBit = startBit_[kk];
413            int size = 1 << numberBits;
414            int start = 1 << startBit;
415            int mask = start * (size - 1);
416            int level = (i & mask) >> startBit;
417            if (level < lower2[kk]) {
418              good = false;
419              break;
420            }
421          }
422          if (good) {
423            break;
424          }
425        }
426      }
427      delete[] lower2;
428      target_ = i;
429    }
430  }
431  delete[] rhs;
432  if (allowableSize && size_ > allowableSize) {
433    COIN_DETAIL_PRINT(printf("Too large - need %d entries x 8 bytes\n", size_));
434    return -1; // too big
435  } else {
436    return algorithm_;
437  }
438}
439
440// Resets stuff if model changes
441void CbcFathomDynamicProgramming::resetModel(CbcModel *model)
442{
443  model_ = model;
444  type_ = checkPossible();
445}
446int CbcFathomDynamicProgramming::fathom(double *&betterSolution)
447{
448  int returnCode = 0;
449  int type = checkPossible(maximumSizeAllowed_);
450  assert(type != -1);
451  if (type == -2) {
452    // infeasible (so complete search done)
453    return 1;
454  }
455  if (algorithm_ >= 0) {
456    OsiSolverInterface *solver = model_->solver();
457    const double *lower = solver->getColLower();
458    const double *upper = solver->getColUpper();
459    const double *objective = solver->getObjCoefficients();
460    double direction = solver->getObjSense();
461    const CoinPackedMatrix *matrix = solver->getMatrixByCol();
462    // Column copy
463    const double *element = matrix->getElements();
464    const int *row = matrix->getIndices();
465    const CoinBigIndex *columnStart = matrix->getVectorStarts();
466    const int *columnLength = matrix->getVectorLengths();
467    const double *rowLower = solver->getRowLower();
468    const double *rowUpper = solver->getRowUpper();
469    int numberRows = model_->getNumRows();
470
471    int numberColumns = solver->getNumCols();
472    double offset;
473    solver->getDblParam(OsiObjOffset, offset);
474    double fixedObj = -offset;
475    int i;
476    // may be possible
477    double bestAtTarget = COIN_DBL_MAX;
478    for (i = 0; i < numberColumns; i++) {
479      if (size_ > 10000000 && (i % 100) == 0)
480        COIN_DETAIL_PRINT(printf("column %d\n", i));
481      double lowerValue = lower[i];
482      assert(lowerValue == floor(lowerValue));
483      double cost = direction * objective[i];
484      fixedObj += lowerValue * cost;
485      int gap = static_cast<int>(upper[i] - lowerValue);
486      CoinBigIndex start = columnStart[i];
487      tryColumn(columnLength[i], row + start, element + start, cost, gap);
488      if (cost_[target_] < bestAtTarget) {
489        if (model_->messageHandler()->logLevel() > 1)
490          printf("At column %d new best objective of %g\n", i, cost_[target_]);
491        bestAtTarget = cost_[target_];
492      }
493    }
494    returnCode = 1;
495    int needed = 0;
496    double bestValue = COIN_DBL_MAX;
497    int iBest = -1;
498    if (algorithm_ == 0) {
499      int numberActive = 0;
500      for (i = 0; i < numberRows; i++) {
501        int newRow = lookup_[i];
502        if (newRow >= 0) {
503          if (rowLower[i] == rowUpper[i]) {
504            needed += 1 << numberActive;
505            numberActive++;
506          }
507        }
508      }
509      for (i = 0; i < size_; i++) {
510        if ((i & needed) == needed) {
511          // this one will do
512          if (cost_[i] < bestValue) {
513            bestValue = cost_[i];
514            iBest = i;
515          }
516        }
517      }
518    } else {
519      int *lower = new int[numberActive_];
520      for (i = 0; i < numberRows; i++) {
521        int newRow = lookup_[i];
522        if (newRow >= 0) {
523          int gap = static_cast<int>(rowUpper[i] - CoinMax(0.0, rowLower[i]));
524          lower[newRow] = rhs_[newRow] - gap;
525          int numberBits = numberBits_[newRow];
526          int startBit = startBit_[newRow];
527          if (numberBits == 1 && !gap) {
528            needed |= 1 << startBit;
529          }
530        }
531      }
532      for (i = 0; i < size_; i++) {
533        if ((i & needed) == needed) {
534          // this one may do
535          bool good = true;
536          for (int kk = 0; kk < numberActive_; kk++) {
537            int numberBits = numberBits_[kk];
538            int startBit = startBit_[kk];
539            int size = 1 << numberBits;
540            int start = 1 << startBit;
541            int mask = start * (size - 1);
542            int level = (i & mask) >> startBit;
543            if (level < lower[kk]) {
544              good = false;
545              break;
546            }
547          }
548          if (good && cost_[i] < bestValue) {
549            bestValue = cost_[i];
550            iBest = i;
551          }
552        }
553      }
554      delete[] lower;
555    }
556    if (bestValue < COIN_DBL_MAX) {
557      bestValue += fixedObj;
558      if (model_->messageHandler()->logLevel() > 1)
559        printf("Can get solution of %g\n", bestValue);
560      if (bestValue < model_->getMinimizationObjValue()) {
561        // set up solution
562        betterSolution = new double[numberColumns];
563        memcpy(betterSolution, lower, numberColumns * sizeof(double));
564        while (iBest > 0) {
565          int n = decodeBitPattern(iBest - back_[iBest], indices_, numberRows);
566          // Search for cheapest
567          double bestCost = COIN_DBL_MAX;
568          int iColumn = -1;
569          for (i = 0; i < numberColumns; i++) {
570            if (n == columnLength[i]) {
571              bool good = true;
572              for (CoinBigIndex j = columnStart[i];
573                   j < columnStart[i] + columnLength[i]; j++) {
574                int iRow = row[j];
575                double value = element[j];
576                int iValue = static_cast<int>(value);
577                if (iValue != indices_[iRow]) {
578                  good = false;
579                  break;
580                }
581              }
582              if (good && objective[i] < bestCost && betterSolution[i] < upper[i]) {
583                bestCost = objective[i];
584                iColumn = i;
585              }
586            }
587          }
588          assert(iColumn >= 0);
589          betterSolution[iColumn]++;
590          assert(betterSolution[iColumn] <= upper[iColumn]);
591          iBest = back_[iBest];
592        }
593      }
594      // paranoid check
595      double *rowActivity = new double[numberRows];
596      memset(rowActivity, 0, numberRows * sizeof(double));
597      for (i = 0; i < numberColumns; i++) {
598        CoinBigIndex j;
599        double value = betterSolution[i];
600        if (value) {
601          for (j = columnStart[i];
602               j < columnStart[i] + columnLength[i]; j++) {
603            int iRow = row[j];
604            rowActivity[iRow] += value * element[j];
605          }
606        }
607      }
608      // check was feasible
609      bool feasible = true;
610      for (i = 0; i < numberRows; i++) {
611        if (rowActivity[i] < rowLower[i]) {
612          if (rowActivity[i] < rowLower[i] - 1.0e-8)
613            feasible = false;
614        } else if (rowActivity[i] > rowUpper[i]) {
615          if (rowActivity[i] > rowUpper[i] + 1.0e-8)
616            feasible = false;
617        }
618      }
619      if (feasible) {
620        if (model_->messageHandler()->logLevel() > 0)
621          printf("** good solution of %g by dynamic programming\n", bestValue);
622      }
623      delete[] rowActivity;
624    }
625    gutsOfDelete();
626  }
627  return returnCode;
628}
629/* Tries a column
630   returns true if was used in making any changes.
631*/
632bool CbcFathomDynamicProgramming::tryColumn(int numberElements, const int *rows,
633  const double *coefficients, double cost,
634  int upper)
635{
636  bool touched = false;
637  int n = 0;
638  if (algorithm_ == 0) {
639    for (int j = 0; j < numberElements; j++) {
640      int iRow = rows[j];
641      double value = coefficients[j];
642      int newRow = lookup_[iRow];
643      if (newRow < 0 || value > rhs_[newRow]) {
644        n = 0;
645        break; //can't use
646      } else {
647        indices_[n++] = newRow;
648      }
649    }
650    if (n && upper) {
651      touched = addOneColumn0(n, indices_, cost);
652    }
653  } else {
654    for (int j = 0; j < numberElements; j++) {
655      int iRow = rows[j];
656      double value = coefficients[j];
657      int iValue = static_cast<int>(value);
658      int newRow = lookup_[iRow];
659      if (newRow < 0 || iValue > rhs_[newRow]) {
660        n = 0;
661        break; //can't use
662      } else {
663        coefficients_[n] = iValue;
664        indices_[n++] = newRow;
665        if (upper * iValue > rhs_[newRow]) {
666          upper = rhs_[newRow] / iValue;
667        }
668      }
669    }
670    if (n) {
671      if (algorithm_ == 1) {
672        for (int k = 1; k <= upper; k++) {
673          bool t = addOneColumn1(n, indices_, coefficients_, cost);
674          if (t)
675            touched = true;
676        }
677      } else {
678        CoinSort_2(indices_, indices_ + n, coefficients_);
679        for (int k = 1; k <= upper; k++) {
680          bool t = addOneColumn1A(n, indices_, coefficients_, cost);
681          if (t)
682            touched = true;
683        }
684      }
685    }
686  }
687  return touched;
688}
689/* Adds one column if type 0,
690   returns true if was used in making any changes
691*/
692bool CbcFathomDynamicProgramming::addOneColumn0(int numberElements, const int *rows,
693  double cost)
694{
695  // build up mask
696  int mask = 0;
697  int i;
698  for (i = 0; i < numberElements; i++) {
699    int iRow = rows[i];
700    mask |= 1 << iRow;
701  }
702  bitPattern_ = mask;
703  i = size_ - 1 - mask;
704  bool touched = false;
705  while (i >= 0) {
706    int kMask = i & mask;
707    if (kMask == 0) {
708      double thisCost = cost_[i];
709      if (thisCost != COIN_DBL_MAX) {
710        // possible
711        double newCost = thisCost + cost;
712        int next = i + mask;
713        if (cost_[next] > newCost) {
714          cost_[next] = newCost;
715          back_[next] = i;
716          touched = true;
717        }
718      }
719      i--;
720    } else {
721      // we can skip some
722      int k = (i & ~mask);
723#ifdef CBC_DEBUG
724      for (int j = i - 1; j > k; j--) {
725        int jMask = j & mask;
726        assert(jMask != 0);
727      }
728#endif
729      i = k;
730    }
731  }
732  return touched;
733}
734/* Adds one attempt of one column of type 1,
735   returns true if was used in making any changes.
736   At present the user has to call it once for each possible value
737*/
738bool CbcFathomDynamicProgramming::addOneColumn1(int numberElements, const int *rows,
739  const int *coefficients, double cost)
740{
741  /* build up masks.
742       a) mask for 1 rhs
743       b) mask for addition
744       c) mask so adding will overflow
745       d) individual masks
746    */
747  int mask1 = 0;
748  int maskAdd = 0;
749  int mask2 = 0;
750  int i;
751  int n2 = 0;
752  int mask[40];
753  int adjust[40];
754  assert(numberElements <= 40);
755  for (i = 0; i < numberElements; i++) {
756    int iRow = rows[i];
757    int numberBits = numberBits_[iRow];
758    int startBit = startBit_[iRow];
759    if (numberBits == 1) {
760      mask1 |= 1 << startBit;
761      maskAdd |= 1 << startBit;
762      mask2 |= 1 << startBit;
763    } else {
764      int value = coefficients[i];
765      int size = 1 << numberBits;
766      int start = 1 << startBit;
767      assert(value < size);
768      maskAdd |= start * value;
769      int gap = size - rhs_[iRow] - 1;
770      assert(gap >= 0);
771      int hi2 = rhs_[iRow] - value;
772      if (hi2 < size - 1)
773        hi2++;
774      adjust[n2] = start * hi2;
775      mask2 += start * gap;
776      mask[n2++] = start * (size - 1);
777    }
778  }
779  bitPattern_ = maskAdd;
780  i = size_ - 1 - maskAdd;
781  bool touched = false;
782  while (i >= 0) {
783    int kMask = i & mask1;
784    if (kMask == 0) {
785      bool good = true;
786      for (int kk = n2 - 1; kk >= 0; kk--) {
787        int iMask = mask[kk];
788        int jMask = iMask & mask2;
789        int kkMask = iMask & i;
790        kkMask += jMask;
791        if (kkMask > iMask) {
792          // we can skip some
793          int k = (i & ~iMask);
794          k |= adjust[kk];
795#ifdef CBC_DEBUG
796          for (int j = i - 1; j > k; j--) {
797            int jMask = j & mask1;
798            if (jMask == 0) {
799              bool good = true;
800              for (int kk = n2 - 1; kk >= 0; kk--) {
801                int iMask = mask[kk];
802                int jMask = iMask & mask2;
803                int kkMask = iMask & i;
804                kkMask += jMask;
805                if (kkMask > iMask) {
806                  good = false;
807                  break;
808                }
809              }
810              assert(!good);
811            }
812          }
813#endif
814          i = k;
815          good = false;
816          break;
817        }
818      }
819      if (good) {
820        double thisCost = cost_[i];
821        if (thisCost != COIN_DBL_MAX) {
822          // possible
823          double newCost = thisCost + cost;
824          int next = i + maskAdd;
825          if (cost_[next] > newCost) {
826            cost_[next] = newCost;
827            back_[next] = i;
828            touched = true;
829          }
830        }
831      }
832      i--;
833    } else {
834      // we can skip some
835      // we can skip some
836      int k = (i & ~mask1);
837#ifdef CBC_DEBUG
838      for (int j = i - 1; j > k; j--) {
839        int jMask = j & mask1;
840        assert(jMask != 0);
841      }
842#endif
843      i = k;
844    }
845  }
846  return touched;
847}
848/* Adds one attempt of one column of type 1,
849   returns true if was used in making any changes.
850   At present the user has to call it once for each possible value
851   This version is when there are enough 1 rhs to do faster
852*/
853bool CbcFathomDynamicProgramming::addOneColumn1A(int numberElements, const int *rows,
854  const int *coefficients, double cost)
855{
856  /* build up masks.
857       a) mask for 1 rhs
858       b) mask for addition
859       c) mask so adding will overflow
860       d) mask for non 1 rhs
861    */
862  int maskA = 0;
863  int maskAdd = 0;
864  int maskC = 0;
865  int maskD = 0;
866  int i;
867  for (i = 0; i < numberElements; i++) {
868    int iRow = rows[i];
869    int numberBits = numberBits_[iRow];
870    int startBit = startBit_[iRow];
871    if (numberBits == 1) {
872      maskA |= 1 << startBit;
873      maskAdd |= 1 << startBit;
874    } else {
875      int value = coefficients[i];
876      int size = 1 << numberBits;
877      int start = 1 << startBit;
878      assert(value < size);
879      maskAdd |= start * value;
880      int gap = size - rhs_[iRow] + value - 1;
881      assert(gap > 0 && gap <= size - 1);
882      maskC |= start * gap;
883      maskD |= start * (size - 1);
884    }
885  }
886  bitPattern_ = maskAdd;
887  int maskDiff = maskD - maskC;
888  i = size_ - 1 - maskAdd;
889  bool touched = false;
890  if (!maskD) {
891    // Just ones
892    while (i >= 0) {
893      int kMask = i & maskA;
894      if (kMask == 0) {
895        double thisCost = cost_[i];
896        if (thisCost != COIN_DBL_MAX) {
897          // possible
898          double newCost = thisCost + cost;
899          int next = i + maskAdd;
900          if (cost_[next] > newCost) {
901            cost_[next] = newCost;
902            back_[next] = i;
903            touched = true;
904          }
905        }
906        i--;
907      } else {
908        // we can skip some
909        int k = (i & ~maskA);
910        i = k;
911      }
912    }
913  } else {
914    // More general
915    while (i >= 0) {
916      int kMask = i & maskA;
917      if (kMask == 0) {
918        int added = i & maskD; // just bits belonging to non 1 rhs
919        added += maskC; // will overflow mask if bad
920        added &= (~maskD);
921        if (added == 0) {
922          double thisCost = cost_[i];
923          if (thisCost != COIN_DBL_MAX) {
924            // possible
925            double newCost = thisCost + cost;
926            int next = i + maskAdd;
927            if (cost_[next] > newCost) {
928              cost_[next] = newCost;
929              back_[next] = i;
930              touched = true;
931            }
932          }
933          i--;
934        } else {
935          // we can skip some
936          int k = i & ~maskD; // clear all
937          // Put back enough - but only below where we are
938          int kk = (numberNonOne_ << 1) - 2;
939          assert(rhs_[kk] > 1);
940          int iMask = 0;
941          for (; kk >= 0; kk -= 2) {
942            iMask = 1 << startBit_[kk + 1];
943            if ((added & iMask) != 0) {
944              iMask--;
945              break;
946            }
947          }
948          assert(kk >= 0);
949          iMask &= maskDiff;
950          k |= iMask;
951          assert(k < i);
952          i = k;
953        }
954      } else {
955        // we can skip some
956        int k = (i & ~maskA);
957        i = k;
958      }
959    }
960  }
961  return touched;
962}
963// update model
964void CbcFathomDynamicProgramming::setModel(CbcModel *model)
965{
966  model_ = model;
967  type_ = checkPossible();
968}
969// Gets bit pattern from original column
970int CbcFathomDynamicProgramming::bitPattern(int numberElements, const int *rows,
971  const int *coefficients)
972{
973  int i;
974  int mask = 0;
975  switch (algorithm_) {
976    // just ones
977  case 0:
978    for (i = 0; i < numberElements; i++) {
979      int iRow = rows[i];
980      iRow = lookup_[iRow];
981      if (iRow >= 0)
982        mask |= 1 << iRow;
983    }
984    break;
985    //
986  case 1:
987  case 2:
988    for (i = 0; i < numberElements; i++) {
989      int iRow = rows[i];
990      iRow = lookup_[iRow];
991      if (iRow >= 0) {
992        int startBit = startBit_[iRow];
993        int value = coefficients[i];
994        int start = 1 << startBit;
995        mask |= start * value;
996      }
997    }
998    break;
999  }
1000  return mask;
1001}
1002// Fills in original column (dense) from bit pattern
1003int CbcFathomDynamicProgramming::decodeBitPattern(int bitPattern,
1004  int *values,
1005  int numberRows)
1006{
1007  int i;
1008  int n = 0;
1009  switch (algorithm_) {
1010    // just ones
1011  case 0:
1012    for (i = 0; i < numberRows; i++) {
1013      values[i] = 0;
1014      int iRow = lookup_[i];
1015      if (iRow >= 0) {
1016        if ((bitPattern & (1 << iRow)) != 0) {
1017          values[i] = 1;
1018          n++;
1019        }
1020      }
1021    }
1022    break;
1023    //
1024  case 1:
1025  case 2:
1026    for (i = 0; i < numberRows; i++) {
1027      values[i] = 0;
1028      int iRow = lookup_[i];
1029      if (iRow >= 0) {
1030        int startBit = startBit_[iRow];
1031        int numberBits = numberBits_[iRow];
1032        int iValue = bitPattern >> startBit;
1033        iValue &= ((1 << numberBits) - 1);
1034        if (iValue) {
1035          values[i] = iValue;
1036          n++;
1037        }
1038      }
1039    }
1040    break;
1041  }
1042  return n;
1043}
Note: See TracBrowser for help on using the repository browser.