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

Last change on this file since 2094 was 2094, checked in by forrest, 3 years ago

for memory leaks and heuristics and some experimental stuff

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