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

Last change on this file since 1424 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

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