Ignore:
Timestamp:
Feb 8, 2007 12:45:15 PM (13 years ago)
Author:
forrest
Message:

for nonlinear

File:
1 edited

Legend:

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

    r525 r530  
    3838#include "OsiRowCutDebugger.hpp"
    3939#include "OsiChooseVariable.hpp"
     40//#define CLP_MALLOC_STATISTICS
     41#ifdef CLP_MALLOC_STATISTICS
     42#include <exception>
     43#include <new>
     44static double malloc_times=0.0;
     45static double malloc_total=0.0;
     46static int malloc_amount[]={0,32,128,256,1024,4096,16384,65536,262144,INT_MAX};
     47static int malloc_n=10;
     48double malloc_counts[10]={0,0,0,0,0,0,0,0,0,0};
     49void * operator new (size_t size)
     50{
     51  malloc_times ++;
     52  malloc_total += size;
     53  int i;
     54  for (i=0;i<malloc_n;i++) {
     55    if ((int) size<=malloc_amount[i]) {
     56      malloc_counts[i]++;
     57      break;
     58    }
     59  }
     60  void * p =malloc(size);
     61  return p;
     62}
     63void operator delete (void *p)
     64{
     65  free(p);
     66}
     67static voif malloc_stats()
     68{
     69  double average = malloc_total/malloc_times;
     70  printf("count %g bytes %g - average %g\n",malloc_times,malloc_total,average);
     71  for (int i=0;i<malloc_n;i++)
     72    printf("%g ",malloc_counts[i]);
     73  printf("\n");
     74}
     75#endif
    4076#ifdef DMALLOC
    4177#include "dmalloc.h"
     
    5692#include "OsiRowCut.hpp"
    5793#include "OsiColCut.hpp"
     94#ifndef COIN_HAS_LINK
    5895#ifdef COIN_HAS_ASL
    5996#define COIN_HAS_LINK
    6097#endif
    61 #ifndef COIN_DEVELOP
    62 #undef COIN_HAS_LINK
    6398#endif
    6499#ifdef COIN_HAS_LINK
     
    426461  }
    427462}
     463/*  Returns OsiSolverInterface (User should delete)
     464    On entry numberKnapsack is maximum number of Total entries
     465*/
     466static OsiSolverInterface * 
     467expandKnapsack(CoinModel & model, int * whichColumn, int * knapsackStart,
     468               int * knapsackRow, int &numberKnapsack,
     469               CglStored & stored, int logLevel,
     470               int fixedPriority)
     471{
     472  int maxTotal = numberKnapsack;
     473  // load from coin model
     474  OsiSolverLink *si = new OsiSolverLink();
     475  OsiSolverInterface * finalModel=NULL;
     476  si->setDefaultMeshSize(0.001);
     477  // need some relative granularity
     478  si->setDefaultBound(100.0);
     479  si->setDefaultMeshSize(0.01);
     480  si->setDefaultBound(1000.0);
     481  si->setIntegerPriority(1000);
     482  si->setBiLinearPriority(10000);
     483  si->load(model,true,logLevel);
     484  // get priorities
     485  const int * priorities=model.priorities();
     486  int numberColumns = model.numberColumns();
     487  int SOSPriority=10000;
     488  if (priorities) {
     489    OsiObject ** objects = si->objects();
     490    int numberObjects = si->numberObjects();
     491    for (int iObj = 0;iObj<numberObjects;iObj++) {
     492      int iColumn = objects[iObj]->columnNumber();
     493      if (iColumn>=0&&iColumn<numberColumns) {
     494        OsiSimpleInteger * obj =
     495          dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
     496        assert (obj);
     497        int iPriority = priorities[iColumn];
     498        if (iPriority>0)
     499          objects[iObj]->setPriority(iPriority);
     500      }
     501    }
     502    if (fixedPriority>0) {
     503      si->setFixedPriority(fixedPriority);
     504      SOSPriority=fixedPriority+1;
     505    }
     506  }
     507  CoinModel coinModel=*si->coinModel();
     508  assert(coinModel.numberRows()>0);
     509  int numberRows = coinModel.numberRows();
     510  // Mark variables
     511  int * whichKnapsack = new int [numberColumns];
     512  int iRow,iColumn;
     513  for (iColumn=0;iColumn<numberColumns;iColumn++)
     514    whichKnapsack[iColumn]=-1;
     515  int kRow;
     516  bool badModel=false;
     517  // analyze
     518  if (logLevel>1) {
     519    for (iRow=0;iRow<numberRows;iRow++) {
     520      /* Just obvious one at first
     521         positive non unit coefficients
     522         all integer
     523         positive rowUpper
     524         for now - linear (but further down in code may use nonlinear)
     525         column bounds should be tight
     526      */
     527      //double lower = coinModel.getRowLower(iRow);
     528      double upper = coinModel.getRowUpper(iRow);
     529      if (upper<1.0e10) {
     530        CoinModelLink triple=coinModel.firstInRow(iRow);
     531        bool possible=true;
     532        int n=0;
     533        int n1=0;
     534        while (triple.column()>=0) {
     535          int iColumn = triple.column();
     536          const char *  el = coinModel.getElementAsString(iRow,iColumn);
     537          if (!strcmp("Numeric",el)) {
     538            if (coinModel.columnLower(iColumn)==coinModel.columnUpper(iColumn)) {
     539              triple=coinModel.next(triple);
     540              continue; // fixed
     541            }
     542            double value=coinModel.getElement(iRow,iColumn);
     543            if (value<0.0) {
     544              possible=false;
     545            } else {
     546              n++;
     547              if (value==1.0)
     548                n1++;
     549              if (coinModel.columnLower(iColumn)<0.0)
     550                possible=false;
     551              if (!coinModel.isInteger(iColumn))
     552                possible=false;
     553              if (whichKnapsack[iColumn]>=0)
     554                possible=false;
     555            }
     556          } else {
     557            possible=false; // non linear
     558          }
     559          triple=coinModel.next(triple);
     560        }
     561        if (n-n1>1&&possible) {
     562          double lower = coinModel.getRowLower(iRow);
     563          double upper = coinModel.getRowUpper(iRow);
     564          CoinModelLink triple=coinModel.firstInRow(iRow);
     565          while (triple.column()>=0) {
     566            int iColumn = triple.column();
     567            lower -= coinModel.columnLower(iColumn)*triple.value();
     568            upper -= coinModel.columnLower(iColumn)*triple.value();
     569            triple=coinModel.next(triple);
     570          }
     571          printf("%d is possible %g <=",iRow,lower);
     572          // print
     573          triple=coinModel.firstInRow(iRow);
     574          while (triple.column()>=0) {
     575            int iColumn = triple.column();
     576            if (coinModel.columnLower(iColumn)!=coinModel.columnUpper(iColumn))
     577              printf(" (%d,el %g up %g)",iColumn,triple.value(),
     578                     coinModel.columnUpper(iColumn)-coinModel.columnLower(iColumn));
     579            triple=coinModel.next(triple);
     580          }
     581          printf(" <= %g\n",upper);
     582        }
     583      }
     584    }
     585  }
     586  numberKnapsack=0;
     587  for (kRow=0;kRow<numberRows;kRow++) {
     588    iRow=kRow;
     589    /* Just obvious one at first
     590       positive non unit coefficients
     591       all integer
     592       positive rowUpper
     593       for now - linear (but further down in code may use nonlinear)
     594       column bounds should be tight
     595    */
     596    //double lower = coinModel.getRowLower(iRow);
     597    double upper = coinModel.getRowUpper(iRow);
     598    if (upper<1.0e10) {
     599      CoinModelLink triple=coinModel.firstInRow(iRow);
     600      bool possible=true;
     601      int n=0;
     602      int n1=0;
     603      while (triple.column()>=0) {
     604        int iColumn = triple.column();
     605        const char *  el = coinModel.getElementAsString(iRow,iColumn);
     606        if (!strcmp("Numeric",el)) {
     607          if (coinModel.columnLower(iColumn)==coinModel.columnUpper(iColumn)) {
     608            triple=coinModel.next(triple);
     609            continue; // fixed
     610          }
     611          double value=coinModel.getElement(iRow,iColumn);
     612          if (value<0.0) {
     613            possible=false;
     614          } else {
     615            n++;
     616            if (value==1.0)
     617              n1++;
     618            if (coinModel.columnLower(iColumn)<0.0)
     619              possible=false;
     620            if (!coinModel.isInteger(iColumn))
     621              possible=false;
     622            if (whichKnapsack[iColumn]>=0)
     623              possible=false;
     624          }
     625        } else {
     626          possible=false; // non linear
     627        }
     628        triple=coinModel.next(triple);
     629      }
     630      if (n-n1>1&&possible) {
     631        // try
     632        CoinModelLink triple=coinModel.firstInRow(iRow);
     633        while (triple.column()>=0) {
     634          int iColumn = triple.column();
     635          if (coinModel.columnLower(iColumn)!=coinModel.columnUpper(iColumn))
     636            whichKnapsack[iColumn]=numberKnapsack;
     637          triple=coinModel.next(triple);
     638        }
     639        knapsackRow[numberKnapsack++]=iRow;
     640      }
     641    }
     642  }
     643  if (logLevel>0)
     644    printf("%d out of %d candidate rows are possible\n",numberKnapsack,numberRows);
     645  // Check whether we can get rid of nonlinearities
     646  /* mark rows
     647     -2 in knapsack and other variables
     648     -1 not involved
     649     n only in knapsack n
     650  */
     651  int * markRow = new int [numberRows];
     652  for (iRow=0;iRow<numberRows;iRow++)
     653    markRow[iRow]=-1;
     654  int canDo=1; // OK and linear
     655  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     656    CoinModelLink triple=coinModel.firstInColumn(iColumn);
     657    int iKnapsack = whichKnapsack[iColumn];
     658    bool linear=true;
     659    // See if quadratic objective
     660    const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
     661    if (strcmp(expr,"Numeric")) {
     662      linear=false;
     663    }
     664    while (triple.row()>=0) {
     665      int iRow = triple.row();
     666      if (iKnapsack>=0) {
     667        if (markRow[iRow]==-1) {
     668          markRow[iRow]=iKnapsack;
     669        } else if (markRow[iRow]!=iKnapsack) {
     670          markRow[iRow]=-2;
     671        }
     672      }
     673      const char * expr = coinModel.getElementAsString(iRow,iColumn);
     674      if (strcmp(expr,"Numeric")) {
     675        linear=false;
     676      }
     677      triple=coinModel.next(triple);
     678    }
     679    if (!linear) {
     680      if (whichKnapsack[iColumn]<0) {
     681        canDo=0;
     682        break;
     683      } else {
     684        canDo=2;
     685      }
     686    }
     687  }
     688  int * markKnapsack = NULL;
     689  double * coefficient = NULL;
     690  double * linear = NULL;
     691  int * whichRow = NULL;
     692  int * lookupRow = NULL;
     693  badModel=(canDo==0);
     694  if (numberKnapsack&&canDo) {
     695    /* double check - OK if
     696       no nonlinear
     697       nonlinear only on columns in knapsack
     698       nonlinear only on columns in knapsack * ONE other - same for all in knapsack
     699       AND that is only row connected to knapsack
     700       (theoretically could split knapsack if two other and small numbers)
     701       also ONE could be ONE expression - not just a variable
     702    */
     703    int iKnapsack;
     704    markKnapsack = new int [numberKnapsack];
     705    coefficient = new double [numberKnapsack];
     706    linear = new double [numberColumns];
     707    for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++)
     708      markKnapsack[iKnapsack]=-1;
     709    if (canDo==2) {
     710      for (iRow=-1;iRow<numberRows;iRow++) {
     711        int numberOdd;
     712        CoinPackedMatrix * row = coinModel.quadraticRow(iRow,linear,numberOdd);
     713        if (row) {
     714          // see if valid
     715          const double * element = row->getElements();
     716          const int * column = row->getIndices();
     717          const CoinBigIndex * columnStart = row->getVectorStarts();
     718          const int * columnLength = row->getVectorLengths();
     719          int numberLook = row->getNumCols();
     720          for (int i=0;i<numberLook;i++) {
     721            int iKnapsack=whichKnapsack[i];
     722            if (iKnapsack<0) {
     723              // might be able to swap - but for now can't have knapsack in
     724              for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     725                int iColumn = column[j];
     726                if (whichKnapsack[iColumn]>=0) {
     727                  canDo=0; // no good
     728                  badModel=true;
     729                  break;
     730                }
     731              }
     732            } else {
     733              // OK if in same knapsack - or maybe just one
     734              int marked=markKnapsack[iKnapsack];
     735              for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     736                int iColumn = column[j];
     737                if (whichKnapsack[iColumn]!=iKnapsack&&whichKnapsack[iColumn]>=0) {
     738                  canDo=0; // no good
     739                  badModel=true;
     740                  break;
     741                } else if (marked==-1) {
     742                  markKnapsack[iKnapsack]=iColumn;
     743                  marked=iColumn;
     744                  coefficient[iKnapsack]=element[j];
     745                  coinModel.associateElement(coinModel.columnName(iColumn),1.0);
     746                } else if (marked!=iColumn) {
     747                  badModel=true;
     748                  canDo=0; // no good
     749                  break;
     750                } else {
     751                  // could manage with different coefficients - but for now ...
     752                  assert(coefficient[iKnapsack]==element[j]);
     753                }
     754              }
     755            }
     756          }
     757          delete row;
     758        }
     759      }
     760    }
     761    if (canDo) {
     762      // for any rows which are cuts
     763      whichRow = new int [numberRows];
     764      lookupRow = new int [numberRows];
     765      bool someNonlinear=false;
     766      double maxCoefficient=1.0;
     767      for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
     768        if (markKnapsack[iKnapsack]>=0) {
     769          someNonlinear=true;
     770          int iColumn = markKnapsack[iKnapsack];
     771          maxCoefficient = CoinMax(maxCoefficient,fabs(coefficient[iKnapsack]*coinModel.columnUpper(iColumn)));
     772        }
     773      }
     774      if (someNonlinear) {
     775        // associate all columns to stop possible error messages
     776        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     777          coinModel.associateElement(coinModel.columnName(iColumn),1.0);
     778        }
     779      }
     780      ClpSimplex tempModel;
     781      tempModel.loadProblem(coinModel);
     782      // Create final model - first without knapsacks
     783      int nCol=0;
     784      int nRow=0;
     785      for (iRow=0;iRow<numberRows;iRow++) {
     786        if (markRow[iRow]<0) {
     787          lookupRow[iRow]=nRow;
     788          whichRow[nRow++]=iRow;
     789        } else {
     790          lookupRow[iRow]=-1;
     791        }
     792      }
     793      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     794        if (whichKnapsack[iColumn]<0)
     795          whichColumn[nCol++]=iColumn;
     796      }
     797      ClpSimplex finalModelX(&tempModel,nRow,whichRow,nCol,whichColumn,false,false,false);
     798      OsiClpSolverInterface finalModelY(&finalModelX,true);
     799      finalModel = finalModelY.clone();
     800      finalModelY.releaseClp();
     801      // Put back priorities
     802      const int * priorities=model.priorities();
     803      if (priorities) {
     804        finalModel->findIntegers(false);
     805        OsiObject ** objects = finalModel->objects();
     806        int numberObjects = finalModel->numberObjects();
     807        for (int iObj = 0;iObj<numberObjects;iObj++) {
     808          int iColumn = objects[iObj]->columnNumber();
     809          if (iColumn>=0&&iColumn<nCol) {
     810            OsiSimpleInteger * obj =
     811              dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
     812            assert (obj);
     813            int iPriority = priorities[whichColumn[iColumn]];
     814            if (iPriority>0)
     815              objects[iObj]->setPriority(iPriority);
     816          }
     817        }
     818      }
     819      for (iRow=0;iRow<numberRows;iRow++) {
     820        whichRow[iRow]=iRow;
     821      }
     822      int numberOther=finalModel->getNumCols();
     823      int nLargest=0;
     824      int nelLargest=0;
     825      int nTotal=0;
     826      for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
     827        iRow = knapsackRow[iKnapsack];
     828        int nCreate = maxTotal;
     829        int nelCreate=coinModel.expandKnapsack(iRow,nCreate,NULL,NULL,NULL,NULL);
     830        if (nelCreate<0)
     831          badModel=true;
     832        nTotal+=nCreate;
     833        nLargest = CoinMax(nLargest,nCreate);
     834        nelLargest = CoinMax(nelLargest,nelCreate);
     835      }
     836      if (nTotal>maxTotal)
     837        badModel=true;
     838      if (!badModel) {
     839        // Now arrays for building
     840        nelLargest = CoinMax(nelLargest,nLargest)+1;
     841        double * buildObj = new double [nLargest];
     842        double * buildElement = new double [nelLargest];
     843        int * buildStart = new int[nLargest+1];
     844        int * buildRow = new int[nelLargest];
     845        OsiObject ** object = new OsiObject * [numberKnapsack];
     846        int nSOS=0;
     847        for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
     848          knapsackStart[iKnapsack]=finalModel->getNumCols();
     849          iRow = knapsackRow[iKnapsack];
     850          int nCreate = 10000;
     851          coinModel.expandKnapsack(iRow,nCreate,buildObj,buildStart,buildRow,buildElement);
     852          // Redo row numbers
     853          for (iColumn=0;iColumn<nCreate;iColumn++) {
     854            for (int j=buildStart[iColumn];j<buildStart[iColumn+1];j++) {
     855              int jRow=buildRow[j];
     856              jRow=lookupRow[jRow];
     857              assert (jRow>=0&&jRow<nRow);
     858              buildRow[j]=jRow;
     859            }
     860          }
     861          finalModel->addCols(nCreate,buildStart,buildRow,buildElement,NULL,NULL,buildObj);
     862          int numberFinal=finalModel->getNumCols();
     863          for (iColumn=numberOther;iColumn<numberFinal;iColumn++) {
     864            if (markKnapsack[iKnapsack]<0) {
     865              finalModel->setColUpper(iColumn,maxCoefficient);
     866              finalModel->setInteger(iColumn);
     867            } else {
     868              finalModel->setColUpper(iColumn,maxCoefficient+1.0);
     869              finalModel->setInteger(iColumn);
     870            }
     871            buildRow[iColumn-numberOther]=iColumn;
     872            buildElement[iColumn-numberOther]=1.0;
     873          }
     874          if (markKnapsack[iKnapsack]<0) {
     875            // convexity row
     876            finalModel->addRow(numberFinal-numberOther,buildRow,buildElement,1.0,1.0);
     877          } else {
     878            int iColumn = markKnapsack[iKnapsack];
     879            int n=numberFinal-numberOther;
     880            buildRow[n]=iColumn;
     881            buildElement[n++]=coefficient[iKnapsack];
     882            // convexity row (sort of)
     883            finalModel->addRow(n,buildRow,buildElement,0.0,0.0);
     884            OsiSOS * sosObject = new OsiSOS(finalModel,n-1,buildRow,NULL,1);
     885            sosObject->setPriority(iKnapsack+SOSPriority);
     886            // Say not integral even if is (switch off heuristics)
     887            sosObject->setIntegerValued(false);
     888            object[nSOS++]=sosObject;
     889          }
     890          numberOther=numberFinal;
     891        }
     892        finalModel->addObjects(nSOS,object);
     893        for (iKnapsack=0;iKnapsack<nSOS;iKnapsack++)
     894          delete object[iKnapsack];
     895        delete [] object;
     896        // Can we move any rows to cuts
     897        const int * cutMarker = coinModel.cutMarker();
     898        if (cutMarker) {
     899          // Row copy
     900          const CoinPackedMatrix * matrixByRow = finalModel->getMatrixByRow();
     901          const double * elementByRow = matrixByRow->getElements();
     902          const int * column = matrixByRow->getIndices();
     903          const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
     904          const int * rowLength = matrixByRow->getVectorLengths();
     905         
     906          const double * rowLower = finalModel->getRowLower();
     907          const double * rowUpper = finalModel->getRowUpper();
     908          int nDelete=0;
     909          for (iRow=0;iRow<numberRows;iRow++) {
     910            if (cutMarker[iRow]&&lookupRow[iRow]>=0) {
     911              int jRow=lookupRow[iRow];
     912              whichRow[nDelete++]=jRow;
     913              int start = rowStart[jRow];
     914              stored.addCut(rowLower[jRow],rowUpper[jRow],
     915                            rowLength[jRow],column+start,elementByRow+start);
     916            }
     917          }
     918          finalModel->deleteRows(nDelete,whichRow);
     919        }
     920        knapsackStart[numberKnapsack]=finalModel->getNumCols();
     921        delete [] buildObj;
     922        delete [] buildElement;
     923        delete [] buildStart;
     924        delete [] buildRow;
     925        finalModel->writeMps("full");
     926      }
     927    }
     928  }
     929  delete [] whichKnapsack;
     930  delete [] markRow;
     931  delete [] markKnapsack;
     932  delete [] coefficient;
     933  delete [] linear;
     934  delete [] whichRow;
     935  delete [] lookupRow;
     936  delete si;
     937  si=NULL;
     938  if (!badModel) {
     939    return finalModel;
     940  } else {
     941    delete finalModel;
     942    return NULL;
     943  }
     944}
     945// Fills in original solution (coinModel length)
     946static void
     947afterKnapsack(const CoinModel & coinModel2, const int * whichColumn, const int * knapsackStart,
     948                   const int * knapsackRow, int numberKnapsack,
     949                   const double * knapsackSolution, double * solution, int logLevel)
     950{
     951  CoinModel coinModel = coinModel2;
     952  int numberColumns = coinModel.numberColumns();
     953  int iColumn;
     954  // associate all columns to stop possible error messages
     955  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     956    coinModel.associateElement(coinModel.columnName(iColumn),1.0);
     957  }
     958  CoinZeroN(solution,numberColumns);
     959  int nCol=knapsackStart[0];
     960  for (iColumn=0;iColumn<nCol;iColumn++) {
     961    int jColumn = whichColumn[iColumn];
     962    solution[jColumn]=knapsackSolution[iColumn];
     963  }
     964  int * buildRow = new int [numberColumns]; // wild overkill
     965  double * buildElement = new double [numberColumns];
     966  int iKnapsack;
     967  for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
     968    int k=-1;
     969    double value=0.0;
     970    for (iColumn=knapsackStart[iKnapsack];iColumn<knapsackStart[iKnapsack+1];iColumn++) {
     971      if (knapsackSolution[iColumn]>1.0e-5) {
     972        if (k>=0) {
     973          printf("Two nonzero values for knapsack %d at (%d,%g) and (%d,%g)\n",iKnapsack,
     974                 k,knapsackSolution[k],iColumn,knapsackSolution[iColumn]);
     975          abort();
     976        }
     977        k=iColumn;
     978        value=floor(knapsackSolution[iColumn]+0.5);
     979        assert (fabs(value-knapsackSolution[iColumn])<1.0e-5);
     980      }
     981    }
     982    if (k>=0) {
     983      int iRow = knapsackRow[iKnapsack];
     984      int nCreate = 10000;
     985      int nel=coinModel.expandKnapsack(iRow,nCreate,NULL,NULL,buildRow,buildElement,k-knapsackStart[iKnapsack]);
     986      assert (nel);
     987      if (logLevel>0)
     988        printf("expanded column %d in knapsack %d has %d nonzero entries:\n",
     989               k-knapsackStart[iKnapsack],iKnapsack,nel);
     990      for (int i=0;i<nel;i++) {
     991        int jColumn = buildRow[i];
     992        double value = buildElement[i];
     993        if (logLevel>0)
     994          printf("%d - original %d has value %g\n",i,jColumn,value);
     995        solution[jColumn]=value;
     996      }
     997    }
     998  }
     999  delete [] buildRow;
     1000  delete [] buildElement;
     1001#if 0
     1002  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1003    if (solution[iColumn]>1.0e-5&&coinModel.isInteger(iColumn))
     1004      printf("%d %g\n",iColumn,solution[iColumn]);
     1005  }
     1006#endif
     1007}
    4281008static int outDupRow(OsiSolverInterface * solver)
    4291009{
     
    6061186    char * sosType = NULL;
    6071187    double * sosReference = NULL;
     1188    int * cut=NULL;
    6081189    int * sosPriority=NULL;
    6091190#ifdef COIN_HAS_ASL
    6101191    ampl_info info;
     1192    CglStored storedAmpl;
    6111193    CoinModel * coinModel = NULL;
     1194    CoinModel saveCoinModel;
     1195    int * whichColumn = NULL;
     1196    int * knapsackStart=NULL;
     1197    int * knapsackRow=NULL;
     1198    int numberKnapsack=0;
    6121199    memset(&info,0,sizeof(info));
    6131200    if (argc>2&&!strcmp(argv[2],"-AMPL")) {
     
    6531240                          info.columnLower,info.columnUpper,info.objective,
    6541241                          info.rowLower,info.rowUpper);
     1242      // take off cuts if ampl wants that
     1243      if (info.cut) {
     1244        int numberRows = info.numberRows;
     1245        int * whichRow = new int [numberRows];
     1246        // Row copy
     1247        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
     1248        const double * elementByRow = matrixByRow->getElements();
     1249        const int * column = matrixByRow->getIndices();
     1250        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
     1251        const int * rowLength = matrixByRow->getVectorLengths();
     1252       
     1253        const double * rowLower = solver->getRowLower();
     1254        const double * rowUpper = solver->getRowUpper();
     1255        int nDelete=0;
     1256        for (int iRow=0;iRow<numberRows;iRow++) {
     1257          if (info.cut[iRow]) {
     1258            whichRow[nDelete++]=iRow;
     1259            int start = rowStart[iRow];
     1260            storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
     1261                          rowLength[iRow],column+start,elementByRow+start);
     1262          }
     1263        }
     1264        solver->deleteRows(nDelete,whichRow);
     1265        delete [] whichRow;
     1266      }
    6551267#ifdef COIN_HAS_LINK
    6561268      } else {
     1269        // save
     1270        saveCoinModel = *coinModel;
    6571271        // load from coin model
    6581272        OsiSolverLink solver1;
     
    18202434              }
    18212435              if (!miplib) {
     2436                if (!preSolve) {
     2437                  model.solver()->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
     2438                  model.solver()->setHintParam(OsiDoPresolveInResolve,false,OsiHintTry);
     2439                }
    18222440                model.initialSolve();
    18232441                OsiSolverInterface * solver = model.solver();
     
    19702588                preProcess=0;
    19712589                useStrategy=true;
     2590                // empty out any cuts
     2591                if (storedAmpl.sizeRowCuts()) {
     2592                  printf("Emptying ampl stored cuts as internal preprocessing\n");
     2593                  CglStored temp;
     2594                  storedAmpl=temp;
     2595                }
    19722596              }
    19732597              if (preProcess&&type==BAB) {
     
    21272751              }
    21282752              int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
     2753#ifdef COIN_HAS_LINK
     2754              // If linked then see if expansion wanted
     2755              {
     2756                OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
     2757                if (solver3) {
     2758                  int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000;
     2759                  if (options) {
     2760                    /*
     2761                      1 - force mini branch and bound
     2762                      2 - set priorities high on continuous
     2763                      4 - try adding OA cuts
     2764                      8 - try doing quadratic linearization
     2765                      16 - try expanding knapsacks
     2766                    */
     2767                    if ((options&16)) {
     2768                      int numberColumns = saveCoinModel.numberColumns();
     2769                      int numberRows = saveCoinModel.numberRows();
     2770                      whichColumn = new int[numberColumns];
     2771                      knapsackStart=new int[numberRows+1];
     2772                      knapsackRow=new int[numberRows];
     2773                      numberKnapsack=10000;
     2774                      OsiSolverInterface * solver = expandKnapsack(saveCoinModel,whichColumn,knapsackStart,
     2775                                                                   knapsackRow,numberKnapsack,
     2776                                                                   storedAmpl,2,slpValue);
     2777                      if (solver) {
     2778                        clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
     2779                        assert (clpSolver);
     2780                        lpSolver = clpSolver->getModelPtr();
     2781                        babModel->assignSolver(solver);
     2782                        testOsiOptions=0;
     2783                      } else {
     2784                        numberKnapsack=0;
     2785                        delete [] whichColumn;
     2786                        delete [] knapsackStart;
     2787                        delete [] knapsackRow;
     2788                        whichColumn=NULL;
     2789                        knapsackStart=NULL;
     2790                        knapsackRow=NULL;
     2791                      }
     2792                    }
     2793                  }
     2794                }
     2795              }
     2796#endif
    21292797              if (useCosts&&testOsiOptions<0) {
    21302798                int numberColumns = babModel->getNumCols();
     
    28353503                  free(sosReference);
    28363504                  sosReference=NULL;
     3505                  free(cut);
     3506                  cut=NULL;
    28373507                  free(sosPriority);
    28383508                  sosPriority=NULL;
     
    29183588                          4 - try adding OA cuts
    29193589                          8 - try doing quadratic linearization
     3590                          16 - try expanding knapsacks
    29203591                        */
    29213592                        if ((options&2))
     
    30113682                  osiclp->setSolveOptions(clpSolve);
    30123683                  osiclp->setHintParam(OsiDoDualInResolve,false);
     3684                  // switch off row copy
     3685                  osiclp->getModelPtr()->setSpecialOptions(osiclp->getModelPtr()->specialOptions()|256);
     3686                }
     3687                if (storedAmpl.sizeRowCuts()) {
     3688                  if (preProcess) {
     3689                    const int * originalColumns = process.originalColumns();
     3690                    int numberColumns = babModel->getNumCols();
     3691                    int * newColumn = new int[numberOriginalColumns];
     3692                    int i;
     3693                    for (i=0;i<numberOriginalColumns;i++)
     3694                      newColumn[i]=-1;
     3695                    for (i=0;i<numberColumns;i++) {
     3696                      int iColumn = originalColumns[i];
     3697                      newColumn[iColumn]=i;
     3698                    }
     3699                    int * buildColumn = new int[numberColumns];
     3700                    // Build up valid cuts
     3701                    int nBad=0;
     3702                    int nCuts = storedAmpl.sizeRowCuts();
     3703                    CglStored newCuts;
     3704                    for (i=0;i<nCuts;i++) {
     3705                      const OsiRowCut * cut = storedAmpl.rowCutPointer(i);
     3706                      double lb = cut->lb();
     3707                      double ub = cut->ub();
     3708                      int n=cut->row().getNumElements();
     3709                      const int * column = cut->row().getIndices();
     3710                      const double * element = cut->row().getElements();
     3711                      bool bad=false;
     3712                      for (int i=0;i<n;i++) {
     3713                        int iColumn = column[i];
     3714                        iColumn = newColumn[iColumn];
     3715                        if (iColumn>=0) {
     3716                          buildColumn[i]=iColumn;
     3717                        } else {
     3718                          bad=true;
     3719                          break;
     3720                        }
     3721                      }
     3722                      if (!bad) {
     3723                        newCuts.addCut(lb,ub,n,buildColumn,element);
     3724                      } else {
     3725                        nBad++;
     3726                      }
     3727                    }
     3728                    storedAmpl=newCuts;
     3729                    if (nBad)
     3730                      printf("%d cuts dropped\n",nBad);
     3731                    delete [] newColumn;
     3732                    delete [] buildColumn;
     3733                  }
     3734                  if (storedAmpl.sizeRowCuts())
     3735                    babModel->addCutGenerator(&storedAmpl,1,"AmplStored");
    30133736                }
    30143737                babModel->branchAndBound(statistics);
     
    30763799              if (type==STRENGTHEN&&strengthenedModel)
    30773800                clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel);
     3801              else if (usingAmpl)
     3802                clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
    30783803              lpSolver = clpSolver->getModelPtr();
    30793804              if (numberChanged) {
     
    31063831#ifdef COIN_HAS_ASL
    31073832                if (usingAmpl) {
     3833                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
     3834                  lpSolver = clpSolver->getModelPtr();
    31083835                  double value = babModel->getObjValue()*lpSolver->getObjSense();
    31093836                  char buf[300];
     
    31483875                  if (bestSolution) {
    31493876                    free(info.primalSolution);
    3150                     info.primalSolution = (double *) malloc(n*sizeof(double));
    3151                     CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
    3152                     int numberRows = lpSolver->numberRows();
    3153                     free(info.dualSolution);
    3154                     info.dualSolution = (double *) malloc(numberRows*sizeof(double));
    3155                     CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
     3877                    if (!numberKnapsack) {
     3878                      info.primalSolution = (double *) malloc(n*sizeof(double));
     3879                      CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
     3880                      int numberRows = lpSolver->numberRows();
     3881                      free(info.dualSolution);
     3882                      info.dualSolution = (double *) malloc(numberRows*sizeof(double));
     3883                      CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
     3884                    } else {
     3885                      // expanded knapsack
     3886                      info.dualSolution=NULL;
     3887                      int numberColumns = saveCoinModel.numberColumns();
     3888                      info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
     3889                      // Fills in original solution (coinModel length)
     3890                      afterKnapsack(saveCoinModel,  whichColumn,  knapsackStart,
     3891                                    knapsackRow,  numberKnapsack,
     3892                                    lpSolver->primalColumnSolution(), info.primalSolution,1);
     3893                    }
    31563894                  } else {
    31573895                    info.primalSolution=NULL;
     
    31983936                free(sosReference);
    31993937                sosReference=NULL;
     3938                free(cut);
     3939                cut=NULL;
    32003940                free(sosPriority);
    32013941                sosPriority=NULL;
Note: See TracChangeset for help on using the changeset viewer.