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

Last change on this file since 1359 was 1359, checked in by forrest, 12 years ago

out printf

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