source: stable/2.9/Cbc/src/CbcHeuristicPivotAndFix.cpp @ 2238

/* $Id: CbcHeuristicPivotAndFix.cpp 2094 2014-11-18 11:15:36Z unxusr $ */
// 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 "OsiSolverInterface.hpp"
#include "CbcModel.hpp"
#include "CbcMessage.hpp"
#include "CbcHeuristicPivotAndFix.hpp"
#include "OsiClpSolverInterface.hpp"
#include  "CoinTime.hpp"

//#define FORNOW

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

// Constructor with model - assumed before cuts

CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix(CbcModel & model)
        : CbcHeuristic(model)
{
}

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

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

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

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

// Resets stuff if model changes
void
CbcHeuristicPivotAndFix::resetModel(CbcModel * /*model*/)
{
    //CbcHeuristic::resetModel(model);
}
/*
  Comments needed
  Returns 1 if solution, 0 if not */
int
CbcHeuristicPivotAndFix::solution(double & /*solutionValue*/,
                                  double * /*betterSolution*/)
{

    numCouldRun_++; // Todo: Ask JJHF what this for.
    std::cout << "Entering Pivot-and-Fix Heuristic" << std::endl;

#ifdef HEURISTIC_INFORM
    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
           heuristicName(),numRuns_,numCouldRun_,when_);
#endif
#ifdef FORNOW
    std::cout << "Lucky you! You're in the Pivot-and-Fix Heuristic" << std::endl;
    // The struct should be moved to member data



    typedef struct {
        int numberSolutions;
        int maximumSolutions;
        int numberColumns;
        double ** solution;
        int * numberUnsatisfied;
    } clpSolution;

    double start = CoinCpuTime();

    OsiClpSolverInterface * clpSolverOriginal
    = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
    assert (clpSolverOriginal);
    OsiClpSolverInterface *clpSolver(clpSolverOriginal);

    ClpSimplex * simplex = clpSolver->getModelPtr();

    // Initialize the structure holding the solutions
    clpSolution solutions;
    // Set typeStruct field of ClpTrustedData struct to one.
    // This tells Clp it's "Mahdi!"
    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();
    simplex->primal();

    // -------------------------------------------------
    // Get the problem information
    // - get the number of cols and rows
    int numCols = clpSolver->getNumCols();
    int numRows = clpSolver->getNumRows();

    // - get the right hand side of the rows
    const double * rhs = clpSolver->getRightHandSide();

    // - find the integer variables
    bool * varClassInt = new bool[numCols];
    int numInt = 0;
    for (int i = 0; i < numCols; i++) {
        if (clpSolver->isContinuous(i))
            varClassInt[i] = 0;
        else {
            varClassInt[i] = 1;
            numInt++;
        }
    }

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

    // -Get the objective coefficients
    const double *objCoefficients = clpSolver->getObjCoefficients();
    double *originalObjCoeff = new double [numCols];
    for (int i = 0; i < numCols; i++)
        originalObjCoeff[i] = objCoefficients[i];

    // -Get the matrix of the problem
    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];
        }
    }

    // The newObj is the randomly perturbed constraint used to find new
    // corner points
    double * newObj = new double [numCols];

    // Set the random seed
    srand ( time(NULL) + 1);
    int randNum;

    // We're going to add a new row to the LP formulation
    // after finding each new solution.
    // Adding a new row requires the new elements and the new indices.
    // The elements are original objective function coefficients.
    // The indicies are the (dense) columns indices stored in addRowIndex.
    // The rhs is the value of the new solution stored in solutionValue.
    int * addRowIndex = new int[numCols];
    for (int i = 0; i < numCols; i++)
        addRowIndex[i] = i;

    // The number of feasible solutions found by the PF heuristic.
    // This controls the return code of the solution() method.
    int numFeasibles = 0;

    // Shuffle the rows
    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 + (rand() % (numRows - i));
        index[i] = index[randNumTemp];
        index[randNumTemp] = temp;
    }

    // In the clpSolution struct, we store a lot of column solutions.
    // For each perturb objective, we store the solution from each
    // iteration of the LP solve.
    // For each perturb objective, we look at the collection of
    // solutions to do something extremly intelligent :-)
    // We could (and should..and will :-) wipe out the block of
    // solutions when we're done with them. But for now, we just move on
    // and store the next block of solutions for the next (perturbed)
    // objective.
    // The variable startIndex tells us where the new block begins.
    int startIndex = 0;

    // At most "fixThreshold" number of integer variables can be unsatisfied
    // for calling smallBranchAndBound().
    // The PF Heuristic only fixes fixThreshold number of variables to
    // their integer values. Not more. Not less. The reason is to give
    // the smallBB some opportunity to find better solutions. If we fix
    // everything it might be too many (leading the heuristic to come up
    // with infeasibility rather than a useful result).
    // (This is an important paramater. And it is dynamically set.)
    double fixThreshold;
    /*
        if(numInt > 400)
        fixThreshold = 17*sqrt(numInt);
        if(numInt<=400 && numInt>100)
        fixThreshold = 5*sqrt(numInt);
        if(numInt<=100)
        fixThreshold = 4*sqrt(numInt);
    */
    // Initialize fixThreshold based on the number of integer
    // variables
    if (numInt <= 100)
        fixThreshold = .35 * numInt;
    if (numInt > 100 && numInt < 1000)
        fixThreshold = .85 * numInt;
    if (numInt >= 1000)
        fixThreshold = .1 * numInt;

    // Whenever the dynamic system for changing fixThreshold
    // kicks in, it changes the parameter by the
    // fixThresholdChange amount.
    // (The 25% should be member data and tuned. Another paper!)
    double fixThresholdChange = 0.25 * fixThreshold;

    // maxNode is the maximum number of nodes we allow smallBB to
    // search. It's initialized to 400 and changed dynamically.
    // The 400 should be member data, if we become virtuous.
    int maxNode = 400;

    // We control the decision to change maxNode through the boolean
    // variable  changeMaxNode. The boolean variable is initialized to
    // true and gets set to false under a condition (and is never true
    // again.)
    // It's flipped off and stays off (in the current incarnation of PF)
    bool changeMaxNode = 1;

    // The sumReturnCode is used for the dynamic system that sets
    // fixThreshold and changeMaxNode.
    //
    // We track what's happening in sumReturnCode. There are 8 switches.
    // The first 5 switches corresponds to a return code for smallBB.
    //
    // We want to know how many times we consecutively get the same
    // return code.
    //
    // If "good" return codes are happening often enough, we're happy.
    //
    // If a "bad" returncodes happen consecutively, we want to
    // change something.
    //
    // The switch 5 is the number of times PF didn't call smallBB
    // becuase the number of integer variables that took integer values
    // was less than fixThreshold.
    //
    // The swicth 6 was added for a brilliant idea...to be announced
    // later (another paper!)
    //
    // The switch 7 is the one that changes the max node. Read the
    // code. (Todo: Verbalize the brilliant idea for the masses.)
    //
    int sumReturnCode[8];
    /*
      sumReturnCode[0] ~ -1 --> problem too big for smallBB
      sumReturnCode[1] ~ 0  --> smallBB not finshed and no soln
      sumReturnCode[2] ~ 1  --> smallBB not finshed and there is a soln
      sumReturnCode[3] ~ 2  --> smallBB finished and no soln
      sumReturnCode[4] ~ 3  --> smallBB finished and there is a soln
      sumReturnCode[5] ~ didn't call smallBranchAndBound too few to fix
      sumReturnCode[6] ~ didn't call smallBranchAndBound too many unsatisfied
      sumReturnCode[7] ~ the same as sumReturnCode[1] but becomes zero just if the returnCode is not 0
    */

    for (int i = 0; i < 8; i++)
        sumReturnCode[i] = 0;
    int * colIndex = new int[numCols];
    for (int i = 0; i < numCols; i++)
        colIndex[i] = i;
    double cutoff = COIN_DBL_MAX;
    bool didMiniBB;

    // Main loop
    for (int i = 0; i < numRows; i++) {
        // track the number of mini-bb for the dynamic threshold setting
        didMiniBB = 0;

        for (int k = startIndex; k < solutions.numberSolutions; k++)
            //if the point has 0 unsatisfied variables; make sure it is
            //feasible. Check integer feasiblity and constraints.
            if (solutions.numberUnsatisfied[k] == 0) {
                double feasibility = 1;
                //check integer feasibility
                for (int icol = 0; icol < numCols; icol++) {
                    double closest = floor(solutions.solution[k][icol] + 0.5);
                    if (varClassInt[icol] && (fabs(solutions.solution[k][icol] - closest) > 1e-6)) {
                        feasibility = 0;
                        break;
                    }
                }
                //check if the solution satisfies the constraints
                for (int irow = 0; irow < numRows; irow++) {
                    double lhs = 0;
                    for (int j = 0; j < numCols; j++)
                        lhs += matrix[irow][j] * solutions.solution[k][j];
                    if (rowSense[irow] == 'L' && lhs > rhs[irow] + 1e-6) {
                        feasibility = 0;
                        break;
                    }
                    if (rowSense[irow] == 'G' && lhs < rhs[irow] - 1e-6) {
                        feasibility = 0;
                        break;
                    }
                    if (rowSense[irow] == 'E' && (lhs - rhs[irow] > 1e-6 || lhs - rhs[irow] < -1e-6)) {
                        feasibility = 0;
                        break;
                    }
                }

                //if feasible, find the objective value and set the cutoff
                // for the smallBB and add a new constraint to the LP
                // (and update the best solution found so far for the
                // return arguments)
                if (feasibility) {
                    double objectiveValue = 0;
                    for (int j = 0; j < numCols; j++)
                        objectiveValue += solutions.solution[k][j] * originalObjCoeff[j];
                    cutoff = objectiveValue;
                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);

                    // Todo: pick up the best solution in the block (not
                    // the last).
                    solutionValue = objectiveValue;
                    for (int m = 0; m < numCols; m++)
                        betterSolution[m] = solutions.solution[k][m];
                    numFeasibles++;
                }
            }

        // Go through the block of solution and decide if to call smallBB
        for (int k = startIndex; k < solutions.numberSolutions; k++) {
            if (solutions.numberUnsatisfied[k] <= fixThreshold) {
                // get new copy
                OsiSolverInterface * newSolver;
                newSolver = new OsiClpSolverInterface(*clpSolver);
                newSolver->setObjSense(1);
                newSolver->setObjective(originalObjCoeff);
                int numberColumns = newSolver->getNumCols();
                int numFixed = 0;

                // Fix the first fixThreshold number of integer vars
                // that are satisfied
                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
                    if (newSolver->isInteger(iColumn)) {
                        double value = solutions.solution[k][iColumn];
                        double intValue = floor(value + 0.5);
                        if (fabs(value - intValue) < 1.0e-5) {
                            newSolver->setColLower(iColumn, intValue);
                            newSolver->setColUpper(iColumn, intValue);
                            numFixed++;
                            if (numFixed > numInt - fixThreshold)
                                break;
                        }
                    }
                }
                COIN_DETAIL_PRINT(printf("numFixed: %d\n", numFixed));
                COIN_DETAIL_PRINT(printf("fixThreshold: %f\n", fixThreshold));
                COIN_DETAIL_PRINT(printf("numInt: %d\n", numInt));
                double *newSolution = new double[numCols];
                double newSolutionValue;

                // Call smallBB on the modified problem
                int returnCode = smallBranchAndBound(newSolver, maxNode, newSolution,
                                                     newSolutionValue, cutoff, "mini");

                // If smallBB found a solution, update the better
                // solution and solutionValue (we gave smallBB our
                // cutoff, so it only finds improving solutions)
                if (returnCode == 1 || returnCode == 3) {
                    numFeasibles ++;
                    solutionValue = newSolutionValue;
                    for (int m = 0; m < numCols; m++)
                        betterSolution[m] = newSolution[m];
                    COIN_DETAIL_PRINT(printf("cutoff: %f\n", newSolutionValue));
                    COIN_DETAIL_PRINT(printf("time: %.2lf\n", CoinCpuTime() - start));
                }
                didMiniBB = 1;
                COIN_DETAIL_PRINT(printf("returnCode: %d\n", returnCode));

                //Update sumReturnCode array
                for (int iRC = 0; iRC < 6; iRC++) {
                    if (iRC == returnCode + 1)
                        sumReturnCode[iRC]++;
                    else
                        sumReturnCode[iRC] = 0;
                }
                if (returnCode != 0)
                    sumReturnCode[7] = 0;
                else
                    sumReturnCode[7]++;
                if (returnCode == 1 || returnCode == 3) {
                    cutoff = newSolutionValue;
                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);
                    COIN_DETAIL_PRINT(printf("******************\n\n*****************\n"));
                }
                break;
            }
        }

        if (!didMiniBB && solutions.numberSolutions - startIndex > 0) {
            sumReturnCode[5]++;
            for (int iRC = 0; iRC < 5; iRC++)
                sumReturnCode[iRC] = 0;
        }

        //Change "fixThreshold" if needed
        // using the data we've recorded in sumReturnCode
        if (sumReturnCode[1] >= 3)
            fixThreshold -= fixThresholdChange;
        if (sumReturnCode[7] >= 3 && changeMaxNode) {
            maxNode *= 5;
            changeMaxNode = 0;
        }
        if (sumReturnCode[3] >= 3 && fixThreshold < 0.95 * numInt)
            fixThreshold += fixThresholdChange;
        if (sumReturnCode[5] >= 4)
            fixThreshold += fixThresholdChange;
        if (sumReturnCode[0] > 3)
            fixThreshold -= fixThresholdChange;

        startIndex = solutions.numberSolutions;

        //Check if the maximum iterations limit is reached
        // rlh: Ask John how this is working with the change to trustedUserPtr.
        if (solutions.numberSolutions > 20000)
            break;

        // The first time in this loop PF solves orig LP.

        //Generate the random objective function
        randNum = rand() % 10 + 1;
        randNum = fmod(randNum, 2);
        for (int j = 0; j < numCols; j++) {
            if (randNum == 1)
                if (fabs(matrix[index[i]][j]) < 1e-6)
                    newObj[j] = 0.1;
                else
                    newObj[j] = matrix[index[i]][j] * 1.1;
            else if (fabs(matrix[index[i]][j]) < 1e-6)
                newObj[j] = -0.1;
            else
                newObj[j] = matrix[index[i]][j] * 0.9;
        }
        clpSolver->setObjective(newObj);
        if (rowSense[i] == 'L')
            clpSolver->setObjSense(-1);
        else
            // Todo #1: We don't need to solve the LPs to optimality.
            // We just need corner points.
            // There's a problem in stopping Clp that needs to be looked
            // into. So for now, we solve optimality.
            clpSolver->setObjSense(1);
        //        simplex->setMaximumIterations(100);
        clpSolver->getModelPtr()->primal(1);
        //        simplex->setMaximumIterations(100000);
#ifdef COIN_DETAIL
        printf("cutoff: %f\n", cutoff);
        printf("time: %.2f\n", CoinCpuTime() - start);
        for (int iRC = 0; iRC < 8; iRC++)
            printf("%d ", sumReturnCode[iRC]);
        printf("\nfixThreshold: %f\n", fixThreshold);
        printf("numInt: %d\n", numInt);
        printf("\n---------------------------------------------------------------- %d\n", i);
#endif

        //temp:
        if (i > 3) break;

    }

    COIN_DETAIL_PRINT(printf("Best Feasible Found: %f\n", cutoff));
    COIN_DETAIL_PRINT(printf("Total time: %.2f\n", CoinCpuTime() - start));

    if (numFeasibles == 0) {
        return 0;
    }



    // We found something better
    std::cout << "See you soon! You're leaving the Pivot-and-Fix Heuristic" << std::endl;
    std::cout << std::endl;

    return 1;
#endif

    return 0;

}




// update model
void CbcHeuristicPivotAndFix::setModel(CbcModel * )
{
    // probably same as resetModel
}