source: branches/sandbox/Cbc/src/CbcHeuristicLocal.cpp @ 1315

Last change on this file since 1315 was 1315, checked in by forrest, 10 years ago

final changes before cleaning

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.4 KB
Line 
1/* $Id: CbcHeuristicLocal.cpp 1315 2009-11-28 11:09:56Z forrest $ */
2// Copyright (C) 2002, 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 "CbcMessage.hpp"
16#include "CbcHeuristicLocal.hpp"
17#include "CbcBranchActual.hpp"
18#include "CbcStrategy.hpp"
19#include "CglPreProcess.hpp"
20
21// Default Constructor
22CbcHeuristicLocal::CbcHeuristicLocal()
23        : CbcHeuristic()
24{
25    numberSolutions_ = 0;
26    swap_ = 0;
27    used_ = NULL;
28}
29
30// Constructor with model - assumed before cuts
31
32CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model)
33        : CbcHeuristic(model)
34{
35    numberSolutions_ = 0;
36    swap_ = 0;
37    // Get a copy of original matrix
38    assert(model.solver());
39    if (model.solver()->getNumRows()) {
40        matrix_ = *model.solver()->getMatrixByCol();
41    }
42    int numberColumns = model.solver()->getNumCols();
43    used_ = new int[numberColumns];
44    memset(used_, 0, numberColumns*sizeof(int));
45}
46
47// Destructor
48CbcHeuristicLocal::~CbcHeuristicLocal ()
49{
50    delete [] used_;
51}
52
53// Clone
54CbcHeuristic *
55CbcHeuristicLocal::clone() const
56{
57    return new CbcHeuristicLocal(*this);
58}
59// Create C++ lines to get to current state
60void
61CbcHeuristicLocal::generateCpp( FILE * fp)
62{
63    CbcHeuristicLocal other;
64    fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
65    fprintf(fp, "3  CbcHeuristicLocal heuristicLocal(*cbcModel);\n");
66    CbcHeuristic::generateCpp(fp, "heuristicLocal");
67    if (swap_ != other.swap_)
68        fprintf(fp, "3  heuristicLocal.setSearchType(%d);\n", swap_);
69    else
70        fprintf(fp, "4  heuristicLocal.setSearchType(%d);\n", swap_);
71    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicLocal);\n");
72}
73
74// Copy constructor
75CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs)
76        :
77        CbcHeuristic(rhs),
78        matrix_(rhs.matrix_),
79        numberSolutions_(rhs.numberSolutions_),
80        swap_(rhs.swap_)
81{
82    if (model_ && rhs.used_) {
83        int numberColumns = model_->solver()->getNumCols();
84        used_ = CoinCopyOfArray(rhs.used_, numberColumns);
85    } else {
86        used_ = NULL;
87    }
88}
89
90// Assignment operator
91CbcHeuristicLocal &
92CbcHeuristicLocal::operator=( const CbcHeuristicLocal & rhs)
93{
94    if (this != &rhs) {
95        CbcHeuristic::operator=(rhs);
96        matrix_ = rhs.matrix_;
97        numberSolutions_ = rhs.numberSolutions_;
98        swap_ = rhs.swap_;
99        delete [] used_;
100        if (model_ && rhs.used_) {
101            int numberColumns = model_->solver()->getNumCols();
102            used_ = CoinCopyOfArray(rhs.used_, numberColumns);
103        } else {
104            used_ = NULL;
105        }
106    }
107    return *this;
108}
109
110// Resets stuff if model changes
111void
112CbcHeuristicLocal::resetModel(CbcModel * /*model*/)
113{
114    //CbcHeuristic::resetModel(model);
115    delete [] used_;
116    if (model_ && used_) {
117        int numberColumns = model_->solver()->getNumCols();
118        used_ = new int[numberColumns];
119        memset(used_, 0, numberColumns*sizeof(int));
120    } else {
121        used_ = NULL;
122    }
123}
124// This version fixes stuff and does IP
125int
126CbcHeuristicLocal::solutionFix(double & objectiveValue,
127                               double * newSolution,
128                               const int * /*keep*/)
129{
130    numCouldRun_++;
131    // See if to do
132    if (!when() || (when() == 1 && model_->phase() != 1))
133        return 0; // switched off
134    // Don't do if it was this heuristic which found solution!
135    if (this == model_->lastHeuristic())
136        return 0;
137    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
138    const double * colLower = newSolver->getColLower();
139    //const double * colUpper = newSolver->getColUpper();
140
141    int numberIntegers = model_->numberIntegers();
142    const int * integerVariable = model_->integerVariable();
143
144    int i;
145    int nFix = 0;
146    for (i = 0; i < numberIntegers; i++) {
147        int iColumn = integerVariable[i];
148        const OsiObject * object = model_->object(i);
149        // get original bounds
150        double originalLower;
151        double originalUpper;
152        getIntegerInformation( object, originalLower, originalUpper);
153        newSolver->setColLower(iColumn, CoinMax(colLower[iColumn], originalLower));
154        if (!used_[iColumn]) {
155            newSolver->setColUpper(iColumn, colLower[iColumn]);
156            nFix++;
157        }
158    }
159    int returnCode = 0;
160#ifdef CLP_INVESTIGATE2
161    printf("Fixing %d out of %d (%d continuous)\n",
162           nFix, numberIntegers, newSolver->getNumCols() - numberIntegers);
163#endif
164    if (nFix*10 <= numberIntegers) {
165        // see if we can fix more
166        int * which = new int [2*(numberIntegers-nFix)];
167        int * sort = which + (numberIntegers - nFix);
168        int n = 0;
169        for (i = 0; i < numberIntegers; i++) {
170            int iColumn = integerVariable[i];
171            if (used_[iColumn]) {
172                which[n] = iColumn;
173                sort[n++] = used_[iColumn];
174            }
175        }
176        CoinSort_2(sort, sort + n, which);
177        // only half fixed in total
178        n = CoinMin(n, numberIntegers / 2 - nFix);
179        int allow = CoinMax(numberSolutions_ - 2, sort[0]);
180        int nFix2 = 0;
181        for (i = 0; i < n; i++) {
182            int iColumn = integerVariable[i];
183            if (used_[iColumn] <= allow) {
184                newSolver->setColUpper(iColumn, colLower[iColumn]);
185                nFix2++;
186            } else {
187                break;
188            }
189        }
190        delete [] which;
191        nFix += nFix2;
192        printf("Number fixed increased from %d to %d\n",
193               nFix - nFix2, nFix);
194    }
195    if (nFix*10 > numberIntegers) {
196        returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue,
197                                         objectiveValue, "CbcHeuristicLocal");
198        if (returnCode < 0) {
199            returnCode = 0; // returned on size
200            int numberColumns = newSolver->getNumCols();
201            int numberContinuous = numberColumns - numberIntegers;
202            if (numberContinuous > 2*numberIntegers &&
203                    nFix*10 < numberColumns) {
204#define LOCAL_FIX_CONTINUOUS
205#ifdef LOCAL_FIX_CONTINUOUS
206                //const double * colUpper = newSolver->getColUpper();
207                const double * colLower = newSolver->getColLower();
208                int nAtLb = 0;
209                //double sumDj=0.0;
210                const double * dj = newSolver->getReducedCost();
211                double direction = newSolver->getObjSense();
212                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
213                    if (!newSolver->isInteger(iColumn)) {
214                        if (!used_[iColumn]) {
215                            //double djValue = dj[iColumn]*direction;
216                            nAtLb++;
217                            //sumDj += djValue;
218                        }
219                    }
220                }
221                if (nAtLb) {
222                    // fix some continuous
223                    double * sort = new double[nAtLb];
224                    int * which = new int [nAtLb];
225                    //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
226                    int nFix2 = 0;
227                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
228                        if (!newSolver->isInteger(iColumn)) {
229                            if (!used_[iColumn]) {
230                                double djValue = dj[iColumn] * direction;
231                                if (djValue > 1.0e-6) {
232                                    sort[nFix2] = -djValue;
233                                    which[nFix2++] = iColumn;
234                                }
235                            }
236                        }
237                    }
238                    CoinSort_2(sort, sort + nFix2, which);
239                    int divisor = 2;
240                    nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
241                    for (int i = 0; i < nFix2; i++) {
242                        int iColumn = which[i];
243                        newSolver->setColUpper(iColumn, colLower[iColumn]);
244                    }
245                    delete [] sort;
246                    delete [] which;
247#ifdef CLP_INVESTIGATE2
248                    printf("%d integers have zero value, and %d continuous fixed at lb\n",
249                           nFix, nFix2);
250#endif
251                    returnCode = smallBranchAndBound(newSolver,
252                                                     numberNodes_, newSolution,
253                                                     objectiveValue,
254                                                     objectiveValue, "CbcHeuristicLocal");
255                    if (returnCode < 0)
256                        returnCode = 0; // returned on size
257                }
258#endif
259            }
260        }
261    }
262    if ((returnCode&2) != 0) {
263        // could add cut
264        returnCode &= ~2;
265    }
266
267    delete newSolver;
268    return returnCode;
269}
270/*
271  First tries setting a variable to better value.  If feasible then
272  tries setting others.  If not feasible then tries swaps
273  Returns 1 if solution, 0 if not */
274int
275CbcHeuristicLocal::solution(double & solutionValue,
276                            double * betterSolution)
277{
278
279    numCouldRun_++;
280    if (numberSolutions_ == model_->getSolutionCount())
281        return 0;
282    numberSolutions_ = model_->getSolutionCount();
283    if ((model_->getNumCols() > 100000 && model_->getNumCols() >
284            10*model_->getNumRows()) || numberSolutions_ <= 1)
285        return 0; // probably not worth it
286    // worth trying
287
288    OsiSolverInterface * solver = model_->solver();
289    const double * rowLower = solver->getRowLower();
290    const double * rowUpper = solver->getRowUpper();
291    const double * solution = model_->bestSolution();
292    if (!solution)
293        return 0; // No solution found yet
294    const double * objective = solver->getObjCoefficients();
295    double primalTolerance;
296    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
297
298    int numberRows = matrix_.getNumRows();
299
300    int numberIntegers = model_->numberIntegers();
301    const int * integerVariable = model_->integerVariable();
302
303    int i;
304    double direction = solver->getObjSense();
305    double newSolutionValue = model_->getObjValue() * direction;
306    int returnCode = 0;
307    numRuns_++;
308    // Column copy
309    const double * element = matrix_.getElements();
310    const int * row = matrix_.getIndices();
311    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
312    const int * columnLength = matrix_.getVectorLengths();
313
314    // Get solution array for heuristic solution
315    int numberColumns = solver->getNumCols();
316    double * newSolution = new double [numberColumns];
317    memcpy(newSolution, solution, numberColumns*sizeof(double));
318#ifdef LOCAL_FIX_CONTINUOUS
319    // mark continuous used
320    const double * columnLower = solver->getColLower();
321    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
322        if (!solver->isInteger(iColumn)) {
323            if (solution[iColumn] > columnLower[iColumn] + 1.0e-8)
324                used_[iColumn] = numberSolutions_;
325        }
326    }
327#endif
328
329    // way is 1 if down possible, 2 if up possible, 3 if both possible
330    char * way = new char[numberIntegers];
331    // corrected costs
332    double * cost = new double[numberIntegers];
333    // for array to mark infeasible rows after iColumn branch
334    char * mark = new char[numberRows];
335    memset(mark, 0, numberRows);
336    // space to save values so we don't introduce rounding errors
337    double * save = new double[numberRows];
338
339    // clean solution
340    for (i = 0; i < numberIntegers; i++) {
341        int iColumn = integerVariable[i];
342        const OsiObject * object = model_->object(i);
343        // get original bounds
344        double originalLower;
345        double originalUpper;
346        getIntegerInformation( object, originalLower, originalUpper);
347        double value = newSolution[iColumn];
348        if (value < originalLower) {
349            value = originalLower;
350            newSolution[iColumn] = value;
351        } else if (value > originalUpper) {
352            value = originalUpper;
353            newSolution[iColumn] = value;
354        }
355        double nearest = floor(value + 0.5);
356        //assert(fabs(value-nearest)<10.0*primalTolerance);
357        value = nearest;
358        newSolution[iColumn] = nearest;
359        // if away from lower bound mark that fact
360        if (nearest > originalLower) {
361            used_[iColumn] = numberSolutions_;
362        }
363        cost[i] = direction * objective[iColumn];
364        int iway = 0;
365
366        if (value > originalLower + 0.5)
367            iway = 1;
368        if (value < originalUpper - 0.5)
369            iway |= 2;
370        way[i] = static_cast<char>(iway);
371    }
372    // get row activities
373    double * rowActivity = new double[numberRows];
374    memset(rowActivity, 0, numberRows*sizeof(double));
375
376    for (i = 0; i < numberColumns; i++) {
377        int j;
378        double value = newSolution[i];
379        if (value) {
380            for (j = columnStart[i];
381                    j < columnStart[i] + columnLength[i]; j++) {
382                int iRow = row[j];
383                rowActivity[iRow] += value * element[j];
384            }
385        }
386    }
387    // check was feasible - if not adjust (cleaning may move)
388    // if very infeasible then give up
389    bool tryHeuristic = true;
390    for (i = 0; i < numberRows; i++) {
391        if (rowActivity[i] < rowLower[i]) {
392            if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance)
393                tryHeuristic = false;
394            rowActivity[i] = rowLower[i];
395        } else if (rowActivity[i] > rowUpper[i]) {
396            if (rowActivity[i] < rowUpper[i] + 10.0*primalTolerance)
397                tryHeuristic = false;
398            rowActivity[i] = rowUpper[i];
399        }
400    }
401    // Switch off if may take too long
402    if (model_->getNumCols() > 10000 && model_->getNumCols() >
403            10*model_->getNumRows())
404        tryHeuristic = false;
405    if (tryHeuristic) {
406
407        // best change in objective
408        double bestChange = 0.0;
409
410        for (i = 0; i < numberIntegers; i++) {
411            int iColumn = integerVariable[i];
412
413            double objectiveCoefficient = cost[i];
414            int k;
415            int j;
416            int goodK = -1;
417            int wayK = -1, wayI = -1;
418            if ((way[i]&1) != 0) {
419                int numberInfeasible = 0;
420                // save row activities and adjust
421                for (j = columnStart[iColumn];
422                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
423                    int iRow = row[j];
424                    save[iRow] = rowActivity[iRow];
425                    rowActivity[iRow] -= element[j];
426                    if (rowActivity[iRow] < rowLower[iRow] - primalTolerance ||
427                            rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
428                        // mark row
429                        mark[iRow] = 1;
430                        numberInfeasible++;
431                    }
432                }
433                // try down
434                for (k = i + 1; k < numberIntegers; k++) {
435                    if ((way[k]&1) != 0) {
436                        // try down
437                        if (-objectiveCoefficient - cost[k] < bestChange) {
438                            // see if feasible down
439                            bool good = true;
440                            int numberMarked = 0;
441                            int kColumn = integerVariable[k];
442                            for (j = columnStart[kColumn];
443                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
444                                int iRow = row[j];
445                                double newValue = rowActivity[iRow] - element[j];
446                                if (newValue < rowLower[iRow] - primalTolerance ||
447                                        newValue > rowUpper[iRow] + primalTolerance) {
448                                    good = false;
449                                    break;
450                                } else if (mark[iRow]) {
451                                    // made feasible
452                                    numberMarked++;
453                                }
454                            }
455                            if (good && numberMarked == numberInfeasible) {
456                                // better solution
457                                goodK = k;
458                                wayK = -1;
459                                wayI = -1;
460                                bestChange = -objectiveCoefficient - cost[k];
461                            }
462                        }
463                    }
464                    if ((way[k]&2) != 0) {
465                        // try up
466                        if (-objectiveCoefficient + cost[k] < bestChange) {
467                            // see if feasible up
468                            bool good = true;
469                            int numberMarked = 0;
470                            int kColumn = integerVariable[k];
471                            for (j = columnStart[kColumn];
472                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
473                                int iRow = row[j];
474                                double newValue = rowActivity[iRow] + element[j];
475                                if (newValue < rowLower[iRow] - primalTolerance ||
476                                        newValue > rowUpper[iRow] + primalTolerance) {
477                                    good = false;
478                                    break;
479                                } else if (mark[iRow]) {
480                                    // made feasible
481                                    numberMarked++;
482                                }
483                            }
484                            if (good && numberMarked == numberInfeasible) {
485                                // better solution
486                                goodK = k;
487                                wayK = 1;
488                                wayI = -1;
489                                bestChange = -objectiveCoefficient + cost[k];
490                            }
491                        }
492                    }
493                }
494                // restore row activities
495                for (j = columnStart[iColumn];
496                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
497                    int iRow = row[j];
498                    rowActivity[iRow] = save[iRow];
499                    mark[iRow] = 0;
500                }
501            }
502            if ((way[i]&2) != 0) {
503                int numberInfeasible = 0;
504                // save row activities and adjust
505                for (j = columnStart[iColumn];
506                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
507                    int iRow = row[j];
508                    save[iRow] = rowActivity[iRow];
509                    rowActivity[iRow] += element[j];
510                    if (rowActivity[iRow] < rowLower[iRow] - primalTolerance ||
511                            rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
512                        // mark row
513                        mark[iRow] = 1;
514                        numberInfeasible++;
515                    }
516                }
517                // try up
518                for (k = i + 1; k < numberIntegers; k++) {
519                    if ((way[k]&1) != 0) {
520                        // try down
521                        if (objectiveCoefficient - cost[k] < bestChange) {
522                            // see if feasible down
523                            bool good = true;
524                            int numberMarked = 0;
525                            int kColumn = integerVariable[k];
526                            for (j = columnStart[kColumn];
527                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
528                                int iRow = row[j];
529                                double newValue = rowActivity[iRow] - element[j];
530                                if (newValue < rowLower[iRow] - primalTolerance ||
531                                        newValue > rowUpper[iRow] + primalTolerance) {
532                                    good = false;
533                                    break;
534                                } else if (mark[iRow]) {
535                                    // made feasible
536                                    numberMarked++;
537                                }
538                            }
539                            if (good && numberMarked == numberInfeasible) {
540                                // better solution
541                                goodK = k;
542                                wayK = -1;
543                                wayI = 1;
544                                bestChange = objectiveCoefficient - cost[k];
545                            }
546                        }
547                    }
548                    if ((way[k]&2) != 0) {
549                        // try up
550                        if (objectiveCoefficient + cost[k] < bestChange) {
551                            // see if feasible up
552                            bool good = true;
553                            int numberMarked = 0;
554                            int kColumn = integerVariable[k];
555                            for (j = columnStart[kColumn];
556                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
557                                int iRow = row[j];
558                                double newValue = rowActivity[iRow] + element[j];
559                                if (newValue < rowLower[iRow] - primalTolerance ||
560                                        newValue > rowUpper[iRow] + primalTolerance) {
561                                    good = false;
562                                    break;
563                                } else if (mark[iRow]) {
564                                    // made feasible
565                                    numberMarked++;
566                                }
567                            }
568                            if (good && numberMarked == numberInfeasible) {
569                                // better solution
570                                goodK = k;
571                                wayK = 1;
572                                wayI = 1;
573                                bestChange = objectiveCoefficient + cost[k];
574                            }
575                        }
576                    }
577                }
578                // restore row activities
579                for (j = columnStart[iColumn];
580                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
581                    int iRow = row[j];
582                    rowActivity[iRow] = save[iRow];
583                    mark[iRow] = 0;
584                }
585            }
586            if (goodK >= 0) {
587                // we found something - update solution
588                for (j = columnStart[iColumn];
589                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
590                    int iRow = row[j];
591                    rowActivity[iRow]  += wayI * element[j];
592                }
593                newSolution[iColumn] += wayI;
594                int kColumn = integerVariable[goodK];
595                for (j = columnStart[kColumn];
596                        j < columnStart[kColumn] + columnLength[kColumn]; j++) {
597                    int iRow = row[j];
598                    rowActivity[iRow]  += wayK * element[j];
599                }
600                newSolution[kColumn] += wayK;
601                // See if k can go further ?
602                const OsiObject * object = model_->object(goodK);
603                // get original bounds
604                double originalLower;
605                double originalUpper;
606                getIntegerInformation( object, originalLower, originalUpper);
607
608                double value = newSolution[kColumn];
609                int iway = 0;
610
611                if (value > originalLower + 0.5)
612                    iway = 1;
613                if (value < originalUpper - 0.5)
614                    iway |= 2;
615                way[goodK] = static_cast<char>(iway);
616            }
617        }
618        if (bestChange + newSolutionValue < solutionValue) {
619            // paranoid check
620            memset(rowActivity, 0, numberRows*sizeof(double));
621
622            for (i = 0; i < numberColumns; i++) {
623                int j;
624                double value = newSolution[i];
625                if (value) {
626                    for (j = columnStart[i];
627                            j < columnStart[i] + columnLength[i]; j++) {
628                        int iRow = row[j];
629                        rowActivity[iRow] += value * element[j];
630                    }
631                }
632            }
633            int numberBad = 0;
634            double sumBad = 0.0;
635            // check was approximately feasible
636            for (i = 0; i < numberRows; i++) {
637                if (rowActivity[i] < rowLower[i]) {
638                    sumBad += rowLower[i] - rowActivity[i];
639                    if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance)
640                        numberBad++;
641                } else if (rowActivity[i] > rowUpper[i]) {
642                    sumBad += rowUpper[i] - rowActivity[i];
643                    if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance)
644                        numberBad++;
645                }
646            }
647            if (!numberBad) {
648                for (i = 0; i < numberIntegers; i++) {
649                    int iColumn = integerVariable[i];
650                    const OsiObject * object = model_->object(i);
651                    // get original bounds
652                    double originalLower;
653                    double originalUpper;
654                    getIntegerInformation( object, originalLower, originalUpper);
655
656                    double value = newSolution[iColumn];
657                    // if away from lower bound mark that fact
658                    if (value > originalLower) {
659                        used_[iColumn] = numberSolutions_;
660                    }
661                }
662                // new solution
663                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
664                CoinWarmStartBasis * basis =
665                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
666                if (basis) {
667                    model_->setBestSolutionBasis(* basis);
668                    delete basis;
669                }
670                returnCode = 1;
671                solutionValue = newSolutionValue + bestChange;
672            } else {
673                // bad solution - should not happen so debug if see message
674                printf("Local search got bad solution with %d infeasibilities summing to %g\n",
675                       numberBad, sumBad);
676            }
677        }
678    }
679    delete [] newSolution;
680    delete [] rowActivity;
681    delete [] way;
682    delete [] cost;
683    delete [] save;
684    delete [] mark;
685    if (numberSolutions_ > 1 && swap_ == 1) {
686        // try merge
687        int returnCode2 = solutionFix( solutionValue, betterSolution, NULL);
688        if (returnCode2)
689            returnCode = 1;
690    }
691    return returnCode;
692}
693// update model
694void CbcHeuristicLocal::setModel(CbcModel * model)
695{
696    model_ = model;
697    // Get a copy of original matrix
698    assert(model_->solver());
699    if (model_->solver()->getNumRows()) {
700        matrix_ = *model_->solver()->getMatrixByCol();
701    }
702    delete [] used_;
703    int numberColumns = model->solver()->getNumCols();
704    used_ = new int[numberColumns];
705    memset(used_, 0, numberColumns*sizeof(int));
706}
707// Default Constructor
708CbcHeuristicNaive::CbcHeuristicNaive()
709        : CbcHeuristic()
710{
711    large_ = 1.0e6;
712}
713
714// Constructor with model - assumed before cuts
715
716CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model)
717        : CbcHeuristic(model)
718{
719    large_ = 1.0e6;
720}
721
722// Destructor
723CbcHeuristicNaive::~CbcHeuristicNaive ()
724{
725}
726
727// Clone
728CbcHeuristic *
729CbcHeuristicNaive::clone() const
730{
731    return new CbcHeuristicNaive(*this);
732}
733// Create C++ lines to get to current state
734void
735CbcHeuristicNaive::generateCpp( FILE * fp)
736{
737    CbcHeuristicNaive other;
738    fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
739    fprintf(fp, "3  CbcHeuristicNaive naive(*cbcModel);\n");
740    CbcHeuristic::generateCpp(fp, "naive");
741    if (large_ != other.large_)
742        fprintf(fp, "3  naive.setLarge(%g);\n", large_);
743    else
744        fprintf(fp, "4  naive.setLarge(%g);\n", large_);
745    fprintf(fp, "3  cbcModel->addHeuristic(&naive);\n");
746}
747
748// Copy constructor
749CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs)
750        :
751        CbcHeuristic(rhs),
752        large_(rhs.large_)
753{
754}
755
756// Assignment operator
757CbcHeuristicNaive &
758CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs)
759{
760    if (this != &rhs) {
761        CbcHeuristic::operator=(rhs);
762        large_ = rhs.large_;
763    }
764    return *this;
765}
766
767// Resets stuff if model changes
768void
769CbcHeuristicNaive::resetModel(CbcModel * model)
770{
771    CbcHeuristic::resetModel(model);
772}
773int
774CbcHeuristicNaive::solution(double & solutionValue,
775                            double * betterSolution)
776{
777    numCouldRun_++;
778    // See if to do
779    bool atRoot = model_->getNodeCount() == 0;
780    int passNumber = model_->getCurrentPassNumber();
781    if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber != 1)
782        return 0; // switched off
783    // Don't do if it was this heuristic which found solution!
784    if (this == model_->lastHeuristic())
785        return 0;
786    numRuns_++;
787    double cutoff;
788    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
789    double direction = model_->solver()->getObjSense();
790    cutoff *= direction;
791    cutoff = CoinMin(cutoff, solutionValue);
792    OsiSolverInterface * solver = model_->continuousSolver();
793    if (!solver)
794        solver = model_->solver();
795    const double * colLower = solver->getColLower();
796    const double * colUpper = solver->getColUpper();
797    const double * objective = solver->getObjCoefficients();
798
799    int numberColumns = model_->getNumCols();
800    int numberIntegers = model_->numberIntegers();
801    const int * integerVariable = model_->integerVariable();
802
803    int i;
804    bool solutionFound = false;
805    CoinWarmStartBasis saveBasis;
806    CoinWarmStartBasis * basis =
807        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
808    if (basis) {
809        saveBasis = * basis;
810        delete basis;
811    }
812    // First just fix all integers as close to zero as possible
813    OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
814    for (i = 0; i < numberIntegers; i++) {
815        int iColumn = integerVariable[i];
816        double lower = colLower[iColumn];
817        double upper = colUpper[iColumn];
818        double value;
819        if (lower > 0.0)
820            value = lower;
821        else if (upper < 0.0)
822            value = upper;
823        else
824            value = 0.0;
825        newSolver->setColLower(iColumn, value);
826        newSolver->setColUpper(iColumn, value);
827    }
828    newSolver->initialSolve();
829    if (newSolver->isProvenOptimal()) {
830        double solValue = newSolver->getObjValue() * direction ;
831        if (solValue < cutoff) {
832            // we have a solution
833            solutionFound = true;
834            solutionValue = solValue;
835            memcpy(betterSolution, newSolver->getColSolution(),
836                   numberColumns*sizeof(double));
837            printf("Naive fixing close to zero gave solution of %g\n", solutionValue);
838            cutoff = solValue - model_->getCutoffIncrement();
839        }
840    }
841    // Now fix all integers as close to zero if zero or large cost
842    int nFix = 0;
843    for (i = 0; i < numberIntegers; i++) {
844        int iColumn = integerVariable[i];
845        double lower = colLower[iColumn];
846        double upper = colUpper[iColumn];
847        double value;
848        if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) {
849            nFix++;
850            if (lower > 0.0)
851                value = lower;
852            else if (upper < 0.0)
853                value = upper;
854            else
855                value = 0.0;
856            newSolver->setColLower(iColumn, value);
857            newSolver->setColUpper(iColumn, value);
858        } else {
859            // set back to original
860            newSolver->setColLower(iColumn, lower);
861            newSolver->setColUpper(iColumn, upper);
862        }
863    }
864    const double * solution = solver->getColSolution();
865    if (nFix) {
866        newSolver->setWarmStart(&saveBasis);
867        newSolver->setColSolution(solution);
868        newSolver->initialSolve();
869        if (newSolver->isProvenOptimal()) {
870            double solValue = newSolver->getObjValue() * direction ;
871            if (solValue < cutoff) {
872                // try branch and bound
873                double * newSolution = new double [numberColumns];
874                printf("%d fixed after fixing costs\n", nFix);
875                int returnCode = smallBranchAndBound(newSolver,
876                                                     numberNodes_, newSolution,
877                                                     solutionValue,
878                                                     solutionValue, "CbcHeuristicNaive1");
879                if (returnCode < 0)
880                    returnCode = 0; // returned on size
881                if ((returnCode&2) != 0) {
882                    // could add cut
883                    returnCode &= ~2;
884                }
885                if (returnCode == 1) {
886                    // solution
887                    solutionFound = true;
888                    memcpy(betterSolution, newSolution,
889                           numberColumns*sizeof(double));
890                    printf("Naive fixing zeros gave solution of %g\n", solutionValue);
891                    cutoff = solutionValue - model_->getCutoffIncrement();
892                }
893                delete [] newSolution;
894            }
895        }
896    }
897#if 1
898    newSolver->setObjSense(-direction); // maximize
899    newSolver->setWarmStart(&saveBasis);
900    newSolver->setColSolution(solution);
901    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
902        double value = solution[iColumn];
903        double lower = colLower[iColumn];
904        double upper = colUpper[iColumn];
905        double newLower;
906        double newUpper;
907        if (newSolver->isInteger(iColumn)) {
908            newLower = CoinMax(lower, floor(value) - 2.0);
909            newUpper = CoinMin(upper, ceil(value) + 2.0);
910        } else {
911            newLower = CoinMax(lower, value - 1.0e5);
912            newUpper = CoinMin(upper, value + 1.0e-5);
913        }
914        newSolver->setColLower(iColumn, newLower);
915        newSolver->setColUpper(iColumn, newUpper);
916    }
917    newSolver->initialSolve();
918    if (newSolver->isProvenOptimal()) {
919        double solValue = newSolver->getObjValue() * direction ;
920        if (solValue < cutoff) {
921            nFix = 0;
922            newSolver->setObjSense(direction); // correct direction
923            //const double * thisSolution = newSolver->getColSolution();
924            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
925                double value = solution[iColumn];
926                double lower = colLower[iColumn];
927                double upper = colUpper[iColumn];
928                double newLower = lower;
929                double newUpper = upper;
930                if (newSolver->isInteger(iColumn)) {
931                    if (value < lower + 1.0e-6) {
932                        nFix++;
933                        newUpper = lower;
934                    } else if (value > upper - 1.0e-6) {
935                        nFix++;
936                        newLower = upper;
937                    } else {
938                        newLower = CoinMax(lower, floor(value) - 2.0);
939                        newUpper = CoinMin(upper, ceil(value) + 2.0);
940                    }
941                }
942                newSolver->setColLower(iColumn, newLower);
943                newSolver->setColUpper(iColumn, newUpper);
944            }
945            // try branch and bound
946            double * newSolution = new double [numberColumns];
947            printf("%d fixed after maximizing\n", nFix);
948            int returnCode = smallBranchAndBound(newSolver,
949                                                 numberNodes_, newSolution,
950                                                 solutionValue,
951                                                 solutionValue, "CbcHeuristicNaive1");
952            if (returnCode < 0)
953                returnCode = 0; // returned on size
954            if ((returnCode&2) != 0) {
955                // could add cut
956                returnCode &= ~2;
957            }
958            if (returnCode == 1) {
959                // solution
960                solutionFound = true;
961                memcpy(betterSolution, newSolution,
962                       numberColumns*sizeof(double));
963                printf("Naive maximizing gave solution of %g\n", solutionValue);
964                cutoff = solutionValue - model_->getCutoffIncrement();
965            }
966            delete [] newSolution;
967        }
968    }
969#endif
970    delete newSolver;
971    return solutionFound ? 1 : 0;
972}
973// update model
974void CbcHeuristicNaive::setModel(CbcModel * model)
975{
976    model_ = model;
977}
978// Default Constructor
979CbcHeuristicCrossover::CbcHeuristicCrossover()
980        : CbcHeuristic(),
981        numberSolutions_(0),
982        useNumber_(3)
983{
984    setWhen(1);
985}
986
987// Constructor with model - assumed before cuts
988
989CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model)
990        : CbcHeuristic(model),
991        numberSolutions_(0),
992        useNumber_(3)
993{
994    setWhen(1);
995    for (int i = 0; i < 10; i++)
996        random_[i] = model.randomNumberGenerator()->randomDouble();
997}
998
999// Destructor
1000CbcHeuristicCrossover::~CbcHeuristicCrossover ()
1001{
1002}
1003
1004// Clone
1005CbcHeuristic *
1006CbcHeuristicCrossover::clone() const
1007{
1008    return new CbcHeuristicCrossover(*this);
1009}
1010// Create C++ lines to get to current state
1011void
1012CbcHeuristicCrossover::generateCpp( FILE * fp)
1013{
1014    CbcHeuristicCrossover other;
1015    fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
1016    fprintf(fp, "3  CbcHeuristicCrossover crossover(*cbcModel);\n");
1017    CbcHeuristic::generateCpp(fp, "crossover");
1018    if (useNumber_ != other.useNumber_)
1019        fprintf(fp, "3  crossover.setNumberSolutions(%d);\n", useNumber_);
1020    else
1021        fprintf(fp, "4  crossover.setNumberSolutions(%d);\n", useNumber_);
1022    fprintf(fp, "3  cbcModel->addHeuristic(&crossover);\n");
1023}
1024
1025// Copy constructor
1026CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs)
1027        :
1028        CbcHeuristic(rhs),
1029        attempts_(rhs.attempts_),
1030        numberSolutions_(rhs.numberSolutions_),
1031        useNumber_(rhs.useNumber_)
1032{
1033    memcpy(random_, rhs.random_, 10*sizeof(double));
1034}
1035
1036// Assignment operator
1037CbcHeuristicCrossover &
1038CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs)
1039{
1040    if (this != &rhs) {
1041        CbcHeuristic::operator=(rhs);
1042        useNumber_ = rhs.useNumber_;
1043        attempts_ = rhs.attempts_;
1044        numberSolutions_ = rhs.numberSolutions_;
1045        memcpy(random_, rhs.random_, 10*sizeof(double));
1046    }
1047    return *this;
1048}
1049
1050// Resets stuff if model changes
1051void
1052CbcHeuristicCrossover::resetModel(CbcModel * model)
1053{
1054    CbcHeuristic::resetModel(model);
1055}
1056int
1057CbcHeuristicCrossover::solution(double & solutionValue,
1058                                double * betterSolution)
1059{
1060    if (when_ == 0)
1061        return 0;
1062    numCouldRun_++;
1063    bool useBest = (numberSolutions_ != model_->getSolutionCount());
1064    if (!useBest && (when_ % 10) == 1)
1065        return 0;
1066    numberSolutions_ = model_->getSolutionCount();
1067    OsiSolverInterface * continuousSolver = model_->continuousSolver();
1068    int useNumber = CoinMin(model_->numberSavedSolutions(), useNumber_);
1069    if (useNumber < 2 || !continuousSolver)
1070        return 0;
1071    // Fix later
1072    if (!useBest)
1073        abort();
1074    numRuns_++;
1075    double cutoff;
1076    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1077    double direction = model_->solver()->getObjSense();
1078    cutoff *= direction;
1079    cutoff = CoinMin(cutoff, solutionValue);
1080    OsiSolverInterface * solver = cloneBut(2);
1081    // But reset bounds
1082    solver->setColLower(continuousSolver->getColLower());
1083    solver->setColUpper(continuousSolver->getColUpper());
1084    int numberColumns = solver->getNumCols();
1085    // Fixed
1086    double * fixed = new double [numberColumns];
1087    for (int i = 0; i < numberColumns; i++)
1088        fixed[i] = -COIN_DBL_MAX;
1089    int whichSolution[10];
1090    for (int i = 0; i < useNumber; i++)
1091        whichSolution[i] = i;
1092    for (int i = 0; i < useNumber; i++) {
1093        int k = whichSolution[i];
1094        const double * solution = model_->savedSolution(k);
1095        for (int j = 0; j < numberColumns; j++) {
1096            if (solver->isInteger(j)) {
1097                if (fixed[j] == -COIN_DBL_MAX)
1098                    fixed[j] = floor(solution[j] + 0.5);
1099                else if (fabs(fixed[j] - solution[j]) > 1.0e-7)
1100                    fixed[j] = COIN_DBL_MAX;
1101            }
1102        }
1103    }
1104    const double * colLower = solver->getColLower();
1105    for (int i = 0; i < numberColumns; i++) {
1106        if (solver->isInteger(i)) {
1107            double value = fixed[i];
1108            if (value != COIN_DBL_MAX) {
1109                if (when_ < 10) {
1110                    solver->setColLower(i, value);
1111                    solver->setColUpper(i, value);
1112                } else if (value == colLower[i]) {
1113                    solver->setColUpper(i, value);
1114                }
1115            }
1116        }
1117    }
1118    int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution,
1119                                         solutionValue,
1120                                         solutionValue, "CbcHeuristicCrossover");
1121    if (returnCode < 0)
1122        returnCode = 0; // returned on size
1123    if ((returnCode&2) != 0) {
1124        // could add cut
1125        returnCode &= ~2;
1126    }
1127
1128    delete solver;
1129    return returnCode;
1130}
1131// update model
1132void CbcHeuristicCrossover::setModel(CbcModel * model)
1133{
1134    model_ = model;
1135    if (model) {
1136        for (int i = 0; i < 10; i++)
1137            random_[i] = model->randomNumberGenerator()->randomDouble();
1138    }
1139}
1140
1141
Note: See TracBrowser for help on using the repository browser.