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

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

Change to EPL license notice.

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