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

Last change on this file since 2094 was 2094, checked in by forrest, 4 years ago

for memory leaks and heuristics and some experimental stuff

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