source: trunk/Cbc/src/CbcHeuristicGreedy.cpp @ 1589

Last change on this file since 1589 was 1589, checked in by forrest, 9 years ago

allow SOS greedy with 0-1

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.1 KB
Line 
1/* $Id: CbcHeuristicGreedy.cpp 1589 2011-01-28 09:16:51Z forrest $ */
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
54CbcHeuristicGreedyCover::gutsOfConstructor(CbcModel * model)
55{
56    model_ = model;
57    // Get a copy of original matrix
58    assert(model->solver());
59    if (model->solver()->getNumRows()) {
60        matrix_ = *model->solver()->getMatrixByCol();
61    }
62    originalNumberRows_ = model->solver()->getNumRows();
63}
64// Create C++ lines to get to current state
65void
66CbcHeuristicGreedyCover::generateCpp( FILE * fp)
67{
68    CbcHeuristicGreedyCover other;
69    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
70    fprintf(fp, "3  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n");
71    CbcHeuristic::generateCpp(fp, "heuristicGreedyCover");
72    if (algorithm_ != other.algorithm_)
73        fprintf(fp, "3  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
74    else
75        fprintf(fp, "4  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
76    if (numberTimes_ != other.numberTimes_)
77        fprintf(fp, "3  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
78    else
79        fprintf(fp, "4  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
80    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyCover);\n");
81}
82
83// Copy constructor
84CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(const CbcHeuristicGreedyCover & rhs)
85        :
86        CbcHeuristic(rhs),
87        matrix_(rhs.matrix_),
88        originalNumberRows_(rhs.originalNumberRows_),
89        algorithm_(rhs.algorithm_),
90        numberTimes_(rhs.numberTimes_)
91{
92}
93
94// Assignment operator
95CbcHeuristicGreedyCover &
96CbcHeuristicGreedyCover::operator=( const CbcHeuristicGreedyCover & rhs)
97{
98    if (this != &rhs) {
99        CbcHeuristic::operator=(rhs);
100        matrix_ = rhs.matrix_;
101        originalNumberRows_ = rhs.originalNumberRows_;
102        algorithm_ = rhs.algorithm_;
103        numberTimes_ = rhs.numberTimes_;
104    }
105    return *this;
106}
107// Returns 1 if solution, 0 if not
108int
109CbcHeuristicGreedyCover::solution(double & solutionValue,
110                                  double * betterSolution)
111{
112    numCouldRun_++;
113    if (!model_)
114        return 0;
115    // See if to do
116    if (!when() || (when() == 1 && model_->phase() != 1))
117        return 0; // switched off
118    if (model_->getNodeCount() > numberTimes_)
119        return 0;
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 (solver->isInteger(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 (solver->isInteger(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 (solver->isInteger(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 &&
284                                element[j]*step + rowActivity[iRow] >= rowLower[iRow]) {
285                            step = (rowLower[iRow] - rowActivity[iRow]) / element[j];;
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 &&
299                                    newActivity >= rowLower[iRow] - 1.0e-12) {
300                                sum += element[j];
301                            }
302                        }
303                        assert (sum > 0.0);
304                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
305                        if (ratio < bestRatio) {
306                            bestRatio = ratio;
307                            bestColumn = iColumn;
308                            bestStepSize = step;
309                        }
310                    }
311                }
312            }
313        }
314        if (bestColumn < 0)
315            break; // we have finished
316        // Increase chosen column
317        newSolution[bestColumn] += bestStepSize;
318        double cost = direction * objective[bestColumn];
319        newSolutionValue += bestStepSize * cost;
320        for (CoinBigIndex j = columnStart[bestColumn];
321                j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
322            int iRow = row[j];
323            rowActivity[iRow] += bestStepSize * element[j];
324        }
325    }
326    delete [] which;
327    if (newSolutionValue < solutionValue) {
328        // check feasible
329        memset(rowActivity, 0, numberRows*sizeof(double));
330        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
331            CoinBigIndex j;
332            double value = newSolution[iColumn];
333            if (value) {
334                for (j = columnStart[iColumn];
335                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
336                    int iRow = row[j];
337                    rowActivity[iRow] += value * element[j];
338                }
339            }
340        }
341        // check was approximately feasible
342        bool feasible = true;
343        for (iRow = 0; iRow < numberRows; iRow++) {
344            if (rowActivity[iRow] < rowLower[iRow]) {
345                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
346                    feasible = false;
347            }
348        }
349        if (feasible) {
350            // new solution
351            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
352            solutionValue = newSolutionValue;
353            //printf("** Solution of %g found by rounding\n",newSolutionValue);
354            returnCode = 1;
355        } else {
356            // Can easily happen
357            //printf("Debug CbcHeuristicGreedyCover giving bad solution\n");
358        }
359    }
360    delete [] newSolution;
361    delete [] rowActivity;
362    return returnCode;
363}
364// update model
365void CbcHeuristicGreedyCover::setModel(CbcModel * model)
366{
367    gutsOfConstructor(model);
368    validate();
369}
370// Resets stuff if model changes
371void
372CbcHeuristicGreedyCover::resetModel(CbcModel * model)
373{
374    gutsOfConstructor(model);
375}
376// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
377void
378CbcHeuristicGreedyCover::validate()
379{
380    if (model_ && when() < 10) {
381        if (model_->numberIntegers() !=
382                model_->numberObjects() && (model_->numberObjects() ||
383                                            (model_->specialOptions()&1024) == 0)) {
384            int numberOdd = 0;
385            for (int i = 0; i < model_->numberObjects(); i++) {
386                if (!model_->object(i)->canDoHeuristics())
387                    numberOdd++;
388            }
389            if (numberOdd)
390                setWhen(0);
391        }
392        // Only works if costs positive, coefficients positive and all rows G
393        OsiSolverInterface * solver = model_->solver();
394        const double * columnLower = solver->getColLower();
395        const double * rowUpper = solver->getRowUpper();
396        const double * objective = solver->getObjCoefficients();
397        double direction = solver->getObjSense();
398
399        int numberRows = solver->getNumRows();
400        // Column copy
401        const double * element = matrix_.getElements();
402        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
403        const int * columnLength = matrix_.getVectorLengths();
404        bool good = true;
405        for (int iRow = 0; iRow < numberRows; iRow++) {
406            if (rowUpper[iRow] < 1.0e30)
407                good = false;
408        }
409        int numberColumns = solver->getNumCols();
410        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
411            if (objective[iColumn]*direction < 0.0)
412                good = false;
413            if (columnLower[iColumn] < 0.0)
414                good = false;
415            CoinBigIndex j;
416            for (j = columnStart[iColumn];
417                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
418                if (element[j] < 0.0)
419                    good = false;
420            }
421        }
422        if (!good)
423            setWhen(0); // switch off
424    }
425}
426// Default Constructor
427CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality()
428        : CbcHeuristic()
429{
430    // matrix  will automatically be empty
431    fraction_ = 1.0; // no branch and bound
432    originalNumberRows_ = 0;
433    algorithm_ = 0;
434    numberTimes_ = 100;
435    whereFrom_ = 1;
436}
437
438// Constructor from model
439CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel & model)
440        : CbcHeuristic(model)
441{
442    // Get a copy of original matrix
443    gutsOfConstructor(&model);
444    fraction_ = 1.0; // no branch and bound
445    algorithm_ = 0;
446    numberTimes_ = 100;
447    whereFrom_ = 1;
448}
449
450// Destructor
451CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality ()
452{
453}
454
455// Clone
456CbcHeuristic *
457CbcHeuristicGreedyEquality::clone() const
458{
459    return new CbcHeuristicGreedyEquality(*this);
460}
461// Guts of constructor from a CbcModel
462void
463CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel * model)
464{
465    model_ = model;
466    // Get a copy of original matrix
467    assert(model->solver());
468    if (model->solver()->getNumRows()) {
469        matrix_ = *model->solver()->getMatrixByCol();
470    }
471    originalNumberRows_ = model->solver()->getNumRows();
472}
473// Create C++ lines to get to current state
474void
475CbcHeuristicGreedyEquality::generateCpp( FILE * fp)
476{
477    CbcHeuristicGreedyEquality other;
478    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
479    fprintf(fp, "3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
480    CbcHeuristic::generateCpp(fp, "heuristicGreedyEquality");
481    if (algorithm_ != other.algorithm_)
482        fprintf(fp, "3  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
483    else
484        fprintf(fp, "4  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
485    if (fraction_ != other.fraction_)
486        fprintf(fp, "3  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
487    else
488        fprintf(fp, "4  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
489    if (numberTimes_ != other.numberTimes_)
490        fprintf(fp, "3  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
491    else
492        fprintf(fp, "4  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
493    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
494}
495
496// Copy constructor
497CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality & rhs)
498        :
499        CbcHeuristic(rhs),
500        matrix_(rhs.matrix_),
501        fraction_(rhs.fraction_),
502        originalNumberRows_(rhs.originalNumberRows_),
503        algorithm_(rhs.algorithm_),
504        numberTimes_(rhs.numberTimes_)
505{
506}
507
508// Assignment operator
509CbcHeuristicGreedyEquality &
510CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality & rhs)
511{
512    if (this != &rhs) {
513        CbcHeuristic::operator=(rhs);
514        matrix_ = rhs.matrix_;
515        fraction_ = rhs.fraction_;
516        originalNumberRows_ = rhs.originalNumberRows_;
517        algorithm_ = rhs.algorithm_;
518        numberTimes_ = rhs.numberTimes_;
519    }
520    return *this;
521}
522// Returns 1 if solution, 0 if not
523int
524CbcHeuristicGreedyEquality::solution(double & solutionValue,
525                                     double * betterSolution)
526{
527    numCouldRun_++;
528    if (!model_)
529        return 0;
530    // See if to do
531    if (!when() || (when() == 1 && model_->phase() != 1))
532        return 0; // switched off
533    if (model_->getNodeCount() > numberTimes_)
534        return 0;
535    // See if at root node
536    bool atRoot = model_->getNodeCount() == 0;
537    int passNumber = model_->getCurrentPassNumber();
538    if (atRoot && passNumber != 1)
539        return 0;
540    OsiSolverInterface * solver = model_->solver();
541    const double * columnLower = solver->getColLower();
542    const double * columnUpper = solver->getColUpper();
543    // And original upper bounds in case we want to use them
544    const double * originalUpper = model_->continuousSolver()->getColUpper();
545    // But not if algorithm says so
546    if ((algorithm_ % 10) == 0)
547        originalUpper = columnUpper;
548    const double * rowLower = solver->getRowLower();
549    const double * rowUpper = solver->getRowUpper();
550    const double * solution = solver->getColSolution();
551    const double * objective = solver->getObjCoefficients();
552    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
553    double primalTolerance;
554    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
555
556    // This is number of rows when matrix was passed in
557    int numberRows = originalNumberRows_;
558    if (!numberRows)
559        return 0; // switched off
560    numRuns_++;
561
562    assert (numberRows == matrix_.getNumRows());
563    int iRow, iColumn;
564    double direction = solver->getObjSense();
565    double offset;
566    solver->getDblParam(OsiObjOffset, offset);
567    double newSolutionValue = -offset;
568    int returnCode = 0;
569
570    // Column copy
571    const double * element = matrix_.getElements();
572    const int * row = matrix_.getIndices();
573    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
574    const int * columnLength = matrix_.getVectorLengths();
575
576    // Get solution array for heuristic solution
577    int numberColumns = solver->getNumCols();
578    double * newSolution = new double [numberColumns];
579    double * rowActivity = new double[numberRows];
580    memset(rowActivity, 0, numberRows*sizeof(double));
581    double rhsNeeded = 0;
582    for (iRow = 0; iRow < numberRows; iRow++)
583        rhsNeeded += rowUpper[iRow];
584    rhsNeeded *= fraction_;
585    bool allOnes = true;
586    // Get rounded down solution
587    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
588        CoinBigIndex j;
589        double value = solution[iColumn];
590        if (solver->isInteger(iColumn)) {
591            // Round down integer
592            if (fabs(floor(value + 0.5) - value) < integerTolerance) {
593                value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
594            } else {
595                value = CoinMax(floor(value), columnLower[iColumn]);
596            }
597        }
598        // make sure clean
599        value = CoinMin(value, columnUpper[iColumn]);
600        value = CoinMax(value, columnLower[iColumn]);
601        newSolution[iColumn] = value;
602        double cost = direction * objective[iColumn];
603        newSolutionValue += value * cost;
604        for (j = columnStart[iColumn];
605                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
606            int iRow = row[j];
607            rowActivity[iRow] += value * element[j];
608            rhsNeeded -= value * element[j];
609            if (element[j] != 1.0)
610                allOnes = false;
611        }
612    }
613    // See if we round up
614    bool roundup = ((algorithm_ % 100) != 0);
615    if (roundup && allOnes) {
616        // Get rounded up solution
617        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
618            CoinBigIndex j;
619            double value = solution[iColumn];
620            if (solver->isInteger(iColumn)) {
621                // but round up if no activity
622                if (roundup && value >= 0.6 && !newSolution[iColumn]) {
623                    bool choose = true;
624                    for (j = columnStart[iColumn];
625                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
626                        int iRow = row[j];
627                        if (rowActivity[iRow]) {
628                            choose = false;
629                            break;
630                        }
631                    }
632                    if (choose) {
633                        newSolution[iColumn] = 1.0;
634                        double cost = direction * objective[iColumn];
635                        newSolutionValue += cost;
636                        for (j = columnStart[iColumn];
637                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
638                            int iRow = row[j];
639                            rowActivity[iRow] += 1.0;
640                            rhsNeeded -= 1.0;
641                        }
642                    }
643                }
644            }
645        }
646    }
647    // Get initial list
648    int * which = new int [numberColumns];
649    for (iColumn = 0; iColumn < numberColumns; iColumn++)
650        which[iColumn] = iColumn;
651    int numberLook = numberColumns;
652    // See if we want to perturb more
653    double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
654    // Keep going round until a solution
655    while (true) {
656        // Get column with best ratio
657        int bestColumn = -1;
658        double bestRatio = COIN_DBL_MAX;
659        double bestStepSize = 0.0;
660        int newNumber = 0;
661        for (int jColumn = 0; jColumn < numberLook; jColumn++) {
662            int iColumn = which[jColumn];
663            CoinBigIndex j;
664            double value = newSolution[iColumn];
665            double cost = direction * objective[iColumn];
666            if (solver->isInteger(iColumn)) {
667                // use current upper or original upper
668                if (value + 0.9999 < originalUpper[iColumn]) {
669                    double movement = 1.0;
670                    double sum = 0.0;
671                    for (j = columnStart[iColumn];
672                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
673                        int iRow = row[j];
674                        double gap = rowUpper[iRow] - rowActivity[iRow];
675                        double elementValue = allOnes ? 1.0 : element[j];
676                        sum += elementValue;
677                        if (movement*elementValue > gap) {
678                            movement = gap / elementValue;
679                        }
680                    }
681                    if (movement > 0.999999) {
682                        // add to next time
683                        which[newNumber++] = iColumn;
684                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
685                        // If at root
686                        if (atRoot) {
687                            if (fraction_ == 1.0)
688                                ratio = iColumn; // choose first
689                            else
690                                ratio = - solution[iColumn]; // choose largest
691                        }
692                        if (ratio < bestRatio) {
693                            bestRatio = ratio;
694                            bestColumn = iColumn;
695                            bestStepSize = 1.0;
696                        }
697                    }
698                }
699            } else {
700                // continuous
701                if (value < columnUpper[iColumn]) {
702                    double movement = 1.0e50;
703                    double sum = 0.0;
704                    for (j = columnStart[iColumn];
705                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
706                        int iRow = row[j];
707                        if (element[j]*movement + rowActivity[iRow] > rowUpper[iRow]) {
708                            movement = (rowUpper[iRow] - rowActivity[iRow]) / element[j];;
709                        }
710                        sum += element[j];
711                    }
712                    // now ratio
713                    if (movement > 1.0e-7) {
714                        // add to next time
715                        which[newNumber++] = iColumn;
716                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
717                        if (ratio < bestRatio) {
718                            bestRatio = ratio;
719                            bestColumn = iColumn;
720                            bestStepSize = movement;
721                        }
722                    }
723                }
724            }
725        }
726        if (bestColumn < 0)
727            break; // we have finished
728        // Increase chosen column
729        newSolution[bestColumn] += bestStepSize;
730        double cost = direction * objective[bestColumn];
731        newSolutionValue += bestStepSize * cost;
732        for (CoinBigIndex j = columnStart[bestColumn];
733                j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
734            int iRow = row[j];
735            rowActivity[iRow] += bestStepSize * element[j];
736            rhsNeeded -= bestStepSize * element[j];
737        }
738        if (rhsNeeded < 1.0e-8)
739            break;
740    }
741    delete [] which;
742    if (fraction_ < 1.0 && rhsNeeded < 1.0e-8 && newSolutionValue < solutionValue) {
743        // do branch and cut
744        // fix all nonzero
745        OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
746        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
747            if (newSolver->isInteger(iColumn))
748                newSolver->setColLower(iColumn, newSolution[iColumn]);
749        }
750        int returnCode = smallBranchAndBound(newSolver, 200, newSolution, newSolutionValue,
751                                             solutionValue, "CbcHeuristicGreedy");
752        if (returnCode < 0)
753            returnCode = 0; // returned on size
754        if ((returnCode&2) != 0) {
755            // could add cut
756            returnCode &= ~2;
757        }
758        rhsNeeded = 1.0 - returnCode;
759        delete newSolver;
760    }
761    if (newSolutionValue < solutionValue && rhsNeeded < 1.0e-8) {
762        // check feasible
763        memset(rowActivity, 0, numberRows*sizeof(double));
764        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
765            CoinBigIndex j;
766            double value = newSolution[iColumn];
767            if (value) {
768                for (j = columnStart[iColumn];
769                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
770                    int iRow = row[j];
771                    rowActivity[iRow] += value * element[j];
772                }
773            }
774        }
775        // check was approximately feasible
776        bool feasible = true;
777        for (iRow = 0; iRow < numberRows; iRow++) {
778            if (rowActivity[iRow] < rowLower[iRow]) {
779                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
780                    feasible = false;
781            }
782        }
783        if (feasible) {
784            // new solution
785            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
786            solutionValue = newSolutionValue;
787            returnCode = 1;
788        }
789    }
790    delete [] newSolution;
791    delete [] rowActivity;
792    if (atRoot && fraction_ == 1.0) {
793        // try quick search
794        fraction_ = 0.4;
795        int newCode = this->solution(solutionValue, betterSolution);
796        if (newCode)
797            returnCode = 1;
798        fraction_ = 1.0;
799    }
800    return returnCode;
801}
802// update model
803void CbcHeuristicGreedyEquality::setModel(CbcModel * model)
804{
805    gutsOfConstructor(model);
806    validate();
807}
808// Resets stuff if model changes
809void
810CbcHeuristicGreedyEquality::resetModel(CbcModel * model)
811{
812    gutsOfConstructor(model);
813}
814// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
815void
816CbcHeuristicGreedyEquality::validate()
817{
818    if (model_ && when() < 10) {
819        if (model_->numberIntegers() !=
820                model_->numberObjects())
821            setWhen(0);
822        // Only works if costs positive, coefficients positive and all rows E or L
823        // And if values are integer
824        OsiSolverInterface * solver = model_->solver();
825        const double * columnLower = solver->getColLower();
826        const double * rowUpper = solver->getRowUpper();
827        const double * rowLower = solver->getRowLower();
828        const double * objective = solver->getObjCoefficients();
829        double direction = solver->getObjSense();
830
831        int numberRows = solver->getNumRows();
832        // Column copy
833        const double * element = matrix_.getElements();
834        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
835        const int * columnLength = matrix_.getVectorLengths();
836        bool good = true;
837        for (int iRow = 0; iRow < numberRows; iRow++) {
838            if (rowUpper[iRow] > 1.0e30)
839                good = false;
840            if (rowLower[iRow] > 0.0 && rowLower[iRow] != rowUpper[iRow])
841                good = false;
842            if (floor(rowUpper[iRow] + 0.5) != rowUpper[iRow])
843                good = false;
844        }
845        int numberColumns = solver->getNumCols();
846        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
847            if (objective[iColumn]*direction < 0.0)
848                good = false;
849            if (columnLower[iColumn] < 0.0)
850                good = false;
851            CoinBigIndex j;
852            for (j = columnStart[iColumn];
853                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
854                if (element[j] < 0.0)
855                    good = false;
856                if (floor(element[j] + 0.5) != element[j])
857                    good = false;
858            }
859        }
860        if (!good)
861            setWhen(0); // switch off
862    }
863}
864
865// Default Constructor
866CbcHeuristicGreedySOS::CbcHeuristicGreedySOS()
867        : CbcHeuristic()
868{
869    originalRhs_ = NULL;
870    // matrix  will automatically be empty
871    originalNumberRows_ = 0;
872    algorithm_ = 0;
873    numberTimes_ = 100;
874}
875
876// Constructor from model
877CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(CbcModel & model)
878        : CbcHeuristic(model)
879{
880    gutsOfConstructor(&model);
881    algorithm_ = 2;
882    numberTimes_ = 100;
883    whereFrom_ = 1;
884}
885
886// Destructor
887CbcHeuristicGreedySOS::~CbcHeuristicGreedySOS ()
888{
889  delete [] originalRhs_;
890}
891
892// Clone
893CbcHeuristic *
894CbcHeuristicGreedySOS::clone() const
895{
896    return new CbcHeuristicGreedySOS(*this);
897}
898// Guts of constructor from a CbcModel
899void
900CbcHeuristicGreedySOS::gutsOfConstructor(CbcModel * model)
901{
902    model_ = model;
903    // Get a copy of original matrix
904    assert(model->solver());
905    if (model->solver()->getNumRows()) {
906        matrix_ = *model->solver()->getMatrixByCol();
907    }
908    originalNumberRows_ = model->solver()->getNumRows();
909    originalRhs_ = new double [originalNumberRows_];
910}
911// Create C++ lines to get to current state
912void
913CbcHeuristicGreedySOS::generateCpp( FILE * fp)
914{
915    CbcHeuristicGreedySOS other;
916    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
917    fprintf(fp, "3  CbcHeuristicGreedySOS heuristicGreedySOS(*cbcModel);\n");
918    CbcHeuristic::generateCpp(fp, "heuristicGreedySOS");
919    if (algorithm_ != other.algorithm_)
920        fprintf(fp, "3  heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_);
921    else
922        fprintf(fp, "4  heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_);
923    if (numberTimes_ != other.numberTimes_)
924        fprintf(fp, "3  heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_);
925    else
926        fprintf(fp, "4  heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_);
927    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedySOS);\n");
928}
929
930// Copy constructor
931CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(const CbcHeuristicGreedySOS & rhs)
932        :
933        CbcHeuristic(rhs),
934        matrix_(rhs.matrix_),
935        originalNumberRows_(rhs.originalNumberRows_),
936        algorithm_(rhs.algorithm_),
937        numberTimes_(rhs.numberTimes_)
938{
939  originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_);
940}
941
942// Assignment operator
943CbcHeuristicGreedySOS &
944CbcHeuristicGreedySOS::operator=( const CbcHeuristicGreedySOS & rhs)
945{
946    if (this != &rhs) {
947        CbcHeuristic::operator=(rhs);
948        matrix_ = rhs.matrix_;
949        originalNumberRows_ = rhs.originalNumberRows_;
950        algorithm_ = rhs.algorithm_;
951        numberTimes_ = rhs.numberTimes_;
952        delete [] originalRhs_;
953        originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_);
954    }
955    return *this;
956}
957// Returns 1 if solution, 0 if not
958int
959CbcHeuristicGreedySOS::solution(double & solutionValue,
960                                  double * betterSolution)
961{
962    numCouldRun_++;
963    if (!model_)
964        return 0;
965    // See if to do
966    if (!when() || (when() == 1 && model_->phase() != 1))
967        return 0; // switched off
968    if (model_->getNodeCount() > numberTimes_)
969        return 0;
970    // See if at root node
971    bool atRoot = model_->getNodeCount() == 0;
972    int passNumber = model_->getCurrentPassNumber();
973    if (atRoot && passNumber != 1)
974        return 0;
975    OsiSolverInterface * solver = model_->solver();
976    int numberColumns = solver->getNumCols();
977    // This is number of rows when matrix was passed in
978    int numberRows = originalNumberRows_;
979    if (!numberRows)
980        return 0; // switched off
981
982    const double * columnLower = solver->getColLower();
983    const double * columnUpper = solver->getColUpper();
984    // modified rhs
985    double * rhs = CoinCopyOfArray(originalRhs_,numberRows);
986    // Column copy
987    const double * element = matrix_.getElements();
988    const int * row = matrix_.getIndices();
989    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
990    const int * columnLength = matrix_.getVectorLengths();
991    int * sosRow = new int [numberColumns];
992    // If bit set then use current
993    if ((algorithm_&1)!=0) {
994      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
995      element = matrix->getElements();
996      row = matrix->getIndices();
997      columnStart = matrix->getVectorStarts();
998      columnLength = matrix->getVectorLengths();
999      //rhs = new double [numberRows];
1000      const double * rowLower = solver->getRowLower();
1001      const double * rowUpper = solver->getRowUpper();
1002      bool good = true;
1003      for (int iRow = 0; iRow < numberRows; iRow++) {
1004        if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1005          // SOS
1006          rhs[iRow]=-1.0;
1007        } else if (rowLower[iRow] > 0.0) {
1008          good = false;
1009        } else if (rowUpper[iRow] < 0.0) {
1010          good = false;
1011        } else {
1012          rhs[iRow]=rowUpper[iRow];
1013        }
1014      }
1015      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1016        if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1017          good = false;
1018        CoinBigIndex j;
1019        int nSOS=0;
1020        int iSOS=-1;
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) {
1027            if (element[j] != 1.0)
1028              good = false;
1029            iSOS=iRow;
1030            nSOS++;
1031          }
1032        }
1033        if (nSOS>1||!solver->isBinary(iColumn))
1034          good = false;
1035        sosRow[iColumn] = iSOS;
1036      }
1037      if (!good) {
1038        delete [] rhs;
1039        delete [] sosRow;
1040        setWhen(0); // switch off
1041        return 0;
1042      }
1043    } else {
1044      abort(); // not allowed yet
1045    }
1046    const double * solution = solver->getColSolution();
1047    const double * objective = solver->getObjCoefficients();
1048    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1049    double primalTolerance;
1050    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1051
1052    numRuns_++;
1053    assert (numberRows == matrix_.getNumRows());
1054    int iRow, iColumn;
1055    double direction = solver->getObjSense();
1056    double * slackCost = new double [numberRows];
1057    double * modifiedCost = CoinCopyOfArray(objective,numberColumns);
1058    for (int iRow = 0;iRow < numberRows; iRow++)
1059      slackCost[iRow]=1.0e30;
1060    // Take off cost of gub slack
1061    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1062      int iRow = sosRow[iColumn];
1063      if (columnLength[iColumn] == 1&&iRow>=0) {
1064        // SOS slack
1065        double cost = direction*objective[iColumn];
1066        assert (rhs[iRow]<0.0);
1067        slackCost[iRow]=CoinMin(slackCost[iRow],cost);
1068      }
1069    }
1070    double offset2 = 0.0;
1071    char * sos = new char [numberRows];
1072    for (int iRow = 0;iRow < numberRows; iRow++) {
1073      sos[iRow]=0;
1074      if (rhs[iRow]<0.0) {
1075        sos[iRow]=1;
1076        rhs[iRow]=1.0;
1077      }
1078      if( slackCost[iRow] == 1.0e30) {
1079        slackCost[iRow]=0.0;
1080      } else {
1081        offset2 += slackCost[iRow];
1082        sos[iRow] = 2;
1083      }
1084    }
1085    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1086      double cost = direction * modifiedCost[iColumn];
1087      CoinBigIndex j;
1088      for (j = columnStart[iColumn];
1089           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1090        int iRow = row[j];
1091        if (sos[iRow]) {
1092          cost -= slackCost[iRow];
1093        }
1094      }
1095      modifiedCost[iColumn] = cost;
1096    }
1097    delete [] slackCost;
1098    double offset;
1099    solver->getDblParam(OsiObjOffset, offset);
1100    double newSolutionValue = -offset+offset2;
1101    int returnCode = 0;
1102
1103
1104    // Get solution array for heuristic solution
1105    double * newSolution = new double [numberColumns];
1106    double * rowActivity = new double[numberRows];
1107    memset(rowActivity, 0, numberRows*sizeof(double));
1108    if ((algorithm_&(2|4))==0) {
1109      // get solution as small as possible
1110      for (iColumn = 0; iColumn < numberColumns; iColumn++) 
1111        newSolution[iColumn] = columnLower[iColumn];
1112    } else {
1113      // Get rounded down solution
1114      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1115        double value = solution[iColumn];
1116        // Round down integer
1117        if (fabs(floor(value + 0.5) - value) < integerTolerance) {
1118          value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
1119        } else {
1120          value = CoinMax(floor(value), columnLower[iColumn]);
1121        }
1122        // make sure clean
1123        value = CoinMin(value, columnUpper[iColumn]);
1124        value = CoinMax(value, columnLower[iColumn]);
1125        newSolution[iColumn] = value;
1126      }
1127    }
1128    // get row activity
1129    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1130      CoinBigIndex j;
1131      double value = newSolution[iColumn];
1132      for (j = columnStart[iColumn];
1133           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1134        int iRow = row[j];
1135        rowActivity[iRow] += value * element[j];
1136      }
1137    }
1138    double * contribution = new double [numberColumns];
1139    int * which = new int [numberColumns];
1140    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1141      CoinBigIndex j;
1142      double value = newSolution[iColumn];
1143      double cost =  modifiedCost[iColumn];
1144      double forSort = 0.0;
1145      bool hasSlack=false;
1146      bool willFit=true;
1147      newSolutionValue += value * cost;
1148      for (j = columnStart[iColumn];
1149           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1150        int iRow = row[j];
1151        if (sos[iRow] == 2) 
1152          hasSlack = true;
1153        forSort += element[j];
1154        double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1155        if (gap<element[j])
1156          willFit = false;
1157      }
1158      bool isSlack = hasSlack && (columnLength[iColumn]==1);
1159      if ((algorithm_&4)!=0) 
1160        forSort=1.0;
1161      // Use smallest cost if will fit
1162      if (willFit /*&& hasSlack*/ && 
1163          value == 0.0 && columnUpper[iColumn]) {
1164        if (hasSlack) {
1165          if (cost>0.0) {
1166            forSort = 2.0e30;
1167          } else if (cost==0.0) {
1168            if (!isSlack)
1169              forSort = 1.0e29;
1170            else
1171              forSort = 1.0e28;
1172          } else {
1173            forSort = cost/forSort;
1174          }
1175        } else {
1176          forSort = cost/forSort;
1177        }
1178      } else {
1179        // put at end
1180        forSort = 1.0e30;
1181      }
1182      which[iColumn]=iColumn;
1183      contribution[iColumn]= forSort;
1184    }
1185    CoinSort_2(contribution,contribution+numberColumns,which);
1186    // Go through columns
1187    for (int jColumn = 0; jColumn < numberColumns; jColumn++) {
1188      int iColumn = which[jColumn];
1189      double value = newSolution[iColumn];
1190      if (value)
1191        continue;
1192      bool possible = true;
1193      CoinBigIndex j;
1194      for (j = columnStart[iColumn];
1195           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1196        int iRow = row[j];
1197        if (sos[iRow]&&rowActivity[iRow]) {
1198          possible = false;
1199        } else {
1200          double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1201          if (gap<element[j])
1202            possible = false;
1203        }
1204      }
1205      if (possible) {
1206        // Increase chosen column
1207        newSolution[iColumn] = 1.0;
1208        double cost =  modifiedCost[iColumn];
1209        newSolutionValue += cost;
1210        for (CoinBigIndex j = columnStart[iColumn];
1211             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1212          int iRow = row[j];
1213          rowActivity[iRow] += element[j];
1214        }
1215      }
1216    }
1217    delete [] sos;
1218    if (newSolutionValue < solutionValue) {
1219        // check feasible
1220      const double * rowLower = solver->getRowLower();
1221      const double * rowUpper = solver->getRowUpper();
1222      memset(rowActivity, 0, numberRows*sizeof(double));
1223        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1224            CoinBigIndex j;
1225            double value = newSolution[iColumn];
1226            if (value) {
1227                for (j = columnStart[iColumn];
1228                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1229                    int iRow = row[j];
1230                    rowActivity[iRow] += value * element[j];
1231                }
1232            }
1233        }
1234        // check was approximately feasible
1235        bool feasible = true;
1236        for (iRow = 0; iRow < numberRows; iRow++) {
1237            if (rowActivity[iRow] < rowLower[iRow]) {
1238                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
1239                    feasible = false;
1240            } else if (rowActivity[iRow] > rowUpper[iRow]) {
1241                if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance)
1242                    feasible = false;
1243            }
1244        }
1245        if (feasible) {
1246            // new solution
1247            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1248            solutionValue = newSolutionValue;
1249            //printf("** Solution of %g found by rounding\n",newSolutionValue);
1250            returnCode = 1;
1251        } else {
1252            // Can easily happen
1253            //printf("Debug CbcHeuristicGreedySOS giving bad solution\n");
1254        }
1255    }
1256    delete [] sosRow;
1257    delete [] newSolution;
1258    delete [] rowActivity;
1259    delete [] modifiedCost;
1260    delete [] contribution;
1261    delete [] which;
1262    delete [] rhs;
1263    return returnCode;
1264}
1265// update model
1266void CbcHeuristicGreedySOS::setModel(CbcModel * model)
1267{
1268    delete [] originalRhs_;
1269    gutsOfConstructor(model);
1270    validate();
1271}
1272// Resets stuff if model changes
1273void
1274CbcHeuristicGreedySOS::resetModel(CbcModel * model)
1275{
1276    delete [] originalRhs_;
1277    gutsOfConstructor(model);
1278}
1279// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1280void
1281CbcHeuristicGreedySOS::validate()
1282{
1283    if (model_ && when() < 10) {
1284        if (model_->numberIntegers() !=
1285                model_->numberObjects() && (model_->numberObjects() ||
1286                                            (model_->specialOptions()&1024) == 0)) {
1287            int numberOdd = 0;
1288            for (int i = 0; i < model_->numberObjects(); i++) {
1289                if (!model_->object(i)->canDoHeuristics())
1290                    numberOdd++;
1291            }
1292            if (numberOdd)
1293                setWhen(0);
1294        }
1295        // Only works if coefficients positive and all rows L or SOS
1296        OsiSolverInterface * solver = model_->solver();
1297        const double * columnUpper = solver->getColUpper();
1298        const double * columnLower = solver->getColLower();
1299        const double * rowLower = solver->getRowLower();
1300        const double * rowUpper = solver->getRowUpper();
1301
1302        int numberRows = solver->getNumRows();
1303        // Column copy
1304        const double * element = matrix_.getElements();
1305        const int * row = matrix_.getIndices();
1306        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1307        const int * columnLength = matrix_.getVectorLengths();
1308        bool good = true;
1309        assert (originalRhs_);
1310        for (int iRow = 0; iRow < numberRows; iRow++) {
1311          if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1312            // SOS
1313            originalRhs_[iRow]=-1.0;
1314          } else if (rowLower[iRow] > 0.0) {
1315                good = false;
1316          } else if (rowUpper[iRow] < 0.0) {
1317                good = false;
1318          } else {
1319            originalRhs_[iRow]=rowUpper[iRow];
1320          }
1321        }
1322        int numberColumns = solver->getNumCols();
1323        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1324            if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1325                good = false;
1326            CoinBigIndex j;
1327            int nSOS=0;
1328            for (j = columnStart[iColumn];
1329                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1330                if (element[j] < 0.0)
1331                    good = false;
1332                int iRow = row[j];
1333                if (originalRhs_[iRow]==-1.0) {
1334                  if (element[j] != 1.0)
1335                    good = false;
1336                  nSOS++;
1337                }
1338            }
1339            if (nSOS>1||!solver->isBinary(iColumn))
1340              good = false;
1341        }
1342        if (!good)
1343            setWhen(0); // switch off
1344    }
1345}
Note: See TracBrowser for help on using the repository browser.