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

Last change on this file since 2076 was 1957, checked in by forrest, 6 years ago

minor fixes and take out printout

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