File: CbcHeuristicRandRound.cpp

package info (click to toggle)
coinor-cbc 2.9.9+repack1-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 7,848 kB
  • ctags: 5,787
  • sloc: cpp: 104,337; sh: 8,921; xml: 2,950; makefile: 520; ansic: 491; awk: 197
file content (522 lines) | stat: -rw-r--r-- 17,786 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
/* $Id: CbcHeuristicRandRound.cpp 2094 2014-11-18 11:15:36Z forrest $ */
// Copyright (C) 2008, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <cfloat>
#include <vector>

#include "CoinHelperFunctions.hpp"
#include "OsiSolverInterface.hpp"
#include "CbcModel.hpp"
#include "CbcMessage.hpp"
#include "CbcHeuristicRandRound.hpp"
#include "OsiClpSolverInterface.hpp"
#include  "CoinTime.hpp"

static inline int intRand(const int range)
{
    return static_cast<int> (floor(CoinDrand48() * range));
}

// Default Constructor
CbcHeuristicRandRound::CbcHeuristicRandRound()
        : CbcHeuristic()
{
}

// Constructor with model - assumed before cuts
CbcHeuristicRandRound::CbcHeuristicRandRound(CbcModel & model)
        : CbcHeuristic(model)
{
    model_ = &model;
    setWhen(1);
}

// Destructor
CbcHeuristicRandRound::~CbcHeuristicRandRound ()
{
}

// Clone
CbcHeuristic *
CbcHeuristicRandRound::clone() const
{
    return new CbcHeuristicRandRound(*this);
}

// Create C++ lines to get to current state
void
CbcHeuristicRandRound::generateCpp( FILE * fp)
{
    CbcHeuristicRandRound other;
    fprintf(fp, "0#include \"CbcHeuristicRandRound.hpp\"\n");
    fprintf(fp, "3  CbcHeuristicRandRound heuristicPFX(*cbcModel);\n");
    CbcHeuristic::generateCpp(fp, "heuristicPFX");
    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicPFX);\n");
}

// Copy constructor
CbcHeuristicRandRound::CbcHeuristicRandRound(const CbcHeuristicRandRound & rhs)
        :
        CbcHeuristic(rhs)
{
}

// Assignment operator
CbcHeuristicRandRound &
CbcHeuristicRandRound::operator=( const CbcHeuristicRandRound & rhs)
{
    if (this != &rhs) {
        CbcHeuristic::operator=(rhs);
    }
    return *this;
}

// Resets stuff if model changes
void
CbcHeuristicRandRound::resetModel(CbcModel * model)
{
    CbcHeuristic::resetModel(model);
}

/*
  Randomized Rounding Heuristic
  Returns 1 if solution, 0 if not
*/
int
CbcHeuristicRandRound::solution(double & solutionValue,
                                double * betterSolution)
{
    // rlh: Todo: Memory Cleanup

    //  std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;

    setWhen(1);  // setWhen(1) didn't have the effect I expected (e.g., run once).

    // Run only once.
    //
    //    See if at root node
    bool atRoot = model_->getNodeCount() == 0;
    int passNumber = model_->getCurrentPassNumber();
    //    Just do once
    if (!atRoot || passNumber > 1) {
        // std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
        return 0;
    }

    std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
    typedef struct {
        int numberSolutions;
        int maximumSolutions;
        int numberColumns;
        double ** solution;
        int * numberUnsatisfied;
    } clpSolution;

    double start = CoinCpuTime();
    numCouldRun_++; //
#ifdef HEURISTIC_INFORM
    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
	   heuristicName(),numRuns_,numCouldRun_,when_);
#endif
    // Todo: Ask JJHF what "number of times
    // the heuristic could run" means.

    OsiSolverInterface * solver = model_->solver()->clone();
    double primalTolerance ;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
    OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *> (solver);
    assert (clpSolver);
    ClpSimplex * simplex = clpSolver->getModelPtr();

    // Initialize the structure holding the solutions for the Simplex iterations
    clpSolution solutions;
    // Set typeStruct field of ClpTrustedData struct to 1 to indicate
    // desired behavior for  RandRound heuristic (which is what?)
    ClpTrustedData trustedSolutions;
    trustedSolutions.typeStruct = 1;
    trustedSolutions.data = &solutions;
    solutions.numberSolutions = 0;
    solutions.maximumSolutions = 0;
    solutions.numberColumns = simplex->numberColumns();
    solutions.solution = NULL;
    solutions.numberUnsatisfied = NULL;
    simplex->setTrustedUserPointer(&trustedSolutions);

    // Solve from all slack to get some points
    simplex->allSlackBasis();

    // Calling primal() invalidates pointers to some rim vectors,
    // like...row sense (!)
    simplex->primal();

    // 1. Okay - so a workaround would be to copy the data I want BEFORE
    // calling primal.
    // 2. Another approach is to ask the simplex solvers NOT to mess up my
    // rims.
    // 3. See freeCachedResults() for what is getting
    // deleted. Everything else points into the structure.
    // ...or use collower and colupper rather than rowsense.
    // ..store address of where one of these

    // Store the basic problem information
    // -Get the number of columns, rows and rhs vector
    int numCols = clpSolver->getNumCols();
    int numRows = clpSolver->getNumRows();

    // Find the integer variables (use columnType(?))
    // One if not continuous, that is binary or general integer)
    // columnType() = 0 continuous
    //              = 1 binary
    //              = 2 general integer
    bool * varClassInt = new bool[numCols];
    const char* columnType = clpSolver->columnType();
    int numGenInt = 0;
    for (int i = 0; i < numCols; i++) {
        if (clpSolver->isContinuous(i))
            varClassInt[i] = 0;
        else
            varClassInt[i] = 1;
        if (columnType[i] == 2) numGenInt++;
    }

    // Heuristic is for problems with general integer variables.
    // If there are none, quit.
    if (numGenInt++ < 1) {
        delete [] varClassInt ;
        std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
        return 0;
    }


    // -Get the rows sense
    const char * rowSense;
    rowSense = clpSolver->getRowSense();

    // -Get the objective coefficients
    double *originalObjCoeff = CoinCopyOfArray(clpSolver->getObjCoefficients(), numCols);

    // -Get the matrix of the problem
    // rlh: look at using sparse representation
    double ** matrix = new double * [numRows];
    for (int i = 0; i < numRows; i++) {
        matrix[i] = new double[numCols];
        for (int j = 0; j < numCols; j++)
            matrix[i][j] = 0;
    }

    const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
    const double * matrixElements = matrixByRow->getElements();
    const int * matrixIndices = matrixByRow->getIndices();
    const int * matrixStarts = matrixByRow->getVectorStarts();
    for (int j = 0; j < numRows; j++) {
        for (int i = matrixStarts[j]; i < matrixStarts[j+1]; i++) {
            matrix[j][matrixIndices[i]] = matrixElements[i];
        }
    }

    double * newObj = new double [numCols];
    srand ( static_cast<unsigned int>(time(NULL) + 1));
    int randNum;

    // Shuffle the rows:
    // Put the rows in a random order
    // so that the optimal solution is a different corner point than the
    // starting point.
    int * index = new int [numRows];
    for (int i = 0; i < numRows; i++)
        index[i] = i;
    for (int i = 0; i < numRows; i++) {
        int temp = index[i];
        int randNumTemp = i + intRand(numRows - i);
        index[i] = index[randNumTemp];
        index[randNumTemp] = temp;
    }

    // Start finding corner points by iteratively doing the following:
    // - contruct a randomly tilted objective
    // - solve
    for (int i = 0; i < numRows; i++) {
        // TODO: that 10,000 could be a param in the member data
        if (solutions.numberSolutions  > 10000)
            break;
        randNum = intRand(2);
        for (int j = 0; j < numCols; j++) {
            // for row i and column j vary the coefficient "a bit"
            if (randNum == 1)
                // if the element is zero, then set the new obj
                // coefficient to 0.1 (i.e., round up)
                if (fabs(matrix[index[i]][j]) < primalTolerance)
                    newObj[j] = 0.1;
                else
                    // if the element is nonzero, then increase the new obj
                    // coefficient "a bit"
                    newObj[j] = matrix[index[i]][j] * 1.1;
            else
                // if randnum is 2, then
                // if the element is zero, then set the new obj coeffient
                // to NEGATIVE 0.1 (i.e., round down)
                if (fabs(matrix[index[i]][j]) < primalTolerance)
                    newObj[j] = -0.1;
                else
                    // if the element is nonzero, then DEcrease the new obj coeffienct "a bit"
                    newObj[j] = matrix[index[i]][j] * 0.9;
        }
        // Use the new "tilted" objective
        clpSolver->setObjective(newObj);

        // Based on the row sense, we decide whether to max or min
        if (rowSense[i] == 'L')
            clpSolver->setObjSense(-1);
        else
            clpSolver->setObjSense(1);

        // Solve with primal simplex
        simplex->primal(1);
        // rlh+ll: This was the original code. But we already have the
        // model pointer (it's in simplex). And, calling getModelPtr()
        // invalidates the cached data in the OsiClpSolverInterface
        // object, which means our precious rowsens is lost. So let's
        // not use the line below...
        /******* clpSolver->getModelPtr()->primal(1); */
        printf("---------------------------------------------------------------- %d\n", i);
    }
    // Iteratively do this process until...
    // either you reach the max number of corner points (aka 10K)
    // or all the rows have been used as an objective.

    // Look at solutions
    int numberSolutions = solutions.numberSolutions;
    //const char * integerInfo = simplex->integerInformation();
    //const double * columnLower = simplex->columnLower();
    //const double * columnUpper = simplex->columnUpper();
    printf("there are %d solutions\n", numberSolutions);

    // Up to here we have all the corner points
    // Now we need to do the random walks and roundings

    double ** cornerPoints = new double * [numberSolutions];
    for (int j = 0; j < numberSolutions; j++)
        cornerPoints[j] = solutions.solution[j];

    bool feasibility = 1;
    // rlh: use some COIN max instead of 1e30 (?)
    double bestObj = 1e30;
    std::vector< std::vector <double> > feasibles;
    int numFeasibles = 0;

    // Check the feasibility of the corner points
    int numCornerPoints = numberSolutions;

    const double * rhs = clpSolver->getRightHandSide();
    // rlh: row sense hasn't changed. why a fresh copy?
    // Delete next line.
    rowSense = clpSolver->getRowSense();

    for (int i = 0; i < numCornerPoints; i++) {
        //get the objective value for this this point
        double objValue = 0;
        for (int k = 0; k < numCols; k++)
            objValue += cornerPoints[i][k] * originalObjCoeff[k];

        if (objValue < bestObj) {
            // check integer feasibility
            feasibility = 1;
            for (int j = 0; j < numCols; j++) {
                if (varClassInt[j]) {
                    double closest = floor(cornerPoints[i][j] + 0.5);
                    if (fabs(cornerPoints[i][j] - closest) > primalTolerance) {
                        feasibility = 0;
                        break;
                    }
                }
            }
            // check all constraints satisfied
            if (feasibility) {
                for (int irow = 0; irow < numRows; irow++) {
                    double lhs = 0;
                    for (int j = 0; j < numCols; j++) {
                        lhs += matrix[irow][j] * cornerPoints[i][j];
                    }
                    if (rowSense[irow] == 'L' && lhs > rhs[irow] + primalTolerance) {
                        feasibility = 0;
                        break;
                    }
                    if (rowSense[irow] == 'G' && lhs < rhs[irow] - primalTolerance) {
                        feasibility = 0;
                        break;
                    }
                    if (rowSense[irow] == 'E' && (lhs - rhs[irow] > primalTolerance || lhs - rhs[irow] < -primalTolerance)) {
                        feasibility = 0;
                        break;
                    }
                }
            }

            if (feasibility) {
                numFeasibles++;
                feasibles.push_back(std::vector <double> (numCols));
                for (int k = 0; k < numCols; k++)
                    feasibles[numFeasibles-1][k] = cornerPoints[i][k];
                printf("obj: %f\n", objValue);
                if (objValue < bestObj)
                    bestObj = objValue;
            }
        }
    }
    int numFeasibleCorners;
    numFeasibleCorners = numFeasibles;
    //find the center of gravity of the corner points as the first random point
    double * rp = new double[numCols];
    for (int i = 0; i < numCols; i++) {
        rp[i] = 0;
        for (int j = 0; j < numCornerPoints; j++) {
            rp[i] += cornerPoints[j][i];
        }
        rp[i] = rp[i] / numCornerPoints;
    }

    //-------------------------------------------
    //main loop:
    // -generate the next random point
    // -round the random point
    // -check the feasibility of the random point
    //-------------------------------------------

    srand ( static_cast<unsigned int>(time(NULL) + 1));
    int numRandomPoints = 0;
    while (numRandomPoints < 50000) {
        numRandomPoints++;
        //generate the next random point
        int randomIndex = intRand(numCornerPoints);
        double random = CoinDrand48();
        for (int i = 0; i < numCols; i++) {
            rp[i] = (random * (cornerPoints[randomIndex][i] - rp[i])) + rp[i];
        }

        //CRISP ROUNDING
        //round the random point just generated
        double * roundRp = new double[numCols];
        for (int i = 0; i < numCols; i++) {
            roundRp[i] = rp[i];
            if (varClassInt[i]) {
                if (rp[i] >= 0) {
                    if (fmod(rp[i], 1) > 0.5)
                        roundRp[i] = floor(rp[i]) + 1;
                    else
                        roundRp[i] = floor(rp[i]);
                } else {
                    if (fabs(fmod(rp[i], 1)) > 0.5)
                        roundRp[i] = floor(rp[i]);
                    else
                        roundRp[i] = floor(rp[i]) + 1;

                }
            }
        }


        //SOFT ROUNDING
        // Look at original files for the "how to" on soft rounding;
        // Soft rounding omitted here.

        //Check the feasibility of the rounded random point
        // -Check the feasibility
        // -Get the rows sense
        rowSense = clpSolver->getRowSense();
        rhs = clpSolver->getRightHandSide();

        //get the objective value for this feasible point
        double objValue = 0;
        for (int i = 0; i < numCols; i++)
            objValue += roundRp[i] * originalObjCoeff[i];

        if (objValue < bestObj) {
            feasibility = 1;
            for (int i = 0; i < numRows; i++) {
                double lhs = 0;
                for (int j = 0; j < numCols; j++) {
                    lhs += matrix[i][j] * roundRp[j];
                }
                if (rowSense[i] == 'L' && lhs > rhs[i] + primalTolerance) {
                    feasibility = 0;
                    break;
                }
                if (rowSense[i] == 'G' && lhs < rhs[i] - primalTolerance) {
                    feasibility = 0;
                    break;
                }
                if (rowSense[i] == 'E' && (lhs - rhs[i] > primalTolerance || lhs - rhs[i] < -primalTolerance)) {
                    feasibility = 0;
                    break;
                }
            }
            if (feasibility) {
                printf("Feasible Found.\n");
                printf("%.2f\n", CoinCpuTime() - start);
                numFeasibles++;
                feasibles.push_back(std::vector <double> (numCols));
                for (int i = 0; i < numCols; i++)
                    feasibles[numFeasibles-1][i] = roundRp[i];
                printf("obj: %f\n", objValue);
                if (objValue < bestObj)
                    bestObj = objValue;
            }
        }
        delete [] roundRp;
    }
    printf("Number of Feasible Corners: %d\n", numFeasibleCorners);
    printf("Number of Feasibles Found: %d\n", numFeasibles);
    if (numFeasibles > 0)
        printf("Best Objective: %f\n", bestObj);
    printf("time: %.2f\n", CoinCpuTime() - start);

    if (numFeasibles == 0) {
        // cleanup
        delete [] varClassInt;
        for (int i = 0; i < numRows; i++)
            delete matrix[i];
        delete [] matrix;
        delete [] newObj;
        delete [] index;
        for (int i = 0; i < numberSolutions; i++)
            delete cornerPoints[i];
        delete [] cornerPoints;
        delete [] rp;
        return 0;
    }

    // We found something better
    solutionValue = bestObj;
    for (int k = 0; k < numCols; k++) {
        betterSolution[k] =  feasibles[numFeasibles-1][k];
    }
    delete [] varClassInt;
    for (int i = 0; i < numRows; i++)
        delete matrix[i];
    delete [] matrix;
    delete [] newObj;
    delete [] index;
    for (int i = 0; i < numberSolutions; i++)
        delete cornerPoints[i];
    delete [] cornerPoints;
    delete [] rp;
    std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
    return 1;

}
// update model
void CbcHeuristicRandRound::setModel(CbcModel * model)
{
    CbcHeuristic::setModel(model);
}