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

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

add some heuristic variants

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