source: branches/devel/Cbc/examples/CbcSolverLink.cpp @ 466

Last change on this file since 466 was 466, checked in by forrest, 13 years ago

for nonlinear

File size: 36.1 KB
Line 
1// Copyright (C) 2006, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <cassert>
5
6#include "CoinTime.hpp"
7
8#include "CoinHelperFunctions.hpp"
9#include "CoinIndexedVector.hpp"
10#include "CoinMpsIO.hpp"
11#include "CoinModel.hpp"
12#include "ClpSimplex.hpp"
13#include "CbcSolverLink.hpp"
14#include "CbcModel.hpp"
15#include "ClpPackedMatrix.hpp"
16#include "CbcCutGenerator.hpp"
17#include "CbcStrategy.hpp"
18#include "CglPreProcess.hpp"
19#include "CoinTime.hpp"
20#include "CglProbing.hpp"
21#include "CglProbing.hpp"
22#include "CglGomory.hpp"
23#include "CglKnapsackCover.hpp"
24#include "CglRedSplit.hpp"
25#include "CglClique.hpp"
26#include "CglMixedIntegerRounding2.hpp"
27#include "CglFlowCover.hpp"
28#include "CglTwomir.hpp"
29#include "CbcHeuristicFPump.hpp"
30#include "CbcHeuristic.hpp"
31#include "CbcHeuristicLocal.hpp"
32#include "CbcHeuristicGreedy.hpp"
33#include "CbcHeuristicGreedy.hpp"
34#include "CbcCompareActual.hpp"
35#include "CbcBranchLink.hpp"
36//#############################################################################
37// Solve methods
38//#############################################################################
39void CbcSolverLink::initialSolve()
40{
41  specialOptions_ =0;
42  modelPtr_->setWhatsChanged(0);
43  ClpMatrixBase * save = NULL;
44  if (numberVariables_) {
45    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
46    // update all bounds before coefficients
47    for (int i=0;i<numberVariables_;i++ ) {
48      info_[i].updateBounds(modelPtr_);
49    }
50    for (int i=0;i<numberVariables_;i++ ) {
51      info_[i].updateCoefficients(modelPtr_,temp);
52    }
53    temp->removeGaps(1.0e-14);
54    save = modelPtr_->clpMatrix();
55    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
56    assert (clpMatrix);
57    if (save->getNumRows()>temp->getNumRows()) {
58      // add in cuts
59      int numberRows = temp->getNumRows();
60      int * which = new int[numberRows];
61      for (int i=0;i<numberRows;i++)
62        which[i]=i;
63      save->deleteRows(numberRows,which);
64      delete [] which;
65      temp->bottomAppendPackedMatrix(*clpMatrix->matrix());
66    }
67    modelPtr_->replaceMatrix(temp);
68  }
69  OsiClpSolverInterface::initialSolve();
70  if (save) {
71    delete save;
72  }
73}
74
75//-----------------------------------------------------------------------------
76void CbcSolverLink::resolve()
77{
78  specialOptions_ =0;
79  modelPtr_->setWhatsChanged(0);
80  ClpMatrixBase * save = NULL;
81  bool allFixed=false;
82  bool feasible=true;
83  if (numberVariables_) {
84    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
85    allFixed=true;
86    //bool best=true;
87    const double * lower = modelPtr_->columnLower();
88    const double * upper = modelPtr_->columnUpper();
89    // update all bounds before coefficients
90    for (int i=0;i<numberVariables_;i++ ) {
91      info_[i].updateBounds(modelPtr_);
92      int iColumn = info_[i].variable();
93      double lo = lower[iColumn];
94      double up = upper[iColumn];
95      if (up>lo)
96        allFixed=false;
97      else if (up<lo)
98        feasible=false;
99    }
100    for (int i=0;i<numberVariables_;i++ ) {
101      info_[i].updateCoefficients(modelPtr_,temp);
102    }
103    temp->removeGaps(1.0e-14);
104    save = modelPtr_->clpMatrix();
105    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
106    assert (clpMatrix);
107    if (save->getNumRows()>temp->getNumRows()) {
108      // add in cuts
109      int numberRows = temp->getNumRows();
110      int * which = new int[numberRows];
111      for (int i=0;i<numberRows;i++)
112        which[i]=i;
113      CoinPackedMatrix * mat = clpMatrix->matrix();
114      // for debug
115      //mat = new CoinPackedMatrix(*mat);
116      mat->deleteRows(numberRows,which);
117      delete [] which;
118      temp->bottomAppendPackedMatrix(*mat);
119      temp->removeGaps(1.0e-14);
120    }
121    if (0) {
122      const CoinPackedMatrix * matrix = modelPtr_->matrix();
123      int numberColumns = matrix->getNumCols();
124      assert (numberColumns==temp->getNumCols());
125      const double * element1 = temp->getMutableElements();
126      const int * row1 = temp->getIndices();
127      const CoinBigIndex * columnStart1 = temp->getVectorStarts();
128      const int * columnLength1 = temp->getVectorLengths();
129      const double * element2 = matrix->getMutableElements();
130      const int * row2 = matrix->getIndices();
131      const CoinBigIndex * columnStart2 = matrix->getVectorStarts();
132      const int * columnLength2 = matrix->getVectorLengths();
133      for (int i=0;i<numberColumns;i++) {
134        assert (columnLength2[i]==columnLength1[i]);
135        int offset = columnStart2[i]-columnStart1[i];
136        for (int j=columnStart1[i];j<columnStart1[i]+columnLength1[i];j++) {
137          assert (row1[j]==row2[j+offset]);
138          assert (element1[j]==element2[j+offset]);
139        }
140      }
141    }
142    modelPtr_->replaceMatrix(temp);
143  }
144  if (!feasible)
145    allFixed=false;
146  if ((specialOptions2_&1)!=0)
147    allFixed=false;
148  int returnCode=-1;
149  if (feasible) {
150    // may do lots of work
151    returnCode=fathom(allFixed);
152  }
153  if (returnCode>=0) {
154    if (returnCode==0)
155      OsiClpSolverInterface::resolve();
156    if (!allFixed&&(specialOptions2_&1)==0) {
157      const double * solution = getColSolution();
158      bool satisfied=true;
159      for (int i=0;i<numberVariables_;i++) {
160        int iColumn = info_[i].variable();
161        double value = solution[iColumn];
162        if (fabs(value-floor(value+0.5))>0.0001)
163          satisfied=false;
164      }
165      //if (satisfied)
166      //printf("satisfied but not fixed\n");
167    }
168  } else {
169    modelPtr_->setProblemStatus(1);
170    modelPtr_->setObjectiveValue(COIN_DBL_MAX);
171  }
172  if (save) {
173    delete save;
174  }
175}
176
177//#############################################################################
178// Constructors, destructors clone and assignment
179//#############################################################################
180
181//-------------------------------------------------------------------
182// Default Constructor
183//-------------------------------------------------------------------
184CbcSolverLink::CbcSolverLink ()
185  : OsiClpSolverInterface()
186{
187  gutsOfDestructor(true);
188}
189/* returns
190   sequence of nonlinear or
191   -1 numeric
192   -2 not found
193   -3 too many terms
194*/
195static int getVariable(const CoinModel & model, char * expression,
196                       int & linear)
197{
198  int non=-1;
199  linear=-1;
200  if (strcmp(expression,"Numeric")) {
201    // function
202    char * first = strchr(expression,'*');
203    int numberColumns = model.numberColumns();
204    int j;
205    if (first) {
206      *first='\0';
207      for (j=0;j<numberColumns;j++) {
208        if (!strcmp(expression,model.columnName(j))) {
209          linear=j;
210          memmove(expression,first+1,strlen(first+1)+1);
211          break;
212        }
213      }
214    }
215    // find nonlinear
216    for (j=0;j<numberColumns;j++) {
217      const char * name = model.columnName(j);
218      first = strstr(expression,name);
219      if (first) {
220        if (first!=expression&&isalnum(*(first-1)))
221          continue; // not real match
222        first += strlen(name);
223        if (!isalnum(*first)) {
224          // match
225          non=j;
226          // but check no others
227          j++;
228          for (;j<numberColumns;j++) {
229            const char * name = model.columnName(j);
230            first = strstr(expression,name);
231            if (first) {
232              if (isalnum(*(first-1)))
233                continue; // not real match
234              first += strlen(name);
235              if (!isalnum(*first)) {
236                // match - ouch
237                non=-3;
238                break;
239              }
240            }
241          }
242          break;
243        }
244      }
245    }
246    if (non==-1)
247      non=-2;
248  } 
249  return non;
250}
251/* This creates from a coinModel object
252
253   if errors.then number of sets is -1
254     
255   This creates linked ordered sets information.  It assumes -
256
257   for product terms syntax is yy*f(zz)
258   also just f(zz) is allowed
259   and even a constant
260
261   modelObject not const as may be changed as part of process.
262*/
263CbcSolverLink::CbcSolverLink ( CoinModel & coinModel)
264  : OsiClpSolverInterface()
265{
266  gutsOfDestructor(true);
267  gdb(coinModel);
268}
269void CbcSolverLink::gdb ( CoinModel & coinModel)
270{
271  // first check and set up arrays
272  int numberColumns = coinModel.numberColumns();
273  int numberRows = coinModel.numberRows();
274  // List of nonlinear entries
275  int * which = new int[numberColumns];
276  numberSets_=0;
277  numberNonlinear_=0;
278  int iColumn;
279  int numberErrors=0;
280  for (iColumn=0;iColumn<numberColumns;iColumn++) {
281    CoinModelLink triple=coinModel.firstInColumn(iColumn);
282    bool linear=true;
283    int n=0;
284    while (triple.row()>=0) {
285      int iRow = triple.row();
286      const char * expr = coinModel.getElementAsString(iRow,iColumn);
287      if (strcmp(expr,"Numeric")) {
288        linear=false;
289        // associate to avoid errors - much better to use current solution values
290        coinModel.associateElement(expr,0.0);
291      }
292      triple=coinModel.next(triple);
293      n++;
294    }
295    if (!linear) {
296      which[numberSets_++]=iColumn;
297      numberNonlinear_ += n;
298    }
299  }
300  // return if nothing
301  simplexModel_ = new ClpSimplex(*modelPtr_);
302  simplexModel_->loadProblem(coinModel);
303  if (!numberSets_) {
304    delete [] which;
305    return;
306  } else {
307    simplexModel_->deleteColumns(numberSets_,which);
308    coinModel_ = coinModel;
309    which_ = CoinCopyOfArray(which,numberSets_);
310    delete [] which;
311    if (!modelPtr_->numberColumns()) {
312      delete modelPtr_;
313      modelPtr_ = new ClpSimplex(*simplexModel_);
314      // synchronize integers
315      int numberColumns=simplexModel_->numberColumns();
316      for (int i=0;i<numberColumns;i++) {
317        if (simplexModel_->isInteger(i))
318          this->setInteger(i);
319      }
320    }
321  }
322  // set starts
323  startSet_ = new int[numberSets_+1];
324  definitionRow_ = new int [numberNonlinear_];
325  definitionRowL_ = new int [numberNonlinear_];
326  functionColumn_ = new int [numberSets_];
327  definitionColumnL_ = new int [numberNonlinear_];
328  usageRow_ = new int [numberNonlinear_];
329  convexityRow_ = new int[numberNonlinear_];
330
331  startSet_[0]=0;
332  int nonLinear=0;
333  int nRow=numberRows;
334  int iSet;
335  char * mentioned=new char[numberColumns];
336  memset(mentioned,0,numberColumns);
337  for (iSet=0;iSet<numberSets_;iSet++) {
338    iColumn = which_[iSet];
339    // first get constants
340    CoinModelLink triple=coinModel.firstInColumn(iColumn);
341    int functionColumn=-1;
342    while (triple.row()>=0) {
343      int iRow = triple.row();
344      const char * expr = coinModel.getElementAsString(iRow,iColumn);
345      int linear;
346      assert (strlen(expr)<1000);
347      char expr2 [1000];
348      strcpy(expr2,expr);
349      int non = getVariable(coinModel,expr2,linear);
350      if (non==-1) {
351        // add constant
352        definitionRowL_[nonLinear]=-2;
353        definitionRow_[nonLinear]=-2;
354        definitionColumnL_[nonLinear]=-2;
355        usageRow_[nonLinear]=iRow;
356      } else if (non<0) {
357        // error
358        printf("bad expression %s for column %s and row %s\n",
359               expr,coinModel.getColumnName(iColumn),coinModel.getRowName(iRow));
360        numberErrors++;
361      } else {
362        if (functionColumn==-1) {
363          functionColumn=non;
364        } else if (functionColumn!=non) {
365          printf("second function of %s when already got %s - expression %s\n",
366                 coinModel.getColumnName(non),coinModel.getColumnName(functionColumn),
367                 expr);
368          numberErrors++;
369        }
370      }
371      triple=coinModel.next(triple);
372    }
373    if (!numberErrors) {
374      triple=coinModel.firstInColumn(iColumn);
375      functionColumn_[iSet]=functionColumn;
376      while (triple.row()>=0) {
377        int iRow = triple.row();
378        const char * expr = coinModel.getElementAsString(iRow,iColumn);
379        int linear;
380        assert (strlen(expr)<1000);
381        char expr2 [1000];
382        strcpy(expr2,expr);
383        int non = getVariable(coinModel,expr2,linear);
384        if (non>=0) {
385          mentioned[non]=1;
386          if(linear>=0)
387            mentioned[linear]=1;
388          // definition row
389          double elementValue;
390          definitionRow_[nonLinear]=nRow;
391          elementValue = -1.0;
392          simplexModel_->addRow(1,&functionColumn,&elementValue,0.0,0.0);
393          nRow++;
394          // convexity row
395          convexityRow_[nonLinear]=nRow;
396          simplexModel_->addRow(0,NULL,NULL,1.0,1.0);
397          nRow++;
398          // other definition row
399          if (linear>=0) {
400            definitionRowL_[nonLinear]=nRow;
401            definitionColumnL_[nonLinear]=linear;
402            elementValue = -1.0;
403            simplexModel_->addRow(1,&linear,&elementValue,0.0,0.0);
404            nRow++;
405          } else {
406            definitionRowL_[nonLinear]=-1;
407            definitionColumnL_[nonLinear]=-1;
408          }
409          usageRow_[nonLinear]=iRow;
410          coinModel_.setElement(iRow,functionColumn,expr2);
411          nonLinear++;
412        }
413        triple=coinModel.next(triple);
414      }
415    }
416    startSet_[iSet+1]=nonLinear;
417  }
418  if (numberErrors) {
419    // errors
420    gutsOfDestructor();
421    numberSets_=-1;
422  } else {
423    int numberRows=coinModel_.numberRows();
424    for (int i=0;i<numberColumns;i++) {
425      if (mentioned[i])
426        coinModel_.setElement(numberRows,i,coinModel_.getColumnName(i));
427    }
428  }
429  delete [] mentioned;
430}
431
432//-------------------------------------------------------------------
433// Clone
434//-------------------------------------------------------------------
435OsiSolverInterface * 
436CbcSolverLink::clone(bool copyData) const
437{
438  assert (copyData);
439  return new CbcSolverLink(*this);
440}
441
442
443//-------------------------------------------------------------------
444// Copy constructor
445//-------------------------------------------------------------------
446CbcSolverLink::CbcSolverLink (
447                  const CbcSolverLink & rhs)
448  : OsiClpSolverInterface(rhs)
449{
450  gutsOfDestructor(true);
451  gutsOfCopy(rhs);
452}
453
454//-------------------------------------------------------------------
455// Destructor
456//-------------------------------------------------------------------
457CbcSolverLink::~CbcSolverLink ()
458{
459  gutsOfDestructor();
460}
461
462//-------------------------------------------------------------------
463// Assignment operator
464//-------------------------------------------------------------------
465CbcSolverLink &
466CbcSolverLink::operator=(const CbcSolverLink& rhs)
467{
468  if (this != &rhs) { 
469    gutsOfDestructor();
470    OsiClpSolverInterface::operator=(rhs);
471    gutsOfCopy(rhs);
472  }
473  return *this;
474}
475void 
476CbcSolverLink::gutsOfDestructor(bool justNullify)
477{
478  if (!justNullify) {
479    delete matrix_;
480    delete [] info_;
481    delete [] bestSolution_;
482    delete [] startSet_;
483    delete [] definitionRow_;
484    delete [] definitionRowL_;
485    delete [] functionColumn_;
486    delete [] definitionColumnL_;
487    delete [] usageRow_;
488    delete [] convexityRow_;
489    delete [] granularityTop_;
490    delete [] granularityBottom_;
491    delete [] which_;
492   
493    delete simplexModel_;
494    assert (savedModel_==NULL);
495  } 
496  matrix_ = NULL;
497  info_ = NULL;
498  bestSolution_ = NULL;
499  startSet_ = NULL;
500  definitionRow_ = NULL;
501  definitionRowL_ = NULL;
502  definitionColumnL_ = NULL;
503  functionColumn_ = NULL; 
504  usageRow_ = NULL;
505  convexityRow_ = NULL;
506  simplexModel_ = NULL;
507  savedModel_=NULL;
508  granularityTop_ = NULL;
509  granularityBottom_ = NULL; 
510  which_ = NULL;
511
512  model_=NULL;
513  numberVariables_ = 0;
514  bestObjectiveValue_ =1.0e100;
515  specialOptions2_ = 0;
516  numberSets_ = 0;
517  numberNonlinear_ = 0;
518}
519void 
520CbcSolverLink::gutsOfCopy(const CbcSolverLink & rhs)
521{
522  model_ = rhs.model_;
523  coinModel_ = rhs.coinModel_;
524  numberVariables_ = rhs.numberVariables_;
525  bestObjectiveValue_ = rhs.bestObjectiveValue_;
526  specialOptions2_ = rhs.specialOptions2_;
527  if (numberVariables_) { 
528    matrix_ = new CoinPackedMatrix(*rhs.matrix_);
529    info_ = new CbcLinkedBound [numberVariables_];
530    for (int i=0;i<numberVariables_;i++) {
531      info_[i] = CbcLinkedBound(rhs.info_[i]);
532    }
533    if (rhs.bestSolution_) {
534      bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->getNumCols());
535    } else {
536      bestSolution_=NULL;
537    }
538  }
539  numberSets_ = rhs.numberSets_;
540  numberNonlinear_ = rhs.numberNonlinear_;
541  if (numberSets_>0) {
542    startSet_ = CoinCopyOfArray(rhs.startSet_,numberSets_+1);
543    definitionRow_ = CoinCopyOfArray(rhs.definitionRow_,numberNonlinear_);
544    definitionRowL_ = CoinCopyOfArray(rhs.definitionRowL_,numberNonlinear_);
545    functionColumn_ = CoinCopyOfArray(rhs.functionColumn_,numberSets_);
546    definitionColumnL_ = CoinCopyOfArray(rhs.definitionColumnL_,numberNonlinear_);
547    convexityRow_ = CoinCopyOfArray(rhs.convexityRow_,numberNonlinear_);
548    usageRow_ = CoinCopyOfArray(rhs.usageRow_,numberNonlinear_);
549    simplexModel_ = new ClpSimplex(*rhs.simplexModel_);
550    assert (rhs.savedModel_==NULL);
551    granularityTop_ = CoinCopyOfArray(rhs.granularityTop_,numberSets_);
552    granularityBottom_ = CoinCopyOfArray(rhs.granularityBottom_,numberSets_);
553    which_ = CoinCopyOfArray(rhs.which_,numberSets_);
554  }
555}
556//-------------------------------------------------------------------
557// Real initializer
558//-------------------------------------------------------------------
559void
560CbcSolverLink::initialize(CbcModel * model, const CoinPackedMatrix * matrix,
561                          int numberVariables, const int * which,
562                          const int *starts, const int * positionL, const int * positionU)
563{
564  model_=model;
565  numberVariables_ = numberVariables;
566  matrix_ = new CoinPackedMatrix(*matrix);
567  info_ = new CbcLinkedBound [numberVariables_];
568  // get multipliers from matrix
569  int n = starts[numberVariables_];
570  double * multipliers = new double[n];
571  int i;
572  const double * elements = matrix_->getElements();
573  for ( i=0;i<n;i++) {
574    int position = positionL[i];
575    double value = elements[position];
576    multipliers[i]=value;
577    position = positionU[i];
578    value = elements[position];
579    assert(multipliers[i]==value);
580  }
581  for ( i=0;i<numberVariables_;i++) {
582    info_[i] = CbcLinkedBound(model,which[i],starts[i+1]-starts[i],
583                              positionL+starts[i],
584                              positionU+starts[i],
585                              multipliers+starts[i]);
586  }
587  delete [] multipliers;
588}
589void
590CbcSolverLink::initialize ( CbcModel * model, bool top)
591{
592  model_=model;
593  int numberExtraColumns=0;
594  int numberExtraElements=0;
595  int iSet;
596  int iColumn;
597  double * granularity = CoinCopyOfArray(top ? granularityTop_ : granularityBottom_, numberSets_);
598  if (!granularity) {
599    printf("no granularity given\n");
600    abort();
601  }
602  int * numberInSet = new int[numberSets_];
603  int * setType = new int[numberSets_];
604  assert (simplexModel_);
605  // Use Current bounds
606  const double * lower = modelPtr_->columnLower();
607  const double * upper = modelPtr_->columnUpper();
608  // To point to interesting variables
609  int numberColumns = simplexModel_->numberColumns();
610  int * back = new int[numberColumns];
611  CoinZeroN(back,numberColumns);
612  int largestSet=0;
613  for (iSet=0;iSet<numberSets_;iSet++) {
614    int zColumn = functionColumn_[iSet];
615    double lo = lower[zColumn];
616    double up = upper[zColumn];
617    assert (up-lo<1.0e10);
618    double g = granularity[iSet];
619    assert (g);
620    int number;
621    if (g>0) {
622      number = (int) floor((up-lo)/g+0.5);
623      if (number<3)
624        number = (int) ceil((up-lo)/g);
625      setType[iSet]=2;
626    } else {
627      g = -g;
628      // must be exact
629      number = (int) floor((up-lo)/g+1.0e-6*g);
630      if (fabs(up-lo+number*g)>1.0e-6) {
631        printf("granularity of %g for set %d does not match bounds %g and %g\n",
632               -g,iSet,lo,up);
633        abort();
634      }
635      setType[iSet]=1;
636    }
637    granularity[iSet] = (up-lo) / ((double) number);
638    // and end point
639    number++;
640    numberInSet[iSet]=number;
641    int nElements=0;
642    int nLink=0;
643    int i;
644    // allow zeros to stay
645    for (i=startSet_[iSet];i<startSet_[iSet+1];i++) {
646      if (definitionRow_[i]>=0) {
647        nElements += 2*number;
648        nLink++;
649      }
650      if (definitionRowL_[i]>=0)
651        nElements += 2*number;
652      // output
653      nElements += 2*number;
654      // convexity
655      nElements += 2*number;
656      int iColumn = definitionColumnL_[i];
657      if (iColumn>=0)
658        back[iColumn]=1;
659    }
660    int functionColumn = functionColumn_[iSet];
661    back[functionColumn]=1;
662    number *= 2*nLink;
663    numberExtraColumns+=number;
664    largestSet = CoinMax(largestSet,number);
665    numberExtraElements+=number*nElements;
666  }
667  CbcObject ** objects = new CbcObject * [numberSets_];
668  double * newLower = new double[numberExtraColumns];
669  double * newUpper = new double[numberExtraColumns];
670  double * newObjective = new double[numberExtraColumns];
671  for (iColumn=0;iColumn<numberExtraColumns;iColumn++) {
672    newLower[iColumn]=0.0;
673    newUpper[iColumn]=COIN_DBL_MAX;
674    newObjective[iColumn]=0.0;
675  }
676  int * start = new int[numberExtraColumns+1];
677  int * row = new int[numberExtraElements];
678  double * element = new double[numberExtraElements];
679  // For objects
680  double * weight = new double[largestSet];
681  int * sequence = new int [largestSet];
682  // get copy
683  CbcSolverLink * newSolver = new CbcSolverLink(*this);
684  // replace model
685  newSolver->freeCachedResults();
686  delete newSolver->modelPtr_;
687  newSolver->modelPtr_=new ClpSimplex(*simplexModel_);
688  ClpSimplex * newClp = newSolver->modelPtr_;
689  // correct bounds
690  memcpy(newClp->columnLower(),lower,numberColumns*sizeof(double));
691  memcpy(newClp->columnUpper(),upper,numberColumns*sizeof(double));
692  // space for bound information
693  numberVariables_ = 0;
694  delete [] info_;
695  for (iColumn=0;iColumn<numberColumns;iColumn++) {
696    if (back[iColumn]) {
697      back[iColumn]=numberVariables_;
698      numberVariables_++;
699    } else {
700      back[iColumn]=-1;
701    }
702  }
703  info_ = new CbcLinkedBound [numberVariables_];
704  for (iColumn=0;iColumn<numberColumns;iColumn++) {
705    int iVar = back[iColumn];
706    if (iVar>=0) 
707      info_[iVar] = CbcLinkedBound(NULL,iColumn,0,NULL,NULL,NULL);
708  }
709  // First run just puts in junk and does objects
710  start[0]=0;
711  int nel=0;
712  int n=0;
713  for (iSet=0;iSet<numberSets_;iSet++) {
714    int functionColumn = functionColumn_[iSet];
715    int iSeq=0;
716    int iNon;
717    for (iNon =startSet_[iSet];iNon<startSet_[iSet+1];iNon++) {
718      if (definitionRowL_[iNon]>=0)
719        break;
720    }
721    int startNon=iNon;
722    int number = numberInSet[iSet];
723    double lo = lower[functionColumn];
724    double up = upper[functionColumn];
725    double g = granularity[iSet];
726    double valueZ=lo;
727    for (int i=0;i<number;i++) {
728      for (int iNon =startNon;iNon<startSet_[iSet+1];iNon++) {
729        // for Link
730        weight[iSeq]=valueZ;
731        sequence[iSeq++]=n+numberColumns;
732        for (int i=0;i<2;i++) {
733          if (iNon==startNon) {
734            // linear terms
735            for (int iNon =startSet_[iSet];iNon<startNon;iNon++) {
736              row[nel]=usageRow_[iNon];
737              element[nel++]=1.0;
738            }
739          }
740          // convexity
741          row[nel]=convexityRow_[iNon];
742          element[nel++]=1.0;
743          // definition of z
744          assert (definitionRow_[iNon]>=0);
745          row[nel]=definitionRow_[iNon];
746          element[nel++]=1.0;
747          if (definitionRowL_[iNon]) {
748            row[nel]=definitionRowL_[iNon];
749            element[nel++]=1.0;
750          }
751          row[nel]=usageRow_[iNon];
752          element[nel++]=1.0;
753          n++;
754          start[n]=nel;
755        }
756      }
757      valueZ += g;
758      if (fabs(valueZ-up)<1.0e-7)
759        valueZ=up;
760    }
761    int nLink = iSeq/number;
762    objects[iSet]=new CbcLink(NULL,number,nLink,setType[iSet],sequence,weight,iSet);
763    objects[iSet]->setPriority(1000);
764    objects[iSet]->setPriority(10000);
765   
766  }
767  newClp->addColumns(n,newLower,newUpper,newObjective,
768                     start,row,element);
769  delete [] newLower;
770  delete [] newUpper;
771  delete [] newObjective;
772  delete [] start;
773  delete [] row;
774  delete [] element;
775  lower = newClp->columnLower();
776  upper = newClp->columnUpper();
777  // pack down
778  CoinPackedMatrix * matrix = newClp->matrix();
779  matrix->removeGaps();
780  row = matrix->getMutableIndices();
781  start = matrix->getMutableVectorStarts();
782  element = matrix->getMutableElements();
783  // position
784  n = numberColumns;
785  nel = start[n];
786  // Now do real work
787  for (iSet=0;iSet<numberSets_;iSet++) {
788    int functionColumn = functionColumn_[iSet];
789    const char * name = coinModel_.getColumnName(functionColumn);
790    int iSeq=0;
791    int iNon;
792    for (iNon =startSet_[iSet];iNon<startSet_[iSet+1];iNon++) {
793      if (definitionRowL_[iNon]>=0)
794        break;
795    }
796    int startNon=iNon;
797    int number = numberInSet[iSet];
798    double lo = lower[functionColumn];
799    double up = upper[functionColumn];
800    double g = granularity[iSet];
801    double valueZ=lo;
802    for (int i=0;i<number;i++) {
803      for (int iNon =startNon;iNon<startSet_[iSet+1];iNon++) {
804        // for Link
805        weight[iSeq]=valueZ;
806        coinModel_.associateElement(name,valueZ);
807        sequence[iSeq++]=n+numberColumns;
808        for (int i=0;i<2;i++) {
809          if (iNon==startNon) {
810            // linear terms
811            for (int iNon =startSet_[iSet];iNon<startNon;iNon++) {
812              row[nel]=usageRow_[iNon];
813              element[nel++]=coinModel_.getElement(usageRow_[iNon],functionColumn)*valueZ;
814            }
815          }
816          // convexity
817          row[nel]=convexityRow_[iNon];
818          element[nel++]=1.0;
819          // definition of z
820          assert (definitionRow_[iNon]>=0);
821          row[nel]=definitionRow_[iNon];
822          element[nel++]=valueZ;
823          if (definitionRowL_[iNon]) {
824            row[nel]=definitionRowL_[iNon];
825            int jColumn = definitionColumnL_[iNon];
826            int jVar = back[jColumn];
827            assert (jVar>=0);
828            double bound = (!i) ? lower[jColumn] :upper[jColumn];
829            info_[jVar].addCoefficientModifier((i!=0),nel,1.0);
830            element[nel++]=bound;
831            row[nel]=usageRow_[iNon];
832#if 0
833            CoinModel dummy;
834            dummy.setElement(0,0,coinModel_.getColumnName(jColumn));
835            dummy.setElement(0,1,coinModel_.getElementAsString(row[nel],functionColumn));
836            double associated[2];
837            associated[0]=valueZ;
838            associated[1]=dummy.unsetValue();
839            dummy.computeAssociated(associated);
840            double value = associated[1];
841#endif
842            double value = getFunctionValueFromString(coinModel_.getElementAsString(row[nel],functionColumn),
843                                               name,valueZ);
844            info_[jVar].addCoefficientModifier((i!=0),nel,value);
845            element[nel++]=bound*value;
846          } else {
847            // just z
848            row[nel]=usageRow_[iNon];
849            CoinModel dummy;
850            dummy.setElement(0,0,name);
851            dummy.setElement(0,1,coinModel_.getElementAsString(row[nel],functionColumn));
852            double associated[2];
853            associated[0]=valueZ;
854            associated[1]=dummy.unsetValue();
855            dummy.computeAssociated(associated);
856            double value = associated[1];
857            element[nel++]=value;
858            element[nel++]=1.0;
859          }
860          n++;
861          start[n]=nel;
862        }
863      }
864      valueZ += g;
865      if (fabs(valueZ-up)<1.0e-7)
866        valueZ=up;
867    }
868  }
869  matrix_ = new CoinPackedMatrix(*matrix);
870  delete [] weight;
871  delete [] sequence;
872  delete [] numberInSet;
873  delete [] setType;
874  delete [] granularity;
875  // now load up model
876  assert (!savedModel_);
877  ClpSimplex * saved = modelPtr_;
878  modelPtr_=newClp;
879  OsiSolverInterface * solver = this->clone();
880  savedModel_ = saved;
881  model_->assignSolver(solver);
882  for (iSet=0;iSet<numberSets_;iSet++) {
883    objects[iSet]->setModel(model_);
884  }
885  model_->addObjects(numberSets_,objects);
886  for (iSet=0;iSet<numberSets_;iSet++) {
887    delete objects[iSet];
888  }
889  delete [] objects;
890}
891// revert to previous solver
892void 
893CbcSolverLink::revert()
894{
895  delete modelPtr_;
896  modelPtr_=savedModel_;
897  savedModel_=NULL;
898}
899// Add a bound modifier
900void 
901CbcSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected, 
902                                double multiplier)
903{
904  bool found=false;
905  int i;
906  for ( i=0;i<numberVariables_;i++) {
907    if (info_[i].variable()==whichVariable) {
908      found=true;
909      break;
910    }
911  }
912  if (!found) {
913    // add in
914    CbcLinkedBound * temp = new CbcLinkedBound [numberVariables_+1];
915    for (int i=0;i<numberVariables_;i++) 
916      temp[i]= info_[i];
917    delete [] info_;
918    info_=temp;
919    info_[numberVariables_++] = CbcLinkedBound(model_,whichVariable,0,NULL,NULL,NULL);
920  }
921  info_[i].addBoundModifier(upperBoundAffected, useUpperBound,whichVariableAffected,multiplier);
922}
923// Sets granularity
924void 
925CbcSolverLink::setGranularity(double top, double bottom, const double * topG, const double * bottomG)
926{
927  delete [] granularityTop_;
928  delete [] granularityBottom_;
929  granularityTop_ = NULL;
930  granularityBottom_ = NULL;
931  if (numberSets_>0) {
932    if (topG) {
933      granularityTop_ = CoinCopyOfArray(topG,numberSets_);
934    } else {
935      granularityTop_ = new double [numberSets_];
936      for (int i=0;i<numberSets_;i++) {
937        granularityTop_[i]=top;
938      }
939    }
940    if (bottomG) {
941      granularityBottom_ = CoinCopyOfArray(bottomG,numberSets_);
942    } else {
943      granularityBottom_ = new double [numberSets_];
944      for (int i=0;i<numberSets_;i++) {
945        granularityBottom_[i]=bottom;
946      }
947    }
948  }
949}
950//#############################################################################
951// Constructors, destructors  and assignment
952//#############################################################################
953
954//-------------------------------------------------------------------
955// Default Constructor
956//-------------------------------------------------------------------
957CbcLinkedBound::CbcLinkedBound ()
958{
959  model_ = NULL;
960  variable_ = -1;
961  numberAffected_ = 0;
962  maximumAffected_ = numberAffected_;
963  affected_ = NULL;
964  useObject_=-1;
965  numberSets_=0;
966  whichSets_ = NULL;
967}
968// Useful Constructor
969CbcLinkedBound::CbcLinkedBound(CbcModel * model, int variable,
970                               int numberAffected, const int * positionL, 
971                               const int * positionU, const double * multiplier)
972{
973  model_ = model;
974  variable_ = variable;
975  numberAffected_ = 2*numberAffected;
976  maximumAffected_ = numberAffected_;
977  numberSets_ = 0;
978  whichSets_ = NULL;
979  useObject_=-1;
980  if (numberAffected_) { 
981    affected_ = new boundElementAction[numberAffected_];
982    int n=0;
983    for (int i=0;i<numberAffected;i++) {
984      // LB
985      boundElementAction action;
986      action.affect=2;
987      action.ubUsed=0;
988      action.type=0;
989      action.affected=positionL[i];
990      action.multiplier=multiplier[i];
991      affected_[n++]=action;
992      // UB
993      action.affect=2;
994      action.ubUsed=1;
995      action.type=0;
996      action.affected=positionU[i];
997      action.multiplier=multiplier[i];
998      affected_[n++]=action;
999    }
1000  } else {
1001    affected_ = NULL;
1002  }
1003}
1004
1005//-------------------------------------------------------------------
1006// Copy constructor
1007//-------------------------------------------------------------------
1008CbcLinkedBound::CbcLinkedBound (
1009                  const CbcLinkedBound & rhs)
1010{
1011  model_ = rhs.model_;
1012  variable_ = rhs.variable_;
1013  numberAffected_ = rhs.numberAffected_;
1014  maximumAffected_ = rhs.maximumAffected_;
1015  numberSets_ = rhs.numberSets_;
1016  useObject_=rhs.useObject_;
1017  if (numberAffected_) { 
1018    affected_ = new boundElementAction[maximumAffected_];
1019    memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
1020    whichSets_ = CoinCopyOfArray(rhs.whichSets_,numberSets_);
1021  } else {
1022    affected_ = NULL;
1023    whichSets_ = NULL;
1024  }
1025}
1026
1027//-------------------------------------------------------------------
1028// Destructor
1029//-------------------------------------------------------------------
1030CbcLinkedBound::~CbcLinkedBound ()
1031{
1032  delete [] affected_;
1033  delete [] whichSets_;
1034}
1035
1036//-------------------------------------------------------------------
1037// Assignment operator
1038//-------------------------------------------------------------------
1039CbcLinkedBound &
1040CbcLinkedBound::operator=(const CbcLinkedBound& rhs)
1041{
1042  if (this != &rhs) { 
1043    delete [] affected_;
1044    delete [] whichSets_;
1045    model_ = rhs.model_;
1046    variable_ = rhs.variable_;
1047    numberAffected_ = rhs.numberAffected_;
1048    maximumAffected_ = rhs.maximumAffected_;
1049    numberSets_ = rhs.numberSets_;
1050    useObject_=rhs.useObject_;
1051    if (numberAffected_) { 
1052      affected_ = new boundElementAction[maximumAffected_];
1053      memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
1054      whichSets_ = CoinCopyOfArray(rhs.whichSets_,numberSets_);
1055    } else {
1056      affected_ = NULL;
1057      whichSets_ = NULL;
1058    }
1059  }
1060  return *this;
1061}
1062// Add a bound modifier
1063void 
1064CbcLinkedBound::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, 
1065                                 double multiplier)
1066{
1067  if (numberAffected_==maximumAffected_) {
1068    maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
1069    boundElementAction * temp = new boundElementAction[maximumAffected_];
1070    memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
1071    delete [] affected_;
1072    affected_ = temp;
1073  }
1074  boundElementAction action;
1075  action.affect=upperBoundAffected ? 1 : 0;
1076  action.ubUsed=useUpperBound ? 1 : 0;
1077  action.type=2;
1078  action.affected=whichVariable;
1079  action.multiplier=multiplier;
1080  affected_[numberAffected_++]=action;
1081 
1082}
1083// Update other bounds
1084void 
1085CbcLinkedBound::updateBounds(ClpSimplex * solver)
1086{
1087  double * lower = solver->columnLower();
1088  double * upper = solver->columnUpper();
1089  double lo = lower[variable_];
1090  double up = upper[variable_];
1091  // printf("bounds for %d are %g and %g\n",variable_,lo,up);
1092  for (int j=0;j<numberAffected_;j++) {
1093    if (affected_[j].affect<2) {
1094      double multiplier = affected_[j].multiplier;
1095      assert (affected_[j].type==2);
1096      int iColumn = affected_[j].affected;
1097      double useValue = (affected_[j].ubUsed) ? up : lo;
1098      if (affected_[j].affect==0) 
1099        lower[iColumn] = CoinMin(upper[iColumn],CoinMax(lower[iColumn],multiplier*useValue));
1100      else
1101        upper[iColumn] = CoinMax(lower[iColumn],CoinMin(upper[iColumn],multiplier*useValue));
1102    }
1103  }
1104}
1105// Add an element modifier
1106void 
1107CbcLinkedBound::addCoefficientModifier(bool useUpperBound, int position, 
1108                                 double multiplier)
1109{
1110  if (numberAffected_==maximumAffected_) {
1111    maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
1112    boundElementAction * temp = new boundElementAction[maximumAffected_];
1113    memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
1114    delete [] affected_;
1115    affected_ = temp;
1116  }
1117  boundElementAction action;
1118  action.affect=2;
1119  action.ubUsed=useUpperBound ? 1 : 0;
1120  action.type=0;
1121  action.affected=position;
1122  action.multiplier=multiplier;
1123  affected_[numberAffected_++]=action;
1124 
1125}
1126// Update coefficients
1127void 
1128CbcLinkedBound::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
1129{
1130  double * lower = solver->columnLower();
1131  double * upper = solver->columnUpper();
1132  double * element = matrix->getMutableElements();
1133  double lo = lower[variable_];
1134  double up = upper[variable_];
1135  // printf("bounds for %d are %g and %g\n",variable_,lo,up);
1136  for (int j=0;j<numberAffected_;j++) {
1137    if (affected_[j].affect==2) {
1138      double multiplier = affected_[j].multiplier;
1139      assert (affected_[j].type==0);
1140      int position = affected_[j].affected;
1141      //double old = element[position];
1142      if (affected_[j].ubUsed)
1143        element[position] = multiplier*up;
1144      else
1145        element[position] = multiplier*lo;
1146      //if ( old != element[position])
1147      //printf("change at %d from %g to %g\n",position,old,element[position]);
1148    }
1149  }
1150}
1151// Fix variables in LOS
1152void 
1153CbcLinkedBound::fixLOS(ClpSimplex * solver)
1154{
1155  abort();
1156}
1157// Update this variables bounds from reference row
1158void 
1159CbcLinkedBound::updateBoundsForThis(ClpSimplex * solver)
1160{
1161  if (useObject_>=0) {
1162    const OsiObject * object = model_->object(useObject_);
1163    const CbcLink * link = dynamic_cast<const CbcLink *> (object);
1164    assert (link);
1165    double * lower = solver->columnLower();
1166    double * upper = solver->columnUpper();
1167    double lo = lower[variable_];
1168    double up = upper[variable_];
1169    int numberMembers = link->numberMembers();
1170    int numberLinks = link->numberLinks();
1171    const int * which = link->which();
1172    const double * weights = link->weights();
1173    int i;
1174    double loWeight=-COIN_DBL_MAX;
1175    double upWeight=COIN_DBL_MAX;
1176    int base=0;
1177    for (i=0;i<numberMembers;i++) {
1178      int nFixed=0;
1179      for (int j=0;j<numberLinks;j++) {
1180        int iColumn = which[j+base];
1181        if (!upper[iColumn])
1182          nFixed++;
1183      }
1184      if (nFixed) {
1185        assert (nFixed==numberLinks);
1186        if (loWeight!=-COIN_DBL_MAX) {
1187          assert(i>0);
1188          upWeight=weights[i-1];
1189          break;
1190        }
1191      } else if (loWeight==-COIN_DBL_MAX) {
1192        loWeight=weights[i];
1193      }
1194    }
1195    upper[variable_] = CoinMin(up,upWeight);
1196    lower[variable_] = CoinMax(lo,loWeight);
1197    assert (upper[variable_]>=lower[variable_]);
1198  }
1199}
1200// Default Constructor
1201CbcHeuristicDynamic2::CbcHeuristicDynamic2() 
1202  :CbcHeuristic()
1203{
1204}
1205
1206// Constructor from model
1207CbcHeuristicDynamic2::CbcHeuristicDynamic2(CbcModel & model)
1208  :CbcHeuristic(model)
1209{
1210}
1211
1212// Destructor
1213CbcHeuristicDynamic2::~CbcHeuristicDynamic2 ()
1214{
1215}
1216
1217// Clone
1218CbcHeuristic *
1219CbcHeuristicDynamic2::clone() const
1220{
1221  return new CbcHeuristicDynamic2(*this);
1222}
1223
1224// Copy constructor
1225CbcHeuristicDynamic2::CbcHeuristicDynamic2(const CbcHeuristicDynamic2 & rhs)
1226:
1227  CbcHeuristic(rhs)
1228{
1229}
1230
1231// Returns 1 if solution, 0 if not
1232int
1233CbcHeuristicDynamic2::solution(double & solutionValue,
1234                         double * betterSolution)
1235{
1236  if (!model_)
1237    return 0;
1238  CbcSolverLink * clpSolver
1239    = dynamic_cast<CbcSolverLink *> (model_->solver());
1240  assert (clpSolver); 
1241  double newSolutionValue = clpSolver->bestObjectiveValue();
1242  const double * solution = clpSolver->bestSolution();
1243  if (newSolutionValue<solutionValue&&solution) {
1244    int numberColumns = clpSolver->getNumCols();
1245    // new solution
1246    memcpy(betterSolution,solution,numberColumns*sizeof(double));
1247    solutionValue = newSolutionValue;
1248    return 1;
1249  } else {
1250    return 0;
1251  }
1252}
1253// update model
1254void CbcHeuristicDynamic2::setModel(CbcModel * model)
1255{
1256  model_ = model;
1257}
1258// Resets stuff if model changes
1259void 
1260CbcHeuristicDynamic2::resetModel(CbcModel * model)
1261{
1262  model_ = model;
1263}
Note: See TracBrowser for help on using the repository browser.