source: trunk/Alps/examples/Abc/AbcTreeNode.cpp @ 277

Last change on this file since 277 was 277, checked in by andreasw, 13 years ago

first working version with autotools

File size: 35.9 KB
Line 
1/*===========================================================================*
2 * This file is part of the Abstract Library for Parallel Search (ALPS).     *
3 *                                                                           *
4 * ALPS is distributed under the Common Public License as part of the        *
5 * COIN-OR repository (http://www.coin-or.org).                              *
6 *                                                                           *
7 * Authors: Yan Xu, SAS Institute Inc.                                       *
8 *          Ted Ralphs, Lehigh University                                    *
9 *          Laszlo Ladanyi, IBM T.J. Watson Research Center                  *
10 *          Matthew Saltzman, Clemson University                             *
11 *                                                                           *
12 *                                                                           *
13 * Copyright (C) 2001-2004, International Business Machines                  *
14 * Corporation, Lehigh University, Yan Xu, Ted Ralphs, Matthew Salzman and   *
15 * others. All Rights Reserved.                                              *
16 *===========================================================================*/
17
18#include <cassert>
19#include <iostream>
20#include <utility>
21#include <cmath>
22
23#include "CoinUtility.hpp"
24#include "OsiSolverInterface.hpp"
25#include "OsiClpSolverInterface.hpp"
26#include "OsiRowCut.hpp"
27#include "OsiColCut.hpp"
28#include "OsiRowCutDebugger.hpp"
29#include "OsiCuts.hpp"
30
31#include "AlpsKnowledge.h"
32#include "AlpsEnumProcessT.h"
33
34#include "AbcBranchActual.h"
35#include "AbcMessage.h"
36#include "AbcParams.h"
37#include "AbcTreeNode.h"
38#include "AbcSolution.h"
39
40//#############################################################################
41
42/// Infeasibility - large is 0.5
43static double checkInfeasibility(AbcModel* model, const int columnNumber, 
44                                 int & preferredWay, double & otherWay)
45{
46    OsiSolverInterface * solver = model->solver();
47    const double * solution = model->currentSolution();
48    const double * lower = solver->getColLower();
49    const double * upper = solver->getColUpper();
50    double value = solution[columnNumber];
51    value = max(value, lower[columnNumber]);
52    value = min(value, upper[columnNumber]);
53    /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
54    solution[columnNumber_],upper[columnNumber_]);*/
55    double nearest = floor(value + 0.5);
56    double integerTolerance = 
57        model->getDblParam(AbcModel::AbcIntegerTolerance);
58    if (nearest > value) 
59        preferredWay = 1;
60    else
61        preferredWay = -1;
62    otherWay = 1.0 - fabs(value - nearest);
63    if (fabs(value - nearest) <= integerTolerance) 
64        return 0.0;
65    else
66        return fabs(value - nearest);
67}
68
69//#############################################################################
70
71AlpsTreeNode*
72AbcTreeNode::createNewTreeNode(AlpsNodeDesc *&desc) const
73{
74    // Create a new tree node
75    AbcNodeDesc* d = dynamic_cast<AbcNodeDesc*>(desc);
76    AbcTreeNode* node = new AbcTreeNode(d);
77    desc = 0;
78    return(node);
79}
80
81//#############################################################################
82
83int
84AbcTreeNode::process(bool isRoot, bool rampUp)
85{
86    bool betterSolution = false;
87    double bestValue = getKnowledgeBroker()->getIncumbentValue();
88    double parentObjValue = getObjValue();
89    double primalTolerance = 1.0e-7;
90    bool cutDuringRampup;
91    AlpsProcessType myType = getKnowledgeBroker()->getProcType();
92
93    AbcModel *model=dynamic_cast<AbcModel *>(getKnowledgeBroker()->getModel());
94
95    cutDuringRampup = true;
96
97    if (rampUp) {
98        if (myType == AlpsProcessTypeHub || myType == AlpsProcessTypeMaster){
99            cutDuringRampup=model->AbcPar()->entry(AbcParams::cutDuringRampup);
100        }
101    }
102   
103
104    // Check if can quit early. Becare the root.
105    // REMOVE
106    //if (parentObjValue < 1.0e99) bestValue = 3983;
107   
108    if (parentObjValue - primalTolerance > bestValue) {
109        setStatus(AlpsNodeStatusFathomed);
110        return false;
111    }
112
113    int i = -1;
114    AbcNodeDesc *desc = dynamic_cast<AbcNodeDesc*>(desc_);
115    AbcModel *m = dynamic_cast<AbcModel*>(desc->getModel());
116
117    // Update cutoff if recv a better solution from other process
118    if (bestValue < m->getCutoff()) {
119        double cutoff = bestValue - 
120            m->getDblParam(AbcModel::AbcCutoffIncrement);
121        m->setCutoff(cutoff);
122    }
123
124    //-------------------------------------------------------------------------
125    // Report search status
126    int nodeInterval = model->AbcPar()->entry(AbcParams::statusInterval);
127
128    assert(nodeInterval > 0);
129
130    if(m->getNodeCount() % nodeInterval == 0){
131        const int nodeLeft = getKnowledgeBroker()->getNumNodesLeft();
132        m->messageHandler()->message(ABC_STATUS, m->messages())
133            << getKnowledgeBroker()->getProcRank() << (m->getNodeCount()) 
134            << nodeLeft <<  m->getCutoff() << parentObjValue << CoinMessageEol;
135    }
136     
137    //-------------------------------------------------------------------------
138    // Load and solve the lp relaxation.
139    // (1) LP infeasible
140    //     a. set status to be fathom
141    // (2) LP feasible
142    //     a. MILP feasible. Check whether need update incumbent.
143    //     b. LP feasible but not MIP feasible. Check whether can be
144    //        fathomed, if not, choose a branch variable
145    //-------------------------------------------------------------------------
146    const int numCols = m->solver()->getNumCols();
147    const double *lbs = desc->lowerBounds();
148    const double *ubs = desc->upperBounds();
149
150    for(i = 0; i < numCols; ++i) {
151        m->solver()->setColBounds(i, lbs[i], ubs[i]);
152    }
153     
154    //-------------------------------------------------------------
155    OsiCuts cuts;
156    int maximumCutPassesAtRoot = m->getMaximumCutPassesAtRoot();
157    int numberOldActiveCuts = 0;
158    int numberNewCuts = 0;
159    int maximumWhich = 1000;
160    int *whichGenerator = new int[maximumWhich];
161    int heurFound = -1; 
162 
163    bool feasible = m->solveWithCuts(cuts, maximumCutPassesAtRoot,
164                                     NULL, numberOldActiveCuts, 
165                                     numberNewCuts, maximumWhich, 
166                                     whichGenerator, 
167                                     cutDuringRampup,
168                                     heurFound);
169   
170    m->setCurrentNumberCuts(m->currentNumberCuts() + numberNewCuts);
171   
172    // if ( (heurFound >= 0) && (m->getObjValue() < bestValue) ) {
173    if (heurFound >= 0) {
174        betterSolution = true;
175        AbcSolution* sol = new AbcSolution(numCols, 
176                                           m->bestSolution(),
177                                           m->getObjValue());
178        getKnowledgeBroker()->addKnowledge(ALPS_SOLUTION, sol, 
179                                           m->getObjValue());
180    }
181   
182
183    //-------------------------------------------------------------
184     
185    if (!feasible) {
186        setStatus(AlpsNodeStatusFathomed);
187        setObjValue(-ALPS_OBJ_MAX);       // Remove it as soon as possilbe
188    }
189    else {
190        double val = (m->getCurrentObjValue()) * (m->getObjSense());
191        double xS = desc->getBranchedOnValue();
192        int bDir = desc->getBranchedDir();
193        int bInd = desc->getBranchedOn();
194       
195        /* Update pseudo: not for the root */
196        if ((m->numberStrong() == 0) && (parent_ != 0) && (bInd >= 0)) {
197            /* FIXME: not sure which bInd < 0 */
198            // std::cout << "bInd=" << bInd << "; bVal=" << xS << std::endl;
199            (m->getPseudoList()[bInd])->update(bDir, quality_, val, xS);
200        }
201       
202        int numberInfeasibleIntegers;
203        if(m->feasibleSolution(numberInfeasibleIntegers)) { // MIP feasible
204            if (val < bestValue) {
205                betterSolution = true;
206                m->setBestSolution(ABC_BRANCHSOL,
207                                   val, m->getColSolution());
208                AbcSolution* ksol = new AbcSolution(numCols, 
209                                                    m->getColSolution(), 
210                                                    val);
211                getKnowledgeBroker()->addKnowledge(ALPS_SOLUTION, ksol, val); 
212            }
213            setStatus(AlpsNodeStatusFathomed);
214        }
215        else {
216            if (val < bestValue) {
217                bool strongFound = false;
218                int action = -1;
219                while (action == -1) { 
220                    if(getKnowledgeBroker()->getProcRank() == -1) {
221                        std::cout << "*** I AM RANK ONE: before choose:action = " << action
222                                  << std::endl;
223                    }
224                    action = chooseBranch(m, strongFound);
225                    if ( (strongFound) && 
226                         (m->getObjValue() < bestValue) ) {
227                        betterSolution = true;
228                        AbcSolution *sol = new AbcSolution(numCols, 
229                                                           m->bestSolution(),
230                                                           m->getObjValue());
231                        getKnowledgeBroker()->addKnowledge(ALPS_SOLUTION, sol, 
232                                                           m->getObjValue());
233                       
234                    }   
235                    if (action == -1) { 
236                        feasible = m->resolve();
237                        //resolved = true ;
238#if defined(ABC_DEBUG_MORE)
239                        printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
240                               m->solver()->getObjValue(),
241                               m->solver()->getNumRows());
242#endif
243                        if (!feasible) action = -2; 
244                    }
245                    if (action == -2) { 
246                        feasible = false; 
247                    } 
248                    if(getKnowledgeBroker()->getProcRank() == -1) {
249                        std::cout << "*** I AM RANK ONE: action = " << action
250                                  << std::endl;
251                    }
252                }
253                assert(action != -1);
254
255
256                if (action >= 0) {
257                    const double * newLbs = m->getColLower();
258                    const double * newUbs = m->getColUpper();
259                    desc->setLowerBounds(newLbs, numCols);
260                    desc->setUpperBounds(newUbs, numCols);
261#if defined(ABC_DEBUG_MORE)
262                    std::cout << "SetPregnant: branchedOn = " << branchedOn_
263                              << "; index = " << index_ << std::endl;
264#endif             
265                    setStatus(AlpsNodeStatusPregnant);
266                }
267                else if (action == -2) {
268                    setStatus(AlpsNodeStatusFathomed);
269                }
270                else {
271                    throw CoinError("No branch object found", "process", 
272                                    "AbcTreeNode");
273                }
274               
275            }
276            else {
277                setStatus(AlpsNodeStatusFathomed);
278            }
279        }
280        quality_ = val;
281    }
282
283    m->takeOffCuts();
284
285    delete [] whichGenerator;
286   
287    return betterSolution;
288}
289
290//#############################################################################
291
292std::vector< CoinTriple<AlpsNodeDesc*, AlpsNodeStatus, double> >
293AbcTreeNode::branch()
294{
295    AbcNodeDesc* desc = 
296        dynamic_cast<AbcNodeDesc*>(desc_);
297    AbcModel* m = dynamic_cast<AbcModel*>(desc->getModel());
298
299    double* oldLbs = desc->lowerBounds();
300    double* oldUbs = desc->upperBounds();
301    const int numCols = m->getNumCols();
302    assert(oldLbs && oldUbs && numCols);
303
304    if ( (branchedOn_ < 0) || (branchedOn_ >= numCols) ) {
305        std::cout << "AbcError: branchedOn_ = "<< branchedOn_ << "; numCols = "
306                  << numCols << "; index_ = " << index_ << std::endl;
307        throw CoinError("branch index is out of range", 
308                        "branch", "AbcTreeNode");
309    }
310
311    double* newLbs = new double[numCols];
312    double* newUbs = new double[numCols];
313    std::copy(oldLbs, oldLbs + numCols, newLbs);
314    std::copy(oldUbs, oldUbs + numCols, newUbs);
315    //memcpy(newLbs, oldLbs, sizeof(double) * numCols);
316    //memcpy(newUbs, oldUbs, sizeof(double) * numCols);
317
318    std::vector< CoinTriple<AlpsNodeDesc*, AlpsNodeStatus, double> > newNodes;
319
320    double objVal;
321    // if (getParent() == 0) { //ROOT
322        objVal = getObjValue();
323        //}
324        //else {
325        //objVal = std::min(getParent()->getObjValue(), getObjValue());
326        //}
327   
328    // Branch down
329    newLbs[branchedOn_] = oldLbs[branchedOn_];
330    newUbs[branchedOn_] = floor(branchedOnVal_);//floor(branchedOnVal_+1.0e-5);
331    AbcNodeDesc* child;
332    assert(branchedOn_ >= 0);
333    child = new AbcNodeDesc(m, newLbs, newUbs);
334    child->setBranchedOn(branchedOn_);
335    child->setBranchedOnValue(branchedOnVal_);
336    child->setBranchedDir(-1);
337    newNodes.push_back(CoinMakeTriple(static_cast<AlpsNodeDesc *>(child),
338                                      AlpsNodeStatusCandidate,
339                                      objVal));
340   
341    // Branch up
342    newUbs[branchedOn_] = oldUbs[branchedOn_];
343    newLbs[branchedOn_] = ceil(branchedOnVal_);//ceil(branchedOnVal_ - 1.0e-5);
344    child = 0;
345    child = new AbcNodeDesc(m, newLbs, newUbs);
346    child->setBranchedOn(branchedOn_);
347    child->setBranchedOnValue(branchedOnVal_);
348    child->setBranchedDir(1);
349    newNodes.push_back(CoinMakeTriple(static_cast<AlpsNodeDesc *>(child),
350                                      AlpsNodeStatusCandidate,
351                                      objVal));
352    if (newLbs != 0) {
353        delete [] newLbs;
354        newLbs = 0;
355    }
356    if (newUbs != 0) {
357        delete [] newUbs;
358        newUbs = 0;
359    }
360    return newNodes;
361}
362
363//#############################################################################
364#if 0
365// Right now just choose the variable farest to be integral
366int
367AbcTreeNode::chooseBranch(AbcModel* model, bool& strongFound)
368{
369    strongFound = false;
370    int i = -1;
371    const int numInts = model->numberIntegers();
372    const int* intVar = model->integerVariable();
373    const double* sol = model->solver()->getColSolution();
374    double maxDist = model->getDblParam(AbcModel::AbcIntegerTolerance);
375    double value = 0.0;
376    double nearest = 0.0;
377    double distance = 0.0;
378
379    branchedOn_ = -1;
380
381    for (i = 0; i < numInts; ++i) {
382        value = sol[intVar[i]];
383        double nearest = floor(value + 0.5);
384        double distance = fabs(value - nearest);
385       
386        if (distance > maxDist) {
387            maxDist = distance;
388            branchedOn_ = intVar[i];
389            branchedOnVal_ = value;
390        }
391    }
392
393    if (branchedOn_ == -1) {
394        return -2;  // Should not happend
395    }
396    else {
397        return 2;   // Find a branch variable
398    }
399}
400#endif
401
402//#############################################################################
403//#############################################################################
404
405AlpsEncoded*
406AbcTreeNode::encode() const 
407{
408#if defined(ABC_DEBUG_MORE)
409    std::cout << "AbcTreeNode::encode()--start to encode" << std::endl;
410#endif
411    AlpsEncoded* encoded = new AlpsEncoded("ALPS_NODE");
412    AbcNodeDesc* desc = dynamic_cast<AbcNodeDesc*>(desc_);
413    AbcModel* model = dynamic_cast<AbcModel*>(desc->getModel());
414
415    int numCols = model->getNumCols();
416    assert(numCols);
417   
418    const double * lb = desc->lowerBounds();
419    const double * ub = desc->upperBounds();
420
421    encoded->writeRep(explicit_);
422    encoded->writeRep(numCols);
423    encoded->writeRep(lb, numCols);
424    encoded->writeRep(ub, numCols);
425    encoded->writeRep(index_);
426    encoded->writeRep(depth_);
427    encoded->writeRep(quality_);
428    encoded->writeRep(parentIndex_);
429    encoded->writeRep(numChildren_);
430    encoded->writeRep(status_);
431    encoded->writeRep(sentMark_);
432    encoded->writeRep(branchedOn_);
433    encoded->writeRep(branchedOnVal_);
434
435#if defined(ABC_DEBUG_MORE)
436    std::cout << "numCols = " << numCols << "; ";
437    for (int i = 0; i < numCols; ++i) {
438        std::cout << "[" << lb[i] << "," << ub[i] <<"], ";
439    }
440    std::cout << std::endl;
441    std::cout << "index_ = " << index_ << "; ";
442    std::cout << "depth_ = " << depth_ << "; ";
443    std::cout << "quality_ = " << quality_ << "; ";
444    std::cout << "parentIndex_ = " << parentIndex_ << "; ";
445    std::cout << "numChildren_ = " << numChildren_ << std::endl;
446#endif
447
448    return encoded;
449}
450
451AlpsKnowledge* 
452AbcTreeNode::decode(AlpsEncoded& encoded) const 
453{
454    int expli;
455    int numCols;
456    double* lb;
457    double* ub;
458    int index;
459    int depth;
460    double objValue;
461    int parentIndex;
462    int numChildren;
463    AlpsNodeStatus     nodeStatus;
464    int sentMark;
465    int branchedOn;
466    double branchedOnVal;
467   
468    encoded.readRep(expli);  // Check whether has full or diff
469    encoded.readRep(numCols);
470    encoded.readRep(lb, numCols);
471    encoded.readRep(ub, numCols);
472    encoded.readRep(index);
473    encoded.readRep(depth);
474    encoded.readRep(objValue);
475    encoded.readRep(parentIndex);
476    encoded.readRep(numChildren);
477    encoded.readRep(nodeStatus);
478    encoded.readRep(sentMark);
479    encoded.readRep(branchedOn);
480    encoded.readRep(branchedOnVal);
481   
482
483    if (nodeStatus == AlpsNodeStatusPregnant) {
484        assert(branchedOn >= 0 && branchedOn < numCols);
485    }
486       
487    AbcNodeDesc* nodeDesc = new AbcNodeDesc(dynamic_cast<AbcModel*>
488                                            (desc_->getModel()), lb, ub);
489 
490    //  nodeDesc->setModel(getKnowledgeBroker()->getDataPool()->getModel());
491    AbcTreeNode* treeNode = new AbcTreeNode(nodeDesc);
492    treeNode->setIndex(index);
493    treeNode->setDepth(depth);
494    treeNode->setObjValue(objValue);
495    treeNode->setParentIndex(parentIndex);
496    treeNode->setParent(0);
497    treeNode->setNumChildren(numChildren);
498    treeNode->setStatus(nodeStatus);
499    treeNode->setSentMark(sentMark);
500    treeNode->setBranchedOn(branchedOn);
501    treeNode->setBranchedOnValue(branchedOnVal);
502
503    if(!lb) {
504        delete [] lb;
505        lb = 0;
506    }
507    if(!ub) {
508        delete [] ub;
509        ub = 0;
510    }
511   
512#if defined(ABC_DEBUG_MORE)
513    std::cout << "numCols = " << numCols << "; ";
514    std::cout << "index = " << index << "; ";
515    std::cout << "depth = " << depth << "; ";
516    std::cout << "objValue = " << objValue<<"; ";
517    std::cout << "parentIndex = " << parentIndex << "; ";
518    std::cout << "status = " << nodeStatus << "; ";
519    std::cout << "branchedOn = " << branchedOn << "; ";
520    std::cout << "branchedOnVal = " << branchedOnVal << "; ";
521    std::cout << "numChildren = " << numChildren << std::endl;
522#endif
523
524    return treeNode;
525}
526
527//#############################################################################
528// This function tests each integer using strong branching and selects
529// the one with the least objective degradation. If strong branching is
530// disabled, the most infeasible integer is choosen
531// Returns:     2  Find a branching integer
532//             -1  Monotone
533//             -2  Infeasible
534int AbcTreeNode::chooseBranch(AbcModel *model, bool& strongFound)
535{ 
536    int numberIntegerInfeasibilities;
537    int anyAction = 0;
538       
539    strongFound = false;
540
541    // Make sure we are talking about the same solver
542    bool feasible = model->resolve(); 
543    if (!feasible) {
544        anyAction = -2;
545        return anyAction;
546    }
547
548    OsiSolverInterface * solver = model->solver();     
549    model->feasibleSolution(numberIntegerInfeasibilities);
550    if (numberIntegerInfeasibilities <= 0) {
551        double objectiveValue = 
552            solver->getObjSense() * solver->getObjValue();
553        strongFound = model->setBestSolution(ABC_STRONGSOL,
554                                             objectiveValue,
555                                             solver->getColSolution());
556        anyAction = -2;
557        return anyAction;
558    }
559   
560    double saveObjectiveValue = solver->getObjValue();
561    double objectiveValue = solver->getObjSense() * saveObjectiveValue;
562   
563    if(getKnowledgeBroker()->getProcRank() == -1) {
564        std::cout << "*** I AM RANK FIVE: obj = "<< objectiveValue
565                  << ", cutoff = "<< model->getCutoff() << std::endl;
566    }
567
568    if (objectiveValue >=  model->getCutoff()) {
569        anyAction = -2;
570        return anyAction;
571    }
572
573    const double * lower = solver->getColLower();
574    const double * upper = solver->getColUpper();
575
576    double integerTolerance = 
577        model->getDblParam(AbcModel::AbcIntegerTolerance);
578    int i, j;
579    bool beforeSolution = model->getSolutionCount()==0;
580    int numberStrong = model->numberStrong();
581    int numberObjects = model->numberIntegers();
582    int maximumStrong = max(min(model->numberStrong(), numberObjects), 1);
583    int numberColumns = model->getNumCols();
584    double * saveUpper = new double[numberColumns];
585    double * saveLower = new double[numberColumns];
586
587    // Save solution in case heuristics need good solution later
588    double * saveSolution = new double[numberColumns];
589    memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
590
591    // Get a branching decision object.
592    AbcBranchDecision * decision = model->branchingMethod();
593    if (!decision) {
594        decision = new AbcBranchDefaultDecision();
595    }
596   
597    typedef struct {
598        int possibleBranch;   // branching integer variable
599        double upMovement;    // cost going up (and initial away from feasible)
600        double downMovement;  // cost going down
601        int numIntInfeasUp;   // without odd ones
602        int numObjInfeasUp;   // just odd ones
603        bool finishedUp;      // true if solver finished
604        int numItersUp;       // number of iterations in solver
605        int numIntInfeasDown; // without odd ones
606        int numObjInfeasDown; // just odd ones
607        bool finishedDown;    // true if solver finished
608        int numItersDown;     // number of iterations in solver
609        int objectNumber;     // Which object it is
610    } Strong;
611
612    Strong * choice = new Strong[maximumStrong];
613    for (i = 0; i < numberColumns; ++i) {
614        saveLower[i] = lower[i];
615        saveUpper[i] = upper[i];
616    }
617
618    //assert(model->getForcePriority() < 0);     // For hot start
619
620    double saveOtherWay = 0.0;        // just for non-strong branching
621    numberUnsatisfied_ = 0;
622    int bestPriority = INT_MAX;
623
624    // Scan for branching objects that indicate infeasibility.
625    // Choose the best maximumStrong candidates, using priority as the
626    // first criteria, then integer infeasibility.
627    // The algorithm is to fill the choice array with a set of good
628    // candidates (by infeasibility) with priority bestPriority.  Finding
629    // a candidate with priority better (less) than bestPriority flushes
630    // the choice array. (This serves as initialization when the first
631    // candidate is found.)  When a candidate is added,
632    // it replaces the candidate with the smallest infeasibility (tracked by
633    // iSmallest)
634    int iSmallest;
635    double mostAway;
636    for (i = 0; i < maximumStrong; ++i) choice[i].possibleBranch = -1;
637    numberStrong = 0;
638    for (i = 0; i < numberObjects; ++i) {
639        const int object = model->integerVariable()[i];
640        int preferredWay;
641        double otherWay;
642        double infeasibility = checkInfeasibility(model, object, 
643                                                  preferredWay, otherWay);
644        if (infeasibility > integerTolerance) {
645            if (model->numberStrong() == 0) {  /* pseudo cost branching */
646                model->getPseudoIndices()[numberUnsatisfied_] = object;
647            }
648            ++numberUnsatisfied_;
649            int priorityLevel = model->priority(i);
650            if (priorityLevel < bestPriority) {
651                int j;
652                iSmallest = 0;
653                for (j = 0; j < maximumStrong; ++j) {
654                    choice[j].upMovement = 0.0;
655                    choice[j].possibleBranch = -1;
656                }
657                bestPriority = priorityLevel;
658                mostAway = integerTolerance;
659                numberStrong = 0;
660            } 
661            else if (priorityLevel > bestPriority) {
662                continue;
663            }
664     
665            // Check for suitability based on infeasibility.
666            if (infeasibility > mostAway) {
667                choice[iSmallest].upMovement = infeasibility;
668                choice[iSmallest].downMovement = otherWay;
669                saveOtherWay = otherWay;
670                choice[iSmallest].possibleBranch = object;
671                numberStrong = max(numberStrong, iSmallest + 1);
672                choice[iSmallest].objectNumber = i;
673                int j;
674                iSmallest = -1;
675                mostAway = 1.0e50;
676                for (j = 0; j < maximumStrong; ++j) {
677                    if (choice[j].upMovement < mostAway) {
678                        mostAway = choice[j].upMovement;
679                        iSmallest = j;
680                    }
681                }
682            }
683        }
684    }
685
686#if 0
687    if (numberStrong == 0) {  // FIXME: Should not happen
688        bool fea = model->feasibleSolution(numberIntegerInfeasibilities);
689        if (fea) {
690            double objectiveValue =
691                solver->getObjSense() * solver->getObjValue();
692            strongFound = model->setBestSolution(ABC_STRONGSOL,
693                                                 objectiveValue,
694                                                 solver->getColSolution());
695        }
696        //abort();
697        if (!model->branchingMethod()) delete decision;
698        delete [] choice;
699        delete [] saveLower;
700        delete [] saveUpper;
701        // restore solution
702        solver->setColSolution(saveSolution);
703        delete [] saveSolution;
704        anyAction = -2;
705        return anyAction;
706    }
707#endif   
708
709    // Some solvers can do the strong branching calculations faster if
710    // they do them all at once.  At present only Clp does for ordinary
711    // integers but I think this coding would be easy to modify
712    bool allNormal = true; // to say if we can do fast strong branching
713    for (i = 0; i < numberStrong; ++i) {
714        choice[i].numIntInfeasUp = numberUnsatisfied_;
715        choice[i].numIntInfeasDown = numberUnsatisfied_;
716    }
717
718    // Setup for strong branching involves saving the current basis (for
719    // restoration afterwards) and setting up for hot starts.
720    if (model->numberStrong() >0) {
721        // set true to say look at all even if some fixed (experiment)
722        bool solveAll = false; 
723        // Save basis
724        CoinWarmStart * ws = solver->getWarmStart();
725        // save limit
726        int saveLimit;
727        solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
728        if (beforeSolution)  // go to end
729            solver->setIntParam(OsiMaxNumIterationHotStart, 10000); 
730
731        // If we are doing all strong branching in one go then we create
732        // new arrays to store information. If clp NULL then doing old way.
733        // Going down -
734        //  outputSolution[2*i] is final solution.
735        //  outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown)
736        //  outputStuff[2*i+numberStrong] is number iterations
737        //  On entry newUpper[i] is new upper bound, on exit obj change
738        // Going up -
739        //  outputSolution[2*i+1] is final solution.
740        //  outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
741        //  outputStuff[2*i+1+numberStrong] is number iterations
742        //  On entry newLower[i] is new lower bound, on exit obj change
743        OsiClpSolverInterface * osiclp = 
744            dynamic_cast< OsiClpSolverInterface*>(solver);
745        ClpSimplex * clp=NULL;
746        double * newLower = NULL;
747        double * newUpper = NULL;
748        double ** outputSolution = NULL;
749        int * outputStuff = NULL;
750        int saveLogLevel;
751
752        //allNormal=false;
753        if (osiclp && allNormal) {
754            clp = osiclp->getModelPtr();
755            saveLogLevel = clp->logLevel();
756            int saveMaxIts = clp->maximumIterations();
757            clp->setMaximumIterations(saveLimit);
758            clp->setLogLevel(0);
759            newLower = new double[numberStrong];
760            newUpper = new double[numberStrong];
761            outputSolution = new double * [2 * numberStrong];
762            outputStuff = new int [4 * numberStrong];
763            int * which = new int [numberStrong];
764            for (i = 0; i < numberStrong; ++i) {
765                int iObject = choice[i].objectNumber;
766                int iSequence = model->integerVariable()[iObject];
767                newLower[i] = ceil(saveSolution[iSequence]);
768                newUpper[i] = floor(saveSolution[iSequence]);
769                which[i] = iSequence;
770                outputSolution[2*i] = new double [numberColumns];
771                outputSolution[2*i+1] = new double [numberColumns];
772            }
773            clp->strongBranching(numberStrong, which,
774                                 newLower, newUpper, 
775                                 outputSolution, outputStuff,
776                                 outputStuff + 2 * numberStrong,
777                                 !solveAll, false);
778            clp->setMaximumIterations(saveMaxIts);
779            delete [] which;
780        } 
781        else { // Doing normal way
782            solver->markHotStart();
783        }
784
785        //---------------------------------------------------------------------
786        // Open a loop to do the strong branching LPs. For each candidate
787        // variable, solve an LP with the variable forced down, then up.
788        // If a direction turns out to be infeasible or monotonic (i.e.,
789        // over the dual objective cutoff), force the objective change to
790        // be big (1.0e100). If we determine the problem is infeasible,
791        // or find a monotone variable, escape the loop.
792
793        for (i = 0; i < numberStrong; ++i) { 
794            double objectiveChange ;
795            double newObjectiveValue = 1.0e100;
796            // status is 0 finished, 1 infeasible and other
797            int iStatus;
798            int colInd = choice[i].possibleBranch;
799
800            // Try the down direction first.
801            if (!clp) {
802                solver->setColUpper(colInd, floor(saveSolution[colInd]));
803                solver->solveFromHotStart();
804                // restore bounds
805                for (int j = 0; j < numberColumns; ++j) {
806                    if (saveLower[j] != lower[j])
807                        solver->setColLower(j, saveLower[j]);
808                    if (saveUpper[j] != upper[j])
809                        solver->setColUpper(j, saveUpper[j]);
810                }
811
812                if (solver->isProvenOptimal())
813                    iStatus = 0; // optimal
814                else if (solver->isIterationLimitReached()
815                         &&!solver->isDualObjectiveLimitReached())
816                    iStatus = 2; // unknown
817                else
818                    iStatus = 1; // infeasible
819                newObjectiveValue = 
820                    solver->getObjSense() * solver->getObjValue();
821                choice[i].numItersDown = solver->getIterationCount();
822                //objectiveChange = newObjectiveValue - objectiveValue ;
823            } 
824            else {
825                iStatus = outputStuff[2*i];
826                choice[i].numItersDown = outputStuff[2*numberStrong + 2*i];
827                newObjectiveValue = objectiveValue + newUpper[i];
828                solver->setColSolution(outputSolution[2*i]);
829            }
830            objectiveChange = newObjectiveValue - objectiveValue;
831            if (!iStatus) {
832                choice[i].finishedDown = true;
833                if (newObjectiveValue > model->getCutoff()) {
834                    objectiveChange = 1.0e100;          // discard it
835                } 
836                else {
837                    // See if integer solution
838                    if (model->feasibleSolution(choice[i].numIntInfeasDown)){
839                        strongFound = 
840                            model->setBestSolution(ABC_STRONGSOL,
841                                                   newObjectiveValue,
842                                                   solver->getColSolution());
843                        if (newObjectiveValue > model->getCutoff())
844                            objectiveChange = 1.0e100;  // discard it
845                    }
846                }
847            } 
848            else if (iStatus == 1) {
849                objectiveChange = 1.0e100;
850            } 
851            else {              // Can't say much as we did not finish
852                choice[i].finishedDown = false ;
853            }
854            choice[i].downMovement = objectiveChange ;
855           
856            // repeat the whole exercise, forcing the variable up
857            if (!clp) {
858                solver->setColLower(colInd, ceil(saveSolution[colInd]));
859                solver->solveFromHotStart();
860                // restore bounds
861                for (int j = 0; j < numberColumns; j++) {
862                    if (saveLower[j] != lower[j])
863                        solver->setColLower(j, saveLower[j]);
864                    if (saveUpper[j] != upper[j])
865                        solver->setColUpper(j, saveUpper[j]);
866                }
867                if (solver->isProvenOptimal())
868                    iStatus = 0; // optimal
869                else if (solver->isIterationLimitReached()
870                         &&!solver->isDualObjectiveLimitReached())
871                    iStatus = 2; // unknown
872                else
873                    iStatus = 1; // infeasible
874                newObjectiveValue = 
875                    solver->getObjSense() * solver->getObjValue();
876                choice[i].numItersUp = solver->getIterationCount();
877                objectiveChange = newObjectiveValue - objectiveValue ;
878            } 
879            else {
880                iStatus = outputStuff[2 * i + 1];
881                choice[i].numItersUp = outputStuff[2*numberStrong + 2*i + 1];
882                newObjectiveValue = objectiveValue + newLower[i];
883                solver->setColSolution(outputSolution[2*i + 1]);
884            }
885            objectiveChange = newObjectiveValue - objectiveValue;
886            if (!iStatus) {
887                choice[i].finishedUp = true;
888                if (newObjectiveValue > model->getCutoff()) {
889                    objectiveChange = 1.0e100;
890                } 
891                else {
892                    // See if integer solution
893                    if (model->feasibleSolution(choice[i].numIntInfeasUp)){
894                        strongFound = 
895                            model->setBestSolution(ABC_STRONGSOL,
896                                                   newObjectiveValue,
897                                                   solver->getColSolution());
898                        if (newObjectiveValue > model->getCutoff())
899                            objectiveChange = 1.0e100;
900                    }
901                }
902            } 
903            else if (iStatus == 1) {
904                objectiveChange = 1.0e100;
905            } 
906            else {
907                // Can't say much as we did not finish
908                choice[i].finishedUp = false;
909            }
910            choice[i].upMovement = objectiveChange;
911
912            // End of evaluation for this candidate variable. Possibilities
913            // are:
914            // * Both sides below cutoff; this variable is a candidate
915            //   for branching.
916            // * Both sides infeasible or above the objective cutoff:
917            //   no further action here. Break from the evaluation loop and
918            //   assume the node will be purged by the caller.
919            // * One side below cutoff: Install the branch (i.e., fix the
920            //   variable). Break from the evaluation loop and assume
921            //   the node will be reoptimised by the caller.
922
923            AbcNodeDesc* desc = dynamic_cast<AbcNodeDesc*>(desc_);
924            int monoIndex;         
925
926            if (choice[i].upMovement < 1.0e100) {
927                if(choice[i].downMovement < 1.0e100) {
928                    anyAction = 2;   // Both
929                } 
930                else {               // Only up
931                    anyAction = -1; 
932                    monoIndex = choice[i].possibleBranch;
933                    solver->setColLower(monoIndex,
934                                        ceil(saveSolution[monoIndex]));
935                    desc->setLowerBound(monoIndex, 
936                                        ceil(saveSolution[monoIndex]));
937                    break;
938                }
939            } 
940            else {
941                if(choice[i].downMovement < 1.0e100) {
942                    anyAction = -1;  // Only down
943                    monoIndex = choice[i].possibleBranch;
944                    solver->setColUpper(monoIndex,
945                                        floor(saveSolution[monoIndex]));
946                    desc->setUpperBound(monoIndex, 
947                                        floor(saveSolution[monoIndex]));
948                    break;
949                } 
950                else {                // neither
951                    anyAction = -2;
952                    break;
953                }
954            }
955        }
956       
957        if (!clp) {
958            solver->unmarkHotStart();
959        } else {
960            clp->setLogLevel(saveLogLevel);
961            delete [] newLower;
962            delete [] newUpper;
963            delete [] outputStuff;
964            for (int i = 0; i < 2*numberStrong; ++i)
965                delete [] outputSolution[i];
966            delete [] outputSolution;
967        }
968
969        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
970        solver->setWarmStart(ws);       // restore basis
971        delete ws;
972
973        if(getKnowledgeBroker()->getProcRank() == -1) {
974            std::cout << "*** I AM RANK ONE: chooseBranch():anyAction = " 
975                      << anyAction << "numberStrong = "<< numberStrong
976                      << std::endl;
977        }
978
979        // Sift through the candidates for the best one.
980        // QUERY: Setting numberNodes looks to be a distributed noop.
981        // numberNodes is local to this code block. Perhaps should be
982        // numberNodes_ from model?
983        // Unclear what this calculation is doing.
984        if (anyAction > 0) {
985            int numberNodes = model->getNodeCount();
986            // get average cost per iteration and assume stopped ones
987            // would stop after 50% more iterations at average cost??? !!! ???
988            double averageCostPerIteration = 0.0;
989            double totalNumberIterations = 1.0;
990            int smallestNumberInfeasibilities = INT_MAX;
991            for (i = 0; i < numberStrong; ++i) {
992                totalNumberIterations += choice[i].numItersDown +
993                    choice[i].numItersUp;
994                averageCostPerIteration += choice[i].downMovement +
995                    choice[i].upMovement;
996                smallestNumberInfeasibilities = 
997                    min(min(choice[i].numIntInfeasDown,
998                            choice[i].numIntInfeasUp),
999                        smallestNumberInfeasibilities);
1000            }
1001            if (smallestNumberInfeasibilities >= numberIntegerInfeasibilities)
1002                numberNodes = 1000000; // switch off search for better solution
1003            numberNodes = 1000000;     // switch off anyway
1004            averageCostPerIteration /= totalNumberIterations;
1005            // all feasible - choose best bet
1006
1007            // New method does all at once so it can be more sophisticated
1008            // in deciding how to balance actions.
1009            // But it does need arrays
1010            double * changeUp = new double [numberStrong];
1011            int * numberInfeasibilitiesUp = new int [numberStrong];
1012            double * changeDown = new double [numberStrong];
1013            int * numberInfeasibilitiesDown = new int [numberStrong];
1014            int * objects = new int [numberStrong];
1015
1016            for (i = 0; i < numberStrong; ++i) {
1017
1018                int iColumn = choice[i].possibleBranch;
1019
1020                model->messageHandler()->message(ABC_STRONG,model->messages())
1021                    << i << iColumn
1022                    << choice[i].downMovement << choice[i].numIntInfeasDown
1023                    << choice[i].upMovement << choice[i].numIntInfeasUp
1024                    << choice[i].possibleBranch
1025                    << CoinMessageEol;
1026
1027                changeUp[i] = choice[i].upMovement;
1028                numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1029                changeDown[i] = choice[i].downMovement;
1030                numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1031                objects[i] = choice[i].possibleBranch;
1032            }
1033            int whichObject = 
1034                decision->bestBranch(model,
1035                                     objects, numberStrong,
1036                                     numberUnsatisfied_, changeUp,
1037                                     numberInfeasibilitiesUp, changeDown,
1038                                     numberInfeasibilitiesDown,objectiveValue);
1039
1040            // move branching object and make sure it will not be deleted
1041            if (whichObject >= 0) {
1042                branchedOn_ = objects[whichObject];
1043                branchedOnVal_= saveSolution[branchedOn_];
1044                assert(branchedOn_ >= 0 && branchedOn_ < numberColumns);
1045
1046#if 0           // FIXME
1047                std::cout << "AbcStrongBranch: col index : ";
1048                for (i = 0; i < numberStrong; ++i) {
1049                    std::cout<< objects[i] << " ";
1050                }
1051                std::cout <<"\nAbcStrongBranch: branchedOn_ = " << branchedOn_
1052                          << "; branchedOnVal_ = " << branchedOnVal_
1053                          << std::endl;
1054#endif
1055            }
1056            else {
1057                std::cout << "AbcError: whichObject = " << whichObject
1058                          << "; numberStrong = " << numberStrong << std::endl;
1059                throw CoinError("bestBranch find nothing", 
1060                                "chooseBranch", 
1061                                "AbcTreeNode");
1062            }
1063
1064            delete [] changeUp;
1065            delete [] numberInfeasibilitiesUp;
1066            delete [] changeDown;
1067            delete [] numberInfeasibilitiesDown;
1068            delete [] objects;
1069        }
1070       
1071        if (anyAction == 0){
1072            throw CoinError("Can't find candidates", 
1073                            "chooseBranch", 
1074                            "AbcTreeNode");
1075        }
1076    }
1077    else {  // pseudocost branching
1078#if 1
1079        int object;
1080        int mostInd = -1;
1081        double mostDeg = -1.0;
1082        double deg, upCost, downCost, fraction;
1083        AbcPseudocost *pseudoC;
1084       
1085        for(j = 0; j < numberUnsatisfied_; ++j) {
1086            //for(j = 0; j < numberStrong; ++j) {
1087            object = model->getPseudoIndices()[j];
1088            //object = choice[j].possibleBranch;
1089            fraction = saveSolution[object] - floor(saveSolution[object]);
1090            pseudoC = model->getPseudoList()[object];
1091            upCost = pseudoC->upCost_ * (1.0 - fraction);
1092            downCost = pseudoC->downCost_ * fraction;
1093            deg = 4.0 * std::min(upCost, downCost) +
1094                1.0 * std::max(upCost, downCost);
1095            if (deg > mostDeg) {
1096                mostDeg = deg;
1097                mostInd = object;
1098            }
1099        }
1100       
1101        branchedOn_ = mostInd;
1102#else
1103        branchedOn_ = choice[0].possibleBranch;
1104#endif 
1105        branchedOnVal_= saveSolution[branchedOn_];
1106        assert( (branchedOn_ >= 0) && (branchedOn_ < numberColumns));
1107    }
1108
1109    if (!model->branchingMethod()) delete decision;
1110
1111    delete [] choice;
1112    delete [] saveLower;
1113    delete [] saveUpper;
1114
1115    // restore solution
1116    solver->setColSolution(saveSolution);
1117    delete [] saveSolution;
1118
1119    return anyAction;
1120}
1121
Note: See TracBrowser for help on using the repository browser.