source: stable/2.7/Cbc/src/CbcHeuristicGreedy.cpp @ 1797

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

fix for overlapping sos
  • 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 1797 2012-10-11 09:43:41Z 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          // Only in one sos
1122          break;
1123        }
1124      }
1125      modifiedCost[iColumn] = cost;
1126    }
1127    delete [] slackCost;
1128    double offset;
1129    solver->getDblParam(OsiObjOffset, offset);
1130    double newSolutionValue = -offset+offset2;
1131    int returnCode = 0;
1132
1133
1134    // Get solution array for heuristic solution
1135    double * newSolution = new double [numberColumns];
1136    double * rowActivity = new double[numberRows];
1137    double * contribution = new double [numberColumns];
1138    int * which = new int [numberColumns];
1139    double * newSolution0 = new double [numberColumns];
1140    if ((algorithm_&(2|4))==0) {
1141      // get solution as small as possible
1142      for (iColumn = 0; iColumn < numberColumns; iColumn++) 
1143        newSolution0[iColumn] = columnLower[iColumn];
1144    } else {
1145      // Get rounded down solution
1146      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1147        double value = solution[iColumn];
1148        // Round down integer
1149        if (fabs(floor(value + 0.5) - value) < integerTolerance) {
1150          value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
1151        } else {
1152          value = CoinMax(floor(value), columnLower[iColumn]);
1153        }
1154        // make sure clean
1155        value = CoinMin(value, columnUpper[iColumn]);
1156        value = CoinMax(value, columnLower[iColumn]);
1157        newSolution0[iColumn] = value;
1158      }
1159    }
1160    double * rowWeight = new double [numberRows];
1161    for (int i=0;i<numberRows;i++)
1162      rowWeight[i]=1.0;
1163    double costBias = 0.0;
1164    int nPass = ((algorithm_&4)!=0) ? 1 : 10;
1165    for (int iPass=0;iPass<nPass;iPass++) {
1166      newSolutionValue = -offset+offset2;
1167      memcpy(newSolution,newSolution0,numberColumns*sizeof(double));
1168      // get row activity
1169      memset(rowActivity, 0, numberRows*sizeof(double));
1170      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1171        CoinBigIndex j;
1172        double value = newSolution[iColumn];
1173        for (j = columnStart[iColumn];
1174             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1175          int iRow = row[j];
1176          rowActivity[iRow] += value * element[j];
1177        }
1178      }
1179      if (!rowWeight) {
1180        rowWeight = CoinCopyOfArray(rowActivity,numberRows);
1181      }
1182      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1183        CoinBigIndex j;
1184        double value = newSolution[iColumn];
1185        double cost =  modifiedCost[iColumn];
1186        double forSort = 1.0e-24;
1187        bool hasSlack=false;
1188        bool willFit=true;
1189        bool gRow=false;
1190        newSolutionValue += value * cost;
1191        cost += 1.0e-12;
1192        for (j = columnStart[iColumn];
1193             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1194          int iRow = row[j];
1195          int type = sos[iRow];
1196          double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1197          switch (type) {
1198          case -1:
1199            // G row
1200            gRow = true;
1201#if 0
1202            if (rhs[iRow]>rowWeight[iRow]||(algorithm_&(2|4))!=0)
1203              forSort += element[j];
1204            else
1205              forSort += 0.1*element[j];
1206#else
1207            forSort += rowWeight[iRow]*element[j];
1208#endif
1209            break;
1210          case 0:
1211            // L row
1212            if (gap<element[j]) {
1213              willFit = false;
1214            } else {
1215              forSort += element[j];
1216            }
1217            break;
1218          case 1:
1219            // SOS without slack
1220            if (gap<element[j]) {
1221              willFit = false;
1222            }
1223            break;
1224          case 2:
1225            // SOS with slack
1226            hasSlack = true;
1227            if (gap<element[j]) {
1228              willFit = false;
1229            }
1230            break;
1231          }
1232        }
1233        bool isSlack = hasSlack && (columnLength[iColumn]==1);
1234        if (forSort<1.0e-24)
1235          forSort = 1.0e-12;
1236        if ((algorithm_&4)!=0 && forSort > 1.0e-24) 
1237          forSort=1.0;
1238        // Use smallest cost if will fit
1239        if (willFit && (hasSlack||gRow) && 
1240            value == 0.0 && columnUpper[iColumn]) {
1241          if (hasSlack && !gRow) {
1242            if (cost>1.0e-12) {
1243              forSort = 2.0e30;
1244            } else if (cost==1.0e-12) {
1245              if (!isSlack)
1246                forSort = 1.0e29;
1247              else
1248                forSort = 1.0e28;
1249            } else {
1250              forSort = cost/forSort;
1251            }
1252          } else {
1253            if (!gRow||true)
1254              forSort = (cost+costBias)/forSort;
1255            else
1256              forSort = 1.0e-12/forSort;
1257          }
1258        } else {
1259          // put at end
1260          forSort = 1.0e30;
1261        }
1262        which[iColumn]=iColumn;
1263        contribution[iColumn]= forSort;
1264      }
1265      CoinSort_2(contribution,contribution+numberColumns,which);
1266      // Go through columns
1267      int nAdded=0;
1268      int nSlacks=0;
1269      for (int jColumn = 0; jColumn < numberColumns; jColumn++) {
1270        if (contribution[jColumn]>=1.0e30)
1271          break;
1272        int iColumn = which[jColumn];
1273        double value = newSolution[iColumn];
1274        if (value)
1275          continue;
1276        bool possible = true;
1277        CoinBigIndex j;
1278        for (j = columnStart[iColumn];
1279             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1280          int iRow = row[j];
1281          if (sos[iRow]>0&&rowActivity[iRow]) {
1282            possible = false;
1283          } else {
1284            double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1285            if (gap<element[j]&&sos[iRow]>=0)
1286              possible = false;
1287          }
1288        }
1289        if (possible) {
1290          //#define REPORT 1
1291#ifdef REPORT
1292          if ((nAdded%1000)==0) {
1293            double gap=0.0;
1294            for (int i=0;i<numberRows;i++) {
1295              if (rowUpper[i]>1.0e20)
1296                gap += CoinMax(rowLower[i]-rowActivity[i],0.0);
1297            }
1298            if (gap)
1299              printf("after %d added gap %g - %d slacks\n",
1300                     nAdded,gap,nSlacks);
1301          }
1302#endif
1303          nAdded++;
1304          if (columnLength[iColumn]==1)
1305            nSlacks++;
1306          // Increase chosen column
1307          newSolution[iColumn] = 1.0;
1308          double cost =  modifiedCost[iColumn];
1309          newSolutionValue += cost;
1310          for (CoinBigIndex j = columnStart[iColumn];
1311               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1312            int iRow = row[j];
1313            rowActivity[iRow] += element[j];
1314          }
1315        }
1316      }
1317#ifdef REPORT
1318      {
1319        double under=0.0;
1320        double over=0.0;
1321        double gap = 0.0;
1322        int nUnder=0;
1323        int nOver=0;
1324        int nGap=0;
1325        for (iRow = 0; iRow < numberRows; iRow++) {
1326          if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1327            double value = rowLower[iRow]-rowActivity[iRow];
1328#if REPORT>1
1329            printf("below on %d is %g - activity %g lower %g\n",
1330                   iRow,value,rowActivity[iRow],rowLower[iRow]);
1331#endif
1332            under += value;
1333            nUnder++;
1334          } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) {
1335            double value = rowActivity[iRow]-rowUpper[iRow];
1336#if REPORT>1
1337            printf("above on %d is %g - activity %g upper %g\n",
1338                   iRow,value,rowActivity[iRow],rowUpper[iRow]);
1339#endif
1340            over += value;
1341            nOver++;
1342          } else {
1343            double value = rowActivity[iRow]-rowLower[iRow];
1344            if (value && value < 1.0e20) {
1345#if REPORT>1
1346              printf("gap on %d is %g - activity %g lower %g\n",
1347                     iRow,value,rowActivity[iRow],rowLower[iRow]);
1348#endif
1349              gap += value;
1350              nGap++;
1351            }
1352          }
1353        }
1354        printf("final under %g (%d) - over %g (%d) - free %g (%d) - %d added - solvalue %g\n",
1355               under,nUnder,over,nOver,gap,nGap,nAdded,newSolutionValue);
1356      }
1357#endif
1358      double gap = 0.0;
1359      double over = 0.0;
1360      int nL=0;
1361      int nG=0;
1362      int iUnder=-1;
1363      int nUnder=0;
1364      for (iRow = 0; iRow < numberRows; iRow++) {
1365        if (rowLower[iRow]<-1.0e20)
1366          nL++;
1367        if (rowUpper[iRow]>1.0e20)
1368          nG++;
1369        if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1370          gap += rowLower[iRow]-rowActivity[iRow];
1371          nUnder++;
1372          iUnder=iRow;
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.