Changeset 1033 for branches


Ignore:
Timestamp:
Jun 25, 2007 10:21:30 AM (12 years ago)
Author:
forrest
Message:

update clp/devel

Location:
branches/devel/Clp/src
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Clp/src/CbcOrClpParam.cpp

    r1010 r1033  
    546546    case SPECIALOPTIONS:
    547547      model->setSpecialOptions(value);
     548#ifdef COIN_HAS_CBC
    548549    case THREADS:
    549550      model->setNumberThreads(value);
    550551      break;
     552#endif
    551553    default:
    552554      break;
     
    578580    value=model->specialOptions();
    579581    break;
     582#ifndef COIN_HAS_CBC
    580583  case THREADS:
    581584    value = model->numberThreads();
     585#endif
    582586  default:
    583587    value=intValue_;
     
    784788      model.setSizeMiniTree(value);
    785789      break;
     790    case CUTPASSINTREE:
     791      oldValue=model.getMaximumCutPasses();
     792      model.setMaximumCutPasses(value);
     793      break;
     794    case CUTPASS:
     795      oldValue=model.getMaximumCutPassesAtRoot();
     796      model.setMaximumCutPassesAtRoot(value);
     797      break;
     798#ifdef COIN_HAS_CBC
     799    case THREADS:
     800      oldValue=model.getNumberThreads();
     801      model.setNumberThreads(value);
     802      break;
     803#endif
    786804    default:
    787805      break;
     
    822840    value=model.sizeMiniTree();
    823841    break;
     842  case CUTPASSINTREE:
     843    value=model.getMaximumCutPasses();
     844    break;
     845  case CUTPASS:
     846    value=model.getMaximumCutPassesAtRoot();
     847    break;
     848#ifdef COIN_HAS_CBC
     849  case THREADS:
     850    value = model.getNumberThreads();
     851#endif
    824852  default:
    825853    value=intValue_;
     
    19511979 more lightweight or do more if improvements are still being made.  As Presolve will return\
    19521980 if nothing is being taken out, you should not normally need to use this fine tuning."
    1953      );
     1981     );
     1982#endif
     1983#ifdef COIN_HAS_CBC
     1984  parameters[numberParameters++]=
     1985    CbcOrClpParam("passT!reeCuts","Number of cut passes in tree",
     1986                  -999999,999999,CUTPASSINTREE);
     1987  parameters[numberParameters-1].setLonghelp
     1988    (
     1989     "The default is one pass"
     1990     );
     1991#endif
     1992#ifdef COIN_HAS_CLP
    19541993  parameters[numberParameters++]=
    19551994    CbcOrClpParam("pertV!alue","Method of perturbation",
     
    22092248See branchAndCut for information on options."
    22102249     );
     2250  parameters[numberParameters++]=
     2251    CbcOrClpParam("residual!CapacityCuts","Whether to use Residual Capacity cuts",
     2252                  "off",RESIDCUTS);
     2253  parameters[numberParameters-1].append("on");
     2254  parameters[numberParameters-1].append("root");
     2255  parameters[numberParameters-1].append("ifmove");
     2256  parameters[numberParameters-1].append("forceOn");
     2257  parameters[numberParameters-1].setLonghelp
     2258    (
     2259     "Residual capacity cuts. \
     2260See branchAndCut for information on options."
     2261     );
    22112262#endif
    22122263#ifdef COIN_HAS_CLP
     
    24462497                  -1,INT_MAX,TESTOSI,false);
    24472498#endif
    2448 #ifdef COIN_HAS_CLP
     2499#ifdef CBC_THREAD
    24492500  parameters[numberParameters++]=
    24502501    CbcOrClpParam("thread!s","Number of threads to try and use",
    2451                   -2,64,THREADS,false);
     2502                  -100,10000,THREADS,false);
     2503  parameters[numberParameters-1].setLonghelp
     2504    (
     2505     "To use multiple threads, set threads to number wanted.  It may be better \
     2506to use one or two more than number of cpus available.  If 100+n then n threads and \
     2507threads used in sub-trees, if 200+n use threads for root cuts, 300+n - both."
     2508     );
    24522509#endif
    24532510#ifdef COIN_HAS_CBC
  • branches/devel/Clp/src/CbcOrClpParam.hpp

    r918 r1033  
    6060    MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    6161    OUTPUTFORMAT,SLPVALUE,PRESOLVEOPTIONS,PRINTOPTIONS,SPECIALOPTIONS,
    62     SUBSTITUTION,DUALIZE,VERBOSE,THREADS,CPP,CUTPASS,PROCESSTUNE,
     62    SUBSTITUTION,DUALIZE,VERBOSE,CPP,PROCESSTUNE,
    6363
    6464    STRONGBRANCHING=151,CUTDEPTH, MAXNODES,NUMBERBEFORE,NUMBERANALYZE,
    6565    NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,MAXSOLS,
    66     FPUMPTUNE,TESTOSI,EXTRA1,EXTRA2,EXTRA3,EXTRA4,
     66    FPUMPTUNE,TESTOSI,EXTRA1,EXTRA2,EXTRA3,EXTRA4,CUTPASSINTREE,
     67    THREADS,CUTPASS,
    6768#ifdef COIN_HAS_CBC
    6869    LOGLEVEL ,
     
    7778    ROUNDING,SOLVER,CLIQUECUTS,COSTSTRATEGY,FLOWCUTS,MIXEDCUTS,
    7879    TWOMIRCUTS,PREPROCESS,FPUMP,GREEDY,COMBINE,LOCALTREE,USESOLUTION,SOS,
    79     LANDPCUTS,RINS,
     80    LANDPCUTS,RINS,RESIDCUTS,
    8081   
    8182    DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,EITHERSIMPLEX,
  • branches/devel/Clp/src/ClpConstraint.cpp

    r1005 r1033  
    6161  return *this;
    6262}
     63// Constraint function value
     64double
     65ClpConstraint::functionValue (const ClpSimplex * model,
     66                              const double * solution,
     67                              bool useScaling,
     68                              bool refresh) const
     69{
     70  double offset;
     71  double value;
     72  int n = model->numberColumns();
     73  double * grad = new double [n];
     74  gradient(model,solution,grad,value,offset,useScaling,refresh);
     75  delete [] grad;
     76  return value;
     77}
    6378
  • branches/devel/Clp/src/ClpConstraint.hpp

    r1005 r1033  
    2525      If refresh is false then uses last solution
    2626      Uses model for scaling
    27       Returns non-zero if gradient udefined at current solution
     27      Returns non-zero if gradient undefined at current solution
    2828  */
    2929  virtual int gradient(const ClpSimplex * model,
     
    3434                       bool useScaling=false,
    3535                       bool refresh=true) const =0;
     36  /// Constraint function value
     37  virtual double functionValue (const ClpSimplex * model,
     38                        const double * solution,
     39                        bool useScaling=false,
     40                        bool refresh=true) const ;
    3641  /// Resize constraint
    3742  virtual void resize(int newNumberColumns) = 0;
     
    8287  virtual int numberCoefficients() const = 0;
    8388 
    84   /// Constraint function value
     89  /// Stored constraint function value
    8590  inline double functionValue () const
    8691  { return functionValue_;};
     
    8994  inline double offset () const
    9095  { return offset_;};
     96  /// Say we have new primal solution - so may need to recompute
     97  virtual void newXValues() {};
    9198  //@}
    9299
  • branches/devel/Clp/src/ClpConstraintLinear.cpp

    r1005 r1033  
    134134{
    135135  if (numberColumns_!=newNumberColumns) {
     136#ifndef NDEBUG
    136137    int lastColumn = column_[numberCoefficients_-1];
     138#endif
    137139    assert (newNumberColumns>lastColumn);
    138140    delete [] lastGradient_;
  • branches/devel/Clp/src/ClpConstraintQuadratic.cpp

    r1007 r1033  
    189189  if (numberColumns_!=newNumberColumns) {
    190190    abort();
     191#ifndef NDEBUG
    191192    int lastColumn = column_[numberCoefficients_-1];
     193#endif
    192194    assert (newNumberColumns>lastColumn);
    193195    delete [] lastGradient_;
  • branches/devel/Clp/src/ClpLinearObjective.cpp

    r754 r1033  
    198198    return 0.0;
    199199  }
     200}
     201// Return objective value (without any ClpModel offset) (model may be NULL)
     202double
     203ClpLinearObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
     204{
     205  const double * cost = objective_;
     206  if (model&&model->costRegion())
     207    cost = model->costRegion();
     208  double currentObj=0.0;
     209  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     210    currentObj += cost[iColumn]*solution[iColumn];
     211  }
     212  return currentObj;
    200213}
    201214//-------------------------------------------------------------------
  • branches/devel/Clp/src/ClpLinearObjective.hpp

    r754 r1033  
    4545                            double & predictedObj,
    4646                            double & thetaObj);
     47  /// Return objective value (without any ClpModel offset) (model may be NULL)
     48  virtual double objectiveValue(const ClpSimplex * model, const double * solution) const ;
    4749  /// Resize objective
    4850  virtual void resize(int newNumberColumns) ;
  • branches/devel/Clp/src/ClpMessage.cpp

    r918 r1033  
    100100  {CLP_BAD_STRING_VALUES,3005,1,"%d string elements had no values associated with them"},
    101101  {CLP_CRUNCH_STATS,61,2,"Crunch %d (%d) rows, %d (%d) columns and %d (%d) elements"},
     102  {CLP_GENERAL,1000,1,"%s"},
    102103  {CLP_DUMMY_END,999999,0,""}
    103104};
  • branches/devel/Clp/src/ClpMessage.hpp

    r754 r1033  
    9999  CLP_BAD_STRING_VALUES,
    100100  CLP_CRUNCH_STATS,
     101  CLP_GENERAL,
    101102  CLP_DUMMY_END
    102103};
  • branches/devel/Clp/src/ClpModel.cpp

    r996 r1033  
    291291  gutsOfLoadModel(numrows, numcols,
    292292                  collb, colub, obj, rowlb, rowub, rowObjective);
    293   CoinPackedMatrix matrix(true,numrows,numcols,start[numcols],
     293  int numberElements = start ? start[numcols] : 0;
     294  CoinPackedMatrix matrix(true,numrows,numrows ? numcols : 0,numberElements,
    294295                              value,index,start,NULL);
    295296  matrix_ = new ClpPackedMatrix(matrix);
     
    322323ClpModel::loadProblem (  CoinModel & modelObject,bool tryPlusMinusOne)
    323324{
    324   if (modelObject.numberElements()==0)
     325  if (modelObject.numberColumns()==0&&modelObject.numberRows()==0)
    325326    return 0;
    326327  int numberErrors = 0;
     
    991992  CoinPackedMatrix matrix2;
    992993  matrix_=new ClpPackedMatrix(matrix2);
     994}
     995/* Really clean up matrix.
     996   a) eliminate all duplicate AND small elements in matrix
     997   b) remove all gaps and set extraGap_ and extraMajor_ to 0.0
     998   c) reallocate arrays and make max lengths equal to lengths
     999   d) orders elements
     1000   returns number of elements eliminated or -1 if not ClpMatrix
     1001*/
     1002int
     1003ClpModel::cleanMatrix(double threshold)
     1004{
     1005  ClpPackedMatrix * matrix = (dynamic_cast< ClpPackedMatrix*>(matrix_));
     1006  if (matrix) {
     1007    return matrix->getPackedMatrix()->cleanMatrix(threshold);
     1008  } else {
     1009    return -1;
     1010  }
    9931011}
    9941012// Resizes
  • branches/devel/Clp/src/ClpModel.hpp

    r989 r1033  
    249249  /// Create empty ClpPackedMatrix
    250250  void createEmptyMatrix();
     251  /** Really clean up matrix (if ClpPackedMatrix).
     252      a) eliminate all duplicate AND small elements in matrix
     253      b) remove all gaps and set extraGap_ and extraMajor_ to 0.0
     254      c) reallocate arrays and make max lengths equal to lengths
     255      d) orders elements
     256      returns number of elements eliminated or -1 if not ClpPackedMatrix
     257  */
     258  int cleanMatrix(double threshold=1.0e-20);
    251259#ifndef CLP_NO_STD
    252260  /// Drops names - makes lengthnames 0 and names empty
  • branches/devel/Clp/src/ClpNonLinearCost.cpp

    r1018 r1033  
    7979  feasibleCost_=0.0;
    8080  infeasibilityWeight_ = -1.0;
     81  double * cost = model_->costRegion();
     82  // check if all 0
     83  int iSequence;
     84  bool allZero=true;
     85  for (iSequence=0;iSequence<numberTotal1;iSequence++) {
     86    if (cost[iSequence]) {
     87      allZero=false;
     88      break;
     89    }
     90  }
     91  if (allZero)
     92    model_->setInfeasibilityCost(1.0);
    8193  double infeasibilityCost = model_->infeasibilityCost();
    8294  sumInfeasibilities_=0.0;
     
    94106  infeasible_=NULL;
    95107
    96   int iSequence;
    97108  double * upper = model_->upperRegion();
    98109  double * lower = model_->lowerRegion();
    99   double * cost = model_->costRegion();
    100110
    101111  // See how we are storing things
     
    610620            lowerValue = lower_[iRange+1];
    611621            if (value-lowerValue<-primalTolerance) {
    612               value = lowerValue-value;
     622              value = lowerValue-value-primalTolerance;
    613623#ifndef NDEBUG
    614624              if(value>1.0e15)
     
    627637            upperValue = lower_[iRange];
    628638            if (value-upperValue>primalTolerance) {
    629               value = value-upperValue;
     639              value = value-upperValue-primalTolerance;
    630640#ifndef NDEBUG
    631641              if(value>1.0e15)
     
    812822            // below
    813823            newWhere=CLP_BELOW_LOWER;
    814             double infeasibility = lowerValue-value;
     824            double infeasibility = lowerValue-value-primalTolerance;
    815825            sumInfeasibilities_ += infeasibility;
    816826            largestInfeasibility_ = CoinMax(largestInfeasibility_,infeasibility);
     
    822832          // above
    823833          newWhere = CLP_ABOVE_UPPER;
    824           double infeasibility = value-upperValue;
     834          double infeasibility = value-upperValue-primalTolerance;
    825835          sumInfeasibilities_ += infeasibility;
    826836          largestInfeasibility_ = CoinMax(largestInfeasibility_,infeasibility);
  • branches/devel/Clp/src/ClpObjective.hpp

    r754 r1033  
    4848                            double & predictedObj,
    4949                            double & thetaObj)=0;
     50  /// Return objective value (without any ClpModel offset) (model may be NULL)
     51  virtual double objectiveValue(const ClpSimplex * model, const double * solution) const = 0;
    5052  /// Resize objective
    5153  virtual void resize(int newNumberColumns) = 0;
     
    5860   */
    5961  virtual int markNonlinear(char * which);
     62  /// Say we have new primal solution - so may need to recompute
     63  virtual void newXValues() {};
    6064  //@}
    6165 
  • branches/devel/Clp/src/ClpQuadraticObjective.cpp

    r920 r1033  
    1111// Constructors / Destructor / Assignment
    1212//#############################################################################
    13 
    1413//-------------------------------------------------------------------
    1514// Default Constructor
     
    263262    delete quadraticObjective_;
    264263    quadraticObjective_ = NULL;
     264    delete [] objective_;
     265    delete [] gradient_;
    265266    ClpObjective::operator=(rhs);
    266267    numberColumns_=rhs.numberColumns_;
     
    894895  return CoinMin(theta,maximumTheta);
    895896}
     897// Return objective value (without any ClpModel offset) (model may be NULL)
     898double
     899ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
     900{
     901  bool scaling=false;
     902  if (model&&(model->rowScale()||
     903              model->objectiveScale()!=1.0))
     904    scaling=true;
     905  const double * cost = NULL;
     906  if (model)
     907    cost = model->costRegion();
     908  if (!cost) {
     909    // not in solve
     910    cost = objective_;
     911    scaling=false;
     912  }
     913  double linearCost =0.0;
     914  int numberColumns = model->numberColumns();
     915  int numberTotal = numberColumns;
     916  double currentObj=0.0;
     917  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
     918    linearCost += cost[iColumn]*solution[iColumn];
     919  }
     920  if (!activated_||!quadraticObjective_) {
     921    currentObj=linearCost;
     922    return currentObj;
     923  }
     924  assert (model);
     925  const int * columnQuadratic = quadraticObjective_->getIndices();
     926  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     927  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     928  const double * quadraticElement = quadraticObjective_->getElements();
     929  double c=0.0;
     930  if (!scaling) {
     931    if (!fullMatrix_) {
     932      int iColumn;
     933      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     934        double valueI = solution[iColumn];
     935        CoinBigIndex j;
     936        for (j=columnQuadraticStart[iColumn];
     937             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     938          int jColumn = columnQuadratic[j];
     939          double valueJ = solution[jColumn];
     940          double elementValue = quadraticElement[j];
     941          if (iColumn!=jColumn) {
     942            c += valueI*valueJ*elementValue;
     943          } else {
     944            c += 0.5*valueI*valueI*elementValue;
     945          }
     946        }
     947      }
     948    } else {
     949      // full matrix stored
     950      int iColumn;
     951      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     952        double valueI = solution[iColumn];
     953        CoinBigIndex j;
     954        for (j=columnQuadraticStart[iColumn];
     955             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     956          int jColumn = columnQuadratic[j];
     957          double valueJ = solution[jColumn];
     958          double elementValue = quadraticElement[j];
     959          valueJ *= elementValue;
     960          c += valueI*valueJ;
     961        }
     962      }
     963      c *= 0.5;
     964    }
     965  } else {
     966    // scaling
     967    // for now only if half
     968    assert (!fullMatrix_);
     969    const double * columnScale = model->columnScale();
     970    double direction = model->objectiveScale();
     971    // direction is actually scale out not scale in
     972    if (direction)
     973      direction = 1.0/direction;
     974    if (!columnScale) {
     975      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     976        double valueI = solution[iColumn];
     977        CoinBigIndex j;
     978        for (j=columnQuadraticStart[iColumn];
     979             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     980          int jColumn = columnQuadratic[j];
     981          double valueJ = solution[jColumn];
     982          double elementValue = quadraticElement[j];
     983          elementValue *= direction;
     984          if (iColumn!=jColumn) {
     985            c += valueI*valueJ*elementValue;
     986          } else {
     987            c += 0.5*valueI*valueI*elementValue;
     988          }
     989        }
     990      }
     991    } else {
     992      // scaling
     993      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     994        double valueI = solution[iColumn];
     995        double scaleI = columnScale[iColumn]*direction;
     996        CoinBigIndex j;
     997        for (j=columnQuadraticStart[iColumn];
     998             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     999          int jColumn = columnQuadratic[j];
     1000          double valueJ = solution[jColumn];
     1001          double elementValue = quadraticElement[j];
     1002          elementValue *= scaleI*columnScale[jColumn];
     1003          if (iColumn!=jColumn) {
     1004            c += valueI*valueJ*elementValue;
     1005          } else {
     1006            c += 0.5*valueI*valueI*elementValue;
     1007          }
     1008        }
     1009      }
     1010    }
     1011  }
     1012  currentObj = c+linearCost;
     1013  return currentObj;
     1014}
    8961015// Scale objective
    8971016void
  • branches/devel/Clp/src/ClpQuadraticObjective.hpp

    r754 r1033  
    4747                            double & predictedObj,
    4848                            double & thetaObj);
     49  /// Return objective value (without any ClpModel offset) (model may be NULL)
     50  virtual double objectiveValue(const ClpSimplex * model, const double * solution) const ;
    4951  virtual void resize(int newNumberColumns) ;
    5052  /// Delete columns in  objective
  • branches/devel/Clp/src/ClpSimplex.cpp

    r1009 r1033  
    35263526  if (!numberRows||!numberColumns) {
    35273527    numberRows=0;
    3528     numberColumns=0;
     3528    if (objective_->type()<2)
     3529      numberColumns=0;
    35293530  }
    35303531  if (!auxiliaryModel_) {
     
    74857486{
    74867487  model_ = model;
    7487   int i;
    7488   for (i=0;i<CLP_PROGRESS;i++) {
    7489     if (model_->algorithm()>=0)
    7490       objective_[i] = COIN_DBL_MAX;
    7491     else
    7492       objective_[i] = -COIN_DBL_MAX;
    7493     infeasibility_[i] = -1.0; // set to an impossible value
    7494     realInfeasibility_[i] = COIN_DBL_MAX;
    7495     numberInfeasibilities_[i]=-1;
    7496     iterationNumber_[i]=-1;
    7497   }
    7498 #ifdef CLP_PROGRESS_WEIGHT
    7499   for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
    7500     objectiveWeight_[i] = COIN_DBL_MAX;
    7501     infeasibilityWeight_[i] = -1.0; // set to an impossible value
    7502     realInfeasibilityWeight_[i] = COIN_DBL_MAX;
    7503     numberInfeasibilitiesWeight_[i]=-1;
    7504     iterationNumberWeight_[i]=-1;
    7505   }
    7506   drop_ =0.0;
    7507   best_ =0.0;
    7508 #endif
     7488  reset();
    75097489  initialWeight_=0.0;
    7510   for (i=0;i<CLP_CYCLE;i++) {
    7511     //obj_[i]=COIN_DBL_MAX;
    7512     in_[i]=-1;
    7513     out_[i]=-1;
    7514     way_[i]=0;
    7515   }
    7516   numberTimes_ = 0;
    7517   numberBadTimes_ = 0;
    7518   oddState_=0;
    75197490}
    75207491// Assignment operator. This copies the data
     
    77177688  return -1;
    77187689}
     7690// Resets as much as possible
     7691void
     7692ClpSimplexProgress::reset()
     7693{
     7694  int i;
     7695  for (i=0;i<CLP_PROGRESS;i++) {
     7696    if (model_->algorithm()>=0)
     7697      objective_[i] = COIN_DBL_MAX;
     7698    else
     7699      objective_[i] = -COIN_DBL_MAX;
     7700    infeasibility_[i] = -1.0; // set to an impossible value
     7701    realInfeasibility_[i] = COIN_DBL_MAX;
     7702    numberInfeasibilities_[i]=-1;
     7703    iterationNumber_[i]=-1;
     7704  }
     7705#ifdef CLP_PROGRESS_WEIGHT
     7706  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7707    objectiveWeight_[i] = COIN_DBL_MAX;
     7708    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     7709    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     7710    numberInfeasibilitiesWeight_[i]=-1;
     7711    iterationNumberWeight_[i]=-1;
     7712  }
     7713  drop_ =0.0;
     7714  best_ =0.0;
     7715#endif
     7716  for (i=0;i<CLP_CYCLE;i++) {
     7717    //obj_[i]=COIN_DBL_MAX;
     7718    in_[i]=-1;
     7719    out_[i]=-1;
     7720    way_[i]=0;
     7721  }
     7722  numberTimes_ = 0;
     7723  numberBadTimes_ = 0;
     7724  oddState_=0;
     7725}
    77197726// Returns previous objective (if -1) - current if (0)
    77207727double
  • branches/devel/Clp/src/ClpSimplex.hpp

    r998 r1033  
    202202  /// Puts solution back into small model
    203203  void getbackSolution(const ClpSimplex & smallModel,const int * whichRow, const int * whichColumn);
     204  /** Load nonlinear part of problem from AMPL info
     205      Returns 0 if linear
     206      1 if quadratic objective
     207      2 if quadratic constraints
     208      3 if nonlinear objective
     209      4 if nonlinear constraints
     210      -1 on failure
     211  */
     212  int loadNonLinear(void * info, int & numberConstraints,
     213                    ClpConstraint ** & constraints);
    204214  //@}
    205215
     
    12961306  /// Destructor
    12971307   ~ClpSimplexProgress (  );
     1308  /// Resets as much as possible
     1309   void reset();
    12981310  //@}
    12991311
  • branches/devel/Clp/src/ClpSimplexDual.cpp

    r1010 r1033  
    31883188    // so (dualIn_+modification)==theta_*alpha_
    31893189    double modification = theta_*alpha_-dualIn_;
     3190    // But should not move objectivetoo much ??
     3191#define DONT_MOVE_OBJECTIVE
     3192#ifdef DONT_MOVE_OBJECTIVE
     3193    double moveObjective = fabs(modification*solution_[sequenceIn_]);
     3194    double smallMove = CoinMax(fabs(objectiveValue_),1.0e-3);
     3195    if (moveObjective>smallMove) {
     3196      printf("would move objective by %g - original mod %g sol value %g\n",moveObjective,
     3197             modification,solution_[sequenceIn_]);
     3198      modification *= smallMove/moveObjective;
     3199    }
     3200#endif
    31903201    if ((specialOptions_&(2048+4096))!=0) {
    31913202      if ((specialOptions_&16384)!=0) {
  • branches/devel/Clp/src/ClpSimplexNonlinear.cpp

    r1007 r1033  
    27142714  double offset=0.0;
    27152715  double objValue=0.0;
     2716  int exitPass = 2*numberPasses+10;
    27162717  for (iPass=0;iPass<numberPasses;iPass++) {
     2718    exitPass--;
    27172719    // redo objective
    27182720    offset=0.0;
    27192721    objValue=-objectiveOffset;
     2722    // make sure x updated
     2723    trueObjective->newXValues();
    27202724    double theta=-1.0;
    27212725    double maxTheta=COIN_DBL_MAX;
     
    28002804    double theta2 = trueObjective->stepLength(this,saveSolution,changeRegion,maxTheta2,
    28012805                                             objValue,predictedObj,thetaObj);
    2802    
     2806    int lastMoveStatus = goodMove;
    28032807    if (goodMove>=0) {
    28042808      theta = CoinMin(theta2,maxTheta);
     
    28232827#define LOCAL
    28242828#ifdef LOCAL
     2829        bool absolutelyOptimal=false;
    28252830        int saveScaling = scalingFlag();
    28262831        scaling(0);
     
    28882893            directionVector(longArray,spare,rowArray,0,
    28892894                            normFlagged,normUnflagged,numberNonBasic);
    2890             if (normFlagged+normUnflagged<1.0e-8)
     2895            if (normFlagged+normUnflagged<1.0e-8) {
     2896              absolutelyOptimal=true;
    28912897              break;  //looks optimal
     2898            }
    28922899            double djNorm=0.0;
    28932900            int iIndex;
     
    29022909            if (!numberNonBasic) {
    29032910              // looks optimal
     2911              absolutelyOptimal=true;
    29042912              break;
    29052913            }
     
    30093017        memset(rowActivity_,0,numberRows_*sizeof(double));
    30103018        times(1.0,solution,rowActivity_);
    3011         if (theta>0.99999&&theta2<1.9) {
     3019        if (theta>0.99999&&theta2<1.9&&!absolutelyOptimal) {
    30123020          // If big changes then tighten
    30133021          /*  thetaObj is objvalue + a*2*2 +b*2
     
    30283036    memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
    30293037           numberColumns*sizeof(double));
     3038    //printf("offset comp %g orig %g\n",offset,objectiveOffset);
    30303039    // could do faster
    30313040    trueObjective->stepLength(this,solution,changeRegion,0.0,
    30323041                                             objValue,predictedObj,thetaObj);
     3042    printf("offset comp %g orig %g - obj from stepLength %g\n",offset,objectiveOffset,objValue);
    30333043    setDblParam(ClpObjOffset,objectiveOffset+offset);
    30343044    int * temp=last[2];
     
    30563066    double maxDelta=0.0;
    30573067    if (goodMove>=0) {
    3058       if (objValue<=lastObjective+1.0e-15*fabs(lastObjective))
     3068      if (objValue-lastObjective<=1.0e-15*fabs(lastObjective))
    30593069        goodMove=1;
    30603070      else
     
    30993109        trust[jNon] *=0.0001;
    31003110    }
     3111    if(lastMoveStatus==-1&&goodMove==-1)
     3112      goodMove=1; // to force solve
    31013113    if (goodMove>0) {
    31023114      double drop = lastObjective-objValue;
     
    32613273      memcpy(status_,saveStatus,numberRows+numberColumns);
    32623274      iPass--;
     3275      assert (exitPass>0);
    32633276      goodMove=-1;
    32643277    }
     
    33213334{
    33223335  if (!numberConstraints) {
    3323     // no nonlinear part
    3324     return ClpSimplex::primal(0);
     3336    // no nonlinear constraints - may be nonlinear objective
     3337    return primalSLP(numberPasses,deltaTolerance);
    33253338  }
    33263339  // check all matrix for odd rows is empty
     
    33373350  //const double * elementByRow = copy.getElements();
    33383351  int numberArtificials=0;
     3352  // see how many extra we need
     3353  CoinBigIndex numberExtra=0;
    33393354  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
    33403355    ClpConstraint * constraint = constraints[iConstraint];
    33413356    int iRow = constraint->rowNumber();
    3342     if (iRow>=0) {
    3343       if (iRow>=numberRows_)
    3344         numberBad++;
    3345       else if (rowLength[iRow])
    3346         numberBad++;
    3347       if (rowLower_[iRow]>-1.0e20)
    3348         numberArtificials++;
    3349       if (rowUpper_[iRow]<1.0e20)
    3350         numberArtificials++;
    3351     } else {
    3352       // objective
    3353       assert (iRow==-1);
    3354       ClpLinearObjective * linearObj = (dynamic_cast< ClpLinearObjective*>(objective_));
    3355       assert (linearObj);
    3356       // check empty
    3357       const double * obj = objective();
    3358       int iColumn;
    3359       for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3360         if (obj[iColumn])
    3361           numberBad++;
    3362       }
    3363     }
     3357    assert (iRow>=0);
     3358    int n = constraint->numberCoefficients() -rowLength[iRow];
     3359    numberExtra += n;
     3360    if (iRow>=numberRows_)
     3361      numberBad++;
     3362    else if (rowLength[iRow]&&n)
     3363      numberBad++;
     3364    if (rowLower_[iRow]>-1.0e20)
     3365      numberArtificials++;
     3366    if (rowUpper_[iRow]<1.0e20)
     3367      numberArtificials++;
    33643368  }
    33653369  if (numberBad)
    33663370    return numberBad;
     3371  ClpObjective * trueObjective = NULL;
     3372  if (objective_->type()>=2) {
     3373    // Replace objective
     3374    trueObjective = objective_;
     3375    objective_=new ClpLinearObjective(NULL,numberColumns_);
     3376  }
    33673377  ClpSimplex newModel(*this);
     3378  // we can put back true objective
     3379  if (trueObjective) {
     3380    // Replace objective
     3381    delete objective_;
     3382    objective_=trueObjective;
     3383  }
    33683384  int numberColumns2 = numberColumns_;
    33693385  double penalty=1.0e9;
     
    33793395      ClpConstraint * constraint = constraints[iConstraint];
    33803396      int iRow = constraint->rowNumber();
    3381       if (iRow>=0) {
    3382         if (rowLower_[iRow]>-1.0e20) {
    3383           addRow[numberArtificials]=iRow;
    3384           addElement[numberArtificials]=1.0;
    3385           addCost[numberArtificials]=penalty;
    3386           numberArtificials++;
    3387           addStarts[numberArtificials]=numberArtificials;
    3388         }
    3389         if (rowUpper_[iRow]<1.0e20) {
    3390           addRow[numberArtificials]=iRow;
    3391           addElement[numberArtificials]=-1.0;
    3392           addCost[numberArtificials]=penalty;
    3393           numberArtificials++;
    3394           addStarts[numberArtificials]=numberArtificials;
    3395         }
     3397      if (rowLower_[iRow]>-1.0e20) {
     3398        addRow[numberArtificials]=iRow;
     3399        addElement[numberArtificials]=1.0;
     3400        addCost[numberArtificials]=penalty;
     3401        numberArtificials++;
     3402        addStarts[numberArtificials]=numberArtificials;
     3403      }
     3404      if (rowUpper_[iRow]<1.0e20) {
     3405        addRow[numberArtificials]=iRow;
     3406        addElement[numberArtificials]=-1.0;
     3407        addCost[numberArtificials]=penalty;
     3408        numberArtificials++;
     3409        addStarts[numberArtificials]=numberArtificials;
    33963410      }
    33973411    }
     
    34113425    constraint->markNonlinear(mark);
    34123426  }
     3427  if (trueObjective)
     3428    trueObjective->markNonlinear(mark);
    34133429  int iColumn;
    34143430  int numberNonLinearColumns=0;
     
    34263442  rowLength = copy.getVectorLengths();
    34273443  const double * elementByRow = copy.getElements();
    3428   CoinBigIndex numberExtra=0;
    3429   for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
    3430     ClpConstraint * constraint = constraints[iConstraint];
    3431     numberExtra+=constraint->markNonzero(mark);
    3432   }
    34333444  numberExtra +=copy.getNumElements();
    34343445  CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1];
     
    34483459  for (iRow=0;iRow<numberRows_;iRow++) {
    34493460    if (backRow[iRow]<0) {
    3450       // copy normal
     3461      // copy normal 
    34513462      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
    34523463           j++) {
     
    34553466      }
    34563467    } else {
    3457       // all possible
    3458       memset(mark,0,numberColumns_);
    34593468      ClpConstraint * constraint = constraints[backRow[iRow]];
    34603469      assert(iRow == constraint->rowNumber());
    3461       constraint->markNonzero(mark);
    3462       for (int k=0;k<numberColumns_;k++) {
    3463         if (mark[k]) {
    3464           newColumn[numberExtra]=k;
    3465           newElement[numberExtra++] = 1.0;
    3466         }
    3467       }
    3468       // copy artificials
    3469       for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
    3470            j++) {
    3471         newColumn[numberExtra]=column[j];
    3472         newElement[numberExtra++] = elementByRow[j];
     3470      int numberArtificials=0;
     3471      if (rowLower_[iRow]>-1.0e20)
     3472        numberArtificials++;
     3473      if (rowUpper_[iRow]<1.0e20)
     3474        numberArtificials++;
     3475      if (numberArtificials==rowLength[iRow]) {
     3476        // all possible
     3477        memset(mark,0,numberColumns_);
     3478        constraint->markNonzero(mark);
     3479        for (int k=0;k<numberColumns_;k++) {
     3480          if (mark[k]) {
     3481            newColumn[numberExtra]=k;
     3482            newElement[numberExtra++] = 1.0;
     3483          }
     3484        }
     3485        // copy artificials
     3486        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3487             j++) {
     3488          newColumn[numberExtra]=column[j];
     3489          newElement[numberExtra++] = elementByRow[j];
     3490        }
     3491      } else {
     3492        // there already
     3493        // copy
     3494        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3495             j++) {
     3496          newColumn[numberExtra]=column[j];
     3497          assert (elementByRow[j]);
     3498          newElement[numberExtra++] = elementByRow[j];
     3499        }
    34733500      }
    34743501    }
     
    34823509  delete [] newElement;
    34833510  delete [] mark;
     3511  // get feasible
     3512  newModel.primal(1);
     3513  // still infeasible
     3514  if (newModel.numberPrimalInfeasibilities()) {
     3515    delete [] listNonLinearColumn;
     3516    return 0;
     3517  }
    34843518  int numberRows = newModel.numberRows();
    34853519  double * columnLower = newModel.columnLower();
     
    35023536    double upper = columnUpper[iColumn];
    35033537    double lower = columnLower[iColumn];
     3538    if (solution[iColumn]<trueLower[jNon])
     3539      solution[iColumn]=trueLower[jNon];
     3540    else if (solution[iColumn]>trueUpper[jNon])
     3541      solution[iColumn]=trueUpper[jNon];
     3542#if 0
     3543    double large = CoinMax(1000.0,10.0*fabs(solution[iColumn]));
    35043544    if (upper>1.0e10)
    3505       upper = max(1000.0,10*solution[iColumn]);
     3545      upper = solution[iColumn]+large;
     3546    if (lower<-1.0e10)
     3547      lower = solution[iColumn]-large;
     3548#else
     3549    upper = solution[iColumn]+0.5;
     3550    lower = solution[iColumn]-0.5;
     3551#endif
    35063552    //columnUpper[iColumn]=upper;
    35073553    trust[jNon]=0.05*(1.0+upper-lower);
    3508     trueLower[jNon]=lower;
     3554    trueLower[jNon]=columnLower[iColumn];
    35093555    //trueUpper[jNon]=upper;
    35103556    trueUpper[jNon]=columnUpper[iColumn];
     
    35343580      else if (solution[iColumn]>trueUpper[jNon])
    35353581        solution[iColumn]=trueUpper[jNon];
    3536       if (iPass) {
    3537         columnLower[iColumn]=max(solution[iColumn]
    3538                                  -trust[jNon],
    3539                                  trueLower[jNon]);
    3540         if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
    3541           columnLower[iColumn]=SMALL_FIX;
    3542         columnUpper[iColumn]=min(solution[iColumn]
    3543                                  +trust[jNon],
    3544                                  trueUpper[jNon]);
    3545         if (!trueLower[jNon])
    3546           columnUpper[iColumn] = max(columnUpper[iColumn],
    3547                                      columnLower[iColumn]+SMALL_FIX);
    3548         if (!trueLower[jNon]&&tryFix&&
    3549             columnLower[iColumn]==SMALL_FIX&&
    3550             columnUpper[iColumn]<3.0*SMALL_FIX) {
    3551           columnLower[iColumn]=0.0;
    3552           solution[iColumn]=0.0;
    3553           columnUpper[iColumn]=0.0;
    3554           printf("fixing %d\n",iColumn);
    3555         }
    3556       } else {
    3557 #if 0
    3558         // fix nonlinear variables
    3559         columnLower[iColumn]=solution[iColumn];
    3560         columnUpper[iColumn]=solution[iColumn];
    3561 #endif
     3582      columnLower[iColumn]=CoinMax(solution[iColumn]
     3583                                   -trust[jNon],
     3584                                   trueLower[jNon]);
     3585      if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
     3586        columnLower[iColumn]=SMALL_FIX;
     3587      columnUpper[iColumn]=CoinMin(solution[iColumn]
     3588                                   +trust[jNon],
     3589                                   trueUpper[jNon]);
     3590      if (!trueLower[jNon])
     3591        columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
     3592                                       columnLower[iColumn]+SMALL_FIX);
     3593      if (!trueLower[jNon]&&tryFix&&
     3594          columnLower[iColumn]==SMALL_FIX&&
     3595          columnUpper[iColumn]<3.0*SMALL_FIX) {
     3596        columnLower[iColumn]=0.0;
     3597        solution[iColumn]=0.0;
     3598        columnUpper[iColumn]=0.0;
     3599        printf("fixing %d\n",iColumn);
    35623600      }
    35633601    }
     
    35693607    rowStart = newMatrix.getVectorStarts();
    35703608    rowLength = newMatrix.getVectorLengths();
     3609    // make sure x updated
     3610    if (numberConstraints)
     3611      constraints[0]->newXValues();
     3612    else
     3613      trueObjective->newXValues();
    35713614    double * changeableElement = newMatrix.getMutableElements();
     3615    if (trueObjective) {
     3616      memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
     3617             numberColumns_*sizeof(double));
     3618    } else {
     3619      memcpy(objective,objective_->gradient(this,solution,offset,true,2),
     3620             numberColumns_*sizeof(double));
     3621    }
    35723622    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
    35733623      ClpConstraint * constraint = constraints[iConstraint];
    35743624      int iRow = constraint->rowNumber();
    35753625      double functionValue;
    3576       int numberErrors = constraint->gradient(&newModel,solution,gradient,functionValue,offset);
     3626#ifndef NDEBUG
     3627      int numberErrors =
     3628#endif
     3629      constraint->gradient(&newModel,solution,gradient,functionValue,offset);
    35773630      assert (!numberErrors);
    3578       if (iRow>=0) {
    3579         if (numberRows<50)
    3580           printf("For row %d current value is %g (activity %g) , dual is %g - offset %g\n",iRow,functionValue,
    3581                  newModel.primalRowSolution()[iRow],
    3582                  newModel.dualRowSolution()[iRow],offset);
    3583         int numberCoefficients = constraint->numberCoefficients();
    3584         for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+numberCoefficients;j++) {
    3585           int iColumn = column[j];
    3586           changeableElement[j] = gradient[iColumn];
    3587           gradient[iColumn]=0.0;
    3588         }
    3589         for (int k=0;k<numberColumns_;k++)
    3590           assert (!gradient[k]);
    3591         if (rowLower_[iRow]>-1.0e20)
    3592           rowLower[iRow] = rowLower_[iRow] - offset;
    3593         if (rowUpper_[iRow]<1.0e20)
    3594           rowUpper[iRow] = rowUpper_[iRow] - offset;
    3595       } else {
    3596         memcpy(objective,gradient,numberColumns_*sizeof(double));
    3597         printf("objective offset %g\n",offset);
    3598         objectiveOffset2 =objectiveOffset-offset;
    3599         newModel.setObjectiveOffset(objectiveOffset2);
    3600       }
     3631      double dualValue = newModel.dualRowSolution()[iRow];
     3632      int numberCoefficients = constraint->numberCoefficients();
     3633      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+numberCoefficients;j++) {
     3634        int iColumn = column[j];
     3635        changeableElement[j] = gradient[iColumn];
     3636        //objective[iColumn] -= dualValue*gradient[iColumn];
     3637        gradient[iColumn]=0.0;
     3638      }
     3639      for (int k=0;k<numberColumns_;k++)
     3640        assert (!gradient[k]);
     3641      if (rowLower_[iRow]>-1.0e20)
     3642        rowLower[iRow] = rowLower_[iRow] - offset;
     3643      if (rowUpper_[iRow]<1.0e20)
     3644        rowUpper[iRow] = rowUpper_[iRow] - offset;
    36013645    }
    36023646    // Replace matrix
     
    36053649    columnCopy->reverseOrderedCopyOf(newMatrix);
    36063650    newModel.replaceMatrix(columnCopy,true);
     3651    // solve
     3652    newModel.primal(1);
     3653    if (newModel.status()==1) {
     3654      // Infeasible!
     3655      newModel.allSlackBasis();
     3656      newModel.primal();
     3657      newModel.writeMps("infeas.mps");
     3658      assert(!newModel.status());
     3659    }
     3660    double sumInfeas=0.0;
     3661    int numberInfeas=0;
     3662    for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
     3663      if (solution[iColumn]>1.0e-8) {
     3664        numberInfeas++;
     3665        sumInfeas += solution[iColumn];
     3666      }
     3667    }
     3668    printf("%d artificial infeasibilities - summing to %g\n",
     3669           numberInfeas,sumInfeas);
     3670    double infValue=0.0;
    36073671    double objValue=0.0;
    3608     double infValue=0.0;
    3609     {
    3610       // get matrix data pointers
    3611       const int * row = columnCopy->getIndices();
    3612       const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
    3613       const int * columnLength = columnCopy->getVectorLengths();
    3614       const double * elementByColumn = columnCopy->getElements();
    3615       int iRow,iColumn;
    3616       double * accumulate = new double [numberRows];
    3617       //const double * solution = newModel.primalColumnSolution();
    3618       double * partSolution = new double [numberColumns2];
    3619       memcpy(partSolution,solution,numberColumns_*sizeof(double));
    3620       memset(partSolution+numberColumns_,0,
    3621              (numberColumns2-numberColumns_)*sizeof(double));
    3622       memset(accumulate,0,numberRows*sizeof(double));
    3623       //newModel.matrix()->times(partSolution,accumulate);
    3624       for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3625         if (solution[iColumn]) {
    3626           double value=solution[iColumn];
    3627           objValue += objective[iColumn]*solution[iColumn];
    3628           int j;
    3629           for (j=columnStart[iColumn];
    3630                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    3631             int iRow=row[j];
    3632             accumulate[iRow] += value*elementByColumn[j];
    3633           }
    3634         }
    3635       }
    3636       for (iRow=0;iRow<numberRows;iRow++) {
    3637         if (accumulate[iRow]<rowLower[iRow]-1.0e-5) {
    3638           double under = rowLower[iRow]-accumulate[iRow];
    3639           infValue += under;
    3640 #if 0
    3641           printf("Row %d l=%g u=%g a=%g",iRow,rowLower[iRow],
    3642                  rowUpper[iRow],accumulate[iRow]);
    3643           printf(" ** under by %g",under);
    3644           printf("\n");
    3645 #endif
    3646         } else if (accumulate[iRow]>rowUpper[iRow]+1.0e-5) {
    3647           double over = accumulate[iRow]-rowUpper[iRow];
    3648           infValue += over;
    3649 #if 0
    3650           printf("Row %d l=%g u=%g a=%g",iRow,rowLower[iRow],
    3651                  rowUpper[iRow],accumulate[iRow]);
    3652           printf(" ** over by %g",over);
    3653           printf("\n");
    3654 #endif
    3655         }
    3656       }
    3657       if (infValue)
    3658         printf("Sum infeasibilities %g ",infValue);
    3659       if (infValue<0.1)
    3660         infValue=0.0;
    3661       infValue *= penalty;
    3662       if (infValue)
    3663         printf("Infeasible obj %g ",infValue);
    3664       delete [] accumulate; 
    3665       delete [] partSolution;
    3666     }
     3672    // make sure x updated
     3673    if (numberConstraints)
     3674      constraints[0]->newXValues();
     3675    else
     3676      trueObjective->newXValues();
     3677    if (trueObjective) {
     3678      objValue=trueObjective->objectiveValue(this,solution);
     3679      printf("objective offset %g\n",offset);
     3680      objectiveOffset2 =objectiveOffset+offset; // ? sign
     3681      newModel.setObjectiveOffset(objectiveOffset2);
     3682    } else {
     3683      objValue=objective_->objectiveValue(this,solution);
     3684    }
     3685    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3686      ClpConstraint * constraint = constraints[iConstraint];
     3687      int iRow = constraint->rowNumber();
     3688      double functionValue=constraint->functionValue(this,solution);
     3689      double dualValue = newModel.dualRowSolution()[iRow];
     3690      if (numberConstraints<50)
     3691        printf("For row %d current value is %g (activity %g) , dual is %g - offset %g\n",iRow,functionValue,
     3692               newModel.primalRowSolution()[iRow],
     3693               dualValue,offset);
     3694      if (functionValue<rowLower_[iRow]-1.0e-5) {
     3695        double under = rowLower_[iRow]-functionValue;
     3696        infValue += under;
     3697      } else if (functionValue>rowUpper_[iRow]+1.0e-5) {
     3698        double over = functionValue-rowUpper_[iRow];
     3699        infValue += over;
     3700      }
     3701    }
     3702    if (infValue)
     3703      printf("Sum infeasibilities %g ",infValue);
     3704    if (infValue<0.1)
     3705      infValue=0.0;
     3706    infValue *= penalty;
     3707    if (infValue)
     3708      printf("Infeasible obj %g ",infValue);
    36673709    if (objectiveOffset2)
    36683710      printf("offset2 %g\n",objectiveOffset2);
     
    37243766        assert (gap>=0.0);
    37253767        if (gap)
    3726           smallestGap = min(smallestGap,gap);
     3768          smallestGap = CoinMin(smallestGap,gap);
    37273769      }
    37283770      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     
    37313773        assert (gap>=0.0);
    37323774        if (gap) {
    3733           smallestNonLinearGap = min(smallestNonLinearGap,gap);
     3775          smallestNonLinearGap = CoinMin(smallestNonLinearGap,gap);
    37343776          if (gap<1.0e-7&&iPass==1) {
    37353777            printf("Small gap %d %d %g %g %g\n",
     
    37403782          }
    37413783        }
    3742         maxDelta = max(maxDelta,
     3784        maxDelta = CoinMax(maxDelta,
    37433785                       fabs(solution[iColumn]-saveSolution[iColumn]));
    37443786        if (last[0][jNon]*last[1][jNon]<0) {
     
    37573799          }
    37583800        }
    3759         smallestTrust = min(smallestTrust,trust[jNon]);
     3801        smallestTrust = CoinMin(smallestTrust,trust[jNon]);
    37603802      }
    37613803      std::cout<<"largest delta is "<<maxDelta
     
    37733815            tryFix=true;
    37743816            for (jNon=0;jNon<numberNonLinearColumns;jNon++)
    3775               trust[jNon] = max(trust[jNon],1.0e-1);
     3817              trust[jNon] = CoinMax(trust[jNon],1.0e-1);
    37763818            numberZeroPasses=0;
    37773819          }
     
    38013843        if (dj<-1.0e-6) {
    38023844          double drop = -dj*(columnUpper[iColumn]-solution[iColumn]);
    3803           //double upper = min(trueUpper[jNon],solution[iColumn]+0.1);
     3845          //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
    38043846          //double drop2 = -dj*(upper-solution[iColumn]);
    38053847#if 0
     
    38183860        } else if (dj>1.0e-6) {
    38193861          double drop = -dj*(columnLower[iColumn]-solution[iColumn]);
    3820           //double lower = max(trueLower[jNon],solution[iColumn]-0.1);
     3862          //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
    38213863          //double drop2 = -dj*(lower-solution[iColumn]);
    38223864#if 0
     
    38603902    writer.writeMps("xx.mps");
    38613903#endif
    3862     // solve
    3863     newModel.primal(1);
    3864     if (newModel.status()==1) {
    3865       // Infeasible!
    3866       newModel.allSlackBasis();
    3867       newModel.primal();
    3868       newModel.writeMps("infeas.mps");
    3869       assert(!newModel.status());
    3870     }
    3871     double sumInfeas=0.0;
    3872     int numberInfeas=0;
    3873     for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
    3874       if (solution[iColumn]>1.0e-8) {
    3875         numberInfeas++;
    3876         sumInfeas += solution[iColumn];
    3877       }
    3878     }
    3879     printf("%d infeasibilities - summing to %g\n",
    3880            numberInfeas,sumInfeas);
    38813904  }
    38823905  delete [] saveSolution;
  • branches/devel/Clp/src/ClpSimplexPrimal.cpp

    r1010 r1033  
    10151015    infeasibilityCost_ = CoinMin(CoinMax(100.0*sumDualInfeasibilities_,1.0e8),1.00000001e10);
    10161016    // reset looping criterion
    1017     *progress = ClpSimplexProgress();
     1017    progress->reset();
    10181018    trueInfeasibility = 1.123456e10;
    10191019  }
     
    10211021    // If infeasibility going up may change weights
    10221022    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
    1023     double lastInf = progress->lastInfeasibility();
    1024     if(lastInf<testValue||trueInfeasibility==1.123456e10) {
     1023    double lastInf = progress->lastInfeasibility(1);
     1024    double lastInf3 = progress->lastInfeasibility(3);
     1025    double thisObj = progress->lastObjective(0);
     1026    double thisInf = progress->lastInfeasibility(0);
     1027    thisObj += infeasibilityCost_*thisInf;
     1028    double lastObj = progress->lastObjective(1);
     1029    lastObj += infeasibilityCost_*lastInf;
     1030    double lastObj3 = progress->lastObjective(3);
     1031    lastObj3 += infeasibilityCost_*lastInf3;
     1032    if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-7
     1033        &&firstFree_<0) {
     1034      if (handler_->logLevel()==63)
     1035        printf("lastobj %g this %g force %d ",lastObj,thisObj,forceFactorization_);
     1036      int maxFactor = factorization_->maximumPivots();
     1037      if (maxFactor>10) {
     1038        if (forceFactorization_<0)
     1039          forceFactorization_= maxFactor;
     1040        forceFactorization_ = CoinMax(1,(forceFactorization_>>2));
     1041        if (handler_->logLevel()==63)
     1042          printf("Reducing factorization frequency to %d\n",forceFactorization_);
     1043      }
     1044    } else if (lastObj3<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj3))-1.0e-7
     1045        &&firstFree_<0) {
     1046      if (handler_->logLevel()==63)
     1047        printf("lastobj3 %g this3 %g `force %d ",lastObj3,thisObj,forceFactorization_);
     1048      int maxFactor = factorization_->maximumPivots();
     1049      if (maxFactor>10) {
     1050        if (forceFactorization_<0)
     1051          forceFactorization_= maxFactor;
     1052        forceFactorization_ = CoinMax(1,(forceFactorization_*2)/3);
     1053        if (handler_->logLevel()==63)
     1054        printf("Reducing factorization frequency to %d\n",forceFactorization_);
     1055      }
     1056    } else if(lastInf<testValue||trueInfeasibility==1.123456e10) {
    10251057      if (infeasibilityCost_<1.0e14) {
    10261058        infeasibilityCost_ *= 1.5;
    10271059        // reset looping criterion
    1028         *progress = ClpSimplexProgress();
    1029         //printf("increasing weight to %g\n",infeasibilityCost_);
     1060        progress->reset();
     1061        if (handler_->logLevel()==63)
     1062          printf("increasing weight to %g\n",infeasibilityCost_);
    10301063        gutsOfSolution(NULL,NULL,ifValuesPass!=0);
    10311064      }
     
    10671100        numberPrimalInfeasibilities_ = ninfeas;
    10681101        sumPrimalInfeasibilities_ = sum;
    1069         bool unflagged = unflag();
     1102#ifdef COIN_DEVELOP
     1103        bool unflagged =
     1104#endif
     1105        unflag();
    10701106#ifdef COIN_DEVELOP
    10711107        if (unflagged&&handler_->logLevel()>0)
     
    11431179          infeasibilityCost_ *= 5.0;
    11441180          // reset looping criterion
    1145           *progress = ClpSimplexProgress();
     1181          progress->reset();
    11461182          changeMade_++; // say change made
    11471183          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
     
    12571293          infeasibilityCost_ = 1.0e13;
    12581294          // reset looping criterion
    1259           *progress = ClpSimplexProgress();
     1295          progress->reset();
    12601296          unPerturb(); // stop any further perturbation
    12611297        }
     
    12641300          infeasibilityCost_ *= 5.0;
    12651301          // reset looping criterion
    1266           *progress = ClpSimplexProgress();
     1302          progress->reset();
    12671303          changeMade_++; // say change made
    12681304          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
     
    13551391    nextSuperBasic(1,NULL);
    13561392  }
    1357 #if 0
    1358   double thisObj = progress->lastObjective(0);
    1359   double lastObj = progress->lastObjective(1);
    1360   if (lastObj<thisObj-1.0e-7*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-8
    1361       &&firstFree_<0) {
    1362     int maxFactor = factorization_->maximumPivots();
    1363     if (maxFactor>10) {
    1364       if (forceFactorization_<0)
    1365         forceFactorization_= maxFactor;
    1366       forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
    1367       printf("Reducing factorization frequency\n");
    1368     }
    1369   }
    1370 #endif
    13711393  // Allow matrices to be sorted etc
    13721394  int fake=-999; // signal sort
  • branches/devel/Clp/src/Idiot.cpp

    r1011 r1033  
    519519      if (iCol>=0) {
    520520        CoinBigIndex j=columnStart[iCol];
     521#ifndef NDEBUG
    521522        int iRow = row[j];
     523#endif
    522524        assert (element[j]>0.0);
    523525        assert (iRow==i);
     
    528530          iCol = nextSlack[iCol];
    529531          CoinBigIndex j=columnStart[iCol];
     532#ifndef NDEBUG
    530533          int iRow = row[j];
     534#endif
    531535          assert (element[j]>0.0);
    532536          assert (iRow==i);
     
    551555      if (iCol>=0) {
    552556        CoinBigIndex j=columnStart[iCol];
     557#ifndef NDEBUG
    553558        int iRow = row[j];
     559#endif
    554560        assert (element[j]<0.0);
    555561        assert (iRow==i);
     
    560566          iCol = nextSlack[iCol];
    561567          CoinBigIndex j=columnStart[iCol];
     568#ifndef NDEBUG
    562569          int iRow = row[j];
     570#endif
    563571          assert (element[j]<0.0);
    564572          assert (iRow==i);
     
    869877            (result.infeas>lastResult.infeas*0.9
    870878             &&result.weighted>lastResult.weighted
    871              -dropEnoughWeighted_*fabs(lastResult.weighted))) {
     879             -dropEnoughWeighted_*CoinMax(fabs(lastResult.weighted),fabs(result.weighted)))) {
    872880          mu*=changeMu;
    873881          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas&&0) {
     
    875883            printf("reasonable infeas now %g\n",reasonableInfeas);
    876884          }
     885          result.weighted=1.0e60;
    877886          nTry=0;
    878887          bestFeasible=1.0e31;
Note: See TracChangeset for help on using the changeset viewer.