source: trunk/CbcFathomDynamicProgramming.cpp @ 36

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

start of CbcFathom? class

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.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// Default Constructor
18CbcFathomDynamicProgramming::CbcFathomDynamicProgramming() 
19  :CbcFathom(),
20   size_(0),
21   type_(-1),
22   cost_(NULL),
23   back_(NULL),
24   id_(NULL),
25   lookup_(NULL),
26   numberActive_(0),
27   startBit_(NULL),
28   numberBits_(NULL),
29   rhs_(NULL)
30{
31
32}
33
34// Constructor from model
35CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(CbcModel & model)
36  :CbcFathom(model),
37   cost_(NULL),
38   back_(NULL),
39   id_(NULL),
40   lookup_(NULL),
41   numberActive_(0),
42   startBit_(NULL),
43   numberBits_(NULL),
44   rhs_(NULL)
45{
46  type_=gutsOfCheckPossible();
47}
48
49// Destructor
50CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming ()
51{
52  gutsOfDelete();
53}
54// Does deleteions
55void 
56CbcFathomDynamicProgramming::gutsOfDelete()
57{
58  delete [] cost_;
59  delete [] back_;
60  delete [] id_;
61  delete [] lookup_;
62  delete [] startBit_;
63  delete [] numberBits_;
64  delete [] rhs_;
65  cost_ = NULL;
66  back_ = NULL;
67  id_ = NULL;
68  lookup_ = NULL;
69  numberActive_ = 0;
70  startBit_ = NULL;
71  numberBits_ = NULL;
72  rhs_ = NULL;
73}
74// Clone
75CbcFathom *
76CbcFathomDynamicProgramming::clone() const
77{
78  return new CbcFathomDynamicProgramming(*this);
79}
80
81// Copy constructor
82CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(const CbcFathomDynamicProgramming & rhs)
83:
84  CbcFathom(rhs),
85  size_(rhs.size_),
86  type_(rhs.type_),
87  cost_(NULL),
88  back_(NULL),
89  id_(NULL),
90  lookup_(NULL),
91  numberActive_(rhs.numberActive_),
92  startBit_(NULL),
93  numberBits_(NULL),
94  rhs_(NULL)
95{
96  if (size_) {
97    cost_=CoinCopyOfArray(rhs.cost_,size_);
98    back_=CoinCopyOfArray(rhs.back_,size_);
99    id_=CoinCopyOfArray(rhs.id_,size_);
100    int numberRows=model_->getNumRows();
101    lookup_=CoinCopyOfArray(rhs.lookup_,numberRows);
102    startBit_=CoinCopyOfArray(rhs.startBit_,numberActive_);
103    numberBits_=CoinCopyOfArray(rhs.numberBits_,numberActive_);
104    rhs_=CoinCopyOfArray(rhs.rhs_,numberActive_);
105  }
106}
107// Returns type
108int 
109CbcFathomDynamicProgramming::gutsOfCheckPossible(int allowableSize)
110{
111  assert(model_->solver());
112  OsiSolverInterface * solver = model_->solver();
113  const CoinPackedMatrix * matrix = solver->getMatrixByCol();
114
115  int numberIntegers = model_->numberIntegers();
116  int numberColumns = solver->getNumCols();
117  size_=0;
118  if (numberIntegers!=numberColumns)
119    return -1; // can't do dynamic programming
120
121  const double * lower = solver->getColLower();
122  const double * upper = solver->getColUpper();
123  const double * rowUpper = solver->getRowUpper();
124
125  int numberRows = model_->getNumRows();
126  int i;
127
128  // First check columns to see if possible
129  double * rhs = new double [numberRows];
130  CoinCopyN(rowUpper,numberRows,rhs);
131
132  // Column copy
133  const double * element = matrix->getElements();
134  const int * row = matrix->getIndices();
135  const CoinBigIndex * columnStart = matrix->getVectorStarts();
136  const int * columnLength = matrix->getVectorLengths();
137  bool bad=false;
138  /* It is just possible that we could say okay as
139     variables may get fixed but seems unlikely */
140  for (i=0;i<numberColumns;i++) {
141    int j;
142    double lowerValue = lower[i];
143    assert (lowerValue==floor(lowerValue));
144    for (j=columnStart[i];
145         j<columnStart[i]+columnLength[i];j++) {
146      int iRow=row[j];
147      double value = element[j];
148      if (upper[i]>lowerValue&&(value<=0.0||value!=floor(value)))
149        bad=true;
150      if (lowerValue)
151        rhs[iRow] -= lowerValue*value;
152    }
153  }
154  // check possible (at present do not allow covering)
155  int numberActive=0;
156  for (i=0;i<numberRows;i++) {
157    if (rhs[i]>1.0e5||fabs(rhs[i]-floor(rhs[i]+0.5))>1.0e-7)
158      bad=true;
159    else if (rhs[i]>0.0)
160      numberActive++;
161  }
162  if (bad) {
163    delete [] rhs;
164    return -1;
165  }
166  // check size of array needed
167  double size=1.0;
168  double check = INT_MAX;
169  for (i=0;i<numberRows;i++) {
170    int n= (int) floor(rhs[i]+0.5);
171    if (n) {
172      n++; // allow for 0,1... n
173      if (numberActive!=1) {
174        // power of 2
175        int iBit=0;
176        int k=n;
177        k &= ~1;
178        while (k) {
179          iBit++;
180          k &= ~(1<<iBit);
181        }
182        // See if exact power
183        if (n!=(1<<iBit)) {
184          // round up to next power of 2
185          n= 1<<(iBit+1);
186        }
187        size *= n;
188        if (size>=check)
189          break;
190      } else {
191        size = n; // just one constraint
192      }
193    }
194  }
195  // set size needed
196  if (size>=check)
197    size_=INT_MAX;
198  else
199    size_=(int) size;
200       
201  int n01=0;
202  int nbadcoeff=0;
203  // See if we can tighten bounds
204  for (i=0;i<numberColumns;i++) {
205    int j;
206    double lowerValue = lower[i];
207    double gap = upper[i]-lowerValue;
208    for (j=columnStart[i];
209         j<columnStart[i]+columnLength[i];j++) {
210      int iRow=row[j];
211      double value = element[j];
212      if (value!=1.0)
213        nbadcoeff++;
214      if (gap*value>rhs[iRow]+1.0e-8)
215        gap = rhs[iRow]/value;
216    }
217    gap=lowerValue+floor(gap+1.0e-7);
218    if (gap<upper[i])
219      solver->setColUpper(i,gap);
220    if (gap<=1.0)
221      n01++;
222  }
223  if (allowableSize&&size_<=allowableSize) {
224    numberActive_=numberActive;
225    cost_ = new double [size_];
226    CoinFillN(cost_,size_,COIN_DBL_MAX);
227    // but do nothing is okay
228    cost_[0]=0.0;
229    back_ = new int[size_];
230    CoinFillN(back_,size_,-1);
231    id_ = new int[size_];
232    CoinFillN(id_,size_,-1);
233    startBit_=new int[numberActive_];
234    numberBits_=new int[numberActive_];
235    lookup_ = new int [numberRows];
236    rhs_ = new int [numberActive_];
237    numberActive=0;
238    int kBit=0;
239    for (i=0;i<numberRows;i++) {
240      int n= (int) floor(rhs[i]+0.5);
241      if (n) {
242        lookup_[i]=numberActive;
243        rhs_[numberActive]=n;
244        startBit_[numberActive]=kBit;
245        n++; // allow for 0,1... n
246        int iBit=0;
247        // power of 2
248        int k=n;
249        k &= ~1;
250        while (k) {
251          iBit++;
252          k &= ~(1<<iBit);
253        }
254        // See if exact power
255        if (n!=(1<<iBit)) {
256          // round up to next power of 2
257          iBit++;
258        }
259        if (numberActive!=1) {
260          n= 1<<iBit;
261          size *= n;
262          if (size>=check)
263            break;
264        } else {
265          size = n; // just one constraint
266        }
267        numberBits_[numberActive++]=iBit;
268        kBit += iBit;
269      } else {
270        lookup_[i]=-1;
271      }
272    }
273  }
274  delete [] rhs;
275  if (n01==numberColumns&&!nbadcoeff)
276    return 0; // easiest
277  else
278    return 1;
279}
280
281// Resets stuff if model changes
282void 
283CbcFathomDynamicProgramming::resetModel(CbcModel * model)
284{
285  model_=model;
286  type_=gutsOfCheckPossible();
287}
288int
289CbcFathomDynamicProgramming::fathom(double * & betterSolution)
290{
291  int returnCode=0;
292  int type=gutsOfCheckPossible(1000000);
293  if (type>=0) {
294    bool gotSolution=false;
295    OsiSolverInterface * solver = model_->solver();
296    const double * lower = solver->getColLower();
297    const double * upper = solver->getColUpper();
298    const double * objective = solver->getObjCoefficients();
299    double direction = solver->getObjSense();
300    const CoinPackedMatrix * matrix = solver->getMatrixByCol();
301    // Column copy
302    const double * element = matrix->getElements();
303    const int * row = matrix->getIndices();
304    const CoinBigIndex * columnStart = matrix->getVectorStarts();
305    const int * columnLength = matrix->getVectorLengths();
306    const double * rowLower = solver->getRowLower();
307    const double * rowUpper = solver->getRowUpper();
308    int numberRows = model_->getNumRows();
309   
310    int numberColumns = solver->getNumCols();
311    // may be possible
312    // for test - just if all 1
313    if (type==0) {
314      int i;
315
316      int * indices = new int [numberActive_];
317      double offset;
318      solver->getDblParam(OsiObjOffset,offset);
319      double fixedObj=-offset;
320      for (i=0;i<numberColumns;i++) {
321        int j;
322        double lowerValue = lower[i];
323        assert (lowerValue==floor(lowerValue));
324        double cost = direction * objective[i];
325        fixedObj += lowerValue*cost;
326        if (lowerValue==upper[i])
327          continue;
328        int n=0;
329        for (j=columnStart[i];
330             j<columnStart[i]+columnLength[i];j++) {
331          int iRow=row[j];
332          double value = element[j];
333          int newRow = lookup_[iRow];
334          if (newRow<0||value>rhs_[newRow]) {
335            n=0;
336            break; //can't use
337          } else {
338            indices[n++]=newRow;
339          }
340        }
341        if (n) {
342          addOneColumn0(i,n,indices,cost);
343        }
344      }
345      int needed=0;
346      int numberActive=0;
347      for (i=0;i<numberRows;i++) {
348        int newRow = lookup_[i];
349        if (newRow>=0) {
350          if (rowLower[i]==rowUpper[i]) {
351            needed += 1<<numberActive;
352          }
353          numberActive++;
354        }
355      }
356      double bestValue=COIN_DBL_MAX;
357      int iBest=-1;
358      for (i=0;i<size_;i++) {
359        if ((i&needed)==needed) {
360          // this one will do
361          if (cost_[i]<bestValue) {
362            bestValue=cost_[i];
363            iBest=i;
364          }
365        }
366      }
367      returnCode=1;
368      if (bestValue<COIN_DBL_MAX) {
369        bestValue += fixedObj;
370        printf("Can get solution of %g\n",bestValue);
371        if (bestValue<model_->getMinimizationObjValue()) {
372          // set up solution
373          betterSolution = new double[numberColumns];
374          memcpy(betterSolution,lower,numberColumns*sizeof(double));
375          while (iBest>0) {
376            int iColumn = id_[iBest];
377            assert (iColumn>=0);
378            betterSolution[iColumn]++;
379            assert (betterSolution[iColumn]<=upper[iColumn]);
380            iBest = back_[iBest];
381          }
382        } else {
383        }
384      }
385      delete [] indices;
386    }
387    gutsOfDelete();
388    if (gotSolution) {
389      int i;
390      // paranoid check
391      double * rowActivity = new double [numberRows];
392      memset(rowActivity,0,numberRows*sizeof(double));
393      for (i=0;i<numberColumns;i++) {
394        int j;
395        double value = betterSolution[i];
396        if (value) {
397          for (j=columnStart[i];
398               j<columnStart[i]+columnLength[i];j++) {
399            int iRow=row[j];
400            rowActivity[iRow] += value*element[j];
401          }
402        }
403      }
404      // check was feasible
405      bool feasible=true;
406      for (i=0;i<numberRows;i++) {
407        if(rowActivity[i]<rowLower[i]) {
408          if (rowActivity[i]<rowLower[i]-1.0e-8)
409            feasible = false;
410        } else if(rowActivity[i]>rowUpper[i]) {
411          if (rowActivity[i]>rowUpper[i]+1.0e-8)
412            feasible = false;
413        }
414      }
415      if (feasible) {
416        if (model_->messageHandler()->logLevel()>1)
417          printf("** good solution by dynamic programming\n");
418      }
419      delete [] rowActivity;
420    }
421  }
422  return returnCode;
423}
424/* Adds one column if type 0,
425   returns true if was used in making any changes
426*/
427bool 
428CbcFathomDynamicProgramming::addOneColumn0(int id,int numberElements, const int * rows,
429                     double cost)
430{
431  // build up mask
432  int mask=0;
433  int i;
434  for (i=0;i<numberElements;i++) {
435    int iRow=rows[i];
436    mask |= 1<<iRow;
437  }
438  i=0;
439  bool touched = false;
440  while (i<size_) {
441    int kMask = i&mask;
442    if (kMask==0) {
443      double thisCost = cost_[i];
444      if (thisCost!=COIN_DBL_MAX) {
445        // possible
446        double newCost=thisCost+cost;
447        int next = i + mask;
448        if (cost_[next]>newCost) {
449          cost_[next]=newCost;
450          back_[next]=i;
451          id_[next]=id;
452          touched=true;
453        }
454      }
455      i++;
456    } else {
457      // we can skip some
458      int k=i;
459      int iBit=0;
460      k &= ~1;
461      while ((k&kMask)!=0) {
462        iBit++;
463        k &= ~(1<<iBit);
464      }
465      // onto next
466      k += 1<<(iBit+1);
467#ifdef CBC_DEBUG
468      for (int j=i+1;j<k;j++) {
469        int jMask = j&mask;
470        assert (jMask!=0);
471      }
472#endif
473      i=k;
474    }
475  }
476  return touched;
477}
478// update model
479void CbcFathomDynamicProgramming::setModel(CbcModel * model)
480{
481  model_ = model;
482  type_=gutsOfCheckPossible();
483}
484
485 
Note: See TracBrowser for help on using the repository browser.