Ignore:
Timestamp:
May 31, 2011 4:10:30 AM (8 years ago)
Author:
forrest
Message:

allow end cuts and lagomory

File:
1 edited

Legend:

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

    r1654 r1656  
    19811981#endif
    19821982    bool feasible;
     1983    numberSolves_ = 0 ;
    19831984    // If NLP then we assume already solved outside branchAndbound
    19841985    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
     
    22232224        if (clpSolver)
    22242225            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
     2226#ifdef CLP_RESOLVE
     2227        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
     2228          resolveClp(clpSolver,0);
     2229        }
     2230#endif
    22252231    }
    22262232#endif
     
    23322338    */
    23332339    numberIterations_ = 0 ;
    2334     numberSolves_ = 0 ;
    23352340    numberNodes_ = 0 ;
    23362341    numberNodes2_ = 0 ;
     
    25612566        feasible = false;
    25622567
     2568#ifdef CLP_RESOLVE
     2569    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
     2570      OsiClpSolverInterface * clpSolver
     2571        = dynamic_cast<OsiClpSolverInterface *> (solver_);
     2572      if (clpSolver)
     2573        resolveClp(clpSolver,0);
     2574    }
     2575#endif
    25632576    // make cut generators less aggressive
    25642577    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
     
    69206933    // Whether to increase minimum drop
    69216934    bool increaseDrop = (moreSpecialOptions_ & 8192) != 0;
     6935    for (int i = 0; i < numberCutGenerators_; i++)
     6936      generator_[i]->setWhetherInMustCallAgainMode(false);
    69226937    /*
    69236938      Begin cut generation loop. Cuts generated during each iteration are
     
    69436958                if (!generator_[i]->mustCallAgain())
    69446959                    generator_[i]->setSwitchedOff(true);
     6960                else
     6961                  generator_[i]->setWhetherInMustCallAgainMode(true);
    69456962            }
    69466963        }
     
    75307547            keepGoing=false;
    75317548        }
     7549        if (numberTries <=0 && feasible && !keepGoing && !parentModel_ && !numberNodes_) {
     7550          for (int i = 0; i < numberCutGenerators_; i++) {
     7551            if (generator_[i]->whetherCallAtEnd()
     7552                &&!generator_[i]->whetherInMustCallAgainMode()) {
     7553              keepGoing= true;
     7554              break;
     7555            }
     7556          }
     7557        }
    75327558    } while (numberTries > 0 || keepGoing) ;
    75337559    /*
     
    83458371        int numberRowCutsAfter = numberRowCutsBefore;
    83468372        int numberColumnCutsAfter = numberColumnCutsBefore;
     8373        /*printf("GEN %d %s switches %d\n",
     8374               i,generator_[i]->cutGeneratorName(),
     8375               generator_[i]->switches());*/
    83478376        bool generate = generator_[i]->normal();
    83488377        // skip if not optimal and should be (maybe a cut generator has fixed variables)
     
    83758404            }
    83768405            if (numberRowCutsBefore < numberRowCutsAfter &&
    8377                     generator_[i]->mustCallAgain() && status >= 0)
     8406                generator_[i]->mustCallAgain() && status >= 0)
     8407              /*printf("%s before %d after %d must %c atend %c off %c endmode %c\n",
     8408                   generator_[i]->cutGeneratorName(),
     8409                   numberRowCutsBefore,numberRowCutsAfter,
     8410                   generator_[i]->mustCallAgain() ? 'Y': 'N',
     8411                   generator_[i]->whetherCallAtEnd() ? 'Y': 'N',
     8412                   generator_[i]->switchedOff() ? 'Y': 'N',
     8413                   generator_[i]->whetherInMustCallAgainMode() ? 'Y': 'N');*/
     8414            if (numberRowCutsBefore < numberRowCutsAfter &&
     8415                generator_[i]->mustCallAgain() && status >= 0)
    83788416                status = 1 ; // say must go round
    83798417            // Check last cut to see if infeasible
     
    1210712145    return solver->isProvenOptimal() ? 1 : 0;
    1210812146}
     12147#ifdef CLP_RESOLVE
     12148// Special purpose resolve
     12149int
     12150CbcModel::resolveClp(OsiClpSolverInterface * clpSolver, int type)
     12151{
     12152    numberSolves_++;
     12153    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     12154    int save = clpSimplex->specialOptions();
     12155    clpSimplex->setSpecialOptions(save | 0x11000000); // say is Cbc (and in branch and bound)
     12156    int save2 = clpSolver->specialOptions();
     12157    clpSolver->resolve();
     12158    if (!numberNodes_) {
     12159      double error = CoinMax(clpSimplex->largestDualError(),
     12160                             clpSimplex->largestPrimalError());
     12161      if (error > 1.0e-2 || !clpSolver->isProvenOptimal()) {
     12162#ifdef CLP_INVESTIGATE
     12163        printf("Problem was %s largest dual error %g largest primal %g - safer cuts\n",
     12164               clpSolver->isProvenOptimal() ? "optimal" : "!infeasible",
     12165               clpSimplex->largestDualError(),
     12166               clpSimplex->largestPrimalError());
     12167#endif
     12168        if (!clpSolver->isProvenOptimal()) {
     12169          clpSolver->setSpecialOptions(save2 | 2048);
     12170          clpSimplex->allSlackBasis(true);
     12171          clpSolver->resolve();
     12172          if (!clpSolver->isProvenOptimal()) {
     12173            bool takeHint;
     12174            OsiHintStrength strength;
     12175            clpSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
     12176            clpSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
     12177            clpSolver->resolve();
     12178            clpSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
     12179          }
     12180        }
     12181        // make cuts safer
     12182        for (int iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
     12183          CglCutGenerator * generator = generator_[iCutGenerator]->generator();
     12184          CglGomory * cgl1 = dynamic_cast<CglGomory *>(generator);
     12185          if (cgl1) {
     12186            cgl1->setLimitAtRoot(cgl1->getLimit());
     12187          }
     12188          CglTwomir * cgl2 = dynamic_cast<CglTwomir *>(generator);
     12189          if (cgl2) {
     12190            generator_[iCutGenerator]->setHowOften(-100);
     12191          }
     12192        }
     12193      }
     12194    }
     12195    clpSolver->setSpecialOptions(save2);
     12196#ifdef CLP_INVESTIGATE
     12197    if (clpSimplex->numberIterations() > 1000)
     12198      printf("node %d took %d iterations\n", numberNodes_, clpSimplex->numberIterations());
     12199#endif
     12200    if (type==0 && clpSolver->isProvenOptimal()) {
     12201      ClpSimplex newModel(*clpSimplex);
     12202      newModel.primal();
     12203      int numberColumns = newModel.numberColumns();
     12204      int numberRows = newModel.numberRows();
     12205      double * obj = new double [numberColumns];
     12206      int * which = new int [numberColumns];
     12207      const double * solution = clpSimplex->primalColumnSolution();
     12208      double rhs=1.0e-8;
     12209      int numberObj=0;
     12210      double integerTolerance = getDblParam(CbcIntegerTolerance) ;
     12211      double * objective = newModel.objective();
     12212      for (int i=0;i<numberColumns;i++) {
     12213        if (objective[i]) {
     12214          rhs += objective[i]*solution[i];
     12215          obj[numberObj]=objective[i];
     12216          which[numberObj++]=i;
     12217          objective[i]=0.0;
     12218        }
     12219      }
     12220      if (numberObj) {
     12221        newModel.addRow(numberObj,which,obj,-COIN_DBL_MAX,rhs);
     12222      }
     12223      delete [] obj;
     12224      delete [] which;
     12225      double * lower = newModel.columnLower();
     12226      double * upper = newModel.columnUpper();
     12227      int numberInf = 0;
     12228      int numberLb = 0;
     12229      int numberUb = 0;
     12230      int numberInt = 0;
     12231      double sumInf=0.0;
     12232      for (int i=0;i<numberIntegers_;i++) {
     12233        int iSequence = integerVariable_[i];
     12234        double value = solution[iSequence];
     12235        value = CoinMax(value, lower[iSequence]);
     12236        value = CoinMin(value, upper[iSequence]);
     12237        double nearest = floor(value + 0.5);
     12238        if (value<lower[iSequence]+integerTolerance) {
     12239          objective[iSequence]=1.0;
     12240          numberLb++;
     12241        } else if (value>upper[iSequence]-integerTolerance) {
     12242          objective[iSequence]=-1.0;
     12243          numberUb++;
     12244        } else if (fabs(value - nearest) <= integerTolerance) {
     12245          // fix??
     12246          lower[iSequence]=nearest;
     12247          upper[iSequence]=nearest;
     12248          numberInt++;
     12249        } else {
     12250          lower[iSequence]=floor(value);
     12251          upper[iSequence]=ceil(value);
     12252          if (value>nearest) {
     12253            objective[iSequence]=1.0;
     12254            sumInf += value-nearest;
     12255          } else {
     12256            objective[iSequence]=-1.0;
     12257            sumInf -= value-nearest;
     12258          }
     12259          numberInf++;
     12260        }
     12261      }
     12262      printf("XX %d inf (sum %g), %d at lb %d at ub %d other integer\n",
     12263             numberInf,sumInf,numberLb,numberUb,numberInt);
     12264      if (numberInf) {
     12265        newModel.primal(1);
     12266        if (!newModel.isProvenOptimal()) {
     12267          printf("not optimal - scaling issue - switch off\n");
     12268          clpSimplex->setSpecialOptions(save);
     12269          if (clpSimplex->status()==4)
     12270            clpSimplex->setProblemStatus(1);
     12271          return clpSolver->isProvenOptimal() ? 1 : 0;
     12272        }
     12273        //newModel.writeMps("bad.mps");
     12274        //assert (newModel.isProvenOptimal());
     12275        printf("%d iterations\n",newModel.numberIterations());
     12276        int numberInf2 = 0;
     12277        int numberLb2 = 0;
     12278        int numberUb2 = 0;
     12279        int numberInt2 = 0;
     12280        double sumInf2=0.0;
     12281        const double * solution = newModel.primalColumnSolution();
     12282        const double * lower = clpSimplex->columnLower();
     12283        const double * upper = clpSimplex->columnUpper();
     12284        for (int i=0;i<numberIntegers_;i++) {
     12285          int iSequence = integerVariable_[i];
     12286          double value = solution[iSequence];
     12287          value = CoinMax(value, lower[iSequence]);
     12288          value = CoinMin(value, upper[iSequence]);
     12289          double nearest = floor(value + 0.5);
     12290          if (value<lower[iSequence]+integerTolerance) {
     12291            numberLb2++;
     12292          } else if (value>upper[iSequence]-integerTolerance) {
     12293            numberUb2++;
     12294          } else if (fabs(value - nearest) <= integerTolerance) {
     12295            numberInt2++;
     12296          } else {
     12297            if (value>nearest) {
     12298              sumInf2 += value-nearest;
     12299            } else {
     12300              sumInf2 -= value-nearest;
     12301            }
     12302            numberInf2++;
     12303          }
     12304        }
     12305        printf("XXX %d inf (sum %g), %d at lb %d at ub %d other integer\n",
     12306               numberInf2,sumInf2,numberLb2,numberUb2,numberInt2);
     12307        if (sumInf2<sumInf*0.95) {
     12308          printf("XXXX suminf reduced from %g (%d) to %g (%d)\n",
     12309                 sumInf,numberInf,sumInf2,numberInf2);
     12310          if (numberObj) {
     12311            newModel.deleteRows(1,&numberRows);
     12312          }
     12313          memcpy(newModel.objective(),
     12314                 clpSimplex->objective(),
     12315                 numberColumns*sizeof(double));
     12316          memcpy(newModel.columnLower(),
     12317                 clpSimplex->columnLower(),
     12318                 numberColumns*sizeof(double));
     12319          memcpy(newModel.columnUpper(),
     12320                 clpSimplex->columnUpper(),
     12321                 numberColumns*sizeof(double));
     12322          newModel.setClpScaledMatrix(NULL);
     12323          newModel.primal(1);
     12324          printf("%d iterations\n",newModel.numberIterations());
     12325          int numberInf3 = 0;
     12326          int numberLb3 = 0;
     12327          int numberUb3 = 0;
     12328          int numberInt3 = 0;
     12329          double sumInf3=0.0;
     12330          const double * solution = newModel.primalColumnSolution();
     12331          const double * lower = clpSimplex->columnLower();
     12332          const double * upper = clpSimplex->columnUpper();
     12333          for (int i=0;i<numberIntegers_;i++) {
     12334            int iSequence = integerVariable_[i];
     12335            double value = solution[iSequence];
     12336            value = CoinMax(value, lower[iSequence]);
     12337            value = CoinMin(value, upper[iSequence]);
     12338            double nearest = floor(value + 0.5);
     12339            if (value<lower[iSequence]+integerTolerance) {
     12340              numberLb3++;
     12341            } else if (value>upper[iSequence]-integerTolerance) {
     12342              numberUb3++;
     12343            } else if (fabs(value - nearest) <= integerTolerance) {
     12344              numberInt3++;
     12345            } else {
     12346              if (value>nearest) {
     12347                sumInf3 += value-nearest;
     12348              } else {
     12349                sumInf3 -= value-nearest;
     12350              }
     12351              numberInf3++;
     12352            }
     12353          }
     12354          printf("XXXXX %d inf (sum %g), %d at lb %d at ub %d other integer\n",
     12355                 numberInf3,sumInf3,numberLb3,numberUb3,numberInt3);
     12356          if (sumInf3<sumInf*0.95) {
     12357            memcpy(clpSimplex->primalColumnSolution(),
     12358                   newModel.primalColumnSolution(),
     12359                   numberColumns*sizeof(double));
     12360            memcpy(clpSimplex->dualColumnSolution(),
     12361                   newModel.dualColumnSolution(),
     12362                   numberColumns*sizeof(double));
     12363            memcpy(clpSimplex->primalRowSolution(),
     12364                   newModel.primalRowSolution(),
     12365                   numberRows*sizeof(double));
     12366            memcpy(clpSimplex->dualRowSolution(),
     12367                   newModel.dualRowSolution(),
     12368                   numberRows*sizeof(double));
     12369            memcpy(clpSimplex->statusArray(),
     12370                   newModel.statusArray(),
     12371                   (numberColumns+numberRows)*sizeof(unsigned char));
     12372            clpSolver->setWarmStart(NULL);
     12373          }
     12374        }
     12375      }
     12376    }
     12377    clpSimplex->setSpecialOptions(save);
     12378    if (clpSimplex->status()==4)
     12379      clpSimplex->setProblemStatus(1);
     12380    return clpSolver->isProvenOptimal() ? 1 : 0;
     12381}
     12382#endif
    1210912383/*!
    1211012384    \todo It'd be really nice if there were an overload for this method that
Note: See TracChangeset for help on using the changeset viewer.