source: trunk/Cbc/src/CbcHeuristicDW.cpp @ 1951

Last change on this file since 1951 was 1951, checked in by forrest, 8 years ago

changes to fix bug on max nodes

File size: 59.9 KB
Line 
1// $Id: CbcHeuristicDW.cpp 1899 2013-04-09 18:12:08Z stefan $
2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6// edwin 12/5/09 carved out of CbcHeuristicRINS
7
8#if defined(_MSC_VER)
9// Turn off compiler warning about long names
10#  pragma warning(disable:4786)
11#endif
12#include <cassert>
13#include <cstdlib>
14#include <cmath>
15#include <cfloat>
16
17#include "OsiSolverInterface.hpp"
18#include "OsiClpSolverInterface.hpp"
19#include "CbcModel.hpp"
20#include "CbcMessage.hpp"
21#include "CbcHeuristicDW.hpp"
22#include "CbcStrategy.hpp"
23#include "ClpPresolve.hpp"
24#include "CglProbing.hpp"
25
26static int dummyCallBack(CbcHeuristicDW * /*heuristic*/, 
27                         CbcModel * /*thisModel*/ , int /*whereFrom*/)
28{
29    return 0;
30}
31
32// Default Constructor
33CbcHeuristicDW::CbcHeuristicDW()
34  : CbcHeuristic()
35{
36  setDefaults();
37}
38
39// Constructor with model - assumed before cuts
40
41CbcHeuristicDW::CbcHeuristicDW(CbcModel & model, int keepContinuous)
42  : CbcHeuristic(model)
43{
44  setDefaults();
45  functionPointer_ = dummyCallBack;
46  assert(model.solver());
47  solver_ = model.solver()->clone();
48  findStructure();
49}
50/* Constructor with model - assumed before cuts
51 */
52CbcHeuristicDW::CbcHeuristicDW (CbcModel & model,
53                                int callBack(CbcHeuristicDW * currentHeuristic,
54                                             CbcModel * thisModel,
55                                             int whereFrom),
56                                int keepContinuous)
57  : CbcHeuristic(model)
58{
59  setDefaults();
60  functionPointer_ = callBack;
61  assert(model.solver());
62  solver_ = model.solver()->clone();
63  findStructure();
64}
65// Set default values
66void 
67CbcHeuristicDW::setDefaults()
68{
69  targetObjective_ = -COIN_DBL_MAX;
70  bestObjective_ = COIN_DBL_MAX;
71  lastObjective_ = COIN_DBL_MAX;
72  fullDWEverySoOften_ = 0;
73  numberPasses_ = 0;
74  howOften_ = 100;
75  decayFactor_ = 0.5;
76  functionPointer_ = NULL;
77  solver_=NULL;
78  dwSolver_=NULL;
79  bestSolution_=NULL;
80  continuousSolution_ = NULL;
81  fixedDj_ = NULL;
82  saveLower_=NULL;
83  saveUpper_=NULL;
84  random_=NULL;
85  affinity_=NULL;
86  weights_=NULL;
87  objectiveDW_=NULL;
88  numberColumnsDW_=NULL;
89  whichRowBlock_=NULL;
90  whichColumnBlock_=NULL;
91  dwBlock_=NULL;
92  backwardRow_=NULL;
93  rowsInBlock_=NULL;
94  columnsInBlock_=NULL;
95  startRowBlock_=NULL;
96  startColumnBlock_=NULL;
97  intsInBlock_=NULL;
98  fingerPrint_=NULL;
99  fullDWEverySoOften_=0;
100  numberPasses_=0;
101  maximumDW_=0;
102  numberDW_=0;
103  numberDWTimes_=0;
104  sizeFingerPrint_=0;
105  numberMasterColumns_=0;
106  numberMasterRows_=0;
107  numberBlocks_=0;
108  keepContinuous_=0;
109  phase_=0;
110  pass_ = 0;
111  nNeededBase_=200;
112  nNodesBase_=500;
113  nNeeded_=nNeededBase_;
114  nNodes_=nNodesBase_;
115  solveState_=0;
116}
117// Guts of copy
118void 
119CbcHeuristicDW::gutsOfCopy(const CbcHeuristicDW & rhs)
120{
121  targetObjective_ = rhs.targetObjective_;
122  bestObjective_ = rhs.bestObjective_;
123  lastObjective_ = rhs.lastObjective_;
124  fullDWEverySoOften_ = rhs.fullDWEverySoOften_;
125  numberPasses_ = rhs.numberPasses_;
126  howOften_ = rhs.howOften_;
127  decayFactor_ = rhs.decayFactor_;
128  fullDWEverySoOften_ = rhs.fullDWEverySoOften_;
129  numberPasses_ = rhs.numberPasses_;
130  maximumDW_ = rhs.maximumDW_;
131  numberDW_ = rhs.numberDW_;
132  numberDWTimes_ = rhs.numberDWTimes_;
133  sizeFingerPrint_ = rhs.sizeFingerPrint_;
134  numberMasterColumns_ = rhs.numberMasterColumns_;
135  numberMasterRows_ = rhs.numberMasterRows_;
136  numberBlocks_ = rhs.numberBlocks_;
137  keepContinuous_ = rhs.keepContinuous_;
138  phase_ = rhs.phase_;
139  pass_ = rhs.pass_;
140  nNeededBase_ = rhs.nNeededBase_;
141  nNodesBase_ = rhs.nNodesBase_;
142  nNeeded_ = rhs.nNeeded_;
143  nNodes_ = rhs.nNodes_;
144  solveState_ = rhs.solveState_;
145  functionPointer_ = rhs.functionPointer_;
146  if (rhs.solver_)
147    solver_ = rhs.solver_->clone();
148  else
149    solver_ = NULL;
150  if (rhs.dwSolver_)
151    dwSolver_ = rhs.dwSolver_->clone();
152  else
153    dwSolver_=NULL;
154  if (rhs.saveLower_) {
155    int numberColumns = solver_->getNumCols();
156    int numberRows = solver_->getNumRows();
157    saveLower_ = CoinCopyOfArray(rhs.saveLower_,numberColumns);
158    saveUpper_ = CoinCopyOfArray(rhs.saveUpper_,numberColumns);
159    whichColumnBlock_ = CoinCopyOfArray(rhs.whichColumnBlock_,numberColumns);
160    columnsInBlock_ = CoinCopyOfArray(rhs.columnsInBlock_,numberColumns);
161    whichRowBlock_ = CoinCopyOfArray(rhs.whichRowBlock_,numberRows);
162    rowsInBlock_ = CoinCopyOfArray(rhs.rowsInBlock_,numberRows);
163    if (rhs.affinity_)
164      affinity_ = CoinCopyOfArray(rhs.affinity_,numberBlocks_*numberBlocks_);
165    else
166      affinity_ = NULL;
167    backwardRow_ = CoinCopyOfArray(rhs.backwardRow_,numberRows);
168    startRowBlock_ = CoinCopyOfArray(rhs.startRowBlock_,numberBlocks_+1);
169    startColumnBlock_ = CoinCopyOfArray(rhs.startColumnBlock_,numberBlocks_+1);
170    intsInBlock_ = CoinCopyOfArray(rhs.intsInBlock_,numberBlocks_);
171  } else {
172    saveLower_=NULL;
173    saveUpper_=NULL;
174    affinity_=NULL;
175    whichRowBlock_=NULL;
176    whichColumnBlock_=NULL;
177    backwardRow_=NULL;
178    rowsInBlock_=NULL;
179    columnsInBlock_=NULL;
180    startRowBlock_=NULL;
181    startColumnBlock_=NULL;
182    intsInBlock_=NULL;
183  }
184  if (rhs.weights_) {
185    assert (maximumDW_);
186    weights_ = CoinCopyOfArray(rhs.weights_,maximumDW_);
187    random_ = CoinCopyOfArray(rhs.random_,numberMasterRows_);
188    dwBlock_ = CoinCopyOfArray(rhs.dwBlock_,maximumDW_);
189    fingerPrint_ = CoinCopyOfArray(rhs.fingerPrint_,
190                                   sizeFingerPrint_*maximumDW_);
191    objectiveDW_ = CoinCopyOfArray(rhs.objectiveDW_,numberDWTimes_);
192    numberColumnsDW_ = CoinCopyOfArray(rhs.numberColumnsDW_,numberDWTimes_);
193  } else {
194    random_=NULL;
195    weights_=NULL;
196    objectiveDW_=NULL;
197    numberColumnsDW_=NULL;
198    dwBlock_=NULL;
199    fingerPrint_=NULL;
200  }
201  if (rhs.bestSolution_) {
202    int numberColumns = solver_->getNumCols();
203    bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,numberColumns);
204  } else {
205    bestSolution_=NULL;
206  }
207  if (rhs.continuousSolution_) {
208    int numberColumns = solver_->getNumCols();
209    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
210  } else {
211    continuousSolution_=NULL;
212  }
213  if (rhs.fixedDj_) {
214    int numberColumns = solver_->getNumCols();
215    fixedDj_ = CoinCopyOfArray(rhs.fixedDj_,numberColumns);
216  } else {
217    fixedDj_=NULL;
218  }
219}
220// Guts of delete
221void 
222CbcHeuristicDW::gutsOfDelete()
223{
224  delete solver_;
225  delete dwSolver_;
226  delete [] bestSolution_;
227  delete [] continuousSolution_;
228  delete [] fixedDj_;
229  delete [] saveLower_;
230  delete [] saveUpper_;
231  delete [] random_;
232  delete [] affinity_;
233  delete [] weights_;
234  delete [] objectiveDW_;
235  delete [] numberColumnsDW_;
236  delete [] whichRowBlock_;
237  delete [] whichColumnBlock_;
238  delete [] dwBlock_;
239  delete [] backwardRow_;
240  delete [] rowsInBlock_;
241  delete [] columnsInBlock_;
242  delete [] startRowBlock_;
243  delete [] startColumnBlock_;
244  delete [] intsInBlock_;
245  delete [] fingerPrint_;
246  //functionPointer_ = NULL;
247  solver_ = NULL;
248  dwSolver_=NULL;
249  bestSolution_=NULL;
250  continuousSolution_ = NULL;
251  fixedDj_ = NULL;
252  saveLower_=NULL;
253  saveUpper_=NULL;
254  random_=NULL;
255  affinity_=NULL;
256  weights_=NULL;
257  objectiveDW_=NULL;
258  numberColumnsDW_ = NULL;
259  whichRowBlock_=NULL;
260  whichColumnBlock_=NULL;
261  dwBlock_=NULL;
262  backwardRow_=NULL;
263  rowsInBlock_=NULL;
264  columnsInBlock_=NULL;
265  startRowBlock_=NULL;
266  startColumnBlock_=NULL;
267  intsInBlock_=NULL;
268  fingerPrint_=NULL;
269  numberBlocks_=0;
270}
271
272// Destructor
273CbcHeuristicDW::~CbcHeuristicDW ()
274{
275  gutsOfDelete();
276}
277
278// Clone
279CbcHeuristic *
280CbcHeuristicDW::clone() const
281{
282  return new CbcHeuristicDW(*this);
283}
284
285// Assignment operator
286CbcHeuristicDW &
287CbcHeuristicDW::operator=( const CbcHeuristicDW & rhs)
288{
289  if (this != &rhs) {
290    CbcHeuristic::operator=(rhs);
291    gutsOfDelete();
292    gutsOfCopy(rhs);
293  }
294  return *this;
295}
296
297// Create C++ lines to get to current state
298void
299CbcHeuristicDW::generateCpp( FILE * fp)
300{
301  abort();
302}
303
304// Copy constructor
305CbcHeuristicDW::CbcHeuristicDW(const CbcHeuristicDW & rhs)
306  :
307  CbcHeuristic(rhs)
308{
309  gutsOfCopy(rhs);
310}
311// Resets stuff if model changes
312void
313CbcHeuristicDW::resetModel(CbcModel * model)
314{
315  if (model_&&numberBlocks_&&
316      model->getNumCols()!=model->getNumCols())
317    abort();
318  model_=model;
319}
320/*
321  First tries setting a variable to better value.  If feasible then
322  tries setting others.  If not feasible then tries swaps
323  Returns 1 if solution, 0 if not */
324int
325  CbcHeuristicDW::solution(double & solutionValue,
326                           double * betterSolution)
327{
328  numCouldRun_++;
329  int returnCode = 0;
330  const double * bestSolutionIn = model_->bestSolution();
331  if (!bestSolutionIn && !bestSolution_)
332    return 0; // No solution found yet
333  if (numberBlocks_<3)
334    return 0; // no point
335  if (bestSolutionIn&&objectiveValue(bestSolutionIn)<bestObjective_-1.0e-5)
336    passInSolution(bestSolutionIn);
337  printf("Before PASSes objective is %g\n",
338         bestObjective_);
339  double startTime = CoinCpuTime();
340  double startTimeElapsed = CoinGetTimeOfDay();
341  CoinWarmStart * basis = NULL;
342  lastObjective_ = COIN_DBL_MAX;
343  int passesToDW = dwSolver_ ? 0 : -1;
344  bool goodSolution=true;
345  int numberColumns = solver_->getNumCols();
346  // For moment just OsiClp
347  OsiClpSolverInterface * solver = dynamic_cast<OsiClpSolverInterface *>
348    (solver_);
349  ClpSimplex * simplex = solver->getModelPtr();
350  double * columnLower = simplex->columnLower();
351  double * columnUpper = simplex->columnUpper();
352  const double * cost = solver->getObjCoefficients();
353  const double * dj = solver->getReducedCost();
354  assert (solver);
355  if (!continuousSolution_) {
356    bool takeHint;
357    OsiHintStrength strength;
358    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
359    solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
360    solver->resolve();
361    solver->setHintParam(OsiDoDualInResolve, takeHint, OsiHintDo);
362    continuousSolution_ = CoinCopyOfArray(solver->getColSolution(),
363                                          numberColumns);
364  }
365  // data arrays
366  // Block order for choosing (after sort)
367  int * whichBlock = new int [7*numberBlocks_];
368  memset(whichBlock,0,7*numberBlocks_*sizeof(int));
369  // Count of number of times block chosen
370  int * doneBlock = whichBlock + numberBlocks_;
371  // Pass at which block last used
372  int * whenBlock = doneBlock+numberBlocks_;
373  // Number of Big Djs' (? artificial costs) in block
374  int * bigDjBlock = whenBlock+numberBlocks_;
375  // Number of times block has helped improve solution
376  int * goodBlock = bigDjBlock+numberBlocks_;
377  int * priorityBlock = goodBlock+numberBlocks_;
378  int * orderBlock = priorityBlock+numberBlocks_;
379  // Mixture of stuff to sort blocks on
380  double * blockSort = new double [4*numberBlocks_];
381  // Reduced cost (d sub j was old notation) contribution
382  double * blockDj = blockSort + numberBlocks_;
383  // Difference between current best and continuous solutions
384  double * blockDiff = blockDj+numberBlocks_;
385  // Difference between current best and continuous solutions (just integers)
386  double * blockDiffInt = blockDiff+numberBlocks_;
387  delete [] fixedDj_;
388  fixedDj_ = CoinCopyOfArray(dj,numberColumns);
389  int numberImproving=0;
390  int * whenBetter = new int [numberPasses_];
391  double * improvement = new double [numberPasses_];
392  // First has number
393  int ** improvingBlocks = new int * [numberPasses_];
394  int numberBlocksUsed = numberBlocks_;
395  // Get basic priority order
396  for (int i=0;i<numberBlocks_;i++) {
397    whenBlock[i]=-intsInBlock_[i];
398    orderBlock[i]=i;
399  }
400  CoinSort_2(whenBlock,whenBlock+numberBlocks_,orderBlock);
401  for (int i=0;i<numberBlocks_;i++) {
402    whenBlock[i]=-100;
403    doneBlock[i]=0;
404    whichBlock[i]=i;
405  }
406  bestObjective_=objectiveValue(bestSolution_);
407  double startTime2 = CoinCpuTime();
408  double startTime2Elapsed = CoinGetTimeOfDay();
409  // make sure all master columns freed up
410  for ( int iColumn=0 ; iColumn<numberColumns ; ++iColumn ) {
411    if (whichColumnBlock_[iColumn]<0) {
412      columnLower[iColumn]=saveLower_[iColumn];
413      columnUpper[iColumn]=saveUpper_[iColumn];
414    }
415  }
416  for (pass_=0;pass_<numberPasses_;pass_++) {
417    double endTime2 = CoinCpuTime();
418    double endTime2Elapsed = CoinGetTimeOfDay();
419#ifndef SCALE_FACTOR
420#define SCALE_FACTOR 1.0
421#endif
422    if (pass_)
423      printf("PASS %d changed objective from %g to %g in %g seconds (%g elapsed) - total %g (%g elapsed) - current needed %d nodes %d - %s\n",
424             pass_,lastObjective_*SCALE_FACTOR,
425             bestObjective_*SCALE_FACTOR,endTime2-startTime2,
426             endTime2Elapsed-startTime2Elapsed,
427             endTime2-startTime,endTime2Elapsed-startTimeElapsed,
428             nNeeded_,nNodes_,
429             lastObjective_>bestObjective_+1.0e-3 ? "improving" : "");
430    if ((pass_%10)==9) {
431      for (int iImp=CoinMax(1,numberImproving-10);iImp<numberImproving;iImp++) {
432        int * blocks = improvingBlocks[iImp];
433        int nBlocks = blocks[0];
434        blocks++;
435        printf("Pass %d improved objective by %g using %d blocks - ",
436               whenBetter[iImp],improvement[iImp],nBlocks);
437        for (int i=0;i<nBlocks;i++)
438          printf("%d ",blocks[i]);
439        printf("\n");
440      }
441      int * count = new int [numberImproving+1];
442      memset(count,0,(numberImproving+1)*sizeof(int));
443      for (int i=0;i<numberBlocks_;i++)
444        count[goodBlock[i]]++;
445      for (int i=0;i<numberImproving;i++) {
446        if (count[i])
447          printf("%d blocks were involved in improvement %d times\n",
448                 count[i],i);
449      }
450      delete [] count;
451    }
452    startTime2 = CoinCpuTime();
453    startTime2Elapsed = CoinGetTimeOfDay();
454    if (model_->getNodeCount()>=model_->getMaximumNodes()||
455        model_->maximumSecondsReached()) {
456      printf("Exiting on time or interrupt\n");
457      break;
458    }
459    if (bestObjective_>=lastObjective_-1.0e-3) {
460      // what now
461      // 0 - fine, 1 can't be better, 2 max node
462      //assert(solveState);
463      if (solveState_<2) {
464              // more in
465        printf("No improvement - think we need more variables ");
466        nNeeded_ += nNeeded_/10;
467        nNeeded_ = CoinMin(nNeeded_,800);
468        nNodes_=nNodesBase_;
469        (*(functionPointer_))(this,NULL,6);
470      } else {
471        // more nodes fewer in
472        printf("No improvement - stopped on nodes - think we need more nodes ");
473        if (phase_) {
474          nNodes_ += nNodes_/5;
475          nNodes_ = CoinMin(nNodes_,1000);
476        }
477        nNeeded_ -= nNeeded_/20;
478        nNeeded_ = CoinMax(nNeeded_,50);
479        (*(functionPointer_))(this,NULL,7);
480      }
481    } else {
482      // improving
483      (*(functionPointer_))(this,NULL,4);
484      solveState_=0;
485      //lastObjective=bestObjective_;
486      if (phase_) {
487        //nNeededBase_ += nNeededBase_/50;
488        //nNodesBase_ += nNodesBase_/50;
489      }
490      nNeeded_ -= nNeeded_/10;
491      nNeeded_=CoinMax(nNeededBase_,nNeeded_);
492      nNodes_=nNodesBase_;
493      (*(functionPointer_))(this,NULL,8);
494    }
495    printf("new needed %d, nodes %d\n",nNeeded_,nNodes_); 
496    for ( int i=0 ; i<numberColumns ; ++i ) {
497      if (solver_->isInteger(i)) {
498        double value = floor(bestSolution_[i]+0.5);
499        columnLower[i] = value;
500        columnUpper[i] = value;
501      } else {
502        columnLower[i] = saveLower_[i];
503        columnUpper[i] = saveUpper_[i];
504      }
505    }
506    if (goodSolution) {
507      int lastNumberDW=numberDW_;
508      solver->setColSolution(bestSolution_);
509      solver->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
510      if (basis) {
511        solver->setWarmStart(basis);
512        delete basis;
513        basis=NULL;
514      }
515      solver->resolve();
516      solver->setHintParam(OsiDoReducePrint, true, OsiHintDo, 0) ;
517      if (solver->getObjValue()<bestObjective_-1.0e-5) {
518        memcpy(bestSolution_,solver->getColSolution(),
519               numberColumns*sizeof(double));
520        bestObjective_ = solver->getObjValue();
521        int * blocks = new int [numberBlocksUsed+1];
522        blocks[0]=numberBlocksUsed;
523        memcpy(blocks+1,whichBlock,numberBlocksUsed*sizeof(int));
524        improvingBlocks[numberImproving]=blocks;
525        whenBetter[numberImproving]=pass_;
526        improvement[numberImproving]=lastObjective_-bestObjective_;
527        numberImproving++;
528        lastObjective_=bestObjective_;
529        if (pass_) {
530          // update good
531          for (int i=0;i<numberBlocksUsed;i++)
532            goodBlock[whichBlock[i]]++;
533        }
534      }
535      goodSolution=false;
536      if (fullDWEverySoOften_>0) {
537        addDW(bestSolution_,numberBlocksUsed,
538                         whichBlock);
539      }
540      if (passesToDW==0) {
541        passesToDW = fullDWEverySoOften_;
542        const double * duals = solver->getRowPrice();
543        double * bestSolution2 = CoinCopyOfArray(bestSolution_,
544                                                 numberColumns);
545        // Column copy
546        const double * element = solver->getMatrixByCol()->getElements();
547        const int * row = solver->getMatrixByCol()->getIndices();
548        const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
549        const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
550        for (int i=0;i<numberBlocks_;i++)
551          whichBlock[i]=i;
552        for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
553          int start=startColumnBlock_[iBlock];
554          int end=startColumnBlock_[iBlock+1];
555          ClpSimplex * tempModel = new
556            ClpSimplex(solver->getModelPtr(),
557                       startRowBlock_[iBlock+1]-startRowBlock_[iBlock],
558                       rowsInBlock_+startRowBlock_[iBlock],
559                       end-start,
560                       columnsInBlock_+startColumnBlock_[iBlock]);
561          tempModel->setLogLevel(0);
562          double * objectiveX = tempModel->objective();
563          double * columnLowerX = tempModel->columnLower();
564          double * columnUpperX = tempModel->columnUpper();
565          for (int i=start;i<end;i++) {
566            int jColumn=i-start;
567            int iColumn=columnsInBlock_[i];
568            columnLowerX[jColumn]=saveLower_[iColumn];
569            columnUpperX[jColumn]=saveUpper_[iColumn];
570            if (solver->isInteger(iColumn))
571              tempModel->setInteger(jColumn);
572            double cost=objectiveX[jColumn];
573            for (CoinBigIndex j=columnStart[iColumn];
574                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
575              int iRow = row[j];
576              double elementValue=element[j];
577              if (backwardRow_[iRow]>=0) {
578                cost -= elementValue * duals[iRow]; 
579              }
580            }
581            objectiveX[jColumn]=cost;
582          }
583          OsiClpSolverInterface solverX(tempModel,true);
584          CbcModel modelX(solverX);
585          modelX.setLogLevel(0);
586          modelX.branchAndBound();
587          const double * bestSolutionX = modelX.bestSolution();
588          assert (bestSolutionX);
589          for (int i=start;i<end;i++) {
590            int jColumn = i-start;
591            int iColumn=columnsInBlock_[i];
592            bestSolution2[iColumn]=bestSolutionX[jColumn];
593          }
594        }
595        addDW(bestSolution2,numberBlocks_,whichBlock);
596        // now try purer DW
597        bool takeHint;
598        OsiHintStrength strength;
599        dwSolver_->getHintParam(OsiDoDualInResolve, takeHint, strength);
600        dwSolver_->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
601        dwSolver_->resolve();
602        dwSolver_->setHintParam(OsiDoDualInResolve, takeHint, OsiHintDo);
603        duals = dwSolver_->getRowPrice();
604        for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
605          int start=startColumnBlock_[iBlock];
606          int end=startColumnBlock_[iBlock+1];
607          ClpSimplex * tempModel = new
608            ClpSimplex(solver->getModelPtr(),
609                       startRowBlock_[iBlock+1]-startRowBlock_[iBlock],
610                       rowsInBlock_+startRowBlock_[iBlock],
611                       end-start,
612                       columnsInBlock_+startColumnBlock_[iBlock]);
613          tempModel->setLogLevel(0);
614          double * objectiveX = tempModel->objective();
615          double * columnLowerX = tempModel->columnLower();
616          double * columnUpperX = tempModel->columnUpper();
617          double convexityDual = duals[numberMasterRows_+iBlock];
618          for (int i=start;i<end;i++) {
619            int jColumn=i-start;
620            int iColumn=columnsInBlock_[i];
621            columnLowerX[jColumn]=saveLower_[iColumn];
622            columnUpperX[jColumn]=CoinMin(saveUpper_[iColumn],1.0e-12);
623            if (solver->isInteger(iColumn))
624              tempModel->setInteger(jColumn);
625            double cost=objectiveX[jColumn];
626            for (CoinBigIndex j=columnStart[iColumn];
627                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
628              int iRow = row[j];
629              double elementValue=element[j];
630              if (backwardRow_[iRow]>=0) {
631                // duals are from dw
632                cost -= elementValue * duals[backwardRow_[iRow]]; 
633              }
634            }
635            objectiveX[jColumn]=cost;
636          }
637          OsiClpSolverInterface solverX(tempModel,true);
638          solverX.initialSolve();
639          double cObj=solverX.getObjValue();
640          CbcModel modelX(solverX);
641          modelX.setLogLevel(0);
642          modelX.branchAndBound();
643          printf("Block %d contobj %g intobj %g convdual %g\n",
644                 iBlock,cObj,modelX.getObjValue(),convexityDual);
645          const double * bestSolutionX = modelX.bestSolution();
646          assert (bestSolutionX);
647          for (int i=start;i<end;i++) {
648            int iColumn=columnsInBlock_[i];
649            bestSolution2[iColumn]=bestSolutionX[i-start];
650          }
651        }
652        addDW(bestSolution2,numberBlocks_,whichBlock);
653        delete [] bestSolution2;
654      }
655      passesToDW--;
656      if (numberDW_>lastNumberDW) {
657        intArray_=NULL;
658        doubleArray_=NULL;
659        (*(functionPointer_))(this,NULL,3);
660      }
661    }
662    for (int i=0;i<numberBlocks_;i++) {
663      blockDj[i]=0.0;
664      blockDiff[i]=0.0;
665      blockDiffInt[i]=0.0;
666      bigDjBlock[i]=0;
667    }
668    for ( int i=0 ; i<numberColumns ; ++i ) {
669      int kBlock = whichColumnBlock_[i];
670      if (kBlock>=0) {
671        columnLower[i]=bestSolution_[i];
672        columnUpper[i]=bestSolution_[i];
673      }
674    }
675    for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
676      for (int i=startColumnBlock_[iBlock];
677           i<startColumnBlock_[iBlock+1];i++) {
678        int iColumn=columnsInBlock_[i];
679        if (continuousSolution_) {
680          blockDiff[iBlock] += 
681            fabs((bestSolution_[iColumn]-continuousSolution_[iColumn]))
682            *(fabs(cost[iColumn])+1.0e-5);
683          if (solver->isInteger(iColumn)) 
684            blockDiffInt[iBlock] += 
685              fabs((bestSolution_[iColumn]-continuousSolution_[iColumn]))
686              *(fabs(cost[iColumn])+1.0e-5);
687        }
688        if (solver->isInteger(iColumn)) {
689          if (bestSolution_[iColumn]<saveUpper_[iColumn]-1.0e-1) {
690            if (fixedDj_[iColumn]<-1.0e-5)
691              blockDj[iBlock]+=fixedDj_[iColumn];
692            if (fixedDj_[iColumn]<-1.0e4)
693              bigDjBlock[iBlock]++;
694          }
695          if (bestSolution_[iColumn]>saveLower_[iColumn]+1.0e-1) {
696            if (fixedDj_[iColumn]>1.0e-5)
697              blockDj[iBlock]-=fixedDj_[iColumn];
698            if (fixedDj_[iColumn]>1.0e4)
699              bigDjBlock[iBlock]++;
700          }
701        }
702      }
703    }
704    // Get average dj and difference
705    int numberInDj=0;
706    double averageDj=1.0e-12;
707    int numberInDiff=0;
708    double averageDiff=1.0e-12;
709    for (int i=0;i<numberBlocks_;i++) {
710      assert (blockDj[i]<=0.0);
711      if (blockDj[i]<0.0) {
712        numberInDj++;
713        averageDj -= blockDj[i];
714      }
715      assert (blockDiff[i]>=0.0);
716      if (blockDiff[i]>0.0) {
717        numberInDiff++;
718        averageDiff += blockDiff[i];
719      }
720    }
721    if (numberInDj)
722      averageDj /= static_cast<double>(numberInDj);
723    if (numberInDiff)
724      averageDiff /= static_cast<double>(numberInDiff);
725    double ratioDiff = averageDj/averageDiff;
726    // downplay
727    ratioDiff *= 1.0e-3;
728    for (int i=0;i<numberBlocks_;i++) {
729      whichBlock[i]=i;
730      assert (intsInBlock_[i]);
731      blockSort[i]=blockDj[i];
732#define TRY_ADJUST 3
733#if TRY_ADJUST == 1
734      blockSort[i] *= CoinDrand48() + 5.0e-1;
735#elif TRY_ADJUST == 2
736      blockSort[i] -= ratioDiff* blockDiff[i];
737#elif TRY_ADJUST == 3
738      blockSort[i] -= ratioDiff* blockDiff[i];
739      blockSort[i] *= 0.1*CoinDrand48() + 0.9;
740#endif
741      if (phase_==99)
742        blockSort[i] -= 2.0*averageDj*goodBlock[i];
743      blockSort[i] /= static_cast<double>(intsInBlock_[i]);
744      //blockSort[i] /= sqrt(static_cast<double>(intsInBlock_[i]));
745      if (doneBlock[i]) {
746        blockSort[i] /= static_cast<double>(doneBlock[i]+1);
747        if (whenBlock[i]>pass_-10)
748          blockSort[i]+= 1.0e2*averageDj;
749      }
750    }
751    CoinSort_2(blockSort,blockSort+numberBlocks_,whichBlock);
752    // allow user to modify
753    intArray_=whichBlock;
754    doubleArray_=blockSort;
755    (*(functionPointer_))(this,NULL,1);
756    int numberBlocksIn=0;
757    for (int iTry=0;iTry<2;iTry++) {
758      int nFreed=0;
759      int nBigDjBlock=0;
760      numberBlocksUsed=0;
761      for (int i=0;i<numberBlocks_;i++) {
762        int iBlock = whichBlock[i];
763        bool skipBlock=false;
764        assert (iBlock>=0);
765        if ((doneBlock[iBlock]&&!phase_)||!blockSort[i]) {
766          //printf("already done block %d - dj %g\n",iBlock,blockDj[i]);
767          skipBlock=true;
768        } else if (bigDjBlock[iBlock]) {
769          nBigDjBlock++;
770          if (nBigDjBlock>20&&!phase_) {
771            skipBlock=true;
772          }
773        }
774        int returnCode = (*(functionPointer_))(this,NULL,5);
775        if (returnCode<0)
776          skipBlock=true;
777        else if (returnCode>0)
778          skipBlock=false;
779        if (skipBlock) {
780          whichBlock[i] -= 1000000;
781        } else {
782          // free up
783          printf("freeing block %d (already freed %d times) - dj %g, diff %g, diffint %g - %d columns (%d integer)\n",
784                 iBlock,doneBlock[iBlock],blockDj[iBlock],
785                 blockDiff[iBlock],blockDiffInt[iBlock],
786                 startColumnBlock_[iBlock+1]-startColumnBlock_[iBlock],
787                 intsInBlock_[iBlock]);
788          numberBlocksIn++;
789          doneBlock[iBlock]++;
790          whenBlock[iBlock]=pass_;
791          nFreed += intsInBlock_[iBlock];
792          for (int j=startColumnBlock_[iBlock];
793           j<startColumnBlock_[iBlock+1];j++) {
794            int iColumn=columnsInBlock_[j];
795            columnLower[iColumn]=saveLower_[iColumn];
796            columnUpper[iColumn]=saveUpper_[iColumn];
797          }
798          if (!numberBlocksUsed) {
799            // re-sort rest using affinity
800            const unsigned short int * aff = affinity_+iBlock*numberBlocks_;
801            for (int j=i+1;j<numberBlocks_;j++) {
802              int jBlock=whichBlock[j];
803              blockSort[j] -= 1.0e2*aff[jBlock];
804            }
805            CoinSort_2(blockSort+i+1,blockSort+numberBlocks_,whichBlock+i+1);
806          }
807          numberBlocksUsed=i+1;
808          if (nFreed>=nNeeded_ && numberBlocksIn>3)
809            break;
810        }
811      }
812      printf("%d big dj blocks found\n",nBigDjBlock);
813      if (nFreed)
814        break;
815      phase_=1; // round again
816      printf("Changing phase\n");
817      for (int i=0;i<numberBlocks_;i++) 
818        whichBlock[i] += 1000000;
819      //nNeeded=500; // allow more
820    }
821    if (numberBlocksIn==numberBlocks_)
822      numberPasses_=1; // no point
823    // pack down blocks
824    int n=0;
825    for (int i=0;i<numberBlocksUsed;i++) {
826      if (whichBlock[i]>=0)
827        whichBlock[n++]=whichBlock[i];
828    }
829    numberBlocksUsed = n;
830    int nFixedInts=0;
831    int nFixedContinuous=0;
832    int nFreeInts=0;
833    int nFreeContinuous=0;
834    int nFreeMaster=0;
835    int nFixedMaster=0;
836    for ( int i=0 ; i<numberColumns ; ++i ) {
837      int kBlock = whichColumnBlock_[i];
838      if (columnUpper[i]>columnLower[i]) {
839        if (kBlock>=0) {
840          if (solver->isInteger(i))
841            nFreeInts++;
842          else
843            nFreeContinuous++;
844        } else {
845          nFreeMaster++;
846        }
847      } else {
848        if (kBlock>=0) {
849          if (solver->isInteger(i))
850            nFixedInts++;
851          else
852            nFixedContinuous++;
853        } else {
854          nFixedMaster++;
855        }
856      }
857    }
858    printf("Fixed %d ints, %d c, %d m - free %d, %d, %d\n",
859           nFixedInts,nFixedContinuous,nFixedMaster,
860           nFreeInts,nFreeContinuous,nFreeMaster);
861    // But free up and then fix again
862    for ( int i=0 ; i<numberColumns ; ++i ) {
863      columnLower[i]=saveLower_[i];
864      columnUpper[i]=saveUpper_[i];
865      int kBlock = whichColumnBlock_[i];
866      if (kBlock>=0&&solver->isInteger(i)&&whenBlock[kBlock]!=pass_) {
867        columnLower[i]=bestSolution_[i];
868        columnUpper[i]=bestSolution_[i];
869      }
870    }
871    bool takeHint;
872    OsiHintStrength strength;
873    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
874    solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
875    solver->messageHandler()->setLogLevel(1) ;
876    solver->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
877    solver->resolve();
878    solver->setHintParam(OsiDoReducePrint, true, OsiHintDo, 0) ;
879    //solver->messageHandler()->setLogLevel(0) ;
880    solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
881    if (solver->getObjValue()>bestObjective_+1.0e-5*(1.0+fabs(bestObjective_))) {
882      // trouble
883      for (int i=0;i<numberBlocks_;i++) {
884        if (whenBlock[i]==pass_) {
885          printf("Block %d free\n",i);
886        }
887      }
888      solver->writeMps("bad","mps");
889      const double * lower = solver->getColLower();
890      const double * upper = solver->getColUpper();
891      printf("best obj %g\n",objectiveValue(bestSolution_));
892      for ( int i=0 ; i<numberColumns ; ++i ) {
893        double value = bestSolution_[i];
894        if (value<lower[i]-1.0e-5||value>upper[i]+1.0e-5)
895          printf("column %d (block %d) %g %g <= %g <= %g %g\n",
896                 i,whichColumnBlock_[i],saveLower_[i],
897                 lower[i],value,upper[i],saveUpper_[i]);
898      }
899      abort();
900    }
901    const double * tempSol = solver->getColSolution();
902    // But free up and then fix again
903    for ( int i=0 ; i<numberColumns ; ++i ) {
904      int kBlock = whichColumnBlock_[i];
905      if (kBlock>=0&&!solver->isInteger(i)&&whenBlock[kBlock]!=pass_) {
906        columnLower[i]=tempSol[i];
907        columnUpper[i]=tempSol[i];
908      }
909    }
910    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
911    solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
912    solver->messageHandler()->setLogLevel(1) ;
913    solver->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
914    solver->resolve();
915    solver->setHintParam(OsiDoReducePrint, true, OsiHintDo, 0) ;
916    //solver->messageHandler()->setLogLevel(0) ;
917    solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
918    //solver->messageHandler()->setLogLevel(0) ;
919    //lp->setLogLevel(0);
920    ClpPresolve pinfo;
921    // fix small infeasibilities
922    pinfo.setPresolveActions(pinfo.presolveActions()|0x4000);
923    int numberPasses = 2; // can change this
924    ClpSimplex * model2 = pinfo.presolvedModel(*simplex, 1.0e-8,
925                                               true, 
926                                               numberPasses, true);
927    if (!model2) {
928      abort();
929    } else {
930      printf("Reduced model has %d rows and %d columns\n",
931             model2->numberRows(),model2->numberColumns());
932      //model2->setLogLevel(0);
933      OsiClpSolverInterface solver2(model2);
934      solver2.setWarmStart(NULL);
935      //solver2.messageHandler()->setLogLevel(0) ;
936      CbcModel model(solver2);
937      model.setMaximumNodes(nNodes_);
938      //model.setMaximumSolutions(2);
939      model.solver()->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
940      model.solver()->setHintParam(OsiDoPresolveInInitial, false, OsiHintDo, 0) ;
941      //model.solver()->setHintParam(OsiDoPresolveInResolve, false, OsiHintDo, 0) ;
942      model.initialSolve();
943      if (bestObjective_ > model.solver()->getObjValue()+1.0e-1) {
944        const double * dj = model.solver()->getReducedCost();
945        int nFix=0;
946        const double * lower = model.solver()->getColLower();
947        const double * upper = model.solver()->getColUpper();
948        const double * solution = model.solver()->getColSolution();
949        double gap = CoinMax(bestObjective_-model.solver()->getObjValue(),
950                                   1.0e-3);
951        int numberColumns2=model.solver()->getNumCols();
952#ifdef HOT_START
953        // Set up hot start
954        const int * originalColumns = pinfo.originalColumns();
955        double * hot = new double[2*numberColumns2];
956        int * hotPriorities = new int[numberColumns2*2];
957        double * hotWeight = hot+numberColumns2; 
958        memset(hot,0,numberColumns2*sizeof(double));
959        int * sort = hotPriorities+numberColumns2;
960        for (int i=0;i<numberColumns2;i++) {
961          hotWeight[i]=COIN_DBL_MAX;
962          sort[i]=i;
963          if (solver2.isInteger(i)) {
964            int iColumn=originalColumns[i];
965            hot[i]=bestSolution_[iColumn];
966            hotWeight[i]=-fabs(fixedDj_[iColumn]);
967            if (bestSolution_[i]<saveLower_[iColumn]+1.0e-6) {
968              if (fixedDj_[i]>0.0) 
969                hotWeight[i]=fixedDj_[iColumn];
970            } else if (bestSolution_[i]>saveUpper_[iColumn]-1.0e-6) {
971              if (fixedDj_[i]<0.0) 
972                hotWeight[i]=-fixedDj_[iColumn];
973            }
974            if (solution[i]<saveLower_[iColumn]+1.0e-6) {
975              if (dj[i]>gap) {
976                solver2.setColUpper(i,saveLower_[iColumn]);
977                nFix++;
978              } 
979            } else if (solution[i]>saveUpper_[iColumn]-1.0e-6) {
980              if (-dj[i]>gap) {
981                solver2.setColLower(i,saveUpper_[iColumn]);
982                nFix++;
983              }
984            }
985          }
986        }
987        CoinSort_2(hotWeight,hotWeight+numberColumns2,sort);
988        for (int i=0;i<numberColumns2;i++) {
989          hotPriorities[sort[i]]=i+1;
990        }
991        model.setHotstartSolution(hot,hotPriorities);
992        delete [] hot;
993        delete [] hotPriorities;
994#endif
995        if (nFix)
996          printf("Fixed another %d integers\n",nFix);
997        {
998          // priorities
999          memset(priorityBlock,0,numberBlocks_*sizeof(int));
1000          for (int i=0;i<numberBlocks_;i++) {
1001            int iBlock=whichBlock[i];
1002            if (iBlock>=0)
1003              priorityBlock[iBlock]=i+1;
1004          }
1005#if 1
1006          assert (numberBlocks_<4000);
1007          // But make sure we do one block before next
1008          for (int i=0;i<numberBlocks_;i++) {
1009            priorityBlock[i] += 4000*orderBlock[i];
1010          }
1011#else
1012          // just do ones with many first
1013          for (int i=0;i<numberBlocks;i++) {
1014            priorityBlock[i]=10000-intsInBlock_[i];
1015          }
1016#endif
1017          int numberColumns2=model.solver()->getNumCols();
1018          const int * original = pinfo.originalColumns();
1019          int * priorities = new int [numberColumns2];
1020          int n=0;
1021          for (int i=0;i<numberColumns2;i++) {
1022            if (solver2.isInteger(i)) {
1023              int iBlock = whichColumnBlock_[original[i]];
1024              priorities[n++]=priorityBlock[iBlock];
1025            }
1026          }
1027          model.passInPriorities(priorities,false);
1028          delete [] priorities;
1029        }
1030        CglProbing probingGen;
1031        probingGen.setUsingObjective(1);
1032        probingGen.setMaxPass(1);
1033        probingGen.setMaxPassRoot(1);
1034        // Number of unsatisfied variables to look at
1035        probingGen.setMaxProbe(10);
1036        probingGen.setMaxProbeRoot(50);
1037        // How far to follow the consequences
1038        probingGen.setMaxLook(10);
1039        probingGen.setMaxLookRoot(50);
1040        probingGen.setMaxLookRoot(10);
1041        // Only look at rows with fewer than this number of elements
1042        probingGen.setMaxElements(200);
1043        probingGen.setMaxElementsRoot(300);
1044        probingGen.setRowCuts(0);
1045        model.addCutGenerator(&probingGen, 5, "Probing");
1046        model.setNumberThreads(model_->getNumberThreads());
1047        model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintDo, 0) ;
1048        model.setPrintingMode(1);
1049        model.setCutoff(lastObjective_-model_->getCutoffIncrement());
1050        model.setLogLevel(1);
1051        intArray_=NULL;
1052        doubleArray_=NULL;
1053        (*(functionPointer_))(this,&model,2);
1054        model.branchAndBound();
1055        printf("After B&B status %d objective %g\n",model.status(),
1056               model.getMinimizationObjValue());
1057        int modelStatus = model.status();
1058        model.solver()->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
1059        if (model.bestSolution() && 
1060            model.getMinimizationObjValue()<lastObjective_) {
1061          const int * original = pinfo.originalColumns();
1062          int n2=solver2.getNumCols();
1063          const double * bestSolution2 = model.bestSolution();
1064          for (int i=0;i<n2;i++) {
1065            int iColumn = original[i];
1066            if (solver2.isInteger(i)) {
1067              double value =floor(bestSolution2[i]+0.5);
1068              if (fabs(bestSolution2[i]-value)>1.0e-5) {
1069                printf("bad %d %g\n",i,bestSolution2[i]);
1070              } else {
1071                solver->setColLower(iColumn,value);
1072                solver->setColUpper(iColumn,value);
1073              }
1074            }
1075          }
1076          pinfo.postsolve(true);
1077          delete model2;
1078          bool takeHint;
1079          OsiHintStrength strength;
1080          solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1081          solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
1082          solver->messageHandler()->setLogLevel(1) ;
1083          solver->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
1084          solver->resolve();
1085          if (basis) 
1086            delete basis;
1087          basis = solver->getWarmStart();
1088          memcpy(fixedDj_,solver->getReducedCost(),
1089                 numberColumns*sizeof(double));
1090          solver->setHintParam(OsiDoReducePrint, true, OsiHintDo, 0) ;
1091          //solver->messageHandler()->setLogLevel(0) ;
1092          solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
1093          //solver->
1094          CbcModel model(*solver);
1095          model.setNumberBeforeTrust(-1);
1096          model.setNumberStrong(0);
1097          model.setMoreSpecialOptions2(57);
1098          model.branchAndBound();
1099          if (model.getMinimizationObjValue()<lastObjective_) {
1100            memcpy(bestSolution_,model.bestSolution(),
1101                   numberColumns*sizeof(double));
1102            bestObjective_ = model.getObjValue();
1103            for (int i=0;i<numberColumns;i++) {
1104              if (simplex->isInteger(i)) {
1105                if (fabs(bestSolution_[i]-floor(bestSolution_[i]+0.5))>1.0e-5) {
1106                  printf("bad after %d %g\n",i,bestSolution_[i]);
1107                }
1108              }
1109            }
1110            goodSolution=true;
1111          }
1112        } else if (modelStatus!=0) {
1113          // stopped on nodes
1114          // 0 - fine, 1 can't be better, 2 max node
1115          solveState_=2;
1116          pinfo.destroyPresolve();
1117          delete model2;
1118        } else {
1119          // complete search
1120          solveState_=1;
1121          pinfo.destroyPresolve();
1122          delete model2;
1123        }
1124      } else {
1125        // can't be better
1126        // 0 - fine, 1 can't be better, 2 max node
1127        solveState_=1;
1128        pinfo.destroyPresolve();
1129        delete model2;
1130      }
1131    }
1132  }
1133  delete [] whichBlock;
1134  delete [] blockSort;
1135  delete [] whenBetter;
1136  delete [] improvement;
1137  for (int i=0;i<numberImproving;i++)
1138    delete [] improvingBlocks[i];
1139  delete [] improvingBlocks;
1140  delete basis;
1141  return returnCode;
1142}
1143// update model
1144void 
1145CbcHeuristicDW::setModel(CbcModel * model)
1146{
1147  if (model!=model_) {
1148    gutsOfDelete();
1149    model_ = model;
1150    assert(model->solver());
1151    solver_ = model->solver()->clone();
1152    findStructure();
1153  }
1154}
1155// Find structure
1156void 
1157CbcHeuristicDW::findStructure()
1158{
1159  int numberRows = solver_->getNumRows();
1160  int numberColumns = solver_->getNumCols();
1161  // look for DW
1162  int * blockStart = new int [3*(numberRows+numberColumns)+1];
1163  int * columnBlock = blockStart+numberRows;
1164  int * nextColumn = columnBlock+numberColumns;
1165  int * blockCount = nextColumn+numberColumns;
1166  int * blockEls = blockCount+numberRows+1;
1167  int * countIntegers = blockEls+numberRows;
1168  memset(countIntegers,0,numberColumns*sizeof(int));
1169  int direction[2]={-1,1};
1170  int bestBreak=-1;
1171  double bestValue=0.0;
1172  int iPass=0;
1173  int halfway=(numberRows+1)/2;
1174  int firstMaster=-1;
1175  int lastMaster=-2;
1176  // Column copy
1177  const CoinPackedMatrix * columnCopy = solver_->getMatrixByCol();
1178  //const double * element = columnCopy->getElements();
1179  const int * row = columnCopy->getIndices();
1180  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
1181  const int * columnLength = columnCopy->getVectorLengths();
1182  // Row copy
1183  const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
1184  const int * column = rowCopy->getIndices();
1185  const int * rowLength = rowCopy->getVectorLengths();
1186  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1187  //const double * elementByRow = rowCopy->getElements();
1188  while (iPass<2) {
1189    int increment=direction[iPass];
1190    int start= increment>0 ? 0 : numberRows-1;
1191    int stop=increment>0 ? numberRows : -1;
1192    int numberBlocks=0;
1193    int thisBestBreak=-1;
1194    double thisBestValue=COIN_DBL_MAX;
1195    int numberRowsDone=0;
1196    int numberMarkedColumns=0;
1197    int maximumBlockSize=0;
1198    for (int i=0;i<numberRows+2*numberColumns;i++) 
1199      blockStart[i]=-1;
1200    for (int i=0;i<numberRows+1;i++)
1201      blockCount[i]=0;
1202    for (int iRow=start;iRow!=stop;iRow+=increment) {
1203      int iBlock = -1;
1204      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
1205        int iColumn=column[j];
1206        int whichColumnBlock=columnBlock[iColumn];
1207        if (whichColumnBlock>=0) {
1208          // column marked
1209          if (iBlock<0) {
1210            // put row in that block
1211            iBlock=whichColumnBlock;
1212          } else if (iBlock!=whichColumnBlock) {
1213            // merge
1214            blockCount[iBlock]+=blockCount[whichColumnBlock];
1215            blockCount[whichColumnBlock]=0;
1216            int jColumn=blockStart[whichColumnBlock];
1217            while (jColumn>=0) {
1218              columnBlock[jColumn]=iBlock;
1219              iColumn=jColumn;
1220              jColumn=nextColumn[jColumn];
1221            }
1222            nextColumn[iColumn]=blockStart[iBlock];
1223            blockStart[iBlock]=blockStart[whichColumnBlock];
1224            blockStart[whichColumnBlock]=-1;
1225          }
1226        }
1227      }
1228      int n=numberMarkedColumns;
1229      if (iBlock<0) {
1230        //new block
1231        if (rowLength[iRow]) {
1232          numberBlocks++;
1233          iBlock=numberBlocks;
1234          int jColumn=column[rowStart[iRow]];
1235          columnBlock[jColumn]=iBlock;
1236          blockStart[iBlock]=jColumn;
1237          numberMarkedColumns++;
1238          for (CoinBigIndex j=rowStart[iRow]+1;j<rowStart[iRow]+rowLength[iRow];j++) {
1239            int iColumn=column[j];
1240            columnBlock[iColumn]=iBlock;
1241            numberMarkedColumns++;
1242            nextColumn[jColumn]=iColumn;
1243            jColumn=iColumn;
1244          }
1245          blockCount[iBlock]=numberMarkedColumns-n;
1246        } else {
1247          // empty
1248          iBlock=numberRows;
1249        }
1250      } else {
1251        // put in existing block
1252        int jColumn=blockStart[iBlock];
1253        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
1254          int iColumn=column[j];
1255          assert (columnBlock[iColumn]<0||columnBlock[iColumn]==iBlock);
1256          if (columnBlock[iColumn]<0) {
1257            columnBlock[iColumn]=iBlock;
1258            numberMarkedColumns++;
1259            nextColumn[iColumn]=jColumn;
1260            jColumn=iColumn;
1261          }
1262        }
1263        blockStart[iBlock]=jColumn;
1264        blockCount[iBlock]+=numberMarkedColumns-n;
1265      }
1266      maximumBlockSize=CoinMax(maximumBlockSize,blockCount[iBlock]);
1267      numberRowsDone++;
1268      if (thisBestValue*numberRowsDone > maximumBlockSize&&numberRowsDone>halfway) { 
1269        thisBestBreak=iRow;
1270        thisBestValue=static_cast<double>(maximumBlockSize)/
1271          static_cast<double>(numberRowsDone);
1272      }
1273    }
1274    if (thisBestBreak==stop)
1275      thisBestValue=COIN_DBL_MAX;
1276    iPass++;
1277    if (iPass==1) {
1278      bestBreak=thisBestBreak;
1279      bestValue=thisBestValue;
1280    } else {
1281      if (bestValue<thisBestValue) {
1282        firstMaster=0;
1283        lastMaster=bestBreak;
1284      } else {
1285        firstMaster=thisBestBreak; // ? +1
1286        lastMaster=numberRows;
1287      }
1288    }
1289  }
1290  numberBlocks_=0;
1291  if (firstMaster<lastMaster) {
1292    printf("%d master rows %d <= < %d\n",lastMaster-firstMaster,
1293           firstMaster,lastMaster);
1294    for (int i=0;i<numberRows+2*numberColumns;i++) 
1295      blockStart[i]=-1;
1296    for (int i=firstMaster;i<lastMaster;i++)
1297      blockStart[i]=-2;
1298    int firstRow=0;
1299    int numberBlocks=-1;
1300    while (true) {
1301      for (;firstRow<numberRows;firstRow++) {
1302        if (blockStart[firstRow]==-1)
1303          break;
1304      }
1305      if (firstRow==numberRows)
1306        break;
1307      int nRows=0;
1308      numberBlocks++;
1309      int numberStack=1;
1310      blockCount[0] = firstRow;
1311      while (numberStack) {
1312        int iRow=blockCount[--numberStack];
1313        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
1314          int iColumn=column[j];
1315          int iBlock=columnBlock[iColumn];
1316          if (iBlock<0) {
1317            columnBlock[iColumn]=numberBlocks;
1318            for (CoinBigIndex k=columnStart[iColumn];
1319                 k<columnStart[iColumn]+columnLength[iColumn];k++) {
1320              int jRow=row[k];
1321              int rowBlock=blockStart[jRow];
1322              if (rowBlock==-1) {
1323                nRows++;
1324                blockStart[jRow]=numberBlocks;
1325                blockCount[numberStack++]=jRow;
1326              }
1327            }
1328          }
1329        }
1330      }
1331      if (!nRows) {
1332        // empty!!
1333        numberBlocks--;
1334      }
1335      firstRow++;
1336    }
1337    // adjust
1338    numberBlocks++;
1339    for (int i=0;i<numberBlocks;i++) {
1340      blockCount[i]=0;
1341      nextColumn[i]=0;
1342    }
1343    int numberEmpty=0;
1344    int numberMaster=0;
1345    memset(blockEls,0,numberBlocks*sizeof(int));
1346    for (int iRow = 0; iRow < numberRows; iRow++) {
1347      int iBlock=blockStart[iRow];
1348      if (iBlock>=0) {
1349        blockCount[iBlock]++;
1350        blockEls[iBlock]+=rowLength[iRow];
1351      } else {
1352        if (iBlock==-2)
1353          numberMaster++;
1354        else
1355          numberEmpty++;
1356      }
1357    }
1358    int numberEmptyColumns=0;
1359    int numberMasterColumns=0;
1360    int numberMasterIntegers=0;
1361    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1362      int iBlock=columnBlock[iColumn];
1363      bool isInteger = (solver_->isInteger(iColumn));
1364      if (iBlock>=0) {
1365        nextColumn[iBlock]++;
1366        if (isInteger)
1367          countIntegers[iBlock]++;
1368      } else {
1369        if (isInteger)
1370          numberMasterIntegers++;
1371        if (columnLength[iColumn])
1372          numberMasterColumns++;
1373        else
1374          numberEmptyColumns++;
1375      }
1376    }
1377    int largestRows=0;
1378    int largestColumns=0;
1379    for (int i=0;i<numberBlocks;i++) {
1380      if (blockCount[i]+nextColumn[i]>largestRows+largestColumns) {
1381        largestRows=blockCount[i];
1382        largestColumns=nextColumn[i];
1383      }
1384    }
1385    bool useful=true;
1386    if (numberMaster>halfway||largestRows*3>numberRows)
1387      useful=false;
1388    printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty, %d integer) out of %d\n",
1389           useful ? "**Useful" : "NoGood",
1390           numberBlocks,largestRows,largestColumns,numberMaster,numberEmpty,numberRows,
1391           numberMasterColumns,numberEmptyColumns,numberMasterIntegers,
1392           numberColumns);
1393    // columnBlock is columnBlock and blockStart is rowBlock
1394    // See if we want to compress
1395    if (!keepContinuous_) {
1396      // use blockEls
1397      int newNumber=0;
1398      for (int i=0;i<numberBlocks;i++) { 
1399        printf("Block %d has %d rows and %d columns (%d elements, %d integers)\n",
1400               i,blockCount[i],nextColumn[i],blockEls[i],countIntegers[i]);
1401        if (countIntegers[i]) {
1402          blockEls[i]=newNumber;
1403          newNumber++;
1404        } else {
1405          blockEls[i]=-1;
1406        }
1407      }
1408      for (int i=0;i<numberRows;i++) {
1409        int iBlock=blockStart[i];
1410        if (iBlock>=0)
1411          blockStart[i]=blockEls[iBlock];
1412      }
1413      for (int i=0;i<numberColumns;i++) {
1414        int iBlock=columnBlock[i];
1415        if (iBlock>=0)
1416          columnBlock[i]=blockEls[iBlock];
1417      }
1418      numberBlocks=newNumber;
1419    }
1420    // now set up structures
1421    numberBlocks_ = numberBlocks;
1422    // so callBack can modify
1423    whichRowBlock_ = blockStart;
1424    whichColumnBlock_ = columnBlock;
1425    intArray_=NULL;
1426    doubleArray_=NULL;
1427    (*(functionPointer_))(this,NULL,0);
1428
1429    saveLower_=CoinCopyOfArray(solver_->getColLower(),numberColumns);
1430    saveUpper_=CoinCopyOfArray(solver_->getColUpper(),numberColumns);
1431    startRowBlock_ = new int [numberBlocks_+2];
1432    backwardRow_ = new int [numberRows];
1433    rowsInBlock_ = new int [numberRows];
1434    whichRowBlock_ = new int [numberRows];
1435    startColumnBlock_ = new int [numberBlocks_+2];
1436    columnsInBlock_ = new int [numberColumns];
1437    whichColumnBlock_ = new int [numberColumns];
1438    intsInBlock_ = new int [numberBlocks_];
1439    // use for counts
1440    memset(rowsInBlock_,0,numberBlocks_*sizeof(int));
1441    numberMasterRows_=0;
1442    for (int i=0;i<numberRows;i++) {
1443      int iBlock=blockStart[i];
1444      if (iBlock>=0) {
1445        rowsInBlock_[iBlock]++;
1446        whichRowBlock_[i]=iBlock;
1447        backwardRow_[i]=-1;
1448      } else {
1449        whichRowBlock_[i]=-1;
1450        backwardRow_[i]=numberMasterRows_;
1451        numberMasterRows_++;
1452      }
1453    }
1454    memset(columnsInBlock_,0,numberBlocks_*sizeof(int));
1455    memset(intsInBlock_,0,numberBlocks_*sizeof(int));
1456    numberMasterColumns_=0;
1457    for (int i=0;i<numberColumns;i++) {
1458      int iBlock=columnBlock[i];
1459      if (iBlock>=0) {
1460        columnsInBlock_[iBlock]++;
1461        whichColumnBlock_[i]=iBlock;
1462        if (solver_->isInteger(i))
1463          intsInBlock_[iBlock]++;
1464      } else {
1465        whichColumnBlock_[i]=-1;
1466        numberMasterColumns_++;
1467      }
1468    }
1469    // starts
1470    int nRow=0;
1471    int nColumn=0;
1472    int maxIntsInBlock=0;
1473    for (int i=0;i<numberBlocks_;i++) {
1474      maxIntsInBlock = CoinMax(maxIntsInBlock,intsInBlock_[i]);
1475      startRowBlock_[i]=nRow;
1476      startColumnBlock_[i]=nColumn;
1477      nRow+=rowsInBlock_[i];
1478      nColumn+=columnsInBlock_[i];
1479    }
1480    // may not be used - but set anyway
1481    sizeFingerPrint_ = (maxIntsInBlock+31)/32;
1482    startRowBlock_[numberBlocks_]=nRow;
1483    startColumnBlock_[numberBlocks_]=nColumn;
1484    startRowBlock_[numberBlocks_+1]=numberRows;
1485    startColumnBlock_[numberBlocks_+1]=numberColumns;
1486    // do lists
1487    for ( int i=0 ; i<numberRows ; ++i ) {
1488      int iBlock = whichRowBlock_[i];
1489      if (iBlock<0)
1490        iBlock=numberBlocks_;
1491      int k=startRowBlock_[iBlock];
1492      startRowBlock_[iBlock]=k+1;
1493      rowsInBlock_[k]=i;
1494    }
1495    for (int i=numberBlocks+1;i>0;i--)
1496      startRowBlock_[i]=startRowBlock_[i-1];
1497    startRowBlock_[0]=0;
1498    for ( int i=0 ; i<numberColumns ; ++i ) {
1499      int iBlock = whichColumnBlock_[i];
1500      if (iBlock<0)
1501        iBlock=numberBlocks_;
1502      int k=startColumnBlock_[iBlock];
1503      startColumnBlock_[iBlock]=k+1;
1504      columnsInBlock_[k]=i;
1505    }
1506    for (int i=numberBlocks+1;i>0;i--)
1507      startColumnBlock_[i]=startColumnBlock_[i-1];
1508    startColumnBlock_[0]=0;
1509    if (numberBlocks_<10000) {
1510      affinity_ = new unsigned short [numberBlocks_*numberBlocks_];
1511      // compute space needed
1512      int * build = new int [numberMasterRows_];
1513      memset(build,0,numberMasterRows_*sizeof(int));
1514      int nSpace=0;
1515      for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
1516        int start=startColumnBlock_[iBlock];
1517        int end=startColumnBlock_[iBlock+1];
1518        for (int i=start;i<end;i++) {
1519          int iColumn=columnsInBlock_[i];
1520          for (CoinBigIndex j=columnStart[iColumn];
1521               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1522            int iRow = row[j];
1523            iRow=backwardRow_[iRow];
1524            if (iRow>=0) 
1525              build[iRow] ++;
1526          }
1527        }
1528        for (int i=0;i<numberMasterRows_;i++) {
1529          int value=build[i];
1530          if (value) {
1531            build[i]=0;
1532            nSpace++;
1533          }
1534        }
1535      }
1536      // get arrays
1537      int * starts = new int [numberBlocks_+1+2*nSpace];
1538      memset(affinity_,0,numberBlocks_*numberBlocks_*sizeof(unsigned short));
1539      // fill arrays
1540      int * rowM = starts+numberBlocks_+1;
1541      int * sumM = rowM+nSpace;
1542      nSpace=0;
1543      starts[0]=0;
1544      for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
1545        int start=startColumnBlock_[iBlock];
1546        int end=startColumnBlock_[iBlock+1];
1547        for (int i=start;i<end;i++) {
1548          int iColumn=columnsInBlock_[i];
1549          for (CoinBigIndex j=columnStart[iColumn];
1550               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1551            int iRow = row[j];
1552            iRow=backwardRow_[iRow];
1553            if (iRow>=0) 
1554              build[iRow] ++;
1555          }
1556        }
1557        for (int i=0;i<numberMasterRows_;i++) {
1558          int value=build[i];
1559          if (value) {
1560            build[i]=0;
1561            sumM[nSpace]=value;
1562            rowM[nSpace++]=i;
1563          }
1564        }
1565        starts[iBlock+1]=nSpace;
1566      }
1567      for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
1568        int startI=starts[iBlock];
1569        int endI=starts[iBlock+1];
1570        if (endI==startI)
1571          continue;
1572        for (int jBlock=iBlock+1;jBlock<numberBlocks_;jBlock++) {
1573          int startJ=starts[jBlock];
1574          int endJ=starts[jBlock+1];
1575          if (endJ==startJ)
1576            continue;
1577          double sum=0.0;
1578          int i=startI;
1579          int j=startJ;
1580          int rowI=rowM[i];
1581          int rowJ=rowM[j];
1582          while (rowI!=COIN_INT_MAX&&rowJ!=COIN_INT_MAX) {
1583            if (rowI<rowJ) {
1584              i++;
1585              if (i<endI)
1586                rowI=rowM[i];
1587              else
1588                rowI=COIN_INT_MAX;
1589            } else if (rowI>rowJ) {
1590              j++;
1591              if (j<endJ)
1592                rowJ=rowM[j];
1593              else
1594                rowJ=COIN_INT_MAX;
1595            } else {
1596              // bias ????????
1597              sum += sumM[i]*sumM[j];
1598              i++;
1599              if (i<endI)
1600                rowI=rowM[i];
1601              else
1602                rowI=COIN_INT_MAX;
1603              j++;
1604              if (j<endJ)
1605                rowJ=rowM[j];
1606              else
1607                rowJ=COIN_INT_MAX;
1608            }
1609          }
1610          if (sum>65535)
1611            sum = 65535;
1612          unsigned short value = static_cast<unsigned short>(sum);
1613          affinity_[iBlock*numberBlocks+jBlock]=value;
1614          affinity_[jBlock*numberBlocks+iBlock]=value;
1615        }
1616      }
1617      // statistics
1618      int nTotalZero=0;
1619      int base=0;
1620      for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
1621        int aff=0;
1622        int nZero=0;
1623        for (int jBlock=0;jBlock<numberBlocks_;jBlock++) {
1624          if (iBlock!=jBlock) {
1625            if (affinity_[base+jBlock])
1626              aff+=affinity_[base+jBlock];
1627            else
1628              nZero++;
1629          }
1630        }         
1631        //printf("Block %d has affinity %d but zero with %d blocks\n",
1632        //     iBlock,aff,nZero);
1633        nTotalZero+=nZero;
1634        base+=numberBlocks;
1635      }
1636      printf("Total not affinity %d - average %g%%\n",
1637             nTotalZero,100.0*(static_cast<double>(nTotalZero)
1638                               /(numberBlocks*numberBlocks)));
1639     
1640      delete [] starts;
1641      delete [] build;
1642    } else {
1643      printf("Too many blocks - no affinity\n");
1644    }
1645    if (fullDWEverySoOften_>0) {
1646      setupDWStructures();
1647    }
1648  }
1649  delete [] blockStart;
1650}
1651// Add DW proposals
1652int 
1653CbcHeuristicDW::addDW(const double * solution,int numberBlocksUsed, 
1654                      const int * whichBlocks)
1655{
1656  if (numberDW_+numberBlocksUsed>maximumDW_) {
1657    // extend
1658    int n = maximumDW_+5*numberBlocks_;
1659    double * weightsT = new double [n];
1660    int * dwBlockT = new int[n];
1661    unsigned int * fingerT = new unsigned int [n*sizeFingerPrint_];
1662    memcpy(weightsT,weights_,numberDW_*sizeof(double));
1663    memcpy(dwBlockT,dwBlock_,numberDW_*sizeof(int));
1664    memcpy(fingerT,fingerPrint_,
1665           numberDW_*sizeFingerPrint_*sizeof(unsigned int));
1666    delete [] weights_;
1667    weights_=weightsT;
1668    delete [] dwBlock_;
1669    dwBlock_=dwBlockT;
1670    delete [] fingerPrint_;
1671    fingerPrint_ = fingerT;
1672    maximumDW_=n;
1673  }
1674  //int numberColumns = solver_->getNumCols();
1675  //int numberRows = solver_->getNumRows();
1676  // get space to add elements
1677#define MAX_ADD 100000
1678  int * startsDW = new int[numberBlocks_+1+MAX_ADD];
1679  int * rowDW = startsDW+numberBlocks_+1;
1680  double * elementDW = new double[MAX_ADD+3*numberBlocks_+numberMasterRows_];
1681  double * newCost = elementDW+MAX_ADD;
1682  double * newLower = newCost+numberBlocks_;
1683  double * newUpper = newLower+numberBlocks_;
1684  double * build = newUpper+numberBlocks_;
1685  memset(build,0,numberMasterRows_*sizeof(double));
1686  int nAdd=0;
1687  int nTotalAdded=0;
1688  int nEls=0;
1689  startsDW[0]=0;
1690  // Column copy
1691  const double * element = solver_->getMatrixByCol()->getElements();
1692  const int * row = solver_->getMatrixByCol()->getIndices();
1693  const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
1694  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
1695  const double * objective = solver_->getObjCoefficients();
1696  for (int jBlock=0;jBlock<numberBlocksUsed;jBlock++) {
1697    double thisWeight=0.0;
1698    double thisWeightC=0.0;
1699    double thisCost=0.0;
1700    int nElInMaster=0;
1701    int nElIntInMaster=0;
1702    int nElIntInMaster1=0;
1703    int iBlock = whichBlocks[jBlock];
1704    int start=startColumnBlock_[iBlock];
1705    int end=startColumnBlock_[iBlock+1];
1706    unsigned int * finger = fingerPrint_+sizeFingerPrint_*(nAdd+numberDW_);
1707    memset(finger,0,sizeFingerPrint_*sizeof(unsigned int));
1708    int iBit=0;
1709    for (int i=start;i<end;i++) {
1710      int iColumn=columnsInBlock_[i];
1711      bool isInteger = solver_->isInteger(iColumn);
1712      double value = solution[iColumn];
1713      if (isInteger) {
1714        if(value>1.0e-6) 
1715          finger[0] |= 1<<iBit;
1716        iBit++;
1717        if (iBit==32) {
1718          finger++;
1719          iBit++;
1720        }
1721      }
1722      thisCost += value*objective[iColumn];
1723      for (CoinBigIndex j=columnStart[iColumn];
1724           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1725        int iRow = row[j];
1726        iRow=backwardRow_[iRow];
1727        if (iRow>=0) {
1728          nElInMaster++;
1729          double elementValue=element[j];
1730          build[iRow]+=value*elementValue;
1731          if (isInteger) {
1732            nElIntInMaster++;
1733            if (value)
1734              nElIntInMaster1++;
1735            thisWeight += random_[iRow]*value*elementValue;
1736          } else {
1737            value = 0.0001*floor(value*10000.0+0.5);
1738          }
1739          thisWeightC += random_[iRow]*value*elementValue;
1740        }
1741      }
1742    }
1743    // see if already in
1744    printf("block %d nel %d nelInt %d nelInt1 %d - weight %g (%g)\n",
1745           iBlock,nElInMaster,nElIntInMaster,nElIntInMaster1,
1746           thisWeight,thisWeightC);
1747    int iProposal;
1748    for (iProposal=0;iProposal<numberDW_;iProposal++) {
1749      if (iBlock==dwBlock_[iProposal]&&weights_[iProposal]==thisWeightC)
1750        break;
1751    }
1752    if (iProposal<numberDW_) {
1753      printf("above looks like duplicate\n");
1754      memset(build,0,numberMasterRows_*sizeof(double));
1755      //iProposal=numberDW;
1756    }
1757    if (iProposal==numberDW_) {
1758      // new
1759      // could build in stages
1760      assert (nEls+numberMasterRows_<MAX_ADD);
1761      for (int i=0;i<numberMasterRows_;i++) {
1762        double value=build[i];
1763        if (value) {
1764          build[i]=0.0;
1765          if (fabs(value)>1.0e-10) {
1766            elementDW[nEls]=value;
1767            rowDW[nEls++]=i;
1768          }
1769        }
1770      }
1771      // convexity
1772      elementDW[nEls]=1.0;
1773      rowDW[nEls++]=numberMasterRows_+iBlock;
1774      weights_[nAdd+numberDW_]=thisWeightC;
1775      dwBlock_[nAdd+numberDW_]=iBlock;
1776      newLower[nAdd]=0.0;
1777      newUpper[nAdd]=1.0;
1778      newCost[nAdd++]=thisCost;
1779      startsDW[nAdd]=nEls;
1780    }
1781    if (nEls+numberMasterRows_>MAX_ADD) {
1782      printf("Adding %d proposals with %d elements - out of room\n",
1783             nAdd,nEls);
1784      dwSolver_->addCols(nAdd,startsDW,rowDW,elementDW,newLower,
1785                        newUpper,newCost);
1786      numberDW_+=nAdd;
1787      nTotalAdded+=nAdd;
1788      nAdd=0;
1789      nEls=0;
1790    }
1791  }
1792  if (nAdd) {
1793    printf("Adding %d proposals with %d elements\n",
1794           nAdd,nEls);
1795    dwSolver_->addCols(nAdd,startsDW,rowDW,elementDW,newLower,
1796                     newUpper,newCost);
1797    nTotalAdded+=nAdd;
1798    numberDW_+=nAdd;
1799  }
1800  delete [] startsDW;
1801  delete [] elementDW;
1802  if (nTotalAdded) {
1803    double * objs = new double[numberDWTimes_+1];
1804    memcpy(objs,objectiveDW_,numberDWTimes_*sizeof(double));
1805    delete [] objectiveDW_;
1806    objectiveDW_=objs;
1807    int * temp = new int[numberDWTimes_+1];
1808    memcpy(temp,numberColumnsDW_,numberDWTimes_*sizeof(int));
1809    delete [] numberColumnsDW_;
1810    numberColumnsDW_=temp;
1811    numberColumnsDW_[numberDWTimes_]= dwSolver_->getNumCols();
1812    objectiveDW_[numberDWTimes_++]= objectiveValue(solution);
1813  }
1814  return nTotalAdded;
1815}
1816// Pass in a solution
1817void 
1818CbcHeuristicDW::passInSolution(const double * solution)
1819{
1820  if (fullDWEverySoOften_>0) {
1821    int * which = new int[numberBlocks_];
1822    for (int i=0;i<numberBlocks_;i++)
1823      which[i]=i;
1824    addDW(solution,numberBlocks_,which);
1825    delete [] which;
1826  }
1827  if (objectiveValue(solution)<bestObjective_-1.0e-5) {
1828    bestObjective_=objectiveValue(solution);
1829    int numberColumns = solver_->getNumCols();
1830    if (!bestSolution_) 
1831      bestSolution_ = new double [numberColumns];
1832    memcpy(bestSolution_,solution,numberColumns*sizeof(double));
1833  }
1834}
1835// Pass in continuous solution
1836void 
1837CbcHeuristicDW::passInContinuousSolution(const double * solution)
1838{
1839  int numberColumns = solver_->getNumCols();
1840  if (!continuousSolution_) 
1841    continuousSolution_ = new double [numberColumns];
1842  memcpy(continuousSolution_,solution,numberColumns*sizeof(double));
1843}
1844// Objective value (could also check validity)
1845double 
1846CbcHeuristicDW::objectiveValue(const double * solution)
1847{
1848  // compute objective value
1849  double objOffset = 0.0;
1850  solver_->getDblParam(OsiObjOffset, objOffset);
1851  double objectiveValue = -objOffset;
1852  int numberColumns = solver_->getNumCols();
1853  const double * objective = solver_->getObjCoefficients();
1854  for (int i=0;i<numberColumns;i++) {
1855    double value = solution[i];
1856    if (solver_->isInteger(i)) {
1857      if (fabs(value-floor(value+0.5))>1.0e-7)
1858        printf("Bad integer value for %d of %g\n",i,value);
1859    }
1860    objectiveValue += objective[i]*value;
1861  }
1862  return objectiveValue;
1863}
1864// Objective value when whichDw created
1865double 
1866CbcHeuristicDW::objectiveValueWhen(int whichDW) const
1867{
1868  if (whichDW>=numberDWTimes_)
1869    return COIN_DBL_MAX;
1870  else
1871    return objectiveDW_[whichDW];
1872}
1873// Number of columns in DW
1874int 
1875CbcHeuristicDW::numberColumnsDW(int whichDW) const
1876{
1877  if (whichDW>=numberDWTimes_)
1878    return COIN_INT_MAX;
1879  else
1880    return numberColumnsDW_[whichDW];
1881}
1882// DW model (user must delete)
1883OsiSolverInterface * 
1884CbcHeuristicDW::DWModel(int whichDW) const
1885{
1886  if (whichDW>=numberDWTimes_)
1887    return NULL;
1888  OsiSolverInterface * newSolver = dwSolver_->clone();
1889  int numberColumns2=newSolver->getNumCols();
1890  int numberColumns=numberColumnsDW_[whichDW];
1891  if (numberColumns<numberColumns2) {
1892    int * del = new int [numberColumns2-numberColumns];
1893    for (int i=numberColumns;i<numberColumns2;i++)
1894      del[i-numberColumns]=i;
1895    newSolver->deleteCols(numberColumns2-numberColumns,del);
1896    delete [] del;
1897  }
1898  // Set all to integer that need setting
1899  for (int i=numberMasterColumns_;i<numberColumns;i++) {
1900    newSolver->setContinuous(i);
1901  }
1902  int numberDW=numberColumns-numberMasterColumns_;
1903  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
1904    bool allSame=true;
1905    unsigned int * finger = fingerPrint_;
1906    unsigned int * fingerTest = NULL;
1907    for (int i=0;i<numberDW;i++) {
1908      if (dwBlock_[i]==iBlock) {
1909        if (fingerTest) {
1910          for (int j=0;j<sizeFingerPrint_;j++) {
1911            if (finger[j]!=fingerTest[j]) {
1912              allSame=false;
1913              break;
1914            }
1915          }
1916          if (!allSame)
1917            break;
1918        } else {
1919          fingerTest = finger;
1920        }
1921      }
1922      finger+=sizeFingerPrint_;
1923    }
1924    if (!allSame) {
1925      // Set all to integer that need setting
1926      for (int i=0;i<numberDW;i++) {
1927        if (iBlock==dwBlock_[i]) {
1928          int iColumn = numberMasterColumns_+i;
1929          newSolver->setInteger(iColumn);
1930        }
1931      }
1932    }
1933  }
1934  //newSolver->writeMps("dw","mps");
1935  return newSolver;
1936}
1937/* DW Proposal actions
1938   fullDWEverySoOften -
1939   0 - off
1940   k - every k times solution gets better
1941*/
1942void
1943CbcHeuristicDW::setProposalActions(int fullDWEverySoOften)
1944{ 
1945  fullDWEverySoOften_=fullDWEverySoOften;
1946  if (fullDWEverySoOften_>0&&!random_)
1947    setupDWStructures();
1948}
1949// Set up DW structure
1950void
1951CbcHeuristicDW::setupDWStructures()
1952{
1953  random_=new double [numberMasterRows_];
1954  for (int i=0;i<numberMasterRows_;i++)
1955    random_[i]=CoinDrand48();
1956  weights_ = new double [numberBlocks_];
1957  dwBlock_ = new int [numberBlocks_];
1958  fingerPrint_ = new unsigned int[numberBlocks_*sizeFingerPrint_];
1959  // create dwSolver
1960  int numberColumns = solver_->getNumCols();
1961  int numberRows = solver_->getNumRows();
1962  int * tempRow = new int [numberRows+numberColumns];
1963  int * tempColumn = tempRow + numberRows;
1964  int numberMasterRows=0;
1965  for (int i=0;i<numberRows;i++) {
1966    int iBlock=whichRowBlock_[i];
1967    if (iBlock<0) 
1968      tempRow[numberMasterRows++]=i;
1969  }
1970  int numberMasterColumns=0;
1971  for (int i=0;i<numberColumns;i++) {
1972    int iBlock=whichColumnBlock_[i];
1973    if (iBlock<0) 
1974      tempColumn[numberMasterColumns++]=i;
1975  }
1976  OsiClpSolverInterface * solver = 
1977    dynamic_cast<OsiClpSolverInterface *>(solver_);
1978  ClpSimplex * tempModel = new ClpSimplex(solver->getModelPtr(),
1979                                          numberMasterRows,tempRow,
1980                                          numberMasterColumns,tempColumn);
1981  // add convexity constraints
1982  double * rhs = new double[numberBlocks_];
1983  for (int i=0;i<numberBlocks_;i++)
1984    rhs[i]=1.0;
1985  tempModel->addRows(numberBlocks_,rhs,rhs,NULL,NULL,NULL);
1986  delete [] rhs;
1987  dwSolver_ = new OsiClpSolverInterface(tempModel,true);
1988  printf("DW model has %d master rows, %d master columns and %d convexity rows\n",
1989         numberMasterRows,numberMasterColumns,numberBlocks_);
1990  // do master integers
1991  for (int i=0;i<numberMasterColumns;i++) {
1992    int iColumn=tempColumn[i];
1993    if (solver->isInteger(iColumn))
1994      dwSolver_->setInteger(i);
1995  }
1996  delete [] tempRow;
1997}
1998
1999
Note: See TracBrowser for help on using the repository browser.