source: stable/2.4/Cbc/src/CbcFathomDynamicProgramming.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.3 KB
Line 
1/* $Id: CbcFathomDynamicProgramming.cpp 1271 2009-11-05 15:57:25Z forrest $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcFathomDynamicProgramming.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "CoinPackedMatrix.hpp"
19#include "CoinSort.hpp"
20// Default Constructor
21CbcFathomDynamicProgramming::CbcFathomDynamicProgramming() 
22  :CbcFathom(),
23   size_(0),
24   type_(-1),
25   cost_(NULL),
26   back_(NULL),
27   lookup_(NULL),
28   indices_(NULL),
29   numberActive_(0),
30   maximumSizeAllowed_(1000000),
31   startBit_(NULL),
32   numberBits_(NULL),
33   rhs_(NULL),
34   coefficients_(NULL),
35   target_(0),
36   numberNonOne_(0),
37   bitPattern_(0),
38   algorithm_(-1)
39{
40
41}
42
43// Constructor from model
44CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(CbcModel & model)
45  :CbcFathom(model),
46   cost_(NULL),
47   back_(NULL),
48   lookup_(NULL),
49   indices_(NULL),
50   numberActive_(0),
51   maximumSizeAllowed_(1000000),
52   startBit_(NULL),
53   numberBits_(NULL),
54   rhs_(NULL),
55   coefficients_(NULL),
56   target_(0),
57   numberNonOne_(0),
58   bitPattern_(0),
59   algorithm_(-1)
60{
61  type_=checkPossible();
62}
63
64// Destructor
65CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming ()
66{
67  gutsOfDelete();
68}
69// Does deleteions
70void 
71CbcFathomDynamicProgramming::gutsOfDelete()
72{
73  delete [] cost_;
74  delete [] back_;
75  delete [] lookup_;
76  delete [] indices_;
77  delete [] startBit_;
78  delete [] numberBits_;
79  delete [] rhs_;
80  delete [] coefficients_;
81  cost_ = NULL;
82  back_ = NULL;
83  lookup_ = NULL;
84  indices_ =NULL;
85  startBit_ = NULL;
86  numberBits_ = NULL;
87  rhs_ = NULL;
88  coefficients_ = NULL;
89}
90// Clone
91CbcFathom *
92CbcFathomDynamicProgramming::clone() const
93{
94  return new CbcFathomDynamicProgramming(*this);
95}
96
97// Copy constructor
98CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(const CbcFathomDynamicProgramming & rhs)
99:
100  CbcFathom(rhs),
101  size_(rhs.size_),
102  type_(rhs.type_),
103  cost_(NULL),
104  back_(NULL),
105  lookup_(NULL),
106  indices_(NULL),
107  numberActive_(rhs.numberActive_),
108  maximumSizeAllowed_(rhs.maximumSizeAllowed_),
109  startBit_(NULL),
110  numberBits_(NULL),
111  rhs_(NULL),
112  coefficients_(NULL),
113  target_(rhs.target_),
114  numberNonOne_(rhs.numberNonOne_),
115  bitPattern_(rhs.bitPattern_),
116  algorithm_(rhs.algorithm_)
117{
118  if (size_) {
119    cost_=CoinCopyOfArray(rhs.cost_,size_);
120    back_=CoinCopyOfArray(rhs.back_,size_);
121    int numberRows=model_->getNumRows();
122    lookup_=CoinCopyOfArray(rhs.lookup_,numberRows);
123    startBit_=CoinCopyOfArray(rhs.startBit_,numberActive_);
124    indices_=CoinCopyOfArray(rhs.indices_,numberActive_);
125    numberBits_=CoinCopyOfArray(rhs.numberBits_,numberActive_);
126    rhs_=CoinCopyOfArray(rhs.rhs_,numberActive_);
127    coefficients_=CoinCopyOfArray(rhs.coefficients_,numberActive_);
128  }
129}
130// Returns type
131int 
132CbcFathomDynamicProgramming::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    int 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    int 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    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 
442CbcFathomDynamicProgramming::resetModel(CbcModel * model)
443{
444  model_=model;
445  type_=checkPossible();
446}
447int
448CbcFathomDynamicProgramming::fathom(double * & betterSolution)
449{
450  int returnCode=0;
451  int type=checkPossible(maximumSizeAllowed_);
452  assert (type!=-1);
453  if (type==-2) {
454    // infeasible (so complete search done)
455    return 1;
456  }
457  if (algorithm_>=0) {
458    OsiSolverInterface * solver = model_->solver();
459    const double * lower = solver->getColLower();
460    const double * upper = solver->getColUpper();
461    const double * objective = solver->getObjCoefficients();
462    double direction = solver->getObjSense();
463    const CoinPackedMatrix * matrix = solver->getMatrixByCol();
464    // Column copy
465    const double * element = matrix->getElements();
466    const int * row = matrix->getIndices();
467    const CoinBigIndex * columnStart = matrix->getVectorStarts();
468    const int * columnLength = matrix->getVectorLengths();
469    const double * rowLower = solver->getRowLower();
470    const double * rowUpper = solver->getRowUpper();
471    int numberRows = model_->getNumRows();
472   
473    int numberColumns = solver->getNumCols();
474    double offset;
475    solver->getDblParam(OsiObjOffset,offset);
476    double fixedObj=-offset;
477    int i;
478    // may be possible
479    double bestAtTarget = COIN_DBL_MAX;
480    for (i=0;i<numberColumns;i++) {
481      if (size_>10000000&&(i%100)==0)
482        printf("column %d\n",i);
483      double lowerValue = lower[i];
484      assert (lowerValue==floor(lowerValue));
485      double cost = direction * objective[i];
486      fixedObj += lowerValue*cost;
487      int gap = static_cast<int> (upper[i]-lowerValue);
488      CoinBigIndex start = columnStart[i];
489      tryColumn(columnLength[i],row+start,element+start,cost,gap);
490      if (cost_[target_]<bestAtTarget) {
491        if (model_->messageHandler()->logLevel()>1)
492          printf("At column %d new best objective of %g\n",i,cost_[target_]);
493        bestAtTarget = cost_[target_];
494      }
495    }
496    returnCode=1;
497    int needed=0;
498    double bestValue=COIN_DBL_MAX;
499    int iBest=-1;
500    if (algorithm_==0) {
501      int numberActive=0;
502      for (i=0;i<numberRows;i++) {
503        int newRow = lookup_[i];
504        if (newRow>=0) {
505          if (rowLower[i]==rowUpper[i]) {
506            needed += 1<<numberActive;
507            numberActive++;
508          }
509        }
510      }
511      for (i=0;i<size_;i++) {
512        if ((i&needed)==needed) {
513          // this one will do
514          if (cost_[i]<bestValue) {
515            bestValue=cost_[i];
516            iBest=i;
517          }
518        }
519      }
520    } else {
521      int * lower = new int[numberActive_];
522      for (i=0;i<numberRows;i++) {
523        int newRow = lookup_[i];
524        if (newRow>=0) {
525          int gap=static_cast<int> (rowUpper[i]-CoinMax(0.0,rowLower[i]));
526          lower[newRow]=rhs_[newRow]-gap;
527          int numberBits = numberBits_[newRow];
528          int startBit = startBit_[newRow];
529          if (numberBits==1&&!gap) {
530            needed |= 1<<startBit;
531          }
532        }
533      }
534      for (i=0;i<size_;i++) {
535        if ((i&needed)==needed) {
536        // this one may do
537          bool good = true;
538          for (int kk=0;kk<numberActive_;kk++) {
539            int numberBits = numberBits_[kk];
540            int startBit = startBit_[kk];
541            int size = 1<<numberBits;
542            int start = 1<<startBit;
543            int mask = start*(size-1);
544            int level=(i&mask)>>startBit;
545            if (level<lower[kk]) {
546              good=false;
547              break;
548            }
549          }
550          if (good&&cost_[i]<bestValue) {
551            bestValue=cost_[i];
552            iBest=i;
553          }
554        }
555      }
556      delete [] lower;
557    }
558    if (bestValue<COIN_DBL_MAX) {
559      bestValue += fixedObj;
560      if (model_->messageHandler()->logLevel()>1)
561        printf("Can get solution of %g\n",bestValue);
562      if (bestValue<model_->getMinimizationObjValue()) {
563        // set up solution
564        betterSolution = new double[numberColumns];
565        memcpy(betterSolution,lower,numberColumns*sizeof(double));
566        while (iBest>0) {
567          int n=decodeBitPattern(iBest-back_[iBest],indices_,numberRows);
568          // Search for cheapest
569          double bestCost=COIN_DBL_MAX;
570          int iColumn=-1;
571          for (i=0;i<numberColumns;i++) {
572            if (n==columnLength[i]) {
573              bool good=true;
574              for (int j=columnStart[i];
575                   j<columnStart[i]+columnLength[i];j++) {
576                int iRow=row[j];
577                double value = element[j];
578                int iValue = static_cast<int> (value);
579                if(iValue!=indices_[iRow]) {
580                  good=false;
581                  break;
582                }
583              }
584              if (good && objective[i]<bestCost&&betterSolution[i]<upper[i]) {
585                bestCost=objective[i];
586                iColumn=i;
587              }
588            }
589          }
590          assert (iColumn>=0);
591          betterSolution[iColumn]++;
592          assert (betterSolution[iColumn]<=upper[iColumn]);
593          iBest = back_[iBest];
594        }
595      }
596      // paranoid check
597      double * rowActivity = new double [numberRows];
598      memset(rowActivity,0,numberRows*sizeof(double));
599      for (i=0;i<numberColumns;i++) {
600        int j;
601        double value = betterSolution[i];
602        if (value) {
603          for (j=columnStart[i];
604               j<columnStart[i]+columnLength[i];j++) {
605            int iRow=row[j];
606            rowActivity[iRow] += value*element[j];
607          }
608        }
609      }
610      // check was feasible
611      bool feasible=true;
612      for (i=0;i<numberRows;i++) {
613        if(rowActivity[i]<rowLower[i]) {
614          if (rowActivity[i]<rowLower[i]-1.0e-8)
615            feasible = false;
616        } else if(rowActivity[i]>rowUpper[i]) {
617          if (rowActivity[i]>rowUpper[i]+1.0e-8)
618            feasible = false;
619        }
620      }
621      if (feasible) {
622        if (model_->messageHandler()->logLevel()>0)
623          printf("** good solution of %g by dynamic programming\n",bestValue);
624      }
625      delete [] rowActivity;
626    }
627    gutsOfDelete();
628  }
629  return returnCode;
630}
631/* Tries a column
632   returns true if was used in making any changes.
633*/
634bool 
635CbcFathomDynamicProgramming::tryColumn(int numberElements, const int * rows,
636                                          const double * coefficients, double cost,
637                                          int upper)
638{
639  bool touched=false;
640  int n=0;
641  if (algorithm_==0) {
642    for (int j=0;j<numberElements;j++) {
643      int iRow=rows[j];
644      double value = coefficients[j];
645      int newRow = lookup_[iRow];
646      if (newRow<0||value>rhs_[newRow]) {
647        n=0;
648        break; //can't use
649      } else {
650        indices_[n++]=newRow;
651      }
652    }
653    if (n&&upper) {
654      touched=addOneColumn0(n,indices_,cost);
655    }
656  } else {
657    for (int j=0;j<numberElements;j++) {
658      int iRow=rows[j];
659      double value = coefficients[j];
660      int iValue = static_cast<int> (value);
661      int newRow = lookup_[iRow];
662      if (newRow<0||iValue>rhs_[newRow]) {
663        n=0;
664        break; //can't use
665      } else {
666        coefficients_[n]=iValue;
667        indices_[n++]=newRow;
668        if (upper*iValue>rhs_[newRow]) {
669          upper = rhs_[newRow]/iValue;
670        }
671      }
672    }
673    if (n) {
674      if (algorithm_==1) {
675        for (int k=1;k<=upper;k++) {
676          bool t = addOneColumn1(n,indices_,coefficients_,cost);
677          if (t)
678            touched=true;
679        }
680      } else {
681        CoinSort_2(indices_,indices_+n,coefficients_);
682        for (int k=1;k<=upper;k++) {
683          bool t = addOneColumn1A(n,indices_,coefficients_,cost);
684          if (t)
685            touched=true;
686        }
687      }
688    }
689  }
690  return touched;
691}
692/* Adds one column if type 0,
693   returns true if was used in making any changes
694*/
695bool 
696CbcFathomDynamicProgramming::addOneColumn0(int numberElements, const int * rows,
697                     double cost)
698{
699  // build up mask
700  int mask=0;
701  int i;
702  for (i=0;i<numberElements;i++) {
703    int iRow=rows[i];
704    mask |= 1<<iRow;
705  }
706  bitPattern_=mask;
707  i=size_-1-mask;
708  bool touched = false;
709  while (i>=0) {
710    int kMask = i&mask;
711    if (kMask==0) {
712      double thisCost = cost_[i];
713      if (thisCost!=COIN_DBL_MAX) {
714        // possible
715        double newCost=thisCost+cost;
716        int next = i + mask;
717        if (cost_[next]>newCost) {
718          cost_[next]=newCost;
719          back_[next]=i;
720          touched=true;
721        }
722      }
723      i--;
724    } else {
725      // we can skip some
726      int k=(i&~mask);
727#ifdef CBC_DEBUG
728      for (int j=i-1;j>k;j--) {
729        int jMask = j&mask;
730        assert (jMask!=0);
731      }
732#endif
733      i=k;
734    }
735  }
736  return touched;
737}
738/* Adds one attempt of one column of type 1,
739   returns true if was used in making any changes.
740   At present the user has to call it once for each possible value
741*/
742bool 
743CbcFathomDynamicProgramming::addOneColumn1(int numberElements, const int * rows,
744                                           const int * coefficients, double cost)
745{
746  /* build up masks.
747     a) mask for 1 rhs
748     b) mask for addition
749     c) mask so adding will overflow
750     d) individual masks
751  */
752  int mask1=0;
753  int maskAdd=0;
754  int mask2=0;
755  int i;
756  int n2=0;
757  int mask[40];
758  int nextbit[40];
759  int adjust[40];
760  assert (numberElements<=40);
761  for (i=0;i<numberElements;i++) {
762    int iRow=rows[i];
763    int numberBits = numberBits_[iRow];
764    int startBit = startBit_[iRow];
765    if (numberBits==1) {
766      mask1 |= 1<<startBit;
767      maskAdd |= 1<<startBit;
768      mask2 |= 1<<startBit;
769    } else {
770      int value=coefficients[i];
771      int size = 1<<numberBits;
772      int start = 1<<startBit;
773      assert (value<size);
774      maskAdd |= start*value;
775      int gap = size-rhs_[iRow]-1;
776      assert (gap>=0);
777      int hi2=rhs_[iRow]-value;
778      if (hi2<size-1)
779        hi2++;
780      adjust[n2] = start*hi2;
781      mask2 += start*gap;
782      nextbit[n2]=startBit+numberBits;
783      mask[n2++] = start*(size-1);
784    }
785  }
786  bitPattern_=maskAdd;
787  i=size_-1-maskAdd;
788  bool touched = false;
789  while (i>=0) {
790    int kMask = i&mask1;
791    if (kMask==0) {
792      bool good=true;
793      for (int kk=n2-1;kk>=0;kk--) {
794        int iMask = mask[kk];
795        int jMask = iMask&mask2;
796        int kkMask = iMask&i;
797        kkMask += jMask;
798        if (kkMask>iMask) {
799          // we can skip some
800          int k=(i&~iMask);
801          k |= adjust[kk]; 
802#ifdef CBC_DEBUG
803          for (int j=i-1;j>k;j--) {
804            int jMask = j&mask1;
805            if (jMask==0) {
806              bool good=true;
807              for (int kk=n2-1;kk>=0;kk--) {
808                int iMask = mask[kk];
809                int jMask = iMask&mask2;
810                int kkMask = iMask&i;
811                kkMask += jMask;
812                if (kkMask>iMask) {
813                  good=false;
814                  break;
815                }
816              }
817              assert (!good);
818            }
819          }
820#endif
821          i=k;
822          good=false;
823          break;
824        }
825      }
826      if (good) {
827        double thisCost = cost_[i];
828        if (thisCost!=COIN_DBL_MAX) {
829          // possible
830          double newCost=thisCost+cost;
831          int next = i + maskAdd;
832          if (cost_[next]>newCost) {
833            cost_[next]=newCost;
834            back_[next]=i;
835            touched=true;
836          }
837        }
838      }
839      i--;
840    } else {
841      // we can skip some
842      // we can skip some
843      int k=(i&~mask1);
844#ifdef CBC_DEBUG
845      for (int j=i-1;j>k;j--) {
846        int jMask = j&mask1;
847        assert (jMask!=0);
848      }
849#endif
850      i=k;
851    }
852  }
853  return touched;
854}
855/* Adds one attempt of one column of type 1,
856   returns true if was used in making any changes.
857   At present the user has to call it once for each possible value
858   This version is when there are enough 1 rhs to do faster
859*/
860bool 
861CbcFathomDynamicProgramming::addOneColumn1A(int numberElements, const int * rows,
862                                           const int * coefficients, double cost)
863{
864  /* build up masks.
865     a) mask for 1 rhs
866     b) mask for addition
867     c) mask so adding will overflow
868     d) mask for non 1 rhs
869  */
870  int maskA=0;
871  int maskAdd=0;
872  int maskC=0;
873  int maskD=0;
874  int i;
875  for (i=0;i<numberElements;i++) {
876    int iRow=rows[i];
877    int numberBits = numberBits_[iRow];
878    int startBit = startBit_[iRow];
879    if (numberBits==1) {
880      maskA |= 1<<startBit;
881      maskAdd |= 1<<startBit;
882    } else {
883      int value=coefficients[i];
884      int size = 1<<numberBits;
885      int start = 1<<startBit;
886      assert (value<size);
887      maskAdd |= start*value;
888      int gap = size-rhs_[iRow]+value-1;
889      assert (gap>0&&gap<=size-1);
890      maskC |= start*gap;
891      maskD |= start*(size-1);
892    }
893  }
894  bitPattern_=maskAdd;
895  int maskDiff = maskD-maskC;
896  i=size_-1-maskAdd;
897  bool touched = false;
898  if (!maskD) {
899    // Just ones
900    while (i>=0) {
901      int kMask = i&maskA;
902      if (kMask==0) {
903        double thisCost = cost_[i];
904        if (thisCost!=COIN_DBL_MAX) {
905          // possible
906          double newCost=thisCost+cost;
907          int next = i + maskAdd;
908          if (cost_[next]>newCost) {
909            cost_[next]=newCost;
910            back_[next]=i;
911            touched=true;
912          }
913        }
914        i--;
915      } else {
916        // we can skip some
917        int k=(i&~maskA);
918        i=k;
919      }
920    }
921  } else {
922    // More general
923    while (i>=0) {
924      int kMask = i&maskA;
925      if (kMask==0) {
926        int added = i & maskD; // just bits belonging to non 1 rhs
927        added += maskC; // will overflow mask if bad
928        added &= (~maskD);
929        if (added == 0) {
930          double thisCost = cost_[i];
931          if (thisCost!=COIN_DBL_MAX) {
932            // possible
933            double newCost=thisCost+cost;
934            int next = i + maskAdd;
935            if (cost_[next]>newCost) {
936              cost_[next]=newCost;
937              back_[next]=i;
938              touched=true;
939            }
940          }
941          i--;
942        } else {
943          // we can skip some
944          int k = i & ~ maskD; // clear all
945          // Put back enough - but only below where we are
946          int kk=(numberNonOne_<<1)-2;
947          assert (rhs_[kk]>1);
948          int iMask=0;
949          for(;kk>=0;kk-=2) {
950            iMask = 1<<startBit_[kk+1];
951            if ((added&iMask)!=0) {
952              iMask--;
953              break;
954            }
955          }
956          assert (kk>=0);
957          iMask &= maskDiff;
958          k |= iMask;
959          assert (k<i);
960          i=k;
961        }
962      } else {
963        // we can skip some
964        int k=(i&~maskA);
965        i=k;
966      }
967    }
968  }
969  return touched;
970}
971// update model
972void CbcFathomDynamicProgramming::setModel(CbcModel * model)
973{
974  model_ = model;
975  type_=checkPossible();
976}
977// Gets bit pattern from original column
978int CbcFathomDynamicProgramming::bitPattern(int numberElements, const int * rows,
979                                            const int * coefficients)
980{
981  int i;
982  int mask=0;
983  switch (algorithm_) {
984    // just ones
985  case 0:
986    for (i=0;i<numberElements;i++) {
987      int iRow=rows[i];
988      iRow = lookup_[iRow];
989      if (iRow>=0)
990        mask |= 1<<iRow;
991    }
992    break;
993    //
994  case 1:
995  case 2:
996    for (i=0;i<numberElements;i++) {
997      int iRow=rows[i];
998      iRow = lookup_[iRow];
999      if (iRow>=0) {
1000        int startBit = startBit_[iRow];
1001        int value=coefficients[i];
1002        int start = 1<<startBit;
1003        mask |= start*value;
1004      }
1005    }
1006    break;
1007  }
1008  return mask;
1009}
1010// Fills in original column (dense) from bit pattern
1011int CbcFathomDynamicProgramming::decodeBitPattern(int bitPattern,
1012                                                   int * values,
1013                                                   int numberRows)
1014{
1015  int i;
1016  int n=0;
1017  switch (algorithm_) {
1018    // just ones
1019  case 0:
1020    for (i=0;i<numberRows;i++) {
1021      values[i]=0;
1022      int iRow = lookup_[i];
1023      if (iRow>=0) {
1024        if ((bitPattern&(1<<iRow))!=0) {
1025          values[i]=1;
1026          n++;
1027        }
1028      }
1029    }
1030    break;
1031    //
1032  case 1:
1033  case 2:
1034    for (i=0;i<numberRows;i++) {
1035      values[i]=0;
1036      int iRow = lookup_[i];
1037      if (iRow>=0) {
1038        int startBit = startBit_[iRow];
1039        int numberBits = numberBits_[iRow];
1040        int iValue = bitPattern>>startBit;
1041        iValue &= ((1<<numberBits)-1);
1042        if (iValue) {
1043          values[i]=iValue;
1044          n++;
1045        }
1046      }
1047    }
1048    break;
1049  }
1050  return n;
1051}
1052
1053 
Note: See TracBrowser for help on using the repository browser.