Changeset 602


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

nonlinear

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

Legend:

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

    r596 r602  
    88#endif
    99#endif
     10#include "CoinTime.hpp"
     11
     12#include "CoinHelperFunctions.hpp"
     13#include "CoinModel.hpp"
     14#include "ClpSimplex.hpp"
     15// returns jColumn (-2 if linear term, -1 if unknown) and coefficient
     16static
     17int decodeBit(char * phrase, char * & nextPhrase, double & coefficient, bool ifFirst, const CoinModel & model)
     18{
     19  char * pos = phrase;
     20  // may be leading - (or +)
     21  char * pos2 = pos;
     22  double value=1.0;
     23  if (*pos2=='-'||*pos2=='+')
     24    pos2++;
     25  // next terminator * or + or -
     26  while (*pos2) {
     27    if (*pos2=='*') {
     28      break;
     29    } else if (*pos2=='-'||*pos2=='+') {
     30      if (pos2==pos||*(pos2-1)!='e')
     31        break;
     32    }
     33    pos2++;
     34  }
     35  // if * must be number otherwise must be name
     36  if (*pos2=='*') {
     37    char * pos3 = pos;
     38    while (pos3!=pos2) {
     39      char x = *pos3;
     40      pos3++;
     41      assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-'||x=='e');
     42    }
     43    char saved = *pos2;
     44    *pos2='\0';
     45    value = atof(pos);
     46    *pos2=saved;
     47    // and down to next
     48    pos2++;
     49    pos=pos2;
     50    while (*pos2) {
     51      if (*pos2=='-'||*pos2=='+')
     52        break;
     53      pos2++;
     54    }
     55  }
     56  char saved = *pos2;
     57  *pos2='\0';
     58  // now name
     59  // might have + or -
     60  if (*pos=='+') {
     61    pos++;
     62  } else if (*pos=='-') {
     63    pos++;
     64    assert (value==1.0);
     65    value = - value;
     66  }
     67  int jColumn = model.column(pos);
     68  // must be column unless first when may be linear term
     69  if (jColumn<0) {
     70    if (ifFirst) {
     71      char * pos3 = pos;
     72      while (pos3!=pos2) {
     73        char x = *pos3;
     74        pos3++;
     75        assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-'||x=='e');
     76      }
     77      assert(*pos2=='\0');
     78      // keep possible -
     79      value = value * atof(pos);
     80      jColumn=-2;
     81    } else {
     82      // bad
     83      *pos2=saved;
     84      printf("bad nonlinear term %s\n",phrase);
     85      abort();
     86    }
     87  }
     88  *pos2=saved;
     89  pos=pos2;
     90  coefficient=value;
     91  nextPhrase = pos;
     92  return jColumn;
     93}
    1094#ifdef COIN_HAS_LINK
    1195#include <cassert>
     
    1599#endif
    16100#include "CbcLinked.hpp"
    17 #include "CoinTime.hpp"
    18 
    19 #include "CoinHelperFunctions.hpp"
    20101#include "CoinIndexedVector.hpp"
    21102#include "CoinMpsIO.hpp"
    22 #include "CoinModel.hpp"
    23 #include "ClpSimplex.hpp"
    24103//#include "OsiSolverLink.hpp"
    25104//#include "OsiBranchLink.hpp"
     
    828907  if (up>maximumValue)
    829908    solver->setColUpper(column,maximumValue);
    830 }
    831 // returns jColumn (-2 if linear term, -1 if unknown) and coefficient
    832 static
    833 int decodeBit(char * phrase, char * & nextPhrase, double & coefficient, bool ifFirst, const CoinModel & model)
    834 {
    835   char * pos = phrase;
    836   // may be leading - (or +)
    837   char * pos2 = pos;
    838   double value=1.0;
    839   if (*pos2=='-'||*pos2=='+')
    840     pos2++;
    841   // next terminator * or + or -
    842   while (*pos2) {
    843     if (*pos2=='*') {
    844       break;
    845     } else if (*pos2=='-'||*pos2=='+') {
    846       if (pos2==pos||*(pos2-1)!='e')
    847         break;
    848     }
    849     pos2++;
    850   }
    851   // if * must be number otherwise must be name
    852   if (*pos2=='*') {
    853     char * pos3 = pos;
    854     while (pos3!=pos2) {
    855       char x = *pos3;
    856       pos3++;
    857       assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-'||x=='e');
    858     }
    859     char saved = *pos2;
    860     *pos2='\0';
    861     value = atof(pos);
    862     *pos2=saved;
    863     // and down to next
    864     pos2++;
    865     pos=pos2;
    866     while (*pos2) {
    867       if (*pos2=='-'||*pos2=='+')
    868         break;
    869       pos2++;
    870     }
    871   }
    872   char saved = *pos2;
    873   *pos2='\0';
    874   // now name
    875   // might have + or -
    876   if (*pos=='+') {
    877     pos++;
    878   } else if (*pos=='-') {
    879     pos++;
    880     assert (value==1.0);
    881     value = - value;
    882   }
    883   int jColumn = model.column(pos);
    884   // must be column unless first when may be linear term
    885   if (jColumn<0) {
    886     if (ifFirst) {
    887       char * pos3 = pos;
    888       while (pos3!=pos2) {
    889         char x = *pos3;
    890         pos3++;
    891         assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-'||x=='e');
    892       }
    893       assert(*pos2=='\0');
    894       // keep possible -
    895       value = value * atof(pos);
    896       jColumn=-2;
    897     } else {
    898       // bad
    899       *pos2=saved;
    900       printf("bad nonlinear term %s\n",phrase);
    901       abort();
    902     }
    903   }
    904   *pos2=saved;
    905   pos=pos2;
    906   coefficient=value;
    907   nextPhrase = pos;
    908   return jColumn;
    909909}
    910910void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds,int logLevel)
     
    64856485}
    64866486#endif
     6487#include "ClpConstraint.hpp"
     6488#include "ClpConstraintLinear.hpp"
     6489/* Return an approximate solution to a CoinModel.
     6490    Lots of bounds may be odd to force a solution.
     6491    mode = 0 just tries to get a continuous solution
     6492*/
     6493ClpSimplex *
     6494approximateSolution(CoinModel & coinModel,
     6495                    int numberPasses, double deltaTolerance,
     6496                    int mode)
     6497{
     6498  // first check and set up arrays
     6499  int numberColumns = coinModel.numberColumns();
     6500  int numberRows = coinModel.numberRows();
     6501  // List of nonlinear rows
     6502  int * which = new int[numberRows+1];
     6503  bool testLinear=true;
     6504  int numberConstraints=0;
     6505  int iColumn;
     6506  // matrix etc will be changed
     6507  CoinModel coinModel2 = coinModel;
     6508  bool linearObjective=true;
     6509  int maximumQuadraticElements=0;
     6510  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     6511    // See if quadratic objective
     6512    const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
     6513    if (strcmp(expr,"Numeric")) {
     6514      linearObjective=false;
     6515      // check if value*x+-value*y....
     6516      assert (strlen(expr)<20000);
     6517      char temp[20000];
     6518      strcpy(temp,expr);
     6519      char * pos = temp;
     6520      bool ifFirst=true;
     6521      double linearTerm=0.0;
     6522      while (*pos) {
     6523        double value;
     6524        int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
     6525        // must be column unless first when may be linear term
     6526        if (jColumn>=0) {
     6527          maximumQuadraticElements++;
     6528        } else if (jColumn==-2) {
     6529          linearTerm = value;
     6530        } else {
     6531          printf("bad nonlinear term %s\n",temp);
     6532          abort();
     6533        }
     6534        ifFirst=false;
     6535      }
     6536    }
     6537  }
     6538  if (!linearObjective) {
     6539    // zero objective
     6540    for (iColumn=0;iColumn<numberColumns;iColumn++)
     6541      coinModel2.setObjective(iColumn,0.0);
     6542    which[numberConstraints++]=-1;
     6543  }
     6544  int iRow;
     6545  for (iRow=0;iRow<numberRows;iRow++) {   
     6546    int numberQuadratic=0;
     6547    bool linear=true;
     6548    CoinModelLink triple=coinModel.firstInRow(iRow);
     6549    while (triple.column()>=0) {
     6550      int iColumn = triple.column();
     6551      const char *  expr = coinModel.getElementAsString(iRow,iColumn);
     6552      if (strcmp("Numeric",expr)) {
     6553        linear=false;
     6554        // check if value*x+-value*y....
     6555        assert (strlen(expr)<20000);
     6556        char temp[20000];
     6557        strcpy(temp,expr);
     6558        char * pos = temp;
     6559        bool ifFirst=true;
     6560        double linearTerm=0.0;
     6561        while (*pos) {
     6562          double value;
     6563          int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
     6564          // must be column unless first when may be linear term
     6565          if (jColumn>=0) {
     6566            numberQuadratic++;
     6567          } else if (jColumn==-2) {
     6568            linearTerm = value;
     6569          } else {
     6570            printf("bad nonlinear term %s\n",temp);
     6571            abort();
     6572          }
     6573          ifFirst=false;
     6574        }
     6575      }
     6576      triple=coinModel.next(triple);
     6577    }
     6578    if (!linear||testLinear) {
     6579      CoinModelLink triple=coinModel.firstInRow(iRow);
     6580      while (triple.column()>=0) {
     6581        int iColumn = triple.column();
     6582        coinModel2.setElement(iRow,iColumn,0.0);
     6583        triple=coinModel.next(triple);
     6584      }
     6585      which[numberConstraints++]=iRow;
     6586      maximumQuadraticElements=CoinMax(maximumQuadraticElements,numberQuadratic);
     6587    }
     6588  }
     6589  ClpSimplex * model = new ClpSimplex();
     6590  // return if nothing
     6591  if (!numberConstraints) {
     6592    delete [] which;
     6593    model->loadProblem(coinModel);
     6594    model->dual();
     6595    return model;
     6596  }
     6597  // space for quadratic
     6598  CoinBigIndex * startQuadratic = new CoinBigIndex [numberColumns+1];
     6599  int * columnQuadratic = new int [maximumQuadraticElements];
     6600  double * elementQuadratic = new double [maximumQuadraticElements];
     6601  ClpConstraint ** constraints = new ClpConstraint * [numberConstraints];
     6602  double * linearTerm = new double [numberColumns];
     6603  int saveNumber=numberConstraints;
     6604  numberConstraints=0;
     6605  if (!linearObjective) {
     6606    int numberQuadratic=0;
     6607    CoinZeroN(linearTerm,numberColumns);
     6608    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     6609      startQuadratic[iColumn] = numberQuadratic;
     6610      // See if quadratic objective
     6611      const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
     6612      if (strcmp(expr,"Numeric")) {
     6613        // value*x*y
     6614        char temp[20000];
     6615        strcpy(temp,expr);
     6616        char * pos = temp;
     6617        bool ifFirst=true;
     6618        while (*pos) {
     6619          double value;
     6620          int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
     6621          // must be column unless first when may be linear term
     6622          if (jColumn>=0) {
     6623            columnQuadratic[numberQuadratic]=jColumn;
     6624            elementQuadratic[numberQuadratic++]=2.0*value; // convention
     6625          } else if (jColumn==-2) {
     6626            linearTerm[iColumn] = value;
     6627          } else {
     6628            printf("bad nonlinear term %s\n",temp);
     6629            abort();
     6630          }
     6631          ifFirst=false;
     6632        }
     6633      }
     6634    }
     6635    startQuadratic[numberColumns] = numberQuadratic;
     6636    // here we create ClpConstraint
     6637    abort();
     6638  }
     6639  int iConstraint;
     6640  for (iConstraint=0;iConstraint<saveNumber;iConstraint++) {
     6641    iRow = which[iConstraint];
     6642    if (iRow>=0) {
     6643      int numberQuadratic=0;
     6644      CoinZeroN(linearTerm,numberColumns);
     6645      CoinModelLink triple=coinModel.firstInRow(iRow);
     6646      while (triple.column()>=0) {
     6647        int iColumn = triple.column();
     6648        const char *  expr = coinModel.getElementAsString(iRow,iColumn);
     6649        if (strcmp("Numeric",expr)) {
     6650          // value*x*y
     6651          char temp[20000];
     6652          strcpy(temp,expr);
     6653          char * pos = temp;
     6654          bool ifFirst=true;
     6655          while (*pos) {
     6656            double value;
     6657            int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
     6658            // must be column unless first when may be linear term
     6659            if (jColumn>=0) {
     6660              columnQuadratic[numberQuadratic]=jColumn;
     6661              elementQuadratic[numberQuadratic++]=2.0*value; // convention
     6662            } else if (jColumn==-2) {
     6663              linearTerm[iColumn] = value;
     6664            } else {
     6665              printf("bad nonlinear term %s\n",temp);
     6666              abort();
     6667            }
     6668            ifFirst=false;
     6669          }
     6670        }
     6671        triple=coinModel.next(triple);
     6672      }
     6673      startQuadratic[numberColumns] = numberQuadratic;
     6674      // here we create ClpConstraint
     6675      if (testLinear) {
     6676        int n=0;
     6677        int * indices = new int[numberColumns];
     6678        for (int j=0;j<numberColumns;j++) {
     6679          if (linearTerm[j]) {
     6680            linearTerm[n]=linearTerm[j];
     6681            indices[n++]=j;
     6682          }
     6683        }
     6684        /// Constructor from constraint
     6685        constraints[numberConstraints++] = new ClpConstraintLinear(iRow, n, numberColumns,
     6686                      indices,linearTerm);
     6687        delete [] indices;
     6688      } else {
     6689        abort();
     6690      }
     6691    }
     6692  }
     6693  delete [] startQuadratic;
     6694  delete [] columnQuadratic;
     6695  delete [] elementQuadratic;
     6696  delete [] linearTerm;
     6697  delete [] which;
     6698  model->loadProblem(coinModel2);
     6699  int returnCode = model->nonlinearSLP(numberConstraints, constraints,
     6700                                       numberPasses,deltaTolerance);
     6701  for (iConstraint=0;iConstraint<saveNumber;iConstraint++)
     6702    delete constraints[iConstraint];
     6703 
     6704  delete [] constraints;
     6705  assert (!returnCode);
     6706  abort();
     6707  return model;
     6708}
  • branches/devel/Cbc/src/CbcLinked.hpp

    r588 r602  
    1313#endif
    1414#endif
     15#include "CoinModel.hpp"
    1516#ifdef COIN_HAS_LINK
    1617#include "OsiClpSolverInterface.hpp"
    17 #include "CoinModel.hpp"
    1818class CbcModel;
    1919class CoinPackedMatrix;
     
    11431143};
    11441144#endif
     1145class ClpSimplex;
     1146/** Return an approximate solution to a CoinModel.
     1147    Lots of bounds may be odd to force a solution.
     1148    mode = 0 just tries to get a continuous solution
     1149*/
     1150ClpSimplex * approximateSolution(CoinModel & coinModel,
     1151                                 int numberPasses, double deltaTolerance,
     1152                                 int mode=0);
    11451153#endif
  • branches/devel/Cbc/src/CbcSolver.cpp

    r601 r602  
    23772377              } else {
    23782378                // special solver
    2379                 double * solution = linkSolver->nonlinearSLP(slpValue,1.0e-5);
     2379                int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
     2380                double * solution = NULL;
     2381                if (testOsiOptions<10) {
     2382                  solution = linkSolver->nonlinearSLP(slpValue,1.0e-5);
     2383                } else {
     2384                  CoinModel coinModel = *linkSolver->coinModel();
     2385                  ClpSimplex * tempModel = approximateSolution(coinModel,slpValue>0 ? slpValue : 20 ,1.0e-5,0);
     2386                  assert (tempModel);
     2387                  solution = CoinCopyOfArray(tempModel->primalColumnSolution(),coinModel.numberColumns());
     2388                  delete tempModel;
     2389                }
    23802390                if (solution) {
    23812391                  memcpy(model2->primalColumnSolution(),solution,
     
    25912601                  if (obj) {
    25922602                    preProcess=0;
    2593                     parameters[whichParam(TESTOSI,numberParameters,parameters)].setIntValue(0);
     2603                    int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
     2604                    parameters[whichParam(TESTOSI,numberParameters,parameters)].setIntValue(CoinMax(0,testOsiOptions));
    25942605                    // create coin model
    25952606                    coinModel = lpSolver->createCoinModel();
  • branches/devel/Cbc/src/Cbc_ampl.cpp

    r595 r602  
    377377  for (i=0;i<saveArgc;i++) {
    378378    if (!strncmp(saveArgv[i],"testosi",7)) {
    379       testOsi = atoi(saveArgv[i+1]);
     379      testOsi = atoi(saveArgv[i]+8);
    380380      break;
    381381    }
    382382  }
    383   if (!(nlvc+nlvo)||testOsi>=10) {
     383  if (!(nlvc+nlvo)&&testOsi<10) {
    384384    /* read linear model*/
    385385    f_read(nl,0);
Note: See TracChangeset for help on using the changeset viewer.