source: trunk/Cbc/src/CbcHeuristicLocal.cpp @ 1573

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.6 KB
Line 
1/* $Id: CbcHeuristicLocal.cpp 1573 2011-01-05 01:12:36Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10#include <cassert>
11#include <cstdlib>
12#include <cmath>
13#include <cfloat>
14
15#include "OsiSolverInterface.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcHeuristicLocal.hpp"
19#include "CbcBranchActual.hpp"
20#include "CbcStrategy.hpp"
21#include "CglPreProcess.hpp"
22
23// Default Constructor
24CbcHeuristicLocal::CbcHeuristicLocal()
25        : CbcHeuristic()
26{
27    numberSolutions_ = 0;
28    swap_ = 0;
29    used_ = NULL;
30}
31
32// Constructor with model - assumed before cuts
33
34CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model)
35        : CbcHeuristic(model)
36{
37    numberSolutions_ = 0;
38    swap_ = 0;
39    // Get a copy of original matrix
40    assert(model.solver());
41    if (model.solver()->getNumRows()) {
42        matrix_ = *model.solver()->getMatrixByCol();
43    }
44    int numberColumns = model.solver()->getNumCols();
45    used_ = new int[numberColumns];
46    memset(used_, 0, numberColumns*sizeof(int));
47}
48
49// Destructor
50CbcHeuristicLocal::~CbcHeuristicLocal ()
51{
52    delete [] used_;
53}
54
55// Clone
56CbcHeuristic *
57CbcHeuristicLocal::clone() const
58{
59    return new CbcHeuristicLocal(*this);
60}
61// Create C++ lines to get to current state
62void
63CbcHeuristicLocal::generateCpp( FILE * fp)
64{
65    CbcHeuristicLocal other;
66    fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
67    fprintf(fp, "3  CbcHeuristicLocal heuristicLocal(*cbcModel);\n");
68    CbcHeuristic::generateCpp(fp, "heuristicLocal");
69    if (swap_ != other.swap_)
70        fprintf(fp, "3  heuristicLocal.setSearchType(%d);\n", swap_);
71    else
72        fprintf(fp, "4  heuristicLocal.setSearchType(%d);\n", swap_);
73    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicLocal);\n");
74}
75
76// Copy constructor
77CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs)
78        :
79        CbcHeuristic(rhs),
80        matrix_(rhs.matrix_),
81        numberSolutions_(rhs.numberSolutions_),
82        swap_(rhs.swap_)
83{
84    if (model_ && rhs.used_) {
85        int numberColumns = model_->solver()->getNumCols();
86        used_ = CoinCopyOfArray(rhs.used_, numberColumns);
87    } else {
88        used_ = NULL;
89    }
90}
91
92// Assignment operator
93CbcHeuristicLocal &
94CbcHeuristicLocal::operator=( const CbcHeuristicLocal & rhs)
95{
96    if (this != &rhs) {
97        CbcHeuristic::operator=(rhs);
98        matrix_ = rhs.matrix_;
99        numberSolutions_ = rhs.numberSolutions_;
100        swap_ = rhs.swap_;
101        delete [] used_;
102        if (model_ && rhs.used_) {
103            int numberColumns = model_->solver()->getNumCols();
104            used_ = CoinCopyOfArray(rhs.used_, numberColumns);
105        } else {
106            used_ = NULL;
107        }
108    }
109    return *this;
110}
111
112// Resets stuff if model changes
113void
114CbcHeuristicLocal::resetModel(CbcModel * /*model*/)
115{
116    //CbcHeuristic::resetModel(model);
117    delete [] used_;
118    if (model_ && used_) {
119        int numberColumns = model_->solver()->getNumCols();
120        used_ = new int[numberColumns];
121        memset(used_, 0, numberColumns*sizeof(int));
122    } else {
123        used_ = NULL;
124    }
125}
126/*
127  Run a mini-BaB search after fixing all variables not marked as used by
128  solution(). (See comments there for semantics.)
129
130  Return values are:
131    1: smallBranchAndBound found a solution
132    0: everything else
133
134  The degree of overload as return codes from smallBranchAndBound are folded
135  into 0 is such that it's impossible to distinguish return codes that really
136  require attention from a simple `nothing of interest'.
137*/
138// This version fixes stuff and does IP
139int
140CbcHeuristicLocal::solutionFix(double & objectiveValue,
141                               double * newSolution,
142                               const int * /*keep*/)
143{
144/*
145  If when is set to off (0), or set to root (1) and we're not at the root,
146  return. If this heuristic discovered the current solution, don't continue.
147*/
148
149    numCouldRun_++;
150    // See if to do
151    if (!when() || (when() == 1 && model_->phase() != 1))
152        return 0; // switched off
153    // Don't do if it was this heuristic which found solution!
154    if (this == model_->lastHeuristic())
155        return 0;
156/*
157  Load up a new solver with the solution.
158
159  Why continuousSolver(), as opposed to solver()?
160*/
161    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
162    const double * colLower = newSolver->getColLower();
163    //const double * colUpper = newSolver->getColUpper();
164
165    int numberIntegers = model_->numberIntegers();
166    const int * integerVariable = model_->integerVariable();
167/*
168  The net effect here is that anything that hasn't moved from its lower bound
169  will be fixed at lower bound.
170
171  See comments in solution() w.r.t. asymmetric treatment of upper and lower
172  bounds.
173*/
174
175    int i;
176    int nFix = 0;
177    for (i = 0; i < numberIntegers; i++) {
178        int iColumn = integerVariable[i];
179        const OsiObject * object = model_->object(i);
180        // get original bounds
181        double originalLower;
182        double originalUpper;
183        getIntegerInformation( object, originalLower, originalUpper);
184        newSolver->setColLower(iColumn, CoinMax(colLower[iColumn], originalLower));
185        if (!used_[iColumn]) {
186            newSolver->setColUpper(iColumn, colLower[iColumn]);
187            nFix++;
188        }
189    }
190/*
191  Try a `small' branch-and-bound search. The notion here is that we've fixed a
192  lot of variables and reduced the amount of `free' problem to a point where a
193  small BaB search will suffice to fully explore the remaining problem. This
194  routine will execute integer presolve, then call branchAndBound to do the
195  actual search.
196*/
197    int returnCode = 0;
198#ifdef CLP_INVESTIGATE2
199    printf("Fixing %d out of %d (%d continuous)\n",
200           nFix, numberIntegers, newSolver->getNumCols() - numberIntegers);
201#endif
202    if (nFix*10 <= numberIntegers) {
203        // see if we can fix more
204        int * which = new int [2*(numberIntegers-nFix)];
205        int * sort = which + (numberIntegers - nFix);
206        int n = 0;
207        for (i = 0; i < numberIntegers; i++) {
208            int iColumn = integerVariable[i];
209            if (used_[iColumn]) {
210                which[n] = iColumn;
211                sort[n++] = used_[iColumn];
212            }
213        }
214        CoinSort_2(sort, sort + n, which);
215        // only half fixed in total
216        n = CoinMin(n, numberIntegers / 2 - nFix);
217        int allow = CoinMax(numberSolutions_ - 2, sort[0]);
218        int nFix2 = 0;
219        for (i = 0; i < n; i++) {
220            int iColumn = integerVariable[i];
221            if (used_[iColumn] <= allow) {
222                newSolver->setColUpper(iColumn, colLower[iColumn]);
223                nFix2++;
224            } else {
225                break;
226            }
227        }
228        delete [] which;
229        nFix += nFix2;
230#ifdef CLP_INVESTIGATE2
231        printf("Number fixed increased from %d to %d\n",
232               nFix - nFix2, nFix);
233#endif
234    }
235    if (nFix*10 > numberIntegers) {
236        returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue,
237                                         objectiveValue, "CbcHeuristicLocal");
238 /*
239  -2 is return due to user event, and -1 is overloaded with what look to be
240  two contradictory meanings.
241*/
242       if (returnCode < 0) {
243            returnCode = 0; // returned on size
244            int numberColumns = newSolver->getNumCols();
245            int numberContinuous = numberColumns - numberIntegers;
246            if (numberContinuous > 2*numberIntegers &&
247                    nFix*10 < numberColumns) {
248#define LOCAL_FIX_CONTINUOUS
249#ifdef LOCAL_FIX_CONTINUOUS
250                //const double * colUpper = newSolver->getColUpper();
251                const double * colLower = newSolver->getColLower();
252                int nAtLb = 0;
253                //double sumDj=0.0;
254                const double * dj = newSolver->getReducedCost();
255                double direction = newSolver->getObjSense();
256                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
257                    if (!newSolver->isInteger(iColumn)) {
258                        if (!used_[iColumn]) {
259                            //double djValue = dj[iColumn]*direction;
260                            nAtLb++;
261                            //sumDj += djValue;
262                        }
263                    }
264                }
265                if (nAtLb) {
266                    // fix some continuous
267                    double * sort = new double[nAtLb];
268                    int * which = new int [nAtLb];
269                    //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
270                    int nFix2 = 0;
271                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
272                        if (!newSolver->isInteger(iColumn)) {
273                            if (!used_[iColumn]) {
274                                double djValue = dj[iColumn] * direction;
275                                if (djValue > 1.0e-6) {
276                                    sort[nFix2] = -djValue;
277                                    which[nFix2++] = iColumn;
278                                }
279                            }
280                        }
281                    }
282                    CoinSort_2(sort, sort + nFix2, which);
283                    int divisor = 2;
284                    nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
285                    for (int i = 0; i < nFix2; i++) {
286                        int iColumn = which[i];
287                        newSolver->setColUpper(iColumn, colLower[iColumn]);
288                    }
289                    delete [] sort;
290                    delete [] which;
291#ifdef CLP_INVESTIGATE2
292                    printf("%d integers have zero value, and %d continuous fixed at lb\n",
293                           nFix, nFix2);
294#endif
295                    returnCode = smallBranchAndBound(newSolver,
296                                                     numberNodes_, newSolution,
297                                                     objectiveValue,
298                                                     objectiveValue, "CbcHeuristicLocal");
299                    if (returnCode < 0)
300                        returnCode = 0; // returned on size
301                }
302#endif
303            }
304        }
305    }
306/*
307  If the result is complete exploration with a solution (3) or proven
308  infeasibility (2), we could generate a cut (the AI folks would call it a
309  nogood) to prevent us from going down this route in the future.
310*/
311    if ((returnCode&2) != 0) {
312        // could add cut
313        returnCode &= ~2;
314    }
315
316    delete newSolver;
317    return returnCode;
318}
319/*
320  First tries setting a variable to better value.  If feasible then
321  tries setting others.  If not feasible then tries swaps
322  Returns 1 if solution, 0 if not
323  The main body of this routine implements an O((q^2)/2) brute force search
324  around the current solution, for q = number of integer variables. Call this
325  the inc/dec heuristic.  For each integer variable x<i>, first decrement the
326  value. Then, for integer variables x<i+1>, ..., x<q-1>, try increment and
327  decrement. If one of these permutations produces a better solution,
328  remember it.  Then repeat, with x<i> incremented. If we find a better
329  solution, update our notion of current solution and continue.
330
331  The net effect is a greedy walk: As each improving pair is found, the
332  current solution is updated and the search continues from this updated
333  solution.
334
335  Way down at the end, we call solutionFix, which will create a drastically
336  restricted problem based on variables marked as used, then do mini-BaC on
337  the restricted problem. This can occur even if we don't try the inc/dec
338  heuristic. This would be more obvious if the inc/dec heuristic were broken
339  out as a separate routine and solutionFix had a name that reflected where
340  it was headed.
341
342  The return code of 0 is grossly overloaded, because it maps to a return
343  code of 0 from solutionFix, which is itself grossly overloaded. See
344  comments in solutionFix and in CbcHeuristic::smallBranchAndBound.
345  */
346int
347CbcHeuristicLocal::solution(double & solutionValue,
348                            double * betterSolution)
349{
350/*
351  Execute only if a new solution has been discovered since the last time we
352  were called.
353*/
354
355    numCouldRun_++;
356    if (numberSolutions_ == model_->getSolutionCount())
357        return 0;
358    numberSolutions_ = model_->getSolutionCount();
359/*
360  Exclude long (column), thin (row) systems.
361
362  Given the n^2 nature of the search, more than 100,000 columns could get
363  expensive. But I don't yet see the rationale for the second part of the
364  condition (cols > 10*rows). And cost is proportional to number of integer
365  variables --- shouldn't we use that?
366
367  Why wait until we have more than one solution?
368*/
369    if ((model_->getNumCols() > 100000 && model_->getNumCols() >
370            10*model_->getNumRows()) || numberSolutions_ <= 1)
371        return 0; // probably not worth it
372    // worth trying
373
374    OsiSolverInterface * solver = model_->solver();
375    const double * rowLower = solver->getRowLower();
376    const double * rowUpper = solver->getRowUpper();
377    const double * solution = model_->bestSolution();
378/*
379  Shouldn't this test be redundant if we've already checked that
380  numberSolutions_ > 1? Stronger: shouldn't this be an assertion?
381*/
382    if (!solution)
383        return 0; // No solution found yet
384    const double * objective = solver->getObjCoefficients();
385    double primalTolerance;
386    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
387
388    int numberRows = matrix_.getNumRows();
389
390    int numberIntegers = model_->numberIntegers();
391    const int * integerVariable = model_->integerVariable();
392
393    int i;
394    double direction = solver->getObjSense();
395    double newSolutionValue = model_->getObjValue() * direction;
396    int returnCode = 0;
397    numRuns_++;
398    // Column copy
399    const double * element = matrix_.getElements();
400    const int * row = matrix_.getIndices();
401    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
402    const int * columnLength = matrix_.getVectorLengths();
403
404    // Get solution array for heuristic solution
405    int numberColumns = solver->getNumCols();
406    double * newSolution = new double [numberColumns];
407    memcpy(newSolution, solution, numberColumns*sizeof(double));
408#ifdef LOCAL_FIX_CONTINUOUS
409    // mark continuous used
410    const double * columnLower = solver->getColLower();
411    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
412        if (!solver->isInteger(iColumn)) {
413            if (solution[iColumn] > columnLower[iColumn] + 1.0e-8)
414                used_[iColumn] = numberSolutions_;
415        }
416    }
417#endif
418
419    // way is 1 if down possible, 2 if up possible, 3 if both possible
420    char * way = new char[numberIntegers];
421    // corrected costs
422    double * cost = new double[numberIntegers];
423    // for array to mark infeasible rows after iColumn branch
424    char * mark = new char[numberRows];
425    memset(mark, 0, numberRows);
426    // space to save values so we don't introduce rounding errors
427    double * save = new double[numberRows];
428/*
429  Force variables within their original bounds, then to the nearest integer.
430  Overall, we seem to be prepared to cope with noninteger bounds. Is this
431  necessary? Seems like we'd be better off to force the bounds to integrality
432  as part of preprocessing.  More generally, why do we need to do this? This
433  solution should have been cleaned and checked when it was accepted as a
434  solution!
435
436  Once the value is set, decide whether we can move up or down.
437
438  The only place that used_ is used is in solutionFix; if a variable is not
439  flagged as used, it will be fixed (at lower bound). Why the asymmetric
440  treatment? This makes some sense for binary variables (for which there are
441  only two options). But for general integer variables, why not make a similar
442  test against the original upper bound?
443*/
444
445    // clean solution
446    for (i = 0; i < numberIntegers; i++) {
447        int iColumn = integerVariable[i];
448        const OsiObject * object = model_->object(i);
449        // get original bounds
450        double originalLower;
451        double originalUpper;
452        getIntegerInformation( object, originalLower, originalUpper);
453        double value = newSolution[iColumn];
454        if (value < originalLower) {
455            value = originalLower;
456            newSolution[iColumn] = value;
457        } else if (value > originalUpper) {
458            value = originalUpper;
459            newSolution[iColumn] = value;
460        }
461        double nearest = floor(value + 0.5);
462        //assert(fabs(value-nearest)<10.0*primalTolerance);
463        value = nearest;
464        newSolution[iColumn] = nearest;
465        // if away from lower bound mark that fact
466        if (nearest > originalLower) {
467            used_[iColumn] = numberSolutions_;
468        }
469        cost[i] = direction * objective[iColumn];
470/*
471  Given previous computation we're checking that value is at least 1 away
472  from the original bounds.
473*/
474        int iway = 0;
475
476        if (value > originalLower + 0.5)
477            iway = 1;
478        if (value < originalUpper - 0.5)
479            iway |= 2;
480        way[i] = static_cast<char>(iway);
481    }
482/*
483  Calculate lhs of each constraint for groomed solution.
484*/
485    // get row activities
486    double * rowActivity = new double[numberRows];
487    memset(rowActivity, 0, numberRows*sizeof(double));
488
489    for (i = 0; i < numberColumns; i++) {
490        int j;
491        double value = newSolution[i];
492        if (value) {
493            for (j = columnStart[i];
494                    j < columnStart[i] + columnLength[i]; j++) {
495                int iRow = row[j];
496                rowActivity[iRow] += value * element[j];
497            }
498        }
499    }
500/*
501  Check that constraints are satisfied. For small infeasibility, force the
502  activity within bound. Again, why is this necessary if the current solution
503  was accepted as a valid solution?
504
505  Why are we scanning past the first unacceptable constraint?
506*/
507    // check was feasible - if not adjust (cleaning may move)
508    // if very infeasible then give up
509    bool tryHeuristic = true;
510    for (i = 0; i < numberRows; i++) {
511        if (rowActivity[i] < rowLower[i]) {
512            if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance)
513                tryHeuristic = false;
514            rowActivity[i] = rowLower[i];
515        } else if (rowActivity[i] > rowUpper[i]) {
516            if (rowActivity[i] < rowUpper[i] + 10.0*primalTolerance)
517                tryHeuristic = false;
518            rowActivity[i] = rowUpper[i];
519        }
520    }
521/*
522  This bit of code is not quite totally redundant: it'll bail at 10,000
523  instead of 100,000. Potentially we can do a lot of work to get here, only
524  to abandon it.
525*/
526    // Switch off if may take too long
527    if (model_->getNumCols() > 10000 && model_->getNumCols() >
528            10*model_->getNumRows())
529        tryHeuristic = false;
530/*
531  Try the inc/dec heuristic?
532*/
533    if (tryHeuristic) {
534
535        // best change in objective
536        double bestChange = 0.0;
537/*
538  Outer loop to walk integer variables. Call the current variable x<i>. At the
539  end of this loop, bestChange will contain the best (negative) change in the
540  objective for any single pair.
541
542  The trouble is, we're limited to monotonically increasing improvement.
543  Suppose we discover an improvement of 10 for some pair. If, later in the
544  search, we discover an improvement of 9 for some other pair, we will not use
545  it. That seems wasteful.
546*/
547
548        for (i = 0; i < numberIntegers; i++) {
549            int iColumn = integerVariable[i];
550
551            double objectiveCoefficient = cost[i];
552            int k;
553            int j;
554            int goodK = -1;
555            int wayK = -1, wayI = -1;
556/*
557  Try decrementing x<i>.
558*/
559            if ((way[i]&1) != 0) {
560                int numberInfeasible = 0;
561/*
562  Adjust row activities where x<i> has a nonzero coefficient. Save the old
563  values for restoration. Mark any rows that become infeasible as a result
564  of the decrement.
565*/
566                // save row activities and adjust
567                for (j = columnStart[iColumn];
568                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
569                    int iRow = row[j];
570                    save[iRow] = rowActivity[iRow];
571                    rowActivity[iRow] -= element[j];
572                    if (rowActivity[iRow] < rowLower[iRow] - primalTolerance ||
573                            rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
574                        // mark row
575                        mark[iRow] = 1;
576                        numberInfeasible++;
577                    }
578                }
579  /*
580  Run through the remaining integer variables. Try increment and decrement on
581  each one. If the potential objective change is better than anything we've
582  seen so far, do a full evaluation of x<k> in that direction.  If we can
583  repair all infeasibilities introduced by pushing x<i> down, we have a
584  winner. Remember the best variable, and the direction for x<i> and x<k>.
585*/
586              // try down
587                for (k = i + 1; k < numberIntegers; k++) {
588                    if ((way[k]&1) != 0) {
589                        // try down
590                        if (-objectiveCoefficient - cost[k] < bestChange) {
591                            // see if feasible down
592                            bool good = true;
593                            int numberMarked = 0;
594                            int kColumn = integerVariable[k];
595                            for (j = columnStart[kColumn];
596                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
597                                int iRow = row[j];
598                                double newValue = rowActivity[iRow] - element[j];
599                                if (newValue < rowLower[iRow] - primalTolerance ||
600                                        newValue > rowUpper[iRow] + primalTolerance) {
601                                    good = false;
602                                    break;
603                                } else if (mark[iRow]) {
604                                    // made feasible
605                                    numberMarked++;
606                                }
607                            }
608                            if (good && numberMarked == numberInfeasible) {
609                                // better solution
610                                goodK = k;
611                                wayK = -1;
612                                wayI = -1;
613                                bestChange = -objectiveCoefficient - cost[k];
614                            }
615                        }
616                    }
617                    if ((way[k]&2) != 0) {
618                        // try up
619                        if (-objectiveCoefficient + cost[k] < bestChange) {
620                            // see if feasible up
621                            bool good = true;
622                            int numberMarked = 0;
623                            int kColumn = integerVariable[k];
624                            for (j = columnStart[kColumn];
625                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
626                                int iRow = row[j];
627                                double newValue = rowActivity[iRow] + element[j];
628                                if (newValue < rowLower[iRow] - primalTolerance ||
629                                        newValue > rowUpper[iRow] + primalTolerance) {
630                                    good = false;
631                                    break;
632                                } else if (mark[iRow]) {
633                                    // made feasible
634                                    numberMarked++;
635                                }
636                            }
637                            if (good && numberMarked == numberInfeasible) {
638                                // better solution
639                                goodK = k;
640                                wayK = 1;
641                                wayI = -1;
642                                bestChange = -objectiveCoefficient + cost[k];
643                            }
644                        }
645                    }
646                }
647/*
648  Remove effect of decrementing x<i> by restoring original lhs values.
649*/
650                // restore row activities
651                for (j = columnStart[iColumn];
652                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
653                    int iRow = row[j];
654                    rowActivity[iRow] = save[iRow];
655                    mark[iRow] = 0;
656                }
657            }
658/*
659  Try to increment x<i>. Actions as for decrement.
660*/
661            if ((way[i]&2) != 0) {
662                int numberInfeasible = 0;
663                // save row activities and adjust
664                for (j = columnStart[iColumn];
665                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
666                    int iRow = row[j];
667                    save[iRow] = rowActivity[iRow];
668                    rowActivity[iRow] += element[j];
669                    if (rowActivity[iRow] < rowLower[iRow] - primalTolerance ||
670                            rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
671                        // mark row
672                        mark[iRow] = 1;
673                        numberInfeasible++;
674                    }
675                }
676                // try up
677                for (k = i + 1; k < numberIntegers; k++) {
678                    if ((way[k]&1) != 0) {
679                        // try down
680                        if (objectiveCoefficient - cost[k] < bestChange) {
681                            // see if feasible down
682                            bool good = true;
683                            int numberMarked = 0;
684                            int kColumn = integerVariable[k];
685                            for (j = columnStart[kColumn];
686                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
687                                int iRow = row[j];
688                                double newValue = rowActivity[iRow] - element[j];
689                                if (newValue < rowLower[iRow] - primalTolerance ||
690                                        newValue > rowUpper[iRow] + primalTolerance) {
691                                    good = false;
692                                    break;
693                                } else if (mark[iRow]) {
694                                    // made feasible
695                                    numberMarked++;
696                                }
697                            }
698                            if (good && numberMarked == numberInfeasible) {
699                                // better solution
700                                goodK = k;
701                                wayK = -1;
702                                wayI = 1;
703                                bestChange = objectiveCoefficient - cost[k];
704                            }
705                        }
706                    }
707                    if ((way[k]&2) != 0) {
708                        // try up
709                        if (objectiveCoefficient + cost[k] < bestChange) {
710                            // see if feasible up
711                            bool good = true;
712                            int numberMarked = 0;
713                            int kColumn = integerVariable[k];
714                            for (j = columnStart[kColumn];
715                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
716                                int iRow = row[j];
717                                double newValue = rowActivity[iRow] + element[j];
718                                if (newValue < rowLower[iRow] - primalTolerance ||
719                                        newValue > rowUpper[iRow] + primalTolerance) {
720                                    good = false;
721                                    break;
722                                } else if (mark[iRow]) {
723                                    // made feasible
724                                    numberMarked++;
725                                }
726                            }
727                            if (good && numberMarked == numberInfeasible) {
728                                // better solution
729                                goodK = k;
730                                wayK = 1;
731                                wayI = 1;
732                                bestChange = objectiveCoefficient + cost[k];
733                            }
734                        }
735                    }
736                }
737                // restore row activities
738                for (j = columnStart[iColumn];
739                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
740                    int iRow = row[j];
741                    rowActivity[iRow] = save[iRow];
742                    mark[iRow] = 0;
743                }
744            }
745/*
746  We've found a pair x<i> and x<k> which produce a better solution. Update our
747  notion of current solution to match.
748
749  Why does this not update newSolutionValue?
750*/
751            if (goodK >= 0) {
752                // we found something - update solution
753                for (j = columnStart[iColumn];
754                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
755                    int iRow = row[j];
756                    rowActivity[iRow]  += wayI * element[j];
757                }
758                newSolution[iColumn] += wayI;
759                int kColumn = integerVariable[goodK];
760                for (j = columnStart[kColumn];
761                        j < columnStart[kColumn] + columnLength[kColumn]; j++) {
762                    int iRow = row[j];
763                    rowActivity[iRow]  += wayK * element[j];
764                }
765                newSolution[kColumn] += wayK;
766/*
767  Adjust motion range for x<k>. We may have banged up against a bound with that
768  last move.
769*/
770               // See if k can go further ?
771                const OsiObject * object = model_->object(goodK);
772                // get original bounds
773                double originalLower;
774                double originalUpper;
775                getIntegerInformation( object, originalLower, originalUpper);
776
777                double value = newSolution[kColumn];
778                int iway = 0;
779
780                if (value > originalLower + 0.5)
781                    iway = 1;
782                if (value < originalUpper - 0.5)
783                    iway |= 2;
784                way[goodK] = static_cast<char>(iway);
785            }
786        }
787/*
788  End of loop to try increment/decrement of integer variables.
789
790  newSolutionValue does not necessarily match the current newSolution, and
791  bestChange simply reflects the best single change. Still, that's sufficient
792  to indicate that there's been at least one change. Check that we really do
793  have a valid solution.
794*/
795        if (bestChange + newSolutionValue < solutionValue) {
796            // paranoid check
797            memset(rowActivity, 0, numberRows*sizeof(double));
798
799            for (i = 0; i < numberColumns; i++) {
800                int j;
801                double value = newSolution[i];
802                if (value) {
803                    for (j = columnStart[i];
804                            j < columnStart[i] + columnLength[i]; j++) {
805                        int iRow = row[j];
806                        rowActivity[iRow] += value * element[j];
807                    }
808                }
809            }
810            int numberBad = 0;
811            double sumBad = 0.0;
812            // check was approximately feasible
813            for (i = 0; i < numberRows; i++) {
814                if (rowActivity[i] < rowLower[i]) {
815                    sumBad += rowLower[i] - rowActivity[i];
816                    if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance)
817                        numberBad++;
818                } else if (rowActivity[i] > rowUpper[i]) {
819                    sumBad += rowUpper[i] - rowActivity[i];
820                    if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance)
821                        numberBad++;
822                }
823            }
824            if (!numberBad) {
825                for (i = 0; i < numberIntegers; i++) {
826                    int iColumn = integerVariable[i];
827                    const OsiObject * object = model_->object(i);
828                    // get original bounds
829                    double originalLower;
830                    double originalUpper;
831                    getIntegerInformation( object, originalLower, originalUpper);
832
833                    double value = newSolution[iColumn];
834                    // if away from lower bound mark that fact
835                    if (value > originalLower) {
836                        used_[iColumn] = numberSolutions_;
837                    }
838                }
839/*
840  Copy the solution to the array returned to the client. Grab a basis from
841  the solver (which, if it exists, is almost certainly infeasible, but it
842  should be ok for a dual start). The value returned as solutionValue is
843  conservative because of handling of newSolutionValue and bestChange, as
844  described above.
845*/
846                // new solution
847                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
848                CoinWarmStartBasis * basis =
849                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
850                if (basis) {
851                    model_->setBestSolutionBasis(* basis);
852                    delete basis;
853                }
854                returnCode = 1;
855                solutionValue = newSolutionValue + bestChange;
856            } else {
857                // bad solution - should not happen so debug if see message
858                printf("Local search got bad solution with %d infeasibilities summing to %g\n",
859                       numberBad, sumBad);
860            }
861        }
862    }
863/*
864  We're done. Clean up.
865*/
866    delete [] newSolution;
867    delete [] rowActivity;
868    delete [] way;
869    delete [] cost;
870    delete [] save;
871    delete [] mark;
872/*
873  Do we want to try swapping values between solutions?
874  swap_ is set elsewhere; it's not adjusted during heuristic execution.
875
876  Again, redundant test. We shouldn't be here if numberSolutions_ = 1.
877*/
878    if (numberSolutions_ > 1 && swap_ == 1) {
879        // try merge
880        int returnCode2 = solutionFix( solutionValue, betterSolution, NULL);
881        if (returnCode2)
882            returnCode = 1;
883    }
884    return returnCode;
885}
886// update model
887void CbcHeuristicLocal::setModel(CbcModel * model)
888{
889    model_ = model;
890    // Get a copy of original matrix
891    assert(model_->solver());
892    if (model_->solver()->getNumRows()) {
893        matrix_ = *model_->solver()->getMatrixByCol();
894    }
895    delete [] used_;
896    int numberColumns = model->solver()->getNumCols();
897    used_ = new int[numberColumns];
898    memset(used_, 0, numberColumns*sizeof(int));
899}
900// Default Constructor
901CbcHeuristicNaive::CbcHeuristicNaive()
902        : CbcHeuristic()
903{
904    large_ = 1.0e6;
905}
906
907// Constructor with model - assumed before cuts
908
909CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model)
910        : CbcHeuristic(model)
911{
912    large_ = 1.0e6;
913}
914
915// Destructor
916CbcHeuristicNaive::~CbcHeuristicNaive ()
917{
918}
919
920// Clone
921CbcHeuristic *
922CbcHeuristicNaive::clone() const
923{
924    return new CbcHeuristicNaive(*this);
925}
926// Create C++ lines to get to current state
927void
928CbcHeuristicNaive::generateCpp( FILE * fp)
929{
930    CbcHeuristicNaive other;
931    fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
932    fprintf(fp, "3  CbcHeuristicNaive naive(*cbcModel);\n");
933    CbcHeuristic::generateCpp(fp, "naive");
934    if (large_ != other.large_)
935        fprintf(fp, "3  naive.setLarge(%g);\n", large_);
936    else
937        fprintf(fp, "4  naive.setLarge(%g);\n", large_);
938    fprintf(fp, "3  cbcModel->addHeuristic(&naive);\n");
939}
940
941// Copy constructor
942CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs)
943        :
944        CbcHeuristic(rhs),
945        large_(rhs.large_)
946{
947}
948
949// Assignment operator
950CbcHeuristicNaive &
951CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs)
952{
953    if (this != &rhs) {
954        CbcHeuristic::operator=(rhs);
955        large_ = rhs.large_;
956    }
957    return *this;
958}
959
960// Resets stuff if model changes
961void
962CbcHeuristicNaive::resetModel(CbcModel * model)
963{
964    CbcHeuristic::resetModel(model);
965}
966int
967CbcHeuristicNaive::solution(double & solutionValue,
968                            double * betterSolution)
969{
970    numCouldRun_++;
971    // See if to do
972    bool atRoot = model_->getNodeCount() == 0;
973    int passNumber = model_->getCurrentPassNumber();
974    if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber != 1)
975        return 0; // switched off
976    // Don't do if it was this heuristic which found solution!
977    if (this == model_->lastHeuristic())
978        return 0;
979    numRuns_++;
980    double cutoff;
981    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
982    double direction = model_->solver()->getObjSense();
983    cutoff *= direction;
984    cutoff = CoinMin(cutoff, solutionValue);
985    OsiSolverInterface * solver = model_->continuousSolver();
986    if (!solver)
987        solver = model_->solver();
988    const double * colLower = solver->getColLower();
989    const double * colUpper = solver->getColUpper();
990    const double * objective = solver->getObjCoefficients();
991
992    int numberColumns = model_->getNumCols();
993    int numberIntegers = model_->numberIntegers();
994    const int * integerVariable = model_->integerVariable();
995
996    int i;
997    bool solutionFound = false;
998    CoinWarmStartBasis saveBasis;
999    CoinWarmStartBasis * basis =
1000        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1001    if (basis) {
1002        saveBasis = * basis;
1003        delete basis;
1004    }
1005    // First just fix all integers as close to zero as possible
1006    OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
1007    for (i = 0; i < numberIntegers; i++) {
1008        int iColumn = integerVariable[i];
1009        double lower = colLower[iColumn];
1010        double upper = colUpper[iColumn];
1011        double value;
1012        if (lower > 0.0)
1013            value = lower;
1014        else if (upper < 0.0)
1015            value = upper;
1016        else
1017            value = 0.0;
1018        newSolver->setColLower(iColumn, value);
1019        newSolver->setColUpper(iColumn, value);
1020    }
1021    newSolver->initialSolve();
1022    if (newSolver->isProvenOptimal()) {
1023        double solValue = newSolver->getObjValue() * direction ;
1024        if (solValue < cutoff) {
1025            // we have a solution
1026            solutionFound = true;
1027            solutionValue = solValue;
1028            memcpy(betterSolution, newSolver->getColSolution(),
1029                   numberColumns*sizeof(double));
1030            printf("Naive fixing close to zero gave solution of %g\n", solutionValue);
1031            cutoff = solValue - model_->getCutoffIncrement();
1032        }
1033    }
1034    // Now fix all integers as close to zero if zero or large cost
1035    int nFix = 0;
1036    for (i = 0; i < numberIntegers; i++) {
1037        int iColumn = integerVariable[i];
1038        double lower = colLower[iColumn];
1039        double upper = colUpper[iColumn];
1040        double value;
1041        if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) {
1042            nFix++;
1043            if (lower > 0.0)
1044                value = lower;
1045            else if (upper < 0.0)
1046                value = upper;
1047            else
1048                value = 0.0;
1049            newSolver->setColLower(iColumn, value);
1050            newSolver->setColUpper(iColumn, value);
1051        } else {
1052            // set back to original
1053            newSolver->setColLower(iColumn, lower);
1054            newSolver->setColUpper(iColumn, upper);
1055        }
1056    }
1057    const double * solution = solver->getColSolution();
1058    if (nFix) {
1059        newSolver->setWarmStart(&saveBasis);
1060        newSolver->setColSolution(solution);
1061        newSolver->initialSolve();
1062        if (newSolver->isProvenOptimal()) {
1063            double solValue = newSolver->getObjValue() * direction ;
1064            if (solValue < cutoff) {
1065                // try branch and bound
1066                double * newSolution = new double [numberColumns];
1067                printf("%d fixed after fixing costs\n", nFix);
1068                int returnCode = smallBranchAndBound(newSolver,
1069                                                     numberNodes_, newSolution,
1070                                                     solutionValue,
1071                                                     solutionValue, "CbcHeuristicNaive1");
1072                if (returnCode < 0)
1073                    returnCode = 0; // returned on size
1074                if ((returnCode&2) != 0) {
1075                    // could add cut
1076                    returnCode &= ~2;
1077                }
1078                if (returnCode == 1) {
1079                    // solution
1080                    solutionFound = true;
1081                    memcpy(betterSolution, newSolution,
1082                           numberColumns*sizeof(double));
1083                    printf("Naive fixing zeros gave solution of %g\n", solutionValue);
1084                    cutoff = solutionValue - model_->getCutoffIncrement();
1085                }
1086                delete [] newSolution;
1087            }
1088        }
1089    }
1090#if 1
1091    newSolver->setObjSense(-direction); // maximize
1092    newSolver->setWarmStart(&saveBasis);
1093    newSolver->setColSolution(solution);
1094    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1095        double value = solution[iColumn];
1096        double lower = colLower[iColumn];
1097        double upper = colUpper[iColumn];
1098        double newLower;
1099        double newUpper;
1100        if (newSolver->isInteger(iColumn)) {
1101            newLower = CoinMax(lower, floor(value) - 2.0);
1102            newUpper = CoinMin(upper, ceil(value) + 2.0);
1103        } else {
1104            newLower = CoinMax(lower, value - 1.0e5);
1105            newUpper = CoinMin(upper, value + 1.0e-5);
1106        }
1107        newSolver->setColLower(iColumn, newLower);
1108        newSolver->setColUpper(iColumn, newUpper);
1109    }
1110    newSolver->initialSolve();
1111    if (newSolver->isProvenOptimal()) {
1112        double solValue = newSolver->getObjValue() * direction ;
1113        if (solValue < cutoff) {
1114            nFix = 0;
1115            newSolver->setObjSense(direction); // correct direction
1116            //const double * thisSolution = newSolver->getColSolution();
1117            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1118                double value = solution[iColumn];
1119                double lower = colLower[iColumn];
1120                double upper = colUpper[iColumn];
1121                double newLower = lower;
1122                double newUpper = upper;
1123                if (newSolver->isInteger(iColumn)) {
1124                    if (value < lower + 1.0e-6) {
1125                        nFix++;
1126                        newUpper = lower;
1127                    } else if (value > upper - 1.0e-6) {
1128                        nFix++;
1129                        newLower = upper;
1130                    } else {
1131                        newLower = CoinMax(lower, floor(value) - 2.0);
1132                        newUpper = CoinMin(upper, ceil(value) + 2.0);
1133                    }
1134                }
1135                newSolver->setColLower(iColumn, newLower);
1136                newSolver->setColUpper(iColumn, newUpper);
1137            }
1138            // try branch and bound
1139            double * newSolution = new double [numberColumns];
1140            printf("%d fixed after maximizing\n", nFix);
1141            int returnCode = smallBranchAndBound(newSolver,
1142                                                 numberNodes_, newSolution,
1143                                                 solutionValue,
1144                                                 solutionValue, "CbcHeuristicNaive1");
1145            if (returnCode < 0)
1146                returnCode = 0; // returned on size
1147            if ((returnCode&2) != 0) {
1148                // could add cut
1149                returnCode &= ~2;
1150            }
1151            if (returnCode == 1) {
1152                // solution
1153                solutionFound = true;
1154                memcpy(betterSolution, newSolution,
1155                       numberColumns*sizeof(double));
1156                printf("Naive maximizing gave solution of %g\n", solutionValue);
1157                cutoff = solutionValue - model_->getCutoffIncrement();
1158            }
1159            delete [] newSolution;
1160        }
1161    }
1162#endif
1163    delete newSolver;
1164    return solutionFound ? 1 : 0;
1165}
1166// update model
1167void CbcHeuristicNaive::setModel(CbcModel * model)
1168{
1169    model_ = model;
1170}
1171// Default Constructor
1172CbcHeuristicCrossover::CbcHeuristicCrossover()
1173        : CbcHeuristic(),
1174        numberSolutions_(0),
1175        useNumber_(3)
1176{
1177    setWhen(1);
1178}
1179
1180// Constructor with model - assumed before cuts
1181
1182CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model)
1183        : CbcHeuristic(model),
1184        numberSolutions_(0),
1185        useNumber_(3)
1186{
1187    setWhen(1);
1188    for (int i = 0; i < 10; i++)
1189        random_[i] = model.randomNumberGenerator()->randomDouble();
1190}
1191
1192// Destructor
1193CbcHeuristicCrossover::~CbcHeuristicCrossover ()
1194{
1195}
1196
1197// Clone
1198CbcHeuristic *
1199CbcHeuristicCrossover::clone() const
1200{
1201    return new CbcHeuristicCrossover(*this);
1202}
1203// Create C++ lines to get to current state
1204void
1205CbcHeuristicCrossover::generateCpp( FILE * fp)
1206{
1207    CbcHeuristicCrossover other;
1208    fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
1209    fprintf(fp, "3  CbcHeuristicCrossover crossover(*cbcModel);\n");
1210    CbcHeuristic::generateCpp(fp, "crossover");
1211    if (useNumber_ != other.useNumber_)
1212        fprintf(fp, "3  crossover.setNumberSolutions(%d);\n", useNumber_);
1213    else
1214        fprintf(fp, "4  crossover.setNumberSolutions(%d);\n", useNumber_);
1215    fprintf(fp, "3  cbcModel->addHeuristic(&crossover);\n");
1216}
1217
1218// Copy constructor
1219CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs)
1220        :
1221        CbcHeuristic(rhs),
1222        attempts_(rhs.attempts_),
1223        numberSolutions_(rhs.numberSolutions_),
1224        useNumber_(rhs.useNumber_)
1225{
1226    memcpy(random_, rhs.random_, 10*sizeof(double));
1227}
1228
1229// Assignment operator
1230CbcHeuristicCrossover &
1231CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs)
1232{
1233    if (this != &rhs) {
1234        CbcHeuristic::operator=(rhs);
1235        useNumber_ = rhs.useNumber_;
1236        attempts_ = rhs.attempts_;
1237        numberSolutions_ = rhs.numberSolutions_;
1238        memcpy(random_, rhs.random_, 10*sizeof(double));
1239    }
1240    return *this;
1241}
1242
1243// Resets stuff if model changes
1244void
1245CbcHeuristicCrossover::resetModel(CbcModel * model)
1246{
1247    CbcHeuristic::resetModel(model);
1248}
1249int
1250CbcHeuristicCrossover::solution(double & solutionValue,
1251                                double * betterSolution)
1252{
1253    if (when_ == 0)
1254        return 0;
1255    numCouldRun_++;
1256    bool useBest = (numberSolutions_ != model_->getSolutionCount());
1257    if (!useBest && (when_ % 10) == 1)
1258        return 0;
1259    numberSolutions_ = model_->getSolutionCount();
1260    OsiSolverInterface * continuousSolver = model_->continuousSolver();
1261    int useNumber = CoinMin(model_->numberSavedSolutions(), useNumber_);
1262    if (useNumber < 2 || !continuousSolver)
1263        return 0;
1264    // Fix later
1265    if (!useBest)
1266        abort();
1267    numRuns_++;
1268    double cutoff;
1269    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1270    double direction = model_->solver()->getObjSense();
1271    cutoff *= direction;
1272    cutoff = CoinMin(cutoff, solutionValue);
1273    OsiSolverInterface * solver = cloneBut(2);
1274    // But reset bounds
1275    solver->setColLower(continuousSolver->getColLower());
1276    solver->setColUpper(continuousSolver->getColUpper());
1277    int numberColumns = solver->getNumCols();
1278    // Fixed
1279    double * fixed = new double [numberColumns];
1280    for (int i = 0; i < numberColumns; i++)
1281        fixed[i] = -COIN_DBL_MAX;
1282    int whichSolution[10];
1283    for (int i = 0; i < useNumber; i++)
1284        whichSolution[i] = i;
1285    for (int i = 0; i < useNumber; i++) {
1286        int k = whichSolution[i];
1287        const double * solution = model_->savedSolution(k);
1288        for (int j = 0; j < numberColumns; j++) {
1289            if (solver->isInteger(j)) {
1290                if (fixed[j] == -COIN_DBL_MAX)
1291                    fixed[j] = floor(solution[j] + 0.5);
1292                else if (fabs(fixed[j] - solution[j]) > 1.0e-7)
1293                    fixed[j] = COIN_DBL_MAX;
1294            }
1295        }
1296    }
1297    const double * colLower = solver->getColLower();
1298    for (int i = 0; i < numberColumns; i++) {
1299        if (solver->isInteger(i)) {
1300            double value = fixed[i];
1301            if (value != COIN_DBL_MAX) {
1302                if (when_ < 10) {
1303                    solver->setColLower(i, value);
1304                    solver->setColUpper(i, value);
1305                } else if (value == colLower[i]) {
1306                    solver->setColUpper(i, value);
1307                }
1308            }
1309        }
1310    }
1311    int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution,
1312                                         solutionValue,
1313                                         solutionValue, "CbcHeuristicCrossover");
1314    if (returnCode < 0)
1315        returnCode = 0; // returned on size
1316    if ((returnCode&2) != 0) {
1317        // could add cut
1318        returnCode &= ~2;
1319    }
1320
1321    delete solver;
1322    return returnCode;
1323}
1324// update model
1325void CbcHeuristicCrossover::setModel(CbcModel * model)
1326{
1327    model_ = model;
1328    if (model) {
1329        for (int i = 0; i < 10; i++)
1330            random_[i] = model->randomNumberGenerator()->randomDouble();
1331    }
1332}
1333
1334
Note: See TracBrowser for help on using the repository browser.