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

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

Added comments from Lou

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