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

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

adding a dubious heuristic

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