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

Last change on this file since 1433 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 32.0 KB
Line 
1/* $Id: CbcHeuristicGreedy.cpp 1286 2009-11-09 23:33:07Z tkr $ */
2// Copyright (C) 2005, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcStrategy.hpp"
16#include "CbcHeuristicGreedy.hpp"
17#include "CoinSort.hpp"
18#include "CglPreProcess.hpp"
19// Default Constructor
20CbcHeuristicGreedyCover::CbcHeuristicGreedyCover()
21        : CbcHeuristic()
22{
23    // matrix  will automatically be empty
24    originalNumberRows_ = 0;
25    algorithm_ = 0;
26    numberTimes_ = 100;
27}
28
29// Constructor from model
30CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(CbcModel & model)
31        : CbcHeuristic(model)
32{
33    gutsOfConstructor(&model);
34    algorithm_ = 0;
35    numberTimes_ = 100;
36    whereFrom_ = 1;
37}
38
39// Destructor
40CbcHeuristicGreedyCover::~CbcHeuristicGreedyCover ()
41{
42}
43
44// Clone
45CbcHeuristic *
46CbcHeuristicGreedyCover::clone() const
47{
48    return new CbcHeuristicGreedyCover(*this);
49}
50// Guts of constructor from a CbcModel
51void
52CbcHeuristicGreedyCover::gutsOfConstructor(CbcModel * model)
53{
54    model_ = model;
55    // Get a copy of original matrix
56    assert(model->solver());
57    if (model->solver()->getNumRows()) {
58        matrix_ = *model->solver()->getMatrixByCol();
59    }
60    originalNumberRows_ = model->solver()->getNumRows();
61}
62// Create C++ lines to get to current state
63void
64CbcHeuristicGreedyCover::generateCpp( FILE * fp)
65{
66    CbcHeuristicGreedyCover other;
67    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
68    fprintf(fp, "3  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n");
69    CbcHeuristic::generateCpp(fp, "heuristicGreedyCover");
70    if (algorithm_ != other.algorithm_)
71        fprintf(fp, "3  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
72    else
73        fprintf(fp, "4  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
74    if (numberTimes_ != other.numberTimes_)
75        fprintf(fp, "3  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
76    else
77        fprintf(fp, "4  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
78    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyCover);\n");
79}
80
81// Copy constructor
82CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(const CbcHeuristicGreedyCover & rhs)
83        :
84        CbcHeuristic(rhs),
85        matrix_(rhs.matrix_),
86        originalNumberRows_(rhs.originalNumberRows_),
87        algorithm_(rhs.algorithm_),
88        numberTimes_(rhs.numberTimes_)
89{
90}
91
92// Assignment operator
93CbcHeuristicGreedyCover &
94CbcHeuristicGreedyCover::operator=( const CbcHeuristicGreedyCover & rhs)
95{
96    if (this != &rhs) {
97        CbcHeuristic::operator=(rhs);
98        matrix_ = rhs.matrix_;
99        originalNumberRows_ = rhs.originalNumberRows_;
100        algorithm_ = rhs.algorithm_;
101        numberTimes_ = rhs.numberTimes_;
102    }
103    return *this;
104}
105// Returns 1 if solution, 0 if not
106int
107CbcHeuristicGreedyCover::solution(double & solutionValue,
108                                  double * betterSolution)
109{
110    numCouldRun_++;
111    if (!model_)
112        return 0;
113    // See if to do
114    if (!when() || (when() == 1 && model_->phase() != 1))
115        return 0; // switched off
116    if (model_->getNodeCount() > numberTimes_)
117        return 0;
118    // See if at root node
119    bool atRoot = model_->getNodeCount() == 0;
120    int passNumber = model_->getCurrentPassNumber();
121    if (atRoot && passNumber != 1)
122        return 0;
123    OsiSolverInterface * solver = model_->solver();
124    const double * columnLower = solver->getColLower();
125    const double * columnUpper = solver->getColUpper();
126    // And original upper bounds in case we want to use them
127    const double * originalUpper = model_->continuousSolver()->getColUpper();
128    // But not if algorithm says so
129    if ((algorithm_ % 10) == 0)
130        originalUpper = columnUpper;
131    const double * rowLower = solver->getRowLower();
132    const double * solution = solver->getColSolution();
133    const double * objective = solver->getObjCoefficients();
134    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
135    double primalTolerance;
136    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
137
138    // This is number of rows when matrix was passed in
139    int numberRows = originalNumberRows_;
140    if (!numberRows)
141        return 0; // switched off
142
143    numRuns_++;
144    assert (numberRows == matrix_.getNumRows());
145    int iRow, iColumn;
146    double direction = solver->getObjSense();
147    double offset;
148    solver->getDblParam(OsiObjOffset, offset);
149    double newSolutionValue = -offset;
150    int returnCode = 0;
151
152    // Column copy
153    const double * element = matrix_.getElements();
154    const int * row = matrix_.getIndices();
155    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
156    const int * columnLength = matrix_.getVectorLengths();
157
158    // Get solution array for heuristic solution
159    int numberColumns = solver->getNumCols();
160    double * newSolution = new double [numberColumns];
161    double * rowActivity = new double[numberRows];
162    memset(rowActivity, 0, numberRows*sizeof(double));
163    bool allOnes = true;
164    // Get rounded down solution
165    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
166        CoinBigIndex j;
167        double value = solution[iColumn];
168        if (solver->isInteger(iColumn)) {
169            // Round down integer
170            if (fabs(floor(value + 0.5) - value) < integerTolerance) {
171                value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
172            } else {
173                value = CoinMax(floor(value), columnLower[iColumn]);
174            }
175        }
176        // make sure clean
177        value = CoinMin(value, columnUpper[iColumn]);
178        value = CoinMax(value, columnLower[iColumn]);
179        newSolution[iColumn] = value;
180        double cost = direction * objective[iColumn];
181        newSolutionValue += value * cost;
182        for (j = columnStart[iColumn];
183                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
184            int iRow = row[j];
185            rowActivity[iRow] += value * element[j];
186            if (element[j] != 1.0)
187                allOnes = false;
188        }
189    }
190    // See if we round up
191    bool roundup = ((algorithm_ % 100) != 0);
192    if (roundup && allOnes) {
193        // Get rounded up solution
194        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
195            CoinBigIndex j;
196            double value = solution[iColumn];
197            if (solver->isInteger(iColumn)) {
198                // but round up if no activity
199                if (roundup && value >= 0.499999 && !newSolution[iColumn]) {
200                    bool choose = true;
201                    for (j = columnStart[iColumn];
202                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
203                        int iRow = row[j];
204                        if (rowActivity[iRow]) {
205                            choose = false;
206                            break;
207                        }
208                    }
209                    if (choose) {
210                        newSolution[iColumn] = 1.0;
211                        double cost = direction * objective[iColumn];
212                        newSolutionValue += cost;
213                        for (j = columnStart[iColumn];
214                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
215                            int iRow = row[j];
216                            rowActivity[iRow] += 1.0;
217                        }
218                    }
219                }
220            }
221        }
222    }
223    // Get initial list
224    int * which = new int [numberColumns];
225    for (iColumn = 0; iColumn < numberColumns; iColumn++)
226        which[iColumn] = iColumn;
227    int numberLook = numberColumns;
228    // See if we want to perturb more
229    double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
230    // Keep going round until a solution
231    while (true) {
232        // Get column with best ratio
233        int bestColumn = -1;
234        double bestRatio = COIN_DBL_MAX;
235        double bestStepSize = 0.0;
236        int newNumber = 0;
237        for (int jColumn = 0; jColumn < numberLook; jColumn++) {
238            int iColumn = which[jColumn];
239            CoinBigIndex j;
240            double value = newSolution[iColumn];
241            double cost = direction * objective[iColumn];
242            if (solver->isInteger(iColumn)) {
243                // use current upper or original upper
244                if (value + 0.99 < originalUpper[iColumn]) {
245                    double sum = 0.0;
246                    int numberExact = 0;
247                    for (j = columnStart[iColumn];
248                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
249                        int iRow = row[j];
250                        double gap = rowLower[iRow] - rowActivity[iRow];
251                        double elementValue = allOnes ? 1.0 : element[j];
252                        if (gap > 1.0e-7) {
253                            sum += CoinMin(elementValue, gap);
254                            if (fabs(elementValue - gap) < 1.0e-7)
255                                numberExact++;
256                        }
257                    }
258                    // could bias if exact
259                    if (sum > 0.0) {
260                        // add to next time
261                        which[newNumber++] = iColumn;
262                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
263                        // If at root choose first
264                        if (atRoot)
265                            ratio = iColumn;
266                        if (ratio < bestRatio) {
267                            bestRatio = ratio;
268                            bestColumn = iColumn;
269                            bestStepSize = 1.0;
270                        }
271                    }
272                }
273            } else {
274                // continuous
275                if (value < columnUpper[iColumn]) {
276                    // Go through twice - first to get step length
277                    double step = 1.0e50;
278                    for (j = columnStart[iColumn];
279                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
280                        int iRow = row[j];
281                        if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 &&
282                                element[j]*step + rowActivity[iRow] >= rowLower[iRow]) {
283                            step = (rowLower[iRow] - rowActivity[iRow]) / element[j];;
284                        }
285                    }
286                    // now ratio
287                    if (step < 1.0e50) {
288                        // add to next time
289                        which[newNumber++] = iColumn;
290                        assert (step > 0.0);
291                        double sum = 0.0;
292                        for (j = columnStart[iColumn];
293                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
294                            int iRow = row[j];
295                            double newActivity = element[j] * step + rowActivity[iRow];
296                            if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 &&
297                                    newActivity >= rowLower[iRow] - 1.0e-12) {
298                                sum += element[j];
299                            }
300                        }
301                        assert (sum > 0.0);
302                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
303                        if (ratio < bestRatio) {
304                            bestRatio = ratio;
305                            bestColumn = iColumn;
306                            bestStepSize = step;
307                        }
308                    }
309                }
310            }
311        }
312        if (bestColumn < 0)
313            break; // we have finished
314        // Increase chosen column
315        newSolution[bestColumn] += bestStepSize;
316        double cost = direction * objective[bestColumn];
317        newSolutionValue += bestStepSize * cost;
318        for (CoinBigIndex j = columnStart[bestColumn];
319                j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
320            int iRow = row[j];
321            rowActivity[iRow] += bestStepSize * element[j];
322        }
323    }
324    delete [] which;
325    if (newSolutionValue < solutionValue) {
326        // check feasible
327        memset(rowActivity, 0, numberRows*sizeof(double));
328        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
329            CoinBigIndex j;
330            double value = newSolution[iColumn];
331            if (value) {
332                for (j = columnStart[iColumn];
333                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
334                    int iRow = row[j];
335                    rowActivity[iRow] += value * element[j];
336                }
337            }
338        }
339        // check was approximately feasible
340        bool feasible = true;
341        for (iRow = 0; iRow < numberRows; iRow++) {
342            if (rowActivity[iRow] < rowLower[iRow]) {
343                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
344                    feasible = false;
345            }
346        }
347        if (feasible) {
348            // new solution
349            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
350            solutionValue = newSolutionValue;
351            //printf("** Solution of %g found by rounding\n",newSolutionValue);
352            returnCode = 1;
353        } else {
354            // Can easily happen
355            //printf("Debug CbcHeuristicGreedyCover giving bad solution\n");
356        }
357    }
358    delete [] newSolution;
359    delete [] rowActivity;
360    return returnCode;
361}
362// update model
363void CbcHeuristicGreedyCover::setModel(CbcModel * model)
364{
365    gutsOfConstructor(model);
366    validate();
367}
368// Resets stuff if model changes
369void
370CbcHeuristicGreedyCover::resetModel(CbcModel * model)
371{
372    gutsOfConstructor(model);
373}
374// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
375void
376CbcHeuristicGreedyCover::validate()
377{
378    if (model_ && when() < 10) {
379        if (model_->numberIntegers() !=
380                model_->numberObjects() && (model_->numberObjects() ||
381                                            (model_->specialOptions()&1024) == 0)) {
382            int numberOdd = 0;
383            for (int i = 0; i < model_->numberObjects(); i++) {
384                if (!model_->object(i)->canDoHeuristics())
385                    numberOdd++;
386            }
387            if (numberOdd)
388                setWhen(0);
389        }
390        // Only works if costs positive, coefficients positive and all rows G
391        OsiSolverInterface * solver = model_->solver();
392        const double * columnLower = solver->getColLower();
393        const double * rowUpper = solver->getRowUpper();
394        const double * objective = solver->getObjCoefficients();
395        double direction = solver->getObjSense();
396
397        int numberRows = solver->getNumRows();
398        // Column copy
399        const double * element = matrix_.getElements();
400        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
401        const int * columnLength = matrix_.getVectorLengths();
402        bool good = true;
403        for (int iRow = 0; iRow < numberRows; iRow++) {
404            if (rowUpper[iRow] < 1.0e30)
405                good = false;
406        }
407        int numberColumns = solver->getNumCols();
408        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
409            if (objective[iColumn]*direction < 0.0)
410                good = false;
411            if (columnLower[iColumn] < 0.0)
412                good = false;
413            CoinBigIndex j;
414            for (j = columnStart[iColumn];
415                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
416                if (element[j] < 0.0)
417                    good = false;
418            }
419        }
420        if (!good)
421            setWhen(0); // switch off
422    }
423}
424// Default Constructor
425CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality()
426        : CbcHeuristic()
427{
428    // matrix  will automatically be empty
429    fraction_ = 1.0; // no branch and bound
430    originalNumberRows_ = 0;
431    algorithm_ = 0;
432    numberTimes_ = 100;
433    whereFrom_ = 1;
434}
435
436// Constructor from model
437CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel & model)
438        : CbcHeuristic(model)
439{
440    // Get a copy of original matrix
441    gutsOfConstructor(&model);
442    fraction_ = 1.0; // no branch and bound
443    algorithm_ = 0;
444    numberTimes_ = 100;
445    whereFrom_ = 1;
446}
447
448// Destructor
449CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality ()
450{
451}
452
453// Clone
454CbcHeuristic *
455CbcHeuristicGreedyEquality::clone() const
456{
457    return new CbcHeuristicGreedyEquality(*this);
458}
459// Guts of constructor from a CbcModel
460void
461CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel * model)
462{
463    model_ = model;
464    // Get a copy of original matrix
465    assert(model->solver());
466    if (model->solver()->getNumRows()) {
467        matrix_ = *model->solver()->getMatrixByCol();
468    }
469    originalNumberRows_ = model->solver()->getNumRows();
470}
471// Create C++ lines to get to current state
472void
473CbcHeuristicGreedyEquality::generateCpp( FILE * fp)
474{
475    CbcHeuristicGreedyEquality other;
476    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
477    fprintf(fp, "3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
478    CbcHeuristic::generateCpp(fp, "heuristicGreedyEquality");
479    if (algorithm_ != other.algorithm_)
480        fprintf(fp, "3  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
481    else
482        fprintf(fp, "4  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
483    if (fraction_ != other.fraction_)
484        fprintf(fp, "3  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
485    else
486        fprintf(fp, "4  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
487    if (numberTimes_ != other.numberTimes_)
488        fprintf(fp, "3  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
489    else
490        fprintf(fp, "4  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
491    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
492}
493
494// Copy constructor
495CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality & rhs)
496        :
497        CbcHeuristic(rhs),
498        matrix_(rhs.matrix_),
499        fraction_(rhs.fraction_),
500        originalNumberRows_(rhs.originalNumberRows_),
501        algorithm_(rhs.algorithm_),
502        numberTimes_(rhs.numberTimes_)
503{
504}
505
506// Assignment operator
507CbcHeuristicGreedyEquality &
508CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality & rhs)
509{
510    if (this != &rhs) {
511        CbcHeuristic::operator=(rhs);
512        matrix_ = rhs.matrix_;
513        fraction_ = rhs.fraction_;
514        originalNumberRows_ = rhs.originalNumberRows_;
515        algorithm_ = rhs.algorithm_;
516        numberTimes_ = rhs.numberTimes_;
517    }
518    return *this;
519}
520// Returns 1 if solution, 0 if not
521int
522CbcHeuristicGreedyEquality::solution(double & solutionValue,
523                                     double * betterSolution)
524{
525    numCouldRun_++;
526    if (!model_)
527        return 0;
528    // See if to do
529    if (!when() || (when() == 1 && model_->phase() != 1))
530        return 0; // switched off
531    if (model_->getNodeCount() > numberTimes_)
532        return 0;
533    // See if at root node
534    bool atRoot = model_->getNodeCount() == 0;
535    int passNumber = model_->getCurrentPassNumber();
536    if (atRoot && passNumber != 1)
537        return 0;
538    OsiSolverInterface * solver = model_->solver();
539    const double * columnLower = solver->getColLower();
540    const double * columnUpper = solver->getColUpper();
541    // And original upper bounds in case we want to use them
542    const double * originalUpper = model_->continuousSolver()->getColUpper();
543    // But not if algorithm says so
544    if ((algorithm_ % 10) == 0)
545        originalUpper = columnUpper;
546    const double * rowLower = solver->getRowLower();
547    const double * rowUpper = solver->getRowUpper();
548    const double * solution = solver->getColSolution();
549    const double * objective = solver->getObjCoefficients();
550    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
551    double primalTolerance;
552    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
553
554    // This is number of rows when matrix was passed in
555    int numberRows = originalNumberRows_;
556    if (!numberRows)
557        return 0; // switched off
558    numRuns_++;
559
560    assert (numberRows == matrix_.getNumRows());
561    int iRow, iColumn;
562    double direction = solver->getObjSense();
563    double offset;
564    solver->getDblParam(OsiObjOffset, offset);
565    double newSolutionValue = -offset;
566    int returnCode = 0;
567
568    // Column copy
569    const double * element = matrix_.getElements();
570    const int * row = matrix_.getIndices();
571    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
572    const int * columnLength = matrix_.getVectorLengths();
573
574    // Get solution array for heuristic solution
575    int numberColumns = solver->getNumCols();
576    double * newSolution = new double [numberColumns];
577    double * rowActivity = new double[numberRows];
578    memset(rowActivity, 0, numberRows*sizeof(double));
579    double rhsNeeded = 0;
580    for (iRow = 0; iRow < numberRows; iRow++)
581        rhsNeeded += rowUpper[iRow];
582    rhsNeeded *= fraction_;
583    bool allOnes = true;
584    // Get rounded down solution
585    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
586        CoinBigIndex j;
587        double value = solution[iColumn];
588        if (solver->isInteger(iColumn)) {
589            // Round down integer
590            if (fabs(floor(value + 0.5) - value) < integerTolerance) {
591                value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
592            } else {
593                value = CoinMax(floor(value), columnLower[iColumn]);
594            }
595        }
596        // make sure clean
597        value = CoinMin(value, columnUpper[iColumn]);
598        value = CoinMax(value, columnLower[iColumn]);
599        newSolution[iColumn] = value;
600        double cost = direction * objective[iColumn];
601        newSolutionValue += value * cost;
602        for (j = columnStart[iColumn];
603                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
604            int iRow = row[j];
605            rowActivity[iRow] += value * element[j];
606            rhsNeeded -= value * element[j];
607            if (element[j] != 1.0)
608                allOnes = false;
609        }
610    }
611    // See if we round up
612    bool roundup = ((algorithm_ % 100) != 0);
613    if (roundup && allOnes) {
614        // Get rounded up solution
615        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
616            CoinBigIndex j;
617            double value = solution[iColumn];
618            if (solver->isInteger(iColumn)) {
619                // but round up if no activity
620                if (roundup && value >= 0.6 && !newSolution[iColumn]) {
621                    bool choose = true;
622                    for (j = columnStart[iColumn];
623                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
624                        int iRow = row[j];
625                        if (rowActivity[iRow]) {
626                            choose = false;
627                            break;
628                        }
629                    }
630                    if (choose) {
631                        newSolution[iColumn] = 1.0;
632                        double cost = direction * objective[iColumn];
633                        newSolutionValue += cost;
634                        for (j = columnStart[iColumn];
635                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
636                            int iRow = row[j];
637                            rowActivity[iRow] += 1.0;
638                            rhsNeeded -= 1.0;
639                        }
640                    }
641                }
642            }
643        }
644    }
645    // Get initial list
646    int * which = new int [numberColumns];
647    for (iColumn = 0; iColumn < numberColumns; iColumn++)
648        which[iColumn] = iColumn;
649    int numberLook = numberColumns;
650    // See if we want to perturb more
651    double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
652    // Keep going round until a solution
653    while (true) {
654        // Get column with best ratio
655        int bestColumn = -1;
656        double bestRatio = COIN_DBL_MAX;
657        double bestStepSize = 0.0;
658        int newNumber = 0;
659        for (int jColumn = 0; jColumn < numberLook; jColumn++) {
660            int iColumn = which[jColumn];
661            CoinBigIndex j;
662            double value = newSolution[iColumn];
663            double cost = direction * objective[iColumn];
664            if (solver->isInteger(iColumn)) {
665                // use current upper or original upper
666                if (value + 0.9999 < originalUpper[iColumn]) {
667                    double movement = 1.0;
668                    double sum = 0.0;
669                    for (j = columnStart[iColumn];
670                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
671                        int iRow = row[j];
672                        double gap = rowUpper[iRow] - rowActivity[iRow];
673                        double elementValue = allOnes ? 1.0 : element[j];
674                        sum += elementValue;
675                        if (movement*elementValue > gap) {
676                            movement = gap / elementValue;
677                        }
678                    }
679                    if (movement > 0.999999) {
680                        // add to next time
681                        which[newNumber++] = iColumn;
682                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
683                        // If at root
684                        if (atRoot) {
685                            if (fraction_ == 1.0)
686                                ratio = iColumn; // choose first
687                            else
688                                ratio = - solution[iColumn]; // choose largest
689                        }
690                        if (ratio < bestRatio) {
691                            bestRatio = ratio;
692                            bestColumn = iColumn;
693                            bestStepSize = 1.0;
694                        }
695                    }
696                }
697            } else {
698                // continuous
699                if (value < columnUpper[iColumn]) {
700                    double movement = 1.0e50;
701                    double sum = 0.0;
702                    for (j = columnStart[iColumn];
703                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
704                        int iRow = row[j];
705                        if (element[j]*movement + rowActivity[iRow] > rowUpper[iRow]) {
706                            movement = (rowUpper[iRow] - rowActivity[iRow]) / element[j];;
707                        }
708                        sum += element[j];
709                    }
710                    // now ratio
711                    if (movement > 1.0e-7) {
712                        // add to next time
713                        which[newNumber++] = iColumn;
714                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
715                        if (ratio < bestRatio) {
716                            bestRatio = ratio;
717                            bestColumn = iColumn;
718                            bestStepSize = movement;
719                        }
720                    }
721                }
722            }
723        }
724        if (bestColumn < 0)
725            break; // we have finished
726        // Increase chosen column
727        newSolution[bestColumn] += bestStepSize;
728        double cost = direction * objective[bestColumn];
729        newSolutionValue += bestStepSize * cost;
730        for (CoinBigIndex j = columnStart[bestColumn];
731                j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
732            int iRow = row[j];
733            rowActivity[iRow] += bestStepSize * element[j];
734            rhsNeeded -= bestStepSize * element[j];
735        }
736        if (rhsNeeded < 1.0e-8)
737            break;
738    }
739    delete [] which;
740    if (fraction_ < 1.0 && rhsNeeded < 1.0e-8 && newSolutionValue < solutionValue) {
741        // do branch and cut
742        // fix all nonzero
743        OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
744        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
745            if (newSolver->isInteger(iColumn))
746                newSolver->setColLower(iColumn, newSolution[iColumn]);
747        }
748        int returnCode = smallBranchAndBound(newSolver, 200, newSolution, newSolutionValue,
749                                             solutionValue, "CbcHeuristicGreedy");
750        if (returnCode < 0)
751            returnCode = 0; // returned on size
752        if ((returnCode&2) != 0) {
753            // could add cut
754            returnCode &= ~2;
755        }
756        rhsNeeded = 1.0 - returnCode;
757        delete newSolver;
758    }
759    if (newSolutionValue < solutionValue && rhsNeeded < 1.0e-8) {
760        // check feasible
761        memset(rowActivity, 0, numberRows*sizeof(double));
762        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
763            CoinBigIndex j;
764            double value = newSolution[iColumn];
765            if (value) {
766                for (j = columnStart[iColumn];
767                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
768                    int iRow = row[j];
769                    rowActivity[iRow] += value * element[j];
770                }
771            }
772        }
773        // check was approximately feasible
774        bool feasible = true;
775        for (iRow = 0; iRow < numberRows; iRow++) {
776            if (rowActivity[iRow] < rowLower[iRow]) {
777                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
778                    feasible = false;
779            }
780        }
781        if (feasible) {
782            // new solution
783            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
784            solutionValue = newSolutionValue;
785            returnCode = 1;
786        }
787    }
788    delete [] newSolution;
789    delete [] rowActivity;
790    if (atRoot && fraction_ == 1.0) {
791        // try quick search
792        fraction_ = 0.4;
793        int newCode = this->solution(solutionValue, betterSolution);
794        if (newCode)
795            returnCode = 1;
796        fraction_ = 1.0;
797    }
798    return returnCode;
799}
800// update model
801void CbcHeuristicGreedyEquality::setModel(CbcModel * model)
802{
803    gutsOfConstructor(model);
804    validate();
805}
806// Resets stuff if model changes
807void
808CbcHeuristicGreedyEquality::resetModel(CbcModel * model)
809{
810    gutsOfConstructor(model);
811}
812// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
813void
814CbcHeuristicGreedyEquality::validate()
815{
816    if (model_ && when() < 10) {
817        if (model_->numberIntegers() !=
818                model_->numberObjects())
819            setWhen(0);
820        // Only works if costs positive, coefficients positive and all rows E or L
821        // And if values are integer
822        OsiSolverInterface * solver = model_->solver();
823        const double * columnLower = solver->getColLower();
824        const double * rowUpper = solver->getRowUpper();
825        const double * rowLower = solver->getRowLower();
826        const double * objective = solver->getObjCoefficients();
827        double direction = solver->getObjSense();
828
829        int numberRows = solver->getNumRows();
830        // Column copy
831        const double * element = matrix_.getElements();
832        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
833        const int * columnLength = matrix_.getVectorLengths();
834        bool good = true;
835        for (int iRow = 0; iRow < numberRows; iRow++) {
836            if (rowUpper[iRow] > 1.0e30)
837                good = false;
838            if (rowLower[iRow] > 0.0 && rowLower[iRow] != rowUpper[iRow])
839                good = false;
840            if (floor(rowUpper[iRow] + 0.5) != rowUpper[iRow])
841                good = false;
842        }
843        int numberColumns = solver->getNumCols();
844        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
845            if (objective[iColumn]*direction < 0.0)
846                good = false;
847            if (columnLower[iColumn] < 0.0)
848                good = false;
849            CoinBigIndex j;
850            for (j = columnStart[iColumn];
851                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
852                if (element[j] < 0.0)
853                    good = false;
854                if (floor(element[j] + 0.5) != element[j])
855                    good = false;
856            }
857        }
858        if (!good)
859            setWhen(0); // switch off
860    }
861}
862
863
Note: See TracBrowser for help on using the repository browser.