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

Last change on this file since 904 was 904, checked in by ladanyi, 11 years ago

include cstdlib before cmath to get things to compile on AIX with xlC

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