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

Last change on this file since 1586 was 1586, checked in by forrest, 10 years ago

fix minor bug
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.9 KB
Line 
1/* $Id: CbcHeuristicGreedy.cpp 1586 2011-01-21 17:57:01Z 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    // If bit set then use current
992    if ((algorithm_&1)!=0) {
993      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
994      element = matrix->getElements();
995      row = matrix->getIndices();
996      columnStart = matrix->getVectorStarts();
997      columnLength = matrix->getVectorLengths();
998      rhs = new double [numberRows];
999      const double * rowLower = solver->getRowLower();
1000      const double * rowUpper = solver->getRowUpper();
1001      bool good = true;
1002      for (int iRow = 0; iRow < numberRows; iRow++) {
1003        if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1004          // SOS
1005          rhs[iRow]=-1.0;
1006        } else if (rowLower[iRow] > 0.0) {
1007          good = false;
1008        } else if (rowUpper[iRow] < 0.0) {
1009          good = false;
1010        } else {
1011          rhs[iRow]=rowUpper[iRow];
1012        }
1013      }
1014      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1015        if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1016          good = false;
1017        CoinBigIndex j;
1018        int nSOS=0;
1019        if (!solver->isInteger(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) {
1027            if (element[j] != 1.0)
1028              good = false;
1029            nSOS++;
1030          }
1031        }
1032        if (nSOS!=1)
1033          good = false;
1034      }
1035      if (!good) {
1036        delete [] rhs;
1037        setWhen(0); // switch off
1038        return 0;
1039      }
1040    }
1041    const double * solution = solver->getColSolution();
1042    const double * objective = solver->getObjCoefficients();
1043    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1044    double primalTolerance;
1045    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1046
1047    numRuns_++;
1048    assert (numberRows == matrix_.getNumRows());
1049    int iRow, iColumn;
1050    double direction = solver->getObjSense();
1051    double * slackCost = new double [numberRows];
1052    double * modifiedCost = CoinCopyOfArray(objective,numberColumns);
1053    for (int iRow = 0;iRow < numberRows; iRow++)
1054      slackCost[iRow]=1.0e30;
1055    // Take off cost of gub slack
1056    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1057      if (columnLength[iColumn] == 1) {
1058        // SOS slack
1059        double cost = direction*objective[iColumn];
1060        int iRow = row[columnStart[iColumn]];
1061        assert (rhs[iRow]<0.0);
1062        slackCost[iRow]=CoinMin(slackCost[iRow],cost);
1063      }
1064    }
1065    double offset2 = 0.0;
1066    char * sos = new char [numberRows];
1067    for (int iRow = 0;iRow < numberRows; iRow++) {
1068      sos[iRow]=0;
1069      if (rhs[iRow]<0.0) {
1070        sos[iRow]=1;
1071        rhs[iRow]=1.0;
1072      }
1073      if( slackCost[iRow] == 1.0e30) {
1074        slackCost[iRow]=0.0;
1075      } else {
1076        offset2 += slackCost[iRow];
1077        sos[iRow] = 2;
1078      }
1079    }
1080    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1081      double cost = direction * modifiedCost[iColumn];
1082      CoinBigIndex j;
1083      for (j = columnStart[iColumn];
1084           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1085        int iRow = row[j];
1086        if (sos[iRow]) {
1087          cost -= slackCost[iRow];
1088        }
1089      }
1090      modifiedCost[iColumn] = cost;
1091    }
1092    delete [] slackCost;
1093    double offset;
1094    solver->getDblParam(OsiObjOffset, offset);
1095    double newSolutionValue = -offset+offset2;
1096    int returnCode = 0;
1097
1098
1099    // Get solution array for heuristic solution
1100    double * newSolution = new double [numberColumns];
1101    double * rowActivity = new double[numberRows];
1102    memset(rowActivity, 0, numberRows*sizeof(double));
1103    if ((algorithm_&(2|4))==0) {
1104      // get solution as small as possible
1105      for (iColumn = 0; iColumn < numberColumns; iColumn++) 
1106        newSolution[iColumn] = columnLower[iColumn];
1107    } else {
1108      // Get rounded down solution
1109      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1110        double value = solution[iColumn];
1111        // Round down integer
1112        if (fabs(floor(value + 0.5) - value) < integerTolerance) {
1113          value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
1114        } else {
1115          value = CoinMax(floor(value), columnLower[iColumn]);
1116        }
1117        // make sure clean
1118        value = CoinMin(value, columnUpper[iColumn]);
1119        value = CoinMax(value, columnLower[iColumn]);
1120        newSolution[iColumn] = value;
1121      }
1122    }
1123    // get row activity
1124    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1125      CoinBigIndex j;
1126      double value = newSolution[iColumn];
1127      for (j = columnStart[iColumn];
1128           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1129        int iRow = row[j];
1130        rowActivity[iRow] += value * element[j];
1131      }
1132    }
1133    double * contribution = new double [numberColumns];
1134    int * which = new int [numberColumns];
1135    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1136      CoinBigIndex j;
1137      double value = newSolution[iColumn];
1138      double cost =  modifiedCost[iColumn];
1139      double forSort = 0.0;
1140      bool hasSlack=false;
1141      bool willFit=true;
1142      newSolutionValue += value * cost;
1143      for (j = columnStart[iColumn];
1144           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1145        int iRow = row[j];
1146        if (sos[iRow] == 2) 
1147          hasSlack = true;
1148        forSort += element[j];
1149        double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1150        if (gap<element[j])
1151          willFit = false;
1152      }
1153      bool isSlack = hasSlack && (columnLength[iColumn]==1);
1154      if ((algorithm_&4)!=0) 
1155        forSort=1.0;
1156      // Use smallest cost if will fit
1157      if (willFit && hasSlack && 
1158          value == 0.0 && columnUpper[iColumn]) {
1159        if (hasSlack) {
1160          if (cost>0.0) {
1161            forSort = 2.0e30;
1162          } else if (cost==0.0) {
1163            if (!isSlack)
1164              forSort = 1.0e29;
1165            else
1166              forSort = 1.0e28;
1167          } else {
1168            forSort = cost/forSort;
1169          }
1170        } else {
1171          forSort = cost/forSort;
1172        }
1173      } else {
1174        // put at end
1175        forSort = 1.0e30;
1176      }
1177      which[iColumn]=iColumn;
1178      contribution[iColumn]= forSort;
1179    }
1180    CoinSort_2(contribution,contribution+numberColumns,which);
1181    // Go through columns
1182    for (int jColumn = 0; jColumn < numberColumns; jColumn++) {
1183      int iColumn = which[jColumn];
1184      double value = newSolution[iColumn];
1185      if (value)
1186        continue;
1187      bool possible = true;
1188      CoinBigIndex j;
1189      for (j = columnStart[iColumn];
1190           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1191        int iRow = row[j];
1192        if (sos[iRow]&&rowActivity[iRow]) {
1193          possible = false;
1194        } else {
1195          double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1196          if (gap<element[j])
1197            possible = false;
1198        }
1199      }
1200      if (possible) {
1201        // Increase chosen column
1202        newSolution[iColumn] = 1.0;
1203        double cost =  modifiedCost[iColumn];
1204        newSolutionValue += cost;
1205        for (CoinBigIndex j = columnStart[iColumn];
1206             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1207          int iRow = row[j];
1208          rowActivity[iRow] += element[j];
1209        }
1210      }
1211    }
1212    delete [] sos;
1213    if (newSolutionValue < solutionValue) {
1214        // check feasible
1215      const double * rowLower = solver->getRowLower();
1216      const double * rowUpper = solver->getRowUpper();
1217      memset(rowActivity, 0, numberRows*sizeof(double));
1218        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1219            CoinBigIndex j;
1220            double value = newSolution[iColumn];
1221            if (value) {
1222                for (j = columnStart[iColumn];
1223                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1224                    int iRow = row[j];
1225                    rowActivity[iRow] += value * element[j];
1226                }
1227            }
1228        }
1229        // check was approximately feasible
1230        bool feasible = true;
1231        for (iRow = 0; iRow < numberRows; iRow++) {
1232            if (rowActivity[iRow] < rowLower[iRow]) {
1233                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
1234                    feasible = false;
1235            } else if (rowActivity[iRow] > rowUpper[iRow]) {
1236                if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance)
1237                    feasible = false;
1238            }
1239        }
1240        if (feasible) {
1241            // new solution
1242            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1243            solutionValue = newSolutionValue;
1244            //printf("** Solution of %g found by rounding\n",newSolutionValue);
1245            returnCode = 1;
1246        } else {
1247            // Can easily happen
1248            //printf("Debug CbcHeuristicGreedySOS giving bad solution\n");
1249        }
1250    }
1251    delete [] newSolution;
1252    delete [] rowActivity;
1253    delete [] modifiedCost;
1254    delete [] contribution;
1255    delete [] which;
1256    delete [] rhs;
1257    return returnCode;
1258}
1259// update model
1260void CbcHeuristicGreedySOS::setModel(CbcModel * model)
1261{
1262    delete [] originalRhs_;
1263    gutsOfConstructor(model);
1264    validate();
1265}
1266// Resets stuff if model changes
1267void
1268CbcHeuristicGreedySOS::resetModel(CbcModel * model)
1269{
1270    delete [] originalRhs_;
1271    gutsOfConstructor(model);
1272}
1273// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1274void
1275CbcHeuristicGreedySOS::validate()
1276{
1277    if (model_ && when() < 10) {
1278        if (model_->numberIntegers() !=
1279                model_->numberObjects() && (model_->numberObjects() ||
1280                                            (model_->specialOptions()&1024) == 0)) {
1281            int numberOdd = 0;
1282            for (int i = 0; i < model_->numberObjects(); i++) {
1283                if (!model_->object(i)->canDoHeuristics())
1284                    numberOdd++;
1285            }
1286            if (numberOdd)
1287                setWhen(0);
1288        }
1289        // Only works if coefficients positive and all rows L or SOS
1290        OsiSolverInterface * solver = model_->solver();
1291        const double * columnUpper = solver->getColUpper();
1292        const double * columnLower = solver->getColLower();
1293        const double * rowLower = solver->getRowLower();
1294        const double * rowUpper = solver->getRowUpper();
1295
1296        int numberRows = solver->getNumRows();
1297        // Column copy
1298        const double * element = matrix_.getElements();
1299        const int * row = matrix_.getIndices();
1300        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1301        const int * columnLength = matrix_.getVectorLengths();
1302        bool good = true;
1303        assert (originalRhs_);
1304        for (int iRow = 0; iRow < numberRows; iRow++) {
1305          if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1306            // SOS
1307            originalRhs_[iRow]=-1.0;
1308          } else if (rowLower[iRow] > 0.0) {
1309                good = false;
1310          } else if (rowUpper[iRow] < 0.0) {
1311                good = false;
1312          } else {
1313            originalRhs_[iRow]=rowUpper[iRow];
1314          }
1315        }
1316        int numberColumns = solver->getNumCols();
1317        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1318            if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1319                good = false;
1320            CoinBigIndex j;
1321            int nSOS=0;
1322            if (!solver->isInteger(iColumn))
1323              good = false;
1324            for (j = columnStart[iColumn];
1325                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1326                if (element[j] < 0.0)
1327                    good = false;
1328                int iRow = row[j];
1329                if (originalRhs_[iRow]==-1.0) {
1330                  if (element[j] != 1.0)
1331                    good = false;
1332                  nSOS++;
1333                }
1334            }
1335            if (nSOS!=1)
1336              good = false;
1337        }
1338        if (!good)
1339            setWhen(0); // switch off
1340    }
1341}
Note: See TracBrowser for help on using the repository browser.