source: trunk/Cbc/src/CbcHeuristicDW.cpp

Last change on this file was 2467, checked in by unxusr, 10 months ago

spaces after angles

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