source: stable/2.8/Cbc/src/CbcHeuristicGreedy.cpp @ 2128

Last change on this file since 2128 was 1888, checked in by stefan, 6 years ago

sync with trunk rev 1887

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