source: branches/sandbox/Cbc/src/CbcFathomDynamicProgramming.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 9 years ago

Changed formatting using AStyle -A4 -p

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.2 KB
Line 
1/* $Id: CbcFathomDynamicProgramming.cpp 1286 2009-11-09 23:33:07Z EdwinStraver $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcFathomDynamicProgramming.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "CoinPackedMatrix.hpp"
19#include "CoinSort.hpp"
20// Default Constructor
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)
39{
40
41}
42
43// Constructor from model
44CbcFathomDynamicProgramming::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
65CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming ()
66{
67    gutsOfDelete();
68}
69// Does deleteions
70void
71CbcFathomDynamicProgramming::gutsOfDelete()
72{
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;
89}
90// Clone
91CbcFathom *
92CbcFathomDynamicProgramming::clone() const
93{
94    return new CbcFathomDynamicProgramming(*this);
95}
96
97// Copy constructor
98CbcFathomDynamicProgramming::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    }
129}
130// Returns type
131int
132CbcFathomDynamicProgramming::checkPossible(int allowableSize)
133{
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
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;
630}
631/* Tries a column
632   returns true if was used in making any changes.
633*/
634bool
635CbcFathomDynamicProgramming::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;
691}
692/* Adds one column if type 0,
693   returns true if was used in making any changes
694*/
695bool
696CbcFathomDynamicProgramming::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);
727#ifdef CBC_DEBUG
728            for (int j = i - 1; j > k; j--) {
729                int jMask = j & mask;
730                assert (jMask != 0);
731            }
732#endif
733            i = k;
734        }
735    }
736    return touched;
737}
738/* Adds one attempt of one column of type 1,
739   returns true if was used in making any changes.
740   At present the user has to call it once for each possible value
741*/
742bool
743CbcFathomDynamicProgramming::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];
802#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                    }
820#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);
844#ifdef CBC_DEBUG
845            for (int j = i - 1; j > k; j--) {
846                int jMask = j & mask1;
847                assert (jMask != 0);
848            }
849#endif
850            i = k;
851        }
852    }
853    return touched;
854}
855/* Adds one attempt of one column of type 1,
856   returns true if was used in making any changes.
857   At present the user has to call it once for each possible value
858   This version is when there are enough 1 rhs to do faster
859*/
860bool
861CbcFathomDynamicProgramming::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;
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        }
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;
970}
971// update model
972void CbcFathomDynamicProgramming::setModel(CbcModel * model)
973{
974    model_ = model;
975    type_ = checkPossible();
976}
977// Gets bit pattern from original column
978int 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;
1009}
1010// Fills in original column (dense) from bit pattern
1011int 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
Note: See TracBrowser for help on using the repository browser.