Ignore:
Timestamp:
Jul 21, 2013 5:05:45 AM (6 years ago)
Author:
forrest
Message:

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcModel.cpp

    r1924 r1943  
    24532453        }
    24542454    }
     2455#ifdef SWITCH_VARIABLES
     2456    // see if any switching variables
     2457    if (numberIntegers_<solver_->getNumCols())
     2458      findSwitching();
     2459#endif
    24552460    /*
    24562461      Run heuristics at the root. This is the only opportunity to run FPump; it
     
    34063411    numberFixedAtRoot_ = 0;
    34073412    numberFixedNow_ = 0;
     3413    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
     3414#ifdef COIN_HAS_CLP
     3415      OsiClpSolverInterface * clpSolver
     3416        = dynamic_cast<OsiClpSolverInterface *> (solver_);
     3417      if (clpSolver) {
     3418        if (getCutoff()>1.0e20) {
     3419          printf("Zapping costs\n");
     3420          int numberColumns=solver_->getNumCols();
     3421          double * zeroCost = new double [numberColumns];
     3422          // could make small random
     3423          memset(zeroCost,0,numberColumns*sizeof(double));
     3424          solver_->setObjective(zeroCost);
     3425          double objValue = solver_->getObjValue();
     3426          solver_->setDblParam(OsiObjOffset,-objValue);
     3427          clpSolver->getModelPtr()->setObjectiveValue(objValue);
     3428          delete [] zeroCost;
     3429        } else {
     3430          moreSpecialOptions2_ &= ~2;
     3431        }
     3432      } else {
     3433#endif
     3434          moreSpecialOptions2_ &= ~2;
     3435#ifdef COIN_HAS_CLP
     3436      }
     3437#endif
     3438    }
    34083439    int numberIterationsAtContinuous = numberIterations_;
    34093440    //solverCharacteristics_->setSolver(solver_);
     
    38113842    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
    38123843    // Save address of root node as we don't want to delete it
    3813     // initialize for print out
    3814     int lastDepth = 0;
    3815     int lastUnsatisfied = 0;
    3816     if (newNode)
    3817         lastUnsatisfied = newNode->numberUnsatisfied();
    38183844    /*
    38193845      The common case is that the lp relaxation is feasible but doesn't satisfy
     
    43564382                messageHandler()->message(CBC_STATUS2, messages())
    43574383                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
    4358                 << lastDepth << lastUnsatisfied << numberIterations_
     4384                << tree_->lastDepth() << tree_->lastUnsatisfied()
     4385                << tree_->lastObjective() << numberIterations_
    43594386                << getCurrentSeconds()
    43604387                << CoinMessageEol ;
     
    43624389                messageHandler()->message(CBC_STATUS2, messages())
    43634390                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
    4364                 << lastDepth << lastUnsatisfied << numberIterations_
     4391                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
    43654392                << getCurrentSeconds()
    43664393                << CoinMessageEol ;
     
    43684395                messageHandler()->message(CBC_STATUS3, messages())
    43694396                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
    4370                 << lastDepth << lastUnsatisfied << numberIterations_ << numberExtraIterations_
     4397                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
    43714398                << getCurrentSeconds()
    43724399                << CoinMessageEol ;
     
    44464473        } else {
    44474474            // Deterministic parallel
    4448             if (tree_->size() < CoinMax(numberThreads_, 8) && !goneParallel) {
     4475          if ((tree_->size() < CoinMax(numberThreads_, 8)||
     4476               hotstartSolution_) && !goneParallel) {
    44494477                node = tree_->bestNode(cutoff) ;
    44504478                // Possible one on tree worse than cutoff
     
    51015129        specialOptions_(0),
    51025130        moreSpecialOptions_(0),
     5131        moreSpecialOptions2_(0),
    51035132        topOfTree_(NULL),
    51045133        subTreeModel_(NULL),
     
    52655294        specialOptions_(0),
    52665295        moreSpecialOptions_(0),
     5296        moreSpecialOptions2_(0),
    52675297        topOfTree_(NULL),
    52685298        subTreeModel_(NULL),
     
    55535583        specialOptions_(rhs.specialOptions_),
    55545584        moreSpecialOptions_(rhs.moreSpecialOptions_),
     5585        moreSpecialOptions2_(rhs.moreSpecialOptions2_),
    55555586        topOfTree_(NULL),
    55565587        subTreeModel_(rhs.subTreeModel_),
     
    59075938        specialOptions_ = rhs.specialOptions_;
    59085939        moreSpecialOptions_ = rhs.moreSpecialOptions_;
     5940        moreSpecialOptions2_ = rhs.moreSpecialOptions2_;
    59095941        subTreeModel_ = rhs.subTreeModel_;
    59105942        heuristicModel_ = NULL;
     
    63636395    specialOptions_ = rhs.specialOptions_;
    63646396    moreSpecialOptions_ = rhs.moreSpecialOptions_;
     6397    moreSpecialOptions2_ = rhs.moreSpecialOptions2_;
    63656398    numberStrong_ = rhs.numberStrong_;
    63666399    numberBeforeTrust_ = rhs.numberBeforeTrust_;
     
    72607293    }
    72617294    numberDJFixed_ += numberFixed - numberTightened;
     7295#ifdef SWITCH_VARIABLES
     7296    if (numberFixed)
     7297      fixAssociated(NULL,0);
     7298#endif
    72627299    return numberFixed;
    72637300}
     
    1099411031            // treat as if will cost what it says up
    1099511032            double upCost = costValue;
     11033#ifndef BRANCH_BREAKEVEN
     11034#define BRANCH_BREAKEVEN 0.3
     11035#else
     11036            preferredWay=1;
     11037#endif
    1099611038            // and balance at breakeven of 0.3
    10997             double downCost = (0.7 * upCost) / 0.3;
     11039            double downCost = ((1.0-BRANCH_BREAKEVEN) * upCost) / BRANCH_BREAKEVEN;
    1099811040            if (obj1a) {
    1099911041                upCost = obj1a->upPseudoCost();
     
    1103111073        branchingMethod_ = new CbcBranchDynamicDecision();
    1103211074    }
     11075#ifdef SWITCH_VARIABLES
     11076    // see if any switching variables
     11077    if (numberIntegers_<solver_->getNumCols())
     11078      findSwitching();
     11079#endif
    1103311080    synchronizeNumberBeforeTrust();
    1103411081}
     11082#ifdef SWITCH_VARIABLES
     11083// Convert Dynamic to Switching
     11084int
     11085CbcModel::findSwitching()
     11086{
     11087  if ((moreSpecialOptions2_&1)==0)
     11088    return 0;
     11089  const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
     11090  const int * column = rowCopy->getIndices();
     11091  const int * rowLength = rowCopy->getVectorLengths();
     11092  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     11093  const double * rowLower = solver_->getRowLower();
     11094  const double * rowUpper = solver_->getRowUpper();
     11095  const double * columnLower = solver_->getColLower();
     11096  const double * columnUpper = solver_->getColUpper();
     11097  const double * element = rowCopy->getElements();
     11098  //const double * element = solver_->getMatrixByCol()->getElements();
     11099  const int * row = solver_->getMatrixByCol()->getIndices();
     11100  const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     11101  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     11102  int numberRows = solver_->getNumRows();
     11103  int numberColumns = solver_->getNumCols();
     11104  int * sort = new int[2*numberRows+2+numberColumns];
     11105  int * whichRow = sort+numberRows+1;
     11106  int * marked = whichRow+numberRows+1;
     11107  memset(marked,0,numberColumns*sizeof(int));
     11108  int nnSwitch=0;
     11109  int nnSwitchTotal=0;
     11110  int n2Switch=0;
     11111  double largeRatio1=1000.0;
     11112  double largeRatio2=100.0;
     11113  for (int i=0;i<numberIntegers_;i++) {
     11114    int iColumn = integerVariable_[i];
     11115    if (columnLower[iColumn]||columnUpper[iColumn]!=1.0)
     11116      continue;
     11117    if (!dynamic_cast <CbcSimpleInteger *> (object_[i]))
     11118      continue;
     11119    int nAdd=0;
     11120    bool takeThis=false;
     11121    CoinBigIndex start = columnStart[iColumn];
     11122    CoinBigIndex end = start + columnLength[iColumn];
     11123    for (CoinBigIndex j = start; j < end; j++) {
     11124      int iRow = row[j];
     11125      if (rowLength[iRow]!=2) {
     11126        continue;
     11127      }
     11128      // for now just 0.0 in rhs
     11129      if (!rowLower[iRow]) {
     11130        if (rowUpper[iRow]!=COIN_DBL_MAX)
     11131          continue;
     11132      } else if (rowLower[iRow]!=-COIN_DBL_MAX) {
     11133        continue;
     11134      } else if (rowUpper[iRow]) {
     11135        continue;
     11136      }
     11137      CoinBigIndex k = rowStart[iRow];
     11138      double bValue, cValue;
     11139      int cColumn;
     11140      if (column[k]==iColumn) {
     11141        bValue=element[k];
     11142        cValue=element[k+1];
     11143        cColumn=column[k+1];
     11144      } else {
     11145        bValue=element[k+1];
     11146        cValue=element[k];
     11147        cColumn=column[k];
     11148      }
     11149      if (solver_->isInteger(cColumn))
     11150        continue;
     11151      if (columnLower[cColumn]<0.0)
     11152        continue;
     11153      if (bValue*cValue>0.0)
     11154        continue;
     11155      if (fabs(bValue)>largeRatio1*fabs(cValue))
     11156        takeThis=true;
     11157      // add to list
     11158      whichRow[nAdd]=iRow;
     11159      sort[nAdd++]=cColumn;
     11160    }
     11161    if (nAdd) {
     11162      n2Switch++;
     11163      CoinSort_2(sort,sort+nAdd,whichRow);
     11164      int last=sort[0];
     11165      for (int k=1;k<nAdd;k++) {
     11166        if (sort[k]==last)
     11167          takeThis=true;
     11168        else
     11169          last=sort[k];
     11170      }
     11171      if (takeThis) {
     11172        int last=sort[0];
     11173        marked[last]++;
     11174        for (int k=1;k<nAdd;k++) {
     11175          if (sort[k]!=last) {
     11176            last=sort[k];
     11177            marked[last]++;
     11178          }
     11179        }
     11180        //printf("Column %d has %d other columns\n",iColumn,nAdd);
     11181        sort[nAdd]=COIN_INT_MAX;
     11182        whichRow[nAdd]=COIN_INT_MAX;
     11183        CbcSimpleIntegerDynamicPseudoCost * thisOne =
     11184          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *> (object_[i]);
     11185        if (thisOne) {
     11186          assert(iColumn == thisOne->columnNumber());
     11187          object_[i]=new CbcSwitchingBinary(thisOne,nAdd,sort,whichRow);
     11188          delete thisOne;
     11189        } else {
     11190          CbcSimpleInteger * thisOne =
     11191            dynamic_cast <CbcSimpleInteger *> (object_[i]);
     11192          assert (thisOne);
     11193          assert(iColumn == thisOne->columnNumber());
     11194          CbcSimpleIntegerDynamicPseudoCost tempObj(this,iColumn,0.1);
     11195          object_[i]=new CbcSwitchingBinary(&tempObj,nAdd,sort,whichRow);
     11196          delete thisOne;
     11197        }
     11198      }
     11199    }
     11200    // see if there is an interesting row
     11201    for (CoinBigIndex j = start; j < end; j++) {
     11202      int iRow = row[j];
     11203      // for now just 0.0 in rhs
     11204      if (!rowLower[iRow]) {
     11205        if (rowUpper[iRow]!=COIN_DBL_MAX)
     11206          continue;
     11207      } else if (rowLower[iRow]!=-COIN_DBL_MAX) {
     11208        continue;
     11209      } else if (rowUpper[iRow]) {
     11210        continue;
     11211      }
     11212      int nOther=0;
     11213      double bEl=0.0;
     11214      double cMax=-COIN_DBL_MAX;
     11215      double cMin=COIN_DBL_MAX;
     11216      for (CoinBigIndex k = rowStart[iRow];
     11217           k<rowStart[iRow]+rowLength[iRow];k++) {
     11218        int jColumn = column[k];
     11219        if (jColumn==iColumn) {
     11220          bEl=element[k];
     11221        } else {
     11222          sort[nOther++]=jColumn;
     11223          if (solver_->isInteger(jColumn)) {
     11224            cMin=-1.0;
     11225            cMax=1.0;
     11226            break;
     11227          } else {
     11228            cMax=CoinMax(cMax,element[k]);
     11229            cMin=CoinMin(cMin,element[k]);
     11230            if (columnLower[jColumn]<0.0) {
     11231              cMin=-1.0;
     11232              cMax=1.0;
     11233              break;
     11234            }
     11235          }
     11236        }
     11237      }
     11238      double largestC = CoinMax(fabs(cMin),fabs(cMax));
     11239      if (((cMin>0.0&&bEl<0.0&&!rowUpper[iRow])||
     11240           (cMin<0.0&&bEl>0.0&&!rowLower[iRow]))&&cMin*cMax>0.0&&
     11241          fabs(bEl)>largeRatio2*largestC) {
     11242        // forces to zero
     11243        CbcSwitchingBinary * object =
     11244          dynamic_cast <CbcSwitchingBinary *> (object_[i]);
     11245        if (!object) {
     11246          // create empty one
     11247          CbcSimpleIntegerDynamicPseudoCost * thisOne =
     11248            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *> (object_[i]);
     11249          if (thisOne) {
     11250            assert(iColumn == thisOne->columnNumber());
     11251            object=new CbcSwitchingBinary(thisOne,0,sort,whichRow);
     11252            delete thisOne;
     11253          } else {
     11254            CbcSimpleInteger * thisOne =
     11255              dynamic_cast <CbcSimpleInteger *> (object_[i]);
     11256            assert (thisOne);
     11257            assert(iColumn == thisOne->columnNumber());
     11258            CbcSimpleIntegerDynamicPseudoCost tempObj(this,iColumn,0.1);
     11259            object=new CbcSwitchingBinary(&tempObj,0,sort,whichRow);
     11260            delete thisOne;
     11261          }
     11262          object_[i]=object;
     11263        }
     11264        object->addZeroSwitches(nOther,sort);
     11265        nnSwitch++;
     11266        nnSwitchTotal+=nOther;
     11267      }
     11268    }
     11269  }
     11270  if (n2Switch+nnSwitch) {
     11271    if (handler_->logLevel()>2)
     11272      printf("%d two switch variables - %d multi (total multi %d)\n",
     11273             n2Switch,nnSwitch,nnSwitchTotal);
     11274    memset(whichRow,0,(numberRows+1)*sizeof(int));
     11275    for (int i=0;i<numberColumns;i++) {
     11276      whichRow[marked[i]]++;
     11277    }
     11278    if (handler_->logLevel()>2) {
     11279      for (int i=0;i<numberRows+1;i++) {
     11280        if (whichRow[i])
     11281          printf("%d variables have %d switches\n",whichRow[i],i);
     11282      }
     11283    }
     11284  }
     11285  delete [] sort;
     11286  // say switches exist
     11287  if (n2Switch+nnSwitch)
     11288    moreSpecialOptions2_|=4;
     11289  return n2Switch+nnSwitch;
     11290}
     11291// Fix associated variables
     11292int
     11293CbcModel::fixAssociated(OsiSolverInterface * solver,int cleanBasis)
     11294{
     11295  int nChanged=0;
     11296  if ((moreSpecialOptions2_&4)!=0) {
     11297    int n=-1;
     11298    while (n) {
     11299      n=0;
     11300      for (int i=0;i<numberObjects_;i++) {
     11301        CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11302        if (object) {
     11303          n += object->setAssociatedBounds(solver,cleanBasis);
     11304        }
     11305      }
     11306      nChanged+=n;
     11307    }
     11308  }
     11309  return nChanged;
     11310}
     11311/* Debug associated variables
     11312   printLevel - 1 summary if bad on fixed
     11313                2 summary if bad on satisfied
     11314                3 for individuals
     11315 */
     11316int
     11317CbcModel::checkAssociated(const OsiSolverInterface * solver,
     11318                          const double * solution,int printLevel)
     11319{
     11320  int nBad=0;
     11321  int nBadFixed=0;
     11322  if ((moreSpecialOptions2_&4)!=0) {
     11323    int nAt0=0;
     11324    int nAt1=0;
     11325    int nBetween=0;
     11326    for (int i=0;i<numberObjects_;i++) {
     11327      CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11328      if (object) {
     11329        int state[3];
     11330        nBad += object->checkAssociatedBounds(solver,solution,printLevel,state,
     11331                                              nBadFixed);
     11332        if (state[0]==0)
     11333          nBetween++;
     11334        else if (state[0]==-1)
     11335          nAt0++;
     11336        else
     11337          nAt1++;
     11338      }
     11339    }
     11340    if (handler_->logLevel()>2) {
     11341      if (printLevel>1||(printLevel==1&&nBadFixed)) {
     11342        printf("%d switches, %d at 0, %d at 1, %d between - %d bad values (%d when fixed)\n",
     11343               nBetween+nAt0+nAt1,nAt0,nAt1,nBetween,nBad,nBadFixed);
     11344        if (nBadFixed && printLevel!=3)
     11345          checkAssociated(solver,solution,3);
     11346      }
     11347    }
     11348  }
     11349  return nBad;
     11350}
     11351#endif
    1103511352// Set numberBeforeTrust in all objects
    1103611353void
     
    1139011707        memcpy(saveUpper, getColUpper(), numberColumns*sizeof(double));
    1139111708        memcpy(saveLower, getColLower(), numberColumns*sizeof(double));
     11709        //#define CLP_INVESTIGATE4
     11710#ifdef CLP_INVESTIGATE4
     11711        {
     11712          int nBad=checkAssociated(solver_,solver_->getColSolution(),1);
     11713          if (nBad)
     11714            checkAssociated(solver_,solver_->getColSolution(),3);
     11715          double largestInfeasibility = 0.0;
     11716          double primalTolerance ;
     11717          double offset;
     11718          solver_->getDblParam(OsiObjOffset, offset);
     11719          solver_->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     11720          const double *objective = getObjCoefficients() ;
     11721          const double * rowLower = solver_->getRowLower() ;
     11722          const double * rowUpper = solver_->getRowUpper() ;
     11723          const double * columnLower = solver_->getColLower() ;
     11724          const double * columnUpper = solver_->getColUpper() ;
     11725          int numberRows = solver_->getNumRows() ;
     11726          double *rowActivity = new double[numberRows] ;
     11727          memset(rowActivity, 0, numberRows*sizeof(double)) ;
     11728          double *rowSum = new double[numberRows] ;
     11729          memset(rowSum, 0, numberRows*sizeof(double)) ;
     11730          int * marked = new int [numberColumns];
     11731          for (int i=0;i<numberColumns;i++)
     11732            marked[i]=-1;
     11733          for (int i=0;i<numberIntegers_;i++)
     11734            marked[integerVariable_[i]]=-2;
     11735          if ((moreSpecialOptions2_&4)!=0) {
     11736            for (int i=0;i<numberObjects_;i++) {
     11737              CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11738              if (object) {
     11739                int iColumn = object->columnNumber();
     11740                const int * other = object->otherVariable();
     11741                marked[iColumn]=-3-other[0];
     11742                int n=object->numberOther();
     11743                for (int k=0;k<n;k++)
     11744                  marked[other[k]]=iColumn;
     11745              }
     11746            }
     11747          }
     11748          const double * element = solver_->getMatrixByCol()->getElements();
     11749          const int * row = solver_->getMatrixByCol()->getIndices();
     11750          const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     11751          const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     11752          const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
     11753          const int * column = rowCopy->getIndices();
     11754          const int * rowLength = rowCopy->getVectorLengths();
     11755          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     11756          const double * elementByRow = rowCopy->getElements();
     11757          double objValue=-offset;
     11758          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     11759            double value = solution[iColumn];
     11760            objValue += value*objective[iColumn];
     11761            if (value>columnUpper[iColumn]) {
     11762              if (value-columnUpper[iColumn]>1.0e-8)
     11763                printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     11764              value=columnUpper[iColumn];
     11765            } else if (value<columnLower[iColumn]) {
     11766              if (value-columnLower[iColumn]<-1.0e-8)
     11767                printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     11768              value=columnLower[iColumn];
     11769            }
     11770            if (value) {
     11771              CoinBigIndex start = columnStart[iColumn];
     11772              CoinBigIndex end = start + columnLength[iColumn];
     11773              for (CoinBigIndex j = start; j < end; j++) {
     11774                int iRow = row[j];
     11775                rowActivity[iRow] += value * element[j];
     11776                rowSum[iRow] += fabs(value * element[j]);
     11777              }
     11778            }
     11779          }
     11780          for (int i = 0 ; i < numberRows ; i++) {
     11781#if 0 //def CLP_INVESTIGATE
     11782            double inf;
     11783            inf = rowLower[i] - rowActivity[i];
     11784            if (inf > primalTolerance)
     11785              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11786                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11787            inf = rowActivity[i] - rowUpper[i];
     11788            if (inf > primalTolerance)
     11789              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11790                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11791#endif
     11792            double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     11793                                           rowLower[i]-rowActivity[i]);
     11794            // but allow for errors
     11795            double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     11796            if (infeasibility>largestInfeasibility*factor) {
     11797              largestInfeasibility = infeasibility/factor;
     11798              printf("Ainf of %g on row %d sum %g scaled %g\n",
     11799                     infeasibility,i,rowSum[i],largestInfeasibility);
     11800              if (infeasibility>1.0e10) {
     11801                for (CoinBigIndex j=rowStart[i];
     11802                     j<rowStart[i]+rowLength[i];j++) {
     11803                  printf("col %d element %g marked %d\n",
     11804                         column[j],elementByRow[j],marked[column[j]]);
     11805                }
     11806              }
     11807            }
     11808          }
     11809          delete [] rowActivity ;
     11810          delete [] rowSum;
     11811          delete [] marked;
     11812          if (largestInfeasibility > 10.0*primalTolerance)
     11813            printf("Alargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     11814          else
     11815            printf("Afeasible (%g) - obj %g\n", largestInfeasibility,objValue);
     11816        }
     11817#endif
    1139211818        // point to useful information
    1139311819        OsiBranchingInformation usefulInfo = usefulInformation();
     
    1140211828        for (i = 0; i < numberObjects_; i++)
    1140311829            object_[i]->feasibleRegion(solver_, &usefulInfo);
     11830#ifdef CLP_INVESTIGATE4
     11831        {
     11832          int nBad=checkAssociated(solver_,solver_->getColSolution(),1);
     11833          if (nBad)
     11834            checkAssociated(solver_,solver_->getColSolution(),3);
     11835          double largestInfeasibility = 0.0;
     11836          double primalTolerance ;
     11837          double offset;
     11838          solver_->getDblParam(OsiObjOffset, offset);
     11839          solver_->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     11840          const double *objective = getObjCoefficients() ;
     11841          const double * rowLower = solver_->getRowLower() ;
     11842          const double * rowUpper = solver_->getRowUpper() ;
     11843          const double * columnLower = solver_->getColLower() ;
     11844          const double * columnUpper = solver_->getColUpper() ;
     11845          int numberRows = solver_->getNumRows() ;
     11846          double *rowActivity = new double[numberRows] ;
     11847          memset(rowActivity, 0, numberRows*sizeof(double)) ;
     11848          double *rowSum = new double[numberRows] ;
     11849          memset(rowSum, 0, numberRows*sizeof(double)) ;
     11850          int * marked = new int [numberColumns];
     11851          for (int i=0;i<numberColumns;i++)
     11852            marked[i]=-1;
     11853          for (int i=0;i<numberIntegers_;i++)
     11854            marked[integerVariable_[i]]=-2;
     11855          if ((moreSpecialOptions2_&4)!=0) {
     11856            for (int i=0;i<numberObjects_;i++) {
     11857              CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11858              if (object) {
     11859                int iColumn = object->columnNumber();
     11860                const int * other = object->otherVariable();
     11861                marked[iColumn]=-3-other[0];
     11862                int n=object->numberOther();
     11863                for (int k=0;k<n;k++)
     11864                  marked[other[k]]=iColumn;
     11865              }
     11866            }
     11867          }
     11868          const double * element = solver_->getMatrixByCol()->getElements();
     11869          const int * row = solver_->getMatrixByCol()->getIndices();
     11870          const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     11871          const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     11872          const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
     11873          const int * column = rowCopy->getIndices();
     11874          const int * rowLength = rowCopy->getVectorLengths();
     11875          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     11876          const double * elementByRow = rowCopy->getElements();
     11877          double objValue=-offset;
     11878          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     11879            double value = solution[iColumn];
     11880            objValue += value*objective[iColumn];
     11881            if (value>columnUpper[iColumn]) {
     11882              if (value-columnUpper[iColumn]>1.0e-8)
     11883                printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     11884              value=columnUpper[iColumn];
     11885            } else if (value<columnLower[iColumn]) {
     11886              if (value-columnLower[iColumn]<-1.0e-8)
     11887                printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     11888              value=columnLower[iColumn];
     11889            }
     11890            if (value) {
     11891              CoinBigIndex start = columnStart[iColumn];
     11892              CoinBigIndex end = start + columnLength[iColumn];
     11893              for (CoinBigIndex j = start; j < end; j++) {
     11894                int iRow = row[j];
     11895                rowActivity[iRow] += value * element[j];
     11896                rowSum[iRow] += fabs(value * element[j]);
     11897              }
     11898            }
     11899          }
     11900          for (int i = 0 ; i < numberRows ; i++) {
     11901#if 0 //def CLP_INVESTIGATE
     11902            double inf;
     11903            inf = rowLower[i] - rowActivity[i];
     11904            if (inf > primalTolerance)
     11905              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11906                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11907            inf = rowActivity[i] - rowUpper[i];
     11908            if (inf > primalTolerance)
     11909              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11910                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11911#endif
     11912            double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     11913                                           rowLower[i]-rowActivity[i]);
     11914            // but allow for errors
     11915            double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     11916            if (infeasibility>largestInfeasibility*factor) {
     11917              largestInfeasibility = infeasibility/factor;
     11918              printf("inf of %g on row %d sum %g scaled %g\n",
     11919                     infeasibility,i,rowSum[i],largestInfeasibility);
     11920              if (infeasibility>1.0e10) {
     11921                for (CoinBigIndex j=rowStart[i];
     11922                     j<rowStart[i]+rowLength[i];j++) {
     11923                  printf("col %d element %g marked %d\n",
     11924                         column[j],elementByRow[j],marked[column[j]]);
     11925                }
     11926              }
     11927            }
     11928          }
     11929          delete [] rowActivity ;
     11930          delete [] rowSum;
     11931          delete [] marked;
     11932          if (largestInfeasibility > 10.0*primalTolerance)
     11933            printf("Largest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     11934          else
     11935            printf("Feasible (%g) - obj %g\n", largestInfeasibility,objValue);
     11936        }
     11937#endif
    1140411938        // If relaxed then leave bounds on basic variables
    1140511939        if (fixVariables == -1 && (specialOptions_&16) == 0) {
     
    1142211956        }
    1142311957        // We can switch off check
    11424         if ((specialOptions_&4) == 0) {
     11958        if ((specialOptions_&4) == 0 && (moreSpecialOptions2_&10) != 8) {
    1142511959            if ((specialOptions_&2) == 0 && solverCharacteristics_->warmStart()) {
    1142611960                /*
     
    1144811982            solver_->setHintParam(OsiDoDualInInitial, true, OsiHintTry);
    1144911983            solver_->initialSolve();
     11984#ifdef SWITCH_VARIABLES
     11985            if (solver_->isProvenOptimal()) {
     11986              int nBad=checkAssociated(solver_,solver_->getColSolution(),1);
     11987              if (nBad)
     11988                checkAssociated(solver_,solver_->getColSolution(),3);
     11989            }
     11990#endif
    1145011991#ifdef JJF_ZERO
    1145111992            if (solver_->isProvenOptimal()) {
    11452                 solver_->writeMps("feasible");
     11993              solver_->writeMpsNative("feasible.mps",NULL,NULL,2);
     11994#ifdef COIN_HAS_CLP
     11995              OsiClpSolverInterface * clpSolver
     11996                = dynamic_cast<OsiClpSolverInterface *> (solver_);
     11997              if (clpSolver ) {
     11998                clpSolver->getModelPtr()->writeBasis("feasible.bas",true);
     11999              }
     12000#endif
    1145312001                printf("XXXXXXXXXXXX - saving feasible\n");
    1145412002            }
     
    1162312171                const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
    1162412172                const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
    11625                 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    11626                     double value = solution[iColumn];
     12173                double offset;
     12174                solver_->getDblParam(OsiObjOffset, offset);
     12175                double objValue=-offset;
     12176                const double *objective = getObjCoefficients() ;
     12177                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     12178                  double value = solution[iColumn];
     12179                  objValue += value*objective[iColumn];
    1162712180                    if (value) {
    1162812181                        CoinBigIndex start = columnStart[iColumn];
     
    1164712200                               i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
    1164812201#endif
    11649                     largestInfeasibility = CoinMax(largestInfeasibility,
    11650                                                    rowLower[i] - rowActivity[i]);
    11651                     largestInfeasibility = CoinMax(largestInfeasibility,
    11652                                                    rowActivity[i] - rowUpper[i]);
     12202                    double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     12203                                                   rowLower[i]-rowActivity[i]);
     12204                    // but allow for errors
     12205                    double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     12206                    if (infeasibility>largestInfeasibility*factor) {
     12207                      largestInfeasibility = infeasibility/factor;
     12208                      //printf("inf of %g on row %d sum %g scaled %g\n",
     12209                      //     infeasibility,i,rowSum[i],largestInfeasibility);
     12210                    }
    1165312211                }
    1165412212                delete [] rowActivity ;
    1165512213                delete [] rowSum;
     12214                if (handler_->logLevel()>2) {
     12215                  if (largestInfeasibility > 10.0*primalTolerance)
     12216                    printf("BLargest infeasibility is %g - obj %g (%g)\n", largestInfeasibility,objValue,objectiveValue);
     12217                  else
     12218                    printf("BFeasible (%g) - obj %g %g\n", largestInfeasibility,objValue,objectiveValue);
     12219                }
     12220                //if (fabs(objValue-objectiveValue)>1.0e-7*fabs(objectiveValue)) {
     12221                //printf("Bad obj values\n");
     12222                objectiveValue = objValue;
     12223                //}
    1165612224#ifdef CLP_INVESTIGATE
    1165712225                if (largestInfeasibility > 10.0*primalTolerance)
     
    1186712435        assert(basis != NULL);
    1186812436        objectiveValue = checkSolution(cutoff, solution, fixVariables, objectiveValue);
    11869         if (saveObjectiveValue + 1.0e-3 < objectiveValue) {
     12437        if (cutoff>1.0e40&&objectiveValue<1.0e10)
     12438          saveObjectiveValue = objectiveValue; // take anyway
     12439        if (saveObjectiveValue + 1.0e-3 +1.0e-7*fabs(saveObjectiveValue)
     12440            < objectiveValue) {
    1187012441#if COIN_DEVELOP>1
    1187112442            printf("First try at solution had objective %.16g, rechecked as %.16g\n",
     
    1199612567                    cutoff -= 2.0e-5;
    1199712568            }
     12569            if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
     12570              // put back objective
     12571              solver_->setObjective(continuousSolver_->getObjCoefficients());
     12572              double offset;
     12573              continuousSolver_->getDblParam(OsiObjOffset,offset);
     12574              solver_->setDblParam(OsiObjOffset,offset);
     12575              moreSpecialOptions2_ &= ~2;
     12576            }
    1199812577            // This is not correct - that way cutoff can go up if maximization
    1199912578            //double direction = solver_->getObjSense();
     
    1254513124                    } else {
    1254613125                        // relax a bit
    12547                         value = CoinMin(saveUpper, value + 1.0e-5 * (fabs(saveUpper) + 1));
     13126                        value = CoinMin(saveUpper, value + 1.0e-8 * (fabs(saveUpper) + 1));
    1254813127                    }
    1254913128                    if (value - saveLower < 1.0e-7)
     
    1255913138                    } else {
    1256013139                        // relax a bit
    12561                         value = CoinMax(saveLower, value - 1.0e-5 * (fabs(saveLower) + 1));
     13140                        value = CoinMax(saveLower, value - 1.0e-8 * (fabs(saveLower) + 1));
    1256213141                    }
    1256313142                    if (saveUpper - value < 1.0e-7)
     
    1260713186                                    newLower = CoinMax(lower[jColumn],
    1260813187                                                       newLower
    12609                                                        - 1.0e-5 * (fabs(lower[jColumn]) + 1));
     13188                                                       - 1.0e-8 * (fabs(lower[jColumn]) + 1));
    1261013189                                    newUpper = CoinMin(upper[jColumn],
    1261113190                                                       newUpper
    12612                                                        + 1.0e-5 * (fabs(upper[jColumn]) + 1));
     13191                                                       + 1.0e-8 * (fabs(upper[jColumn]) + 1));
    1261313192                                }
    1261413193                                solver->setColLower(jColumn, newLower);
     
    1314513724    if (probingInfo_ && currentDepth_ > 0) {
    1314613725        int nFix = probingInfo_->fixColumns(*solver);
     13726#ifdef SWITCH_VARIABLES
     13727        if (nFix>0)
     13728          fixAssociated(solver_,0);
     13729#endif
    1314713730        if (nFix < 0) {
    1314813731#ifdef COIN_HAS_CLP
     
    1324213825#else
    1324313826    solver->resolve();
     13827#endif
     13828#ifdef SWITCH_VARIABLES
     13829    if (solver_->isProvenOptimal()) {
     13830      int nBad=checkAssociated(solver_,solver_->getColSolution(),0);
     13831      if (nBad)
     13832        checkAssociated(solver_,solver_->getColSolution(),1);
     13833    }
    1324413834#endif
    1324513835    return solver->isProvenOptimal() ? 1 : 0;
     
    1473215322    int save2 = maximumDepth_;
    1473315323    int retCode = addCuts(node, lastws, numberFixedNow_ > numberFixedAtRoot_);
     15324#ifdef SWITCH_VARIABLES
     15325    fixAssociated(solver_,0);
     15326#endif
    1473415327    //if (save1<maximumNumberCuts_) {
    1473515328    // increased
     
    1571416307                            }
    1571516308                        } else if (ifSol < 0)   { // just returning an estimate
    15716                             estValue = CoinMin(heurValue, estValue) ;
     16309                            estValue = heurValue; //CoinMin(heurValue, estValue) ;
    1571716310                            heurValue = saveValue ;
    1571816311                        }
     
    1604716640        }
    1604816641    }
    16049 #ifdef CLP_INVESTIGATE4
     16642#ifdef CLP_INVESTIGATE5
    1605016643    if (priority)
    1605116644        printf("Before fathom %d not trusted out of %d\n",
     
    1631416907        OsiClpSolverInterface * clpSolver
    1631516908        = dynamic_cast<OsiClpSolverInterface *> (solver_);
    16316         if (clpSolver && numberNodes_ >= numberNodes && numberNodes_ < 2*numberNodes) {
     16909        if (clpSolver && numberNodes_ >= numberNodes && numberNodes_ < 2*numberNodes  && clpSolver->getNumRows() < 10000) {
    1631716910            if (numberIterations_ < (numberSolves_ + numberNodes_)*10) {
    1631816911                //if (numberIterations_<numberNodes_*20) {
Note: See TracChangeset for help on using the changeset viewer.