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

Last change on this file since 1750 was 1641, checked in by forrest, 9 years ago

out some printf statements

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