source: branches/devel/Cbc/examples/ClpDynamicInterface.cpp @ 409

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

now that I understand .am

File size: 35.6 KB
Line 
1// Copyright (C) 2006, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <cassert>
5
6#include "CoinHelperFunctions.hpp"
7//#include "CoinIndexedVector.hpp"
8//#include "ClpSimplex.hpp"
9#include "ClpDynamicInterface.hpp"
10#include "OsiClpSolverInterface.hpp"
11#include "CbcModel.hpp"
12#include "CglProbing.hpp"
13#include "CglKnapsackCover.hpp"
14#include "CbcStrategy.hpp"
15#include "OsiOpbdpSolve.hpp"
16//#############################################################################
17// Solve methods
18//#############################################################################
19void ClpDynamicInterface::initialSolve()
20{
21  OsiClpSolverInterface::initialSolve();
22}
23
24//-----------------------------------------------------------------------------
25void ClpDynamicInterface::resolve()
26{
27  if (cbcModel_) {
28    ClpDynamicInterface * solver
29      = dynamic_cast<ClpDynamicInterface *> (cbcModel_->continuousSolver());
30    if (solver==this) {
31      // use parent
32      OsiClpSolverInterface::resolve();
33      return;
34    }
35  }
36
37  printf("%d proposals\n",proposals_.getNumCols());
38  // build up model
39  ClpSimplex * model = new ClpSimplex(*staticModel_);
40  // Included are costed slacks (numberBlocks_ of them)
41  int numberCommon = staticModel_->numberColumns();
42  // make sure artificials can come in
43  int j;
44  for (j=numberCommon-numberArtificials_;j<numberCommon;j++) {
45    model->setColumnUpper(j,COIN_DBL_MAX);
46  }
47  numberProposalsAdded_=0;
48  setBasis(dynamicBasis_,model);
49  int numberColumns = modelPtr_->numberColumns();
50  const double * columnLower = modelPtr_->columnLower();
51  const double * columnUpper = modelPtr_->columnUpper();
52  // make sure bounds are good
53  int * backward = backward_[numberBlocks_];
54  double * columnLower2 = model->columnLower();
55  double * columnUpper2 = model->columnUpper();
56  int iColumn;
57  for (iColumn=0;iColumn<numberCommon-numberArtificials_;iColumn++) {
58    int jColumn = backward[iColumn];
59    columnLower2[iColumn]=columnLower[jColumn];
60    columnUpper2[iColumn]=columnUpper[jColumn];
61  }
62  double * solution = modelPtr_->primalColumnSolution();
63  int numberProposals=proposals_.getNumCols();
64  model->setLogLevel(0);
65  // here solve and generate
66  // save cutoff
67  double cutoff = modelPtr_->dualObjectiveLimit();
68  model->setDualObjectiveLimit(1.0e50);
69  model->dual();
70  // Big trouble if not feasible
71  assert (!model->problemStatus());
72  // add in existing valid ones
73  {
74    int * which = new int[numberProposals];
75    for (int i=0;i<numberProposals;i++)
76      which[i]=i;
77    addProposals(model,numberProposals,which,false);
78    delete [] which;
79  }
80  model->primal();
81  // generate
82  double lastObjective = model->objectiveValue();
83  double sumDj=0.0;
84  #define NPASS 5
85  #define MAXPASS 100
86  int maxPass=MAXPASS; // but no more than this
87  for (int iPass=0;iPass<NPASS;iPass++) {
88    maxPass--;
89    double * saveObjective = new double[numberColumns];
90    int numberRows2 = staticModel_->numberRows();
91    int numberMasterRows = numberRows2-numberBlocks_;
92    // we need big arrays for storing proposals
93    int size = modelPtr_->numberColumns()+numberRows2+3;
94    int * rows = new int[size];
95    double * elements = new double[size];
96    sumDj=0.0;
97    double * pi = new double[numberRows2];
98    for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
99      CoinMemcpyN(model->dualRowSolution(),numberRows2,pi);;
100      double * solution2 = model->primalColumnSolution();
101      int numberBad=0;
102      int i;
103      for ( i=numberCommon-numberArtificials_;i<numberCommon;i++) {
104        if (fabs(solution2[i])>1.0e-8) {
105          numberBad++;
106        }
107      }
108      int numberColumns2 = subProblem_[iBlock].numberColumns();
109      double * objective = subProblem_[iBlock].objective();
110      CoinMemcpyN(objective,numberColumns2,saveObjective);
111      masterRow_[iBlock].transposeTimes(pi,objective);
112      // and correct objective
113      for ( i=0;i<numberColumns2;i++) {
114        if (!numberBad)
115          objective[i] = saveObjective[i] - objective[i];
116        else
117          objective[i] = - objective[i];
118      }
119      // Now solve IP by some means
120      OsiClpSolverInterface solver (&subProblem_[iBlock],false);
121      // fix variables and clean objective
122      const int * backward = backward_[iBlock];
123      int nInt=0;
124      double * obj = solver.getModelPtr()->objective();
125      double maxObj=0.0;
126      for ( i=0;i<numberColumns2;i++) {
127        maxObj = CoinMax(maxObj,fabs(obj[i]));
128        int iColumn = backward[i];
129        if (solver.isInteger(i)&&!columnLower[i]&&columnUpper[i]==1.0)
130          nInt++;
131        solver.setColLower(i,columnLower[iColumn]);
132        solver.setColUpper(i,columnUpper[iColumn]);
133      }
134      if (maxObj>1.0e7) {
135        double scaleFactor = 1.0e7/maxObj;
136        for ( i=0;i<numberColumns2;i++) 
137          obj[i] *= scaleFactor;
138      }
139      if (nInt==numberColumns2&&!iPass&&!cbcModel_->getNodeCount())
140        printf("pure 0-1 problem\n");
141      //#define ENUMERATE
142#ifndef ENUMERATE
143      CbcModel small(solver);
144      // Normal type strategy
145      CbcStrategyDefault strategy(true,5,5);
146      small.setStrategy(strategy);
147      small.setLogLevel(0);
148      small.branchAndBound();
149      if (small.getMinimizationObjValue()<1.0e50) {
150        double * solution = small.bestSolution();
151        // clean
152        for (i=0;i<numberColumns2;i++) {
153          if (solver.isInteger(i))
154            solution[i]=floor(solution[i]+0.5);
155        }
156#else
157        int rcode=solveOpbdp(&solver);
158        if (rcode>0) {
159          const double * solution = solver.getColSolution();
160#endif
161        double objValue =0.0;
162        // proposal
163        masterRow_[iBlock].times(solution,elements);
164        for (i=0;i<numberColumns2;i++) {
165          objValue += solution[i]*saveObjective[i];
166        }
167        // See if good dj and pack down
168        int number=0;
169        double dj = objValue;
170        double smallest=1.0e100;
171        double largest=0.0;
172        for (i=0;i<numberMasterRows;i++) {
173          double value = elements[i];
174          if (fabs(value)>1.0e-15) {
175            dj -= pi[i]*value;
176            smallest = min(smallest,fabs(value));
177            largest = max(largest,fabs(value));
178            rows[number]=i;
179            elements[number++]=value;
180          }
181        }
182        // and convexity
183        dj -= pi[numberMasterRows+iBlock];
184        rows[number]=numberMasterRows+iBlock;
185        elements[number++]=1.0;
186        // if elements large then scale?
187        //if (largest>1.0e8||smallest<1.0e-8)
188        if (dj<0.0)
189          sumDj += dj;
190        if (dj<-1.0e-6) {
191          printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
192                 iBlock,smallest,largest,dj);
193          // take
194          // objective
195          int row = numberMasterRows+numberBlocks_;
196          rows[number]=row;
197          elements[number++]=objValue;
198          // what went into it
199          row++;
200          for (i=0;i<numberColumns2;i++) {
201            if (fabs(solution[i])>1.0e-15) {
202              rows[number]=row+i;
203              elements[number++]=solution[i];
204            }
205          }
206          // add
207          proposals_.appendCol(number,rows,elements);
208          const double * solution2 = model->primalColumnSolution();
209          // See if any artificials in
210          bool artificialsIn=false;
211          for (int j=numberCommon-numberArtificials_;j<numberCommon;j++) {
212            if (fabs(solution2[j])>1.0e-8||model->getStatus(j)==ClpSimplex::basic) {
213              artificialsIn=true;
214              break;
215            }
216          }
217          if ((maxPass>MAXPASS-5||artificialsIn)&&!cbcModel_->getNodeCount()) {
218            // add in one at a time in case of symmetry
219            int n=proposals_.getNumCols();
220            int * which = new int[n-numberProposals];
221            for (int i=numberProposals;i<n;i++)
222              which[i-numberProposals]=i;
223            addProposals(model,n-numberProposals,which,false);
224            delete [] which;
225            numberProposals=n;
226            //model->setLogLevel(1);
227            model->primal(1);
228            assert (!model->problemStatus());
229            double * solution2 = model->primalColumnSolution();
230            // Take out any artificials if we can
231            for (int j=numberCommon-numberArtificials_;j<numberCommon;j++) {
232              if (fabs(solution2[j])<1.0e-8&&model->getStatus(j)==ClpSimplex::basic) {
233                solution2[j]=0.0;
234                model->setStatus(j,ClpSimplex::atLowerBound);
235                model->setColumnUpper(j,0.0);
236              }
237            }
238          }
239        }
240      }
241      CoinMemcpyN(saveObjective,numberColumns2,objective);
242    }
243    delete [] pi;
244    delete [] rows;
245    delete [] elements;
246    delete [] saveObjective;
247    if ((sumDj>-1.0e-5||iPass==NPASS-1)) {
248      // exit
249      break;
250    } else {
251      int n=proposals_.getNumCols();
252      int * which = new int[n-numberProposals];
253      for (int i=numberProposals;i<n;i++)
254        which[i-numberProposals]=i;
255      addProposals(model,n-numberProposals,which,false);
256      delete [] which;
257      numberProposals=n;
258      model->setLogLevel(1);
259      // Take out any artificials if we can
260      double * solution2 = model->primalColumnSolution();
261      for (int j=numberCommon-numberArtificials_;j<numberCommon;j++) {
262        if (fabs(solution2[j])<1.0e-8&&model->getStatus(j)==ClpSimplex::basic) {
263          solution2[j]=0.0;
264          model->setStatus(j,ClpSimplex::atLowerBound);
265          model->setColumnUpper(j,0.0);
266        }
267      }
268      model->primal(1);
269      // Check if we have valid solution
270      int status=setSolution(model);
271      if (status<0)
272        printf("*** we have a solution of %g\n",model->objectiveValue());
273      if (lastObjective<model->objectiveValue()+1.0e-5&&cbcModel_->getNodeCount())
274        iPass=NPASS-2; // just generate and exit
275      else if (!cbcModel_->getNodeCount()&&status==1)
276        iPass=0; // not feasible - make sure we keep going
277      else if (!cbcModel_->getNodeCount()&&(sumDj<-1.0e-3&&maxPass>0))
278        iPass=0; // not optimal - make sure we keep going
279      lastObjective = model->objectiveValue();
280    }
281  }
282  printf("DW Objective value is %g - best possible %g\n",model->objectiveValue(),
283         model->objectiveValue()+sumDj);
284  dynamicBasis_ = getBasis(model,proposalsAdded_,numberCommon);
285  int numberRows = modelPtr_->numberRows();
286  int problemStatus = model->problemStatus();
287  CoinAssert (model->problemStatus()||model->objectiveValue()<1.0e50);
288  if (model->objectiveValue()+sumDj>cutoff) {
289    // say infeasible
290    problemStatus=1;
291  }
292  modelPtr_->setObjectiveValue(model->objectiveValue()+sumDj);
293  modelPtr_->setSumDualInfeasibilities(model->sumDualInfeasibilities());
294  modelPtr_->setNumberDualInfeasibilities(model->numberDualInfeasibilities());
295  modelPtr_->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities());
296  modelPtr_->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities());
297  if (!problemStatus) {
298    problemStatus=setSolution(model);
299    if (problemStatus<0) {
300      printf("*** We have a valid solution of %g\n",model->objectiveValue());
301      problemStatus=0;
302    }
303  }
304  //  If at root node try BAB
305  if (!cbcModel_->getNodeCount()) {
306    // build up model
307    ClpSimplex * model = new ClpSimplex(*staticModel_);
308    // Included are costed slacks (numberBlocks_ of them)
309    int numberCommon = staticModel_->numberColumns();
310    numberProposalsAdded_=0;
311    const double * columnLower = modelPtr_->columnLower();
312    const double * columnUpper = modelPtr_->columnUpper();
313    // make sure bounds are good
314    int * backward = backward_[numberBlocks_];
315    double * columnLower2 = model->columnLower();
316    double * columnUpper2 = model->columnUpper();
317    int iColumn;
318    for (iColumn=0;iColumn<numberCommon-numberArtificials_;iColumn++) {
319      int jColumn = backward[iColumn];
320      columnLower2[iColumn]=columnLower[jColumn];
321      columnUpper2[iColumn]=columnUpper[jColumn];
322    }
323    int numberProposals=proposals_.getNumCols();
324    // add in existing valid ones
325    {
326      int * which = new int[numberProposals];
327      int i;
328      for ( i=0;i<numberProposals;i++)
329        which[i]=i;
330      addProposals(model,numberProposals,which,false);
331      // make integer
332      for ( i=numberCommon;i<model->numberColumns();i++) {
333        model->setColLower(i,0.0);
334        model->setColUpper(i,1.0);
335        model->setInteger(i);
336      }
337      delete [] which;
338    }
339    OsiClpSolverInterface * clpSolver1 = new OsiClpSolverInterface(model,false);
340    CbcModel small(*clpSolver1);
341    clpSolver1->releaseClp();
342    delete clpSolver1;
343    small.setLogLevel(0);
344    CbcStrategyDefault strategy(true,5,5);
345    small.setStrategy(strategy);
346    small.branchAndBound();
347    if (small.getMinimizationObjValue()<1.0e50) {
348      // save solution round this
349      double * saveSolution = CoinCopyOfArray(modelPtr_->primalColumnSolution(),
350                                         modelPtr_->numberColumns());
351      //const double * solution = small.bestSolution();
352      OsiSolverInterface * solver1 = small.solver();
353      OsiClpSolverInterface * clpSolver
354        = dynamic_cast<OsiClpSolverInterface *> (solver1);
355      assert (clpSolver);
356      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
357      int status=setSolution(clpSimplex);
358      if (status<0)
359        printf("*** we have a solution of %g\n",clpSimplex->objectiveValue());
360      // copy back solution
361      CoinMemcpyN(saveSolution,
362                  modelPtr_->numberColumns(),
363                  modelPtr_->primalColumnSolution());
364      delete [] saveSolution ;
365    }
366    delete model;
367  }
368  modelPtr_->setProblemStatus(problemStatus);
369  // Don't allow reduced cost fixing?
370  CoinZeroN(modelPtr_->dualColumnSolution(),numberColumns);
371  memset(modelPtr_->primalRowSolution(),0,numberRows*sizeof(double));
372  modelPtr_->clpMatrix()->times(1.0,solution,modelPtr_->primalRowSolution());
373  modelPtr_->setNumberIterations(model->numberIterations());
374  delete model;
375}
376
377//#############################################################################
378// Constructors, destructors clone and assignment
379//#############################################################################
380
381//-------------------------------------------------------------------
382// Default Constructor
383//-------------------------------------------------------------------
384ClpDynamicInterface::ClpDynamicInterface ()
385  : OsiClpSolverInterface()
386{
387  numberBlocks_=0;
388  numberArtificials_ = 0;
389  staticModel_=NULL;
390  subProblem_=NULL;
391  masterRow_=NULL;
392  numberProposalsAdded_ = 0;
393  maximumProposalsAdded_ =0;
394  proposalsAdded_ = NULL;
395  backward_=NULL;
396  allGenerated_ = NULL;
397  cbcModel_ = NULL;
398  bestObjectiveValue_ =1.0e100;
399  bestSolution_=NULL;
400}
401
402//-------------------------------------------------------------------
403// Clone
404//-------------------------------------------------------------------
405OsiSolverInterface * 
406ClpDynamicInterface::clone(bool CopyData) const
407{
408  if (CopyData) {
409    return new ClpDynamicInterface(*this);
410  } else {
411    printf("warning ClpDynamicInterface clone with copyData false\n");
412    return new ClpDynamicInterface();
413  }
414}
415
416
417//-------------------------------------------------------------------
418// Copy constructor
419//-------------------------------------------------------------------
420ClpDynamicInterface::ClpDynamicInterface (
421                  const ClpDynamicInterface & rhs)
422  : OsiClpSolverInterface(rhs)
423{
424  dynamicBasis_ = rhs.dynamicBasis_;
425  numberBlocks_ = rhs.numberBlocks_;
426  numberArtificials_ = rhs.numberArtificials_;
427  numberProposalsAdded_ = rhs.numberProposalsAdded_;
428  maximumProposalsAdded_ =rhs.maximumProposalsAdded_;
429  bestObjectiveValue_ = rhs.bestObjectiveValue_;
430  cbcModel_ = rhs.cbcModel_;
431  if (rhs.staticModel_) {
432    staticModel_=new ClpSimplex(*rhs.staticModel_);
433    subProblem_=new ClpSimplex[numberBlocks_];
434    masterRow_ = new CoinPackedMatrix[numberBlocks_];
435    if (numberBlocks_) {
436      backward_ = new int * [numberBlocks_+1];
437      backward_[numberBlocks_]=CoinCopyOfArray(rhs.backward_[numberBlocks_],
438                                               staticModel_->numberColumns()-numberArtificials_);
439      for (int i=0;i<numberBlocks_;i++) {
440        subProblem_[i]=rhs.subProblem_[i];
441        masterRow_[i]=rhs.masterRow_[i];
442        backward_[i]=CoinCopyOfArray(rhs.backward_[i],subProblem_[i].numberColumns());
443      }
444    } else {
445      backward_=NULL;
446    }
447    proposalsAdded_ = CoinCopyOfArray(rhs.proposalsAdded_,maximumProposalsAdded_);;
448    proposals_=rhs.proposals_;
449    allGenerated_ = CoinCopyOfArray(rhs.allGenerated_,numberBlocks_);
450    if (rhs.bestSolution_) {
451      bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->getNumCols());
452    } else {
453      bestSolution_=NULL;
454    }
455  } else {
456    staticModel_=NULL;
457    subProblem_ = NULL;
458    masterRow_=NULL;
459    proposalsAdded_=NULL;
460    backward_=NULL;
461    allGenerated_ = NULL;
462    bestSolution_=NULL;
463  }
464}
465
466//-------------------------------------------------------------------
467// Destructor
468//-------------------------------------------------------------------
469ClpDynamicInterface::~ClpDynamicInterface ()
470{
471  delete staticModel_;
472  delete [] subProblem_;
473  delete [] masterRow_;
474  if (numberBlocks_) {
475    for (int i=0;i<numberBlocks_+1;i++) {
476      delete [] backward_[i];
477    }
478  }
479  delete [] proposalsAdded_;
480  delete [] backward_;
481  delete [] allGenerated_;
482  delete [] bestSolution_;
483}
484
485//-------------------------------------------------------------------
486// Assignment operator
487//-------------------------------------------------------------------
488ClpDynamicInterface &
489ClpDynamicInterface::operator=(const ClpDynamicInterface& rhs)
490{
491  if (this != &rhs) { 
492    OsiClpSolverInterface::operator=(rhs);
493    numberBlocks_ = rhs.numberBlocks_;
494    numberArtificials_ = rhs.numberArtificials_;
495    numberProposalsAdded_ = rhs.numberProposalsAdded_;
496    maximumProposalsAdded_ =rhs.maximumProposalsAdded_;
497    bestObjectiveValue_ = rhs.bestObjectiveValue_;
498    delete staticModel_;
499    delete [] subProblem_;
500    delete [] masterRow_;
501    if (numberBlocks_) {
502      for (int i=0;i<numberBlocks_+1;i++) {
503        delete [] backward_[i];
504      }
505    }
506    delete [] proposalsAdded_;
507    delete [] backward_;
508    delete [] allGenerated_;
509    delete [] bestSolution_;
510    dynamicBasis_ = rhs.dynamicBasis_;
511    proposalsAdded_ = CoinCopyOfArray(rhs.proposalsAdded_,maximumProposalsAdded_);;
512    cbcModel_ = rhs.cbcModel_;
513    if (rhs.staticModel_) {
514      staticModel_=new ClpSimplex(*rhs.staticModel_);
515      subProblem_=new ClpSimplex[numberBlocks_];
516      masterRow_ = new CoinPackedMatrix[numberBlocks_];
517      if (numberBlocks_) {
518        backward_ = new int * [numberBlocks_+1];
519        backward_[numberBlocks_]=CoinCopyOfArray(rhs.backward_[numberBlocks_],
520                                                 staticModel_->numberColumns()-numberArtificials_);
521        for (int i=0;i<numberBlocks_;i++) {
522          subProblem_[i]=rhs.subProblem_[i];
523          masterRow_[i]=rhs.masterRow_[i];
524          backward_[i]=CoinCopyOfArray(rhs.backward_[i],subProblem_[i].numberColumns());
525        }
526      } else {
527        backward_=NULL;
528      }
529      proposals_=rhs.proposals_;
530      allGenerated_ = CoinCopyOfArray(rhs.allGenerated_,numberBlocks_);
531      if (rhs.bestSolution_) {
532        bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->getNumCols());
533      } else {
534        bestSolution_=NULL;
535      }
536    } else {
537      staticModel_=NULL;
538      subProblem_ = NULL;
539      masterRow_=NULL;
540      backward_=NULL;
541      proposals_=CoinPackedMatrix();
542      allGenerated_ = NULL;
543      bestSolution_ = NULL;
544    }
545  }
546  return *this;
547}
548//-------------------------------------------------------------------
549// Real initializer
550//-------------------------------------------------------------------
551void
552ClpDynamicInterface::initialize (int * rowBlock, CbcModel * cbcModel)
553{
554  //
555  delete staticModel_;
556  delete [] subProblem_;
557  delete [] masterRow_;
558  cbcModel_ = cbcModel;
559  int numberColumns = modelPtr_->numberColumns();
560  int numberRows = modelPtr_->numberRows();
561  CoinPackedMatrix * matrix = modelPtr_->matrix();
562  int numberElements = matrix->getNumElements(); 
563  // get row copy
564  CoinPackedMatrix rowCopy = *matrix;
565  rowCopy.reverseOrdering();
566  const int * rowByColumn = matrix->getIndices();
567  const int * columnLength = matrix->getVectorLengths();
568  const CoinBigIndex * columnStart = matrix->getVectorStarts();
569  //const double * elementByColumn = matrix->getElements();
570  const int * column = rowCopy.getIndices();
571  const int * rowLength = rowCopy.getVectorLengths();
572  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
573  //const double * elementByRow = rowCopy.getElements();
574  numberBlocks_ = 0;
575  int * stack = new int [numberRows];
576  // to say if column looked at
577  int * columnBlock = new int[numberColumns];
578  int iRow,iColumn;
579  for (iColumn=0;iColumn<numberColumns;iColumn++)
580    columnBlock[iColumn]=-2;
581  for (iColumn=0;iColumn<numberColumns;iColumn++) {
582    int kstart = columnStart[iColumn];
583    int kend = columnStart[iColumn]+columnLength[iColumn];
584    if (columnBlock[iColumn]==-2) {
585      // column not allocated
586      int j;
587      int nstack=0;
588      for (j=kstart;j<kend;j++) {
589        int iRow= rowByColumn[j];
590        if (rowBlock[iRow]!=-1) {
591          assert(rowBlock[iRow]==-2);
592          rowBlock[iRow]=numberBlocks_; // mark
593          stack[nstack++] = iRow;
594        }
595      }
596      if (nstack) {
597        // new block - put all connected in
598        numberBlocks_++;
599        columnBlock[iColumn]=numberBlocks_-1;
600        while (nstack) {
601          int iRow = stack[--nstack];
602          int k;
603          for (k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {
604            int iColumn = column[k];
605            int kkstart = columnStart[iColumn];
606            int kkend = kkstart + columnLength[iColumn];
607            if (columnBlock[iColumn]==-2) {
608              columnBlock[iColumn]=numberBlocks_-1; // mark
609              // column not allocated
610              int jj;
611              for (jj=kkstart;jj<kkend;jj++) {
612                int jRow= rowByColumn[jj];
613                if (rowBlock[jRow]==-2) {
614                  rowBlock[jRow]=numberBlocks_-1;
615                  stack[nstack++]=jRow;
616                }
617              }
618            } else {
619              assert (columnBlock[iColumn]==numberBlocks_-1);
620            }
621          }
622        }
623      } else {
624        // Only in master
625        columnBlock[iColumn]=-1;
626      }
627    }
628  }
629  printf("%d blocks found\n",numberBlocks_);
630  delete [] stack;
631  int i;
632  subProblem_ = new ClpSimplex[numberBlocks_];
633  masterRow_ = new CoinPackedMatrix[numberBlocks_];
634  int * which = new int[numberColumns];
635  int * whichRow = new int[numberRows];
636  int * whichMaster = new int[numberRows];
637  CoinBigIndex * starts = new int [numberColumns+1];
638  int * row = new int[numberElements];
639  double * element = new double[numberElements];
640  // Master
641  int numberRows2;
642  int numberColumns2;
643  numberRows2=0;
644  numberColumns2=0;
645  for (i=0;i<numberRows;i++) {
646    if (rowBlock[i]==-1)
647      whichMaster[numberRows2++]=i;
648  }
649  // and do backward pointers
650  backward_ = new int * [numberBlocks_+1];
651  int numberMasterRows = numberRows2;
652  for (iColumn=0;iColumn<numberColumns;iColumn++) {
653    if (columnBlock[iColumn]==-1)
654      which[numberColumns2++]=iColumn;
655  }
656  backward_[numberBlocks_]=CoinCopyOfArray(which,numberColumns2);
657  staticModel_ = new ClpSimplex(modelPtr_,numberRows2,whichMaster,
658                                numberColumns2,which,true,false,false);
659  printf("Master has %d rows and %d columns\n",numberRows2,numberColumns2);
660  // Add convexity rows
661  for (i=0;i<numberBlocks_;i++) {
662    staticModel_->addRow(0,NULL,NULL,1.0,1.0);
663  }
664  // Add costed slacks to make sure feasible
665  // *** We could do better than this
666  const double * rowLower = staticModel_->rowLower();
667  const double * rowUpper = staticModel_->rowUpper();
668  numberArtificials_ = 0;
669  for (i=0;i<numberRows2;i++) {
670    if (rowLower[i]>-1.0e20) {
671      numberArtificials_ ++;
672      double value = 1.0;
673      staticModel_->addColumn(1,&i,&value,0.0,COIN_DBL_MAX,1.0e10);
674    }
675    if (rowUpper[i]<1.0e20) {
676      numberArtificials_ ++;
677      double value = -1.0;
678      staticModel_->addColumn(1,&i,&value,0.0,COIN_DBL_MAX,1.0e10);
679    }
680  }
681  for (i=0;i<numberBlocks_;i++) {
682    int row=i+numberMasterRows;
683    double value=1.0;
684    staticModel_->addColumn(1,&row,&value,0.0,COIN_DBL_MAX,1.0e10);
685  }
686  numberArtificials_ += numberBlocks_;
687  // Blocks
688  for (i=0;i<numberBlocks_;i++) {
689    numberRows2=0;
690    numberColumns2=0;
691    for (iRow=0;iRow<numberRows;iRow++) {
692      if (rowBlock[iRow]==i)
693        whichRow[numberRows2++]=iRow;
694    }
695    for (iColumn=0;iColumn<numberColumns;iColumn++) {
696    if (columnBlock[iColumn]==i)
697      which[numberColumns2++]=iColumn;
698    }
699    backward_[i]=CoinCopyOfArray(which,numberColumns2);
700    subProblem_[i] = ClpSimplex(modelPtr_,numberRows2,whichRow,
701                                numberColumns2,which,true,false,false);
702    printf("Subproblem %d has %d rows and %d columns\n",
703           i,numberRows2,numberColumns2);
704    masterRow_[i] = CoinPackedMatrix(*matrix,numberMasterRows,whichMaster,
705                                     numberColumns2,which);
706  }
707  delete [] starts;
708  delete [] row;
709  delete [] element;
710  delete [] which;
711  delete [] whichRow;
712  delete [] whichMaster;
713  delete [] columnBlock;
714  // Could see if we can generate all upfront
715  allGenerated_ = new int[numberBlocks_];
716  CoinZeroN(allGenerated_,numberBlocks_);
717}
718// Warm start
719CoinWarmStartBasisDynamic
720ClpDynamicInterface::getBasis(ClpSimplex * model, const int * whichColumns,
721                              int numberCommon) const
722{
723  int iRow,iColumn;
724  int numberRows = model->numberRows();
725  int numberColumns = model->numberColumns();
726  CoinWarmStartBasisDynamic basis;
727  int * dynamic = new int [numberColumns];
728  // compute number of columns
729  int numberColumns2=numberCommon;
730  const double * columnLower = model->columnLower();
731  for (iColumn=numberCommon;iColumn<numberColumns;iColumn++) {
732    int iStatus = model->getColumnStatus(iColumn);
733    if (columnLower[iColumn]||iStatus==1) 
734      numberColumns2++;
735  }
736  basis.setSize(numberColumns2,numberRows);
737  if (model->statusExists()) {
738    // Flip slacks
739    int lookupA[]={0,1,3,2,0,2};
740    for (iRow=0;iRow<numberRows;iRow++) {
741      int iStatus = model->getRowStatus(iRow);
742      iStatus = lookupA[iStatus];
743      basis.setArtifStatus(iRow,(CoinWarmStartBasis::Status) iStatus);
744    }
745    int numberColumns2=0;
746    int numberDynamic=0;
747    const double * columnLower = model->columnLower();
748    int lookupS[]={0,1,2,3,0,3};
749    for (iColumn=0;iColumn<numberColumns;iColumn++) {
750      int iStatus = model->getColumnStatus(iColumn);
751      if (columnLower[iColumn]||iStatus==1||iColumn<numberCommon) {
752        iStatus = lookupS[iStatus];
753        basis.setStructStatus(numberColumns2,(CoinWarmStartBasis::Status) iStatus);
754        if (iColumn>=numberCommon) 
755          dynamic[numberDynamic++]=whichColumns[iColumn-numberCommon];
756        numberColumns2++;
757      }
758    }
759    basis.setDynamicVariables(numberDynamic,dynamic);
760    delete [] dynamic;
761    basis.setNumCommonVariables(numberCommon);
762  }
763  return basis;
764}
765// Sets up basis
766void 
767ClpDynamicInterface::setBasis ( const CoinWarmStartBasisDynamic & basis,
768                                  ClpSimplex * model)
769{
770  // Return if null basis
771  if (!dynamicBasis_.getNumStructural())
772    return;
773  int numberCommon = basis.getNumCommonVariables();
774  assert (numberCommon==model->numberColumns());
775  int numberDynamic = basis.getNumDynamicVariables();
776  const int * dynamic = basis.getDynamicVariables();
777  // add in proposals
778  addProposals(model, numberDynamic, dynamic,true);
779  // Do rest of basis stuff
780  OsiClpSolverInterface::setBasis(basis,model);
781}
782// Adds proposals to model
783void 
784ClpDynamicInterface::addProposals(ClpSimplex * model, int number,
785                                  const int * which, bool addEvenIfFixed)
786{
787  const double * columnLower = modelPtr_->columnLower();
788  const double * columnUpper = modelPtr_->columnUpper();
789  double * cost = new double [number];
790  double * lower = new double [number];
791  double * upper = new double [number];
792  int numberProposals = proposals_.getNumCols();
793  // mark ones in
794  char * mark = new char[numberProposals];
795  memset(mark,0,numberProposals);
796  int i;
797  // re-allocate added if necessary
798  if (number+numberProposalsAdded_>maximumProposalsAdded_) {
799    int * temp = proposalsAdded_;
800    maximumProposalsAdded_ = 
801      CoinMax(2*maximumProposalsAdded_+100,number+numberProposalsAdded_);
802    proposalsAdded_ = new int[maximumProposalsAdded_];
803    CoinZeroN(proposalsAdded_+numberProposalsAdded_,
804              maximumProposalsAdded_-numberProposalsAdded_);
805    CoinMemcpyN(temp,numberProposalsAdded_,proposalsAdded_);
806    delete [] temp;
807  }
808  for (i=0;i<numberProposalsAdded_;i++) 
809    mark[proposalsAdded_[i]]=1;
810  int numberRows2 = model->numberRows();
811  int numberMasterRows = numberRows2-numberBlocks_;
812  int objRow=numberRows2;
813  int jColumn;
814  // Get counts of variables fixed to 1
815  int * numberAtOne = new int[numberBlocks_];
816  for (int i=0;i<numberBlocks_;i++) {
817    int * backward = backward_[i];
818    int n=subProblem_[i].numberColumns();
819    int count=0;
820    for (int j=0;j<n;j++) {
821      int iColumn = backward[j];
822      if (columnLower[iColumn])
823        count++;
824    }
825    numberAtOne[i]=count;
826  }
827  int numberElements = proposals_.getNumElements();
828  CoinBigIndex * starts = new int [number+1];
829  int * row = new int[numberElements];
830  double * element = new double[numberElements];
831  const int * rowByColumn = proposals_.getIndices();
832  const int * columnLength = proposals_.getVectorLengths();
833  const CoinBigIndex * columnStart = proposals_.getVectorStarts();
834  const double * elementByColumn = proposals_.getElements();
835  starts[0]=0;
836  CoinBigIndex numberElements2=0;
837  int numberColumns2=0;
838  for (jColumn=0;jColumn<number;jColumn++) {
839    int iColumn = which[jColumn];
840    // skip if already in
841    if (mark[iColumn])
842      continue;
843    double costValue=0.0;
844    int iBlock=-1;
845    int * backward = NULL;
846    bool out=false;
847    int numberOne=0;
848    for (CoinBigIndex j=columnStart[iColumn];
849         j<columnStart[iColumn]+columnLength[iColumn];j++) {
850      int iRow = rowByColumn[j];
851      if (iRow<objRow) {
852        row[numberElements2]=iRow;
853        element[numberElements2++]=elementByColumn[j];
854        if (iRow>=numberMasterRows) {
855          iBlock = iRow-numberMasterRows;
856          backward = backward_[iBlock];
857        }
858      } else if (iRow==objRow) {
859        costValue = elementByColumn[j];
860      } else {
861        assert(backward);
862        if (iBlock==8&&proposals_.getNumCols()==1493)
863          printf("block 8 row %d\n",iRow);
864        int kColumn =backward[iRow-objRow-1];
865        // only works for 0-1 at present
866        if (!columnUpper[kColumn]) {
867          out=true;
868          break;
869        } else if (columnLower[kColumn]) {
870          numberOne++;
871        }
872      }
873    }
874    // If does not match state of problem - discard
875    if (numberOne!=numberAtOne[iBlock])
876      out=true;
877    cost[numberColumns2]=costValue;
878    lower[numberColumns2]=0.0;
879    upper[numberColumns2]= out ? 0.0 : COIN_DBL_MAX;
880    starts[numberColumns2+1]=numberElements2;
881    // see if we take
882    if (!out||addEvenIfFixed) {
883      numberColumns2++;
884      mark[iColumn]=1;
885      proposalsAdded_[numberProposalsAdded_++]=iColumn;
886    } else {
887      // take off elements
888      numberElements2=starts[numberColumns2];
889    }
890  }
891  model->addColumns(numberColumns2,lower,upper,cost,
892                   starts,row,element);
893  delete [] starts;
894  delete [] row;
895  delete [] element;
896  delete [] cost;
897  delete [] lower;
898  delete [] upper;
899  delete [] numberAtOne;
900  delete [] mark;
901}
902/* Creates modelPtr_ solution.
903   Returns 1 if not feasible, -1 if integer solution
904   0 otherwise */
905int 
906ClpDynamicInterface::setSolution(ClpSimplex * model)
907{
908  int numberColumns = modelPtr_->numberColumns();
909  const double * columnLower = modelPtr_->columnLower();
910  const double * columnUpper = modelPtr_->columnUpper();
911  double * solution = modelPtr_->primalColumnSolution();
912  int iColumn;
913  for (iColumn=0;iColumn<numberColumns;iColumn++) {
914    // make solution sensible
915    solution[iColumn]=0.0;
916  }
917  const double * solution2 = model->primalColumnSolution();
918  int problemStatus=0;
919  int numberCommon = staticModel_->numberColumns();
920  // copy back common (without costed slacks)
921  int * backward = backward_[numberBlocks_];
922  for (iColumn=0;iColumn<numberCommon-numberArtificials_;iColumn++) {
923    solution[backward[iColumn]]= solution2[iColumn];
924  }
925  // infeasible if artificials in
926  for (;iColumn<numberCommon;iColumn++) {
927    if (fabs(solution2[iColumn])>1.0e-5)
928      problemStatus=1;
929  }
930  const int * rowByColumn = proposals_.getIndices();
931  const int * columnLength = proposals_.getVectorLengths();
932  const CoinBigIndex * columnStart = proposals_.getVectorStarts();
933  const double * elementByColumn = proposals_.getElements();
934  int numberMasterRows = model->numberRows()-numberBlocks_;
935  int objRow=model->numberRows();
936  int numberColumns2=model->numberColumns()-numberCommon;
937  for (int jColumn=0;jColumn<numberColumns2;jColumn++) {
938    int iColumn = proposalsAdded_[jColumn];
939    int iBlock=-1;
940    int * backward = NULL;
941    double value = solution2[jColumn+numberCommon];
942    if (value) {
943      for (CoinBigIndex j=columnStart[iColumn];
944           j<columnStart[iColumn]+columnLength[iColumn];j++) {
945        int iRow = rowByColumn[j];
946        if (iRow<objRow&&iRow>=numberMasterRows) {
947          iBlock = iRow-numberMasterRows;
948          backward = backward_[iBlock];
949        } else if (iRow>objRow) {
950          assert(backward);
951          int kColumn =backward[iRow-objRow-1];
952          solution[kColumn] += value*elementByColumn[j];
953        }
954      }
955    }
956  }
957  if (problemStatus==0) {
958    problemStatus=-1;
959    for (iColumn=0;iColumn<numberColumns;iColumn++) {
960      assert (solution[iColumn]>=columnLower[iColumn]-1.0e-7);
961      assert (solution[iColumn]<=columnUpper[iColumn]+1.0e-7);
962      if (modelPtr_->isInteger(iColumn)&&fabs(solution[iColumn]-floor(solution[iColumn]+0.5))>1.0e-7)
963        problemStatus=0;
964    }
965    if (problemStatus&&model->objectiveValue()<bestObjectiveValue_) {
966      delete [] bestSolution_;
967      bestSolution_ = CoinCopyOfArray(solution,modelPtr_->getNumCols());
968      bestObjectiveValue_ = model->objectiveValue();
969      int numberRows = modelPtr_->numberRows();
970      double * rowActivity = modelPtr_->primalRowSolution();
971      CoinZeroN(rowActivity,numberRows);
972      modelPtr_->times(1.0,bestSolution_,rowActivity);
973      const double * rowLower = modelPtr_->rowLower();
974      const double * rowUpper = modelPtr_->rowUpper();
975      bool good=true;
976      for (int i=0;i<numberRows;i++) {
977        if(rowActivity[i]<rowLower[i]-1.0e-5) {
978          printf("%d %g %g %g\n",i,rowLower[i],rowActivity[i],rowUpper[i]);
979          good=false;
980        } else if(rowActivity[i]>rowUpper[i]+1.0e-5) {
981          printf("%d %g %g %g\n",i,rowLower[i],rowActivity[i],rowUpper[i]);
982          good=false;
983        }
984      }
985      assert (good);
986    }
987  }
988  return problemStatus;
989}
990// Default Constructor
991CbcHeuristicDynamic::CbcHeuristicDynamic() 
992  :CbcHeuristic()
993{
994}
995
996// Constructor from model
997CbcHeuristicDynamic::CbcHeuristicDynamic(CbcModel & model)
998  :CbcHeuristic(model)
999{
1000}
1001
1002// Destructor
1003CbcHeuristicDynamic::~CbcHeuristicDynamic ()
1004{
1005}
1006
1007// Clone
1008CbcHeuristic *
1009CbcHeuristicDynamic::clone() const
1010{
1011  return new CbcHeuristicDynamic(*this);
1012}
1013
1014// Copy constructor
1015CbcHeuristicDynamic::CbcHeuristicDynamic(const CbcHeuristicDynamic & rhs)
1016:
1017  CbcHeuristic(rhs)
1018{
1019}
1020
1021// Returns 1 if solution, 0 if not
1022int
1023CbcHeuristicDynamic::solution(double & solutionValue,
1024                         double * betterSolution)
1025{
1026  if (!model_)
1027    return 0;
1028  ClpDynamicInterface * clpSolver
1029    = dynamic_cast<ClpDynamicInterface *> (model_->solver());
1030  assert (clpSolver); 
1031  double newSolutionValue = clpSolver->bestObjectiveValue();
1032  const double * solution = clpSolver->bestSolution();
1033  if (newSolutionValue<solutionValue&&solution) {
1034    int numberColumns = clpSolver->getNumCols();
1035    // new solution
1036    memcpy(betterSolution,solution,numberColumns*sizeof(double));
1037    solutionValue = newSolutionValue;
1038    return 1;
1039  } else {
1040    return 0;
1041  }
1042}
1043// update model
1044void CbcHeuristicDynamic::setModel(CbcModel * model)
1045{
1046  model_ = model;
1047}
1048// Resets stuff if model changes
1049void 
1050CbcHeuristicDynamic::resetModel(CbcModel * model)
1051{
1052  model_ = model;
1053}
Note: See TracBrowser for help on using the repository browser.