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

Last change on this file since 1886 was 1886, checked in by stefan, 5 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.2 KB
Line 
1/*
2  $Id: CbcFathomDynamicProgramming.cpp 1886 2013-04-06 18:13:35Z stefan $
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 adjust[40];
763    assert (numberElements <= 40);
764    for (i = 0; i < numberElements; i++) {
765        int iRow = rows[i];
766        int numberBits = numberBits_[iRow];
767        int startBit = startBit_[iRow];
768        if (numberBits == 1) {
769            mask1 |= 1 << startBit;
770            maskAdd |= 1 << startBit;
771            mask2 |= 1 << startBit;
772        } else {
773            int value = coefficients[i];
774            int size = 1 << numberBits;
775            int start = 1 << startBit;
776            assert (value < size);
777            maskAdd |= start * value;
778            int gap = size - rhs_[iRow] - 1;
779            assert (gap >= 0);
780            int hi2 = rhs_[iRow] - value;
781            if (hi2 < size - 1)
782                hi2++;
783            adjust[n2] = start * hi2;
784            mask2 += start * gap;
785            mask[n2++] = start * (size - 1);
786        }
787    }
788    bitPattern_ = maskAdd;
789    i = size_ - 1 - maskAdd;
790    bool touched = false;
791    while (i >= 0) {
792        int kMask = i & mask1;
793        if (kMask == 0) {
794            bool good = true;
795            for (int kk = n2 - 1; kk >= 0; kk--) {
796                int iMask = mask[kk];
797                int jMask = iMask & mask2;
798                int kkMask = iMask & i;
799                kkMask += jMask;
800                if (kkMask > iMask) {
801                    // we can skip some
802                    int k = (i&~iMask);
803                    k |= adjust[kk];
804#ifdef CBC_DEBUG
805                    for (int j = i - 1; j > k; j--) {
806                        int jMask = j & mask1;
807                        if (jMask == 0) {
808                            bool good = true;
809                            for (int kk = n2 - 1; kk >= 0; kk--) {
810                                int iMask = mask[kk];
811                                int jMask = iMask & mask2;
812                                int kkMask = iMask & i;
813                                kkMask += jMask;
814                                if (kkMask > iMask) {
815                                    good = false;
816                                    break;
817                                }
818                            }
819                            assert (!good);
820                        }
821                    }
822#endif
823                    i = k;
824                    good = false;
825                    break;
826                }
827            }
828            if (good) {
829                double thisCost = cost_[i];
830                if (thisCost != COIN_DBL_MAX) {
831                    // possible
832                    double newCost = thisCost + cost;
833                    int next = i + maskAdd;
834                    if (cost_[next] > newCost) {
835                        cost_[next] = newCost;
836                        back_[next] = i;
837                        touched = true;
838                    }
839                }
840            }
841            i--;
842        } else {
843            // we can skip some
844            // we can skip some
845            int k = (i&~mask1);
846#ifdef CBC_DEBUG
847            for (int j = i - 1; j > k; j--) {
848                int jMask = j & mask1;
849                assert (jMask != 0);
850            }
851#endif
852            i = k;
853        }
854    }
855    return touched;
856}
857/* Adds one attempt of one column of type 1,
858   returns true if was used in making any changes.
859   At present the user has to call it once for each possible value
860   This version is when there are enough 1 rhs to do faster
861*/
862bool
863CbcFathomDynamicProgramming::addOneColumn1A(int numberElements, const int * rows,
864        const int * coefficients, double cost)
865{
866    /* build up masks.
867       a) mask for 1 rhs
868       b) mask for addition
869       c) mask so adding will overflow
870       d) mask for non 1 rhs
871    */
872    int maskA = 0;
873    int maskAdd = 0;
874    int maskC = 0;
875    int maskD = 0;
876    int i;
877    for (i = 0; i < numberElements; i++) {
878        int iRow = rows[i];
879        int numberBits = numberBits_[iRow];
880        int startBit = startBit_[iRow];
881        if (numberBits == 1) {
882            maskA |= 1 << startBit;
883            maskAdd |= 1 << startBit;
884        } else {
885            int value = coefficients[i];
886            int size = 1 << numberBits;
887            int start = 1 << startBit;
888            assert (value < size);
889            maskAdd |= start * value;
890            int gap = size - rhs_[iRow] + value - 1;
891            assert (gap > 0 && gap <= size - 1);
892            maskC |= start * gap;
893            maskD |= start * (size - 1);
894        }
895    }
896    bitPattern_ = maskAdd;
897    int maskDiff = maskD - maskC;
898    i = size_ - 1 - maskAdd;
899    bool touched = false;
900    if (!maskD) {
901        // Just ones
902        while (i >= 0) {
903            int kMask = i & maskA;
904            if (kMask == 0) {
905                double thisCost = cost_[i];
906                if (thisCost != COIN_DBL_MAX) {
907                    // possible
908                    double newCost = thisCost + cost;
909                    int next = i + maskAdd;
910                    if (cost_[next] > newCost) {
911                        cost_[next] = newCost;
912                        back_[next] = i;
913                        touched = true;
914                    }
915                }
916                i--;
917            } else {
918                // we can skip some
919                int k = (i&~maskA);
920                i = k;
921            }
922        }
923    } else {
924        // More general
925        while (i >= 0) {
926            int kMask = i & maskA;
927            if (kMask == 0) {
928                int added = i & maskD; // just bits belonging to non 1 rhs
929                added += maskC; // will overflow mask if bad
930                added &= (~maskD);
931                if (added == 0) {
932                    double thisCost = cost_[i];
933                    if (thisCost != COIN_DBL_MAX) {
934                        // possible
935                        double newCost = thisCost + cost;
936                        int next = i + maskAdd;
937                        if (cost_[next] > newCost) {
938                            cost_[next] = newCost;
939                            back_[next] = i;
940                            touched = true;
941                        }
942                    }
943                    i--;
944                } else {
945                    // we can skip some
946                    int k = i & ~ maskD; // clear all
947                    // Put back enough - but only below where we are
948                    int kk = (numberNonOne_ << 1) - 2;
949                    assert (rhs_[kk] > 1);
950                    int iMask = 0;
951                    for (; kk >= 0; kk -= 2) {
952                        iMask = 1 << startBit_[kk+1];
953                        if ((added&iMask) != 0) {
954                            iMask--;
955                            break;
956                        }
957                    }
958                    assert (kk >= 0);
959                    iMask &= maskDiff;
960                    k |= iMask;
961                    assert (k < i);
962                    i = k;
963                }
964            } else {
965                // we can skip some
966                int k = (i&~maskA);
967                i = k;
968            }
969        }
970    }
971    return touched;
972}
973// update model
974void CbcFathomDynamicProgramming::setModel(CbcModel * model)
975{
976    model_ = model;
977    type_ = checkPossible();
978}
979// Gets bit pattern from original column
980int CbcFathomDynamicProgramming::bitPattern(int numberElements, const int * rows,
981        const int * coefficients)
982{
983    int i;
984    int mask = 0;
985    switch (algorithm_) {
986        // just ones
987    case 0:
988        for (i = 0; i < numberElements; i++) {
989            int iRow = rows[i];
990            iRow = lookup_[iRow];
991            if (iRow >= 0)
992                mask |= 1 << iRow;
993        }
994        break;
995        //
996    case 1:
997    case 2:
998        for (i = 0; i < numberElements; i++) {
999            int iRow = rows[i];
1000            iRow = lookup_[iRow];
1001            if (iRow >= 0) {
1002                int startBit = startBit_[iRow];
1003                int value = coefficients[i];
1004                int start = 1 << startBit;
1005                mask |= start * value;
1006            }
1007        }
1008        break;
1009    }
1010    return mask;
1011}
1012// Fills in original column (dense) from bit pattern
1013int CbcFathomDynamicProgramming::decodeBitPattern(int bitPattern,
1014        int * values,
1015        int numberRows)
1016{
1017    int i;
1018    int n = 0;
1019    switch (algorithm_) {
1020        // just ones
1021    case 0:
1022        for (i = 0; i < numberRows; i++) {
1023            values[i] = 0;
1024            int iRow = lookup_[i];
1025            if (iRow >= 0) {
1026                if ((bitPattern&(1 << iRow)) != 0) {
1027                    values[i] = 1;
1028                    n++;
1029                }
1030            }
1031        }
1032        break;
1033        //
1034    case 1:
1035    case 2:
1036        for (i = 0; i < numberRows; i++) {
1037            values[i] = 0;
1038            int iRow = lookup_[i];
1039            if (iRow >= 0) {
1040                int startBit = startBit_[iRow];
1041                int numberBits = numberBits_[iRow];
1042                int iValue = bitPattern >> startBit;
1043                iValue &= ((1 << numberBits) - 1);
1044                if (iValue) {
1045                    values[i] = iValue;
1046                    n++;
1047                }
1048            }
1049        }
1050        break;
1051    }
1052    return n;
1053}
1054
1055
Note: See TracBrowser for help on using the repository browser.