source: trunk/Cbc/src/CbcHeuristicRandRound.cpp @ 1514

Last change on this file since 1514 was 1337, checked in by bjarni, 10 years ago

Add static_cast<unsigned int> for time(NULL) to avoid MSVC warnings when used as argument for srand()

File size: 17.1 KB
Line 
1/* $Id: CbcHeuristicRandRound.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 "CoinHelperFunctions.hpp"
15#include "OsiSolverInterface.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcHeuristicRandRound.hpp"
19#include "OsiClpSolverInterface.hpp"
20#include  "CoinTime.hpp"
21
22static inline int intRand(const int range)
23{
24    return static_cast<int> (floor(CoinDrand48() * range));
25}
26
27// Default Constructor
28CbcHeuristicRandRound::CbcHeuristicRandRound()
29        : CbcHeuristic()
30{
31}
32
33// Constructor with model - assumed before cuts
34CbcHeuristicRandRound::CbcHeuristicRandRound(CbcModel & model)
35        : CbcHeuristic(model)
36{
37    model_ = &model;
38    setWhen(1);
39}
40
41// Destructor
42CbcHeuristicRandRound::~CbcHeuristicRandRound ()
43{
44}
45
46// Clone
47CbcHeuristic *
48CbcHeuristicRandRound::clone() const
49{
50    return new CbcHeuristicRandRound(*this);
51}
52
53// Create C++ lines to get to current state
54void
55CbcHeuristicRandRound::generateCpp( FILE * fp)
56{
57    CbcHeuristicRandRound other;
58    fprintf(fp, "0#include \"CbcHeuristicRandRound.hpp\"\n");
59    fprintf(fp, "3  CbcHeuristicRandRound heuristicPFX(*cbcModel);\n");
60    CbcHeuristic::generateCpp(fp, "heuristicPFX");
61    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicPFX);\n");
62}
63
64// Copy constructor
65CbcHeuristicRandRound::CbcHeuristicRandRound(const CbcHeuristicRandRound & rhs)
66        :
67        CbcHeuristic(rhs)
68{
69}
70
71// Assignment operator
72CbcHeuristicRandRound &
73CbcHeuristicRandRound::operator=( const CbcHeuristicRandRound & rhs)
74{
75    if (this != &rhs) {
76        CbcHeuristic::operator=(rhs);
77    }
78    return *this;
79}
80
81// Resets stuff if model changes
82void
83CbcHeuristicRandRound::resetModel(CbcModel * model)
84{
85    CbcHeuristic::resetModel(model);
86}
87
88/*
89  Randomized Rounding Heuristic
90  Returns 1 if solution, 0 if not
91*/
92int
93CbcHeuristicRandRound::solution(double & solutionValue,
94                                double * betterSolution)
95{
96    // rlh: Todo: Memory Cleanup
97
98    //  std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
99
100    setWhen(1);  // setWhen(1) didn't have the effect I expected (e.g., run once).
101
102    // Run only once.
103    //
104    //    See if at root node
105    bool atRoot = model_->getNodeCount() == 0;
106    int passNumber = model_->getCurrentPassNumber();
107    //    Just do once
108    if (!atRoot || passNumber != 1) {
109        // std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
110        return 0;
111    }
112
113    std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
114    typedef struct {
115        int numberSolutions;
116        int maximumSolutions;
117        int numberColumns;
118        double ** solution;
119        int * numberUnsatisfied;
120    } clpSolution;
121
122    double start = CoinCpuTime();
123    numCouldRun_++; //
124    // Todo: Ask JJHF what "number of times
125    // the heuristic could run" means.
126
127    OsiSolverInterface * solver = model_->solver()->clone();
128    double primalTolerance ;
129    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
130    OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *> (solver);
131    assert (clpSolver);
132    ClpSimplex * simplex = clpSolver->getModelPtr();
133
134    // Initialize the structure holding the solutions for the Simplex iterations
135    clpSolution solutions;
136    // Set typeStruct field of ClpTrustedData struct to 1 to indicate
137    // desired behavior for  RandRound heuristic (which is what?)
138    ClpTrustedData trustedSolutions;
139    trustedSolutions.typeStruct = 1;
140    trustedSolutions.data = &solutions;
141    solutions.numberSolutions = 0;
142    solutions.maximumSolutions = 0;
143    solutions.numberColumns = simplex->numberColumns();
144    solutions.solution = NULL;
145    solutions.numberUnsatisfied = NULL;
146    simplex->setTrustedUserPointer(&trustedSolutions);
147
148    // Solve from all slack to get some points
149    simplex->allSlackBasis();
150
151    // Calling primal() invalidates pointers to some rim vectors,
152    // like...row sense (!)
153    simplex->primal();
154
155    // 1. Okay - so a workaround would be to copy the data I want BEFORE
156    // calling primal.
157    // 2. Another approach is to ask the simplex solvers NOT to mess up my
158    // rims.
159    // 3. See freeCachedResults() for what is getting
160    // deleted. Everything else points into the structure.
161    // ...or use collower and colupper rather than rowsense.
162    // ..store address of where one of these
163
164    // Store the basic problem information
165    // -Get the number of columns, rows and rhs vector
166    int numCols = clpSolver->getNumCols();
167    int numRows = clpSolver->getNumRows();
168
169    // Find the integer variables (use columnType(?))
170    // One if not continuous, that is binary or general integer)
171    // columnType() = 0 continuous
172    //              = 1 binary
173    //              = 2 general integer
174    bool * varClassInt = new bool[numCols];
175    const char* columnType = clpSolver->columnType();
176    int numGenInt = 0;
177    for (int i = 0; i < numCols; i++) {
178        if (clpSolver->isContinuous(i))
179            varClassInt[i] = 0;
180        else
181            varClassInt[i] = 1;
182        if (columnType[i] == 2) numGenInt++;
183    }
184
185    // Heuristic is for problems with general integer variables.
186    // If there are none, quit.
187    if (numGenInt++ < 1) {
188        delete [] varClassInt ;
189        std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
190        return 0;
191    }
192
193
194    // -Get the rows sense
195    const char * rowSense;
196    rowSense = clpSolver->getRowSense();
197
198    // -Get the objective coefficients
199    double *originalObjCoeff = CoinCopyOfArray(clpSolver->getObjCoefficients(), numCols);
200
201    // -Get the matrix of the problem
202    // rlh: look at using sparse representation
203    double ** matrix = new double * [numRows];
204    for (int i = 0; i < numRows; i++) {
205        matrix[i] = new double[numCols];
206        for (int j = 0; j < numCols; j++)
207            matrix[i][j] = 0;
208    }
209
210    const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
211    const double * matrixElements = matrixByRow->getElements();
212    const int * matrixIndices = matrixByRow->getIndices();
213    const int * matrixStarts = matrixByRow->getVectorStarts();
214    for (int j = 0; j < numRows; j++) {
215        for (int i = matrixStarts[j]; i < matrixStarts[j+1]; i++) {
216            matrix[j][matrixIndices[i]] = matrixElements[i];
217        }
218    }
219
220    double * newObj = new double [numCols];
221    srand ( static_cast<unsigned int>(time(NULL) + 1));
222    int randNum;
223
224    // Shuffle the rows:
225    // Put the rows in a random order
226    // so that the optimal solution is a different corner point than the
227    // starting point.
228    int * index = new int [numRows];
229    for (int i = 0; i < numRows; i++)
230        index[i] = i;
231    for (int i = 0; i < numRows; i++) {
232        int temp = index[i];
233        int randNumTemp = i + intRand(numRows - i);
234        index[i] = index[randNumTemp];
235        index[randNumTemp] = temp;
236    }
237
238    // Start finding corner points by iteratively doing the following:
239    // - contruct a randomly tilted objective
240    // - solve
241    for (int i = 0; i < numRows; i++) {
242        // TODO: that 10,000 could be a param in the member data
243        if (solutions.numberSolutions  > 10000)
244            break;
245        randNum = intRand(2);
246        for (int j = 0; j < numCols; j++) {
247            // for row i and column j vary the coefficient "a bit"
248            if (randNum == 1)
249                // if the element is zero, then set the new obj
250                // coefficient to 0.1 (i.e., round up)
251                if (fabs(matrix[index[i]][j]) < primalTolerance)
252                    newObj[j] = 0.1;
253                else
254                    // if the element is nonzero, then increase the new obj
255                    // coefficient "a bit"
256                    newObj[j] = matrix[index[i]][j] * 1.1;
257            else
258                // if randnum is 2, then
259                // if the element is zero, then set the new obj coeffient
260                // to NEGATIVE 0.1 (i.e., round down)
261                if (fabs(matrix[index[i]][j]) < primalTolerance)
262                    newObj[j] = -0.1;
263                else
264                    // if the element is nonzero, then DEcrease the new obj coeffienct "a bit"
265                    newObj[j] = matrix[index[i]][j] * 0.9;
266        }
267        // Use the new "tilted" objective
268        clpSolver->setObjective(newObj);
269
270        // Based on the row sense, we decide whether to max or min
271        if (rowSense[i] == 'L')
272            clpSolver->setObjSense(-1);
273        else
274            clpSolver->setObjSense(1);
275
276        // Solve with primal simplex
277        simplex->primal(1);
278        // rlh+ll: This was the original code. But we already have the
279        // model pointer (it's in simplex). And, calling getModelPtr()
280        // invalidates the cached data in the OsiClpSolverInterface
281        // object, which means our precious rowsens is lost. So let's
282        // not use the line below...
283        /******* clpSolver->getModelPtr()->primal(1); */
284        printf("---------------------------------------------------------------- %d\n", i);
285    }
286    // Iteratively do this process until...
287    // either you reach the max number of corner points (aka 10K)
288    // or all the rows have been used as an objective.
289
290    // Look at solutions
291    int numberSolutions = solutions.numberSolutions;
292    //const char * integerInfo = simplex->integerInformation();
293    //const double * columnLower = simplex->columnLower();
294    //const double * columnUpper = simplex->columnUpper();
295    printf("there are %d solutions\n", numberSolutions);
296
297    // Up to here we have all the corner points
298    // Now we need to do the random walks and roundings
299
300    double ** cornerPoints = new double * [numberSolutions];
301    for (int j = 0; j < numberSolutions; j++)
302        cornerPoints[j] = solutions.solution[j];
303
304    bool feasibility = 1;
305    // rlh: use some COIN max instead of 1e30 (?)
306    double bestObj = 1e30;
307    std::vector< std::vector <double> > feasibles;
308    int numFeasibles = 0;
309
310    // Check the feasibility of the corner points
311    int numCornerPoints = numberSolutions;
312
313    const double * rhs = clpSolver->getRightHandSide();
314    // rlh: row sense hasn't changed. why a fresh copy?
315    // Delete next line.
316    rowSense = clpSolver->getRowSense();
317
318    for (int i = 0; i < numCornerPoints; i++) {
319        //get the objective value for this this point
320        double objValue = 0;
321        for (int k = 0; k < numCols; k++)
322            objValue += cornerPoints[i][k] * originalObjCoeff[k];
323
324        if (objValue < bestObj) {
325            // check integer feasibility
326            feasibility = 1;
327            for (int j = 0; j < numCols; j++) {
328                if (varClassInt[j]) {
329                    double closest = floor(cornerPoints[i][j] + 0.5);
330                    if (fabs(cornerPoints[i][j] - closest) > primalTolerance) {
331                        feasibility = 0;
332                        break;
333                    }
334                }
335            }
336            // check all constraints satisfied
337            if (feasibility) {
338                for (int irow = 0; irow < numRows; irow++) {
339                    double lhs = 0;
340                    for (int j = 0; j < numCols; j++) {
341                        lhs += matrix[irow][j] * cornerPoints[i][j];
342                    }
343                    if (rowSense[irow] == 'L' && lhs > rhs[irow] + primalTolerance) {
344                        feasibility = 0;
345                        break;
346                    }
347                    if (rowSense[irow] == 'G' && lhs < rhs[irow] - primalTolerance) {
348                        feasibility = 0;
349                        break;
350                    }
351                    if (rowSense[irow] == 'E' && (lhs - rhs[irow] > primalTolerance || lhs - rhs[irow] < -primalTolerance)) {
352                        feasibility = 0;
353                        break;
354                    }
355                }
356            }
357
358            if (feasibility) {
359                numFeasibles++;
360                feasibles.push_back(std::vector <double> (numCols));
361                for (int k = 0; k < numCols; k++)
362                    feasibles[numFeasibles-1][k] = cornerPoints[i][k];
363                printf("obj: %f\n", objValue);
364                if (objValue < bestObj)
365                    bestObj = objValue;
366            }
367        }
368    }
369    int numFeasibleCorners;
370    numFeasibleCorners = numFeasibles;
371    //find the center of gravity of the corner points as the first random point
372    double * rp = new double[numCols];
373    for (int i = 0; i < numCols; i++) {
374        rp[i] = 0;
375        for (int j = 0; j < numCornerPoints; j++) {
376            rp[i] += cornerPoints[j][i];
377        }
378        rp[i] = rp[i] / numCornerPoints;
379    }
380
381    //-------------------------------------------
382    //main loop:
383    // -generate the next random point
384    // -round the random point
385    // -check the feasibility of the random point
386    //-------------------------------------------
387
388    srand ( static_cast<unsigned int>(time(NULL) + 1));
389    int numRandomPoints = 0;
390    while (numRandomPoints < 50000) {
391        numRandomPoints++;
392        //generate the next random point
393        int randomIndex = intRand(numCornerPoints);
394        double random = CoinDrand48();
395        for (int i = 0; i < numCols; i++) {
396            rp[i] = (random * (cornerPoints[randomIndex][i] - rp[i])) + rp[i];
397        }
398
399        //CRISP ROUNDING
400        //round the random point just generated
401        double * roundRp = new double[numCols];
402        for (int i = 0; i < numCols; i++) {
403            roundRp[i] = rp[i];
404            if (varClassInt[i]) {
405                if (rp[i] >= 0) {
406                    if (fmod(rp[i], 1) > 0.5)
407                        roundRp[i] = floor(rp[i]) + 1;
408                    else
409                        roundRp[i] = floor(rp[i]);
410                } else {
411                    if (fabs(fmod(rp[i], 1)) > 0.5)
412                        roundRp[i] = floor(rp[i]);
413                    else
414                        roundRp[i] = floor(rp[i]) + 1;
415
416                }
417            }
418        }
419
420
421        //SOFT ROUNDING
422        // Look at original files for the "how to" on soft rounding;
423        // Soft rounding omitted here.
424
425        //Check the feasibility of the rounded random point
426        // -Check the feasibility
427        // -Get the rows sense
428        rowSense = clpSolver->getRowSense();
429        rhs = clpSolver->getRightHandSide();
430
431        //get the objective value for this feasible point
432        double objValue = 0;
433        for (int i = 0; i < numCols; i++)
434            objValue += roundRp[i] * originalObjCoeff[i];
435
436        if (objValue < bestObj) {
437            feasibility = 1;
438            for (int i = 0; i < numRows; i++) {
439                double lhs = 0;
440                for (int j = 0; j < numCols; j++) {
441                    lhs += matrix[i][j] * roundRp[j];
442                }
443                if (rowSense[i] == 'L' && lhs > rhs[i] + primalTolerance) {
444                    feasibility = 0;
445                    break;
446                }
447                if (rowSense[i] == 'G' && lhs < rhs[i] - primalTolerance) {
448                    feasibility = 0;
449                    break;
450                }
451                if (rowSense[i] == 'E' && (lhs - rhs[i] > primalTolerance || lhs - rhs[i] < -primalTolerance)) {
452                    feasibility = 0;
453                    break;
454                }
455            }
456            if (feasibility) {
457                printf("Feasible Found.\n");
458                printf("%.2f\n", CoinCpuTime() - start);
459                numFeasibles++;
460                feasibles.push_back(std::vector <double> (numCols));
461                for (int i = 0; i < numCols; i++)
462                    feasibles[numFeasibles-1][i] = roundRp[i];
463                printf("obj: %f\n", objValue);
464                if (objValue < bestObj)
465                    bestObj = objValue;
466            }
467        }
468        delete [] roundRp;
469    }
470    printf("Number of Feasible Corners: %d\n", numFeasibleCorners);
471    printf("Number of Feasibles Found: %d\n", numFeasibles);
472    if (numFeasibles > 0)
473        printf("Best Objective: %f\n", bestObj);
474    printf("time: %.2f\n", CoinCpuTime() - start);
475
476    if (numFeasibles == 0) {
477        // cleanup
478        delete [] varClassInt;
479        for (int i = 0; i < numRows; i++)
480            delete matrix[i];
481        delete [] matrix;
482        delete [] newObj;
483        delete [] index;
484        for (int i = 0; i < numberSolutions; i++)
485            delete cornerPoints[i];
486        delete [] cornerPoints;
487        delete [] rp;
488        return 0;
489    }
490
491    // We found something better
492    solutionValue = bestObj;
493    for (int k = 0; k < numCols; k++) {
494        betterSolution[k] =  feasibles[numFeasibles-1][k];
495    }
496    delete [] varClassInt;
497    for (int i = 0; i < numRows; i++)
498        delete matrix[i];
499    delete [] matrix;
500    delete [] newObj;
501    delete [] index;
502    for (int i = 0; i < numberSolutions; i++)
503        delete cornerPoints[i];
504    delete [] cornerPoints;
505    delete [] rp;
506    std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
507    return 1;
508
509}
510// update model
511void CbcHeuristicRandRound::setModel(CbcModel * model)
512{
513    CbcHeuristic::setModel(model);
514}
515
516
Note: See TracBrowser for help on using the repository browser.