source: branches/sandbox/Cbc/src/CbcLotsize.cpp @ 1308

Last change on this file since 1308 was 1308, checked in by EdwinStraver, 10 years ago

Broke up CbcBranchDynamic? and CbcBranchLotsize?.cpp.
Updated spreadsheets.

File size: 20.8 KB
Line 
1//Edwin 11/17/2009 - carved out of CbcBranchLotsize
2#if defined(_MSC_VER)
3// Turn off compiler warning about long names
4#  pragma warning(disable:4786)
5#endif
6#include <cassert>
7#include <cstdlib>
8#include <cmath>
9#include <cfloat>
10
11#include "OsiSolverInterface.hpp"
12#include "CbcModel.hpp"
13#include "CbcMessage.hpp"
14#include "CbcBranchLotsize.hpp"
15#include "CoinSort.hpp"
16#include "CoinError.hpp"
17#include "CbcLotsize.hpp"
18
19/*
20  CBC_PRINT 1 just does sanity checks - no printing
21            2
22*/
23//#define CBC_PRINT 1
24// First/last variable to print info on
25#if CBC_PRINT
26// preset does all - change to x,x to just do x
27static int firstPrint = 0;
28static int lastPrint = 1000000;
29static CbcModel * saveModel = NULL;
30#endif
31// Just for debug (CBC_PRINT defined in CbcBranchLotsize.cpp)
32void
33#if CBC_PRINT
34CbcLotsize::printLotsize(double value, bool condition, int type) const
35#else
36CbcLotsize::printLotsize(double , bool , int ) const
37#endif
38{
39#if CBC_PRINT
40    if (columnNumber_ >= firstPrint && columnNumber_ <= lastPrint) {
41        int printIt = CBC_PRINT - 1;
42        // Get details
43        OsiSolverInterface * solver = saveModel->solver();
44        double currentLower = solver->getColLower()[columnNumber_];
45        double currentUpper = solver->getColUpper()[columnNumber_];
46        int i;
47        // See if in a valid range (with two tolerances)
48        bool inRange = false;
49        bool inRange2 = false;
50        double integerTolerance =
51            model_->getDblParam(CbcModel::CbcIntegerTolerance);
52        // increase if type 2
53        if (type == 2) {
54            integerTolerance *= 100.0;
55            type = 0;
56            printIt = 2; // always print
57        }
58        // bounds should match some bound
59        int rangeL = -1;
60        int rangeU = -1;
61        if (rangeType_ == 1) {
62            for (i = 0; i < numberRanges_; i++) {
63                if (fabs(currentLower - bound_[i]) < 1.0e-12)
64                    rangeL = i;
65                if (fabs(currentUpper - bound_[i]) < 1.0e-12)
66                    rangeU = i;
67                if (fabs(value - bound_[i]) < integerTolerance)
68                    inRange = true;
69                if (fabs(value - bound_[i]) < 1.0e8)
70                    inRange2 = true;
71            }
72        } else {
73            for (i = 0; i < numberRanges_; i++) {
74                if (fabs(currentLower - bound_[2*i]) < 1.0e-12)
75                    rangeL = i;
76                if (fabs(currentUpper - bound_[2*i+1]) < 1.0e-12)
77                    rangeU = i;
78                if (value > bound_[2*i] - integerTolerance &&
79                        value < bound_[2*i+1] + integerTolerance)
80                    inRange = true;
81                if (value > bound_[2*i] - integerTolerance &&
82                        value < bound_[2*i+1] + integerTolerance)
83                    inRange = true;
84            }
85        }
86        assert (rangeL >= 0 && rangeU >= 0);
87        bool abortIt = false;
88        switch (type) {
89            // returning from findRange (fall through to just check)
90        case 0:
91            if (printIt) {
92                printf("findRange returns %s for column %d and value %g",
93                       condition ? "true" : "false", columnNumber_, value);
94                if (printIt > 1)
95                    printf(" LP bounds %g, %g", currentLower, currentUpper);
96                printf("\n");
97            }
98            // Should match
99        case 1:
100            if (inRange != condition) {
101                printIt = 2;
102                abortIt = true;
103            }
104            break;
105            //
106        case 2:
107            break;
108            //
109        case 3:
110            break;
111            //
112        case 4:
113            break;
114        }
115    }
116#endif
117}
118/** Default Constructor
119
120*/
121CbcLotsize::CbcLotsize ()
122        : CbcObject(),
123        columnNumber_(-1),
124        rangeType_(0),
125        numberRanges_(0),
126        largestGap_(0),
127        bound_(NULL),
128        range_(0)
129{
130}
131
132/** Useful constructor
133
134  Loads actual upper & lower bounds for the specified variable.
135*/
136CbcLotsize::CbcLotsize (CbcModel * model,
137                        int iColumn, int numberPoints,
138                        const double * points, bool range)
139        : CbcObject(model)
140{
141#if CBC_PRINT
142    if (!saveModel)
143        saveModel = model;
144#endif
145    assert (numberPoints > 0);
146    columnNumber_ = iColumn ;
147    // and set id so can be used for branching
148    id_ = iColumn;
149    // sort ranges
150    int * sort = new int[numberPoints];
151    double * weight = new double [numberPoints];
152    int i;
153    if (range) {
154        rangeType_ = 2;
155    } else {
156        rangeType_ = 1;
157    }
158    for (i = 0; i < numberPoints; i++) {
159        sort[i] = i;
160        weight[i] = points[i*rangeType_];
161    }
162    CoinSort_2(weight, weight + numberPoints, sort);
163    numberRanges_ = 1;
164    largestGap_ = 0;
165    if (rangeType_ == 1) {
166        bound_ = new double[numberPoints+1];
167        bound_[0] = weight[0];
168        for (i = 1; i < numberPoints; i++) {
169            if (weight[i] != weight[i-1])
170                bound_[numberRanges_++] = weight[i];
171        }
172        // and for safety
173        bound_[numberRanges_] = bound_[numberRanges_-1];
174        for (i = 1; i < numberRanges_; i++) {
175            largestGap_ = CoinMax(largestGap_, bound_[i] - bound_[i-1]);
176        }
177    } else {
178        bound_ = new double[2*numberPoints+2];
179        bound_[0] = points[sort[0] * 2];
180        bound_[1] = points[sort[0] * 2 + 1];
181        double lo = bound_[0];
182        double hi = bound_[1];
183        assert (hi >= lo);
184        for (i = 1; i < numberPoints; i++) {
185            double thisLo = points[sort[i] * 2];
186            double thisHi = points[sort[i] * 2 + 1];
187            assert (thisHi >= thisLo);
188            if (thisLo > hi) {
189                bound_[2*numberRanges_] = thisLo;
190                bound_[2*numberRanges_+1] = thisHi;
191                numberRanges_++;
192                lo = thisLo;
193                hi = thisHi;
194            } else {
195                //overlap
196                hi = CoinMax(hi, thisHi);
197                bound_[2*numberRanges_-1] = hi;
198            }
199        }
200        // and for safety
201        bound_[2*numberRanges_] = bound_[2*numberRanges_-2];
202        bound_[2*numberRanges_+1] = bound_[2*numberRanges_-1];
203        for (i = 1; i < numberRanges_; i++) {
204            largestGap_ = CoinMax(largestGap_, bound_[2*i] - bound_[2*i-1]);
205        }
206    }
207    delete [] sort;
208    delete [] weight;
209    range_ = 0;
210}
211
212// Copy constructor
213CbcLotsize::CbcLotsize ( const CbcLotsize & rhs)
214        : CbcObject(rhs)
215
216{
217    columnNumber_ = rhs.columnNumber_;
218    rangeType_ = rhs.rangeType_;
219    numberRanges_ = rhs.numberRanges_;
220    range_ = rhs.range_;
221    largestGap_ = rhs.largestGap_;
222    if (numberRanges_) {
223        assert (rangeType_ > 0 && rangeType_ < 3);
224        bound_ = new double [(numberRanges_+1)*rangeType_];
225        memcpy(bound_, rhs.bound_, (numberRanges_ + 1)*rangeType_*sizeof(double));
226    } else {
227        bound_ = NULL;
228    }
229}
230
231// Clone
232CbcObject *
233CbcLotsize::clone() const
234{
235    return new CbcLotsize(*this);
236}
237
238// Assignment operator
239CbcLotsize &
240CbcLotsize::operator=( const CbcLotsize & rhs)
241{
242    if (this != &rhs) {
243        CbcObject::operator=(rhs);
244        columnNumber_ = rhs.columnNumber_;
245        rangeType_ = rhs.rangeType_;
246        numberRanges_ = rhs.numberRanges_;
247        largestGap_ = rhs.largestGap_;
248        delete [] bound_;
249        range_ = rhs.range_;
250        if (numberRanges_) {
251            assert (rangeType_ > 0 && rangeType_ < 3);
252            bound_ = new double [(numberRanges_+1)*rangeType_];
253            memcpy(bound_, rhs.bound_, (numberRanges_ + 1)*rangeType_*sizeof(double));
254        } else {
255            bound_ = NULL;
256        }
257    }
258    return *this;
259}
260
261// Destructor
262CbcLotsize::~CbcLotsize ()
263{
264    delete [] bound_;
265}
266/* Finds range of interest so value is feasible in range range_ or infeasible
267   between hi[range_] and lo[range_+1].  Returns true if feasible.
268*/
269bool
270CbcLotsize::findRange(double value) const
271{
272    assert (range_ >= 0 && range_ < numberRanges_ + 1);
273    double integerTolerance =
274        model_->getDblParam(CbcModel::CbcIntegerTolerance);
275    int iLo;
276    int iHi;
277    double infeasibility = 0.0;
278    if (rangeType_ == 1) {
279        if (value < bound_[range_] - integerTolerance) {
280            iLo = 0;
281            iHi = range_ - 1;
282        } else if (value < bound_[range_] + integerTolerance) {
283#if CBC_PRINT
284            printLotsize(value, true, 0);
285#endif
286            return true;
287        } else if (value < bound_[range_+1] - integerTolerance) {
288#ifdef CBC_PRINT
289            printLotsize(value, false, 0);
290#endif
291            return false;
292        } else {
293            iLo = range_ + 1;
294            iHi = numberRanges_ - 1;
295        }
296        // check lo and hi
297        bool found = false;
298        if (value > bound_[iLo] - integerTolerance && value < bound_[iLo+1] + integerTolerance) {
299            range_ = iLo;
300            found = true;
301        } else if (value > bound_[iHi] - integerTolerance && value < bound_[iHi+1] + integerTolerance) {
302            range_ = iHi;
303            found = true;
304        } else {
305            range_ = (iLo + iHi) >> 1;
306        }
307        //points
308        while (!found) {
309            if (value < bound_[range_]) {
310                if (value >= bound_[range_-1]) {
311                    // found
312                    range_--;
313                    break;
314                } else {
315                    iHi = range_;
316                }
317            } else {
318                if (value < bound_[range_+1]) {
319                    // found
320                    break;
321                } else {
322                    iLo = range_;
323                }
324            }
325            range_ = (iLo + iHi) >> 1;
326        }
327        if (value - bound_[range_] <= bound_[range_+1] - value) {
328            infeasibility = value - bound_[range_];
329        } else {
330            infeasibility = bound_[range_+1] - value;
331            if (infeasibility < integerTolerance)
332                range_++;
333        }
334#ifdef CBC_PRINT
335        printLotsize(value, (infeasibility < integerTolerance), 0);
336#endif
337        return (infeasibility < integerTolerance);
338    } else {
339        // ranges
340        if (value < bound_[2*range_] - integerTolerance) {
341            iLo = 0;
342            iHi = range_ - 1;
343        } else if (value < bound_[2*range_+1] + integerTolerance) {
344#ifdef CBC_PRINT
345            printLotsize(value, true, 0);
346#endif
347            return true;
348        } else if (value < bound_[2*range_+2] - integerTolerance) {
349#ifdef CBC_PRINT
350            printLotsize(value, false, 0);
351#endif
352            return false;
353        } else {
354            iLo = range_ + 1;
355            iHi = numberRanges_ - 1;
356        }
357        // check lo and hi
358        bool found = false;
359        if (value > bound_[2*iLo] - integerTolerance && value < bound_[2*iLo+2] - integerTolerance) {
360            range_ = iLo;
361            found = true;
362        } else if (value >= bound_[2*iHi] - integerTolerance) {
363            range_ = iHi;
364            found = true;
365        } else {
366            range_ = (iLo + iHi) >> 1;
367        }
368        //points
369        while (!found) {
370            if (value < bound_[2*range_]) {
371                if (value >= bound_[2*range_-2]) {
372                    // found
373                    range_--;
374                    break;
375                } else {
376                    iHi = range_;
377                }
378            } else {
379                if (value < bound_[2*range_+2]) {
380                    // found
381                    break;
382                } else {
383                    iLo = range_;
384                }
385            }
386            range_ = (iLo + iHi) >> 1;
387        }
388        if (value >= bound_[2*range_] - integerTolerance && value <= bound_[2*range_+1] + integerTolerance)
389            infeasibility = 0.0;
390        else if (value - bound_[2*range_+1] < bound_[2*range_+2] - value) {
391            infeasibility = value - bound_[2*range_+1];
392        } else {
393            infeasibility = bound_[2*range_+2] - value;
394        }
395#ifdef CBC_PRINT
396        printLotsize(value, (infeasibility < integerTolerance), 0);
397#endif
398        return (infeasibility < integerTolerance);
399    }
400}
401/* Returns floor and ceiling
402 */
403void
404CbcLotsize::floorCeiling(double & floorLotsize, double & ceilingLotsize, double value,
405                         double /*tolerance*/) const
406{
407    bool feasible = findRange(value);
408    if (rangeType_ == 1) {
409        floorLotsize = bound_[range_];
410        ceilingLotsize = bound_[range_+1];
411        // may be able to adjust
412        if (feasible && fabs(value - floorLotsize) > fabs(value - ceilingLotsize)) {
413            floorLotsize = bound_[range_+1];
414            ceilingLotsize = bound_[range_+2];
415        }
416    } else {
417        // ranges
418        assert (value >= bound_[2*range_+1]);
419        floorLotsize = bound_[2*range_+1];
420        ceilingLotsize = bound_[2*range_+2];
421    }
422}
423double
424CbcLotsize::infeasibility(const OsiBranchingInformation * /*info*/,
425                          int &preferredWay) const
426{
427    OsiSolverInterface * solver = model_->solver();
428    const double * solution = model_->testSolution();
429    const double * lower = solver->getColLower();
430    const double * upper = solver->getColUpper();
431    double value = solution[columnNumber_];
432    value = CoinMax(value, lower[columnNumber_]);
433    value = CoinMin(value, upper[columnNumber_]);
434    double integerTolerance =
435        model_->getDblParam(CbcModel::CbcIntegerTolerance);
436    /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
437      solution[columnNumber_],upper[columnNumber_]);*/
438    assert (value >= bound_[0] - integerTolerance
439            && value <= bound_[rangeType_*numberRanges_-1] + integerTolerance);
440    double infeasibility = 0.0;
441    bool feasible = findRange(value);
442    if (!feasible) {
443        if (rangeType_ == 1) {
444            if (value - bound_[range_] < bound_[range_+1] - value) {
445                preferredWay = -1;
446                infeasibility = value - bound_[range_];
447            } else {
448                preferredWay = 1;
449                infeasibility = bound_[range_+1] - value;
450            }
451        } else {
452            // ranges
453            if (value - bound_[2*range_+1] < bound_[2*range_+2] - value) {
454                preferredWay = -1;
455                infeasibility = value - bound_[2*range_+1];
456            } else {
457                preferredWay = 1;
458                infeasibility = bound_[2*range_+2] - value;
459            }
460        }
461    } else {
462        // always satisfied
463        preferredWay = -1;
464    }
465    if (infeasibility < integerTolerance)
466        infeasibility = 0.0;
467    else
468        infeasibility /= largestGap_;
469#ifdef CBC_PRINT
470    printLotsize(value, infeasibility, 1);
471#endif
472    return infeasibility;
473}
474/* Column number if single column object -1 otherwise,
475   so returns >= 0
476   Used by heuristics
477*/
478int
479CbcLotsize::columnNumber() const
480{
481    return columnNumber_;
482}
483// This looks at solution and sets bounds to contain solution
484/** More precisely: it first forces the variable within the existing
485    bounds, and then tightens the bounds to make sure the variable is feasible
486*/
487void
488CbcLotsize::feasibleRegion()
489{
490    OsiSolverInterface * solver = model_->solver();
491    const double * lower = solver->getColLower();
492    const double * upper = solver->getColUpper();
493    const double * solution = model_->testSolution();
494    double value = solution[columnNumber_];
495    value = CoinMax(value, lower[columnNumber_]);
496    value = CoinMin(value, upper[columnNumber_]);
497    findRange(value);
498    double nearest;
499    if (rangeType_ == 1) {
500        nearest = bound_[range_];
501        solver->setColLower(columnNumber_, nearest);
502        solver->setColUpper(columnNumber_, nearest);
503    } else {
504        // ranges
505        solver->setColLower(columnNumber_, bound_[2*range_]);
506        solver->setColUpper(columnNumber_, bound_[2*range_+1]);
507        if (value > bound_[2*range_+1])
508            nearest = bound_[2*range_+1];
509        else if (value < bound_[2*range_])
510            nearest = bound_[2*range_];
511        else
512            nearest = value;
513    }
514#ifdef CBC_PRINT
515    // print details
516    printLotsize(value, true, 2);
517#endif
518    // Scaling may have moved it a bit
519    // Lotsizing variables could be a lot larger
520#ifndef NDEBUG
521    double integerTolerance =
522        model_->getDblParam(CbcModel::CbcIntegerTolerance);
523    assert (fabs(value - nearest) <= (100.0 + 10.0*fabs(nearest))*integerTolerance);
524#endif
525}
526CbcBranchingObject *
527CbcLotsize::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
528{
529    //OsiSolverInterface * solver = model_->solver();
530    const double * solution = model_->testSolution();
531    const double * lower = solver->getColLower();
532    const double * upper = solver->getColUpper();
533    double value = solution[columnNumber_];
534    value = CoinMax(value, lower[columnNumber_]);
535    value = CoinMin(value, upper[columnNumber_]);
536    assert (!findRange(value));
537    return new CbcLotsizeBranchingObject(model_, columnNumber_, way,
538                                         value, this);
539}
540
541
542/* Given valid solution (i.e. satisfied) and reduced costs etc
543   returns a branching object which would give a new feasible
544   point in direction reduced cost says would be cheaper.
545   If no feasible point returns null
546*/
547CbcBranchingObject *
548CbcLotsize::preferredNewFeasible() const
549{
550    OsiSolverInterface * solver = model_->solver();
551
552    assert (findRange(model_->testSolution()[columnNumber_]));
553    double dj = solver->getObjSense() * solver->getReducedCost()[columnNumber_];
554    CbcLotsizeBranchingObject * object = NULL;
555    double lo, up;
556    if (dj >= 0.0) {
557        // can we go down
558        if (range_) {
559            // yes
560            if (rangeType_ == 1) {
561                lo = bound_[range_-1];
562                up = bound_[range_-1];
563            } else {
564                lo = bound_[2*range_-2];
565                up = bound_[2*range_-1];
566            }
567            object = new CbcLotsizeBranchingObject(model_, columnNumber_, -1,
568                                                   lo, up);
569        }
570    } else {
571        // can we go up
572        if (range_ < numberRanges_ - 1) {
573            // yes
574            if (rangeType_ == 1) {
575                lo = bound_[range_+1];
576                up = bound_[range_+1];
577            } else {
578                lo = bound_[2*range_+2];
579                up = bound_[2*range_+3];
580            }
581            object = new CbcLotsizeBranchingObject(model_, columnNumber_, -1,
582                                                   lo, up);
583        }
584    }
585    return object;
586}
587
588/* Given valid solution (i.e. satisfied) and reduced costs etc
589   returns a branching object which would give a new feasible
590   point in direction opposite to one reduced cost says would be cheaper.
591   If no feasible point returns null
592*/
593CbcBranchingObject *
594CbcLotsize::notPreferredNewFeasible() const
595{
596    OsiSolverInterface * solver = model_->solver();
597
598#ifndef NDEBUG
599    double value = model_->testSolution()[columnNumber_];
600    double nearest = floor(value + 0.5);
601    double integerTolerance =
602        model_->getDblParam(CbcModel::CbcIntegerTolerance);
603    // Scaling may have moved it a bit
604    // Lotsizing variables could be a lot larger
605    assert (fabs(value - nearest) <= (10.0 + 10.0*fabs(nearest))*integerTolerance);
606#endif
607    double dj = solver->getObjSense() * solver->getReducedCost()[columnNumber_];
608    CbcLotsizeBranchingObject * object = NULL;
609    double lo, up;
610    if (dj <= 0.0) {
611        // can we go down
612        if (range_) {
613            // yes
614            if (rangeType_ == 1) {
615                lo = bound_[range_-1];
616                up = bound_[range_-1];
617            } else {
618                lo = bound_[2*range_-2];
619                up = bound_[2*range_-1];
620            }
621            object = new CbcLotsizeBranchingObject(model_, columnNumber_, -1,
622                                                   lo, up);
623        }
624    } else {
625        // can we go up
626        if (range_ < numberRanges_ - 1) {
627            // yes
628            if (rangeType_ == 1) {
629                lo = bound_[range_+1];
630                up = bound_[range_+1];
631            } else {
632                lo = bound_[2*range_+2];
633                up = bound_[2*range_+3];
634            }
635            object = new CbcLotsizeBranchingObject(model_, columnNumber_, -1,
636                                                   lo, up);
637        }
638    }
639    return object;
640}
641
642/*
643  Bounds may be tightened, so it may be good to be able to refresh the local
644  copy of the original bounds.
645 */
646void
647CbcLotsize::resetBounds(const OsiSolverInterface * /*solver*/)
648{
649}
Note: See TracBrowser for help on using the repository browser.