source: trunk/Couenne/src/branch/CouenneChooseStrong.cpp @ 593

Last change on this file since 593 was 593, checked in by fmargot, 9 years ago

Make sure artificial cutoff from option file works; fix small bug in CouenneChooseStrong?.cpp; Add ifdef FM_ALWAYS_SORT

  • Property svn:keywords set to Author Date Id Revision
File size: 48.3 KB
Line 
1/* $Id: CouenneChooseStrong.cpp 593 2011-06-01 12:16:35Z fmargot $
2 *
3 * Name:    CouenneChooseStrong.cpp
4 * Authors: Andreas Waechter, IBM Corp.
5 *          Pietro Belotti, Lehigh University
6 * Purpose: Strong branching objects for Couenne
7 *
8 * (C) Carnegie-Mellon University, 2006-10.
9 * This file is licensed under the Eclipse Public License (EPL)
10 */
11
12#include "CouenneObject.hpp"
13#include "BonChooseVariable.hpp"
14#include "CouenneChooseStrong.hpp"
15#include "CouenneProblem.hpp"
16#include "CouenneProblemElem.hpp"
17#include "CouenneBranchingObject.hpp"
18
19#include "CouenneRecordBestSol.hpp"
20
21//#define TRACE_STRONG
22//#define TRACE_STRONG2
23//#define FM_MOD
24//#define FM_SORT_STRONG
25//#define FM_ALWAYS_SORT
26//#define OLD_STYLE
27
28using namespace Couenne;
29
30const CouNumber estProdEps = 1e-6;
31
32
33  /// constructor
34  CouenneChooseStrong::CouenneChooseStrong (Bonmin::BabSetupBase &b, CouenneProblem* p, JnlstPtr jnlst) :
35
36  Bonmin::BonChooseVariable (b, b.continuousSolver()),
37  problem_          (p),
38  jnlst_            (jnlst),
39  branchtime_       (0.) {
40
41    std::string s;
42
43    b.options () -> GetStringValue ("pseudocost_mult_lp", s, "couenne.");
44    pseudoUpdateLP_ = (s == "yes");     
45
46    b.options () -> GetStringValue ("estimate_select", s, "couenne.");
47    estimateProduct_ = (s == "product");
48
49    b.options () -> GetStringValue ("trust_strong", s, "couenne.");
50
51    // trust solution from strong branching to provide right LB
52
53    setTrustStrongForSolution (s == "yes");
54    setTrustStrongForBound    (s == "yes");
55
56    minDepthPrint_ = -1;
57  }
58
59  /// copy constructor
60  CouenneChooseStrong::CouenneChooseStrong (const CouenneChooseStrong& rhs) :
61    Bonmin::BonChooseVariable (rhs),
62    problem_          (rhs.problem_),
63    pseudoUpdateLP_   (rhs.pseudoUpdateLP_),
64    estimateProduct_  (rhs.estimateProduct_),
65    jnlst_            (rhs.jnlst_),
66    branchtime_       (rhs.branchtime_),
67    minDepthPrint_    (rhs.minDepthPrint_)
68  {}
69
70  /// destructor
71  CouenneChooseStrong::~CouenneChooseStrong()
72  {if (branchtime_ > 1e-9) jnlst_ -> Printf (J_ERROR, J_BRANCHING, "Strong branching: total time %g\n", branchtime_);}
73
74  /// cloning method
75  OsiChooseVariable *
76  CouenneChooseStrong::clone() const
77  {return new CouenneChooseStrong(*this);}
78
79  /// assignment operator
80  CouenneChooseStrong&
81  CouenneChooseStrong::operator=(const CouenneChooseStrong & rhs)
82  {
83    if (this != &rhs) {
84      Bonmin::BonChooseVariable::operator=(rhs);
85      problem_         = rhs.problem_;
86      pseudoUpdateLP_  = rhs.pseudoUpdateLP_;
87      estimateProduct_ = rhs.estimateProduct_;
88      jnlst_           = rhs.jnlst_;
89      branchtime_      = rhs.branchtime_;
90      minDepthPrint_   = rhs.minDepthPrint_;
91    }
92    return *this;
93  }
94
95
96  /* Choose a variable. Returns:
97
98    -1  Node is infeasible
99     0  Normal termination - we have a candidate
100     1  All look satisfied - no candidate
101     2  We can change the bound on a variable - but we also have a strong branching candidate
102     3  We can change the bound on a variable - but we have a non-strong branching candidate
103     4  We can change the bound on a variable - no other candidates
104
105     We can pick up branch from whichObject() and whichWay()
106     We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
107     If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
108  */
109#ifdef FM_MOD
110/*******************************************************************/
111  int CouenneChooseStrong::chooseVariable (OsiSolverInterface * solver,
112                                           OsiBranchingInformation *info,
113                                           bool fixVariables) {
114
115    /// Note: had to copy code from
116    /// BonChooseVariable::chooseVariable() in order to test product
117    /// thing
118
119    problem_ -> domain () -> push
120      (problem_ -> nVars (),
121       info -> solution_,
122       info -> lower_,
123       info -> upper_);
124
125#ifdef TRACE_STRONG2
126    int pv = -1;
127    if(info->depth_ > minDepthPrint_) {
128      if(pv > -1) {
129        printf("CCS: beg: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
130               pv, solver->getColSolution()[pv], solver->getColLower()[pv], solver->getColUpper()[pv]);
131        printf("CCS: info: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
132               pv, info->solution_[pv], info->lower_[pv], info->upper_[pv]);
133        printf("CCS: problem: lb: %10.4f  ub: %10.4f\n",
134               problem_->Lb(pv), problem_->Ub(pv));
135      }
136    }
137 #endif                                                                         
138
139    int retval;
140
141    //int retval = BonChooseVariable::chooseVariable (solver, info, fixVariables);
142
143    // COPY of Bonmin starts here ////////////////////////////////////////////
144
145    // We assume here that chooseVariable is called once at the very
146    // beginning with fixVariables set to true.  This is then the root
147    // node.
148
149    bool isRoot = isRootNode(info);
150    int numberStrong = numberStrong_;
151    if (isRoot) {
152      numberStrong = CoinMax(numberStrong_, numberStrongRoot_);
153    }
154    if (numberUnsatisfied_) {
155      int cardIndForPseudo = 0, *indForPseudo = new int[numberUnsatisfied_];
156      OsiObject ** object = solver->objects();
157      const double* upTotalChange = pseudoCosts_.upTotalChange();
158      const double* downTotalChange = pseudoCosts_.downTotalChange();
159      const int* upNumber = pseudoCosts_.upNumber();
160      const int* downNumber = pseudoCosts_.downNumber();
161      int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
162      int numberLeft = CoinMin(numberStrong -numberStrongDone_,numberUnsatisfied_);
163      results_.clear();
164      int returnCode = 0;
165      int returnCodeSB = 0;
166      bestObjectIndex_ = -1;
167      bestWhichWay_ = -1;
168      firstForcedObjectIndex_ = -1;
169      firstForcedWhichWay_ =-1;
170      double bestTrusted=-COIN_DBL_MAX;
171      for (int i=0;i<numberLeft;i++) {
172        int iObject = list_[i];
173        if (numberBeforeTrusted == 0||
174            i < minNumberStrongBranch_ ||
175            (
176              only_pseudo_when_trusted_ && number_not_trusted_>0 ) ||
177              ( !isRoot && (upNumber[iObject]<numberBeforeTrusted ||
178                          downNumber[iObject]<numberBeforeTrusted ))||
179              ( isRoot && (!upNumber[iObject] && !downNumber[iObject])) ) {
180         
181#ifdef TRACE_STRONG
182          if(info->depth_ > minDepthPrint_) {
183            printf("Push object %d for strong branch\n", iObject);
184          }
185#endif
186          results_.push_back(Bonmin::HotInfo(solver, info,
187                                object, iObject));
188        }
189        else {
190
191#ifdef TRACE_STRONG
192          if(info->depth_ > minDepthPrint_) {
193            printf("Use pseudo cost for object %d\n", iObject);
194          }
195#endif
196          indForPseudo[cardIndForPseudo] = iObject;
197          cardIndForPseudo++;
198        }
199      }
200
201      int numberFixed=0;
202
203      if (results_.size() > 0) {
204
205        // do strong branching //////////////////////////////////////////////////
206        returnCodeSB = doStrongBranching (solver, info, results_.size(), 1);
207
208        if (bb_log_level_>=3) {
209          const char* stat_msg[] = {"NOTDON", "FEAS", "INFEAS", "NOFINI"};
210          message(SB_HEADER)<<CoinMessageEol;
211          for (unsigned int i = 0; i< results_.size(); i++) {
212            double up_change = results_[i].upChange();
213            double down_change = results_[i].downChange();
214            int up_status = results_[i].upStatus();
215            int down_status = results_[i].downStatus();
216            message(SB_RES)<<(int) i<<stat_msg[down_status+1]<<down_change
217            <<stat_msg[up_status+1]<< up_change<< CoinMessageEol;
218          }
219        }
220
221        if (returnCodeSB >= 0 && returnCodeSB <= 2) { // 1, 2: some can be fixed
222          if(returnCodeSB > 0) {
223            returnCode = 4; // no bestObject yet
224          }
225          else {
226            returnCode = 0;
227          }
228          const double prec = problem_->getFeasTol();
229
230          for (unsigned int i=0;i < results_.size (); i++) {
231
232            if((results_[i].upStatus() < 0) || (results_[i].downStatus() < 0))
233              continue;
234
235
236            if((results_[i].upStatus() < 0) || (results_[i].downStatus() < 0)) {
237              continue;
238            }
239
240            int iObject = results_[i].whichObject();
241
242            ///////////////////////////////////////////////////////////////////
243
244            double upEstimate;
245
246            if (results_[i].upStatus()!=1) {
247              assert (results_[i].upStatus()>=0);
248              upEstimate = results_[i].upChange();
249            }
250            else {
251              // infeasible - just say expensive
252              if (info->cutoff_<1.0e50)
253                upEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
254              else
255                upEstimate = 2.0*fabs(info->objectiveValue_);
256              if (firstForcedObjectIndex_ <0) {
257                // first fixed variable
258                firstForcedObjectIndex_ = iObject;
259                firstForcedWhichWay_ =0;
260              }
261
262              numberFixed++;
263              if (fixVariables) {
264                bool needBranch = true;
265                const OsiObject * obj = object[iObject];
266                CouenneObject *co = dynamic_cast <CouenneObject *>(object[iObject]);
267                OsiSimpleInteger *simpl = dynamic_cast <OsiSimpleInteger *>(object[iObject]);
268               int vInd = -1;
269                if(co) {
270                  vInd = co->Reference()->Index();
271                }
272                else {
273                  vInd = obj->columnNumber();
274                }
275                if((vInd >= 0) && (simpl)) {
276                  needBranch = false;
277                  double nearest = floor(info->solution_[vInd] + 0.5);
278                  if(nearest > 0.5) {
279                    nearest = 1 - nearest;
280                  }
281                                                                               
282                  if((nearest > info->integerTolerance_) &&
283                     (solver->getColUpper()[vInd] > solver->getColLower()[vInd]\
284 + prec)) {
285                    needBranch = true;
286                  }
287                }
288                if(needBranch) {
289                  OsiBranchingObject * branch = obj -> createBranch (solver, info, 0);
290                  branch -> branch (solver);
291                  delete branch;
292                }
293              }
294            }
295
296            ///////////////////////////////////////////////////////////////////
297
298            double downEstimate;
299
300            if (results_[i].downStatus()!=1) {
301              assert (results_[i].downStatus()>=0);
302              downEstimate = results_[i].downChange();
303            }
304            else {
305              // infeasible - just say expensive
306              if (info->cutoff_<1.0e50)
307                downEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
308              else
309                downEstimate = 2.0*fabs(info->objectiveValue_);
310              if (firstForcedObjectIndex_ <0) {
311                firstForcedObjectIndex_ = iObject;
312                firstForcedWhichWay_ =1;
313              }
314              numberFixed++;
315              if (fixVariables) {
316                bool needBranch = true;
317                const OsiObject * obj = object[iObject];
318                CouenneObject *co = dynamic_cast <CouenneObject *>(object[iObject]);
319                OsiSimpleInteger *simpl = dynamic_cast <OsiSimpleInteger *>(object[iObject]);
320                int vInd = -1;
321                if(co) {
322                  vInd = co->Reference()->Index();
323                }
324                else {
325                  vInd = obj->columnNumber();
326                }
327                if((vInd >= 0) && (simpl)) {
328                  needBranch = false;
329                  double nearest = floor(info->solution_[vInd] + 0.5);
330                  if(nearest > 0.5) {
331                    nearest = 1 - nearest;
332                  }
333                                                                               
334                  if((nearest > info->integerTolerance_) &&
335                     (solver->getColUpper()[vInd] > solver->getColLower()[vInd]\
336 + prec)) {
337                    needBranch = true;
338                  }
339                }
340                if(needBranch) {
341                  OsiBranchingObject * branch = obj -> createBranch (solver, info, 1);
342                  branch -> branch (solver);
343                  delete branch;
344                }
345              }
346            }
347
348            double
349              MAXMIN_CRITERION = maxminCrit (info),
350              minVal, maxVal, value;
351
352            if (downEstimate < upEstimate) {
353              minVal = downEstimate;
354              maxVal = upEstimate;
355            } else {
356              minVal = upEstimate;
357              maxVal = downEstimate;
358            }
359
360            value = 
361              estimateProduct_ ? 
362              ((estProdEps + minVal) * maxVal) :
363              (       MAXMIN_CRITERION  * minVal + 
364               (1.0 - MAXMIN_CRITERION) * maxVal);
365
366            if (value>bestTrusted) {
367              bestTrusted = value;
368              bestObjectIndex_ = iObject;
369              bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
370              // but override if there is a preferred way
371              const OsiObject * obj = object[iObject];
372              if (obj->preferredWay()>=0&&obj->infeasibility())
373                bestWhichWay_ = obj->preferredWay();
374              if(returnCodeSB) { // 1 or 2
375                returnCode = 2; 
376              }
377            }
378          }
379        }
380        else { // returnCodeSB is -1 or 3
381          if (returnCodeSB == 3) { // max time - just choose one
382            if(bestObjectIndex_ < 0) {
383              bestObjectIndex_ = list_[0];
384              bestWhichWay_ = 0;
385              returnCode = 0;
386            }
387          }
388          else {
389            returnCode = -1;
390          }
391        }
392#ifdef OLD_STYLE
393        if(bestObjectIndex_ < 0) {
394          bestObjectIndex_ = list_[0];
395        }
396        cardIndForPseudo = 0; // do not scan other vars using pseudo costs
397#endif
398      }
399
400      if(returnCodeSB != -1) { 
401                  // if returnCodeSB == -1 (i.e. problem is infeasible)
402                  // no need to scan objects with pseudocosts
403        const double *solverUb = solver->getColUpper();
404        const double *solverLb = solver->getColLower();
405        const double prec = problem_->getFeasTol();
406       
407        // to keep best object skipped on basis of bounds and branching value
408        int bestObjectIndex2 = -1;
409        int bestWhichWay2 = 0;
410        double bestTrusted2 = -COIN_DBL_MAX;
411       
412        for(int ips=0; ips<cardIndForPseudo; ips++) {
413          int iObject = indForPseudo[ips];
414          const OsiObject * obj = object[iObject];
415          double
416            upEstimate       = (upTotalChange[iObject]*obj->upEstimate())/upNumber[iObject],
417            downEstimate     = (downTotalChange[iObject]*obj->downEstimate())/downNumber[iObject],
418            MAXMIN_CRITERION = maxminCrit (info),
419            minVal, maxVal, value;
420         
421          if (downEstimate < upEstimate) {
422            minVal = downEstimate;
423            maxVal = upEstimate;
424          } else {
425            minVal = upEstimate;
426            maxVal = downEstimate;
427          }
428         
429          value = 
430            estimateProduct_ ? 
431            ((estProdEps + minVal) * maxVal) :
432            (       MAXMIN_CRITERION  * minVal + 
433                    (1.0 - MAXMIN_CRITERION) * maxVal);
434         
435         
436          // skip OsiSimpleInteger objects with variable fixed or
437          // branching value outside bounds
438          bool skipIt = false;
439          CouenneObject *co = dynamic_cast <CouenneObject *>(object[iObject]);
440          OsiSimpleInteger *simpl = dynamic_cast <OsiSimpleInteger *>(object[iObject]);
441          int vInd = -1;
442          if(co) {
443            vInd = co->Reference()->Index();
444          }
445          else {
446            vInd = obj->columnNumber();
447          }
448          if(simpl && (vInd >= 0)) {
449            double vUb = solverUb[vInd];
450            double vLb = solverLb[vInd];
451            double vSol = info->solution_[vInd];
452            if((vSol < vLb + prec) || (vSol > vUb - prec) || (vUb-vLb < prec)) {
453              skipIt = true;
454              numberFixed++;
455            }
456          }
457         
458          if(skipIt) {
459            if (value > bestTrusted2) {
460              bestObjectIndex2 = iObject;
461              bestWhichWay2 = upEstimate>downEstimate ? 0 : 1;
462              bestTrusted2 = value;
463            }
464#ifdef TRACE_STRONG
465            if(info->depth_ > minDepthPrint_) {
466              printf("Object %d skip pseudocost\n", iObject);
467            }
468#endif
469          }
470          else {
471            if (value > bestTrusted) {
472              bestObjectIndex_ = iObject;
473              bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
474              bestTrusted = value;
475              returnCode = (returnCode ? 3 : 0); // if returnCode was 2 or 3
476                                                 // it becomes 3
477            }
478#ifdef TRACE_STRONG
479            if(info->depth_ > minDepthPrint_) {
480              printf("Object %d use pseudocost\n", iObject);
481            }
482#endif
483          }
484        }
485   
486        delete[] indForPseudo;
487
488        if((bestObjectIndex_ < 0) && (bestObjectIndex2 >= 0)) {
489          bestObjectIndex_ = bestObjectIndex2;
490          bestWhichWay_ = bestWhichWay2;
491          bestTrusted = bestTrusted2;
492          returnCode = 4;
493        }
494      }
495
496      int objectVarInd = -1;
497      if(bestObjectIndex_ >= 0) {
498        OsiObject * obj = object[bestObjectIndex_];
499        obj->setWhichWay(bestWhichWay_);
500        CouenneObject *co =  dynamic_cast <CouenneObject *>(object[bestObjectIndex_]);
501        if(co) {
502          objectVarInd = co->Reference()->Index();
503        }
504        else {
505          objectVarInd = obj->columnNumber();
506        }
507
508        message(BRANCH_VAR)<<bestObjectIndex_<< bestWhichWay_ <<CoinMessageEol;
509
510#ifdef TRACE_STRONG
511        if(info->depth_ > minDepthPrint_) {
512          if(objectVarInd >= 0) {
513            double vUb = solver->getColUpper()[objectVarInd];
514            double vLb = solver->getColLower()[objectVarInd];
515            double vSolI = info->solution_[objectVarInd];
516            double vSolS = solver->getColSolution()[objectVarInd];
517            printf("Branch on object %d (var: %d): solInfo: %10.4f  SolSolver: %10.4f low: %10.4f  up: %10.4f\n",
518                   bestObjectIndex_, objectVarInd, vSolI, vSolS, vLb, vUb);
519          }
520          else {
521            printf("Branch on object %d (var: -1)\n", bestObjectIndex_);
522          }
523        }
524#endif
525      }
526      message(CHOSEN_VAR)<<bestObjectIndex_<<CoinMessageEol;
527
528      if (numberFixed==numberUnsatisfied_&&numberFixed)
529        returnCode=4;
530
531      if((returnCode == 2) || (returnCode == 3)) {
532        if((objectVarInd > -1) && 
533           (solver->getColUpper()[objectVarInd] < solver->getColLower()[objectVarInd] + problem_->getFeasTol())) {
534          // Can occur: two objects for same var, first scanned object
535          // has both branches feasible and is saved as bestObjectIndex_,
536          // second scanned object fixes the variable
537          returnCode = 4;
538        }
539      }
540      retval = returnCode;
541    }
542    else {
543      retval = 1;
544    }
545
546    // COPY of Bonmin ends here //////////////////////////////////////////////
547
548
549#ifdef TRACE_STRONG2
550    if(info->depth_ > minDepthPrint_) {
551      if(pv > -1) {
552        printf("CCS: end: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
553               pv, solver->getColSolution()[pv], solver->getColLower()[pv], solver->getColUpper()[pv]);
554        printf("CCS: info: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
555               pv, info->solution_[pv], info->lower_[pv], info->upper_[pv]);
556        printf("CCS: problem: lb: %10.4f  ub: %10.4f\n",
557               problem_->Lb(pv), problem_->Ub(pv));
558      }
559      printf("retval: %d\n", retval);
560    }
561#endif
562 
563    problem_ -> domain () -> pop ();
564
565    return retval;
566  }
567/*******************************************************************/
568#else /* not FM_MOD */
569  int CouenneChooseStrong::chooseVariable (OsiSolverInterface * solver,
570                                           OsiBranchingInformation *info,
571                                           bool fixVariables) {
572
573    /// Note: had to copy code from
574    /// BonChooseVariable::chooseVariable() in order to test product
575    /// thing
576
577    problem_ -> domain () -> push
578      (problem_ -> nVars (),
579       info -> solution_,
580       info -> lower_,
581       info -> upper_);
582
583    int retval;
584
585#ifdef TRACE_STRONG2
586    int pv = -1;
587    if(info->depth_ > minDepthPrint_) {
588      if(pv > -1) {
589        printf("CCS: beg: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
590               pv, solver->getColSolution()[pv], solver->getColLower()[pv], solver->getColUpper()[pv]);
591        printf("CCS: info: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
592               pv, info->solution_[pv], info->lower_[pv], info->upper_[pv]);
593        printf("CCS: problem: lb: %10.4f  ub: %10.4f\n",
594               problem_->Lb(pv), problem_->Ub(pv));
595      }
596    }
597#endif
598
599    //int retval = BonChooseVariable::chooseVariable (solver, info, fixVariables);
600
601    // COPY of Bonmin starts here ////////////////////////////////////////////
602
603    // We assume here that chooseVariable is called once at the very
604    // beginning with fixVariables set to true.  This is then the root
605    // node.
606
607    bool isRoot = isRootNode(info);
608    int numberStrong = numberStrong_;
609    if (isRoot) {
610      numberStrong = CoinMax(numberStrong_, numberStrongRoot_);
611    }
612    if (numberUnsatisfied_) {
613      const double* upTotalChange = pseudoCosts_.upTotalChange();
614      const double* downTotalChange = pseudoCosts_.downTotalChange();
615      const int* upNumber = pseudoCosts_.upNumber();
616      const int* downNumber = pseudoCosts_.downNumber();
617      int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
618      int numberLeft = CoinMin(numberStrong -numberStrongDone_,numberUnsatisfied_);
619      results_.clear();
620      int returnCode=0;
621      bestObjectIndex_ = -1;
622      bestWhichWay_ = -1;
623      firstForcedObjectIndex_ = -1;
624      firstForcedWhichWay_ =-1;
625      double bestTrusted=-COIN_DBL_MAX;
626      for (int i=0;i<numberLeft;i++) {
627        int iObject = list_[i];
628        if (numberBeforeTrusted == 0||
629            i < minNumberStrongBranch_ ||
630            (
631              only_pseudo_when_trusted_ && number_not_trusted_>0 ) ||
632              ( !isRoot && (upNumber[iObject]<numberBeforeTrusted ||
633                          downNumber[iObject]<numberBeforeTrusted ))||
634              ( isRoot && (!upNumber[iObject] && !downNumber[iObject])) ) {
635
636#ifdef TRACE_STRONG
637          if(info->depth_ > minDepthPrint_) {
638            printf("Push object %d for strong branch\n", iObject);
639          }
640#endif
641          results_.push_back(Bonmin::HotInfo(solver, info,
642                                             solver->objects(), iObject));
643        }
644        else {
645
646#ifdef TRACE_STRONG
647          if(info->depth_ > minDepthPrint_) {
648            printf("Use pseudo cost for object %d\n", iObject);
649          }
650#endif
651
652          const OsiObject * obj = solver->object(iObject);
653          double
654            upEstimate       = (upTotalChange[iObject]*obj->upEstimate())/upNumber[iObject],
655            downEstimate     = (downTotalChange[iObject]*obj->downEstimate())/downNumber[iObject],
656            MAXMIN_CRITERION = maxminCrit (info),
657            minVal, maxVal, value;
658
659          if (downEstimate < upEstimate) {
660            minVal = downEstimate;
661            maxVal = upEstimate;
662          } else {
663            minVal = upEstimate;
664            maxVal = downEstimate;
665          }
666
667          value = 
668            estimateProduct_ ? 
669            ((estProdEps + minVal) * maxVal) :
670            (       MAXMIN_CRITERION  * minVal + 
671             (1.0 - MAXMIN_CRITERION) * maxVal);
672
673          if (value > bestTrusted) {
674            bestObjectIndex_=iObject;
675            bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
676            bestTrusted = value;
677          }
678        }
679      }
680
681      int numberFixed=0;
682
683      if (results_.size() > 0) {
684
685        // do strong branching //////////////////////////////////////////////////
686        returnCode = doStrongBranching (solver, info, results_.size(), 1);
687
688        if (bb_log_level_>=3) {
689          const char* stat_msg[] = {"NOTDON", "FEAS", "INFEAS", "NOFINI"};
690          message(SB_HEADER)<<CoinMessageEol;
691          for (unsigned int i = 0; i< results_.size(); i++) {
692            double up_change = results_[i].upChange();
693            double down_change = results_[i].downChange();
694            int up_status = results_[i].upStatus();
695            int down_status = results_[i].downStatus();
696            message(SB_RES)<<(int) i<<stat_msg[down_status+1]<<down_change
697            <<stat_msg[up_status+1]<< up_change<< CoinMessageEol;
698          }
699        }
700
701        if (returnCode >= 0 && 
702            returnCode <= 2) {
703
704          if (returnCode)
705            returnCode = (bestObjectIndex_>=0) ? 3 : 4;
706
707          for (unsigned int i=0;i < results_.size (); i++) {
708
709            if((results_[i].upStatus() < 0) || (results_[i].downStatus() < 0))
710              continue;
711
712
713            if((results_[i].upStatus() < 0) || (results_[i].downStatus() < 0)) {
714              continue;
715            }
716
717            int iObject = results_[i].whichObject();
718           
719            ///////////////////////////////////////////////////////////////////
720
721            double upEstimate;
722
723            if (results_[i].upStatus()!=1) {
724              assert (results_[i].upStatus()>=0);
725              upEstimate = results_[i].upChange();
726            }
727            else {
728              // infeasible - just say expensive
729              if (info->cutoff_<1.0e50)
730                upEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
731              else
732                upEstimate = 2.0*fabs(info->objectiveValue_);
733              if (firstForcedObjectIndex_ <0) {
734                // first fixed variable
735                firstForcedObjectIndex_ = iObject;
736                firstForcedWhichWay_ =0;
737              }
738
739              numberFixed++;
740
741              if (fixVariables) {
742
743                const OsiObject * obj = solver->objects()[iObject];
744                OsiBranchingObject * branch = obj -> createBranch (solver, info, 0);
745                branch -> branch (solver);
746                delete branch;
747              }
748            }
749
750            ///////////////////////////////////////////////////////////////////
751
752            double downEstimate;
753
754            if (results_[i].downStatus()!=1) {
755              assert (results_[i].downStatus()>=0);
756              downEstimate = results_[i].downChange();
757            }
758            else {
759              // infeasible - just say expensive
760              if (info->cutoff_<1.0e50)
761                downEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
762              else
763                downEstimate = 2.0*fabs(info->objectiveValue_);
764              if (firstForcedObjectIndex_ <0) {
765                firstForcedObjectIndex_ = iObject;
766                firstForcedWhichWay_ =1;
767              }
768              numberFixed++;
769              if (fixVariables) {
770
771                const OsiObject * obj = solver->objects()[iObject];
772                OsiBranchingObject * branch = obj -> createBranch (solver, info, 1);
773                branch -> branch (solver);
774                delete branch;
775              }
776            }
777
778            double
779              MAXMIN_CRITERION = maxminCrit (info),
780              minVal, maxVal, value;
781
782            if (downEstimate < upEstimate) {
783              minVal = downEstimate;
784              maxVal = upEstimate;
785            } else {
786              minVal = upEstimate;
787              maxVal = downEstimate;
788            }
789
790            value = 
791              estimateProduct_ ? 
792              ((estProdEps + minVal) * maxVal) :
793              (       MAXMIN_CRITERION  * minVal + 
794               (1.0 - MAXMIN_CRITERION) * maxVal);
795
796            if (value>bestTrusted) {
797              bestTrusted = value;
798              bestObjectIndex_ = iObject;
799              bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
800              // but override if there is a preferred way
801              const OsiObject * obj = solver->object(iObject);
802              if (obj->preferredWay()>=0&&obj->infeasibility())
803                bestWhichWay_ = obj->preferredWay();
804              if (returnCode)
805                returnCode = 2;
806            }
807          }
808        }
809        else if (returnCode==3) {
810          // max time - just choose one
811          bestObjectIndex_ = list_[0];
812          bestWhichWay_ = 0;
813          returnCode=0;
814        }
815      }
816      else {
817        bestObjectIndex_=list_[0];
818      }
819   
820      if(bestObjectIndex_ >=0) {
821        OsiObject * obj = solver->objects()[bestObjectIndex_];
822        obj->setWhichWay(bestWhichWay_);
823        message(BRANCH_VAR)<<bestObjectIndex_<< bestWhichWay_ <<CoinMessageEol;
824        CouenneObject *co =  dynamic_cast <CouenneObject *>(solver->objects()[bestObjectIndex_]);
825
826        int objectInd = -1;
827        if(co) {
828          objectInd = co->Reference()->Index();
829        }
830        else {
831          objectInd = obj->columnNumber();
832        }
833
834#ifdef TRACE_STRONG
835        if(info->depth_ > minDepthPrint_) {
836          if(objectInd >= 0) {
837            double vUb = solver->getColUpper()[objectInd];                     
838            double vLb = solver->getColLower()[objectInd];                     
839            double vSolI = info->solution_[objectInd];                         
840            double vSolS = solver->getColSolution()[objectInd];                 
841            printf("Branch on object %d (var: %d): solInfo: %10.4f  SolSolver: \
842%10.4f low: %10.4f  up: %10.4f\n",                                             
843                   bestObjectIndex_, objectInd, vSolI, vSolS, vLb, vUb);       
844          }
845          else {
846            printf("Branch on object %d (var: -1)\n", bestObjectIndex_);           
847          }
848        }
849#endif
850
851      }
852
853      message(CHOSEN_VAR)<<bestObjectIndex_<<CoinMessageEol;
854
855      if (numberFixed==numberUnsatisfied_&&numberFixed)
856        returnCode=4;
857      retval = returnCode;
858    }
859    else {
860      retval = 1;
861    }
862
863    // COPY of Bonmin ends here //////////////////////////////////////////////
864
865#ifdef TRACE_STRONG2
866    if(info->depth_ > minDepthPrint_) {
867      if(pv > -1) {
868        printf("CCS: end: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
869               pv, solver->getColSolution()[pv], solver->getColLower()[pv], solver->getColUpper()[pv]);
870        printf("CCS: info: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
871               pv, info->solution_[pv], info->lower_[pv], info->upper_[pv]);
872        printf("CCS: problem: lb: %10.4f  ub: %10.4f\n",
873               problem_->Lb(pv), problem_->Ub(pv));
874      }                                                   
875      printf("retval: %d\n", retval);
876    }
877#endif
878    problem_ -> domain () -> pop ();
879
880    return retval;
881    }
882#endif /* not FM_MOD */
883
884  // Sets up strong list and clears all if initialize is true.
885  // Returns number of infeasibilities.
886  int CouenneChooseStrong::setupList (OsiBranchingInformation *info, bool initialize) {
887    static bool warned = false;
888
889    initialize = true; // to avoid failed assert in BonChooseVariable::setupList()
890
891    problem_ -> domain () -> push
892      (problem_ -> nVars (),
893       info -> solution_, 
894       info -> lower_, 
895       info -> upper_); // have to alloc+copy
896
897#ifdef COIN_HAS_NTY
898
899    if (problem_ -> orbitalBranching ()) {
900
901      problem_ -> ChangeBounds (info -> lower_,
902                                info -> upper_, 
903                                problem_ -> nVars ());
904   
905      problem_ -> Compute_Symmetry();
906    }
907
908#endif
909
910
911    jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, 
912                      "----------------- (strong) setup list\n");
913
914    if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
915      for (int i=0; i<problem_ -> domain () -> current () -> Dimension (); i++)
916        printf ("%4d %20.4g [%20.4g %20.4g]\n", i,
917                info -> solution_ [i], info -> lower_ [i], info -> upper_ [i]);
918    }
919
920    OsiObject ** object = info->solver_->objects();
921    int numberObjects = info->solver_->numberObjects();
922
923#ifdef TRACE_STRONG
924    if(info->depth_ > minDepthPrint_) {
925      printObjViol(info);
926    }
927#endif
928
929    int retval = gutsOfSetupList(info, initialize);
930
931    if(retval == 0) { // No branching is possible
932      double ckObj = info->objectiveValue_;
933
934#ifdef FM_CHECKNLP2
935      if(!(problem_->checkNLP2(info->solution_, 
936                               info->objectiveValue_, true, // care about obj
937                               false, // do not stop at first viol
938                               true, // checkAll
939                               problem_->getFeasTol()))) {
940                                // false for NOT stopping at first violation
941        if (!warned) {
942          printf("CouenneChooseStrong::setupList(): ### WARNING: checkNLP2() returns infeasible, no branching object selected\n");
943          warned = true;
944        }
945      }
946#else /* not FM_CHECKNLP2 */
947      if(!(problem_->checkNLP(info->solution_, ckObj, true))) {
948        if (!warned) {
949          printf("CouenneChooseStrong::setupList(): ### WARNING: checkNLP() returns infeasible, no branching object selected\n");
950          warned = true;
951        }
952      }
953#endif /* not FM_CHECKNLP2 */
954       
955#ifdef FM_TRACE_OPTSOL
956#ifdef FM_CHECKNLP2
957      problem_->getRecordBestSol()->update();
958#else /* not FM_CHECKNLP2 */
959      problem_->getRecordBestSol()->update(info->solution_, problem_->nVars(),
960                                           ckObj, problem_->getFeasTol());
961#endif /* not FM_CHECKNLP2 */
962#endif
963    }
964
965#ifdef TRACE_STRONG
966    if(info->depth_ > minDepthPrint_) {
967      printf("Strong list: (obj_ind var_ind priority useful)\n");
968      printf("numberStrong: %d  numberStrongRoot: %d  retval: %d\n", 
969             numberStrong_, numberStrongRoot_, retval);
970      for(int i=0; i<retval; i++) {
971        CouenneObject *co =  dynamic_cast <CouenneObject *>(object[list_[i]]);
972        int objectInd = -1;
973        if(co) {
974          objectInd = co->Reference()->Index();
975        }
976        else {
977          objectInd = object[list_[i]]->columnNumber();
978        }
979        printf(" (%d %d %d %6.4f)", list_[i], objectInd, 
980               object[list_[i]]->priority(), useful_[i]);
981      }
982      printf("\n");
983    }
984#endif
985 
986
987    jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, 
988                      "----------------- (strong) setup list done - %d infeasibilities\n", retval);
989
990    problem_ -> domain () -> pop ();
991    return retval;
992  }
993
994/****************************************************************************/
995// Copied from BonChooseVariable.cpp and modified slightly
996  int CouenneChooseStrong::gutsOfSetupList(OsiBranchingInformation *info, 
997                                      bool initialize)
998  {
999    if (numberBeforeTrustedList_ < 0) {
1000      number_not_trusted_ = 1;
1001      printf("CouenneChooseStrong::gutsOfSetupList(): Did not think we were using this; Please double check ...\n");
1002      exit(1);
1003      return OsiChooseVariable::setupList(info, initialize);
1004    }
1005    if (initialize) {
1006      status_=-2;
1007      delete [] goodSolution_;
1008      bestObjectIndex_=-1;
1009      numberStrongDone_=0;
1010      numberStrongIterations_ = 0;
1011      numberStrongFixed_ = 0;
1012      goodSolution_ = NULL;
1013      goodObjectiveValue_ = COIN_DBL_MAX;
1014      number_not_trusted_=0;
1015    }
1016    else {
1017      throw CoinError(CNAME,"setupList","Should not be called with initialize==false");
1018    }
1019    numberOnList_=0;
1020    numberUnsatisfied_=0;
1021    int numberObjects = solver_->numberObjects();
1022    assert (numberObjects);
1023    if (numberObjects>pseudoCosts_.numberObjects()) {
1024      //std::cout<<"Number objects "<<numberObjects<<std::endl;
1025      //AW : How could that ever happen?
1026      //PB : It happens for instance when SOS constraints are added. They are added after the creation of this.
1027      //   assert(false && "Right now, all old content is deleted!");
1028      // redo useful arrays
1029      int saveNumberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
1030      pseudoCosts_.initialize(numberObjects);
1031      pseudoCosts_.setNumberBeforeTrusted(saveNumberBeforeTrusted);
1032    }
1033    double check = -COIN_DBL_MAX;
1034    int checkIndex=0;
1035    int bestPriority=COIN_INT_MAX;
1036    int putOther = numberObjects;
1037    int i;
1038
1039#ifdef FM_SORT_STRONG
1040    int numStr = numberStrong_;
1041    if(isRootNode(info)) {
1042      numStr = numberStrongRoot_;
1043    }
1044    int maximumStrong = CoinMin(numStr, numberObjects) ;
1045    int lastPrio = problem_->getLastPrioSort();
1046    int card_vPriority = 0;
1047    int posEnd_vPriority = numberObjects;
1048    double *vPriority = new double[numberObjects];
1049#else /* not FM_SORT_STRONG */
1050    int maximumStrong = CoinMin(CoinMax(numberStrong_,numberStrongRoot_),
1051        numberObjects) ;
1052    for (i=0;i<numberObjects;i++) {
1053      list_[i]=-1;
1054      useful_[i]=0.0;
1055    }
1056    // We make a second list for most fractional variables
1057    int* list2 = NULL;
1058    double* useful2 = NULL;
1059    double check2 = -COIN_DBL_MAX;
1060    int checkIndex2=0;
1061    int max_most_fra = setup_pseudo_frac_ > 0. ? (int)floor(setup_pseudo_frac_*(double)maximumStrong): 0;
1062    if (setup_pseudo_frac_ > 0.) {
1063      max_most_fra = CoinMax(1, max_most_fra);
1064    }
1065    if (max_most_fra) {
1066      list2 = new int[max_most_fra];
1067      useful2 = new double[max_most_fra];
1068      for (i=0;i<max_most_fra;i++) {
1069        list2[i]=-1;
1070        useful2[i]=0.0;
1071      }
1072    }
1073#endif /* not FM_SORT_STRONG */
1074
1075#ifdef FM_CHECK
1076    const double* upTotalChange = pseudoCosts_.upTotalChange();
1077    const double* downTotalChange = pseudoCosts_.downTotalChange();
1078    int pseudoNum = pseudoCosts_.numberObjects();
1079    for(i=0; i<pseudoNum; i++) {
1080      if(isnan(upTotalChange[i]) || isinf(upTotalChange[i])) {
1081        printf("CouenneChooseStrong::gutsOfSetupList(): upTotalChange[%d]: not a number or infinite\n", i);
1082        exit(1);
1083      }
1084      if(isnan(downTotalChange[i]) || isinf(downTotalChange[i])) {
1085        printf("CouenneChooseStrong::gutsOfSetupList(): downTotalChange[%d]: not a number or infinite\n", i);
1086        exit(1);
1087      }
1088    }
1089#endif
1090
1091    OsiObject ** object = info->solver_->objects();
1092    double upMultiplier, downMultiplier;
1093    computeMultipliers(upMultiplier, downMultiplier);
1094
1095    // Say feasible
1096    bool feasible = true;
1097    const double MAXMIN_CRITERION = maxminCrit(info);
1098
1099    bool firstPass = false; // not important; useful for making two
1100                            // passes, picking different objects
1101    while(numberOnList_ == 0) {
1102      for(i=0;i<numberObjects;i++) {
1103        int way;
1104        double value = object[i]->infeasibility(info, way);
1105        double lbForInfeas = 0.0;
1106        if(value > lbForInfeas) {
1107          numberUnsatisfied_++;
1108          if(value >= 1e50) {
1109            // infeasible
1110            feasible=false;
1111            break;
1112          }
1113          int priorityLevel = object[i]->priority();
1114
1115#ifdef FM_SORT_STRONG
1116          if(priorityLevel > lastPrio) {
1117            posEnd_vPriority--;
1118            vPriority[posEnd_vPriority] = priorityLevel;
1119            list_[posEnd_vPriority] = i;
1120          }
1121          else {
1122            vPriority[card_vPriority] = priorityLevel;
1123            list_[card_vPriority] = i;
1124            card_vPriority++;
1125          }
1126#else /* not FM_SORT_STRONG */
1127          // Better priority? Flush choices.
1128          if(priorityLevel < bestPriority) {
1129            for (int j=maximumStrong-1; j>=0; j--) {
1130              if(list_[j] >= 0) {
1131                int iObject = list_[j];
1132                list_[j]=-1;
1133                useful_[j]=0.0;
1134                list_[--putOther]=iObject;
1135              }
1136            }
1137            maximumStrong = CoinMin(maximumStrong,putOther);
1138            bestPriority = priorityLevel;
1139            check=-COIN_DBL_MAX;
1140            checkIndex=0;
1141            check2=-COIN_DBL_MAX;
1142            checkIndex2=0;
1143            number_not_trusted_=0;
1144            if(max_most_fra > 0) {
1145              for(int j=0; j<max_most_fra; j++) {
1146                list2[j]=-1;
1147                useful2[j]=0.0;
1148              }
1149            }
1150          }
1151          if(priorityLevel == bestPriority) {
1152            // Modify value
1153            double value2;
1154            value = computeUsefulness(MAXMIN_CRITERION,
1155                                      upMultiplier, downMultiplier, value,
1156                                      object[i], i, value2);
1157            if(value > check) {
1158              //add to list
1159              int iObject = list_[checkIndex];
1160              if(iObject >= 0) {
1161                assert (list_[putOther-1]<0);
1162                list_[--putOther]=iObject;  // to end
1163              }
1164              list_[checkIndex]=i;
1165              assert (checkIndex<putOther);
1166              useful_[checkIndex]=value;
1167              // find worst
1168              check=COIN_DBL_MAX;
1169              maximumStrong = CoinMin(maximumStrong,putOther);
1170              for (int j=0; j<maximumStrong; j++) {
1171                if(list_[j]>=0) {
1172                  if (useful_[j]<check) {
1173                    check=useful_[j];
1174                    checkIndex=j;
1175                  }
1176                }
1177                else {
1178                  check=0.0;
1179                  checkIndex = j;
1180                  break;
1181                }
1182              }
1183            }
1184            else {
1185              // to end
1186              assert (list_[putOther-1]<0);
1187              list_[--putOther]=i;
1188              maximumStrong = CoinMin(maximumStrong,putOther);
1189            }
1190            if(max_most_fra > 0 && value2 > check2) {
1191              // add to list of integer infeasibilities
1192              number_not_trusted_++;
1193              list2[checkIndex2]=i;
1194              useful2[checkIndex2]=value2;
1195              // find worst
1196              check2=COIN_DBL_MAX;
1197              for(int j=0; j<max_most_fra; j++) {
1198                if(list2[j] >= 0) {
1199                  if(useful2[j] < check2) {
1200                    check2=useful2[j];
1201                    checkIndex2=j;
1202                  }
1203                }
1204                else {
1205                  check2=0.0;
1206                  checkIndex2 = j;
1207                  break;
1208                }
1209              }
1210            }
1211          }
1212          else {
1213            // worse priority
1214            // to end
1215            assert (list_[putOther-1]<0);
1216            list_[--putOther]=i;
1217            maximumStrong = CoinMin(maximumStrong,putOther);
1218          }
1219#endif /* not FM_SORT_STRONG */
1220        }
1221      }
1222
1223#ifdef FM_SORT_STRONG
1224
1225#ifdef FM_CHECK
1226      if(card_vPriority - posEnd_vPriority + numberObjects != numberUnsatisfied_) {
1227        printf("CouenneChooseStrong::gutsOfSetupList(): ### ERROR: card_vPriority: %d  posEnd_vPriority: %d  numberUnsatisfied: %d numberObjects: %d\n",
1228               card_vPriority, posEnd_vPriority, numberUnsatisfied_, numberObjects);
1229        exit(1);
1230      }
1231#endif
1232
1233      numberOnList_ = 0;
1234      if(feasible) {
1235        int card_smallerThanPrio = card_vPriority;
1236        if(posEnd_vPriority > card_vPriority) {
1237          for(i=posEnd_vPriority; i<numberObjects; i++) {
1238            list_[card_vPriority] = list_[i];
1239            list_[i] = -1;
1240            vPriority[card_vPriority] = vPriority[i]; 
1241            card_vPriority++;
1242          }
1243        }
1244        else {
1245          card_vPriority = numberUnsatisfied_;
1246        }
1247        // correct bounds if card_smallThanPrio >= maximumStrong
1248        int sortFrom = 0;
1249        int sortUpTo = card_smallerThanPrio;
1250        if(card_smallerThanPrio < maximumStrong) {
1251          sortFrom = card_smallerThanPrio;
1252          sortUpTo = card_vPriority;
1253        }
1254        if(card_vPriority > 0) {
1255          numberOnList_ = (card_vPriority < maximumStrong ? card_vPriority : maximumStrong);
1256
1257#ifdef FM_ALWAYS_SORT
1258          bool alwaysSort = true;
1259#else     
1260          bool alwaysSort = false;
1261#endif
1262          if(alwaysSort) {
1263            sortFrom = 0;
1264            sortUpTo = card_vPriority;
1265          }
1266          if((sortUpTo > maximumStrong) || alwaysSort){
1267            // sort list_[card_sortFrom..card_sortUpTo-1] according to priority
1268            CoinSort_2(vPriority + sortFrom, vPriority + sortUpTo, 
1269                       list_ + sortFrom);
1270          }
1271          for(i=0; i<card_vPriority; i++) {
1272            int indObj = list_[i];
1273            double value, value2;
1274            value = computeUsefulness(MAXMIN_CRITERION,
1275                                      upMultiplier, downMultiplier, value,
1276                                      object[indObj], indObj, value2);
1277
1278#ifdef OLD_USEFULLNESS
1279            useful_[i] = -value;
1280#else
1281            if ((sortCrit_ & 1) == 0) {
1282              useful_[i] = -value;
1283            }
1284            else {
1285              useful_[i] = value;
1286            }
1287#endif
1288          }
1289       
1290          if(sortUpTo > maximumStrong) {
1291            // compute from, upto such that
1292            // vPriority[k] == vPriority[maximumStrong] for k in [from..upto-1]
1293            int from = maximumStrong-1, upto = maximumStrong;
1294            int msPrio = vPriority[maximumStrong-1];
1295            problem_->setLastPrioSort(msPrio);
1296            while((from > -1) && (vPriority[from] == msPrio)) {
1297              from--;
1298            }
1299            from++;
1300            while((upto < sortUpTo) && (vPriority[upto] == msPrio)) {
1301              upto++;
1302            }
1303            // sort list[from]..list[upto-1] according to
1304            // useful_[from]..useful_[upto-1]
1305            CoinSort_2(useful_+from, useful_+upto, list_+from);
1306          }
1307        }
1308
1309#ifdef FM_CHECK
1310        // priority of last selected object
1311        double ckPrio = (card_vPriority < numberUnsatisfied_ ?
1312                         vPriority[card_vPriority] : 100000);
1313        double ckUse = (card_vPriority < numberUnsatisfied_ ?
1314                      useful_[card_vPriority] : 100000);
1315        for(i=0; i<card_vPriority; i++) {
1316          int indObj = list_[i];
1317          if(object[indObj]->priority() > ckPrio + 1e-3) {
1318            printf("CouenneChooseStrong::gutsOfSetupList(): ### ERROR: object[%d]->priority(): %d  > ckPrio: %d\n", 
1319                   indObj, object[indObj]->priority(), ckPrio);
1320            exit(1);
1321          }
1322          if(fabs(object[indObj]->priority() - ckPrio) < 1e-3) {
1323            if(useful_[i] > ckUse + 1e-3) {
1324              printf("CouenneChooseStrong::gutsOfSetupList(): ### ERROR: object[%d]->useful: %f  > ckUse: %d\n", 
1325                   indObj, useful_[i], ckUse);
1326              exit(1);
1327            }
1328          }
1329        }
1330        for(i=card_vPriority; i<numberUnsatisfied_; i++) {
1331          int indObj = list_[i];
1332          if(object[indObj]->priority() < ckPrio - 1e-3) {
1333            printf("CouenneChooseStrong::gutsOfSetupList(): ### ERROR: object[%d]->priority(): %d  < ckPrio: %d\n", 
1334                   indObj, object[indObj]->priority(), ckPrio);
1335            exit(1);
1336          }
1337          if(fabs(object[indObj]->priority() - ckPrio) < 1e-3) {
1338            if(useful_[i] < ckUse - 1e-3) {
1339              printf("CouenneChooseStrong::gutsOfSetupList(): ### ERROR: object[%d]->useful: %f  < ckUse: %d\n", 
1340                   indObj, useful_[i], ckUse);
1341              exit(1);
1342            }
1343          }
1344        }
1345#endif
1346       
1347
1348      }
1349      else {
1350        numberUnsatisfied_ = -1;
1351      }
1352#else /* not FM_SORT_STRONG */
1353      // Get list
1354      numberOnList_=0;
1355      if (feasible) {
1356        maximumStrong = CoinMin(maximumStrong,putOther);
1357        for (i=0;i<maximumStrong;i++) {
1358          if (list_[i]>=0) {
1359#ifdef OLD_USEFULLNESS
1360            list_[numberOnList_]=list_[i];
1361            useful_[numberOnList_++]=-useful_[i];
1362           
1363#else
1364            list_[numberOnList_]=list_[i];
1365            if ((sortCrit_ & 1) == 0) {
1366              useful_[numberOnList_++]=-useful_[i];
1367            }
1368            else useful_[numberOnList_++] = useful_[i];
1369#endif
1370            message(CANDIDATE_LIST2)<<numberOnList_-1
1371                                    <<list_[numberOnList_-1]<<numberOnList_-1<<useful_[numberOnList_-1]
1372                                    <<CoinMessageEol;
1373          }
1374        }
1375        if (numberOnList_) {
1376          int tmp_on_list = 0;
1377          if (max_most_fra > 0 && numberOnList_ >= maximumStrong) {
1378            // If we want to force non-trusted in the list, give them huge
1379            // weight here
1380            number_not_trusted_=0;
1381            for (i=0;i<max_most_fra;i++) {
1382              if (list2[i]>=0) {
1383                list2[number_not_trusted_] = list2[i];
1384                useful2[number_not_trusted_++] = useful2[i];
1385                message(CANDIDATE_LIST3)<<number_not_trusted_-1
1386                                        <<list2[number_not_trusted_-1]<<number_not_trusted_-1
1387                                        <<useful2[number_not_trusted_-1]<<CoinMessageEol;
1388              }
1389            }
1390            if (number_not_trusted_) {
1391              CoinSort_2(list_,list_+numberOnList_,useful_);
1392              CoinSort_2(list2,list2+number_not_trusted_,useful2);
1393              int i1=0;
1394              int i2=0;
1395              for (i=0; i<numberObjects; i++) {
1396                bool found1 = (list_[i1]==i);
1397                bool found2 = (list2[i2]==i);
1398                if (found1 && found2) {
1399                  useful_[i1] = -1e150*(1.+useful2[i2]);
1400                  list2[i2] = -1;
1401                }
1402                if (found1) i1++;
1403                if (found2) i2++;
1404                if (i2==max_most_fra) break;
1405              }
1406              for (i=0; i<number_not_trusted_; i++) {
1407                if (list2[i] >= 0) {
1408                  list_[numberOnList_+tmp_on_list] = list2[i];
1409                  useful_[numberOnList_+tmp_on_list] = -1e150*(1.+useful2[i]);
1410                  tmp_on_list++;
1411                }
1412              }
1413            }
1414          }
1415          // Sort
1416          CoinSort_2(useful_,useful_+numberOnList_+tmp_on_list,list_);
1417          // move others
1418          i = numberOnList_;
1419          for (;putOther<numberObjects;putOther++)
1420            list_[i++]=list_[putOther];
1421          assert (i==numberUnsatisfied_);
1422          if (!CoinMax(numberStrong_,numberStrongRoot_))
1423            numberOnList_=0;
1424        }
1425      }
1426      else {
1427        // not feasible
1428        numberUnsatisfied_=-1;
1429      }
1430#endif /* not FM_SORT_STRONG */
1431
1432      if(!firstPass) {
1433        break;
1434      }
1435      firstPass = false;
1436    } /* while(numberOnList_ == 0) */
1437
1438#ifdef TRACE_STRONG
1439      printf("numberStrong_: %d   maximumStrong: %d\n", 
1440             numberStrong_, maximumStrong);
1441#endif
1442
1443#ifdef FM_SORT_STRONG
1444      delete [] vPriority;
1445#else  /* not FM_SORT_STRONG */
1446    delete [] list2;
1447    delete [] useful2;
1448#endif /* not FM_SORT_STRONG */
1449
1450    // Get rid of any shadow prices info
1451    info->defaultDual_ = -1.0; // switch off
1452    delete [] info->usefulRegion_;
1453    delete [] info->indexRegion_;
1454
1455    int way;
1456    if (bb_log_level_>3) {
1457      //for (int i=0; i<Min(numberUnsatisfied_,numberStrong_); i++)
1458      for (int i=0; i<numberOnList_; i++){
1459        message(CANDIDATE_LIST)<<i<< list_[i]<< i<< useful_[i]
1460        <<object[list_[i]]->infeasibility(info,way)
1461        <<CoinMessageEol;
1462      }
1463    }
1464    // return -1 if infeasible to differentiate with numberOnList_==0
1465    // when feasible
1466    if(numberUnsatisfied_ == -1) {
1467      return(-1);
1468    }
1469    return numberOnList_;
1470  }
1471
1472/****************************************************************************/
1473  /// Add list of options to be read from file ////////////////////////////////////////
1474  void CouenneChooseStrong::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
1475
1476    roptions -> AddStringOption6
1477      ("pseudocost_mult",
1478       "Multipliers of pseudocosts for estimating and update estimation of bound",
1479       "interval_br_rev",
1480
1481       "infeasibility", "infeasibility returned by object",
1482
1483       "projectDist",   "distance between current LP point and resulting branches' LP points",
1484
1485       "interval_lp",   "width of the interval between bound and current lp point",
1486       "interval_lp_rev",   "similar to interval_lp, reversed",
1487
1488       "interval_br",   "width of the interval between bound and branching point",
1489       "interval_br_rev",   "similar to interval_br, reversed");
1490
1491    roptions -> AddStringOption2
1492      ("pseudocost_mult_lp",
1493       "Use distance between LP points to update multipliers of pseudocosts " 
1494       "after simulating branching",
1495       "no",
1496       "yes", "",
1497       "no",  "");
1498
1499    roptions -> AddStringOption2
1500      ("estimate_select",
1501       "How the min/max estimates of the subproblems' bounds are used in strong branching",
1502       "normal",
1503       "normal",   "as usual in literature",
1504       "product",  "use their product");
1505
1506    roptions -> AddStringOption2
1507      ("trust_strong",
1508       "Fathom strong branching LPs when their bound is above the cutoff",
1509       "yes",
1510       "yes", "",
1511       "no",  "");
1512  }
1513
1514
1515  // Returns true if solution looks feasible against given objects
1516  bool CouenneChooseStrong::feasibleSolution (const OsiBranchingInformation * info,
1517                                              const double * solution,
1518                                              int numberObjects,
1519                                              const OsiObject ** objects) {
1520
1521    double obj = solution [problem_ -> Obj (0) -> Body () -> Index ()];
1522
1523#ifdef FM_CHECKNLP2
1524    bool isFeas = problem_->checkNLP2(solution, 0, false, true, true, 
1525                                      problem_->getFeasTol());
1526    return isFeas;
1527#else
1528    return problem_ -> checkNLP (solution, obj);
1529#endif
1530  }
1531/****************************************************************************/
1532  void CouenneChooseStrong::printObjViol(OsiBranchingInformation *info) {
1533
1534    OsiObject ** object = info->solver_->objects();
1535    int numberObjects = info->solver_->numberObjects();
1536
1537    printf("CouenneChooseStrong::printObjViol(): Object violations: (obj_ind  var_ind  violation)");
1538    double maxViol = 0;
1539    double minPosViol = 1e50;
1540    for(int i=0; i<numberObjects; i++) {
1541      CouenneObject *co =  dynamic_cast <CouenneObject *>(object[i]);
1542      int indVar = -1;
1543      if(co) {
1544        indVar = co->Reference()->Index();
1545      }
1546      else {
1547        indVar = object[i]->columnNumber();
1548      }
1549      int way;
1550      double value = object[i]->infeasibility(info,way);
1551      maxViol = (value > maxViol ? value : maxViol);
1552      if(value > 0.0) {
1553        printf("(%d %d %f)", i, indVar, value);
1554        minPosViol = (value < minPosViol ? value : minPosViol);
1555      }
1556    }
1557    printf("\nmaxViol: %g  minPosViol: %g\n", maxViol, minPosViol);
1558
1559  } /* printObjViol */
1560
1561
1562//}
Note: See TracBrowser for help on using the repository browser.