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

Last change on this file since 2216 was 2216, checked in by forrest, 5 years ago

fix bug in greedy heuristic
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 61.3 KB
Line 
1/* $Id: CbcHeuristicGreedy.cpp 2216 2015-10-06 16:09:32Z 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    char * sos = new char [numberRows];
999    memset(sos,'a',numberRows);
1000    int nonSOS=0;
1001    // If bit set then use current
1002    if ((algorithm_&1)!=0) {
1003      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
1004      element = matrix->getElements();
1005      row = matrix->getIndices();
1006      columnStart = matrix->getVectorStarts();
1007      columnLength = matrix->getVectorLengths();
1008      //rhs = new double [numberRows];
1009      const double * rowLower = solver->getRowLower();
1010      const double * rowUpper = solver->getRowUpper();
1011      bool good = true;
1012      for (int iRow = 0; iRow < numberRows; iRow++) {
1013        if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1014          // SOS
1015          rhs[iRow]=-1.0;
1016        } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) {
1017          good = false;
1018        } else if (rowUpper[iRow] < 0.0) {
1019          good = false;
1020        } else if (rowUpper[iRow] < 1.0e10) {
1021          rhs[iRow]=rowUpper[iRow];
1022          if (rhs[iRow]<0)
1023            sos[iRow]=0; // can't be SOS
1024        } else {
1025          rhs[iRow]=rowLower[iRow];
1026          if (rhs[iRow]<0)
1027            sos[iRow]=0; // can't be SOS
1028        }
1029      }
1030      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1031        if (!columnLength[iColumn])
1032          continue;
1033        if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1034          good = false;
1035        CoinBigIndex j;
1036        int nSOS=0;
1037        int iSOS=-1;
1038        if (!solver->isInteger(iColumn))
1039          good = false;
1040        for (j = columnStart[iColumn];
1041             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1042          if (element[j] < 0.0)
1043            good = false;
1044          int iRow = row[j];
1045          if (rhs[iRow]==-1.0 && sos[iRow] == 'a') {
1046            if (element[j] != 1.0)
1047              good = false;
1048            iSOS=iRow;
1049            nSOS++;
1050          }
1051        }
1052        if (nSOS>1)
1053          good = false;
1054        else if (!nSOS)
1055          nonSOS++;
1056        sosRow[iColumn] = iSOS;
1057      }
1058      if (!good) {
1059        delete [] sosRow;
1060        delete [] rhs;
1061        delete [] sos;
1062        setWhen(0); // switch off
1063        return 0;
1064      }
1065    } else {
1066      abort(); // not allowed yet
1067    }
1068    const double * solution = solver->getColSolution();
1069    const double * objective = solver->getObjCoefficients();
1070    const double * rowLower = solver->getRowLower();
1071    const double * rowUpper = solver->getRowUpper();
1072    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1073    double primalTolerance;
1074    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1075
1076    numRuns_++;
1077    assert (numberRows == matrix_.getNumRows());
1078    // set up linked list for sets
1079    int * firstGub = new int [numberRows];
1080    int * nextGub = new int [numberColumns];
1081    int iRow, iColumn;
1082    double direction = solver->getObjSense();
1083    double * slackCost = new double [numberRows];
1084    double * modifiedCost = CoinCopyOfArray(objective,numberColumns);
1085    for (int iRow = 0;iRow < numberRows; iRow++) {
1086      slackCost[iRow]=1.0e30;
1087      firstGub[iRow]=-1;
1088    }
1089    // Take off cost of gub slack
1090    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1091      nextGub[iColumn]=-1;
1092      int iRow = sosRow[iColumn];
1093      if (columnLength[iColumn] == 1&&iRow>=0) {
1094        // SOS slack
1095        double cost = direction*objective[iColumn];
1096        assert (rhs[iRow]<0.0);
1097        slackCost[iRow]=CoinMin(slackCost[iRow],cost);
1098      }
1099    }
1100    double offset2 = 0.0;
1101    for (int iRow = 0;iRow < numberRows; iRow++) {
1102      if (sos[iRow]=='a') {
1103        // row is possible
1104        sos[iRow]=0;
1105        if (rhs[iRow]<0.0) {
1106          sos[iRow]=1;
1107          rhs[iRow]=1.0;
1108        } else if (rhs[iRow] != rowUpper[iRow]) {
1109          // G row
1110          sos[iRow]=-1;
1111        }
1112        if( slackCost[iRow] == 1.0e30) {
1113          slackCost[iRow]=0.0;
1114        } else {
1115          offset2 += slackCost[iRow];
1116          sos[iRow] = 2;
1117        }
1118      }
1119    }
1120    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1121      double cost = direction * modifiedCost[iColumn];
1122      CoinBigIndex j;
1123      for (j = columnStart[iColumn];
1124           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1125        int iRow = row[j];
1126        if (sos[iRow]>0) {
1127          cost -= slackCost[iRow];
1128          if (firstGub[iRow]<0) {
1129            firstGub[iRow]=iColumn;
1130          } else {
1131            int jColumn = firstGub[iRow];
1132            while (nextGub[jColumn]>=0)
1133              jColumn=nextGub[jColumn];
1134            nextGub[jColumn]=iColumn;
1135          }
1136          // Only in one sos
1137          break;
1138        }
1139      }
1140      modifiedCost[iColumn] = cost;
1141    }
1142    delete [] slackCost;
1143    double offset;
1144    solver->getDblParam(OsiObjOffset, offset);
1145    double newSolutionValue = -offset+offset2;
1146    int returnCode = 0;
1147
1148
1149    // Get solution array for heuristic solution
1150    double * newSolution = new double [numberColumns];
1151    double * rowActivity = new double[numberRows];
1152    double * contribution = new double [numberColumns];
1153    int * which = new int [numberColumns];
1154    double * newSolution0 = new double [numberColumns];
1155    if ((algorithm_&(2|4))==0) {
1156      // get solution as small as possible
1157      for (iColumn = 0; iColumn < numberColumns; iColumn++) 
1158        newSolution0[iColumn] = columnLower[iColumn];
1159    } else {
1160      // Get rounded down solution
1161      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1162        double value = solution[iColumn];
1163        // Round down integer
1164        if (fabs(floor(value + 0.5) - value) < integerTolerance) {
1165          value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
1166        } else {
1167          value = CoinMax(floor(value), columnLower[iColumn]);
1168        }
1169        // make sure clean
1170        value = CoinMin(value, columnUpper[iColumn]);
1171        value = CoinMax(value, columnLower[iColumn]);
1172        newSolution0[iColumn] = value;
1173      }
1174    }
1175    double * rowWeight = new double [numberRows];
1176    for (int i=0;i<numberRows;i++)
1177      rowWeight[i]=1.0;
1178    double costBias = 0.0;
1179    int nPass = ((algorithm_&4)!=0) ? 1 : 10;
1180    for (int iPass=0;iPass<nPass;iPass++) {
1181      newSolutionValue = -offset+offset2;
1182      memcpy(newSolution,newSolution0,numberColumns*sizeof(double));
1183      // get row activity
1184      memset(rowActivity, 0, numberRows*sizeof(double));
1185      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1186        CoinBigIndex j;
1187        double value = newSolution[iColumn];
1188        for (j = columnStart[iColumn];
1189             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1190          int iRow = row[j];
1191          rowActivity[iRow] += value * element[j];
1192        }
1193      }
1194      if (!rowWeight) {
1195        rowWeight = CoinCopyOfArray(rowActivity,numberRows);
1196      }
1197      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1198        CoinBigIndex j;
1199        double value = newSolution[iColumn];
1200        double cost =  modifiedCost[iColumn];
1201        double forSort = 1.0e-24;
1202        bool hasSlack=false;
1203        bool willFit=true;
1204        bool gRow=false;
1205        newSolutionValue += value * cost;
1206        cost += 1.0e-12;
1207        for (j = columnStart[iColumn];
1208             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1209          int iRow = row[j];
1210          int type = sos[iRow];
1211          double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1212          switch (type) {
1213          case -1:
1214            // G row
1215            gRow = true;
1216#if 0
1217            if (rhs[iRow]>rowWeight[iRow]||(algorithm_&(2|4))!=0)
1218              forSort += element[j];
1219            else
1220              forSort += 0.1*element[j];
1221#else
1222            forSort += rowWeight[iRow]*element[j];
1223#endif
1224            break;
1225          case 0:
1226            // L row
1227            if (gap<element[j]) {
1228              willFit = false;
1229            } else {
1230              forSort += element[j];
1231            }
1232            break;
1233          case 1:
1234            // SOS without slack
1235            if (gap<element[j]) {
1236              willFit = false;
1237            }
1238            break;
1239          case 2:
1240            // SOS with slack
1241            hasSlack = true;
1242            if (gap<element[j]) {
1243              willFit = false;
1244            }
1245            break;
1246          }
1247        }
1248        bool isSlack = hasSlack && (columnLength[iColumn]==1);
1249        if (forSort<1.0e-24)
1250          forSort = 1.0e-12;
1251        if ((algorithm_&4)!=0 && forSort > 1.0e-24) 
1252          forSort=1.0;
1253        // Use smallest cost if will fit
1254        if (willFit && (hasSlack||gRow) && 
1255            value == 0.0 && columnUpper[iColumn]) {
1256          if (hasSlack && !gRow) {
1257            if (cost>1.0e-12) {
1258              forSort = 2.0e30;
1259            } else if (cost==1.0e-12) {
1260              if (!isSlack)
1261                forSort = 1.0e29;
1262              else
1263                forSort = 1.0e28;
1264            } else {
1265              forSort = cost/forSort;
1266            }
1267          } else {
1268            if (!gRow||true)
1269              forSort = (cost+costBias)/forSort;
1270            else
1271              forSort = 1.0e-12/forSort;
1272          }
1273        } else {
1274          // put at end
1275          forSort = 1.0e30;
1276        }
1277        which[iColumn]=iColumn;
1278        contribution[iColumn]= forSort;
1279      }
1280      CoinSort_2(contribution,contribution+numberColumns,which);
1281      // Go through columns
1282      int nAdded=0;
1283      int nSlacks=0;
1284      for (int jColumn = 0; jColumn < numberColumns; jColumn++) {
1285        if (contribution[jColumn]>=1.0e30)
1286          break;
1287        int iColumn = which[jColumn];
1288        double value = newSolution[iColumn];
1289        if (value)
1290          continue;
1291        bool possible = true;
1292        CoinBigIndex j;
1293        for (j = columnStart[iColumn];
1294             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1295          int iRow = row[j];
1296          if (sos[iRow]>0&&rowActivity[iRow]) {
1297            possible = false;
1298          } else {
1299            double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
1300            if (gap<element[j]&&sos[iRow]>=0)
1301              possible = false;
1302          }
1303        }
1304        if (possible) {
1305          //#define REPORT 1
1306#ifdef REPORT
1307          if ((nAdded%1000)==0) {
1308            double gap=0.0;
1309            for (int i=0;i<numberRows;i++) {
1310              if (rowUpper[i]>1.0e20)
1311                gap += CoinMax(rowLower[i]-rowActivity[i],0.0);
1312            }
1313            if (gap)
1314              printf("after %d added gap %g - %d slacks\n",
1315                     nAdded,gap,nSlacks);
1316          }
1317#endif
1318          nAdded++;
1319          if (columnLength[iColumn]==1)
1320            nSlacks++;
1321          // Increase chosen column
1322          newSolution[iColumn] = 1.0;
1323          double cost =  modifiedCost[iColumn];
1324          newSolutionValue += cost;
1325          for (CoinBigIndex j = columnStart[iColumn];
1326               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1327            int iRow = row[j];
1328            rowActivity[iRow] += element[j];
1329          }
1330        }
1331      }
1332#ifdef REPORT
1333      {
1334        double under=0.0;
1335        double over=0.0;
1336        double gap = 0.0;
1337        int nUnder=0;
1338        int nOver=0;
1339        int nGap=0;
1340        for (iRow = 0; iRow < numberRows; iRow++) {
1341          if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1342            double value = rowLower[iRow]-rowActivity[iRow];
1343#if REPORT>1
1344            printf("below on %d is %g - activity %g lower %g\n",
1345                   iRow,value,rowActivity[iRow],rowLower[iRow]);
1346#endif
1347            under += value;
1348            nUnder++;
1349          } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) {
1350            double value = rowActivity[iRow]-rowUpper[iRow];
1351#if REPORT>1
1352            printf("above on %d is %g - activity %g upper %g\n",
1353                   iRow,value,rowActivity[iRow],rowUpper[iRow]);
1354#endif
1355            over += value;
1356            nOver++;
1357          } else {
1358            double value = rowActivity[iRow]-rowLower[iRow];
1359            if (value && value < 1.0e20) {
1360#if REPORT>1
1361              printf("gap on %d is %g - activity %g lower %g\n",
1362                     iRow,value,rowActivity[iRow],rowLower[iRow]);
1363#endif
1364              gap += value;
1365              nGap++;
1366            }
1367          }
1368        }
1369        printf("final under %g (%d) - over %g (%d) - free %g (%d) - %d added - solvalue %g\n",
1370               under,nUnder,over,nOver,gap,nGap,nAdded,newSolutionValue);
1371      }
1372#endif
1373      double gap = 0.0;
1374      double over = 0.0;
1375      int nL=0;
1376      int nG=0;
1377      int nUnder=0;
1378      for (iRow = 0; iRow < numberRows; iRow++) {
1379        if (rowLower[iRow]<-1.0e20)
1380          nL++;
1381        if (rowUpper[iRow]>1.0e20)
1382          nG++;
1383        if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1384          gap += rowLower[iRow]-rowActivity[iRow];
1385          nUnder++;
1386          rowWeight[iRow] *= 1.1;
1387        } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) {
1388          gap += rowActivity[iRow]-rowUpper[iRow];
1389        } else {
1390          over += rowActivity[iRow]-rowLower[iRow];
1391          //rowWeight[iRow] *= 0.9;
1392        }
1393      }
1394      if (nG&&!nL) {
1395        // can we fix
1396        // get list of columns which can go down without making
1397        // things much worse
1398        int nPossible=0;
1399        int nEasyDown=0;
1400        int nSlackDown=0;
1401        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1402          if (newSolution[iColumn]&&
1403              columnUpper[iColumn]>columnLower[iColumn]) {
1404            bool canGoDown=true;
1405            bool under = false;
1406            int iSos=-1;
1407            for (CoinBigIndex j = columnStart[iColumn];
1408                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1409              int iRow = row[j];
1410              if (sos[iRow]<0) {
1411                double over = rowActivity[iRow]-rowLower[iRow];
1412                if (over>=0.0&&element[j]>over+1.0e-12) {
1413                  canGoDown=false;
1414                  break;
1415                } else if (over<0.0) {
1416                  under = true;
1417                }
1418              } else {
1419                iSos=iRow;
1420              }
1421            }
1422            if (canGoDown) { 
1423              if (!under) {
1424                if (iSos>=0) {
1425                  // find cheapest
1426                  double cheapest=modifiedCost[iColumn];
1427                  int iCheapest = -1;
1428                  int jColumn = firstGub[iSos];
1429                  assert (jColumn>=0);
1430                  while (jColumn>=0) {
1431                    if (modifiedCost[jColumn]<cheapest) {
1432                      cheapest=modifiedCost[jColumn];
1433                      iCheapest=jColumn;
1434                    }
1435                    jColumn = nextGub[jColumn];
1436                  }
1437                  if (iCheapest>=0) {
1438                    // Decrease column
1439                    newSolution[iColumn] = 0.0;
1440                    newSolutionValue -= modifiedCost[iColumn];
1441                    for (CoinBigIndex j = columnStart[iColumn];
1442                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1443                      int iRow = row[j];
1444                      rowActivity[iRow] -= element[j];
1445                    }
1446                    // Increase chosen column
1447                    newSolution[iCheapest] = 1.0;
1448                    newSolutionValue += modifiedCost[iCheapest]; 
1449                    for (CoinBigIndex j = columnStart[iCheapest];
1450                         j < columnStart[iCheapest] + columnLength[iCheapest]; j++) {
1451                      int iRow = row[j];
1452                      rowActivity[iRow] += element[j];
1453                    }
1454                    nEasyDown++;
1455                    if (columnLength[iColumn]>1) {
1456                      //printf("%d is easy down\n",iColumn);
1457                    } else {
1458                      nSlackDown++;
1459                    }
1460                  }
1461                } else if (modifiedCost[iColumn]>0.0) {
1462                  // easy down
1463                  // Decrease column
1464                  newSolution[iColumn] = 0.0;
1465                  newSolutionValue -= modifiedCost[iColumn];
1466                  for (CoinBigIndex j = columnStart[iColumn];
1467                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1468                    int iRow = row[j];
1469                    rowActivity[iRow] -= element[j];
1470                  }
1471                  nEasyDown++;
1472                }
1473              } else {
1474                which[nPossible++]=iColumn;
1475              }
1476            }
1477          }
1478        }
1479#ifdef REPORT
1480        printf("%d possible down, %d easy down of which %d are slacks\n",
1481               nPossible,nEasyDown,nSlackDown);
1482#endif
1483        double * needed = new double [numberRows];
1484        for (int i=0;i<numberRows;i++) {
1485          double value = rowLower[i] - rowActivity[i];
1486          if (value<1.0e-8)
1487            value=0.0;
1488          needed[i]=value;
1489        }
1490        if (gap && /*nUnder==1 &&*/ nonSOS) {
1491          double * weight = new double [numberColumns];
1492          int * sort = new int [numberColumns];
1493          // look at ones not in set
1494          int nPossible=0;
1495          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1496            if (!newSolution[iColumn]&&
1497                columnUpper[iColumn]>columnLower[iColumn]) {
1498              int iSos=-1;
1499              double value=0.0;
1500              for (CoinBigIndex j = columnStart[iColumn];
1501                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1502                int iRow = row[j];
1503                if (sos[iRow]<0) {
1504                  if (needed[iRow])
1505                    value += CoinMin(element[j]/needed[iRow],1.0);
1506                } else {
1507                  iSos=iRow;
1508                }
1509              }
1510              if (value && iSos<0) {
1511                weight[nPossible]=-value;
1512                sort[nPossible++]=iColumn;
1513              }
1514            }
1515          }
1516          CoinSort_2(weight,weight+nPossible,sort);
1517          for (int i=0;i<nPossible;i++) {
1518            int iColumn = sort[i];
1519            double helps=0.0;
1520            for (CoinBigIndex j = columnStart[iColumn];
1521                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1522              int iRow = row[j];
1523              if (needed[iRow]) 
1524                helps += CoinMin(needed[iRow],element[j]);
1525            }
1526            if (helps) {
1527              newSolution[iColumn] = 1.0;
1528              newSolutionValue += modifiedCost[iColumn];
1529              for (CoinBigIndex j = columnStart[iColumn];
1530                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1531                int iRow = row[j];
1532                rowActivity[iRow] += element[j];
1533                if (needed[iRow]) {
1534                  needed[iRow] -= element[j];
1535                  if (needed[iRow]<1.0e-8)
1536                    needed[iRow]=0.0;
1537                }
1538              }
1539              gap -= helps;
1540#ifdef REPORT
1541              {
1542                double gap2 = 0.0;
1543                for (iRow = 0; iRow < numberRows; iRow++) {
1544                  if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1545                    gap2 += rowLower[iRow]-rowActivity[iRow];
1546                  }
1547                }
1548                printf("estimated gap (nonsos) %g - computed %g\n",
1549                       gap,gap2);
1550              }
1551#endif
1552              if (gap<1.0e-12)
1553                break;
1554            }
1555          }
1556          delete [] weight;
1557          delete [] sort;
1558        }
1559        if (gap&&nPossible/*&&nUnder==1*/&&true&&model_->bestSolution()) {
1560          double * weight = new double [numberColumns];
1561          int * sort = new int [numberColumns];
1562          // look at ones in sets
1563          const double * goodSolution = model_->bestSolution();
1564          int nPossible=0;
1565          double largestWeight=0.0;
1566          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1567            if (!newSolution[iColumn]&&goodSolution[iColumn]&&
1568                columnUpper[iColumn]>columnLower[iColumn]) {
1569              int iSos=-1;
1570              double value=0.0;
1571              for (CoinBigIndex j = columnStart[iColumn];
1572                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1573                int iRow = row[j];
1574                if (sos[iRow]<0) {
1575                  if (needed[iRow]) 
1576                    value += CoinMin(element[j]/needed[iRow],1.0);
1577                } else {
1578                  iSos=iRow;
1579                }
1580              }
1581              if (value&&iSos>=0) {
1582                // see if value bigger than current
1583                int jColumn = firstGub[iSos];
1584                assert (jColumn>=0);
1585                while (jColumn>=0) {
1586                  if (newSolution[jColumn]) 
1587                    break;
1588                  jColumn = nextGub[jColumn];
1589                }
1590                assert (jColumn>=0);
1591                double value2=0.0;
1592                for (CoinBigIndex j = columnStart[jColumn];
1593                     j < columnStart[jColumn] + columnLength[jColumn]; j++) {
1594                  int iRow = row[j];
1595                  if (needed[iRow]) 
1596                    value2 += CoinMin(element[j]/needed[iRow],1.0);
1597                }
1598                if (value>value2) {
1599                  weight[nPossible]=-(value-value2);
1600                  largestWeight = CoinMax(largestWeight,(value-value2));
1601                  sort[nPossible++]=iColumn;
1602                }
1603              }
1604            }
1605          }
1606          if (nPossible) {
1607            double * temp = new double [numberRows];
1608            int * which2 = new int [numberRows];
1609            memset(temp,0,numberRows*sizeof(double));
1610            // modify so ones just more than gap best
1611            if (largestWeight>gap&&nUnder==1) {
1612              double offset = 4*largestWeight;
1613              for (int i=0;i<nPossible;i++) {
1614                double value = -weight[i];
1615                if (value>gap-1.0e-12)
1616                  weight[i] = -(offset-(value-gap));
1617              }
1618            }
1619            CoinSort_2(weight,weight+nPossible,sort);
1620            for (int i=0;i<nPossible;i++) {
1621              int iColumn = sort[i];
1622              int n=0;
1623              // find jColumn
1624              int iSos=-1;
1625              for (CoinBigIndex j = columnStart[iColumn];
1626                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1627                int iRow = row[j];
1628                temp[iRow]=element[j];
1629                which2[n++]=iRow;
1630                if (sos[iRow]>=0) {
1631                  iSos=iRow;
1632                }
1633              }
1634              int jColumn = firstGub[iSos];
1635              assert (jColumn>=0);
1636              while (jColumn>=0) {
1637                if (newSolution[jColumn]) 
1638                  break;
1639                jColumn = nextGub[jColumn];
1640              }
1641              assert (jColumn>=0);
1642              for (CoinBigIndex j = columnStart[jColumn];
1643                   j < columnStart[jColumn] + columnLength[jColumn]; j++) {
1644                int iRow = row[j];
1645                if (!temp[iRow])
1646                  which2[n++]=iRow;
1647                temp[iRow] -= element[j];
1648              }
1649              double helps = 0.0;
1650              for (int i=0;i<n;i++) {
1651                int iRow = which2[i];
1652                double newValue = rowActivity[iRow]+temp[iRow];
1653                if (temp[iRow]>1.0e-8) {
1654                  if (rowActivity[iRow]<rowLower[iRow]-1.0e-8) {
1655                    helps += CoinMin(temp[iRow],
1656                                     rowLower[iRow]-rowActivity[iRow]);
1657                  }
1658                } else if (temp[iRow]<-1.0e-8) {
1659                  if (newValue<rowLower[iRow]-1.0e-12) { 
1660                    helps -= CoinMin(-temp[iRow],
1661                                     1.0*(rowLower[iRow]-newValue));
1662                  }
1663                }
1664              }
1665              if (helps>0.0) {
1666                newSolution[iColumn]=1.0;
1667                newSolution[jColumn]=0.0;
1668                newSolutionValue += modifiedCost[iColumn]-modifiedCost[jColumn];
1669                for (int i=0;i<n;i++) {
1670                  int iRow = which2[i];
1671                  double newValue = rowActivity[iRow]+temp[iRow];
1672                  rowActivity[iRow] = newValue;
1673                  temp[iRow]=0.0;
1674                }
1675                gap -= helps;
1676#ifdef REPORT
1677                {
1678                  double gap2 = 0.0;
1679                  for (iRow = 0; iRow < numberRows; iRow++) {
1680                    if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1681                      gap2 += rowLower[iRow]-rowActivity[iRow];
1682                    }
1683                  }
1684                  printf("estimated gap %g - computed %g\n",
1685                         gap,gap2);
1686                }
1687#endif
1688                if (gap<1.0e-8)
1689                  break;
1690              } else {
1691                for (int i=0;i<n;i++) 
1692                  temp[which2[i]]=0.0;
1693              }
1694            }
1695            delete [] which2;
1696            delete [] temp;
1697          }
1698          delete [] weight;
1699          delete [] sort;
1700        }
1701        delete [] needed;
1702      }
1703#ifdef REPORT
1704      {
1705        double gap=0.0;
1706        double over = 0.0;
1707        for (iRow = 0; iRow < numberRows; iRow++) {
1708          if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) {
1709            double value = rowLower[iRow]-rowActivity[iRow];
1710#if REPORT>1
1711            printf("below on %d is %g - activity %g lower %g\n",
1712                   iRow,value,rowActivity[iRow],rowLower[iRow]);
1713#endif
1714            gap += value;
1715          } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) {
1716            double value = rowActivity[iRow]-rowUpper[iRow];
1717#if REPORT>1
1718            printf("above on %d is %g - activity %g upper %g\n",
1719                   iRow,value,rowActivity[iRow],rowUpper[iRow]);
1720#endif
1721            gap += value;
1722          } else {
1723            double value = rowActivity[iRow]-rowLower[iRow];
1724            if (value) {
1725#if REPORT>1
1726              printf("over on %d is %g - activity %g lower %g\n",
1727                     iRow,value,rowActivity[iRow],rowLower[iRow]);
1728#endif
1729              over += value;
1730            }
1731          }
1732        }
1733        printf("modified final gap %g - over %g - %d added - solvalue %g\n",
1734               gap,over,nAdded,newSolutionValue);
1735      }
1736#endif
1737      if (!gap) {
1738        break;
1739      } else {
1740        if (iPass==0) {
1741          costBias = 10.0*newSolutionValue/static_cast<double>(nAdded);
1742        } else {
1743          costBias *= 10.0;
1744        }
1745      }
1746    }
1747    delete [] newSolution0;
1748    delete [] rowWeight;
1749    delete [] sos;
1750    delete [] firstGub;
1751    delete [] nextGub;
1752    if (newSolutionValue < solutionValue) {
1753        // check feasible
1754      memset(rowActivity, 0, numberRows*sizeof(double));
1755        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1756            CoinBigIndex j;
1757            double value = newSolution[iColumn];
1758            if (value) {
1759                for (j = columnStart[iColumn];
1760                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1761                    int iRow = row[j];
1762                    rowActivity[iRow] += value * element[j];
1763                }
1764            }
1765        }
1766        // check was approximately feasible
1767        bool feasible = true;
1768        for (iRow = 0; iRow < numberRows; iRow++) {
1769            if (rowActivity[iRow] < rowLower[iRow]) {
1770                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
1771                    feasible = false;
1772            } else if (rowActivity[iRow] > rowUpper[iRow]) {
1773                if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance)
1774                    feasible = false;
1775            }
1776        }
1777        if (feasible) {
1778            // new solution
1779            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1780            solutionValue = newSolutionValue;
1781            //printf("** Solution of %g found by rounding\n",newSolutionValue);
1782            returnCode = 1;
1783        } else {
1784            // Can easily happen
1785            //printf("Debug CbcHeuristicGreedySOS giving bad solution\n");
1786        }
1787    }
1788    delete [] sosRow;
1789    delete [] newSolution;
1790    delete [] rowActivity;
1791    delete [] modifiedCost;
1792    delete [] contribution;
1793    delete [] which;
1794    delete [] rhs;
1795    return returnCode;
1796}
1797// update model
1798void CbcHeuristicGreedySOS::setModel(CbcModel * model)
1799{
1800    delete [] originalRhs_;
1801    gutsOfConstructor(model);
1802    validate();
1803}
1804// Resets stuff if model changes
1805void
1806CbcHeuristicGreedySOS::resetModel(CbcModel * model)
1807{
1808    delete [] originalRhs_;
1809    gutsOfConstructor(model);
1810}
1811// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1812void
1813CbcHeuristicGreedySOS::validate()
1814{
1815    if (model_ && when() < 10) {
1816        if (model_->numberIntegers() !=
1817                model_->numberObjects() && (model_->numberObjects() ||
1818                                            (model_->specialOptions()&1024) == 0)) {
1819            int numberOdd = 0;
1820            for (int i = 0; i < model_->numberObjects(); i++) {
1821                if (!model_->object(i)->canDoHeuristics())
1822                    numberOdd++;
1823            }
1824            if (numberOdd)
1825                setWhen(0);
1826        }
1827        // Only works if coefficients positive and all rows L/G or SOS
1828        OsiSolverInterface * solver = model_->solver();
1829        const double * columnUpper = solver->getColUpper();
1830        const double * columnLower = solver->getColLower();
1831        const double * rowLower = solver->getRowLower();
1832        const double * rowUpper = solver->getRowUpper();
1833
1834        int numberRows = solver->getNumRows();
1835        // Column copy
1836        const double * element = matrix_.getElements();
1837        const int * row = matrix_.getIndices();
1838        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1839        const int * columnLength = matrix_.getVectorLengths();
1840        bool good = true;
1841        assert (originalRhs_);
1842        for (int iRow = 0; iRow < numberRows; iRow++) {
1843          if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
1844            // SOS
1845            originalRhs_[iRow]=-1.0;
1846          } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) {
1847                good = false;
1848          } else if (rowUpper[iRow] < 0.0) {
1849                good = false;
1850          } else if (rowUpper[iRow] < 1.0e10) {
1851            originalRhs_[iRow]=rowUpper[iRow];
1852          } else {
1853            originalRhs_[iRow]=rowLower[iRow];
1854          }
1855        }
1856        int numberColumns = solver->getNumCols();
1857        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1858          if (!columnLength[iColumn])
1859            continue;
1860            if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
1861                good = false;
1862            CoinBigIndex j;
1863            int nSOS=0;
1864            if (!solver->isInteger(iColumn))
1865              good = false;
1866            for (j = columnStart[iColumn];
1867                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1868                if (element[j] < 0.0)
1869                    good = false;
1870                int iRow = row[j];
1871                if (originalRhs_[iRow]==-1.0) {
1872                  if (element[j] != 1.0)
1873                    good = false;
1874                  nSOS++;
1875                }
1876            }
1877            if (nSOS > 1)
1878              good = false;
1879        }
1880        if (!good)
1881            setWhen(0); // switch off
1882    }
1883}
Note: See TracBrowser for help on using the repository browser.