source: trunk/CbcFathomDynamicProgramming.cpp @ 39

Last change on this file since 39 was 39, checked in by forrest, 17 years ago

fixing bug

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