source: trunk/CbcFathomDynamicProgramming.cpp @ 38

Last change on this file since 38 was 38, checked in by forrest, 15 years ago

Slightly faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.8 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 <cmath>
9#include <cfloat>
10
11#include "OsiSolverInterface.hpp"
12#include "CbcModel.hpp"
13#include "CbcMessage.hpp"
14#include "CbcFathomDynamicProgramming.hpp"
15#include "CoinHelperFunctions.hpp"
16#include "CoinPackedMatrix.hpp"
17#include "CoinSort.hpp"
18// Default Constructor
19CbcFathomDynamicProgramming::CbcFathomDynamicProgramming() 
20  :CbcFathom(),
21   size_(0),
22   type_(-1),
23   cost_(NULL),
24   back_(NULL),
25   id_(NULL),
26   lookup_(NULL),
27   numberActive_(0),
28   maximumSizeAllowed_(1000000),
29   startBit_(NULL),
30   numberBits_(NULL),
31   rhs_(NULL),
32   numberNonOne_(0)
33{
34
35}
36
37// Constructor from model
38CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(CbcModel & model)
39  :CbcFathom(model),
40   cost_(NULL),
41   back_(NULL),
42   id_(NULL),
43   lookup_(NULL),
44   numberActive_(0),
45   maximumSizeAllowed_(1000000),
46   startBit_(NULL),
47   numberBits_(NULL),
48   rhs_(NULL),
49   numberNonOne_(0)
50{
51  type_=gutsOfCheckPossible();
52}
53
54// Destructor
55CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming ()
56{
57  gutsOfDelete();
58}
59// Does deleteions
60void 
61CbcFathomDynamicProgramming::gutsOfDelete()
62{
63  delete [] cost_;
64  delete [] back_;
65  delete [] id_;
66  delete [] lookup_;
67  delete [] startBit_;
68  delete [] numberBits_;
69  delete [] rhs_;
70  cost_ = NULL;
71  back_ = NULL;
72  id_ = NULL;
73  lookup_ = NULL;
74  startBit_ = NULL;
75  numberBits_ = NULL;
76  rhs_ = NULL;
77}
78// Clone
79CbcFathom *
80CbcFathomDynamicProgramming::clone() const
81{
82  return new CbcFathomDynamicProgramming(*this);
83}
84
85// Copy constructor
86CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(const CbcFathomDynamicProgramming & rhs)
87:
88  CbcFathom(rhs),
89  size_(rhs.size_),
90  type_(rhs.type_),
91  cost_(NULL),
92  back_(NULL),
93  id_(NULL),
94  lookup_(NULL),
95  numberActive_(rhs.numberActive_),
96  maximumSizeAllowed_(rhs.maximumSizeAllowed_),
97  startBit_(NULL),
98  numberBits_(NULL),
99  rhs_(NULL),
100  numberNonOne_(rhs.numberNonOne_)
101{
102  if (size_) {
103    cost_=CoinCopyOfArray(rhs.cost_,size_);
104    back_=CoinCopyOfArray(rhs.back_,size_);
105    id_=CoinCopyOfArray(rhs.id_,size_);
106    int numberRows=model_->getNumRows();
107    lookup_=CoinCopyOfArray(rhs.lookup_,numberRows);
108    startBit_=CoinCopyOfArray(rhs.startBit_,numberActive_);
109    numberBits_=CoinCopyOfArray(rhs.numberBits_,numberActive_);
110    rhs_=CoinCopyOfArray(rhs.rhs_,numberActive_);
111  }
112}
113// Returns type
114int 
115CbcFathomDynamicProgramming::gutsOfCheckPossible(int allowableSize)
116{
117  assert(model_->solver());
118  OsiSolverInterface * solver = model_->solver();
119  const CoinPackedMatrix * matrix = solver->getMatrixByCol();
120
121  int numberIntegers = model_->numberIntegers();
122  int numberColumns = solver->getNumCols();
123  size_=0;
124  if (numberIntegers!=numberColumns)
125    return -1; // can't do dynamic programming
126
127  const double * lower = solver->getColLower();
128  const double * upper = solver->getColUpper();
129  const double * rowUpper = solver->getRowUpper();
130
131  int numberRows = model_->getNumRows();
132  int i;
133
134  // First check columns to see if possible
135  double * rhs = new double [numberRows];
136  CoinCopyN(rowUpper,numberRows,rhs);
137
138  // Column copy
139  const double * element = matrix->getElements();
140  const int * row = matrix->getIndices();
141  const CoinBigIndex * columnStart = matrix->getVectorStarts();
142  const int * columnLength = matrix->getVectorLengths();
143  bool bad=false;
144  /* It is just possible that we could say okay as
145     variables may get fixed but seems unlikely */
146  for (i=0;i<numberColumns;i++) {
147    int j;
148    double lowerValue = lower[i];
149    assert (lowerValue==floor(lowerValue));
150    for (j=columnStart[i];
151         j<columnStart[i]+columnLength[i];j++) {
152      int iRow=row[j];
153      double value = element[j];
154      if (upper[i]>lowerValue&&(value<=0.0||value!=floor(value)))
155        bad=true;
156      if (lowerValue)
157        rhs[iRow] -= lowerValue*value;
158    }
159  }
160  // check possible (at present do not allow covering)
161  int numberActive=0;
162  for (i=0;i<numberRows;i++) {
163    if (rhs[i]>1.0e5||fabs(rhs[i]-floor(rhs[i]+0.5))>1.0e-7)
164      bad=true;
165    else if (rhs[i]>0.0)
166      numberActive++;
167  }
168  if (bad) {
169    delete [] rhs;
170    return -1;
171  }
172  // check size of array needed
173  double size=1.0;
174  double check = INT_MAX;
175  for (i=0;i<numberRows;i++) {
176    int n= (int) floor(rhs[i]+0.5);
177    if (n) {
178      n++; // allow for 0,1... n
179      if (numberActive!=1) {
180        // power of 2
181        int iBit=0;
182        int k=n;
183        k &= ~1;
184        while (k) {
185          iBit++;
186          k &= ~(1<<iBit);
187        }
188        // See if exact power
189        if (n!=(1<<iBit)) {
190          // round up to next power of 2
191          n= 1<<(iBit+1);
192        }
193        size *= n;
194        if (size>=check)
195          break;
196      } else {
197        size = n; // just one constraint
198      }
199    }
200  }
201  // set size needed
202  if (size>=check)
203    size_=INT_MAX;
204  else
205    size_=(int) size;
206       
207  int n01=0;
208  int nbadcoeff=0;
209  // See if we can tighten bounds
210  for (i=0;i<numberColumns;i++) {
211    int j;
212    double lowerValue = lower[i];
213    double gap = upper[i]-lowerValue;
214    for (j=columnStart[i];
215         j<columnStart[i]+columnLength[i];j++) {
216      int iRow=row[j];
217      double value = element[j];
218      if (value!=1.0)
219        nbadcoeff++;
220      if (gap*value>rhs[iRow]+1.0e-8)
221        gap = rhs[iRow]/value;
222    }
223    gap=lowerValue+floor(gap+1.0e-7);
224    if (gap<upper[i])
225      solver->setColUpper(i,gap);
226    if (gap<=1.0)
227      n01++;
228  }
229  if (allowableSize&&size_<=allowableSize) {
230    numberActive_=numberActive;
231    cost_ = new double [size_];
232    CoinFillN(cost_,size_,COIN_DBL_MAX);
233    // but do nothing is okay
234    cost_[0]=0.0;
235    back_ = new int[size_];
236    CoinFillN(back_,size_,-1);
237    id_ = new int[size_];
238    CoinFillN(id_,size_,-1);
239    startBit_=new int[numberActive_];
240    numberBits_=new int[numberActive_];
241    lookup_ = new int [numberRows];
242    rhs_ = new int [numberActive_];
243    numberActive=0;
244    int kBit=0;
245    for (i=0;i<numberRows;i++) {
246      int n= (int) floor(rhs[i]+0.5);
247      if (n) {
248        lookup_[i]=numberActive;
249        rhs_[numberActive]=n;
250        startBit_[numberActive]=kBit;
251        n++; // allow for 0,1... n
252        int iBit=0;
253        // power of 2
254        int k=n;
255        k &= ~1;
256        while (k) {
257          iBit++;
258          k &= ~(1<<iBit);
259        }
260        // See if exact power
261        if (n!=(1<<iBit)) {
262          // round up to next power of 2
263          iBit++;
264        }
265        if (numberActive!=1) {
266          n= 1<<iBit;
267          size *= n;
268          if (size>=check)
269            break;
270        } else {
271          size = n; // just one constraint
272        }
273        numberBits_[numberActive++]=iBit;
274        kBit += iBit;
275      } else {
276        lookup_[i]=-1;
277      }
278    }
279  }
280  delete [] rhs;
281  if (allowableSize&&size_>allowableSize) {
282    printf("Too large - need %d entries x 16 bytes\n",size_);
283    return -1; // too big
284  }
285  if (n01==numberColumns&&!nbadcoeff)
286    return 0; // easiest
287  else
288    return 1;
289}
290
291// Resets stuff if model changes
292void 
293CbcFathomDynamicProgramming::resetModel(CbcModel * model)
294{
295  model_=model;
296  type_=gutsOfCheckPossible();
297}
298int
299CbcFathomDynamicProgramming::fathom(double * & betterSolution)
300{
301  int returnCode=0;
302  int type=gutsOfCheckPossible(maximumSizeAllowed_);
303  if (type>=0) {
304    bool gotSolution=false;
305    OsiSolverInterface * solver = model_->solver();
306    const double * lower = solver->getColLower();
307    const double * upper = solver->getColUpper();
308    const double * objective = solver->getObjCoefficients();
309    double direction = solver->getObjSense();
310    const CoinPackedMatrix * matrix = solver->getMatrixByCol();
311    // Column copy
312    const double * element = matrix->getElements();
313    const int * row = matrix->getIndices();
314    const CoinBigIndex * columnStart = matrix->getVectorStarts();
315    const int * columnLength = matrix->getVectorLengths();
316    const double * rowLower = solver->getRowLower();
317    const double * rowUpper = solver->getRowUpper();
318    int numberRows = model_->getNumRows();
319   
320    int numberColumns = solver->getNumCols();
321    int * indices = new int [numberActive_];
322    double offset;
323    solver->getDblParam(OsiObjOffset,offset);
324    double fixedObj=-offset;
325    int i;
326    // may be possible
327    if (type==0) {
328      // rhs 1 and coefficients 1
329      for (i=0;i<numberColumns;i++) {
330        int j;
331        double lowerValue = lower[i];
332        assert (lowerValue==floor(lowerValue));
333        double cost = direction * objective[i];
334        fixedObj += lowerValue*cost;
335        if (lowerValue==upper[i])
336          continue;
337        int n=0;
338        for (j=columnStart[i];
339             j<columnStart[i]+columnLength[i];j++) {
340          int iRow=row[j];
341          double value = element[j];
342          int newRow = lookup_[iRow];
343          if (newRow<0||value>rhs_[newRow]) {
344            n=0;
345            break; //can't use
346          } else {
347            indices[n++]=newRow;
348          }
349        }
350        if (n) {
351          addOneColumn0(i,n,indices,cost);
352        }
353      }
354      returnCode=1;
355    } else {
356      // more complex
357      int * coefficients = new int[numberActive_];
358      // If not too many general rhs then we can be more efficient
359      numberNonOne_=0;
360      for (i=0;i<numberActive_;i++) {
361        if (rhs_[i]!=1)
362          numberNonOne_++;
363      }
364      if (numberNonOne_*2<numberActive_) {
365        // put rhs >1 every second
366        int * permute = new int[numberActive_];
367        int * temp = new int[numberActive_];
368        // try different ways
369        int k=0;
370        for (i=0;i<numberRows;i++) {
371          int newRow = lookup_[i];
372          if (newRow>=0&&rhs_[newRow]>1) {
373            permute[newRow]=k;
374            k +=2;
375          }
376        }
377        // adjust so k points to last
378        k -= 2;
379        // and now rest
380        int k1=1;
381        for (i=0;i<numberRows;i++) {
382          int newRow = lookup_[i];
383          if (newRow>=0&&rhs_[newRow]==1) {
384            permute[newRow]=k1;
385            k1++;
386            if (k1<=k)
387              k1++;
388          }
389        }
390        for (i=0;i<numberActive_;i++) {
391          int put = permute[i];
392          temp[put]=rhs_[i];
393        }
394        memcpy(rhs_,temp,numberActive_*sizeof(int));
395        for (i=0;i<numberActive_;i++) {
396          int put = permute[i];
397          temp[put]=numberBits_[i];
398        }
399        memcpy(numberBits_,temp,numberActive_*sizeof(int));
400        k=0;
401        for (i=0;i<numberActive_;i++) {
402          startBit_[i]=k;
403          k+= numberBits_[i];
404        }
405        for (i=0;i<numberRows;i++) {
406          int newRow = lookup_[i];
407          if (newRow>=0)
408            lookup_[i]=permute[newRow];
409        }
410        delete [] permute;
411        delete [] temp;
412        // mark new method
413        type=2;
414      }
415      for (i=0;i<numberColumns;i++) {
416        int j;
417        double lowerValue = lower[i];
418        assert (lowerValue==floor(lowerValue));
419        double cost = direction * objective[i];
420        fixedObj += lowerValue*cost;
421        if (lowerValue==upper[i])
422          continue;
423        int n=0;
424        double gap=upper[i]-lowerValue;
425        for (j=columnStart[i];
426             j<columnStart[i]+columnLength[i];j++) {
427          int iRow=row[j];
428          double value = element[j];
429          int iValue = (int) value;
430          int newRow = lookup_[iRow];
431          if (newRow<0||iValue>rhs_[newRow]) {
432            n=0;
433            break; //can't use
434          } else {
435            coefficients[n]=iValue;
436            indices[n++]=newRow;
437            if (gap*iValue>rhs_[newRow]) {
438              gap = rhs_[newRow]/iValue;
439            }
440          }
441        }
442        int nTotal = (int) gap;
443        if (n) {
444          if (type==1) {
445            for (int k=1;k<=nTotal;k++) {
446              addOneColumn1(i,n,indices,coefficients,cost);
447            }
448          } else {
449            // may be able to take sort out with new method
450            CoinSort_2(indices,indices+n,coefficients);
451            for (int k=1;k<=nTotal;k++) {
452              addOneColumn1A(i,n,indices,coefficients,cost);
453            }
454          }
455        }
456      }
457      delete [] coefficients;
458      returnCode=1;
459    }
460    int needed=0;
461    double bestValue=COIN_DBL_MAX;
462    int iBest=-1;
463    if (type==0) {
464      int numberActive=0;
465      for (i=0;i<numberRows;i++) {
466        int newRow = lookup_[i];
467        if (newRow>=0) {
468          if (rowLower[i]==rowUpper[i]) {
469            needed += 1<<numberActive;
470            numberActive++;
471          }
472        }
473      }
474      for (i=0;i<size_;i++) {
475        if ((i&needed)==needed) {
476          // this one will do
477          if (cost_[i]<bestValue) {
478            bestValue=cost_[i];
479            iBest=i;
480          }
481        }
482      }
483    } else {
484      int * lower = new int[numberActive_];
485      for (i=0;i<numberRows;i++) {
486        int newRow = lookup_[i];
487        if (newRow>=0) {
488          int gap=(int) (rowUpper[i]-CoinMax(0.0,rowLower[i]));
489          lower[newRow]=rhs_[newRow]-gap;
490          int numberBits = numberBits_[newRow];
491          int startBit = startBit_[newRow];
492          if (numberBits==1&&!gap) {
493            needed |= 1<<startBit;
494          }
495        }
496      }
497      for (i=0;i<size_;i++) {
498        if ((i&needed)==needed) {
499        // this one may do
500          bool good = true;
501          for (int kk=0;kk<numberActive_;kk++) {
502            int numberBits = numberBits_[kk];
503            int startBit = startBit_[kk];
504            int size = 1<<numberBits;
505            int start = 1<<startBit;
506            int mask = start*(size-1);
507            int level=(i&mask)>>startBit;
508            if (level<lower[kk]) {
509              good=false;
510              break;
511            }
512          }
513          if (good&&cost_[i]<bestValue) {
514            bestValue=cost_[i];
515            iBest=i;
516          }
517        }
518      }
519      delete [] lower;
520    }
521    if (bestValue<COIN_DBL_MAX) {
522      bestValue += fixedObj;
523      printf("Can get solution of %g\n",bestValue);
524      if (bestValue<model_->getMinimizationObjValue()) {
525        // set up solution
526        betterSolution = new double[numberColumns];
527        memcpy(betterSolution,lower,numberColumns*sizeof(double));
528        while (iBest>0) {
529          int iColumn = id_[iBest];
530          assert (iColumn>=0);
531          betterSolution[iColumn]++;
532          assert (betterSolution[iColumn]<=upper[iColumn]);
533          iBest = back_[iBest];
534        }
535      }
536    }
537    delete [] indices;
538    gutsOfDelete();
539    if (gotSolution) {
540      int i;
541      // paranoid check
542      double * rowActivity = new double [numberRows];
543      memset(rowActivity,0,numberRows*sizeof(double));
544      for (i=0;i<numberColumns;i++) {
545        int j;
546        double value = betterSolution[i];
547        if (value) {
548          for (j=columnStart[i];
549               j<columnStart[i]+columnLength[i];j++) {
550            int iRow=row[j];
551            rowActivity[iRow] += value*element[j];
552          }
553        }
554      }
555      // check was feasible
556      bool feasible=true;
557      for (i=0;i<numberRows;i++) {
558        if(rowActivity[i]<rowLower[i]) {
559          if (rowActivity[i]<rowLower[i]-1.0e-8)
560            feasible = false;
561        } else if(rowActivity[i]>rowUpper[i]) {
562          if (rowActivity[i]>rowUpper[i]+1.0e-8)
563            feasible = false;
564        }
565      }
566      if (feasible) {
567        if (model_->messageHandler()->logLevel()>1)
568          printf("** good solution by dynamic programming\n");
569      }
570      delete [] rowActivity;
571    }
572  }
573  return returnCode;
574}
575/* Adds one column if type 0,
576   returns true if was used in making any changes
577*/
578bool 
579CbcFathomDynamicProgramming::addOneColumn0(int id,int numberElements, const int * rows,
580                     double cost)
581{
582  // build up mask
583  int mask=0;
584  int i;
585  for (i=0;i<numberElements;i++) {
586    int iRow=rows[i];
587    mask |= 1<<iRow;
588  }
589  i=size_-1-mask;
590  bool touched = false;
591  while (i>=0) {
592    int kMask = i&mask;
593    if (kMask==0) {
594      double thisCost = cost_[i];
595      if (thisCost!=COIN_DBL_MAX) {
596        // possible
597        double newCost=thisCost+cost;
598        int next = i + mask;
599        if (cost_[next]>newCost) {
600          cost_[next]=newCost;
601          back_[next]=i;
602          id_[next]=id;
603          touched=true;
604        }
605      }
606      i--;
607    } else {
608      // we can skip some
609      int k=(i&~mask)-1;
610#ifdef CBC_DEBUG
611      for (int j=i-1;j>k;j--) {
612        int jMask = j&mask;
613        assert (jMask!=0);
614      }
615#endif
616      i=k;
617    }
618  }
619  return touched;
620}
621/* Adds one attempt of one column of type 1,
622   returns true if was used in making any changes.
623   At present the user has to call it once for each possible value
624*/
625bool 
626CbcFathomDynamicProgramming::addOneColumn1(int id,int numberElements, const int * rows,
627                                           const int * coefficients, double cost)
628{
629  /* build up masks.
630     a) mask for 1 rhs
631     b) mask for addition
632     c) mask so adding will overflow
633     d) individual masks
634  */
635  int mask1=0;
636  int maskAdd=0;
637  int mask2=0;
638  int i;
639  int n2=0;
640  int mask[40];
641  int nextbit[40];
642  int adjust[40];
643  assert (numberElements<=40);
644  for (i=0;i<numberElements;i++) {
645    int iRow=rows[i];
646    int numberBits = numberBits_[iRow];
647    int startBit = startBit_[iRow];
648    if (numberBits==1) {
649      mask1 |= 1<<startBit;
650      maskAdd |= 1<<startBit;
651      mask2 |= 1<<startBit;
652    } else {
653      int value=coefficients[i];
654      int size = 1<<numberBits;
655      int start = 1<<startBit;
656      assert (value<size);
657      maskAdd |= start*value;
658      int gap = size-rhs_[iRow]-1;
659      assert (gap>=0);
660      int hi2=rhs_[iRow]-value;
661      if (hi2<size-1)
662        hi2++;
663      adjust[n2] = start*hi2;
664      mask2 += start*gap;
665      nextbit[n2]=startBit+numberBits;
666      mask[n2++] = start*(size-1);
667    }
668  }
669  i=size_-1-maskAdd;
670  bool touched = false;
671  while (i>=0) {
672    int kMask = i&mask1;
673    if (kMask==0) {
674      bool good=true;
675      for (int kk=n2-1;kk>=0;kk--) {
676        int iMask = mask[kk];
677        int jMask = iMask&mask2;
678        int kkMask = iMask&i;
679        kkMask += jMask;
680        if (kkMask>iMask) {
681          // we can skip some
682          int k=(i&~iMask);
683          k |= adjust[kk]; 
684#ifdef CBC_DEBUG
685          for (int j=i-1;j>k;j--) {
686            int jMask = j&mask1;
687            if (jMask==0) {
688              bool good=true;
689              for (int kk=n2-1;kk>=0;kk--) {
690                int iMask = mask[kk];
691                int jMask = iMask&mask2;
692                int kkMask = iMask&i;
693                kkMask += jMask;
694                if (kkMask>iMask) {
695                  good=false;
696                  break;
697                }
698              }
699              assert (!good);
700            }
701          }
702#endif
703          i=k;
704          good=false;
705          break;
706        }
707      }
708      if (good) {
709        double thisCost = cost_[i];
710        if (thisCost!=COIN_DBL_MAX) {
711          // possible
712          double newCost=thisCost+cost;
713          int next = i + maskAdd;
714          if (cost_[next]>newCost) {
715            cost_[next]=newCost;
716            back_[next]=i;
717            id_[next]=id;
718            touched=true;
719          }
720        }
721      }
722      i--;
723    } else {
724      // we can skip some
725      // we can skip some
726      int k=(i&~mask1)-1;
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   This version is when there are enough 1 rhs to do faster
742*/
743bool 
744CbcFathomDynamicProgramming::addOneColumn1A(int id,int numberElements, const int * rows,
745                                           const int * coefficients, double cost)
746{
747  /* build up masks.
748     a) mask for 1 rhs
749     b) mask for addition
750     c) mask so adding will overflow
751     d) mask for non 1 rhs
752  */
753  int maskA=0;
754  int maskAdd=0;
755  int maskC=0;
756  int maskD=0;
757  int i;
758  for (i=0;i<numberElements;i++) {
759    int iRow=rows[i];
760    int numberBits = numberBits_[iRow];
761    int startBit = startBit_[iRow];
762    if (numberBits==1) {
763      maskA |= 1<<startBit;
764      maskAdd |= 1<<startBit;
765    } else {
766      int value=coefficients[i];
767      int size = 1<<numberBits;
768      int start = 1<<startBit;
769      assert (value<size);
770      maskAdd |= start*value;
771      int gap = size-rhs_[iRow]-1;
772      assert (gap>=0);
773      maskC |= start*gap;
774      maskD |= start*(size-1);
775    }
776  }
777  int maskDiff = maskD-maskC;
778  i=size_-1-maskAdd;
779  bool touched = false;
780  while (i>=0) {
781    int kMask = i&maskA;
782    if (kMask==0) {
783      int added = i & maskD; // just bits belonging to non 1 rhs
784      added += maskC; // will overflow mask if bad
785      added &= (~maskD);
786      if (added == 0) {
787        double thisCost = cost_[i];
788        if (thisCost!=COIN_DBL_MAX) {
789          // possible
790          double newCost=thisCost+cost;
791          int next = i + maskAdd;
792          if (cost_[next]>newCost) {
793            cost_[next]=newCost;
794            back_[next]=i;
795            id_[next]=id;
796            touched=true;
797          }
798        }
799        i--;
800      } else {
801        // we can skip some
802        int k = i & ~ maskD; // clear all
803        // Put back enough - but only below where we are
804        int kk=(numberNonOne_<<1)-2;
805        assert (rhs_[kk]>1);
806        int iMask=0;
807        for(;kk>=0;kk-=2) {
808          iMask = 1<<startBit_[kk+1];
809          if ((added&iMask)!=0) {
810            iMask--;
811            break;
812          }
813        }
814        assert (kk>=0);
815        iMask &= maskDiff;
816        k |= iMask;
817        assert (k<i);
818        i=k;
819      }
820    } else {
821      // we can skip some
822      // we can skip some
823      int k=(i&~maskA)-1;
824#ifdef CBC_DEBUG
825      for (int j=i-1;j>k;j--) {
826        int jMask = j&mask;
827        assert (jMask!=0);
828      }
829#endif
830      i=k;
831    }
832  }
833  return touched;
834}
835// update model
836void CbcFathomDynamicProgramming::setModel(CbcModel * model)
837{
838  model_ = model;
839  type_=gutsOfCheckPossible();
840}
841
842 
Note: See TracBrowser for help on using the repository browser.