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

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

modifications to heuristics and allow missing out some printout

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