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