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

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

add more twiddles for advanced use of heuristics

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