source: branches/sandbox/Cbc/src/CbcGeneralDepth.cpp @ 1293

Last change on this file since 1293 was 1293, checked in by EdwinStraver, 10 years ago
File size: 14.3 KB
RevLine 
[1293]1// Edwin 11/10/2009-- carved out of CbcBranchActual
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//#define CBC_DEBUG
11
12#include "CoinTypes.hpp"
13#include "OsiSolverInterface.hpp"
14#include "OsiSolverBranch.hpp"
15#include "CbcModel.hpp"
16#include "CbcMessage.hpp"
17#include "CbcGeneralDepth.hpp"
18#include "CbcBranchActual.hpp"
19#include "CoinSort.hpp"
20#include "CoinError.hpp"
21
22#ifdef COIN_HAS_CLP
23#include "OsiClpSolverInterface.hpp"
24#include "CoinWarmStartBasis.hpp"
25#include "ClpNode.hpp"
26#include "CbcBranchDynamic.hpp"
27// Default Constructor
28CbcGeneralDepth::CbcGeneralDepth ()
29        : CbcGeneral(),
30        maximumDepth_(0),
31        maximumNodes_(0),
32        whichSolution_(-1),
33        numberNodes_(0),
34        nodeInfo_(NULL)
35{
36}
37
38// Useful constructor (which are integer indices)
39CbcGeneralDepth::CbcGeneralDepth (CbcModel * model, int maximumDepth)
40        : CbcGeneral(model),
41        maximumDepth_(maximumDepth),
42        maximumNodes_(0),
43        whichSolution_(-1),
44        numberNodes_(0),
45        nodeInfo_(NULL)
46{
47    assert(maximumDepth_ < 1000000);
48    if (maximumDepth_ > 0)
49        maximumNodes_ = (1 << maximumDepth_) + 1 + maximumDepth_;
50    else if (maximumDepth_ < 0)
51        maximumNodes_ = 1 + 1 - maximumDepth_;
52    else
53        maximumNodes_ = 0;
54#define MAX_NODES 100
55    maximumNodes_ = CoinMin(maximumNodes_, 1 + maximumDepth_ + MAX_NODES);
56    if (maximumNodes_) {
57        nodeInfo_ = new ClpNodeStuff();
58        nodeInfo_->maximumNodes_ = maximumNodes_;
59        ClpNodeStuff * info = nodeInfo_;
60        // for reduced costs and duals
61        info->solverOptions_ |= 7;
62        if (maximumDepth_ > 0) {
63            info->nDepth_ = maximumDepth_;
64        } else {
65            info->nDepth_ = - maximumDepth_;
66            info->solverOptions_ |= 32;
67        }
68        ClpNode ** nodeInfo = new ClpNode * [maximumNodes_];
69        for (int i = 0; i < maximumNodes_; i++)
70            nodeInfo[i] = NULL;
71        info->nodeInfo_ = nodeInfo;
72    } else {
73        nodeInfo_ = NULL;
74    }
75}
76
77// Copy constructor
78CbcGeneralDepth::CbcGeneralDepth ( const CbcGeneralDepth & rhs)
79        : CbcGeneral(rhs)
80{
81    maximumDepth_ = rhs.maximumDepth_;
82    maximumNodes_ = rhs.maximumNodes_;
83    whichSolution_ = -1;
84    numberNodes_ = 0;
85    if (maximumNodes_) {
86        assert (rhs.nodeInfo_);
87        nodeInfo_ = new ClpNodeStuff(*rhs.nodeInfo_);
88        nodeInfo_->maximumNodes_ = maximumNodes_;
89        ClpNodeStuff * info = nodeInfo_;
90        if (maximumDepth_ > 0) {
91            info->nDepth_ = maximumDepth_;
92        } else {
93            info->nDepth_ = - maximumDepth_;
94            info->solverOptions_ |= 32;
95        }
96        if (!info->nodeInfo_) {
97            ClpNode ** nodeInfo = new ClpNode * [maximumNodes_];
98            for (int i = 0; i < maximumNodes_; i++)
99                nodeInfo[i] = NULL;
100            info->nodeInfo_ = nodeInfo;
101        }
102    } else {
103        nodeInfo_ = NULL;
104    }
105}
106
107// Clone
108CbcObject *
109CbcGeneralDepth::clone() const
110{
111    return new CbcGeneralDepth(*this);
112}
113
114// Assignment operator
115CbcGeneralDepth &
116CbcGeneralDepth::operator=( const CbcGeneralDepth & rhs)
117{
118    if (this != &rhs) {
119        CbcGeneral::operator=(rhs);
120        delete nodeInfo_;
121        maximumDepth_ = rhs.maximumDepth_;
122        maximumNodes_ = rhs.maximumNodes_;
123        whichSolution_ = -1;
124        numberNodes_ = 0;
125        if (maximumDepth_) {
126            assert (rhs.nodeInfo_);
127            nodeInfo_ = new ClpNodeStuff(*rhs.nodeInfo_);
128            nodeInfo_->maximumNodes_ = maximumNodes_;
129        } else {
130            nodeInfo_ = NULL;
131        }
132    }
133    return *this;
134}
135
136// Destructor
137CbcGeneralDepth::~CbcGeneralDepth ()
138{
139    delete nodeInfo_;
140}
141// Infeasibility - large is 0.5
142double
143CbcGeneralDepth::infeasibility(const OsiBranchingInformation * /*info*/,
144                               int &/*preferredWay*/) const
145{
146    whichSolution_ = -1;
147    // should use genuine OsiBranchingInformation usefulInfo = model_->usefulInformation();
148    // for now assume only called when correct
149    //if (usefulInfo.depth_>=4&&!model_->parentModel()
150    //     &&(usefulInfo.depth_%2)==0) {
151    if (true) {
152        OsiSolverInterface * solver = model_->solver();
153        OsiClpSolverInterface * clpSolver
154        = dynamic_cast<OsiClpSolverInterface *> (solver);
155        if (clpSolver) {
156            ClpNodeStuff * info = nodeInfo_;
157            info->integerTolerance_ = model_->getIntegerTolerance();
158            info->integerIncrement_ = model_->getCutoffIncrement();
159            info->numberBeforeTrust_ = model_->numberBeforeTrust();
160            info->stateOfSearch_ = model_->stateOfSearch();
161            // Compute "small" change in branch
162            int nBranches = model_->getIntParam(CbcModel::CbcNumberBranches);
163            if (nBranches) {
164                double average = model_->getDblParam(CbcModel::CbcSumChange) / static_cast<double>(nBranches);
165                info->smallChange_ =
166                    CoinMax(average * 1.0e-5, model_->getDblParam(CbcModel::CbcSmallestChange));
167                info->smallChange_ = CoinMax(info->smallChange_, 1.0e-8);
168            } else {
169                info->smallChange_ = 1.0e-8;
170            }
171            int numberIntegers = model_->numberIntegers();
172            double * down = new double[numberIntegers];
173            double * up = new double[numberIntegers];
174            int * priority = new int[numberIntegers];
175            int * numberDown = new int[numberIntegers];
176            int * numberUp = new int[numberIntegers];
177            int * numberDownInfeasible = new int[numberIntegers];
178            int * numberUpInfeasible = new int[numberIntegers];
179            model_->fillPseudoCosts(down, up, priority, numberDown, numberUp,
180                                    numberDownInfeasible, numberUpInfeasible);
181            info->fillPseudoCosts(down, up, priority, numberDown, numberUp,
182                                  numberDownInfeasible,
183                                  numberUpInfeasible, numberIntegers);
184            info->presolveType_ = 1;
185            delete [] down;
186            delete [] up;
187            delete [] numberDown;
188            delete [] numberUp;
189            delete [] numberDownInfeasible;
190            delete [] numberUpInfeasible;
191            bool takeHint;
192            OsiHintStrength strength;
193            solver->getHintParam(OsiDoReducePrint, takeHint, strength);
194            ClpSimplex * simplex = clpSolver->getModelPtr();
195            int saveLevel = simplex->logLevel();
196            if (strength != OsiHintIgnore && takeHint && saveLevel == 1)
197                simplex->setLogLevel(0);
198            clpSolver->setBasis();
199            whichSolution_ = simplex->fathomMany(info);
200            //printf("FAT %d nodes, %d iterations\n",
201            //info->numberNodesExplored_,info->numberIterations_);
202            //printf("CbcBranch %d rows, %d columns\n",clpSolver->getNumRows(),
203            //     clpSolver->getNumCols());
204            model_->incrementExtra(info->numberNodesExplored_,
205                                   info->numberIterations_);
206            // update pseudo costs
207            double smallest = 1.0e50;
208            double largest = -1.0;
209            OsiObject ** objects = model_->objects();
210#ifndef NDEBUG
211            const int * integerVariable = model_->integerVariable();
212#endif
213            for (int i = 0; i < numberIntegers; i++) {
214#ifndef NDEBUG
215                CbcSimpleIntegerDynamicPseudoCost * obj =
216                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
217                assert (obj && obj->columnNumber() == integerVariable[i]);
218#else
219                CbcSimpleIntegerDynamicPseudoCost * obj =
220                    static_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
221#endif
222                if (info->numberUp_[i] > 0) {
223                    if (info->downPseudo_[i] > largest)
224                        largest = info->downPseudo_[i];
225                    if (info->downPseudo_[i] < smallest)
226                        smallest = info->downPseudo_[i];
227                    if (info->upPseudo_[i] > largest)
228                        largest = info->upPseudo_[i];
229                    if (info->upPseudo_[i] < smallest)
230                        smallest = info->upPseudo_[i];
231                    obj->updateAfterMini(info->numberDown_[i],
232                                         info->numberDownInfeasible_[i],
233                                         info->downPseudo_[i],
234                                         info->numberUp_[i],
235                                         info->numberUpInfeasible_[i],
236                                         info->upPseudo_[i]);
237                }
238            }
239            //printf("range of costs %g to %g\n",smallest,largest);
240            simplex->setLogLevel(saveLevel);
241            numberNodes_ = info->nNodes_;
242            int numberDo = numberNodes_;
243            if (whichSolution_ >= 0)
244                numberDo--;
245            if (numberDo > 0) {
246                return 0.5;
247            } else {
248                // no solution
249                return COIN_DBL_MAX; // say infeasible
250            }
251        } else {
252            return -1.0;
253        }
254    } else {
255        return -1.0;
256    }
257}
258
259// This looks at solution and sets bounds to contain solution
260void
261CbcGeneralDepth::feasibleRegion()
262{
263    // Other stuff should have done this
264}
265// Redoes data when sequence numbers change
266void
267CbcGeneralDepth::redoSequenceEtc(CbcModel * /*model*/,
268                                 int /*numberColumns*/,
269                                 const int * /*originalColumns*/)
270{
271}
272
273//#define CHECK_PATH
274#ifdef CHECK_PATH
275extern const double * debuggerSolution_Z;
276extern int numberColumns_Z;
277extern int gotGoodNode_Z;
278#endif
279CbcBranchingObject *
280CbcGeneralDepth::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int /*way*/)
281{
282    int numberDo = numberNodes_;
283    if (whichSolution_ >= 0)
284        numberDo--;
285    assert (numberDo > 0);
286    // create object
287    CbcGeneralBranchingObject * branch = new CbcGeneralBranchingObject(model_);
288    // skip solution
289    branch->numberSubProblems_ = numberDo;
290    // If parentBranch_ back in then will have to be 2*
291    branch->numberSubLeft_ = numberDo;
292    branch->setNumberBranches(numberDo);
293    CbcSubProblem * sub = new CbcSubProblem[numberDo];
294    int iProb = 0;
295    branch->subProblems_ = sub;
296    branch->numberRows_ = model_->solver()->getNumRows();
297    int iNode;
298    //OsiSolverInterface * solver = model_->solver();
299    OsiClpSolverInterface * clpSolver
300    = dynamic_cast<OsiClpSolverInterface *> (solver);
301    assert (clpSolver);
302    ClpSimplex * simplex = clpSolver->getModelPtr();
303    int numberColumns = simplex->numberColumns();
304    double * lowerBefore = CoinCopyOfArray(simplex->getColLower(),
305                                           numberColumns);
306    double * upperBefore = CoinCopyOfArray(simplex->getColUpper(),
307                                           numberColumns);
308    ClpNodeStuff * info = nodeInfo_;
309    double * weight = new double[numberNodes_];
310    int * whichNode = new int [numberNodes_];
311    // Sort
312    for (iNode = 0; iNode < numberNodes_; iNode++) {
313        if (iNode != whichSolution_) {
314            double objectiveValue = info->nodeInfo_[iNode]->objectiveValue();
315            double sumInfeasibilities = info->nodeInfo_[iNode]->sumInfeasibilities();
316            int numberInfeasibilities = info->nodeInfo_[iNode]->numberInfeasibilities();
317            double thisWeight = 0.0;
318#if 1
319            // just closest
320            thisWeight = 1.0e9 * numberInfeasibilities;
321            thisWeight += sumInfeasibilities;
322            thisWeight += 1.0e-7 * objectiveValue;
323            // Try estimate
324            thisWeight = info->nodeInfo_[iNode]->estimatedSolution();
325#else
326            thisWeight = 1.0e-3 * numberInfeasibilities;
327            thisWeight += 1.0e-5 * sumInfeasibilities;
328            thisWeight += objectiveValue;
329#endif
330            whichNode[iProb] = iNode;
331            weight[iProb++] = thisWeight;
332        }
333    }
334    assert (iProb == numberDo);
335    CoinSort_2(weight, weight + numberDo, whichNode);
336    for (iProb = 0; iProb < numberDo; iProb++) {
337        iNode = whichNode[iProb];
338        ClpNode * node = info->nodeInfo_[iNode];
339        // move bounds
340        node->applyNode(simplex, 3);
341        // create subproblem
342        sub[iProb] = CbcSubProblem(clpSolver, lowerBefore, upperBefore,
343                                   node->statusArray(), node->depth());
344        sub[iProb].objectiveValue_ = node->objectiveValue();
345        sub[iProb].sumInfeasibilities_ = node->sumInfeasibilities();
346        sub[iProb].numberInfeasibilities_ = node->numberInfeasibilities();
347#ifdef CHECK_PATH
348        if (simplex->numberColumns() == numberColumns_Z) {
349            bool onOptimal = true;
350            const double * columnLower = simplex->columnLower();
351            const double * columnUpper = simplex->columnUpper();
352            for (int i = 0; i < numberColumns_Z; i++) {
353                if (iNode == gotGoodNode_Z)
354                    printf("good %d %d %g %g\n", iNode, i, columnLower[i], columnUpper[i]);
355                if (columnUpper[i] < debuggerSolution_Z[i] || columnLower[i] > debuggerSolution_Z[i] && simplex->isInteger(i)) {
356                    onOptimal = false;
357                    break;
358                }
359            }
360            if (onOptimal) {
361                printf("adding to node %x as %d - objs\n", this, iProb);
362                for (int j = 0; j <= iProb; j++)
363                    printf("%d %g\n", j, sub[j].objectiveValue_);
364            }
365        }
366#endif
367    }
368    delete [] weight;
369    delete [] whichNode;
370    const double * lower = solver->getColLower();
371    const double * upper = solver->getColUpper();
372    // restore bounds
373    for ( int j = 0; j < numberColumns; j++) {
374        if (lowerBefore[j] != lower[j])
375            solver->setColLower(j, lowerBefore[j]);
376        if (upperBefore[j] != upper[j])
377            solver->setColUpper(j, upperBefore[j]);
378    }
379    delete [] upperBefore;
380    delete [] lowerBefore;
381    return branch;
382}
383#endif
Note: See TracBrowser for help on using the repository browser.