source: trunk/Cbc/src/CbcHeuristicPivotAndFix.cpp @ 2094

Last change on this file since 2094 was 2094, checked in by forrest, 4 years ago

for memory leaks and heuristics and some experimental stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.5 KB
Line 
1/* $Id: CbcHeuristicPivotAndFix.cpp 2094 2014-11-18 11:15:36Z forrest $ */
2// Copyright (C) 2008, 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#include <vector>
15
16#include "OsiSolverInterface.hpp"
17#include "CbcModel.hpp"
18#include "CbcMessage.hpp"
19#include "CbcHeuristicPivotAndFix.hpp"
20#include "OsiClpSolverInterface.hpp"
21#include  "CoinTime.hpp"
22
23//#define FORNOW
24
25// Default Constructor
26CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix()
27        : CbcHeuristic()
28{
29}
30
31// Constructor with model - assumed before cuts
32
33CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix(CbcModel & model)
34        : CbcHeuristic(model)
35{
36}
37
38// Destructor
39CbcHeuristicPivotAndFix::~CbcHeuristicPivotAndFix ()
40{
41}
42
43// Clone
44CbcHeuristic *
45CbcHeuristicPivotAndFix::clone() const
46{
47    return new CbcHeuristicPivotAndFix(*this);
48}
49// Create C++ lines to get to current state
50void
51CbcHeuristicPivotAndFix::generateCpp( FILE * fp)
52{
53    CbcHeuristicPivotAndFix other;
54    fprintf(fp, "0#include \"CbcHeuristicPivotAndFix.hpp\"\n");
55    fprintf(fp, "3  CbcHeuristicPivotAndFix heuristicPFX(*cbcModel);\n");
56    CbcHeuristic::generateCpp(fp, "heuristicPFX");
57    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicPFX);\n");
58}
59
60// Copy constructor
61CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix(const CbcHeuristicPivotAndFix & rhs)
62        :
63        CbcHeuristic(rhs)
64{
65}
66
67// Assignment operator
68CbcHeuristicPivotAndFix &
69CbcHeuristicPivotAndFix::operator=( const CbcHeuristicPivotAndFix & rhs)
70{
71    if (this != &rhs) {
72        CbcHeuristic::operator=(rhs);
73    }
74    return *this;
75}
76
77// Resets stuff if model changes
78void
79CbcHeuristicPivotAndFix::resetModel(CbcModel * /*model*/)
80{
81    //CbcHeuristic::resetModel(model);
82}
83/*
84  Comments needed
85  Returns 1 if solution, 0 if not */
86int
87CbcHeuristicPivotAndFix::solution(double & /*solutionValue*/,
88                                  double * /*betterSolution*/)
89{
90
91    numCouldRun_++; // Todo: Ask JJHF what this for.
92    std::cout << "Entering Pivot-and-Fix Heuristic" << std::endl;
93
94#ifdef HEURISTIC_INFORM
95    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
96           heuristicName(),numRuns_,numCouldRun_,when_);
97#endif
98#ifdef FORNOW
99    std::cout << "Lucky you! You're in the Pivot-and-Fix Heuristic" << std::endl;
100    // The struct should be moved to member data
101
102
103
104    typedef struct {
105        int numberSolutions;
106        int maximumSolutions;
107        int numberColumns;
108        double ** solution;
109        int * numberUnsatisfied;
110    } clpSolution;
111
112    double start = CoinCpuTime();
113
114    OsiClpSolverInterface * clpSolverOriginal
115    = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
116    assert (clpSolverOriginal);
117    OsiClpSolverInterface *clpSolver(clpSolverOriginal);
118
119    ClpSimplex * simplex = clpSolver->getModelPtr();
120
121    // Initialize the structure holding the solutions
122    clpSolution solutions;
123    // Set typeStruct field of ClpTrustedData struct to one.
124    // This tells Clp it's "Mahdi!"
125    ClpTrustedData trustedSolutions;
126    trustedSolutions.typeStruct = 1;
127    trustedSolutions.data = &solutions;
128    solutions.numberSolutions = 0;
129    solutions.maximumSolutions = 0;
130    solutions.numberColumns = simplex->numberColumns();
131    solutions.solution = NULL;
132    solutions.numberUnsatisfied = NULL;
133    simplex->setTrustedUserPointer(&trustedSolutions);
134
135    // Solve from all slack to get some points
136    simplex->allSlackBasis();
137    simplex->primal();
138
139    // -------------------------------------------------
140    // Get the problem information
141    // - get the number of cols and rows
142    int numCols = clpSolver->getNumCols();
143    int numRows = clpSolver->getNumRows();
144
145    // - get the right hand side of the rows
146    const double * rhs = clpSolver->getRightHandSide();
147
148    // - find the integer variables
149    bool * varClassInt = new bool[numCols];
150    int numInt = 0;
151    for (int i = 0; i < numCols; i++) {
152        if (clpSolver->isContinuous(i))
153            varClassInt[i] = 0;
154        else {
155            varClassInt[i] = 1;
156            numInt++;
157        }
158    }
159
160    // -Get the rows sense
161    const char * rowSense;
162    rowSense = clpSolver->getRowSense();
163
164    // -Get the objective coefficients
165    const double *objCoefficients = clpSolver->getObjCoefficients();
166    double *originalObjCoeff = new double [numCols];
167    for (int i = 0; i < numCols; i++)
168        originalObjCoeff[i] = objCoefficients[i];
169
170    // -Get the matrix of the problem
171    double ** matrix = new double * [numRows];
172    for (int i = 0; i < numRows; i++) {
173        matrix[i] = new double[numCols];
174        for (int j = 0; j < numCols; j++)
175            matrix[i][j] = 0;
176    }
177    const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
178    const double * matrixElements = matrixByRow->getElements();
179    const int * matrixIndices = matrixByRow->getIndices();
180    const int * matrixStarts = matrixByRow->getVectorStarts();
181    for (int j = 0; j < numRows; j++) {
182        for (int i = matrixStarts[j]; i < matrixStarts[j+1]; i++) {
183            matrix[j][matrixIndices[i]] = matrixElements[i];
184        }
185    }
186
187    // The newObj is the randomly perturbed constraint used to find new
188    // corner points
189    double * newObj = new double [numCols];
190
191    // Set the random seed
192    srand ( time(NULL) + 1);
193    int randNum;
194
195    // We're going to add a new row to the LP formulation
196    // after finding each new solution.
197    // Adding a new row requires the new elements and the new indices.
198    // The elements are original objective function coefficients.
199    // The indicies are the (dense) columns indices stored in addRowIndex.
200    // The rhs is the value of the new solution stored in solutionValue.
201    int * addRowIndex = new int[numCols];
202    for (int i = 0; i < numCols; i++)
203        addRowIndex[i] = i;
204
205    // The number of feasible solutions found by the PF heuristic.
206    // This controls the return code of the solution() method.
207    int numFeasibles = 0;
208
209    // Shuffle the rows
210    int * index = new int [numRows];
211    for (int i = 0; i < numRows; i++)
212        index[i] = i;
213    for (int i = 0; i < numRows; i++) {
214        int temp = index[i];
215        int randNumTemp = i + (rand() % (numRows - i));
216        index[i] = index[randNumTemp];
217        index[randNumTemp] = temp;
218    }
219
220    // In the clpSolution struct, we store a lot of column solutions.
221    // For each perturb objective, we store the solution from each
222    // iteration of the LP solve.
223    // For each perturb objective, we look at the collection of
224    // solutions to do something extremly intelligent :-)
225    // We could (and should..and will :-) wipe out the block of
226    // solutions when we're done with them. But for now, we just move on
227    // and store the next block of solutions for the next (perturbed)
228    // objective.
229    // The variable startIndex tells us where the new block begins.
230    int startIndex = 0;
231
232    // At most "fixThreshold" number of integer variables can be unsatisfied
233    // for calling smallBranchAndBound().
234    // The PF Heuristic only fixes fixThreshold number of variables to
235    // their integer values. Not more. Not less. The reason is to give
236    // the smallBB some opportunity to find better solutions. If we fix
237    // everything it might be too many (leading the heuristic to come up
238    // with infeasibility rather than a useful result).
239    // (This is an important paramater. And it is dynamically set.)
240    double fixThreshold;
241    /*
242        if(numInt > 400)
243        fixThreshold = 17*sqrt(numInt);
244        if(numInt<=400 && numInt>100)
245        fixThreshold = 5*sqrt(numInt);
246        if(numInt<=100)
247        fixThreshold = 4*sqrt(numInt);
248    */
249    // Initialize fixThreshold based on the number of integer
250    // variables
251    if (numInt <= 100)
252        fixThreshold = .35 * numInt;
253    if (numInt > 100 && numInt < 1000)
254        fixThreshold = .85 * numInt;
255    if (numInt >= 1000)
256        fixThreshold = .1 * numInt;
257
258    // Whenever the dynamic system for changing fixThreshold
259    // kicks in, it changes the parameter by the
260    // fixThresholdChange amount.
261    // (The 25% should be member data and tuned. Another paper!)
262    double fixThresholdChange = 0.25 * fixThreshold;
263
264    // maxNode is the maximum number of nodes we allow smallBB to
265    // search. It's initialized to 400 and changed dynamically.
266    // The 400 should be member data, if we become virtuous.
267    int maxNode = 400;
268
269    // We control the decision to change maxNode through the boolean
270    // variable  changeMaxNode. The boolean variable is initialized to
271    // true and gets set to false under a condition (and is never true
272    // again.)
273    // It's flipped off and stays off (in the current incarnation of PF)
274    bool changeMaxNode = 1;
275
276    // The sumReturnCode is used for the dynamic system that sets
277    // fixThreshold and changeMaxNode.
278    //
279    // We track what's happening in sumReturnCode. There are 8 switches.
280    // The first 5 switches corresponds to a return code for smallBB.
281    //
282    // We want to know how many times we consecutively get the same
283    // return code.
284    //
285    // If "good" return codes are happening often enough, we're happy.
286    //
287    // If a "bad" returncodes happen consecutively, we want to
288    // change something.
289    //
290    // The switch 5 is the number of times PF didn't call smallBB
291    // becuase the number of integer variables that took integer values
292    // was less than fixThreshold.
293    //
294    // The swicth 6 was added for a brilliant idea...to be announced
295    // later (another paper!)
296    //
297    // The switch 7 is the one that changes the max node. Read the
298    // code. (Todo: Verbalize the brilliant idea for the masses.)
299    //
300    int sumReturnCode[8];
301    /*
302      sumReturnCode[0] ~ -1 --> problem too big for smallBB
303      sumReturnCode[1] ~ 0  --> smallBB not finshed and no soln
304      sumReturnCode[2] ~ 1  --> smallBB not finshed and there is a soln
305      sumReturnCode[3] ~ 2  --> smallBB finished and no soln
306      sumReturnCode[4] ~ 3  --> smallBB finished and there is a soln
307      sumReturnCode[5] ~ didn't call smallBranchAndBound too few to fix
308      sumReturnCode[6] ~ didn't call smallBranchAndBound too many unsatisfied
309      sumReturnCode[7] ~ the same as sumReturnCode[1] but becomes zero just if the returnCode is not 0
310    */
311
312    for (int i = 0; i < 8; i++)
313        sumReturnCode[i] = 0;
314    int * colIndex = new int[numCols];
315    for (int i = 0; i < numCols; i++)
316        colIndex[i] = i;
317    double cutoff = COIN_DBL_MAX;
318    bool didMiniBB;
319
320    // Main loop
321    for (int i = 0; i < numRows; i++) {
322        // track the number of mini-bb for the dynamic threshold setting
323        didMiniBB = 0;
324
325        for (int k = startIndex; k < solutions.numberSolutions; k++)
326            //if the point has 0 unsatisfied variables; make sure it is
327            //feasible. Check integer feasiblity and constraints.
328            if (solutions.numberUnsatisfied[k] == 0) {
329                double feasibility = 1;
330                //check integer feasibility
331                for (int icol = 0; icol < numCols; icol++) {
332                    double closest = floor(solutions.solution[k][icol] + 0.5);
333                    if (varClassInt[icol] && (fabs(solutions.solution[k][icol] - closest) > 1e-6)) {
334                        feasibility = 0;
335                        break;
336                    }
337                }
338                //check if the solution satisfies the constraints
339                for (int irow = 0; irow < numRows; irow++) {
340                    double lhs = 0;
341                    for (int j = 0; j < numCols; j++)
342                        lhs += matrix[irow][j] * solutions.solution[k][j];
343                    if (rowSense[irow] == 'L' && lhs > rhs[irow] + 1e-6) {
344                        feasibility = 0;
345                        break;
346                    }
347                    if (rowSense[irow] == 'G' && lhs < rhs[irow] - 1e-6) {
348                        feasibility = 0;
349                        break;
350                    }
351                    if (rowSense[irow] == 'E' && (lhs - rhs[irow] > 1e-6 || lhs - rhs[irow] < -1e-6)) {
352                        feasibility = 0;
353                        break;
354                    }
355                }
356
357                //if feasible, find the objective value and set the cutoff
358                // for the smallBB and add a new constraint to the LP
359                // (and update the best solution found so far for the
360                // return arguments)
361                if (feasibility) {
362                    double objectiveValue = 0;
363                    for (int j = 0; j < numCols; j++)
364                        objectiveValue += solutions.solution[k][j] * originalObjCoeff[j];
365                    cutoff = objectiveValue;
366                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);
367
368                    // Todo: pick up the best solution in the block (not
369                    // the last).
370                    solutionValue = objectiveValue;
371                    for (int m = 0; m < numCols; m++)
372                        betterSolution[m] = solutions.solution[k][m];
373                    numFeasibles++;
374                }
375            }
376
377        // Go through the block of solution and decide if to call smallBB
378        for (int k = startIndex; k < solutions.numberSolutions; k++) {
379            if (solutions.numberUnsatisfied[k] <= fixThreshold) {
380                // get new copy
381                OsiSolverInterface * newSolver;
382                newSolver = new OsiClpSolverInterface(*clpSolver);
383                newSolver->setObjSense(1);
384                newSolver->setObjective(originalObjCoeff);
385                int numberColumns = newSolver->getNumCols();
386                int numFixed = 0;
387
388                // Fix the first fixThreshold number of integer vars
389                // that are satisfied
390                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
391                    if (newSolver->isInteger(iColumn)) {
392                        double value = solutions.solution[k][iColumn];
393                        double intValue = floor(value + 0.5);
394                        if (fabs(value - intValue) < 1.0e-5) {
395                            newSolver->setColLower(iColumn, intValue);
396                            newSolver->setColUpper(iColumn, intValue);
397                            numFixed++;
398                            if (numFixed > numInt - fixThreshold)
399                                break;
400                        }
401                    }
402                }
403                COIN_DETAIL_PRINT(printf("numFixed: %d\n", numFixed));
404                COIN_DETAIL_PRINT(printf("fixThreshold: %f\n", fixThreshold));
405                COIN_DETAIL_PRINT(printf("numInt: %d\n", numInt));
406                double *newSolution = new double[numCols];
407                double newSolutionValue;
408
409                // Call smallBB on the modified problem
410                int returnCode = smallBranchAndBound(newSolver, maxNode, newSolution,
411                                                     newSolutionValue, cutoff, "mini");
412
413                // If smallBB found a solution, update the better
414                // solution and solutionValue (we gave smallBB our
415                // cutoff, so it only finds improving solutions)
416                if (returnCode == 1 || returnCode == 3) {
417                    numFeasibles ++;
418                    solutionValue = newSolutionValue;
419                    for (int m = 0; m < numCols; m++)
420                        betterSolution[m] = newSolution[m];
421                    COIN_DETAIL_PRINT(printf("cutoff: %f\n", newSolutionValue));
422                    COIN_DETAIL_PRINT(printf("time: %.2lf\n", CoinCpuTime() - start));
423                }
424                didMiniBB = 1;
425                COIN_DETAIL_PRINT(printf("returnCode: %d\n", returnCode));
426
427                //Update sumReturnCode array
428                for (int iRC = 0; iRC < 6; iRC++) {
429                    if (iRC == returnCode + 1)
430                        sumReturnCode[iRC]++;
431                    else
432                        sumReturnCode[iRC] = 0;
433                }
434                if (returnCode != 0)
435                    sumReturnCode[7] = 0;
436                else
437                    sumReturnCode[7]++;
438                if (returnCode == 1 || returnCode == 3) {
439                    cutoff = newSolutionValue;
440                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);
441                    COIN_DETAIL_PRINT(printf("******************\n\n*****************\n"));
442                }
443                break;
444            }
445        }
446
447        if (!didMiniBB && solutions.numberSolutions - startIndex > 0) {
448            sumReturnCode[5]++;
449            for (int iRC = 0; iRC < 5; iRC++)
450                sumReturnCode[iRC] = 0;
451        }
452
453        //Change "fixThreshold" if needed
454        // using the data we've recorded in sumReturnCode
455        if (sumReturnCode[1] >= 3)
456            fixThreshold -= fixThresholdChange;
457        if (sumReturnCode[7] >= 3 && changeMaxNode) {
458            maxNode *= 5;
459            changeMaxNode = 0;
460        }
461        if (sumReturnCode[3] >= 3 && fixThreshold < 0.95 * numInt)
462            fixThreshold += fixThresholdChange;
463        if (sumReturnCode[5] >= 4)
464            fixThreshold += fixThresholdChange;
465        if (sumReturnCode[0] > 3)
466            fixThreshold -= fixThresholdChange;
467
468        startIndex = solutions.numberSolutions;
469
470        //Check if the maximum iterations limit is reached
471        // rlh: Ask John how this is working with the change to trustedUserPtr.
472        if (solutions.numberSolutions > 20000)
473            break;
474
475        // The first time in this loop PF solves orig LP.
476
477        //Generate the random objective function
478        randNum = rand() % 10 + 1;
479        randNum = fmod(randNum, 2);
480        for (int j = 0; j < numCols; j++) {
481            if (randNum == 1)
482                if (fabs(matrix[index[i]][j]) < 1e-6)
483                    newObj[j] = 0.1;
484                else
485                    newObj[j] = matrix[index[i]][j] * 1.1;
486            else if (fabs(matrix[index[i]][j]) < 1e-6)
487                newObj[j] = -0.1;
488            else
489                newObj[j] = matrix[index[i]][j] * 0.9;
490        }
491        clpSolver->setObjective(newObj);
492        if (rowSense[i] == 'L')
493            clpSolver->setObjSense(-1);
494        else
495            // Todo #1: We don't need to solve the LPs to optimality.
496            // We just need corner points.
497            // There's a problem in stopping Clp that needs to be looked
498            // into. So for now, we solve optimality.
499            clpSolver->setObjSense(1);
500        //        simplex->setMaximumIterations(100);
501        clpSolver->getModelPtr()->primal(1);
502        //        simplex->setMaximumIterations(100000);
503#ifdef COIN_DETAIL
504        printf("cutoff: %f\n", cutoff);
505        printf("time: %.2f\n", CoinCpuTime() - start);
506        for (int iRC = 0; iRC < 8; iRC++)
507            printf("%d ", sumReturnCode[iRC]);
508        printf("\nfixThreshold: %f\n", fixThreshold);
509        printf("numInt: %d\n", numInt);
510        printf("\n---------------------------------------------------------------- %d\n", i);
511#endif
512
513        //temp:
514        if (i > 3) break;
515
516    }
517
518    COIN_DETAIL_PRINT(printf("Best Feasible Found: %f\n", cutoff));
519    COIN_DETAIL_PRINT(printf("Total time: %.2f\n", CoinCpuTime() - start));
520
521    if (numFeasibles == 0) {
522        return 0;
523    }
524
525
526
527    // We found something better
528    std::cout << "See you soon! You're leaving the Pivot-and-Fix Heuristic" << std::endl;
529    std::cout << std::endl;
530
531    return 1;
532#endif
533
534    return 0;
535
536}
537
538
539
540
541// update model
542void CbcHeuristicPivotAndFix::setModel(CbcModel * )
543{
544    // probably same as resetModel
545}
546
547
Note: See TracBrowser for help on using the repository browser.