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

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

correct some printing and allow RINS to start more easily in mini B&B

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