source: trunk/Cbc/src/CbcHeuristicGreedy.cpp

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

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 60.9 KB
Line 
1/* $Id: CbcHeuristicGreedy.cpp 2467 2019-01-03 21:26:29Z unxusr $ */
2// Copyright (C) 2005, 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 "CbcModel.hpp"
17#include "CbcStrategy.hpp"
18#include "CbcHeuristicGreedy.hpp"
19#include "CoinSort.hpp"
20#include "CglPreProcess.hpp"
21// Default Constructor
22CbcHeuristicGreedyCover::CbcHeuristicGreedyCover()
23  : CbcHeuristic()
24{
25  // matrix  will automatically be empty
26  originalNumberRows_ = 0;
27  algorithm_ = 0;
28  numberTimes_ = 100;
29}
30
31// Constructor from model
32CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(CbcModel &model)
33  : CbcHeuristic(model)
34{
35  gutsOfConstructor(&model);
36  algorithm_ = 0;
37  numberTimes_ = 100;
38  whereFrom_ = 1;
39}
40
41// Destructor
42CbcHeuristicGreedyCover::~CbcHeuristicGreedyCover()
43{
44}
45
46// Clone
47CbcHeuristic *
48CbcHeuristicGreedyCover::clone() const
49{
50  return new CbcHeuristicGreedyCover(*this);
51}
52// Guts of constructor from a CbcModel
53void CbcHeuristicGreedyCover::gutsOfConstructor(CbcModel *model)
54{
55  model_ = model;
56  // Get a copy of original matrix
57  assert(model->solver());
58  if (model->solver()->getNumRows()) {
59    matrix_ = *model->solver()->getMatrixByCol();
60  }
61  originalNumberRows_ = model->solver()->getNumRows();
62}
63// Create C++ lines to get to current state
64void CbcHeuristicGreedyCover::generateCpp(FILE *fp)
65{
66  CbcHeuristicGreedyCover other;
67  fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
68  fprintf(fp, "3  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n");
69  CbcHeuristic::generateCpp(fp, "heuristicGreedyCover");
70  if (algorithm_ != other.algorithm_)
71    fprintf(fp, "3  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
72  else
73    fprintf(fp, "4  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
74  if (numberTimes_ != other.numberTimes_)
75    fprintf(fp, "3  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
76  else
77    fprintf(fp, "4  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
78  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyCover);\n");
79}
80
81// Copy constructor
82CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(const CbcHeuristicGreedyCover &rhs)
83  : CbcHeuristic(rhs)
84  , matrix_(rhs.matrix_)
85  , originalNumberRows_(rhs.originalNumberRows_)
86  , algorithm_(rhs.algorithm_)
87  , numberTimes_(rhs.numberTimes_)
88{
89}
90
91// Assignment operator
92CbcHeuristicGreedyCover &
93CbcHeuristicGreedyCover::operator=(const CbcHeuristicGreedyCover &rhs)
94{
95  if (this != &rhs) {
96    CbcHeuristic::operator=(rhs);
97    matrix_ = rhs.matrix_;
98    originalNumberRows_ = rhs.originalNumberRows_;
99    algorithm_ = rhs.algorithm_;
100    numberTimes_ = rhs.numberTimes_;
101  }
102  return *this;
103}
104// Returns 1 if solution, 0 if not
105int CbcHeuristicGreedyCover::solution(double &solutionValue,
106  double *betterSolution)
107{
108  numCouldRun_++;
109  if (!model_)
110    return 0;
111  // See if to do
112  if (!when() || (when() == 1 && model_->phase() != 1))
113    return 0; // switched off
114  if (model_->getNodeCount() > numberTimes_)
115    return 0;
116#ifdef HEURISTIC_INFORM
117  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
118    heuristicName(), numRuns_, numCouldRun_, when_);
119#endif
120  // See if at root node
121  bool atRoot = model_->getNodeCount() == 0;
122  int passNumber = model_->getCurrentPassNumber();
123  if (atRoot && passNumber > 1)
124    return 0;
125  OsiSolverInterface *solver = model_->solver();
126  const double *columnLower = solver->getColLower();
127  const double *columnUpper = solver->getColUpper();
128  // And original upper bounds in case we want to use them
129  const double *originalUpper = model_->continuousSolver()->getColUpper();
130  // But not if algorithm says so
131  if ((algorithm_ % 10) == 0)
132    originalUpper = columnUpper;
133  const double *rowLower = solver->getRowLower();
134  const double *solution = solver->getColSolution();
135  const double *objective = solver->getObjCoefficients();
136  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
137  double primalTolerance;
138  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
139
140  // This is number of rows when matrix was passed in
141  int numberRows = originalNumberRows_;
142  if (!numberRows)
143    return 0; // switched off
144
145  numRuns_++;
146  assert(numberRows == matrix_.getNumRows());
147  int iRow, iColumn;
148  double direction = solver->getObjSense();
149  double offset;
150  solver->getDblParam(OsiObjOffset, offset);
151  double newSolutionValue = -offset;
152  int returnCode = 0;
153
154  // Column copy
155  const double *element = matrix_.getElements();
156  const int *row = matrix_.getIndices();
157  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
158  const int *columnLength = matrix_.getVectorLengths();
159
160  // Get solution array for heuristic solution
161  int numberColumns = solver->getNumCols();
162  double *newSolution = new double[numberColumns];
163  double *rowActivity = new double[numberRows];
164  memset(rowActivity, 0, numberRows * sizeof(double));
165  bool allOnes = true;
166  // Get rounded down solution
167  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
168    CoinBigIndex j;
169    double value = solution[iColumn];
170    if (isHeuristicInteger(solver, iColumn)) {
171      // Round down integer
172      if (fabs(floor(value + 0.5) - value) < integerTolerance) {
173        value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
174      } else {
175        value = CoinMax(floor(value), columnLower[iColumn]);
176      }
177    }
178    // make sure clean
179    value = CoinMin(value, columnUpper[iColumn]);
180    value = CoinMax(value, columnLower[iColumn]);
181    newSolution[iColumn] = value;
182    double cost = direction * objective[iColumn];
183    newSolutionValue += value * cost;
184    for (j = columnStart[iColumn];
185         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
186      int iRow = row[j];
187      rowActivity[iRow] += value * element[j];
188      if (element[j] != 1.0)
189        allOnes = false;
190    }
191  }
192  // See if we round up
193  bool roundup = ((algorithm_ % 100) != 0);
194  if (roundup && allOnes) {
195    // Get rounded up solution
196    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
197      CoinBigIndex j;
198      double value = solution[iColumn];
199      if (isHeuristicInteger(solver, iColumn)) {
200        // but round up if no activity
201        if (roundup && value >= 0.499999 && !newSolution[iColumn]) {
202          bool choose = true;
203          for (j = columnStart[iColumn];
204               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
205            int iRow = row[j];
206            if (rowActivity[iRow]) {
207              choose = false;
208              break;
209            }
210          }
211          if (choose) {
212            newSolution[iColumn] = 1.0;
213            double cost = direction * objective[iColumn];
214            newSolutionValue += cost;
215            for (j = columnStart[iColumn];
216                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
217              int iRow = row[j];
218              rowActivity[iRow] += 1.0;
219            }
220          }
221        }
222      }
223    }
224  }
225  // Get initial list
226  int *which = new int[numberColumns];
227  for (iColumn = 0; iColumn < numberColumns; iColumn++)
228    which[iColumn] = iColumn;
229  int numberLook = numberColumns;
230  // See if we want to perturb more
231  double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
232  // Keep going round until a solution
233  while (true) {
234    // Get column with best ratio
235    int bestColumn = -1;
236    double bestRatio = COIN_DBL_MAX;
237    double bestStepSize = 0.0;
238    int newNumber = 0;
239    for (int jColumn = 0; jColumn < numberLook; jColumn++) {
240      int iColumn = which[jColumn];
241      CoinBigIndex j;
242      double value = newSolution[iColumn];
243      double cost = direction * objective[iColumn];
244      if (isHeuristicInteger(solver, iColumn)) {
245        // use current upper or original upper
246        if (value + 0.99 < originalUpper[iColumn]) {
247          double sum = 0.0;
248          int numberExact = 0;
249          for (j = columnStart[iColumn];
250               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
251            int iRow = row[j];
252            double gap = rowLower[iRow] - rowActivity[iRow];
253            double elementValue = allOnes ? 1.0 : element[j];
254            if (gap > 1.0e-7) {
255              sum += CoinMin(elementValue, gap);
256              if (fabs(elementValue - gap) < 1.0e-7)
257                numberExact++;
258            }
259          }
260          // could bias if exact
261          if (sum > 0.0) {
262            // add to next time
263            which[newNumber++] = iColumn;
264            double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
265            // If at root choose first
266            if (atRoot)
267              ratio = iColumn;
268            if (ratio < bestRatio) {
269              bestRatio = ratio;
270              bestColumn = iColumn;
271              bestStepSize = 1.0;
272            }
273          }
274        }
275      } else {
276        // continuous
277        if (value < columnUpper[iColumn]) {
278          // Go through twice - first to get step length
279          double step = 1.0e50;
280          for (j = columnStart[iColumn];
281               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
282            int iRow = row[j];
283            if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 && element[j] * step + rowActivity[iRow] >= rowLower[iRow]) {
284              step = (rowLower[iRow] - rowActivity[iRow]) / element[j];
285              ;
286            }
287          }
288          // now ratio
289          if (step < 1.0e50) {
290            // add to next time
291            which[newNumber++] = iColumn;
292            assert(step > 0.0);
293            double sum = 0.0;
294            for (j = columnStart[iColumn];
295                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
296              int iRow = row[j];
297              double newActivity = element[j] * step + rowActivity[iRow];
298              if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 && newActivity >= rowLower[iRow] - 1.0e-12) {
299                sum += element[j];
300              }
301            }
302            assert(sum > 0.0);
303            double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
304            if (ratio < bestRatio) {
305              bestRatio = ratio;
306              bestColumn = iColumn;
307              bestStepSize = step;
308            }
309          }
310        }
311      }
312    }
313    if (bestColumn < 0)
314      break; // we have finished
315    // Increase chosen column
316    newSolution[bestColumn] += bestStepSize;
317    double cost = direction * objective[bestColumn];
318    newSolutionValue += bestStepSize * cost;
319    for (CoinBigIndex j = columnStart[bestColumn];
320         j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
321      int iRow = row[j];
322      rowActivity[iRow] += bestStepSize * element[j];
323    }
324  }
325  delete[] which;
326  if (newSolutionValue < solutionValue) {
327    // check feasible
328    memset(rowActivity, 0, numberRows * sizeof(double));
329    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
330      CoinBigIndex j;
331      double value = newSolution[iColumn];
332      if (value) {
333        for (j = columnStart[iColumn];
334             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
335          int iRow = row[j];
336          rowActivity[iRow] += value * element[j];
337        }
338      }
339    }
340    // check was approximately feasible
341    bool feasible = true;
342    for (iRow = 0; iRow < numberRows; iRow++) {
343      if (rowActivity[iRow] < rowLower[iRow]) {
344        if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance)
345          feasible = false;
346      }
347    }
348    if (feasible) {
349      // new solution
350      memcpy(betterSolution, newSolution, numberColumns * sizeof(double));
351      solutionValue = newSolutionValue;
352      //printf("** Solution of %g found by rounding\n",newSolutionValue);
353      returnCode = 1;
354    } else {
355      // Can easily happen
356      //printf("Debug CbcHeuristicGreedyCover giving bad solution\n");
357    }
358  }
359  delete[] newSolution;
360  delete[] rowActivity;
361  return returnCode;
362}
363// update model
364void CbcHeuristicGreedyCover::setModel(CbcModel *model)
365{
366  gutsOfConstructor(model);
367  validate();
368}
369// Resets stuff if model changes
370void CbcHeuristicGreedyCover::resetModel(CbcModel *model)
371{
372  gutsOfConstructor(model);
373}
374// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
375void CbcHeuristicGreedyCover::validate()
376{
377  if (model_ && when() < 10) {
378    if (model_->numberIntegers() != model_->numberObjects() && (model_->numberObjects() || (model_->specialOptions() & 1024) == 0)) {
379      int numberOdd = 0;
380      for (int i = 0; i < model_->numberObjects(); i++) {
381        if (!model_->object(i)->canDoHeuristics())
382          numberOdd++;
383      }
384      if (numberOdd)
385        setWhen(0);
386    }
387    // Only works if costs positive, coefficients positive and all rows G
388    OsiSolverInterface *solver = model_->solver();
389    const double *columnLower = solver->getColLower();
390    const double *rowUpper = solver->getRowUpper();
391    const double *objective = solver->getObjCoefficients();
392    double direction = solver->getObjSense();
393
394    int numberRows = solver->getNumRows();
395    int numberColumns = solver->getNumCols();
396    // Column copy
397    matrix_.setDimensions(numberRows, numberColumns);
398    const double *element = matrix_.getElements();
399    const CoinBigIndex *columnStart = matrix_.getVectorStarts();
400    const int *columnLength = matrix_.getVectorLengths();
401    bool good = true;
402    for (int iRow = 0; iRow < numberRows; iRow++) {
403      if (rowUpper[iRow] < 1.0e30)
404        good = false;
405    }
406    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
407      if (objective[iColumn] * direction < 0.0)
408        good = false;
409      if (columnLower[iColumn] < 0.0)
410        good = false;
411      CoinBigIndex j;
412      for (j = columnStart[iColumn];
413           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
414        if (element[j] < 0.0)
415          good = false;
416      }
417    }
418    if (!good)
419      setWhen(0); // switch off
420  }
421}
422// Default Constructor
423CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality()
424  : CbcHeuristic()
425{
426  // matrix  will automatically be empty
427  fraction_ = 1.0; // no branch and bound
428  originalNumberRows_ = 0;
429  algorithm_ = 0;
430  numberTimes_ = 100;
431  whereFrom_ = 1;
432}
433
434// Constructor from model
435CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel &model)
436  : CbcHeuristic(model)
437{
438  // Get a copy of original matrix
439  gutsOfConstructor(&model);
440  fraction_ = 1.0; // no branch and bound
441  algorithm_ = 0;
442  numberTimes_ = 100;
443  whereFrom_ = 1;
444}
445
446// Destructor
447CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality()
448{
449}
450
451// Clone
452CbcHeuristic *
453CbcHeuristicGreedyEquality::clone() const
454{
455  return new CbcHeuristicGreedyEquality(*this);
456}
457// Guts of constructor from a CbcModel
458void CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel *model)
459{
460  model_ = model;
461  // Get a copy of original matrix
462  assert(model->solver());
463  if (model->solver()->getNumRows()) {
464    matrix_ = *model->solver()->getMatrixByCol();
465  }
466  originalNumberRows_ = model->solver()->getNumRows();
467}
468// Create C++ lines to get to current state
469void CbcHeuristicGreedyEquality::generateCpp(FILE *fp)
470{
471  CbcHeuristicGreedyEquality other;
472  fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
473  fprintf(fp, "3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
474  CbcHeuristic::generateCpp(fp, "heuristicGreedyEquality");
475  if (algorithm_ != other.algorithm_)
476    fprintf(fp, "3  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
477  else
478    fprintf(fp, "4  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
479  if (fraction_ != other.fraction_)
480    fprintf(fp, "3  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
481  else
482    fprintf(fp, "4  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
483  if (numberTimes_ != other.numberTimes_)
484    fprintf(fp, "3  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
485  else
486    fprintf(fp, "4  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
487  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
488}
489
490// Copy constructor
491CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality &rhs)
492  : CbcHeuristic(rhs)
493  , matrix_(rhs.matrix_)
494  , fraction_(rhs.fraction_)
495  , originalNumberRows_(rhs.originalNumberRows_)
496  , algorithm_(rhs.algorithm_)
497  , numberTimes_(rhs.numberTimes_)
498{
499}
500
501// Assignment operator
502CbcHeuristicGreedyEquality &
503CbcHeuristicGreedyEquality::operator=(const CbcHeuristicGreedyEquality &rhs)
504{
505  if (this != &rhs) {
506    CbcHeuristic::operator=(rhs);
507    matrix_ = rhs.matrix_;
508    fraction_ = rhs.fraction_;
509    originalNumberRows_ = rhs.originalNumberRows_;
510    algorithm_ = rhs.algorithm_;
511    numberTimes_ = rhs.numberTimes_;
512  }
513  return *this;
514}
515// Returns 1 if solution, 0 if not
516int CbcHeuristicGreedyEquality::solution(double &solutionValue,
517  double *betterSolution)
518{
519  numCouldRun_++;
520  if (!model_)
521    return 0;
522  // See if to do
523  if (!when() || (when() == 1 && model_->phase() != 1))
524    return 0; // switched off
525  if (model_->getNodeCount() > numberTimes_)
526    return 0;
527  // See if at root node
528  bool atRoot = model_->getNodeCount() == 0;
529  int passNumber = model_->getCurrentPassNumber();
530  if (atRoot && passNumber > 1)
531    return 0;
532  OsiSolverInterface *solver = model_->solver();
533  const double *columnLower = solver->getColLower();
534  const double *columnUpper = solver->getColUpper();
535  // And original upper bounds in case we want to use them
536  const double *originalUpper = model_->continuousSolver()->getColUpper();
537  // But not if algorithm says so
538  if ((algorithm_ % 10) == 0)
539    originalUpper = columnUpper;
540  const double *rowLower = solver->getRowLower();
541  const double *rowUpper = solver->getRowUpper();
542  const double *solution = solver->getColSolution();
543  const double *objective = solver->getObjCoefficients();
544  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
545  double primalTolerance;
546  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
547
548  // This is number of rows when matrix was passed in
549  int numberRows = originalNumberRows_;
550  if (!numberRows)
551    return 0; // switched off
552  numRuns_++;
553
554  assert(numberRows == matrix_.getNumRows());
555  int iRow, iColumn;
556  double direction = solver->getObjSense();
557  double offset;
558  solver->getDblParam(OsiObjOffset, offset);
559  double newSolutionValue = -offset;
560  int returnCode = 0;
561
562  // Column copy
563  const double *element = matrix_.getElements();
564  const int *row = matrix_.getIndices();
565  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
566  const int *columnLength = matrix_.getVectorLengths();
567
568  // Get solution array for heuristic solution
569  int numberColumns = solver->getNumCols();
570  double *newSolution = new double[numberColumns];
571  double *rowActivity = new double[numberRows];
572  memset(rowActivity, 0, numberRows * sizeof(double));
573  double rhsNeeded = 0;
574  for (iRow = 0; iRow < numberRows; iRow++)
575    rhsNeeded += rowUpper[iRow];
576  rhsNeeded *= fraction_;
577  bool allOnes = true;
578  // Get rounded down solution
579  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
580    CoinBigIndex j;
581    double value = solution[iColumn];
582    if (isHeuristicInteger(solver, iColumn)) {
583      // Round down integer
584      if (fabs(floor(value + 0.5) - value) < integerTolerance) {
585        value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
586      } else {
587        value = CoinMax(floor(value), columnLower[iColumn]);
588      }
589    }
590    // make sure clean
591    value = CoinMin(value, columnUpper[iColumn]);
592    value = CoinMax(value, columnLower[iColumn]);
593    newSolution[iColumn] = value;
594    double cost = direction * objective[iColumn];
595    newSolutionValue += value * cost;
596    for (j = columnStart[iColumn];
597         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
598      int iRow = row[j];
599      rowActivity[iRow] += value * element[j];
600      rhsNeeded -= value * element[j];
601      if (element[j] != 1.0)
602        allOnes = false;
603    }
604  }
605  // See if we round up
606  bool roundup = ((algorithm_ % 100) != 0);
607  if (roundup && allOnes) {
608    // Get rounded up solution
609    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
610      CoinBigIndex j;
611      double value = solution[iColumn];
612      if (isHeuristicInteger(solver, iColumn)) {
613        // but round up if no activity
614        if (roundup && value >= 0.6 && !newSolution[iColumn]) {
615          bool choose = true;
616          for (j = columnStart[iColumn];
617               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
618            int iRow = row[j];
619            if (rowActivity[iRow]) {
620              choose = false;
621              break;
622            }
623          }
624          if (choose) {
625            newSolution[iColumn] = 1.0;
626            double cost = direction * objective[iColumn];
627            newSolutionValue += cost;
628            for (j = columnStart[iColumn];
629                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
630              int iRow = row[j];
631              rowActivity[iRow] += 1.0;
632              rhsNeeded -= 1.0;
633            }
634          }
635        }
636      }
637    }
638  }
639  // Get initial list
640  int *which = new int[numberColumns];
641  for (iColumn = 0; iColumn < numberColumns; iColumn++)
642    which[iColumn] = iColumn;
643  int numberLook = numberColumns;
644  // See if we want to perturb more
645  double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
646  // Keep going round until a solution
647  while (true) {
648    // Get column with best ratio
649    int bestColumn = -1;
650    double bestRatio = COIN_DBL_MAX;
651    double bestStepSize = 0.0;
652    int newNumber = 0;
653    for (int jColumn = 0; jColumn < numberLook; jColumn++) {
654      int iColumn = which[jColumn];
655      CoinBigIndex j;
656      double value = newSolution[iColumn];
657      double cost = direction * objective[iColumn];
658      if (isHeuristicInteger(solver, iColumn)) {
659        // use current upper or original upper
660        if (value + 0.9999 < originalUpper[iColumn]) {
661          double movement = 1.0;
662          double sum = 0.0;
663          for (j = columnStart[iColumn];
664               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
665            int iRow = row[j];
666            double gap = rowUpper[iRow] - rowActivity[iRow];
667            double elementValue = allOnes ? 1.0 : element[j];
668            sum += elementValue;
669            if (movement * elementValue > gap) {
670              movement = gap / elementValue;
671            }
672          }
673          if (movement > 0.999999) {
674            // add to next time
675            which[newNumber++] = iColumn;
676            double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
677            // If at root
678            if (atRoot) {
679              if (fraction_ == 1.0)
680                ratio = iColumn; // choose first
681              else
682                ratio = -solution[iColumn]; // choose largest
683            }
684            if (ratio < bestRatio) {
685              bestRatio = ratio;
686              bestColumn = iColumn;
687              bestStepSize = 1.0;
688            }
689          }
690        }
691      } else {
692        // continuous
693        if (value < columnUpper[iColumn]) {
694          double movement = 1.0e50;
695          double sum = 0.0;
696          for (j = columnStart[iColumn];
697               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
698            int iRow = row[j];
699            if (element[j] * movement + rowActivity[iRow] > rowUpper[iRow]) {
700              movement = (rowUpper[iRow] - rowActivity[iRow]) / element[j];
701              ;
702            }
703            sum += element[j];
704          }
705          // now ratio
706          if (movement > 1.0e-7) {
707            // add to next time
708            which[newNumber++] = iColumn;
709            double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
710            if (ratio < bestRatio) {
711              bestRatio = ratio;
712              bestColumn = iColumn;
713              bestStepSize = movement;
714            }
715          }
716        }
717      }
718    }
719    if (bestColumn < 0)
720      break; // we have finished
721    // Increase chosen column
722    newSolution[bestColumn] += bestStepSize;
723    double cost = direction * objective[bestColumn];
724    newSolutionValue += bestStepSize * cost;
725    for (CoinBigIndex j = columnStart[bestColumn];
726         j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
727      int iRow = row[j];
728      rowActivity[iRow] += bestStepSize * element[j];
729      rhsNeeded -= bestStepSize * element[j];
730    }
731    if (rhsNeeded < 1.0e-8)
732      break;
733  }
734  delete[] which;
735  if (fraction_ < 1.0 && rhsNeeded < 1.0e-8 && newSolutionValue < solutionValue) {
736    // do branch and cut
737    // fix all nonzero
738    OsiSolverInterface *newSolver = model_->continuousSolver()->clone();
739    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
740      if (isHeuristicInteger(newSolver, iColumn))
741        newSolver->setColLower(iColumn, newSolution[iColumn]);
742    }
743    int returnCode = smallBranchAndBound(newSolver, 200, newSolution, newSolutionValue,
744      solutionValue, "CbcHeuristicGreedy");
745    if (returnCode < 0)
746      returnCode = 0; // returned on size
747    if ((returnCode & 2) != 0) {
748      // could add cut
749      returnCode &= ~2;
750    }
751    rhsNeeded = 1.0 - returnCode;
752    delete newSolver;
753  }
754  if (newSolutionValue < solutionValue && rhsNeeded < 1.0e-8) {
755    // check feasible
756    memset(rowActivity, 0, numberRows * sizeof(double));
757    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
758      CoinBigIndex j;
759      double value = newSolution[iColumn];
760      if (value) {
761        for (j = columnStart[iColumn];
762             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
763          int iRow = row[j];
764          rowActivity[iRow] += value * element[j];
765        }
766      }
767    }
768    // check was approximately feasible
769    bool feasible = true;
770    for (iRow = 0; iRow < numberRows; iRow++) {
771      if (rowActivity[iRow] < rowLower[iRow]) {
772        if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance)
773          feasible = false;
774      }
775    }
776    if (feasible) {
777      // new solution
778      memcpy(betterSolution, newSolution, numberColumns * sizeof(double));
779      solutionValue = newSolutionValue;
780      returnCode = 1;
781    }
782  }
783  delete[] newSolution;
784  delete[] rowActivity;
785  if (atRoot && fraction_ == 1.0) {
786    // try quick search
787    fraction_ = 0.4;
788    int newCode = this->solution(solutionValue, betterSolution);
789    if (newCode)
790      returnCode = 1;
791    fraction_ = 1.0;
792  }
793  return returnCode;
794}
795// update model
796void CbcHeuristicGreedyEquality::setModel(CbcModel *model)
797{
798  gutsOfConstructor(model);
799  validate();
800}
801// Resets stuff if model changes
802void CbcHeuristicGreedyEquality::resetModel(CbcModel *model)
803{
804  gutsOfConstructor(model);
805}
806// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
807void CbcHeuristicGreedyEquality::validate()
808{
809  if (model_ && when() < 10) {
810    if (model_->numberIntegers() != model_->numberObjects())
811      setWhen(0);
812    // Only works if costs positive, coefficients positive and all rows E or L
813    // And if values are integer
814    OsiSolverInterface *solver = model_->solver();
815    const double *columnLower = solver->getColLower();
816    const double *rowUpper = solver->getRowUpper();
817    const double *rowLower = solver->getRowLower();
818    const double *objective = solver->getObjCoefficients();
819    double direction = solver->getObjSense();
820
821    int numberRows = solver->getNumRows();
822    int numberColumns = solver->getNumCols();
823    matrix_.setDimensions(numberRows, numberColumns);
824    // Column copy
825    const double *element = matrix_.getElements();
826    const CoinBigIndex *columnStart = matrix_.getVectorStarts();
827    const int *columnLength = matrix_.getVectorLengths();
828    bool good = true;
829    for (int iRow = 0; iRow < numberRows; iRow++) {
830      if (rowUpper[iRow] > 1.0e30)
831        good = false;
832      if (rowLower[iRow] > 0.0 && rowLower[iRow] != rowUpper[iRow])
833        good = false;
834      if (floor(rowUpper[iRow] + 0.5) != rowUpper[iRow])
835        good = false;
836    }
837    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
838      if (objective[iColumn] * direction < 0.0)
839        good = false;
840      if (columnLower[iColumn] < 0.0)
841        good = false;
842      CoinBigIndex j;
843      for (j = columnStart[iColumn];
844           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
845        if (element[j] < 0.0)
846          good = false;
847        if (floor(element[j] + 0.5) != element[j])
848          good = false;
849      }
850    }
851    if (!good)
852      setWhen(0); // switch off
853  }
854}
855
856// Default Constructor
857CbcHeuristicGreedySOS::CbcHeuristicGreedySOS()
858  : CbcHeuristic()
859{
860  originalRhs_ = NULL;
861  // matrix  will automatically be empty
862  originalNumberRows_ = 0;
863  algorithm_ = 0;
864  numberTimes_ = 100;
865}
866
867// Constructor from model
868CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(CbcModel &model)
869  : CbcHeuristic(model)
870{
871  gutsOfConstructor(&model);
872  algorithm_ = 2;
873  numberTimes_ = 100;
874  whereFrom_ = 1;
875}
876
877// Destructor
878CbcHeuristicGreedySOS::~CbcHeuristicGreedySOS()
879{
880  delete[] originalRhs_;
881}
882
883// Clone
884CbcHeuristic *
885CbcHeuristicGreedySOS::clone() const
886{
887  return new CbcHeuristicGreedySOS(*this);
888}
889// Guts of constructor from a CbcModel
890void CbcHeuristicGreedySOS::gutsOfConstructor(CbcModel *model)
891{
892  model_ = model;
893  // Get a copy of original matrix
894  assert(model->solver());
895  if (model->solver()->getNumRows()) {
896    matrix_ = *model->solver()->getMatrixByCol();
897  }
898  originalNumberRows_ = model->solver()->getNumRows();
899  originalRhs_ = new double[originalNumberRows_];
900}
901// Create C++ lines to get to current state
902void CbcHeuristicGreedySOS::generateCpp(FILE *fp)
903{
904  CbcHeuristicGreedySOS other;
905  fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
906  fprintf(fp, "3  CbcHeuristicGreedySOS heuristicGreedySOS(*cbcModel);\n");
907  CbcHeuristic::generateCpp(fp, "heuristicGreedySOS");
908  if (algorithm_ != other.algorithm_)
909    fprintf(fp, "3  heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_);
910  else
911    fprintf(fp, "4  heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_);
912  if (numberTimes_ != other.numberTimes_)
913    fprintf(fp, "3  heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_);
914  else
915    fprintf(fp, "4  heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_);
916  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedySOS);\n");
917}
918
919// Copy constructor
920CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(const CbcHeuristicGreedySOS &rhs)
921  : CbcHeuristic(rhs)
922  , matrix_(rhs.matrix_)
923  , originalNumberRows_(rhs.originalNumberRows_)
924  , algorithm_(rhs.algorithm_)
925  , numberTimes_(rhs.numberTimes_)
926{
927  originalRhs_ = CoinCopyOfArray(rhs.originalRhs_, originalNumberRows_);
928}
929
930// Assignment operator
931CbcHeuristicGreedySOS &
932CbcHeuristicGreedySOS::operator=(const CbcHeuristicGreedySOS &rhs)
933{
934  if (this != &rhs) {
935    CbcHeuristic::operator=(rhs);
936    matrix_ = rhs.matrix_;
937    originalNumberRows_ = rhs.originalNumberRows_;
938    algorithm_ = rhs.algorithm_;
939    numberTimes_ = rhs.numberTimes_;
940    delete[] originalRhs_;
941    originalRhs_ = CoinCopyOfArray(rhs.originalRhs_, originalNumberRows_);
942  }
943  return *this;
944}
945// Returns 1 if solution, 0 if not
946int CbcHeuristicGreedySOS::solution(double &solutionValue,
947  double *betterSolution)
948{
949  numCouldRun_++;
950  if (!model_)
951    return 0;
952  // See if to do
953  if (!when() || (when() == 1 && model_->phase() != 1))
954    return 0; // switched off
955  if (model_->getNodeCount() > numberTimes_)
956    return 0;
957  // See if at root node
958  bool atRoot = model_->getNodeCount() == 0;
959  int passNumber = model_->getCurrentPassNumber();
960  if (atRoot && passNumber > 1)
961    return 0;
962  OsiSolverInterface *solver = model_->solver();
963  int numberColumns = solver->getNumCols();
964  // This is number of rows when matrix was passed in
965  int numberRows = originalNumberRows_;
966  if (!numberRows)
967    return 0; // switched off
968
969  const double *columnLower = solver->getColLower();
970  const double *columnUpper = solver->getColUpper();
971  // modified rhs
972  double *rhs = CoinCopyOfArray(originalRhs_, numberRows);
973  // Column copy
974  const double *element = matrix_.getElements();
975  const int *row = matrix_.getIndices();
976  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
977  const int *columnLength = matrix_.getVectorLengths();
978  int *sosRow = new int[numberColumns];
979  char *sos = new char[numberRows];
980  memset(sos, 'a', numberRows);
981  int nonSOS = 0;
982  // If bit set then use current
983  if ((algorithm_ & 1) != 0) {
984    const CoinPackedMatrix *matrix = solver->getMatrixByCol();
985    element = matrix->getElements();
986    row = matrix->getIndices();
987    columnStart = matrix->getVectorStarts();
988    columnLength = matrix->getVectorLengths();
989    //rhs = new double [numberRows];
990    const double *rowLower = solver->getRowLower();
991    const double *rowUpper = solver->getRowUpper();
992    bool good = true;
993    for (int iRow = 0; iRow < numberRows; iRow++) {
994      if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
995        // SOS
996        rhs[iRow] = -1.0;
997      } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) {
998        good = false;
999      } else if (rowUpper[iRow] < 0.0) {
1000        good = false;
1001      } else if (rowUpper[iRow] < 1.0e10) {
1002        rhs[iRow] = rowUpper[iRow];
1003        if (rhs[iRow] < 0)
1004          sos[iRow] = 0; // can't be SOS
1005      } else {
1006        rhs[iRow] = rowLower[iRow];
1007        if (rhs[iRow] < 0)
1008          sos[iRow] = 0; // can't be SOS
1009      }
1010    }
1011    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1012      if (!columnLength[iColumn])
1013        continue;
1014      if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1015        good = false;
1016      CoinBigIndex j;
1017      int nSOS = 0;
1018      int iSOS = -1;
1019      if (!isHeuristicInteger(solver, iColumn))
1020        good = false;
1021      for (j = columnStart[iColumn];
1022           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1023        if (element[j] < 0.0)
1024          good = false;
1025        int iRow = row[j];
1026        if (rhs[iRow] == -1.0 && sos[iRow] == 'a') {
1027          if (element[j] != 1.0)
1028            good = false;
1029          iSOS = iRow;
1030          nSOS++;
1031        }
1032      }
1033      if (nSOS > 1)
1034        good = false;
1035      else if (!nSOS)
1036        nonSOS++;
1037      sosRow[iColumn] = iSOS;
1038    }
1039    if (!good) {
1040      delete[] sosRow;
1041      delete[] rhs;
1042      delete[] sos;
1043      setWhen(0); // switch off
1044      return 0;
1045    }
1046  } else {
1047    abort(); // not allowed yet
1048  }
1049  const double *solution = solver->getColSolution();
1050  const double *objective = solver->getObjCoefficients();
1051  const double *rowLower = solver->getRowLower();
1052  const double *rowUpper = solver->getRowUpper();
1053  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1054  double primalTolerance;
1055  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1056
1057  numRuns_++;
1058  assert(numberRows == matrix_.getNumRows());
1059  // set up linked list for sets
1060  int *firstGub = new int[numberRows];
1061  int *nextGub = new int[numberColumns];
1062  int iRow, iColumn;
1063  double direction = solver->getObjSense();
1064  double *slackCost = new double[numberRows];
1065  double *modifiedCost = CoinCopyOfArray(objective, numberColumns);
1066  for (int iRow = 0; iRow < numberRows; iRow++) {
1067    slackCost[iRow] = 1.0e30;
1068    firstGub[iRow] = -1;
1069  }
1070  // Take off cost of gub slack
1071  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1072    nextGub[iColumn] = -1;
1073    int iRow = sosRow[iColumn];
1074    if (columnLength[iColumn] == 1 && iRow >= 0) {
1075      // SOS slack
1076      double cost = direction * objective[iColumn];
1077      assert(rhs[iRow] < 0.0);
1078      slackCost[iRow] = CoinMin(slackCost[iRow], cost);
1079    }
1080  }
1081  double offset2 = 0.0;
1082  for (int iRow = 0; iRow < numberRows; iRow++) {
1083    if (sos[iRow] == 'a') {
1084      // row is possible
1085      sos[iRow] = 0;
1086      if (rhs[iRow] < 0.0) {
1087        sos[iRow] = 1;
1088        rhs[iRow] = 1.0;
1089      } else if (rhs[iRow] != rowUpper[iRow]) {
1090        // G row
1091        sos[iRow] = -1;
1092      }
1093      if (slackCost[iRow] == 1.0e30) {
1094        slackCost[iRow] = 0.0;
1095      } else {
1096        offset2 += slackCost[iRow];
1097        sos[iRow] = 2;
1098      }
1099    }
1100  }
1101  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1102    double cost = direction * modifiedCost[iColumn];
1103    CoinBigIndex j;
1104    for (j = columnStart[iColumn];
1105         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1106      int iRow = row[j];
1107      if (sos[iRow] > 0) {
1108        cost -= slackCost[iRow];
1109        if (firstGub[iRow] < 0) {
1110          firstGub[iRow] = iColumn;
1111        } else {
1112          int jColumn = firstGub[iRow];
1113          while (nextGub[jColumn] >= 0)
1114            jColumn = nextGub[jColumn];
1115          nextGub[jColumn] = iColumn;
1116        }
1117        // Only in one sos
1118        break;
1119      }
1120    }
1121    modifiedCost[iColumn] = cost;
1122  }
1123  delete[] slackCost;
1124  double offset;
1125  solver->getDblParam(OsiObjOffset, offset);
1126  double newSolutionValue = -offset + offset2;
1127  int returnCode = 0;
1128
1129  // Get solution array for heuristic solution
1130  double *newSolution = new double[numberColumns];
1131  double *rowActivity = new double[numberRows];
1132  double *contribution = new double[numberColumns];
1133  int *which = new int[numberColumns];
1134  double *newSolution0 = new double[numberColumns];
1135  if ((algorithm_ & (2 | 4)) == 0) {
1136    // get solution as small as possible
1137    for (iColumn = 0; iColumn < numberColumns; iColumn++)
1138      newSolution0[iColumn] = columnLower[iColumn];
1139  } else {
1140    // Get rounded down solution
1141    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1142      double value = solution[iColumn];
1143      // Round down integer
1144      if (fabs(floor(value + 0.5) - value) < integerTolerance) {
1145        value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
1146      } else {
1147        value = CoinMax(floor(value), columnLower[iColumn]);
1148      }
1149      // make sure clean
1150      value = CoinMin(value, columnUpper[iColumn]);
1151      value = CoinMax(value, columnLower[iColumn]);
1152      newSolution0[iColumn] = value;
1153    }
1154  }
1155  double *rowWeight = new double[numberRows];
1156  for (int i = 0; i < numberRows; i++)
1157    rowWeight[i] = 1.0;
1158  double costBias = 0.0;
1159  int nPass = ((algorithm_ & 4) != 0) ? 1 : 10;
1160  for (int iPass = 0; iPass < nPass; iPass++) {
1161    newSolutionValue = -offset + offset2;
1162    memcpy(newSolution, newSolution0, numberColumns * sizeof(double));
1163    // get row activity
1164    memset(rowActivity, 0, numberRows * sizeof(double));
1165    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1166      CoinBigIndex j;
1167      double value = newSolution[iColumn];
1168      for (j = columnStart[iColumn];
1169           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1170        int iRow = row[j];
1171        rowActivity[iRow] += value * element[j];
1172      }
1173    }
1174    if (!rowWeight) {
1175      rowWeight = CoinCopyOfArray(rowActivity, numberRows);
1176    }
1177    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1178      CoinBigIndex j;
1179      double value = newSolution[iColumn];
1180      double cost = modifiedCost[iColumn];
1181      double forSort = 1.0e-24;
1182      bool hasSlack = false;
1183      bool willFit = true;
1184      bool gRow = false;
1185      newSolutionValue += value * cost;
1186      cost += 1.0e-12;
1187      for (j = columnStart[iColumn];
1188           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1189        int iRow = row[j];
1190        int type = sos[iRow];
1191        double gap = rhs[iRow] - rowActivity[iRow] + 1.0e-8;
1192        switch (type) {
1193        case -1:
1194          // G row
1195          gRow = true;
1196#if 0
1197            if (rhs[iRow]>rowWeight[iRow]||(algorithm_&(2|4))!=0)
1198              forSort += element[j];
1199            else
1200              forSort += 0.1*element[j];
1201#else
1202          forSort += rowWeight[iRow] * element[j];
1203#endif
1204          break;
1205        case 0:
1206          // L row
1207          if (gap < element[j]) {
1208            willFit = false;
1209          } else {
1210            forSort += element[j];
1211          }
1212          break;
1213        case 1:
1214          // SOS without slack
1215          if (gap < element[j]) {
1216            willFit = false;
1217          }
1218          break;
1219        case 2:
1220          // SOS with slack
1221          hasSlack = true;
1222          if (gap < element[j]) {
1223            willFit = false;
1224          }
1225          break;
1226        }
1227      }
1228      bool isSlack = hasSlack && (columnLength[iColumn] == 1);
1229      if (forSort < 1.0e-24)
1230        forSort = 1.0e-12;
1231      if ((algorithm_ & 4) != 0 && forSort > 1.0e-24)
1232        forSort = 1.0;
1233      // Use smallest cost if will fit
1234      if (willFit && (hasSlack || gRow) && value == 0.0 && columnUpper[iColumn]) {
1235        if (hasSlack && !gRow) {
1236          if (cost > 1.0e-12) {
1237            forSort = 2.0e30;
1238          } else if (cost == 1.0e-12) {
1239            if (!isSlack)
1240              forSort = 1.0e29;
1241            else
1242              forSort = 1.0e28;
1243          } else {
1244            forSort = cost / forSort;
1245          }
1246        } else {
1247          if (!gRow || true)
1248            forSort = (cost + costBias) / forSort;
1249          else
1250            forSort = 1.0e-12 / forSort;
1251        }
1252      } else {
1253        // put at end
1254        forSort = 1.0e30;
1255      }
1256      which[iColumn] = iColumn;
1257      contribution[iColumn] = forSort;
1258    }
1259    CoinSort_2(contribution, contribution + numberColumns, which);
1260    // Go through columns
1261    int nAdded = 0;
1262    int nSlacks = 0;
1263    for (int jColumn = 0; jColumn < numberColumns; jColumn++) {
1264      if (contribution[jColumn] >= 1.0e30)
1265        break;
1266      int iColumn = which[jColumn];
1267      double value = newSolution[iColumn];
1268      if (value)
1269        continue;
1270      bool possible = true;
1271      CoinBigIndex j;
1272      for (j = columnStart[iColumn];
1273           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1274        int iRow = row[j];
1275        if (sos[iRow] > 0 && rowActivity[iRow]) {
1276          possible = false;
1277        } else {
1278          double gap = rhs[iRow] - rowActivity[iRow] + 1.0e-8;
1279          if (gap < element[j] && sos[iRow] >= 0)
1280            possible = false;
1281        }
1282      }
1283      if (possible) {
1284        //#define REPORT 1
1285#ifdef REPORT
1286        if ((nAdded % 1000) == 0) {
1287          double gap = 0.0;
1288          for (int i = 0; i < numberRows; i++) {
1289            if (rowUpper[i] > 1.0e20)
1290              gap += CoinMax(rowLower[i] - rowActivity[i], 0.0);
1291          }
1292          if (gap)
1293            printf("after %d added gap %g - %d slacks\n",
1294              nAdded, gap, nSlacks);
1295        }
1296#endif
1297        nAdded++;
1298        if (columnLength[iColumn] == 1)
1299          nSlacks++;
1300        // Increase chosen column
1301        newSolution[iColumn] = 1.0;
1302        double cost = modifiedCost[iColumn];
1303        newSolutionValue += cost;
1304        for (CoinBigIndex j = columnStart[iColumn];
1305             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1306          int iRow = row[j];
1307          rowActivity[iRow] += element[j];
1308        }
1309      }
1310    }
1311#ifdef REPORT
1312    {
1313      double under = 0.0;
1314      double over = 0.0;
1315      double gap = 0.0;
1316      int nUnder = 0;
1317      int nOver = 0;
1318      int nGap = 0;
1319      for (iRow = 0; iRow < numberRows; iRow++) {
1320        if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance) {
1321          double value = rowLower[iRow] - rowActivity[iRow];
1322#if REPORT > 1
1323          printf("below on %d is %g - activity %g lower %g\n",
1324            iRow, value, rowActivity[iRow], rowLower[iRow]);
1325#endif
1326          under += value;
1327          nUnder++;
1328        } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0 * primalTolerance) {
1329          double value = rowActivity[iRow] - rowUpper[iRow];
1330#if REPORT > 1
1331          printf("above on %d is %g - activity %g upper %g\n",
1332            iRow, value, rowActivity[iRow], rowUpper[iRow]);
1333#endif
1334          over += value;
1335          nOver++;
1336        } else {
1337          double value = rowActivity[iRow] - rowLower[iRow];
1338          if (value && value < 1.0e20) {
1339#if REPORT > 1
1340            printf("gap on %d is %g - activity %g lower %g\n",
1341              iRow, value, rowActivity[iRow], rowLower[iRow]);
1342#endif
1343            gap += value;
1344            nGap++;
1345          }
1346        }
1347      }
1348      printf("final under %g (%d) - over %g (%d) - free %g (%d) - %d added - solvalue %g\n",
1349        under, nUnder, over, nOver, gap, nGap, nAdded, newSolutionValue);
1350    }
1351#endif
1352    double gap = 0.0;
1353    double over = 0.0;
1354    int nL = 0;
1355    int nG = 0;
1356    int nUnder = 0;
1357    for (iRow = 0; iRow < numberRows; iRow++) {
1358      if (rowLower[iRow] < -1.0e20)
1359        nL++;
1360      if (rowUpper[iRow] > 1.0e20)
1361        nG++;
1362      if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance) {
1363        gap += rowLower[iRow] - rowActivity[iRow];
1364        nUnder++;
1365        rowWeight[iRow] *= 1.1;
1366      } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0 * primalTolerance) {
1367        gap += rowActivity[iRow] - rowUpper[iRow];
1368      } else {
1369        over += rowActivity[iRow] - rowLower[iRow];
1370        //rowWeight[iRow] *= 0.9;
1371      }
1372    }
1373    if (nG && !nL) {
1374      // can we fix
1375      // get list of columns which can go down without making
1376      // things much worse
1377      int nPossible = 0;
1378      int nEasyDown = 0;
1379      int nSlackDown = 0;
1380      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1381        if (newSolution[iColumn] && columnUpper[iColumn] > columnLower[iColumn]) {
1382          bool canGoDown = true;
1383          bool under = false;
1384          int iSos = -1;
1385          for (CoinBigIndex j = columnStart[iColumn];
1386               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1387            int iRow = row[j];
1388            if (sos[iRow] < 0) {
1389              double over = rowActivity[iRow] - rowLower[iRow];
1390              if (over >= 0.0 && element[j] > over + 1.0e-12) {
1391                canGoDown = false;
1392                break;
1393              } else if (over < 0.0) {
1394                under = true;
1395              }
1396            } else {
1397              iSos = iRow;
1398            }
1399          }
1400          if (canGoDown) {
1401            if (!under) {
1402              if (iSos >= 0) {
1403                // find cheapest
1404                double cheapest = modifiedCost[iColumn];
1405                int iCheapest = -1;
1406                int jColumn = firstGub[iSos];
1407                assert(jColumn >= 0);
1408                while (jColumn >= 0) {
1409                  if (modifiedCost[jColumn] < cheapest) {
1410                    cheapest = modifiedCost[jColumn];
1411                    iCheapest = jColumn;
1412                  }
1413                  jColumn = nextGub[jColumn];
1414                }
1415                if (iCheapest >= 0) {
1416                  // Decrease column
1417                  newSolution[iColumn] = 0.0;
1418                  newSolutionValue -= modifiedCost[iColumn];
1419                  for (CoinBigIndex j = columnStart[iColumn];
1420                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1421                    int iRow = row[j];
1422                    rowActivity[iRow] -= element[j];
1423                  }
1424                  // Increase chosen column
1425                  newSolution[iCheapest] = 1.0;
1426                  newSolutionValue += modifiedCost[iCheapest];
1427                  for (CoinBigIndex j = columnStart[iCheapest];
1428                       j < columnStart[iCheapest] + columnLength[iCheapest]; j++) {
1429                    int iRow = row[j];
1430                    rowActivity[iRow] += element[j];
1431                  }
1432                  nEasyDown++;
1433                  if (columnLength[iColumn] > 1) {
1434                    //printf("%d is easy down\n",iColumn);
1435                  } else {
1436                    nSlackDown++;
1437                  }
1438                }
1439              } else if (modifiedCost[iColumn] > 0.0) {
1440                // easy down
1441                // Decrease column
1442                newSolution[iColumn] = 0.0;
1443                newSolutionValue -= modifiedCost[iColumn];
1444                for (CoinBigIndex j = columnStart[iColumn];
1445                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1446                  int iRow = row[j];
1447                  rowActivity[iRow] -= element[j];
1448                }
1449                nEasyDown++;
1450              }
1451            } else {
1452              which[nPossible++] = iColumn;
1453            }
1454          }
1455        }
1456      }
1457#ifdef REPORT
1458      printf("%d possible down, %d easy down of which %d are slacks\n",
1459        nPossible, nEasyDown, nSlackDown);
1460#endif
1461      double *needed = new double[numberRows];
1462      for (int i = 0; i < numberRows; i++) {
1463        double value = rowLower[i] - rowActivity[i];
1464        if (value < 1.0e-8)
1465          value = 0.0;
1466        needed[i] = value;
1467      }
1468      if (gap && /*nUnder==1 &&*/ nonSOS) {
1469        double *weight = new double[numberColumns];
1470        int *sort = new int[numberColumns];
1471        // look at ones not in set
1472        int nPossible = 0;
1473        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1474          if (!newSolution[iColumn] && columnUpper[iColumn] > columnLower[iColumn]) {
1475            int iSos = -1;
1476            double value = 0.0;
1477            for (CoinBigIndex j = columnStart[iColumn];
1478                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1479              int iRow = row[j];
1480              if (sos[iRow] < 0) {
1481                if (needed[iRow])
1482                  value += CoinMin(element[j] / needed[iRow], 1.0);
1483              } else {
1484                iSos = iRow;
1485              }
1486            }
1487            if (value && iSos < 0) {
1488              weight[nPossible] = -value;
1489              sort[nPossible++] = iColumn;
1490            }
1491          }
1492        }
1493        CoinSort_2(weight, weight + nPossible, sort);
1494        for (int i = 0; i < nPossible; i++) {
1495          int iColumn = sort[i];
1496          double helps = 0.0;
1497          for (CoinBigIndex j = columnStart[iColumn];
1498               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1499            int iRow = row[j];
1500            if (needed[iRow])
1501              helps += CoinMin(needed[iRow], element[j]);
1502          }
1503          if (helps) {
1504            newSolution[iColumn] = 1.0;
1505            newSolutionValue += modifiedCost[iColumn];
1506            for (CoinBigIndex j = columnStart[iColumn];
1507                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1508              int iRow = row[j];
1509              rowActivity[iRow] += element[j];
1510              if (needed[iRow]) {
1511                needed[iRow] -= element[j];
1512                if (needed[iRow] < 1.0e-8)
1513                  needed[iRow] = 0.0;
1514              }
1515            }
1516            gap -= helps;
1517#ifdef REPORT
1518            {
1519              double gap2 = 0.0;
1520              for (iRow = 0; iRow < numberRows; iRow++) {
1521                if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance) {
1522                  gap2 += rowLower[iRow] - rowActivity[iRow];
1523                }
1524              }
1525              printf("estimated gap (nonsos) %g - computed %g\n",
1526                gap, gap2);
1527            }
1528#endif
1529            if (gap < 1.0e-12)
1530              break;
1531          }
1532        }
1533        delete[] weight;
1534        delete[] sort;
1535      }
1536      if (gap && nPossible /*&&nUnder==1*/ && true && model_->bestSolution()) {
1537        double *weight = new double[numberColumns];
1538        int *sort = new int[numberColumns];
1539        // look at ones in sets
1540        const double *goodSolution = model_->bestSolution();
1541        int nPossible = 0;
1542        double largestWeight = 0.0;
1543        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1544          if (!newSolution[iColumn] && goodSolution[iColumn] && columnUpper[iColumn] > columnLower[iColumn]) {
1545            int iSos = -1;
1546            double value = 0.0;
1547            for (CoinBigIndex j = columnStart[iColumn];
1548                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1549              int iRow = row[j];
1550              if (sos[iRow] < 0) {
1551                if (needed[iRow])
1552                  value += CoinMin(element[j] / needed[iRow], 1.0);
1553              } else {
1554                iSos = iRow;
1555              }
1556            }
1557            if (value && iSos >= 0) {
1558              // see if value bigger than current
1559              int jColumn = firstGub[iSos];
1560              assert(jColumn >= 0);
1561              while (jColumn >= 0) {
1562                if (newSolution[jColumn])
1563                  break;
1564                jColumn = nextGub[jColumn];
1565              }
1566              assert(jColumn >= 0);
1567              double value2 = 0.0;
1568              for (CoinBigIndex j = columnStart[jColumn];
1569                   j < columnStart[jColumn] + columnLength[jColumn]; j++) {
1570                int iRow = row[j];
1571                if (needed[iRow])
1572                  value2 += CoinMin(element[j] / needed[iRow], 1.0);
1573              }
1574              if (value > value2) {
1575                weight[nPossible] = -(value - value2);
1576                largestWeight = CoinMax(largestWeight, (value - value2));
1577                sort[nPossible++] = iColumn;
1578              }
1579            }
1580          }
1581        }
1582        if (nPossible) {
1583          double *temp = new double[numberRows];
1584          int *which2 = new int[numberRows];
1585          memset(temp, 0, numberRows * sizeof(double));
1586          // modify so ones just more than gap best
1587          if (largestWeight > gap && nUnder == 1) {
1588            double offset = 4 * largestWeight;
1589            for (int i = 0; i < nPossible; i++) {
1590              double value = -weight[i];
1591              if (value > gap - 1.0e-12)
1592                weight[i] = -(offset - (value - gap));
1593            }
1594          }
1595          CoinSort_2(weight, weight + nPossible, sort);
1596          for (int i = 0; i < nPossible; i++) {
1597            int iColumn = sort[i];
1598            int n = 0;
1599            // find jColumn
1600            int iSos = -1;
1601            for (CoinBigIndex j = columnStart[iColumn];
1602                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1603              int iRow = row[j];
1604              temp[iRow] = element[j];
1605              which2[n++] = iRow;
1606              if (sos[iRow] >= 0) {
1607                iSos = iRow;
1608              }
1609            }
1610            int jColumn = firstGub[iSos];
1611            assert(jColumn >= 0);
1612            while (jColumn >= 0) {
1613              if (newSolution[jColumn])
1614                break;
1615              jColumn = nextGub[jColumn];
1616            }
1617            assert(jColumn >= 0);
1618            for (CoinBigIndex j = columnStart[jColumn];
1619                 j < columnStart[jColumn] + columnLength[jColumn]; j++) {
1620              int iRow = row[j];
1621              if (!temp[iRow])
1622                which2[n++] = iRow;
1623              temp[iRow] -= element[j];
1624            }
1625            double helps = 0.0;
1626            for (int i = 0; i < n; i++) {
1627              int iRow = which2[i];
1628              double newValue = rowActivity[iRow] + temp[iRow];
1629              if (temp[iRow] > 1.0e-8) {
1630                if (rowActivity[iRow] < rowLower[iRow] - 1.0e-8) {
1631                  helps += CoinMin(temp[iRow],
1632                    rowLower[iRow] - rowActivity[iRow]);
1633                }
1634              } else if (temp[iRow] < -1.0e-8) {
1635                if (newValue < rowLower[iRow] - 1.0e-12) {
1636                  helps -= CoinMin(-temp[iRow],
1637                    1.0 * (rowLower[iRow] - newValue));
1638                }
1639              }
1640            }
1641            if (helps > 0.0) {
1642              newSolution[iColumn] = 1.0;
1643              newSolution[jColumn] = 0.0;
1644              newSolutionValue += modifiedCost[iColumn] - modifiedCost[jColumn];
1645              for (int i = 0; i < n; i++) {
1646                int iRow = which2[i];
1647                double newValue = rowActivity[iRow] + temp[iRow];
1648                rowActivity[iRow] = newValue;
1649                temp[iRow] = 0.0;
1650              }
1651              gap -= helps;
1652#ifdef REPORT
1653              {
1654                double gap2 = 0.0;
1655                for (iRow = 0; iRow < numberRows; iRow++) {
1656                  if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance) {
1657                    gap2 += rowLower[iRow] - rowActivity[iRow];
1658                  }
1659                }
1660                printf("estimated gap %g - computed %g\n",
1661                  gap, gap2);
1662              }
1663#endif
1664              if (gap < 1.0e-8)
1665                break;
1666            } else {
1667              for (int i = 0; i < n; i++)
1668                temp[which2[i]] = 0.0;
1669            }
1670          }
1671          delete[] which2;
1672          delete[] temp;
1673        }
1674        delete[] weight;
1675        delete[] sort;
1676      }
1677      delete[] needed;
1678    }
1679#ifdef REPORT
1680    {
1681      double gap = 0.0;
1682      double over = 0.0;
1683      for (iRow = 0; iRow < numberRows; iRow++) {
1684        if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance) {
1685          double value = rowLower[iRow] - rowActivity[iRow];
1686#if REPORT > 1
1687          printf("below on %d is %g - activity %g lower %g\n",
1688            iRow, value, rowActivity[iRow], rowLower[iRow]);
1689#endif
1690          gap += value;
1691        } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0 * primalTolerance) {
1692          double value = rowActivity[iRow] - rowUpper[iRow];
1693#if REPORT > 1
1694          printf("above on %d is %g - activity %g upper %g\n",
1695            iRow, value, rowActivity[iRow], rowUpper[iRow]);
1696#endif
1697          gap += value;
1698        } else {
1699          double value = rowActivity[iRow] - rowLower[iRow];
1700          if (value) {
1701#if REPORT > 1
1702            printf("over on %d is %g - activity %g lower %g\n",
1703              iRow, value, rowActivity[iRow], rowLower[iRow]);
1704#endif
1705            over += value;
1706          }
1707        }
1708      }
1709      printf("modified final gap %g - over %g - %d added - solvalue %g\n",
1710        gap, over, nAdded, newSolutionValue);
1711    }
1712#endif
1713    if (!gap) {
1714      break;
1715    } else {
1716      if (iPass == 0) {
1717        costBias = 10.0 * newSolutionValue / static_cast< double >(nAdded);
1718      } else {
1719        costBias *= 10.0;
1720      }
1721    }
1722  }
1723  delete[] newSolution0;
1724  delete[] rowWeight;
1725  delete[] sos;
1726  delete[] firstGub;
1727  delete[] nextGub;
1728  if (newSolutionValue < solutionValue) {
1729    // check feasible
1730    memset(rowActivity, 0, numberRows * sizeof(double));
1731    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1732      CoinBigIndex j;
1733      double value = newSolution[iColumn];
1734      if (value) {
1735        for (j = columnStart[iColumn];
1736             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1737          int iRow = row[j];
1738          rowActivity[iRow] += value * element[j];
1739        }
1740      }
1741    }
1742    // check was approximately feasible
1743    bool feasible = true;
1744    for (iRow = 0; iRow < numberRows; iRow++) {
1745      if (rowActivity[iRow] < rowLower[iRow]) {
1746        if (rowActivity[iRow] < rowLower[iRow] - 10.0 * primalTolerance)
1747          feasible = false;
1748      } else if (rowActivity[iRow] > rowUpper[iRow]) {
1749        if (rowActivity[iRow] > rowUpper[iRow] + 10.0 * primalTolerance)
1750          feasible = false;
1751      }
1752    }
1753    if (feasible) {
1754      // new solution
1755      memcpy(betterSolution, newSolution, numberColumns * sizeof(double));
1756      solutionValue = newSolutionValue;
1757      //printf("** Solution of %g found by rounding\n",newSolutionValue);
1758      returnCode = 1;
1759    } else {
1760      // Can easily happen
1761      //printf("Debug CbcHeuristicGreedySOS giving bad solution\n");
1762    }
1763  }
1764  delete[] sosRow;
1765  delete[] newSolution;
1766  delete[] rowActivity;
1767  delete[] modifiedCost;
1768  delete[] contribution;
1769  delete[] which;
1770  delete[] rhs;
1771  return returnCode;
1772}
1773// update model
1774void CbcHeuristicGreedySOS::setModel(CbcModel *model)
1775{
1776  delete[] originalRhs_;
1777  gutsOfConstructor(model);
1778  validate();
1779}
1780// Resets stuff if model changes
1781void CbcHeuristicGreedySOS::resetModel(CbcModel *model)
1782{
1783  delete[] originalRhs_;
1784  gutsOfConstructor(model);
1785}
1786// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1787void CbcHeuristicGreedySOS::validate()
1788{
1789  if (model_ && when() < 10) {
1790    if (model_->numberIntegers() != model_->numberObjects() && (model_->numberObjects() || (model_->specialOptions() & 1024) == 0)) {
1791      int numberOdd = 0;
1792      for (int i = 0; i < model_->numberObjects(); i++) {
1793        if (!model_->object(i)->canDoHeuristics())
1794          numberOdd++;
1795      }
1796      if (numberOdd)
1797        setWhen(0);
1798    }
1799    // Only works if coefficients positive and all rows L/G or SOS
1800    OsiSolverInterface *solver = model_->solver();
1801    const double *columnUpper = solver->getColUpper();
1802    const double *columnLower = solver->getColLower();
1803    const double *rowLower = solver->getRowLower();
1804    const double *rowUpper = solver->getRowUpper();
1805
1806    int numberRows = solver->getNumRows();
1807    // Column copy
1808    const double *element = matrix_.getElements();
1809    const int *row = matrix_.getIndices();
1810    const CoinBigIndex *columnStart = matrix_.getVectorStarts();
1811    const int *columnLength = matrix_.getVectorLengths();
1812    bool good = true;
1813    assert(originalRhs_);
1814    for (int iRow = 0; iRow < numberRows; iRow++) {
1815      if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1816        // SOS
1817        originalRhs_[iRow] = -1.0;
1818      } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) {
1819        good = false;
1820      } else if (rowUpper[iRow] < 0.0) {
1821        good = false;
1822      } else if (rowUpper[iRow] < 1.0e10) {
1823        originalRhs_[iRow] = rowUpper[iRow];
1824      } else {
1825        originalRhs_[iRow] = rowLower[iRow];
1826      }
1827    }
1828    int numberColumns = solver->getNumCols();
1829    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1830      if (!columnLength[iColumn])
1831        continue;
1832      if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1833        good = false;
1834      CoinBigIndex j;
1835      int nSOS = 0;
1836      if (!isHeuristicInteger(solver, iColumn))
1837        good = false;
1838      for (j = columnStart[iColumn];
1839           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1840        if (element[j] < 0.0)
1841          good = false;
1842        int iRow = row[j];
1843        if (originalRhs_[iRow] == -1.0) {
1844          if (element[j] != 1.0)
1845            good = false;
1846          nSOS++;
1847        }
1848      }
1849      if (nSOS > 1)
1850        good = false;
1851    }
1852    if (!good)
1853      setWhen(0); // switch off
1854  }
1855}
1856
1857/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1858*/
Note: See TracBrowser for help on using the repository browser.