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

Last change on this file since 702 was 702, checked in by forrest, 12 years ago

changes for gcc 4.3

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