Changeset 604


Ignore:
Timestamp:
May 10, 2007 5:36:03 AM (12 years ago)
Author:
forrest
Message:

for heuristics and quadratic

Location:
branches/devel/Cbc/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Cbc/src/CbcLinked.cpp

    r602 r604  
    23552355  }
    23562356  return bestObjectiveValue;
     2357}
     2358/* Solves nonlinear problem from CoinModel using SLP - and then tries to get
     2359   heuristic solution
     2360   Returns solution array
     2361*/
     2362double *
     2363OsiSolverLink::heuristicSolution(int numberPasses,double deltaTolerance,int mode)
     2364{
     2365  // get a solution
     2366  CoinModel tempModel = coinModel_;
     2367  ClpSimplex * temp = approximateSolution(tempModel, numberPasses, deltaTolerance);
     2368  int numberColumns = coinModel_.numberColumns();
     2369  double * solution = CoinCopyOfArray(temp->primalColumnSolution(),numberColumns);
     2370  delete temp;
     2371  OsiClpSolverInterface newSolver;
     2372  if (mode==1) {
     2373    // round all with priority < biLinearPriority_
     2374    setFixedPriority(biLinearPriority_);
     2375    // ? should we save and restore coin model
     2376    tempModel = coinModel_;
     2377    // solve modified problem
     2378    for (int iObject =0;iObject<numberObjects_;iObject++) {
     2379      OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
     2380      if (obj&&obj->priority()<biLinearPriority_) {
     2381        int iColumn = obj->columnNumber();
     2382        double value = solution[iColumn];
     2383        value = ceil(value-1.0e-7);
     2384        tempModel.associateElement(coinModel_.columnName(iColumn),value);
     2385      }
     2386    }
     2387    newSolver.loadFromCoinModel(tempModel,true);
     2388    for (int iObject =0;iObject<numberObjects_;iObject++) {
     2389      OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
     2390      if (obj&&obj->priority()<biLinearPriority_) {
     2391        int iColumn = obj->columnNumber();
     2392        double value = solution[iColumn];
     2393        value = ceil(value-1.0e-7);
     2394        newSolver.setColLower(iColumn,value);
     2395        newSolver.setColUpper(iColumn,value);
     2396      }
     2397    }
     2398  }
     2399  CbcModel model(newSolver);
     2400  // Now do requested saves and modifications
     2401  CbcModel * cbcModel = & model;
     2402  OsiSolverInterface * osiModel = model.solver();
     2403  OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
     2404  ClpSimplex * clpModel = osiclpModel->getModelPtr();
     2405  CglProbing probing;
     2406  probing.setMaxProbe(10);
     2407  probing.setMaxLook(10);
     2408  probing.setMaxElements(200);
     2409  probing.setMaxProbeRoot(50);
     2410  probing.setMaxLookRoot(10);
     2411  probing.setRowCuts(3);
     2412  probing.setRowCuts(0);
     2413  probing.setUsingObjective(true);
     2414  cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
     2415 
     2416  CglGomory gomory;
     2417  gomory.setLimitAtRoot(512);
     2418  cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
     2419 
     2420  CglKnapsackCover knapsackCover;
     2421  cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
     2422 
     2423  CglClique clique;
     2424  clique.setStarCliqueReport(false);
     2425  clique.setRowCliqueReport(false);
     2426  clique.setMinViolation(0.1);
     2427  cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
     2428  CglMixedIntegerRounding2 mixedIntegerRounding2;
     2429  cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
     2430 
     2431  CglFlowCover flowCover;
     2432  cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
     2433 
     2434  CglTwomir twomir;
     2435  twomir.setMaxElements(250);
     2436  cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
     2437  cbcModel->cutGenerator(6)->setTiming(true);
     2438 
     2439  CbcHeuristicFPump heuristicFPump(*cbcModel);
     2440  heuristicFPump.setWhen(1);
     2441  heuristicFPump.setMaximumPasses(20);
     2442  heuristicFPump.setDefaultRounding(0.5);
     2443  cbcModel->addHeuristic(&heuristicFPump);
     2444 
     2445  CbcRounding rounding(*cbcModel);
     2446  cbcModel->addHeuristic(&rounding);
     2447 
     2448  CbcHeuristicLocal heuristicLocal(*cbcModel);
     2449  heuristicLocal.setSearchType(1);
     2450  cbcModel->addHeuristic(&heuristicLocal);
     2451 
     2452  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
     2453  cbcModel->addHeuristic(&heuristicGreedyCover);
     2454 
     2455  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
     2456  cbcModel->addHeuristic(&heuristicGreedyEquality);
     2457 
     2458  CbcCompareDefault compare;
     2459  cbcModel->setNodeComparison(compare);
     2460  cbcModel->setNumberBeforeTrust(5);
     2461  cbcModel->setSpecialOptions(2);
     2462  cbcModel->messageHandler()->setLogLevel(1);
     2463  cbcModel->setMaximumCutPassesAtRoot(-100);
     2464  cbcModel->setMaximumCutPasses(1);
     2465  cbcModel->setMinimumDrop(0.05);
     2466  clpModel->setNumberIterations(1);
     2467  // For branchAndBound this may help
     2468  clpModel->defaultFactorizationFrequency();
     2469  clpModel->setDualBound(6.71523e+07);
     2470  clpModel->setPerturbation(50);
     2471  osiclpModel->setSpecialOptions(193);
     2472  osiclpModel->messageHandler()->setLogLevel(0);
     2473  osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
     2474  osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
     2475  // You can save some time by switching off message building
     2476  // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
     2477  // Solve
     2478 
     2479  cbcModel->initialSolve();
     2480  //double cutoff = model_->getCutoff();
     2481  if (!cbcModel_)
     2482    cbcModel->setCutoff(1.0e50);
     2483  else
     2484    cbcModel->setCutoff(cbcModel_->getCutoff());
     2485  // to change exits
     2486  bool isFeasible=false;
     2487  int saveLogLevel=clpModel->logLevel();
     2488  clpModel->setLogLevel(0);
     2489  int returnCode=0;
     2490  if (clpModel->tightenPrimalBounds()!=0) {
     2491    clpModel->setLogLevel(saveLogLevel);
     2492    returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
     2493    clpModel->writeMps("infeas2.mps");
     2494  } else {
     2495    clpModel->setLogLevel(saveLogLevel);
     2496    clpModel->dual();  // clean up
     2497    // compute some things using problem size
     2498    cbcModel->setMinimumDrop(min(5.0e-2,
     2499                                 fabs(cbcModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
     2500    if (cbcModel->getNumCols()<500)
     2501      cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
     2502    else if (cbcModel->getNumCols()<5000)
     2503      cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop
     2504    else
     2505      cbcModel->setMaximumCutPassesAtRoot(20);
     2506    cbcModel->setMaximumCutPasses(1);
     2507    // Hand coded preprocessing
     2508    CglPreProcess process;
     2509    OsiSolverInterface * saveSolver=cbcModel->solver()->clone();
     2510    // Tell solver we are in Branch and Cut
     2511    saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
     2512    // Default set of cut generators
     2513    CglProbing generator1;
     2514    generator1.setUsingObjective(true);
     2515    generator1.setMaxPass(3);
     2516    generator1.setMaxProbeRoot(saveSolver->getNumCols());
     2517    generator1.setMaxElements(100);
     2518    generator1.setMaxLookRoot(50);
     2519    generator1.setRowCuts(3);
     2520    // Add in generators
     2521    process.addCutGenerator(&generator1);
     2522    process.messageHandler()->setLogLevel(cbcModel->logLevel());
     2523    OsiSolverInterface * solver2 =
     2524      process.preProcessNonDefault(*saveSolver,0,10);
     2525    // Tell solver we are not in Branch and Cut
     2526    saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
     2527    if (solver2)
     2528      solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
     2529    if (!solver2) {
     2530      std::cout<<"Pre-processing says infeasible!"<<std::endl;
     2531      delete saveSolver;
     2532      returnCode=-1;
     2533    } else {
     2534      std::cout<<"processed model has "<<solver2->getNumRows()
     2535               <<" rows, "<<solver2->getNumCols()
     2536               <<" and "<<solver2->getNumElements()<<std::endl;
     2537      // we have to keep solver2 so pass clone
     2538      solver2 = solver2->clone();
     2539      //solver2->writeMps("intmodel");
     2540      cbcModel->assignSolver(solver2);
     2541      cbcModel->initialSolve();
     2542      cbcModel->branchAndBound();
     2543      // For best solution
     2544      int numberColumns = newSolver.getNumCols();
     2545      if (cbcModel->getMinimizationObjValue()<1.0e50) {
     2546        // post process
     2547        process.postProcess(*cbcModel->solver());
     2548        // Solution now back in saveSolver
     2549        cbcModel->assignSolver(saveSolver);
     2550        memcpy(cbcModel->bestSolution(),cbcModel->solver()->getColSolution(),
     2551               numberColumns*sizeof(double));
     2552        // put back in original solver
     2553        newSolver.setColSolution(cbcModel->bestSolution());
     2554        isFeasible=true;
     2555      } else {
     2556        delete saveSolver;
     2557      }
     2558    }
     2559  }
     2560  assert (!returnCode);
     2561  abort();
     2562  return solution;
    23572563}
    23582564// Analyze constraints to see which are convex (quadratic)
     
    64876693#include "ClpConstraint.hpp"
    64886694#include "ClpConstraintLinear.hpp"
     6695#include "ClpConstraintQuadratic.hpp"
    64896696/* Return an approximate solution to a CoinModel.
    64906697    Lots of bounds may be odd to force a solution.
     
    65016708  // List of nonlinear rows
    65026709  int * which = new int[numberRows+1];
    6503   bool testLinear=true;
     6710  bool testLinear=false;
    65046711  int numberConstraints=0;
    65056712  int iColumn;
     
    65946801    model->dual();
    65956802    return model;
    6596   } 
     6803  }
    65976804  // space for quadratic
     6805  // allow for linear term
     6806  maximumQuadraticElements += numberColumns;
    65986807  CoinBigIndex * startQuadratic = new CoinBigIndex [numberColumns+1];
    65996808  int * columnQuadratic = new int [maximumQuadraticElements];
     
    66056814  if (!linearObjective) {
    66066815    int numberQuadratic=0;
     6816    int largestColumn=-1;
    66076817    CoinZeroN(linearTerm,numberColumns);
    66086818    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     
    66116821      const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
    66126822      if (strcmp(expr,"Numeric")) {
     6823        largestColumn = CoinMax(largestColumn,iColumn);
    66136824        // value*x*y
    66146825        char temp[20000];
     
    66236834            columnQuadratic[numberQuadratic]=jColumn;
    66246835            elementQuadratic[numberQuadratic++]=2.0*value; // convention
     6836            largestColumn = CoinMax(largestColumn,jColumn);
    66256837          } else if (jColumn==-2) {
    66266838            linearTerm[iColumn] = value;
     6839            // and put in as row -1
     6840            columnQuadratic[numberQuadratic]=-1;
     6841            elementQuadratic[numberQuadratic++]=value;
     6842            largestColumn = CoinMax(largestColumn,iColumn);
    66276843          } else {
    66286844            printf("bad nonlinear term %s\n",temp);
     
    66316847          ifFirst=false;
    66326848        }
     6849      } else {
     6850        // linear part
     6851        linearTerm[iColumn] = coinModel.getColumnObjective(iColumn);
     6852        // and put in as row -1
     6853        columnQuadratic[numberQuadratic]=-1;
     6854        elementQuadratic[numberQuadratic++]=linearTerm[iColumn];
     6855        if (linearTerm[iColumn])
     6856          largestColumn = CoinMax(largestColumn,iColumn);
    66336857      }
    66346858    }
    66356859    startQuadratic[numberColumns] = numberQuadratic;
    66366860    // here we create ClpConstraint
    6637     abort();
     6861    constraints[numberConstraints++] = new ClpConstraintQuadratic(-1, largestColumn+1, numberColumns,
     6862                                                                  startQuadratic,columnQuadratic,elementQuadratic);
    66386863  }
    66396864  int iConstraint;
     
    66426867    if (iRow>=0) {
    66436868      int numberQuadratic=0;
     6869      int lastColumn = -1;
     6870      int largestColumn=-1;
    66446871      CoinZeroN(linearTerm,numberColumns);
    66456872      CoinModelLink triple=coinModel.firstInRow(iRow);
    66466873      while (triple.column()>=0) {
    66476874        int iColumn = triple.column();
     6875        while (lastColumn<iColumn) {
     6876          startQuadratic[lastColumn+1]=numberQuadratic;
     6877          lastColumn++;
     6878        }
    66486879        const char *  expr = coinModel.getElementAsString(iRow,iColumn);
    66496880        if (strcmp("Numeric",expr)) {
     6881          largestColumn = CoinMax(largestColumn,iColumn);
    66506882          // value*x*y
    66516883          char temp[20000];
     
    66606892              columnQuadratic[numberQuadratic]=jColumn;
    66616893              elementQuadratic[numberQuadratic++]=2.0*value; // convention
     6894              largestColumn = CoinMax(largestColumn,jColumn);
    66626895            } else if (jColumn==-2) {
    66636896              linearTerm[iColumn] = value;
     6897              // and put in as row -1
     6898              columnQuadratic[numberQuadratic]=-1;
     6899              elementQuadratic[numberQuadratic++]=value;
     6900              largestColumn = CoinMax(largestColumn,iColumn);
    66646901            } else {
    66656902              printf("bad nonlinear term %s\n",temp);
     
    66686905            ifFirst=false;
    66696906          }
     6907        } else {
     6908          // linear part
     6909          linearTerm[iColumn] = coinModel.getElement(iRow,iColumn);
     6910          // and put in as row -1
     6911          columnQuadratic[numberQuadratic]=-1;
     6912          elementQuadratic[numberQuadratic++]=linearTerm[iColumn];
     6913          if (linearTerm[iColumn])
     6914            largestColumn = CoinMax(largestColumn,iColumn);
    66706915        }
    66716916        triple=coinModel.next(triple);
    66726917      }
    6673       startQuadratic[numberColumns] = numberQuadratic;
     6918      while (lastColumn<numberColumns) {
     6919        startQuadratic[lastColumn+1]=numberQuadratic;
     6920        lastColumn++;
     6921      }
    66746922      // here we create ClpConstraint
    66756923      if (testLinear) {
     
    66876935        delete [] indices;
    66886936      } else {
    6689         abort();
     6937        constraints[numberConstraints++] = new ClpConstraintQuadratic(iRow, largestColumn+1, numberColumns,
     6938                      startQuadratic,columnQuadratic,elementQuadratic);
    66906939      }
    66916940    }
     
    66996948  int returnCode = model->nonlinearSLP(numberConstraints, constraints,
    67006949                                       numberPasses,deltaTolerance);
     6950  // See if any integers
    67016951  for (iConstraint=0;iConstraint<saveNumber;iConstraint++)
    67026952    delete constraints[iConstraint];
     
    67046954  delete [] constraints;
    67056955  assert (!returnCode);
    6706   abort();
    67076956  return model;
    67086957}
  • branches/devel/Cbc/src/CbcLinked.hpp

    r602 r604  
    6060  */
    6161  double linearizedBAB(CglStored * cut) ;
     62  /** Solves nonlinear problem from CoinModel using SLP - and then tries to get
     63      heuristic solution
     64      Returns solution array
     65  */
     66  double * heuristicSolution(int numberPasses,double deltaTolerance,int mode);
    6267  //@}
    6368 
  • branches/devel/Cbc/src/CbcSolver.cpp

    r603 r604  
    22592259              ClpSimplex * model2 = lpSolver;
    22602260              if (dualize) {
     2261                //printf("dualize %d\n",dualize);
    22612262                model2 = ((ClpSimplexOther *) model2)->dualOfModel();
    22622263                printf("Dual of model has %d rows and %d columns\n",
     
    23802381                double * solution = NULL;
    23812382                if (testOsiOptions<10) {
    2382                   solution = linkSolver->nonlinearSLP(slpValue,1.0e-5);
    2383                 } else {
     2383                  solution = linkSolver->nonlinearSLP(slpValue>0 ? slpValue : 20 ,1.0e-5);
     2384                } else if (testOsiOptions==10) {
    23842385                  CoinModel coinModel = *linkSolver->coinModel();
    23852386                  ClpSimplex * tempModel = approximateSolution(coinModel,slpValue>0 ? slpValue : 20 ,1.0e-5,0);
     
    40134014                      //choose.setNumberBeforeTrusted(2000);
    40144015                      //choose.setNumberStrong(20);
     4016                    }
     4017                    // For temporary testing of heuristics
     4018                    //int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
     4019                    if (testOsiOptions>=10) {
     4020                      printf("*** Temp heuristic with mode %d\n",testOsiOptions-10);
     4021                      OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
     4022                      assert (solver3) ;
     4023                      double * solution = solver3->heuristicSolution(slpValue>0 ? slpValue : 20 ,1.0e-5,testOsiOptions-10);
     4024                      assert(solution);
    40154025                    }
    40164026#endif
Note: See TracChangeset for help on using the changeset viewer.