source: branches/devel/Cbc/examples/OsiSolverLink.cpp @ 476

Last change on this file since 476 was 476, checked in by forrest, 14 years ago

osi stuff

File size: 19.7 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 "OsiSolverLink.hpp"
14#include "OsiBranchLink.hpp"
15#include "ClpPackedMatrix.hpp"
16#include "CoinTime.hpp"
17//#############################################################################
18// Solve methods
19//#############################################################################
20void OsiSolverLink::initialSolve()
21{
22  specialOptions_ =0;
23  modelPtr_->setWhatsChanged(0);
24  ClpMatrixBase * save = NULL;
25  if (numberVariables_) {
26    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
27    // update all bounds before coefficients
28    for (int i=0;i<numberVariables_;i++ ) {
29      info_[i].updateBounds(modelPtr_);
30    }
31    updateCoefficients(modelPtr_,temp);
32    temp->removeGaps(1.0e-14);
33    save = modelPtr_->clpMatrix();
34    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
35    assert (clpMatrix);
36    if (save->getNumRows()>temp->getNumRows()) {
37      // add in cuts
38      int numberRows = temp->getNumRows();
39      int * which = new int[numberRows];
40      for (int i=0;i<numberRows;i++)
41        which[i]=i;
42      save->deleteRows(numberRows,which);
43      delete [] which;
44      temp->bottomAppendPackedMatrix(*clpMatrix->matrix());
45    }
46    modelPtr_->replaceMatrix(temp);
47  }
48  OsiClpSolverInterface::initialSolve();
49  if (save) {
50    delete save;
51  }
52}
53//#define WRITE_MATRIX
54#ifdef WRITE_MATRIX
55static int xxxxxx=0;
56#endif
57//-----------------------------------------------------------------------------
58void OsiSolverLink::resolve()
59{
60  specialOptions_ =0;
61  modelPtr_->setWhatsChanged(0);
62  ClpMatrixBase * save = NULL;
63  bool allFixed=false;
64  bool feasible=true;
65  if (numberVariables_) {
66    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
67    allFixed=true;
68    //bool best=true;
69    const double * lower = modelPtr_->columnLower();
70    const double * upper = modelPtr_->columnUpper();
71    // update all bounds before coefficients
72    for (int i=0;i<numberVariables_;i++ ) {
73      info_[i].updateBounds(modelPtr_);
74      int iColumn = info_[i].variable();
75      double lo = lower[iColumn];
76      double up = upper[iColumn];
77      if (up>lo)
78        allFixed=false;
79      else if (up<lo)
80        feasible=false;
81    }
82    updateCoefficients(modelPtr_,temp);
83    temp->removeGaps(1.0e-14);
84    save = modelPtr_->clpMatrix();
85    ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
86    assert (clpMatrix);
87    if (save->getNumRows()>temp->getNumRows()) {
88      // add in cuts
89      int numberRows = temp->getNumRows();
90      int * which = new int[numberRows];
91      for (int i=0;i<numberRows;i++)
92        which[i]=i;
93      CoinPackedMatrix * mat = clpMatrix->matrix();
94      // for debug
95      //mat = new CoinPackedMatrix(*mat);
96      mat->deleteRows(numberRows,which);
97      delete [] which;
98      temp->bottomAppendPackedMatrix(*mat);
99      temp->removeGaps(1.0e-14);
100    }
101    if (0) {
102      const CoinPackedMatrix * matrix = modelPtr_->matrix();
103      int numberColumns = matrix->getNumCols();
104      assert (numberColumns==temp->getNumCols());
105      const double * element1 = temp->getMutableElements();
106      const int * row1 = temp->getIndices();
107      const CoinBigIndex * columnStart1 = temp->getVectorStarts();
108      const int * columnLength1 = temp->getVectorLengths();
109      const double * element2 = matrix->getMutableElements();
110      const int * row2 = matrix->getIndices();
111      const CoinBigIndex * columnStart2 = matrix->getVectorStarts();
112      const int * columnLength2 = matrix->getVectorLengths();
113      for (int i=0;i<numberColumns;i++) {
114        assert (columnLength2[i]==columnLength1[i]);
115        int offset = columnStart2[i]-columnStart1[i];
116        for (int j=columnStart1[i];j<columnStart1[i]+columnLength1[i];j++) {
117          assert (row1[j]==row2[j+offset]);
118          assert (element1[j]==element2[j+offset]);
119        }
120      }
121    }
122    modelPtr_->replaceMatrix(temp);
123  }
124#ifdef WRITE_MATRIX
125  {
126    xxxxxx++;
127    char temp[50];
128    sprintf(temp,"bb%d",xxxxxx);
129    writeMps(temp);
130  }
131#endif
132  if (!feasible)
133    allFixed=false;
134  if ((specialOptions2_&1)!=0)
135    allFixed=false;
136  int returnCode=-1;
137  if (feasible) {
138    // may do lots of work
139    returnCode=fathom(allFixed);
140  }
141  if (returnCode>=0) {
142    if (returnCode==0)
143      OsiClpSolverInterface::resolve();
144    if (!allFixed&&(specialOptions2_&1)==0) {
145      const double * solution = getColSolution();
146      bool satisfied=true;
147      for (int i=0;i<numberVariables_;i++) {
148        int iColumn = info_[i].variable();
149        double value = solution[iColumn];
150        if (fabs(value-floor(value+0.5))>0.0001)
151          satisfied=false;
152      }
153      //if (satisfied)
154      //printf("satisfied but not fixed\n");
155    }
156  } else {
157    modelPtr_->setProblemStatus(1);
158    modelPtr_->setObjectiveValue(COIN_DBL_MAX);
159  }
160  if (save) {
161    delete save;
162  }
163}
164
165//#############################################################################
166// Constructors, destructors clone and assignment
167//#############################################################################
168
169//-------------------------------------------------------------------
170// Default Constructor
171//-------------------------------------------------------------------
172OsiSolverLink::OsiSolverLink ()
173  : OsiClpSolverInterface()
174{
175  gutsOfDestructor(true);
176}
177#if 0
178/* returns
179   sequence of nonlinear or
180   -1 numeric
181   -2 not found
182   -3 too many terms
183*/
184static int getVariable(const CoinModel & model, char * expression,
185                       int & linear)
186{
187  int non=-1;
188  linear=-1;
189  if (strcmp(expression,"Numeric")) {
190    // function
191    char * first = strchr(expression,'*');
192    int numberColumns = model.numberColumns();
193    int j;
194    if (first) {
195      *first='\0';
196      for (j=0;j<numberColumns;j++) {
197        if (!strcmp(expression,model.columnName(j))) {
198          linear=j;
199          memmove(expression,first+1,strlen(first+1)+1);
200          break;
201        }
202      }
203    }
204    // find nonlinear
205    for (j=0;j<numberColumns;j++) {
206      const char * name = model.columnName(j);
207      first = strstr(expression,name);
208      if (first) {
209        if (first!=expression&&isalnum(*(first-1)))
210          continue; // not real match
211        first += strlen(name);
212        if (!isalnum(*first)) {
213          // match
214          non=j;
215          // but check no others
216          j++;
217          for (;j<numberColumns;j++) {
218            const char * name = model.columnName(j);
219            first = strstr(expression,name);
220            if (first) {
221              if (isalnum(*(first-1)))
222                continue; // not real match
223              first += strlen(name);
224              if (!isalnum(*first)) {
225                // match - ouch
226                non=-3;
227                break;
228              }
229            }
230          }
231          break;
232        }
233      }
234    }
235    if (non==-1)
236      non=-2;
237  }
238  return non;
239}
240#endif
241/* This creates from a coinModel object
242
243   if errors.then number of sets is -1
244     
245   This creates linked ordered sets information.  It assumes -
246
247   for product terms syntax is yy*f(zz)
248   also just f(zz) is allowed
249   and even a constant
250
251   modelObject not const as may be changed as part of process.
252*/
253OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
254  : OsiClpSolverInterface()
255{
256  gutsOfDestructor(true);
257  gdb(coinModel);
258}
259void OsiSolverLink::gdb ( CoinModel & coinModel)
260{
261  // first check and set up arrays
262  int numberColumns = coinModel.numberColumns();
263  int numberRows = coinModel.numberRows();
264  // List of nonlinear entries
265  int * which = new int[numberColumns];
266  numberVariables_=0;
267  specialOptions2_=0;
268  int iColumn;
269  int numberErrors=0;
270  for (iColumn=0;iColumn<numberColumns;iColumn++) {
271    CoinModelLink triple=coinModel.firstInColumn(iColumn);
272    bool linear=true;
273    int n=0;
274    while (triple.row()>=0) {
275      int iRow = triple.row();
276      const char * expr = coinModel.getElementAsString(iRow,iColumn);
277      if (strcmp(expr,"Numeric")) {
278        linear=false;
279      }
280      triple=coinModel.next(triple);
281      n++;
282    }
283    if (!linear) {
284      which[numberVariables_++]=iColumn;
285    }
286  }
287  // return if nothing
288  if (!numberVariables_) {
289    delete [] which;
290    return;
291  } else {
292    coinModel_ = coinModel;
293    int nBi=0;
294    int iRow;
295    for (iRow=0;iRow<numberRows;iRow++) {   
296      CoinModelLink triple=coinModel_.firstInRow(iRow);
297      while (triple.column()>=0) {
298        int iColumn = triple.column();
299        const char *  el = coinModel_.getElementAsString(iRow,iColumn);
300        if (strcmp("Numeric",el)) {
301          // check if value*x*y
302          char temp[100];
303          strcpy(temp,el);
304          char * ast = strchr(temp,'*');
305          if (ast) {
306            char * pos = temp;
307            while (pos!=ast) {
308              char x = *pos;
309              pos++;
310              assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-');
311            }
312            // must be column
313            assert (coinModel_.column(ast+1)>=0);
314          } else {
315            // must be column
316            assert (coinModel_.column(temp)>=0);
317          }
318          coinModel.setElement(iRow,iColumn,0.0);
319          nBi++;
320        }
321        triple=coinModel_.next(triple);
322      }
323    }
324    if (!nBi)
325      exit(1);
326    int nInt=0;
327    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
328      if (coinModel_.isInteger(iColumn))
329        nInt++;
330    }
331    printf("There are %d bilinear and %d integers\n",nBi,nInt);
332    loadFromCoinModel(coinModel,true);
333    OsiObject ** objects = new OsiObject * [nBi+nInt];
334    char * marked = new char [numberColumns];
335    memset(marked,0,numberColumns);
336    nBi=nInt;
337    for (iRow=0;iRow<numberRows;iRow++) {   
338      CoinModelLink triple=coinModel_.firstInRow(iRow);
339      while (triple.column()>=0) {
340        int iColumn = triple.column();
341        const char *  el = coinModel_.getElementAsString(iRow,iColumn);
342        if (strcmp("Numeric",el)) {
343          // value*x*y
344          char temp[100];
345          strcpy(temp,el);
346          char * ast = strchr(temp,'*');
347          double value=1.0;
348          if (ast) {
349            *ast='\0';
350            value = atof(temp);
351            ast++;
352          } else {
353            ast=temp;
354          }
355          // other column
356          int jColumn = coinModel_.column(ast);
357          double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
358          if (meshI)
359            marked[iColumn]=1;
360          double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
361          if (meshJ)
362            marked[jColumn]=1;
363          if (!meshJ&&!meshI) {
364            meshI=5.0e-1;
365            meshJ=5.0e-1;
366          }
367          objects[nBi] = new OsiBiLinear(this,iColumn,jColumn,iRow,value,meshI,meshJ,
368                                         nBi-nInt,(const OsiObject **) (objects+nInt));
369          objects[nBi]->setPriority(10000);
370          nBi++;
371        }
372        triple=coinModel_.next(triple);
373      }
374    }
375    nInt=0;
376    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
377      if (coinModel_.isInteger(iColumn)) {
378        objects[nInt] = new OsiSimpleInteger(this,iColumn);
379        if (marked[iColumn])
380          objects[nInt]->setPriority(100000);
381        else
382          objects[nInt]->setPriority(100);
383        nInt++;
384      }
385    }
386    nInt=nBi;
387    delete [] marked;
388    if (numberErrors) {
389      // errors
390      gutsOfDestructor();
391      numberVariables_=-1;
392    } else {
393      addObjects(nInt,objects);
394      int i;
395      for (i=0;i<nInt;i++)
396        delete objects[i];
397      delete [] objects;
398      // Now do dummy bound stuff
399      matrix_ = new CoinPackedMatrix(*getMatrixByCol());
400      info_ = new OsiLinkedBound [numberVariables_];
401      for ( i=0;i<numberVariables_;i++) {
402        info_[i] = OsiLinkedBound(this,which[i],0,NULL,NULL,NULL);
403      }
404    }
405  }
406}
407//-------------------------------------------------------------------
408// Clone
409//-------------------------------------------------------------------
410OsiSolverInterface * 
411OsiSolverLink::clone(bool copyData) const
412{
413  assert (copyData);
414  return new OsiSolverLink(*this);
415}
416
417
418//-------------------------------------------------------------------
419// Copy constructor
420//-------------------------------------------------------------------
421OsiSolverLink::OsiSolverLink (
422                  const OsiSolverLink & rhs)
423  : OsiClpSolverInterface(rhs)
424{
425  gutsOfDestructor(true);
426  gutsOfCopy(rhs);
427  // something odd happens - try this
428  OsiSolverInterface::operator=(rhs);
429}
430
431//-------------------------------------------------------------------
432// Destructor
433//-------------------------------------------------------------------
434OsiSolverLink::~OsiSolverLink ()
435{
436  gutsOfDestructor();
437}
438
439//-------------------------------------------------------------------
440// Assignment operator
441//-------------------------------------------------------------------
442OsiSolverLink &
443OsiSolverLink::operator=(const OsiSolverLink& rhs)
444{
445  if (this != &rhs) { 
446    gutsOfDestructor();
447    OsiClpSolverInterface::operator=(rhs);
448    gutsOfCopy(rhs);
449  }
450  return *this;
451}
452void 
453OsiSolverLink::gutsOfDestructor(bool justNullify)
454{
455  if (!justNullify) {
456    delete matrix_;
457    delete [] info_;
458  } 
459  matrix_ = NULL;
460  info_ = NULL;
461  numberVariables_ = 0;
462  specialOptions2_ = 0;
463}
464void 
465OsiSolverLink::gutsOfCopy(const OsiSolverLink & rhs)
466{
467  coinModel_ = rhs.coinModel_;
468  numberVariables_ = rhs.numberVariables_;
469  specialOptions2_ = rhs.specialOptions2_;
470  if (numberVariables_) { 
471    matrix_ = new CoinPackedMatrix(*rhs.matrix_);
472    info_ = new OsiLinkedBound [numberVariables_];
473    for (int i=0;i<numberVariables_;i++) {
474      info_[i] = OsiLinkedBound(rhs.info_[i]);
475    }
476  }
477}
478// Add a bound modifier
479void 
480OsiSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected, 
481                                double multiplier)
482{
483  bool found=false;
484  int i;
485  for ( i=0;i<numberVariables_;i++) {
486    if (info_[i].variable()==whichVariable) {
487      found=true;
488      break;
489    }
490  }
491  if (!found) {
492    // add in
493    OsiLinkedBound * temp = new OsiLinkedBound [numberVariables_+1];
494    for (int i=0;i<numberVariables_;i++) 
495      temp[i]= info_[i];
496    delete [] info_;
497    info_=temp;
498    info_[numberVariables_++] = OsiLinkedBound(this,whichVariable,0,NULL,NULL,NULL);
499  }
500  info_[i].addBoundModifier(upperBoundAffected, useUpperBound,whichVariableAffected,multiplier);
501}
502// Update coefficients
503void 
504OsiSolverLink::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
505{
506  double * lower = solver->columnLower();
507  double * upper = solver->columnUpper();
508  for (int iObject =0;iObject<numberObjects_;iObject++) {
509    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
510    if (obj) {
511      obj->updateCoefficients(lower,upper,matrix,&basis_);
512    }
513  }
514}
515//#############################################################################
516// Constructors, destructors  and assignment
517//#############################################################################
518
519//-------------------------------------------------------------------
520// Default Constructor
521//-------------------------------------------------------------------
522OsiLinkedBound::OsiLinkedBound ()
523{
524  model_ = NULL;
525  variable_ = -1;
526  numberAffected_ = 0;
527  maximumAffected_ = numberAffected_;
528  affected_ = NULL;
529}
530// Useful Constructor
531OsiLinkedBound::OsiLinkedBound(OsiSolverInterface * model, int variable,
532                               int numberAffected, const int * positionL, 
533                               const int * positionU, const double * multiplier)
534{
535  model_ = model;
536  variable_ = variable;
537  numberAffected_ = 2*numberAffected;
538  maximumAffected_ = numberAffected_;
539  if (numberAffected_) { 
540    affected_ = new boundElementAction[numberAffected_];
541    int n=0;
542    for (int i=0;i<numberAffected;i++) {
543      // LB
544      boundElementAction action;
545      action.affect=2;
546      action.ubUsed=0;
547      action.type=0;
548      action.affected=positionL[i];
549      action.multiplier=multiplier[i];
550      affected_[n++]=action;
551      // UB
552      action.affect=2;
553      action.ubUsed=1;
554      action.type=0;
555      action.affected=positionU[i];
556      action.multiplier=multiplier[i];
557      affected_[n++]=action;
558    }
559  } else {
560    affected_ = NULL;
561  }
562}
563
564//-------------------------------------------------------------------
565// Copy constructor
566//-------------------------------------------------------------------
567OsiLinkedBound::OsiLinkedBound (
568                  const OsiLinkedBound & rhs)
569{
570  model_ = rhs.model_;
571  variable_ = rhs.variable_;
572  numberAffected_ = rhs.numberAffected_;
573  maximumAffected_ = rhs.maximumAffected_;
574  if (numberAffected_) { 
575    affected_ = new boundElementAction[maximumAffected_];
576    memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
577  } else {
578    affected_ = NULL;
579  }
580}
581
582//-------------------------------------------------------------------
583// Destructor
584//-------------------------------------------------------------------
585OsiLinkedBound::~OsiLinkedBound ()
586{
587  delete [] affected_;
588}
589
590//-------------------------------------------------------------------
591// Assignment operator
592//-------------------------------------------------------------------
593OsiLinkedBound &
594OsiLinkedBound::operator=(const OsiLinkedBound& rhs)
595{
596  if (this != &rhs) { 
597    delete [] affected_;
598    model_ = rhs.model_;
599    variable_ = rhs.variable_;
600    numberAffected_ = rhs.numberAffected_;
601    maximumAffected_ = rhs.maximumAffected_;
602    if (numberAffected_) { 
603      affected_ = new boundElementAction[maximumAffected_];
604      memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
605    } else {
606      affected_ = NULL;
607    }
608  }
609  return *this;
610}
611// Add a bound modifier
612void 
613OsiLinkedBound::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, 
614                                 double multiplier)
615{
616  if (numberAffected_==maximumAffected_) {
617    maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
618    boundElementAction * temp = new boundElementAction[maximumAffected_];
619    memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
620    delete [] affected_;
621    affected_ = temp;
622  }
623  boundElementAction action;
624  action.affect=upperBoundAffected ? 1 : 0;
625  action.ubUsed=useUpperBound ? 1 : 0;
626  action.type=2;
627  action.affected=whichVariable;
628  action.multiplier=multiplier;
629  affected_[numberAffected_++]=action;
630 
631}
632// Update other bounds
633void 
634OsiLinkedBound::updateBounds(ClpSimplex * solver)
635{
636  double * lower = solver->columnLower();
637  double * upper = solver->columnUpper();
638  double lo = lower[variable_];
639  double up = upper[variable_];
640  // printf("bounds for %d are %g and %g\n",variable_,lo,up);
641  for (int j=0;j<numberAffected_;j++) {
642    if (affected_[j].affect<2) {
643      double multiplier = affected_[j].multiplier;
644      assert (affected_[j].type==2);
645      int iColumn = affected_[j].affected;
646      double useValue = (affected_[j].ubUsed) ? up : lo;
647      if (affected_[j].affect==0) 
648        lower[iColumn] = CoinMin(upper[iColumn],CoinMax(lower[iColumn],multiplier*useValue));
649      else
650        upper[iColumn] = CoinMax(lower[iColumn],CoinMin(upper[iColumn],multiplier*useValue));
651    }
652  }
653}
654#if 0
655// Add an element modifier
656void
657OsiLinkedBound::addCoefficientModifier(bool useUpperBound, int position,
658                                 double multiplier)
659{
660  if (numberAffected_==maximumAffected_) {
661    maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
662    boundElementAction * temp = new boundElementAction[maximumAffected_];
663    memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
664    delete [] affected_;
665    affected_ = temp;
666  }
667  boundElementAction action;
668  action.affect=2;
669  action.ubUsed=useUpperBound ? 1 : 0;
670  action.type=0;
671  action.affected=position;
672  action.multiplier=multiplier;
673  affected_[numberAffected_++]=action;
674 
675}
676// Update coefficients
677void
678OsiLinkedBound::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
679{
680  double * lower = solver->columnLower();
681  double * upper = solver->columnUpper();
682  double * element = matrix->getMutableElements();
683  double lo = lower[variable_];
684  double up = upper[variable_];
685  // printf("bounds for %d are %g and %g\n",variable_,lo,up);
686  for (int j=0;j<numberAffected_;j++) {
687    if (affected_[j].affect==2) {
688      double multiplier = affected_[j].multiplier;
689      assert (affected_[j].type==0);
690      int position = affected_[j].affected;
691      //double old = element[position];
692      if (affected_[j].ubUsed)
693        element[position] = multiplier*up;
694      else
695        element[position] = multiplier*lo;
696      //if ( old != element[position])
697      //printf("change at %d from %g to %g\n",position,old,element[position]);
698    }
699  }
700}
701#endif
Note: See TracBrowser for help on using the repository browser.