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

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