Changeset 998 for branches


Ignore:
Timestamp:
May 3, 2007 11:45:47 AM (12 years ago)
Author:
forrest
Message:

for nonlinear

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

Legend:

Unmodified
Added
Removed
  • branches/devel/Clp/src/ClpConstraint.hpp

    r997 r998  
    2222 
    2323  /** Fills gradient.  If Linear then solution may be NULL,
    24       also returns value of function
     24      also returns true value of function and offset so we can use x not deltaX in constraint
    2525      If refresh is false then uses last solution
    2626      Uses model for scaling
     
    3131                       double * gradient,
    3232                       double & functionValue ,
     33                       double & offset,
    3334                       bool refresh=true) const =0;
    3435  /// Resize constraint
     
    4243   */
    4344  virtual int markNonlinear(char * which) const = 0;
     45  /** Given a zeroed array sets possible nonzero coefficients to 1.
     46      Returns number of nonzeros
     47   */
     48  virtual int markNonzero(char * which) const = 0;
    4449  //@}
    4550 
  • branches/devel/Clp/src/ClpConstraintLinear.cpp

    r997 r998  
    2828// Useful Constructor
    2929//-------------------------------------------------------------------
    30 ClpConstraintLinear::ClpConstraintLinear (int numberCoefficents ,
     30ClpConstraintLinear::ClpConstraintLinear (int row, int numberCoefficents ,
    3131                                          int numberColumns,
    3232                                          const int * column, const double * coefficient)
     
    3434{
    3535  type_=0;
     36  rowNumber_=row;
    3637  numberColumns_ = numberColumns;
    3738  numberCoefficients_=numberCoefficents;
     
    9293                              const double * solution,
    9394                              double * gradient,
    94                               double & functionValue, bool refresh) const
     95                              double & functionValue,
     96                              double & offset,
     97                              bool refresh) const
    9598{
    9699  if (refresh||!lastGradient_) {
     
    121124  }
    122125  functionValue = functionValue_;
     126  offset=0.0;
    123127  memcpy(gradient,lastGradient_,numberColumns_*sizeof(double));
    124128  return 0;
     
    178182  return 0;
    179183}
     184/* Given a zeroed array sets possible nonzero coefficients to 1.
     185   Returns number of nonzeros
     186*/
     187int
     188ClpConstraintLinear::markNonzero(char * which) const
     189{
     190  for (int i=0;i<numberCoefficients_;i++) {
     191    int iColumn = column_[i];
     192    which[iColumn]=1;
     193  }
     194  return numberCoefficients_;
     195}
  • branches/devel/Clp/src/ClpConstraintLinear.hpp

    r997 r998  
    2121 
    2222  /** Fills gradient.  If Linear then solution may be NULL,
    23       also returns value of function
     23      also returns true value of function and offset so we can use x not deltaX in constraint
    2424      If refresh is false then uses last solution
    2525      Uses model for scaling
     
    3030                       double * gradient,
    3131                       double & functionValue ,
     32                       double & offset,
    3233                       bool refresh=true) const ;
    3334  /// Resize constraint
     
    4142   */
    4243  virtual int markNonlinear(char * which) const ;
     44  /** Given a zeroed array sets possible nonzero coefficients to 1.
     45      Returns number of nonzeros
     46   */
     47  virtual int markNonzero(char * which) const;
    4348  //@}
    4449 
     
    5055 
    5156  /// Constructor from constraint
    52   ClpConstraintLinear(int numberCoefficients, int numberColumns,
     57  ClpConstraintLinear(int row, int numberCoefficients, int numberColumns,
    5358                      const int * column, const double * element);
    5459 
  • branches/devel/Clp/src/ClpSimplex.cpp

    r990 r998  
    48364836{
    48374837  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
     4838}
     4839/* Solves problem with nonlinear constraints using SLP - may be used as crash
     4840   for other algorithms when number of iterations small.
     4841   Also exits if all problematical variables are changing
     4842   less than deltaTolerance
     4843*/
     4844int
     4845ClpSimplex::nonlinearSLP(int numberConstraints, ClpConstraint ** constraints,
     4846                   int numberPasses,double deltaTolerance)
     4847{
     4848  return ((ClpSimplexNonlinear *) this)->primalSLP(numberConstraints,constraints,numberPasses,deltaTolerance);
    48384849}
    48394850// Solves non-linear using reduced gradient
  • branches/devel/Clp/src/ClpSimplex.hpp

    r990 r998  
    2626class CoinWarmStartBasis;
    2727class ClpDisasterHandler;
     28class ClpConstraint;
    2829
    2930/** This solves LPs using the simplex method
     
    249250  */
    250251  int nonlinearSLP(int numberPasses,double deltaTolerance);
     252  /** Solves problem with nonlinear constraints using SLP - may be used as crash
     253      for other algorithms when number of iterations small.
     254      Also exits if all problematical variables are changing
     255      less than deltaTolerance
     256  */
     257  int nonlinearSLP(int numberConstraints, ClpConstraint ** constraints,
     258                   int numberPasses,double deltaTolerance);
    251259  /** Solves using barrier (assumes you have good cholesky factor code).
    252260      Does crossover to simplex if asked*/
  • branches/devel/Clp/src/ClpSimplexNonlinear.cpp

    r913 r998  
    1111#include "ClpNonLinearCost.hpp"
    1212#include "ClpLinearObjective.hpp"
     13#include "ClpConstraint.hpp"
    1314#include "ClpQuadraticObjective.hpp"
    1415#include "CoinPackedMatrix.hpp"
     
    33103311  return 0;
    33113312}
     3313/* Primal algorithm for nonlinear constraints
     3314   Using a semi-trust region approach as for pooling problem
     3315   This is in because I have it lying around
     3316   
     3317*/
     3318int
     3319ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints,
     3320                               int numberPasses, double deltaTolerance)
     3321{
     3322  if (!numberConstraints) {
     3323    // no nonlinear part
     3324    return ClpSimplex::primal(0);
     3325  }
     3326  // check all matrix for odd rows is empty
     3327  int iConstraint;
     3328  int numberBad=0;
     3329  CoinPackedMatrix * columnCopy = matrix();
     3330  // Get a row copy in standard format
     3331  CoinPackedMatrix copy;
     3332  copy.reverseOrderedCopyOf(*columnCopy);
     3333  // get matrix data pointers
     3334  //const int * column = copy.getIndices();
     3335  //const CoinBigIndex * rowStart = copy.getVectorStarts();
     3336  const int * rowLength = copy.getVectorLengths();
     3337  //const double * elementByRow = copy.getElements();
     3338  int numberArtificials=0;
     3339  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3340    ClpConstraint * constraint = constraints[iConstraint];
     3341    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    }
     3364  }
     3365  if (numberBad)
     3366    return numberBad;
     3367  ClpSimplex newModel(*this);
     3368  int numberColumns2 = numberColumns_;
     3369  double penalty=1.0e9;
     3370  if (numberArtificials) {
     3371    numberColumns2 += numberArtificials;
     3372    int * addStarts = new int [numberArtificials+1];
     3373    int * addRow = new int[numberArtificials];
     3374    double * addElement = new double[numberArtificials];
     3375    addStarts[0]=0;
     3376    double * addCost = new double [numberArtificials];
     3377    numberArtificials=0;
     3378    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3379      ClpConstraint * constraint = constraints[iConstraint];
     3380      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        }
     3396      }
     3397    }
     3398    newModel.addColumns(numberArtificials,NULL,NULL,addCost,
     3399                       addStarts,addRow,addElement);
     3400    delete [] addStarts;
     3401    delete [] addRow;
     3402    delete [] addElement;
     3403    delete [] addCost;
     3404  }
     3405  // find nonlinear columns
     3406  int * listNonLinearColumn = new int [numberColumns_];
     3407  char * mark = new char[numberColumns_];
     3408  memset(mark,0,numberColumns_);
     3409  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3410    ClpConstraint * constraint = constraints[iConstraint];
     3411    constraint->markNonlinear(mark);
     3412  }
     3413  int iColumn;
     3414  int numberNonLinearColumns=0;
     3415  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3416    if (mark[iColumn])
     3417      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
     3418  }
     3419  // build row copy of matrix with space for nonzeros
     3420  // Get column copy
     3421  columnCopy = matrix();
     3422  copy.reverseOrderedCopyOf(*columnCopy);
     3423  // get matrix data pointers
     3424  const int * column = copy.getIndices();
     3425  const CoinBigIndex * rowStart = copy.getVectorStarts();
     3426  rowLength = copy.getVectorLengths();
     3427  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  }
     3433  numberExtra +=copy.getNumElements();
     3434  CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1];
     3435  int * newColumn = new int[numberExtra];
     3436  double * newElement = new double[numberExtra];
     3437  newStarts[0]=0;
     3438  int * backRow = new int [numberRows_];
     3439  int iRow;
     3440  for (iRow=0;iRow<numberRows_;iRow++)
     3441    backRow[iRow]=-1;
     3442  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3443    ClpConstraint * constraint = constraints[iConstraint];
     3444    iRow = constraint->rowNumber();
     3445    backRow[iRow]=iConstraint;
     3446  }
     3447  numberExtra=0;
     3448  for (iRow=0;iRow<numberRows_;iRow++) {
     3449    if (backRow[iRow]<0) {
     3450      // copy normal
     3451      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3452           j++) {
     3453        newColumn[numberExtra]=column[j];
     3454        newElement[numberExtra++] = elementByRow[j];
     3455      }
     3456    } else {
     3457      // all possible
     3458      memset(mark,0,numberColumns_);
     3459      ClpConstraint * constraint = constraints[backRow[iRow]];
     3460      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    }
     3469    newStarts[iRow+1]=numberExtra;
     3470  }
     3471  delete [] backRow;
     3472  CoinPackedMatrix saveMatrix(false,numberColumns2,numberRows_,
     3473                              numberExtra,newElement,newColumn,newStarts,NULL,0.0,0.0);
     3474  delete [] newStarts;
     3475  delete [] newColumn;
     3476  delete [] newElement;
     3477  delete [] mark;
     3478  int numberRows = newModel.numberRows();
     3479  double * columnLower = newModel.columnLower();
     3480  double * columnUpper = newModel.columnUpper();
     3481  double * objective = newModel.objective();
     3482  double * rowLower = newModel.rowLower();
     3483  double * rowUpper = newModel.rowUpper();
     3484  double * solution = newModel.primalColumnSolution();
     3485  int jNon;
     3486  int * last[3];
     3487 
     3488  double * trust = new double[numberNonLinearColumns];
     3489  double * trueLower = new double[numberNonLinearColumns];
     3490  double * trueUpper = new double[numberNonLinearColumns];
     3491  double objectiveOffset;
     3492  double objectiveOffset2;
     3493  getDblParam(ClpObjOffset,objectiveOffset);
     3494  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3495    iColumn=listNonLinearColumn[jNon];
     3496    double upper = columnUpper[iColumn];
     3497    double lower = columnLower[iColumn];
     3498    if (upper>1.0e10)
     3499      upper = max(1000.0,10*solution[iColumn]);
     3500    //columnUpper[iColumn]=upper;
     3501    trust[jNon]=0.05*(1.0+upper-lower);
     3502    trueLower[jNon]=lower;
     3503    //trueUpper[jNon]=upper;
     3504    trueUpper[jNon]=columnUpper[iColumn];
     3505  }
     3506  bool tryFix=false;
     3507  int iPass;
     3508  double lastObjective=1.0e31;
     3509  double * saveSolution = new double [numberColumns2+numberRows];
     3510  char * saveStatus = new char [numberColumns2+numberRows];
     3511  double targetDrop=1.0e31;
     3512  //double objectiveOffset2 = objectiveOffset;
     3513  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
     3514  for (iPass=0;iPass<3;iPass++) {
     3515    last[iPass]=new int[numberNonLinearColumns];
     3516    for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3517      last[iPass][jNon]=0;
     3518  }
     3519  int numberZeroPasses=0;
     3520  bool zeroTargetDrop=false;
     3521  double * gradient = new double [numberColumns_];
     3522#define SMALL_FIX 0.0
     3523  for (iPass=0;iPass<numberPasses;iPass++) {
     3524    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3525      iColumn=listNonLinearColumn[jNon];
     3526      if (solution[iColumn]<trueLower[jNon])
     3527        solution[iColumn]=trueLower[jNon];
     3528      else if (solution[iColumn]>trueUpper[jNon])
     3529        solution[iColumn]=trueUpper[jNon];
     3530      if (iPass) {
     3531        columnLower[iColumn]=max(solution[iColumn]
     3532                                 -trust[jNon],
     3533                                 trueLower[jNon]);
     3534        if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
     3535          columnLower[iColumn]=SMALL_FIX;
     3536        columnUpper[iColumn]=min(solution[iColumn]
     3537                                 +trust[jNon],
     3538                                 trueUpper[jNon]);
     3539        if (!trueLower[jNon])
     3540          columnUpper[iColumn] = max(columnUpper[iColumn],
     3541                                     columnLower[iColumn]+SMALL_FIX);
     3542        if (!trueLower[jNon]&&tryFix&&
     3543            columnLower[iColumn]==SMALL_FIX&&
     3544            columnUpper[iColumn]<3.0*SMALL_FIX) {
     3545          columnLower[iColumn]=0.0;
     3546          solution[iColumn]=0.0;
     3547          columnUpper[iColumn]=0.0;
     3548          printf("fixing %d\n",iColumn);
     3549        }
     3550      } else {
     3551#if 1
     3552        // fix nonlinear variables
     3553        columnLower[iColumn]=solution[iColumn];
     3554        columnUpper[iColumn]=solution[iColumn];
     3555#endif
     3556      }
     3557    }
     3558    // redo matrix
     3559    double offset;
     3560    CoinPackedMatrix newMatrix(saveMatrix);
     3561    // get matrix data pointers
     3562    column = newMatrix.getIndices();
     3563    rowStart = newMatrix.getVectorStarts();
     3564    rowLength = newMatrix.getVectorLengths();
     3565    double * changeableElement = newMatrix.getMutableElements();
     3566    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3567      ClpConstraint * constraint = constraints[iConstraint];
     3568      int iRow = constraint->rowNumber();
     3569      double functionValue;
     3570      int numberErrors = constraint->gradient(&newModel,solution,gradient,functionValue,offset);
     3571      assert (!numberErrors);
     3572      if (iRow>=0) {
     3573        if (numberRows<50)
     3574          printf("For row %d current value is %g (activity %g) , dual is %g\n",iRow,offset,
     3575                 newModel.primalRowSolution()[iRow],
     3576                 newModel.dualRowSolution()[iRow]);
     3577        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3578             j++) {
     3579          int iColumn = column[j];
     3580          changeableElement[j] = gradient[iColumn];
     3581          gradient[iColumn]=0.0;
     3582        }
     3583        for (int k=0;k<numberColumns_;k++)
     3584          assert (!gradient[k]);
     3585        if (rowLower_[iRow]>-1.0e20)
     3586          rowLower[iRow] = rowLower_[iRow] - offset;
     3587        if (rowUpper_[iRow]<1.0e20)
     3588          rowUpper[iRow] = rowUpper_[iRow] - offset;
     3589      } else {
     3590        memcpy(objective,gradient,numberColumns_*sizeof(double));
     3591        printf("objective offset %g\n",offset);
     3592        objectiveOffset2 =objectiveOffset-offset;
     3593        newModel.setObjectiveOffset(objectiveOffset2);
     3594      }
     3595    }
     3596    // Replace matrix
     3597    // Get a column copy in standard format
     3598    CoinPackedMatrix * columnCopy = new CoinPackedMatrix();;
     3599    columnCopy->reverseOrderedCopyOf(newMatrix);
     3600    newModel.replaceMatrix(columnCopy,true);
     3601    double objValue=0.0;
     3602    double infValue=0.0;
     3603    {
     3604      // get matrix data pointers
     3605      const int * row = columnCopy->getIndices();
     3606      const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     3607      const int * columnLength = columnCopy->getVectorLengths();
     3608      const double * elementByColumn = columnCopy->getElements();
     3609      int iRow,iColumn;
     3610      double * accumulate = new double [numberRows];
     3611      //const double * solution = newModel.primalColumnSolution();
     3612      double * partSolution = new double [numberColumns2];
     3613      memcpy(partSolution,solution,numberColumns_*sizeof(double));
     3614      memset(partSolution+numberColumns_,0,
     3615             (numberColumns2-numberColumns_)*sizeof(double));
     3616      memset(accumulate,0,numberRows*sizeof(double));
     3617      //newModel.matrix()->times(partSolution,accumulate);
     3618      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3619        if (solution[iColumn]) {
     3620          double value=solution[iColumn];
     3621          objValue += objective[iColumn]*solution[iColumn];
     3622          int j;
     3623          for (j=columnStart[iColumn];
     3624               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     3625            int iRow=row[j];
     3626            accumulate[iRow] += value*elementByColumn[j];
     3627          }
     3628        }
     3629      }
     3630      for (iRow=0;iRow<numberRows;iRow++) {
     3631        if (accumulate[iRow]<rowLower[iRow]-1.0e-5) {
     3632          double under = rowLower[iRow]-accumulate[iRow];
     3633          infValue += under;
     3634#if 0
     3635          printf("Row %d l=%g u=%g a=%g",iRow,rowLower[iRow],
     3636                 rowUpper[iRow],accumulate[iRow]);
     3637          printf(" ** under by %g",under);
     3638          printf("\n");
     3639#endif
     3640        } else if (accumulate[iRow]>rowUpper[iRow]+1.0e-5) {
     3641          double over = accumulate[iRow]-rowUpper[iRow];
     3642          infValue += over;
     3643#if 0
     3644          printf("Row %d l=%g u=%g a=%g",iRow,rowLower[iRow],
     3645                 rowUpper[iRow],accumulate[iRow]);
     3646          printf(" ** over by %g",over);
     3647          printf("\n");
     3648#endif
     3649        }
     3650      }
     3651      if (infValue)
     3652        printf("Sum infeasibilities %g ",infValue);
     3653      if (infValue<0.1)
     3654        infValue=0.0;
     3655      infValue *= penalty;
     3656      if (infValue)
     3657        printf("Infeasible obj %g ",infValue);
     3658      delete [] accumulate; 
     3659      delete [] partSolution;
     3660    }
     3661    if (objectiveOffset2)
     3662      printf("offset2 %g\n",objectiveOffset2);
     3663    objValue = objValue + infValue -objectiveOffset2;
     3664    printf("True objective %g\n",objValue);
     3665    if (iPass) {
     3666      double drop = lastObjective-objValue;
     3667      std::cout<<"True drop was "<<drop<<std::endl;
     3668      if (drop<-0.05*fabs(objValue)-1.0e-4) {
     3669        // pretty bad - go back and halve
     3670        memcpy(solution,saveSolution,numberColumns2*sizeof(double));
     3671        memcpy(newModel.primalRowSolution(),saveSolution+numberColumns2,
     3672               numberRows*sizeof(double));
     3673        memcpy(newModel.statusArray(),saveStatus,
     3674               numberColumns2+numberRows);
     3675        for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3676          if (trust[jNon]>0.1)
     3677            trust[jNon] *= 0.5;
     3678          else
     3679            trust[jNon] *= 0.9;
     3680       
     3681        printf("** Halving trust\n");
     3682        objValue=lastObjective;
     3683        continue;
     3684      } else if ((iPass%25)==-1) {
     3685        for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3686          trust[jNon] *= 2.0;
     3687        printf("** Doubling trust\n");
     3688      }
     3689      int * temp=last[2];
     3690      last[2]=last[1];
     3691      last[1]=last[0];
     3692      last[0]=temp;
     3693      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3694        iColumn=listNonLinearColumn[jNon];
     3695        double change = solution[iColumn]-saveSolution[iColumn];
     3696        if (change<-1.0e-5) {
     3697          if (fabs(change+trust[jNon])<1.0e-5)
     3698            temp[jNon]=-1;
     3699          else
     3700            temp[jNon]=-2;
     3701        } else if(change>1.0e-5) {
     3702          if (fabs(change-trust[jNon])<1.0e-5)
     3703            temp[jNon]=1;
     3704          else
     3705            temp[jNon]=2;
     3706        } else {
     3707          temp[jNon]=0;
     3708        }
     3709      }
     3710      double maxDelta=0.0;
     3711      double smallestTrust=1.0e31;
     3712      double smallestNonLinearGap=1.0e31;
     3713      double smallestGap=1.0e31;
     3714      bool increasing=false;
     3715      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3716        double gap = columnUpper[iColumn]-columnLower[iColumn];
     3717        assert (gap>=0.0);
     3718        if (gap)
     3719          smallestGap = min(smallestGap,gap);
     3720      }
     3721      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3722        iColumn=listNonLinearColumn[jNon];
     3723        double gap = columnUpper[iColumn]-columnLower[iColumn];
     3724        assert (gap>=0.0);
     3725        if (gap) {
     3726          smallestNonLinearGap = min(smallestNonLinearGap,gap);
     3727          if (gap<1.0e-7&&iPass==1) {
     3728            printf("Small gap %d %d %g %g %g\n",
     3729                   jNon,iColumn,columnLower[iColumn],columnUpper[iColumn],
     3730                   gap);
     3731            //trueUpper[jNon]=trueLower[jNon];
     3732            //columnUpper[iColumn]=columnLower[iColumn];
     3733          }
     3734        }
     3735        maxDelta = max(maxDelta,
     3736                       fabs(solution[iColumn]-saveSolution[iColumn]));
     3737        if (last[0][jNon]*last[1][jNon]<0) {
     3738          // halve
     3739          if (trust[jNon]>1.0)
     3740            trust[jNon] *= 0.5;
     3741          else
     3742            trust[jNon] *= 0.7;
     3743        } else {
     3744          // ? only increase if +=1 ?
     3745          if (last[0][jNon]==last[1][jNon]&&
     3746              last[0][jNon]==last[2][jNon]&&
     3747              last[0][jNon]) {
     3748            trust[jNon] *= 1.8;
     3749            increasing=true;
     3750          }
     3751        }
     3752        smallestTrust = min(smallestTrust,trust[jNon]);
     3753      }
     3754      std::cout<<"largest delta is "<<maxDelta
     3755               <<", smallest trust is "<<smallestTrust
     3756               <<", smallest gap is "<<smallestGap
     3757               <<", smallest nonlinear gap is "<<smallestNonLinearGap
     3758               <<std::endl;
     3759      if (maxDelta<deltaTolerance&&!increasing&&iPass>100) {
     3760        numberZeroPasses++;
     3761        if (numberZeroPasses==3) {
     3762          if (tryFix) {
     3763            std::cout<<"Exiting"<<std::endl;
     3764            break;
     3765          } else {
     3766            tryFix=true;
     3767            for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3768              trust[jNon] = max(trust[jNon],1.0e-1);
     3769            numberZeroPasses=0;
     3770          }
     3771        }
     3772      } else {
     3773        numberZeroPasses=0;
     3774      }
     3775    }
     3776    memcpy(saveSolution,solution,numberColumns2*sizeof(double));
     3777    memcpy(saveSolution+numberColumns2,newModel.primalRowSolution(),
     3778           numberRows*sizeof(double));
     3779    memcpy(saveStatus,newModel.statusArray(),
     3780           numberColumns2+numberRows);
     3781   
     3782    targetDrop=0.0;
     3783    if (iPass) {
     3784      // get reduced costs
     3785      const double * pi = newModel.dualRowSolution();
     3786      newModel.matrix()->transposeTimes(pi,
     3787                                        newModel.dualColumnSolution());
     3788      double * r = newModel.dualColumnSolution();
     3789      for (iColumn=0;iColumn<numberColumns_;iColumn++)
     3790        r[iColumn] = objective[iColumn]-r[iColumn];
     3791      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3792        iColumn=listNonLinearColumn[jNon];
     3793        double dj = r[iColumn];
     3794        if (dj<-1.0e-6) {
     3795          double drop = -dj*(columnUpper[iColumn]-solution[iColumn]);
     3796          //double upper = min(trueUpper[jNon],solution[iColumn]+0.1);
     3797          //double drop2 = -dj*(upper-solution[iColumn]);
     3798#if 0
     3799          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2&&iPass>100))
     3800            printf("Big drop %d %g %g %g %g T %g\n",
     3801                   iColumn,columnLower[iColumn],solution[iColumn],
     3802                   columnUpper[iColumn],dj,trueUpper[jNon]);
     3803#endif
     3804          targetDrop += drop;
     3805          if (dj<-1.0e-1&&trust[jNon]<1.0e-3
     3806              &&trueUpper[jNon]-solution[iColumn]>1.0e-2) {
     3807            trust[jNon] *= 1.5;
     3808            //printf("Increasing trust on %d to %g\n",
     3809            //     iColumn,trust[jNon]);
     3810          }
     3811        } else if (dj>1.0e-6) {
     3812          double drop = -dj*(columnLower[iColumn]-solution[iColumn]);
     3813          //double lower = max(trueLower[jNon],solution[iColumn]-0.1);
     3814          //double drop2 = -dj*(lower-solution[iColumn]);
     3815#if 0
     3816          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2))
     3817            printf("Big drop %d %g %g %g %g T %g\n",
     3818                   iColumn,columnLower[iColumn],solution[iColumn],
     3819                   columnUpper[iColumn],dj,trueLower[jNon]);
     3820#endif
     3821          targetDrop += drop;
     3822          if (dj>1.0e-1&&trust[jNon]<1.0e-3
     3823              &&solution[iColumn]-trueLower[jNon]>1.0e-2) {
     3824            trust[jNon] *= 1.5;
     3825            printf("Increasing trust on %d to %g\n",
     3826                   iColumn,trust[jNon]);
     3827          }
     3828        }
     3829      }
     3830    }
     3831    std::cout<<"Pass - "<<iPass
     3832             <<", target drop is "<<targetDrop
     3833             <<std::endl;
     3834    if (iPass>1&&targetDrop<1.0e-5&&zeroTargetDrop)
     3835      break;
     3836    if (iPass>1&&targetDrop<1.0e-5)
     3837      zeroTargetDrop = true;
     3838    else
     3839      zeroTargetDrop = false;
     3840    //if (iPass==5)
     3841    //newModel.setLogLevel(63);
     3842    lastObjective = objValue;
     3843    // take out when ClpPackedMatrix changed
     3844    //newModel.scaling(false);
     3845#if 0
     3846    CoinMpsIO writer;
     3847    writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
     3848                      newModel.getColLower(), newModel.getColUpper(),
     3849                      newModel.getObjCoefficients(),
     3850                      (const char*) 0 /*integrality*/,
     3851                      newModel.getRowLower(), newModel.getRowUpper(),
     3852                      NULL,NULL);
     3853    writer.writeMps("xx.mps");
     3854#endif
     3855    if (iPass<3) {
     3856      newModel.primal(1);
     3857      abort();
     3858    } else {
     3859      abort();
     3860    }
     3861    if (newModel.status()==1) {
     3862      // Infeasible!
     3863      newModel.allSlackBasis();
     3864      newModel.primal();
     3865      assert(!newModel.status());
     3866    }
     3867    double sumInfeas=0.0;
     3868    int numberInfeas=0;
     3869    for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
     3870      if (solution[iColumn]>1.0e-8) {
     3871        numberInfeas++;
     3872        sumInfeas += solution[iColumn];
     3873      }
     3874    }
     3875    printf("%d infeasibilities - summing to %g\n",
     3876           numberInfeas,sumInfeas);
     3877    //model.writeMps("xx");
     3878  }
     3879  delete [] saveSolution;
     3880  delete [] saveStatus;
     3881  for (iPass=0;iPass<3;iPass++)
     3882    delete [] last[iPass];
     3883  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3884    iColumn=listNonLinearColumn[jNon];
     3885    columnLower[iColumn]=trueLower[jNon];
     3886    columnUpper[iColumn]=trueUpper[jNon];
     3887  }
     3888  delete [] trust;
     3889  delete [] trueUpper;
     3890  delete [] trueLower;
     3891  // Simplest way to get true row activity ?
     3892  double * rowActivity = newModel.primalRowSolution();
     3893  for (iRow=0;iRow<numberRows;iRow++) {
     3894    double difference;
     3895    if (fabs(rowLower_[iRow])<fabs(rowUpper_[iRow]))
     3896      difference = rowLower_[iRow]-rowLower[iRow];
     3897    else
     3898      difference = rowUpper_[iRow]-rowUpper[iRow];
     3899    rowLower[iRow]=rowLower_[iRow];
     3900    rowUpper[iRow]=rowUpper_[iRow];
     3901    if (difference) {
     3902      if (numberRows<50)
     3903        printf("For row %d activity changes from %g to %g\n",
     3904               iRow,rowActivity[iRow],rowActivity[iRow]+difference);
     3905      rowActivity[iRow]+= difference;
     3906    }
     3907  }
     3908  delete [] saveSolution;
     3909  delete [] saveStatus;
     3910  for (iPass=0;iPass<3;iPass++)
     3911    delete [] last[iPass];
     3912  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3913    iColumn=listNonLinearColumn[jNon];
     3914    columnLower[iColumn]=trueLower[jNon];
     3915    columnUpper[iColumn]=trueUpper[jNon];
     3916  }
     3917  delete [] trust;
     3918  delete [] trueUpper;
     3919  delete [] trueLower;
     3920  delete [] listNonLinearColumn;
     3921  delete [] gradient;
     3922  printf("solution still in newModel - do objective etc!\n");
     3923  objectiveValue_=newModel.objectiveValue();
     3924  numberIterations_ =newModel.numberIterations();
     3925  problemStatus_ =newModel.problemStatus();
     3926  secondaryStatus_ =newModel.secondaryStatus();
     3927  memcpy(columnActivity_,newModel.primalColumnSolution(),numberColumns_*sizeof(double));
     3928  // should do status region
     3929  CoinZeroN(rowActivity_,numberRows_);
     3930  matrix_->times(1.0,columnActivity_,rowActivity_);
     3931  abort();
     3932  return 0;
     3933}
  • branches/devel/Clp/src/ClpSimplexNonlinear.hpp

    r754 r998  
    1313class ClpNonlinearInfo;
    1414class ClpQuadraticObjective;
     15class ClpConstraint;
    1516
    1617#include "ClpSimplexPrimal.hpp"
     
    4243  */
    4344  int primalSLP(int numberPasses, double deltaTolerance);
     45  /** Primal algorithm for nonlinear constraints
     46      Using a semi-trust region approach as for pooling problem
     47      This is in because I have it lying around
     48
     49  */
     50  int primalSLP(int numberConstraints, ClpConstraint ** constraints,
     51                int numberPasses, double deltaTolerance);
    4452
    4553  /** Creates direction vector.  note longArray is long enough
  • branches/devel/Clp/src/Makefile.in

    r996 r998  
    6969libClp_la_LIBADD =
    7070am_libClp_la_OBJECTS = ClpCholeskyBase.lo ClpCholeskyDense.lo \
    71         ClpCholeskyUfl.lo ClpCholeskyWssmp.lo ClpCholeskyWssmpKKT.lo \
    72         ClpConstraint.lo ClpConstraintLinear.lo Clp_C_Interface.lo \
    73         ClpDualRowDantzig.lo ClpDualRowPivot.lo ClpDualRowSteepest.lo \
    74         ClpDummyMatrix.lo ClpDynamicExampleMatrix.lo \
    75         ClpDynamicMatrix.lo ClpEventHandler.lo ClpFactorization.lo \
    76         ClpGubDynamicMatrix.lo ClpGubMatrix.lo ClpHelperFunctions.lo \
    77         ClpInterior.lo ClpLinearObjective.lo ClpMatrixBase.lo \
    78         ClpMessage.lo ClpModel.lo ClpNetworkBasis.lo \
    79         ClpNetworkMatrix.lo ClpNonLinearCost.lo ClpObjective.lo \
    80         ClpPackedMatrix.lo ClpPlusMinusOneMatrix.lo \
    81         ClpPredictorCorrector.lo ClpPresolve.lo \
    82         ClpPrimalColumnDantzig.lo ClpPrimalColumnPivot.lo \
    83         ClpPrimalColumnSteepest.lo ClpQuadraticObjective.lo \
    84         ClpSimplex.lo ClpSimplexDual.lo ClpSimplexNonlinear.lo \
    85         ClpSimplexOther.lo ClpSimplexPrimal.lo ClpSolve.lo Idiot.lo \
    86         IdiSolve.lo
     71        ClpCholeskyUfl.lo ClpConstraint.lo ClpConstraintLinear.lo \
     72        Clp_C_Interface.lo ClpDualRowDantzig.lo ClpDualRowPivot.lo \
     73        ClpDualRowSteepest.lo ClpDummyMatrix.lo \
     74        ClpDynamicExampleMatrix.lo ClpDynamicMatrix.lo \
     75        ClpEventHandler.lo ClpFactorization.lo ClpGubDynamicMatrix.lo \
     76        ClpGubMatrix.lo ClpHelperFunctions.lo ClpInterior.lo \
     77        ClpLinearObjective.lo ClpMatrixBase.lo ClpMessage.lo \
     78        ClpModel.lo ClpNetworkBasis.lo ClpNetworkMatrix.lo \
     79        ClpNonLinearCost.lo ClpObjective.lo ClpPackedMatrix.lo \
     80        ClpPlusMinusOneMatrix.lo ClpPredictorCorrector.lo \
     81        ClpPresolve.lo ClpPrimalColumnDantzig.lo \
     82        ClpPrimalColumnPivot.lo ClpPrimalColumnSteepest.lo \
     83        ClpQuadraticObjective.lo ClpSimplex.lo ClpSimplexDual.lo \
     84        ClpSimplexNonlinear.lo ClpSimplexOther.lo ClpSimplexPrimal.lo \
     85        ClpSolve.lo Idiot.lo IdiSolve.lo
    8786libClp_la_OBJECTS = $(am_libClp_la_OBJECTS)
    8887binPROGRAMS_INSTALL = $(INSTALL_PROGRAM)
     
    276275        ClpCholeskyDense.cpp ClpCholeskyDense.hpp \
    277276        ClpCholeskyUfl.cpp ClpCholeskyUfl.hpp \
    278         ClpCholeskyWssmp.cpp ClpCholeskyWssmp.hpp \
    279         ClpCholeskyWssmpKKT.cpp ClpCholeskyWssmpKKT.hpp \
    280277        ClpConstraint.cpp ClpConstraint.hpp \
    281278        ClpConstraintLinear.cpp ClpConstraintLinear.hpp \
     
    514511@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyDense.Plo@am__quote@
    515512@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyUfl.Plo@am__quote@
    516 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyWssmp.Plo@am__quote@
    517 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyWssmpKKT.Plo@am__quote@
    518513@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpConstraint.Plo@am__quote@
    519514@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpConstraintLinear.Plo@am__quote@
Note: See TracChangeset for help on using the changeset viewer.