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

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

yet more changes for dubious idea

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