source: stable/2.7/Cbc/src/CbcHeuristicPivotAndFix.cpp @ 1577

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

Change to EPL license notice.

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