Changeset 523 for branches/devel


Ignore:
Timestamp:
Jan 24, 2007 4:52:57 PM (13 years ago)
Author:
forrest
Message:

nonlinear stuff

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

Legend:

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

    r491 r523  
    988988  olb = model_->solver()->getColLower()[iColumn] ;
    989989  oub = model_->solver()->getColUpper()[iColumn] ;
     990#ifdef COIN_DEVELOP
     991  if (olb!=down_[0]||oub!=up_[1]) {
     992    if (way_>0)
     993      printf("branching up on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
     994             iColumn,olb,oub,up_[0],up_[1],down_[0],down_[1]) ;
     995    else
     996      printf("branching down on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
     997             iColumn,olb,oub,down_[0],down_[1],up_[0],up_[1]) ;
     998  }
     999#endif
    9901000  if (way_<0) {
    9911001#ifdef CBC_DEBUG
     
    10121022  }
    10131023  double nlb = model_->solver()->getColLower()[iColumn];
     1024  double nub = model_->solver()->getColUpper()[iColumn];
    10141025  if (nlb<olb) {
    10151026#ifndef NDEBUG
    10161027    printf("bad lb change for column %d from %g to %g\n",iColumn,olb,nlb);
    10171028#endif
    1018     model_->solver()->setColLower(iColumn,olb);
    1019   }
    1020   double nub = model_->solver()->getColUpper()[iColumn];
     1029    model_->solver()->setColLower(iColumn,CoinMin(olb,nub));
     1030    nlb=olb;
     1031  }
    10211032  if (nub>oub) {
    10221033#ifndef NDEBUG
    10231034    printf("bad ub change for column %d from %g to %g\n",iColumn,oub,nub);
    10241035#endif
    1025     model_->solver()->setColUpper(iColumn,oub);
     1036    model_->solver()->setColUpper(iColumn,CoinMax(oub,nlb));
    10261037  }
    10271038#ifndef NDEBUG
  • branches/devel/Cbc/src/CbcLinked.cpp

    r520 r523  
    3131#include "CbcCutGenerator.hpp"
    3232#include "CglStored.hpp"
    33 //#include "CglTemporary.hpp"
     33#include "CglPreProcess.hpp"
     34#include "CglGomory.hpp"
     35#include "CglProbing.hpp"
     36#include "CglKnapsackCover.hpp"
     37#include "CglRedSplit.hpp"
     38#include "CglClique.hpp"
     39#include "CglFlowCover.hpp"
     40#include "CglMixedIntegerRounding2.hpp"
     41#include "CglTwomir.hpp"
     42#include "CglDuplicateRow.hpp"
     43#include "CbcHeuristicFPump.hpp"
     44#include "CbcHeuristic.hpp"
     45#include "CbcHeuristicLocal.hpp"
     46#include "CbcHeuristicGreedy.hpp"
     47#include "ClpLinearObjective.hpp"
     48#include "CbcBranchActual.hpp"
     49#include "CbcCompareActual.hpp"
    3450//#############################################################################
    3551// Solve methods
     
    183199  specialOptions_ =0;
    184200  modelPtr_->setWhatsChanged(0);
    185   bool allFixed=false;
     201  bool allFixed=numberFix_>0;
    186202  bool feasible=true;
    187203  if (numberVariables_) {
    188204    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
    189     allFixed=true;
    190205    //bool best=true;
    191206    const double * lower = modelPtr_->columnLower();
     
    197212      double lo = lower[iColumn];
    198213      double up = upper[iColumn];
    199       if (up>lo)
    200         allFixed=false;
    201       else if (up<lo)
     214      if (up<lo)
    202215        feasible=false;
    203216    }
     
    284297  if (!feasible)
    285298    allFixed=false;
    286   if ((specialOptions2_&1)!=0)
     299  if ((specialOptions2_&1)==0)
    287300    allFixed=false;
    288301  int returnCode=-1;
     
    292305    if(maxIts>10000) {
    293306      // may do lots of work
    294       returnCode=fathom(allFixed);
     307      if ((specialOptions2_&1)!=0) {
     308        // see if fixed
     309        const double * lower = modelPtr_->columnLower();
     310        const double * upper = modelPtr_->columnUpper();
     311        for (int i=0;i<numberFix_;i++ ) {
     312          int iColumn = fixVariables_[i];
     313          double lo = lower[iColumn];
     314          double up = upper[iColumn];
     315          if (up>lo) {
     316            allFixed=false;
     317            break;
     318          }
     319        }
     320        returnCode=allFixed ? fathom(allFixed) : 0;
     321      } else {
     322        returnCode=0;
     323      }
    295324    } else {
    296325      returnCode=0;
     
    300329    if (returnCode==0)
    301330      OsiClpSolverInterface::resolve();
    302     if (!allFixed&&(specialOptions2_&1)==0) {
     331    if (!allFixed&&(specialOptions2_&1)!=0) {
    303332      const double * solution = getColSolution();
    304333      bool satisfied=true;
     
    533562      // ???  - try
    534563      // But skip if strong branching
    535       CbcModel * cbcModel = (modelPtr_->maximumIterations()<10000) ? cbcModel_ : NULL;
     564      CbcModel * cbcModel = (modelPtr_->maximumIterations()>10000) ? cbcModel_ : NULL;
    536565      if ((specialOptions2_&2)!=0) {
    537566        // If model has stored then add cut (if convex)
    538         if (cbcModel&&(specialOptions2_&4)!=0&&quadraticModel_) {
     567        // off until I work out problem with ibell3a
     568        if (cbcModel&&(specialOptions2_&4)!=0&&quadraticModel_&&false) {
    539569          int numberGenerators = cbcModel_->numberCutGenerators();
    540570          int iGenerator;
     
    879909  return jColumn;
    880910}
    881 void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds)
     911void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds,int logLevel)
    882912{
    883913  // first check and set up arrays
     
    10401070      tempModel.deleteRows(nDelete,freeRow);
    10411071      tempModel.setOptimizationDirection(1.0);
     1072      if (logLevel<3)
     1073        tempModel.setLogLevel(0);
    10421074      double * objective = tempModel.objective();
    10431075      CoinZeroN(objective,numberColumns);
     
    10541086            if (coinModel_.isInteger(iColumn))
    10551087              value = ceil(value-0.9e-3);
    1056             printf("lower bound on %d changed from %g to %g\n",iColumn,columnLower[iColumn],value);
     1088            if (logLevel>1)
     1089              printf("lower bound on %d changed from %g to %g\n",iColumn,columnLower[iColumn],value);
    10571090            columnLower[iColumn]=value;
     1091            coinModel_.setColumnLower(iColumn,value);
    10581092          }
    10591093          objective[iColumn]=-1.0;
     
    10631097            if (coinModel_.isInteger(iColumn))
    10641098              value = floor(value+0.9e-3);
    1065             printf("upper bound on %d changed from %g to %g\n",iColumn,columnUpper[iColumn],value);
     1099            if (logLevel>1)
     1100              printf("upper bound on %d changed from %g to %g\n",iColumn,columnUpper[iColumn],value);
    10661101            columnUpper[iColumn]=value;
     1102            coinModel_.setColumnUpper(iColumn,value);
    10671103          }
    10681104          objective[iColumn]=0.0;
     
    14231459  delete [] which;
    14241460}
     1461// Set all biLinear priorities on x-x variables
     1462void
     1463OsiSolverLink::setBiLinearPriorities(int value)
     1464{
     1465  int i;
     1466  for ( i =0;i<numberObjects_;i++) {
     1467    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1468    if (obj) {
     1469      if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
     1470        obj->setPriority(value);
     1471      }
     1472    }
     1473  }
     1474}
     1475// Set all mesh sizes on x-x variables
     1476void
     1477OsiSolverLink::setMeshSizes(double value)
     1478{
     1479  int i;
     1480  for ( i =0;i<numberObjects_;i++) {
     1481    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1482    if (obj) {
     1483      if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
     1484#if 0
     1485        numberContinuous++;
     1486        int xColumn = obj->xColumn();
     1487        double gapX = upper[xColumn]-lower[xColumn];
     1488        int yColumn = obj->yColumn();
     1489        double gapY = upper[yColumn]-lower[yColumn];
     1490        gap = CoinMax(gap,CoinMax(gapX,gapY));
     1491#endif
     1492        obj->setXMeshSize(value);
     1493        obj->setYMeshSize(value);
     1494      }
     1495    }
     1496  }
     1497}
    14251498/* Solves nonlinear problem from CoinModel using SLP - may be used as crash
    14261499   for other algorithms when number of iterations small.
     
    20832156  return CoinCopyOfArray(solution,coinModel.numberColumns());
    20842157}
     2158/* Solve linearized quadratic objective branch and bound.
     2159   Return cutoff and OA cut
     2160*/
     2161double
     2162OsiSolverLink::linearizedBAB(CglStored * cut)
     2163{
     2164  double bestObjectiveValue=COIN_DBL_MAX;
     2165  if (quadraticModel_) {
     2166    ClpSimplex * qp = new ClpSimplex(*quadraticModel_);
     2167    // bounds
     2168    int numberColumns = qp->numberColumns();
     2169    double * lower = qp->columnLower();
     2170    double * upper = qp->columnUpper();
     2171    const double * lower2 = getColLower();
     2172    const double * upper2 = getColUpper();
     2173    for (int i=0;i<numberColumns;i++) {
     2174      lower[i] = CoinMax(lower[i],lower2[i]);
     2175      upper[i] = CoinMin(upper[i],upper2[i]);
     2176    }
     2177    qp->nonlinearSLP(20,1.0e-5);
     2178    qp->primal();
     2179    OsiSolverLinearizedQuadratic solver2(qp);
     2180    const double * solution=NULL;
     2181    // Reduce printout
     2182    solver2.setHintParam(OsiDoReducePrint,true,OsiHintTry);
     2183    CbcModel model2(solver2);
     2184    // Now do requested saves and modifications
     2185    CbcModel * cbcModel = & model2;
     2186    OsiSolverInterface * osiModel = model2.solver();
     2187    OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
     2188    ClpSimplex * clpModel = osiclpModel->getModelPtr();
     2189   
     2190    // Set changed values
     2191   
     2192    CglProbing probing;
     2193    probing.setMaxProbe(10);
     2194    probing.setMaxLook(10);
     2195    probing.setMaxElements(200);
     2196    probing.setMaxProbeRoot(50);
     2197    probing.setMaxLookRoot(10);
     2198    probing.setRowCuts(3);
     2199    probing.setUsingObjective(true);
     2200    cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
     2201    cbcModel->cutGenerator(0)->setTiming(true);
     2202   
     2203    CglGomory gomory;
     2204    gomory.setLimitAtRoot(512);
     2205    cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
     2206    cbcModel->cutGenerator(1)->setTiming(true);
     2207   
     2208    CglKnapsackCover knapsackCover;
     2209    cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
     2210    cbcModel->cutGenerator(2)->setTiming(true);
     2211   
     2212    CglClique clique;
     2213    clique.setStarCliqueReport(false);
     2214    clique.setRowCliqueReport(false);
     2215    clique.setMinViolation(0.1);
     2216    cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
     2217    cbcModel->cutGenerator(3)->setTiming(true);
     2218   
     2219    CglMixedIntegerRounding2 mixedIntegerRounding2;
     2220    cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
     2221    cbcModel->cutGenerator(4)->setTiming(true);
     2222   
     2223    CglFlowCover flowCover;
     2224    cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
     2225    cbcModel->cutGenerator(5)->setTiming(true);
     2226   
     2227    CglTwomir twomir;
     2228    twomir.setMaxElements(250);
     2229    cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
     2230    cbcModel->cutGenerator(6)->setTiming(true);
     2231    // For now - switch off most heuristics (because CglPreProcess is bad with QP)
     2232#if 0   
     2233    CbcHeuristicFPump heuristicFPump(*cbcModel);
     2234    heuristicFPump.setWhen(13);
     2235    heuristicFPump.setMaximumPasses(20);
     2236    heuristicFPump.setMaximumRetries(7);
     2237    heuristicFPump.setAbsoluteIncrement(4332.64);
     2238    cbcModel->addHeuristic(&heuristicFPump);
     2239    heuristicFPump.setInitialWeight(1);
     2240   
     2241    CbcHeuristicLocal heuristicLocal(*cbcModel);
     2242    heuristicLocal.setSearchType(1);
     2243    cbcModel->addHeuristic(&heuristicLocal);
     2244   
     2245    CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
     2246    cbcModel->addHeuristic(&heuristicGreedyCover);
     2247   
     2248    CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
     2249    cbcModel->addHeuristic(&heuristicGreedyEquality);
     2250#endif
     2251   
     2252    CbcRounding rounding(*cbcModel);
     2253    cbcModel->addHeuristic(&rounding);
     2254   
     2255    cbcModel->setNumberBeforeTrust(5);
     2256    cbcModel->setSpecialOptions(2);
     2257    cbcModel->messageHandler()->setLogLevel(1);
     2258    cbcModel->setMaximumCutPassesAtRoot(-100);
     2259    cbcModel->setMaximumCutPasses(1);
     2260    cbcModel->setMinimumDrop(0.05);
     2261    // For branchAndBound this may help
     2262    clpModel->defaultFactorizationFrequency();
     2263    clpModel->setDualBound(1.0001e+08);
     2264    clpModel->setPerturbation(50);
     2265    osiclpModel->setSpecialOptions(193);
     2266    osiclpModel->messageHandler()->setLogLevel(0);
     2267    osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
     2268    osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
     2269    // You can save some time by switching off message building
     2270    // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
     2271   
     2272    // Solve
     2273   
     2274    cbcModel->initialSolve();
     2275    if (clpModel->tightenPrimalBounds()!=0) {
     2276      std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
     2277      delete qp;
     2278      return COIN_DBL_MAX;
     2279    }
     2280    clpModel->dual();  // clean up
     2281    cbcModel->initialSolve();
     2282    cbcModel->branchAndBound();
     2283    OsiSolverLinearizedQuadratic * solver3 = dynamic_cast<OsiSolverLinearizedQuadratic *> (model2.solver());
     2284    assert (solver3);
     2285    solution = solver3->bestSolution();
     2286    bestObjectiveValue = solver3->bestObjectiveValue();
     2287    setBestObjectiveValue(bestObjectiveValue);
     2288    setBestSolution(solution,solver3->getNumCols());
     2289    // if convex
     2290    if ((specialOptions2()&4)!=0) {
     2291      // add OA cut
     2292      double offset;
     2293      double * gradient = new double [numberColumns+1];
     2294      memcpy(gradient,qp->objectiveAsObject()->gradient(qp,solution,offset,true,2),
     2295             numberColumns*sizeof(double));
     2296      double rhs = 0.0;
     2297      int * column = new int[numberColumns+1];
     2298      int n=0;
     2299      for (int i=0;i<numberColumns;i++) {
     2300        double value = gradient[i];
     2301        if (fabs(value)>1.0e-12) {
     2302          gradient[n]=value;
     2303          rhs += value*solution[i];
     2304          column[n++]=i;
     2305        }
     2306      }
     2307      gradient[n]=-1.0;
     2308      column[n++]=numberColumns;
     2309      cut->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
     2310      delete [] gradient;
     2311      delete [] column;
     2312    }
     2313    delete qp;
     2314    printf("obj %g\n",bestObjectiveValue);
     2315  }
     2316  return bestObjectiveValue;
     2317}
    20852318// Analyze constraints to see which are convex (quadratic)
    20862319void
     
    22782511    delete [] convex_;
    22792512    delete [] whichNonLinear_;
     2513    delete [] fixVariables_;
    22802514  }
    22812515  matrix_ = NULL;
     
    22892523  cbcModel_ = NULL;
    22902524  info_ = NULL;
     2525  fixVariables_=NULL;
    22912526  numberVariables_ = 0;
    22922527  specialOptions2_ = 0;
     
    22992534  integerPriority_ = 1000;
    23002535  biLinearPriority_ = 10000;
     2536  numberFix_=0;
    23012537}
    23022538void
     
    23142550  integerPriority_ = rhs.integerPriority_;
    23152551  biLinearPriority_ = rhs.biLinearPriority_;
     2552  numberFix_ = rhs.numberFix_;
    23162553  cbcModel_ = rhs.cbcModel_;
    23172554  if (numberVariables_) {
     
    23462583    quadraticModel_ = NULL;
    23472584  }
     2585  fixVariables_ = CoinCopyOfArray(rhs.fixVariables_,numberFix_);
    23482586}
    23492587// Add a bound modifier
     
    23962634  CoinZeroN(bestSolution_,numberColumnsThis);
    23972635  memcpy(bestSolution_,solution,CoinMin(numberColumns,numberColumnsThis)*sizeof(double));
     2636}
     2637/* Two tier integer problem where when set of variables with priority
     2638   less than this are fixed the problem becomes an easier integer problem
     2639*/
     2640void
     2641OsiSolverLink::setFixedPriority(int priorityValue)
     2642{
     2643  delete [] fixVariables_;
     2644  fixVariables_=NULL;
     2645  numberFix_=0;
     2646  int i;
     2647  for ( i =0;i<numberObjects_;i++) {
     2648    OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[i]);
     2649    if (obj) {
     2650      int iColumn = obj->columnNumber();
     2651      assert (iColumn>=0);
     2652      if (obj->priority()<priorityValue)
     2653        numberFix_++;
     2654    }
     2655  }
     2656  if (numberFix_) {
     2657    specialOptions2_ |= 1;
     2658    fixVariables_ = new int [numberFix_];
     2659    numberFix_=0;
     2660    // need to make sure coinModel_ is correct
     2661    int numberColumns = coinModel_.numberColumns();
     2662    char * highPriority = new char [numberColumns];
     2663    CoinZeroN(highPriority,numberColumns);
     2664    for ( i =0;i<numberObjects_;i++) {
     2665      OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[i]);
     2666      if (obj) {
     2667        int iColumn = obj->columnNumber();
     2668        assert (iColumn>=0);
     2669        if (iColumn<numberColumns) {
     2670          if (obj->priority()<priorityValue) {
     2671            object_[i]=new OsiSimpleFixedInteger(*obj);
     2672            delete obj;
     2673            fixVariables_[numberFix_++]=iColumn;
     2674            highPriority[iColumn]=2;
     2675          } else {
     2676            highPriority[iColumn]=1;
     2677          }
     2678        }
     2679      }
     2680    }
     2681    // All nonlinear terms must involve high priority (as known)
     2682    double * linear = new double[numberColumns];
     2683    int numberRows = coinModel_.numberRows();
     2684    for (int iRow=0;iRow<numberRows;iRow++) {
     2685      CoinPackedMatrix * row = quadraticRow(iRow,linear);
     2686      if (row) {
     2687        // see if valid
     2688        const double * element = row->getElements();
     2689        const int * columnLow = row->getIndices();
     2690        const CoinBigIndex * columnHigh = row->getVectorStarts();
     2691        const int * columnLength = row->getVectorLengths();
     2692        int numberLook = row->getNumCols();
     2693        int canSwap=0;
     2694        for (int i=0;i<numberLook;i++) {
     2695          // this one needs to be available
     2696          int iPriority = highPriority[i];
     2697          for (int j=columnHigh[i];j<columnHigh[i]+columnLength[i];j++) {
     2698            int iColumn = columnLow[j];
     2699            if (highPriority[iColumn]<=1) {
     2700              assert (highPriority[iColumn]==1);
     2701              if (iPriority==1) {
     2702                canSwap=-1; // no good
     2703                break;
     2704              } else {
     2705                canSwap=1;
     2706              }
     2707            }
     2708          }
     2709        }
     2710        if (canSwap) {
     2711          if (canSwap>0) {
     2712            // rewrite row
     2713            /* get triples
     2714               then swap ones needed
     2715               then create packedmatrix
     2716               then replace row
     2717            */
     2718            int numberElements=columnHigh[numberLook];
     2719            int * columnHigh2 = new int [numberElements];
     2720            int * columnLow2 = new int [numberElements];
     2721            double * element2 = new double [numberElements];
     2722            for (int i=0;i<numberLook;i++) {
     2723              // this one needs to be available
     2724              int iPriority = highPriority[i];
     2725              if (iPriority==2) {
     2726                for (int j=columnHigh[i];j<columnHigh[i]+columnLength[i];j++) {
     2727                  columnHigh2[j]=i;
     2728                  columnLow2[j]=columnLow[j];
     2729                  element2[j]=element[j];
     2730                }
     2731              } else {
     2732                for (int j=columnHigh[i];j<columnHigh[i]+columnLength[i];j++) {
     2733                  columnLow2[j]=i;
     2734                  columnHigh2[j]=columnLow[j];
     2735                  element2[j]=element[j];
     2736                }
     2737              }
     2738            }
     2739            delete row;
     2740            row = new CoinPackedMatrix(true,columnHigh2,columnLow2,element2,numberElements);
     2741            delete [] columnHigh2;
     2742            delete [] columnLow2;
     2743            delete [] element2;
     2744            // Now replace row
     2745            replaceQuadraticRow(iRow,linear,row);
     2746            delete row;
     2747          } else {
     2748            delete row;
     2749            printf("Unable to use priority - row %d\n",iRow);
     2750            delete [] fixVariables_;
     2751            fixVariables_ = NULL;
     2752            numberFix_=0;
     2753            break;
     2754          }
     2755        }
     2756      }
     2757    }
     2758    delete [] highPriority;
     2759    delete [] linear;
     2760  }
     2761}
     2762// Replaces a quadratic row
     2763void
     2764OsiSolverLink::replaceQuadraticRow(int rowNumber,const double * linearRow, const CoinPackedMatrix * quadraticPart)
     2765{
     2766  int numberColumns = coinModel_.numberColumns();
     2767  int numberRows = coinModel_.numberRows();
     2768  assert (rowNumber>=0&&rowNumber<numberRows);
     2769  CoinModelLink triple=coinModel_.firstInRow(rowNumber);
     2770  while (triple.column()>=0) {
     2771    int iColumn = triple.column();
     2772    coinModel_.deleteElement(rowNumber,iColumn);
     2773    // triple stale - so start over
     2774    triple=coinModel_.firstInRow(rowNumber);
     2775  }
     2776  const double * element = quadraticPart->getElements();
     2777  const int * column = quadraticPart->getIndices();
     2778  const CoinBigIndex * columnStart = quadraticPart->getVectorStarts();
     2779  const int * columnLength = quadraticPart->getVectorLengths();
     2780  int numberLook = quadraticPart->getNumCols();
     2781  int i;
     2782  for (i=0;i<numberLook;i++) {
     2783    if (!columnLength[i]) {
     2784      // just linear part
     2785      if (linearRow[i])
     2786        coinModel_.setElement(rowNumber,i,linearRow[i]);
     2787    } else {
     2788      char temp[10000];
     2789      int put=0;
     2790      char temp2[30];
     2791      bool first=true;
     2792      if (linearRow[i]) {
     2793        sprintf(temp,"%g",linearRow[i]);
     2794        first=false;
     2795        put = strlen(temp);
     2796      }
     2797      for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     2798        int jColumn = column[j];
     2799        double value = element[j];
     2800        if (value<0.0||first)
     2801          sprintf(temp2,"%g*c%7.7d",value,jColumn);
     2802        else
     2803          sprintf(temp2,"+%g*c%7.7d",value,jColumn);
     2804        int nextPut = put + strlen(temp2);
     2805        assert (nextPut<10000);
     2806        strcpy(temp+put,temp2);
     2807        put = nextPut;
     2808      }
     2809      coinModel_.setElement(rowNumber,i,temp);
     2810    }
     2811  }
     2812  // rest of linear
     2813  for (;i<numberColumns;i++) {
     2814    if (linearRow[i])
     2815      coinModel_.setElement(rowNumber,i,linearRow[i]);
     2816  }
     2817}
     2818// Gets correct form for a quadratic row - user to delete
     2819CoinPackedMatrix *
     2820OsiSolverLink::quadraticRow(int rowNumber,double * linearRow) const
     2821{
     2822  int numberColumns = coinModel_.numberColumns();
     2823  int numberRows = coinModel_.numberRows();
     2824  CoinZeroN(linearRow,numberColumns);
     2825  int numberElements=0;
     2826  assert (rowNumber>=0&&rowNumber<numberRows);
     2827  CoinModelLink triple=coinModel_.firstInRow(rowNumber);
     2828  while (triple.column()>=0) {
     2829    int iColumn = triple.column();
     2830    const char * expr = coinModel_.getElementAsString(rowNumber,iColumn);
     2831    if (strcmp(expr,"Numeric")) {
     2832      // try and see which columns
     2833      assert (strlen(expr)<20000);
     2834      char temp[20000];
     2835      strcpy(temp,expr);
     2836      char * pos = temp;
     2837      bool ifFirst=true;
     2838      while (*pos) {
     2839        double value;
     2840        int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
     2841        // must be column unless first when may be linear term
     2842        if (jColumn>=0) {
     2843          numberElements++;
     2844        } else if (jColumn==-2) {
     2845          linearRow[iColumn]=value;
     2846        } else {
     2847          printf("bad nonlinear term %s\n",temp);
     2848          abort();
     2849        }
     2850        ifFirst=false;
     2851      }
     2852    } else {
     2853      linearRow[iColumn]=coinModel_.getElement(rowNumber,iColumn);
     2854    }
     2855    triple=coinModel_.next(triple);
     2856  }
     2857  if (!numberElements) {
     2858    return NULL;
     2859  } else {
     2860    int * column = new int[numberElements];
     2861    int * column2 = new int[numberElements];
     2862    double * element = new double[numberElements];
     2863    numberElements=0;
     2864    CoinModelLink triple=coinModel_.firstInRow(rowNumber);
     2865    while (triple.column()>=0) {
     2866      int iColumn = triple.column();
     2867      const char * expr = coinModel_.getElementAsString(rowNumber,iColumn);
     2868      if (strcmp(expr,"Numeric")) {
     2869        // try and see which columns
     2870        assert (strlen(expr)<20000);
     2871        char temp[20000];
     2872        strcpy(temp,expr);
     2873        char * pos = temp;
     2874        bool ifFirst=true;
     2875        while (*pos) {
     2876          double value;
     2877          int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
     2878          // must be column unless first when may be linear term
     2879          if (jColumn>=0) {
     2880            column[numberElements]=iColumn;
     2881            column2[numberElements]=jColumn;
     2882            element[numberElements++]=value;
     2883          } else if (jColumn!=-2) {
     2884            printf("bad nonlinear term %s\n",temp);
     2885            abort();
     2886          }
     2887          ifFirst=false;
     2888        }
     2889      }
     2890      triple=coinModel_.next(triple);
     2891    }
     2892    return new CoinPackedMatrix(true,column2,column,element,numberElements);
     2893  }
     2894}
     2895/*
     2896  Problem specific
     2897  Returns -1 if node fathomed and no solution
     2898  0 if did nothing
     2899  1 if node fathomed and solution
     2900  allFixed is true if all LinkedBound variables are fixed
     2901*/
     2902int
     2903OsiSolverLink::fathom(bool allFixed)
     2904{
     2905  int returnCode=0;
     2906  if (allFixed) {
     2907    // all fixed so we can reformulate
     2908    OsiClpSolverInterface newSolver;
     2909    // set values
     2910    const double * lower = modelPtr_->columnLower();
     2911    const double * upper = modelPtr_->columnUpper();
     2912    int i;
     2913    for (i=0;i<numberFix_;i++ ) {
     2914      int iColumn = fixVariables_[i];
     2915      double lo = lower[iColumn];
     2916      double up = upper[iColumn];
     2917      assert (lo==up);
     2918      //printf("column %d fixed to %g\n",iColumn,lo);
     2919      coinModel_.associateElement(coinModel_.columnName(iColumn),lo);
     2920    }
     2921    newSolver.loadFromCoinModel(coinModel_,true);
     2922    for (i=0;i<numberFix_;i++ ) {
     2923      int iColumn = fixVariables_[i];
     2924      newSolver.setColLower(iColumn,lower[iColumn]);
     2925      newSolver.setColUpper(iColumn,lower[iColumn]);
     2926    }
     2927    // see if everything with objective fixed
     2928    const double * objective = modelPtr_->objective();
     2929    int numberColumns = newSolver.getNumCols();
     2930    bool zeroObjective=true;
     2931    double sum=0.0;
     2932    for (i=0;i<numberColumns;i++) {
     2933      if (upper[i]>lower[i]&&objective[i]) {
     2934        zeroObjective=false;
     2935        break;
     2936      } else {
     2937        sum += lower[i]*objective[i];
     2938      }
     2939    }
     2940    //if (fabs(sum-8.3)<1.0e-5)
     2941    //printf("possible\n");
     2942    if (zeroObjective) {
     2943      // randomize objective
     2944      ClpSimplex * clpModel = newSolver.getModelPtr();
     2945      const double * element = clpModel->matrix()->getMutableElements();
     2946      //const int * row = clpModel->matrix()->getIndices();
     2947      const CoinBigIndex * columnStart = clpModel->matrix()->getVectorStarts();
     2948      const int * columnLength = clpModel->matrix()->getVectorLengths();
     2949      double * objective = clpModel->objective();
     2950      for (i=0;i<numberColumns;i++) {
     2951        if (clpModel->isInteger(i)) {
     2952          double value=0.0;
     2953          for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     2954            value += fabs(element[j]);
     2955          }
     2956          objective[i]=value;
     2957        }
     2958      }
     2959    }
     2960    newSolver.writeMps("xx");
     2961    CbcModel model(newSolver);
     2962    // Now do requested saves and modifications
     2963    CbcModel * cbcModel = & model;
     2964    OsiSolverInterface * osiModel = model.solver();
     2965    OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
     2966    ClpSimplex * clpModel = osiclpModel->getModelPtr();
     2967    CglProbing probing;
     2968    probing.setMaxProbe(10);
     2969    probing.setMaxLook(10);
     2970    probing.setMaxElements(200);
     2971    probing.setMaxProbeRoot(50);
     2972    probing.setMaxLookRoot(10);
     2973    probing.setRowCuts(3);
     2974    probing.setRowCuts(0);
     2975    probing.setUsingObjective(true);
     2976    cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
     2977   
     2978    CglGomory gomory;
     2979    gomory.setLimitAtRoot(512);
     2980    cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
     2981   
     2982    CglKnapsackCover knapsackCover;
     2983    cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
     2984 
     2985    CglClique clique;
     2986    clique.setStarCliqueReport(false);
     2987    clique.setRowCliqueReport(false);
     2988    clique.setMinViolation(0.1);
     2989    cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
     2990    CglMixedIntegerRounding2 mixedIntegerRounding2;
     2991    cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
     2992   
     2993    CglFlowCover flowCover;
     2994    cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
     2995   
     2996    CglTwomir twomir;
     2997    twomir.setMaxElements(250);
     2998    cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
     2999    cbcModel->cutGenerator(6)->setTiming(true);
     3000   
     3001    CbcHeuristicFPump heuristicFPump(*cbcModel);
     3002    heuristicFPump.setWhen(1);
     3003    heuristicFPump.setMaximumPasses(20);
     3004    heuristicFPump.setDefaultRounding(0.5);
     3005    cbcModel->addHeuristic(&heuristicFPump);
     3006   
     3007    CbcRounding rounding(*cbcModel);
     3008    cbcModel->addHeuristic(&rounding);
     3009   
     3010    CbcHeuristicLocal heuristicLocal(*cbcModel);
     3011    heuristicLocal.setSearchType(1);
     3012    cbcModel->addHeuristic(&heuristicLocal);
     3013   
     3014    CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
     3015    cbcModel->addHeuristic(&heuristicGreedyCover);
     3016   
     3017    CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
     3018    cbcModel->addHeuristic(&heuristicGreedyEquality);
     3019   
     3020    CbcCompareDefault compare;
     3021    cbcModel->setNodeComparison(compare);
     3022    cbcModel->setNumberBeforeTrust(5);
     3023    cbcModel->setSpecialOptions(2);
     3024    cbcModel->messageHandler()->setLogLevel(1);
     3025    cbcModel->setMaximumCutPassesAtRoot(-100);
     3026    cbcModel->setMaximumCutPasses(1);
     3027    cbcModel->setMinimumDrop(0.05);
     3028    clpModel->setNumberIterations(1);
     3029    // For branchAndBound this may help
     3030    clpModel->defaultFactorizationFrequency();
     3031    clpModel->setDualBound(6.71523e+07);
     3032    clpModel->setPerturbation(50);
     3033    osiclpModel->setSpecialOptions(193);
     3034    osiclpModel->messageHandler()->setLogLevel(0);
     3035    osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
     3036    osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
     3037    // You can save some time by switching off message building
     3038    // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
     3039    // Solve
     3040   
     3041    cbcModel->initialSolve();
     3042    //double cutoff = model_->getCutoff();
     3043    cbcModel->setCutoff(1.0e50);
     3044    // to change exits
     3045    bool isFeasible=false;
     3046    int saveLogLevel=clpModel->logLevel();
     3047    clpModel->setLogLevel(0);
     3048    if (clpModel->tightenPrimalBounds()!=0) {
     3049      clpModel->setLogLevel(saveLogLevel);
     3050      returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
     3051    } else {
     3052      clpModel->setLogLevel(saveLogLevel);
     3053      clpModel->dual();  // clean up
     3054      // compute some things using problem size
     3055      cbcModel->setMinimumDrop(min(5.0e-2,
     3056                                   fabs(cbcModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
     3057      if (cbcModel->getNumCols()<500)
     3058        cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
     3059      else if (cbcModel->getNumCols()<5000)
     3060        cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop
     3061      else
     3062        cbcModel->setMaximumCutPassesAtRoot(20);
     3063      cbcModel->setMaximumCutPasses(1);
     3064      // Hand coded preprocessing
     3065      CglPreProcess process;
     3066      OsiSolverInterface * saveSolver=cbcModel->solver()->clone();
     3067      // Tell solver we are in Branch and Cut
     3068      saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
     3069      // Default set of cut generators
     3070      CglProbing generator1;
     3071      generator1.setUsingObjective(true);
     3072      generator1.setMaxPass(3);
     3073      generator1.setMaxProbeRoot(saveSolver->getNumCols());
     3074      generator1.setMaxElements(100);
     3075      generator1.setMaxLookRoot(50);
     3076      generator1.setRowCuts(3);
     3077      // Add in generators
     3078      process.addCutGenerator(&generator1);
     3079      process.messageHandler()->setLogLevel(cbcModel->logLevel());
     3080      OsiSolverInterface * solver2 =
     3081        process.preProcessNonDefault(*saveSolver,0,10);
     3082      // Tell solver we are not in Branch and Cut
     3083      saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
     3084      if (solver2)
     3085        solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
     3086      if (!solver2) {
     3087        std::cout<<"Pre-processing says infeasible!"<<std::endl;
     3088        delete saveSolver;
     3089        returnCode=-1;
     3090      } else {
     3091        std::cout<<"processed model has "<<solver2->getNumRows()
     3092                 <<" rows, "<<solver2->getNumCols()
     3093                 <<" and "<<solver2->getNumElements()<<std::endl;
     3094        // we have to keep solver2 so pass clone
     3095        solver2 = solver2->clone();
     3096        //solver2->writeMps("intmodel");
     3097        cbcModel->assignSolver(solver2);
     3098        cbcModel->initialSolve();
     3099        if (zeroObjective) {
     3100          cbcModel->setMaximumSolutions(1); // just getting a solution
     3101#if 0
     3102          OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (cbcModel->solver());
     3103          ClpSimplex * clpModel = osiclpModel->getModelPtr();
     3104          const double * element = clpModel->matrix()->getMutableElements();
     3105          //const int * row = clpModel->matrix()->getIndices();
     3106          const CoinBigIndex * columnStart = clpModel->matrix()->getVectorStarts();
     3107          const int * columnLength = clpModel->matrix()->getVectorLengths();
     3108          int n=clpModel->numberColumns();
     3109          int * sort2 = new int[n];
     3110          int * pri = new int[n];
     3111          double * sort = new double[n];
     3112          int i;
     3113          int nint=0;
     3114          for (i=0;i<n;i++) {
     3115            if (clpModel->isInteger(i)) {
     3116              double largest=0.0;
     3117              for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     3118                largest = CoinMax(largest,fabs(element[j]));
     3119              }
     3120              sort2[nint]=nint;
     3121              sort[nint++]=-largest;
     3122            }
     3123          }
     3124          CoinSort_2(sort,sort+nint,sort2);
     3125          int kpri=1;
     3126          double last = sort[0];
     3127          for (i=0;i<nint;i++) {
     3128            if (sort[i]!=last) {
     3129              kpri++;
     3130              last=sort[i];
     3131            }
     3132            pri[sort2[i]]=kpri;
     3133          }
     3134          cbcModel->passInPriorities(pri,false);
     3135          delete [] sort;
     3136          delete [] sort2;
     3137          delete [] pri;
     3138#endif
     3139        }
     3140        cbcModel->branchAndBound();
     3141        // For best solution
     3142        int numberColumns = newSolver.getNumCols();
     3143        if (cbcModel->getMinimizationObjValue()<1.0e50) {
     3144          // post process
     3145          process.postProcess(*cbcModel->solver());
     3146          // Solution now back in saveSolver
     3147          cbcModel->assignSolver(saveSolver);
     3148          memcpy(cbcModel->bestSolution(),cbcModel->solver()->getColSolution(),
     3149                 numberColumns*sizeof(double));
     3150          // put back in original solver
     3151          newSolver.setColSolution(cbcModel->bestSolution());
     3152          isFeasible=true;
     3153        } else {
     3154          delete saveSolver;
     3155        }
     3156      }
     3157      //const double * solution = newSolver.getColSolution();
     3158      if (isFeasible&&cbcModel->getMinimizationObjValue()<1.0e50) {
     3159        int numberColumns = this->getNumCols();
     3160        int i;
     3161        const double * solution = cbcModel->bestSolution();
     3162        int numberColumns2 = newSolver.getNumCols();
     3163        for (i=0;i<numberColumns2;i++) {
     3164          double value = solution[i];
     3165          assert (fabs(value-floor(value+0.5))<0.0001);
     3166          value = floor(value+0.5);
     3167          this->setColLower(i,value);
     3168          this->setColUpper(i,value);
     3169        }
     3170        for (;i<numberColumns;i++) {
     3171          this->setColLower(i,0.0);
     3172          this->setColUpper(i,1.1);
     3173        }
     3174        // but take off cuts
     3175        int numberRows = getNumRows();
     3176        int numberRows2 = cbcModel_->continuousSolver()->getNumRows();
     3177           
     3178        for (i=numberRows2;i<numberRows;i++)
     3179          setRowBounds(i,-COIN_DBL_MAX,COIN_DBL_MAX);
     3180        initialSolve();
     3181        if (!isProvenOptimal())
     3182          getModelPtr()->writeMps("bad.mps");
     3183        if (isProvenOptimal()) {
     3184          delete [] bestSolution_;
     3185          bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols());
     3186          bestObjectiveValue_ = modelPtr_->objectiveValue();
     3187          printf("BB best value %g\n",bestObjectiveValue_);
     3188          returnCode = 1;
     3189        } else {
     3190          printf("*** WHY BAD SOL\n");
     3191          returnCode=-1;
     3192        }
     3193      } else {
     3194        modelPtr_->setProblemStatus(1);
     3195        modelPtr_->setObjectiveValue(COIN_DBL_MAX);
     3196        returnCode = -1;
     3197      }
     3198    }
     3199  }
     3200  return returnCode;
    23983201}
    23993202//#############################################################################
     
    52316034  return *this;
    52326035}
     6036void checkQP(ClpSimplex * model)
     6037{
     6038#if 0
     6039  printf("Checking quadratic model %x\n",model);
     6040  if (model) {
     6041    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
     6042    assert (quadraticObj);
     6043    CoinPackedMatrix * quadraticObjective = quadraticObj->quadraticObjective();
     6044    int numberColumns = quadraticObj->numberColumns();
     6045    const int * columnQuadratic = quadraticObjective->getIndices();
     6046    const CoinBigIndex * columnQuadraticStart = quadraticObjective->getVectorStarts();
     6047    const int * columnQuadraticLength = quadraticObjective->getVectorLengths();
     6048    //const double * quadraticElement = quadraticObjective->getElements();
     6049    for (int i=0;i<numberColumns;i++) {
     6050      for (int j=columnQuadraticStart[i];j<columnQuadraticStart[i]+columnQuadraticLength[i];j++)
     6051        assert (columnQuadratic[j]>=0&&columnQuadratic[j]<1000);
     6052    }
     6053  }
    52336054#endif
     6055}
     6056//#############################################################################
     6057// Solve methods
     6058//#############################################################################
     6059void OsiSolverLinearizedQuadratic::initialSolve()
     6060{
     6061  OsiClpSolverInterface::initialSolve();
     6062  int secondaryStatus = modelPtr_->secondaryStatus();
     6063  if (modelPtr_->status()==0&&(secondaryStatus==2||secondaryStatus==4))
     6064    modelPtr_->cleanup(1);
     6065  if (isProvenOptimal()&&modelPtr_->numberColumns()==quadraticModel_->numberColumns()) {
     6066    // see if qp can get better solution
     6067    const double * solution = modelPtr_->primalColumnSolution();
     6068    int numberColumns = modelPtr_->numberColumns();
     6069    bool satisfied=true;
     6070    for (int i=0;i<numberColumns;i++) {
     6071      if (isInteger(i)) {
     6072        double value = solution[i];
     6073        if (fabs(value-floor(value+0.5))>1.0e-6) {
     6074          satisfied=false;
     6075          break;
     6076        }
     6077      }
     6078    }
     6079    if (satisfied) {
     6080      checkQP(quadraticModel_);
     6081      ClpSimplex qpTemp(*quadraticModel_);
     6082      checkQP(&qpTemp);
     6083      double * lower = qpTemp.columnLower();
     6084      double * upper = qpTemp.columnUpper();
     6085      double * lower2 = modelPtr_->columnLower();
     6086      double * upper2 = modelPtr_->columnUpper();
     6087      for (int i=0;i<numberColumns;i++) {
     6088        if (isInteger(i)) {
     6089          double value = floor(solution[i]+0.5);
     6090          lower[i]=value;
     6091          upper[i]=value;
     6092        } else {
     6093          lower[i]=lower2[i];
     6094          upper[i]=upper2[i];
     6095        }
     6096      }
     6097      //qpTemp.writeMps("bad.mps");
     6098      //modelPtr_->writeMps("bad2.mps");
     6099      //qpTemp.objectiveAsObject()->setActivated(0);
     6100      //qpTemp.primal();
     6101      //qpTemp.objectiveAsObject()->setActivated(1);
     6102      qpTemp.primal();
     6103      //assert (!qpTemp.problemStatus());
     6104      if (qpTemp.objectiveValue()<bestObjectiveValue_&&!qpTemp.problemStatus()) {
     6105        delete [] bestSolution_;
     6106        bestSolution_ = CoinCopyOfArray(qpTemp.primalColumnSolution(),numberColumns);
     6107        bestObjectiveValue_ = qpTemp.objectiveValue();
     6108        printf("better qp objective of %g\n",bestObjectiveValue_);
     6109      }
     6110    }
     6111  }
     6112}
     6113//#############################################################################
     6114// Constructors, destructors clone and assignment
     6115//#############################################################################
     6116//-------------------------------------------------------------------
     6117// Default Constructor
     6118//-------------------------------------------------------------------
     6119OsiSolverLinearizedQuadratic::OsiSolverLinearizedQuadratic ()
     6120  : OsiClpSolverInterface()
     6121{
     6122  bestObjectiveValue_=COIN_DBL_MAX;
     6123  bestSolution_=NULL;
     6124  specialOptions3_=0;
     6125  quadraticModel_=NULL;
     6126}
     6127OsiSolverLinearizedQuadratic::OsiSolverLinearizedQuadratic ( ClpSimplex * quadraticModel)
     6128  : OsiClpSolverInterface(new ClpSimplex(*quadraticModel),true)
     6129{
     6130  bestObjectiveValue_=COIN_DBL_MAX;
     6131  bestSolution_=NULL;
     6132  specialOptions3_=0;
     6133  quadraticModel_=new ClpSimplex(*quadraticModel);
     6134  // linearize
     6135  int numberColumns = modelPtr_->numberColumns();
     6136  const double * solution = modelPtr_->primalColumnSolution();
     6137  // Replace objective
     6138  ClpObjective * trueObjective = modelPtr_->objectiveAsObject();
     6139  ClpObjective * objective=new ClpLinearObjective(NULL,numberColumns);
     6140  modelPtr_->setObjectivePointer(objective);
     6141  double offset;
     6142  double saveOffset = modelPtr_->objectiveOffset();
     6143  memcpy(modelPtr_->objective(),trueObjective->gradient(modelPtr_,solution,offset,true,2),
     6144         numberColumns*sizeof(double));
     6145  modelPtr_->setObjectiveOffset(saveOffset+offset);
     6146  delete trueObjective;
     6147  checkQP(quadraticModel_);
     6148}
     6149//-------------------------------------------------------------------
     6150// Clone
     6151//-------------------------------------------------------------------
     6152OsiSolverInterface *
     6153OsiSolverLinearizedQuadratic::clone(bool copyData) const
     6154{
     6155  assert (copyData);
     6156  return new OsiSolverLinearizedQuadratic(*this);
     6157}
     6158
     6159
     6160//-------------------------------------------------------------------
     6161// Copy constructor
     6162//-------------------------------------------------------------------
     6163OsiSolverLinearizedQuadratic::OsiSolverLinearizedQuadratic (
     6164                  const OsiSolverLinearizedQuadratic & rhs)
     6165  : OsiClpSolverInterface(rhs)
     6166{
     6167  bestObjectiveValue_=rhs.bestObjectiveValue_;
     6168  if (rhs.bestSolution_) {
     6169    bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->numberColumns());
     6170  } else {
     6171    bestSolution_=NULL;
     6172  }
     6173  specialOptions3_=rhs.specialOptions3_;
     6174  if (rhs.quadraticModel_) {
     6175    quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
     6176  } else {
     6177    quadraticModel_=NULL;
     6178  }
     6179  checkQP(rhs.quadraticModel_);
     6180  checkQP(quadraticModel_);
     6181}
     6182
     6183//-------------------------------------------------------------------
     6184// Destructor
     6185//-------------------------------------------------------------------
     6186OsiSolverLinearizedQuadratic::~OsiSolverLinearizedQuadratic ()
     6187{
     6188  delete [] bestSolution_;
     6189  delete quadraticModel_;
     6190}
     6191
     6192//-------------------------------------------------------------------
     6193// Assignment operator
     6194//-------------------------------------------------------------------
     6195OsiSolverLinearizedQuadratic &
     6196OsiSolverLinearizedQuadratic::operator=(const OsiSolverLinearizedQuadratic& rhs)
     6197{
     6198  if (this != &rhs) {
     6199    delete [] bestSolution_;
     6200    delete quadraticModel_;
     6201    OsiClpSolverInterface::operator=(rhs);
     6202    bestObjectiveValue_=rhs.bestObjectiveValue_;
     6203    if (rhs.bestSolution_) {
     6204      bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->numberColumns());
     6205    } else {
     6206      bestSolution_=NULL;
     6207    }
     6208    specialOptions3_=rhs.specialOptions3_;
     6209    if (rhs.quadraticModel_) {
     6210      quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
     6211    } else {
     6212      quadraticModel_=NULL;
     6213    }
     6214    checkQP(rhs.quadraticModel_);
     6215    checkQP(quadraticModel_);
     6216  }
     6217  return *this;
     6218}
     6219#endif
  • branches/devel/Cbc/src/CbcLinked.hpp

    r520 r523  
    1818class OsiLinkedBound;
    1919class OsiObject;
     20class CglStored;
    2021//#############################################################################
    2122
     
    4546     allFixed is true if all LinkedBound variables are fixed
    4647  */
    47   virtual int fathom(bool allFixed) {return 0;};
     48  virtual int fathom(bool allFixed) ;
    4849  /** Solves nonlinear problem from CoinModel using SLP - may be used as crash
    4950      for other algorithms when number of iterations small.
     
    5354  */
    5455  double * nonlinearSLP(int numberPasses,double deltaTolerance);
     56  /** Solve linearized quadratic objective branch and bound.
     57      Return cutoff and OA cut
     58  */
     59  double linearizedBAB(CglStored * cut) ;
    5560  //@}
    5661 
     
    7580  OsiSolverLink(  CoinModel & modelObject);
    7681  // Other way with existing object
    77   void load(  CoinModel & modelObject,bool tightenBounds=false);
     82  void load(  CoinModel & modelObject,bool tightenBounds=false,int logLevel=1);
    7883  /// Clone
    7984  virtual OsiSolverInterface * clone(bool copyData=true) const;
     
    134139  ClpSimplex * quadraticModel() const
    135140  { return quadraticModel_;};
     141  /// Gets correct form for a quadratic row - user to delete
     142  CoinPackedMatrix * quadraticRow(int rowNumber,double * linear) const;
     143  /// Replaces a quadratic row
     144  void replaceQuadraticRow(int rowNumber,const double * linear, const CoinPackedMatrix * quadraticPart);
    136145  /// Default meshSize
    137146  inline double defaultMeshSize() const
     
    162171  inline const CoinModel * coinModel() const
    163172  { return &coinModel_;};
     173  /// Set all biLinear priorities on x-x variables
     174  void setBiLinearPriorities(int value);
     175  /// Set all mesh sizes on x-x variables
     176  void setMeshSizes(double value);
     177  /** Two tier integer problem where when set of variables with priority
     178      less than this are fixed the problem becomes an easier integer problem
     179  */
     180  void setFixedPriority(int priorityValue);
    164181  //@}
    165182 
     
    212229  OsiLinkedBound * info_;
    213230  /**
    214      0 bit (1) - don't do mini B&B
     231     0 bit (1) - call fathom (may do mini B&B)
    215232     1 bit (2) - quadratic only in objective (add OA cuts)
    216233     2 bit (4) - convex
     
    234251  /// Priority for bilinear
    235252  int biLinearPriority_;
     253  /// Number of variables which when fixed help
     254  int numberFix_;
     255  /// list of fixed variables
     256  int * fixVariables_;
    236257  //@}
    237258};
     
    10271048  // Private member data
    10281049};
     1050//#############################################################################
     1051
     1052/**
     1053   
     1054This is to allow the user to replace initialSolve and resolve
     1055*/
     1056
     1057class OsiSolverLinearizedQuadratic : public OsiClpSolverInterface {
     1058 
     1059public:
     1060  //---------------------------------------------------------------------------
     1061  /**@name Solve methods */
     1062  //@{
     1063  /// Solve initial LP relaxation
     1064  virtual void initialSolve();
     1065  //@}
     1066 
     1067 
     1068  /**@name Constructors and destructors */
     1069  //@{
     1070  /// Default Constructor
     1071  OsiSolverLinearizedQuadratic ();
     1072  /// Useful constructor (solution should be good)
     1073  OsiSolverLinearizedQuadratic(  ClpSimplex * quadraticModel);
     1074  /// Clone
     1075  virtual OsiSolverInterface * clone(bool copyData=true) const;
     1076 
     1077  /// Copy constructor
     1078  OsiSolverLinearizedQuadratic (const OsiSolverLinearizedQuadratic &);
     1079 
     1080  /// Assignment operator
     1081  OsiSolverLinearizedQuadratic & operator=(const OsiSolverLinearizedQuadratic& rhs);
     1082 
     1083  /// Destructor
     1084  virtual ~OsiSolverLinearizedQuadratic ();
     1085 
     1086  //@}
     1087 
     1088 
     1089  /**@name Sets and Gets */
     1090  //@{
     1091  /// Objective value of best solution found internally
     1092  inline double bestObjectiveValue() const
     1093  { return bestObjectiveValue_;};
     1094  /// Best solution found internally
     1095  const double * bestSolution() const
     1096  { return bestSolution_;};
     1097  /// Set special options
     1098  inline void setSpecialOptions3(int value)
     1099  { specialOptions3_=value;};
     1100  /// Get special options
     1101  inline int specialOptions3() const
     1102  { return specialOptions3_;};
     1103  /// Copy of quadratic model if one
     1104  ClpSimplex * quadraticModel() const
     1105  { return quadraticModel_;};
     1106  //@}
     1107 
     1108  //---------------------------------------------------------------------------
     1109 
     1110protected:
     1111 
     1112 
     1113  /**@name functions */
     1114  //@{
     1115 
     1116  /**@name Private member data */
     1117  //@{
     1118  /// Objective value of best solution found internally
     1119  double bestObjectiveValue_;
     1120  /// Copy of quadratic model if one
     1121  ClpSimplex * quadraticModel_;
     1122  /// Best solution found internally
     1123  double * bestSolution_;
     1124  /**
     1125     0 bit (1) - don't do mini B&B
     1126     1 bit (2) - quadratic only in objective
     1127  */
     1128  int specialOptions3_;
     1129  //@}
     1130};
    10291131#endif
    10301132#endif
  • branches/devel/Cbc/src/Cbc_ampl.cpp

    r520 r523  
    292292    fileName[0]='\0';
    293293  int saveArgc = argc;
    294   memset(info,0,sizeof(ampl_info));
     294  if (info->numberRows!=-1234567)
     295    memset(info,0,sizeof(ampl_info)); // overwrite unless magic number set
    295296  /* save so can be accessed by decodePhrase */
    296297  saveInfo = info;
     
    449450    if (!strstr(fileName,".nl"))
    450451      strcat(fileName,".nl");
    451     CoinModel * model = new CoinModel(1,fileName);
     452    CoinModel * model = new CoinModel(1,fileName,info);
    452453    if (model->numberRows()>0)
    453454      *coinModel=(void *) model;
     
    466467    info->numberBinary=nbv;
    467468    info->numberIntegers=niv;
    468     if (niv+nbv>0)
     469    if (nbv+niv+nlvbi+nlvci+nlvoi>0)
    469470      mip_stuff(); // get any extra info
    470471  }
     
    619620/* Read a problem from AMPL nl file
    620621 */
    621 CoinModel::CoinModel( int nonLinear, const char * fileName)
     622CoinModel::CoinModel( int nonLinear, const char * fileName,const void * info)
    622623 :  numberRows_(0),
    623624    maximumRows_(0),
     
    671672  }
    672673  if (!status) {
    673     gdb(nonLinear,fileName);
    674   }
    675 }
    676  static real
     674    gdb(nonLinear,fileName,info);
     675  }
     676}
     677#if 0
     678static real
    677679qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
    678680{
     
    692694  return 0.5 * t;
    693695}
    694 
     696#endif
    695697void
    696 CoinModel::gdb( int nonLinear, const char * fileName)
    697 {
     698CoinModel::gdb( int nonLinear, const char * fileName,const void * info)
     699{
     700  const ampl_info * amplInfo = (const ampl_info *) info;
    698701  ograd *og=NULL;
    699702  int i;
     
    10221025              double value = element[k];
    10231026              // ampl gives twice with assumed 0.5
    1024               if (kColumn>j)
     1027              if (kColumn<j)
    10251028                continue;
    10261029              else if (kColumn==j)
     
    10561059              assert (strlen(temp)<1000);
    10571060              setElement(iRow,j,temp);
    1058               printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
     1061              if (amplInfo->logLevel>1)
     1062                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
    10591063            }
    10601064          }
     
    10661070              double value = element[k];
    10671071              // ampl gives twice with assumed 0.5
    1068               if (kColumn>j)
     1072              if (kColumn<j)
    10691073                continue;
    10701074              else if (kColumn==j)
     
    11001104              assert (strlen(temp)<1000);
    11011105              setObjective(j,temp);
    1102               printf("el for objective column c%7.7d is %s\n",j,temp);
     1106              if (amplInfo->logLevel>1)
     1107                printf("el for objective column c%7.7d is %s\n",j,temp);
    11031108            }
    11041109          }
  • branches/devel/Cbc/src/Cbc_ampl.h

    r496 r523  
    3939  char ** arguments;
    4040  char buffer[300];
     41  int logLevel;
    4142} ampl_info;
    4243#ifdef __cplusplus
  • branches/devel/Cbc/src/CoinSolve.cpp

    r521 r523  
    613613    if (argc>2&&!strcmp(argv[2],"-AMPL")) {
    614614      usingAmpl=true;
     615      // see if log in list
     616      noPrinting=true;
     617      for (int i=1;i<argc;i++) {
     618        if (!strncmp(argv[i],"log",3)) {
     619          char * equals = strchr(argv[i],'=');
     620          if (equals&&atoi(equals+1)>0) {
     621            noPrinting=false;
     622            info.logLevel=atoi(equals+1);
     623            // mark so won't be overWritten
     624            info.numberRows=-1234567;
     625            break;
     626          }
     627        }
     628      }
    615629      int returnCode = readAmpl(&info,argc,const_cast<char **>(argv),(void **) (& coinModel));
    616630      if (returnCode)
    617631        return returnCode;
    618632      CbcOrClpRead_mode=2; // so will start with parameters
    619       // see if log in list
    620       noPrinting=true;
     633      // see if log in list (including environment)
    621634      for (int i=1;i<info.numberArguments;i++) {
    622635        if (!strcmp(info.arguments[i],"log")) {
     
    653666        si->setDefaultBound(100.0);
    654667        si->setDefaultMeshSize(0.01);
    655         si->setDefaultBound(100.0);
     668        si->setDefaultBound(1000.0);
    656669        si->setIntegerPriority(1000);
    657670        si->setBiLinearPriority(10000);
    658671        CoinModel * model2 = (CoinModel *) coinModel;
    659         si->load(*model2);
     672        si->load(*model2,true,info.logLevel);
    660673        // redo
    661674        solver = model.solver();
     
    11201133              case TIGHTENFACTOR:
    11211134                tightenFactor=value;
    1122                 defaultSettings=false; // user knows what she is doing
     1135                if(!complicatedInteger)
     1136                  defaultSettings=false; // user knows what she is doing
    11231137                break;
    11241138              default:
     
    18811895                //}
    18821896                // bounds based on continuous
    1883                 if (tightenFactor) {
     1897                if (tightenFactor&&!complicatedInteger) {
    18841898                  if (modelC->tightenPrimalBounds(tightenFactor)!=0) {
    18851899                    std::cout<<"Problem is infeasible!"<<std::endl;
     
    20892103                modelC->dual();
    20902104              }
     2105#if 0
     2106              numberDebugValues=599;
     2107              debugValues = new double[numberDebugValues];
     2108              CoinZeroN(debugValues,numberDebugValues);
     2109              debugValues[3]=1.0;
     2110              debugValues[6]=25.0;
     2111              debugValues[9]=4.0;
     2112              debugValues[26]=4.0;
     2113              debugValues[27]=6.0;
     2114              debugValues[35]=8.0;
     2115              debugValues[53]=21.0;
     2116              debugValues[56]=4.0;
     2117#endif
    20912118              if (debugValues) {
    20922119                // for debug
     
    23562383              // Turn this off if you get problems
    23572384              // Used to be automatically set
    2358               int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue();
     2385              int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()%10000;
    23592386              if (mipOptions!=(128|64|1))
    23602387                printf("mip options %d\n",mipOptions);
     
    28762903                    OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
    28772904                    if (solver3) {
     2905                      if (tightenFactor>0.0) {
     2906                        // set grid size for all continuous bi-linear
     2907                        solver3->setMeshSizes(tightenFactor);
     2908                      }
     2909                      int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000;
     2910                      CglStored stored;
     2911                      if (options) {
     2912                        printf("nlp options %d\n",options);
     2913                        /*
     2914                          1 - force mini branch and bound
     2915                          2 - set priorities high on continuous
     2916                          4 - try adding OA cuts
     2917                          8 - try doing quadratic linearization
     2918                        */
     2919                        if ((options&2))
     2920                          solver3->setBiLinearPriorities(10);
     2921                        if ((options&4))
     2922                          solver3->setSpecialOptions2(solver3->specialOptions2()|(8+4));
     2923                        if ((options&1)!=0&&slpValue>0)
     2924                          solver3->setFixedPriority(slpValue);
     2925                        double cutoff=COIN_DBL_MAX;
     2926                        if ((options&4))
     2927                          cutoff=solver3->linearizedBAB(&stored);
     2928                        if (cutoff<babModel->getCutoff())
     2929                          babModel->setCutoff(cutoff);
     2930                      }
    28782931                      solver3->setCbcModel(babModel);
    2879                       CglStored stored;
    28802932                      babModel->addCutGenerator(&stored,1,"Stored");
    28812933                      CglTemporary temp;
Note: See TracChangeset for help on using the changeset viewer.