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

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

Fix bug in selection of strong branching when requesting non-full strong-branching; add END_OF_HEADER in output file fres.xxx

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