Ignore:
Timestamp:
Nov 9, 2009 6:33:07 PM (10 years ago)
Author:
EdwinStraver
Message:

Changed formatting using AStyle -A4 -p

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/Cbc/src/CbcFathomDynamicProgramming.cpp

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