Ignore:
Timestamp:
Nov 9, 2009 6:33:07 PM (10 years ago)
Author:
EdwinStraver
Message:

Changed formatting using AStyle -A4 -p

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/Cbc/src/CbcHeuristicGreedy.cpp

    r1271 r1286  
    1818#include "CglPreProcess.hpp"
    1919// Default Constructor
    20 CbcHeuristicGreedyCover::CbcHeuristicGreedyCover() 
    21   :CbcHeuristic()
    22 {
    23   // matrix  will automatically be empty
    24   originalNumberRows_=0;
    25   algorithm_=0;
    26   numberTimes_=100;
     20CbcHeuristicGreedyCover::CbcHeuristicGreedyCover()
     21        : CbcHeuristic()
     22{
     23    // matrix  will automatically be empty
     24    originalNumberRows_ = 0;
     25    algorithm_ = 0;
     26    numberTimes_ = 100;
    2727}
    2828
    2929// Constructor from model
    3030CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(CbcModel & model)
    31   :CbcHeuristic(model)
    32 {
    33   gutsOfConstructor(&model);
    34   algorithm_=0;
    35   numberTimes_=100;
    36   whereFrom_=1;
    37 }
    38 
    39 // Destructor 
     31        : CbcHeuristic(model)
     32{
     33    gutsOfConstructor(&model);
     34    algorithm_ = 0;
     35    numberTimes_ = 100;
     36    whereFrom_ = 1;
     37}
     38
     39// Destructor
    4040CbcHeuristicGreedyCover::~CbcHeuristicGreedyCover ()
    4141{
     
    4646CbcHeuristicGreedyCover::clone() const
    4747{
    48   return new CbcHeuristicGreedyCover(*this);
     48    return new CbcHeuristicGreedyCover(*this);
    4949}
    5050// Guts of constructor from a CbcModel
    51 void 
     51void
    5252CbcHeuristicGreedyCover::gutsOfConstructor(CbcModel * model)
    5353{
    54   model_=model;
    55   // Get a copy of original matrix
    56   assert(model->solver());
    57   if (model->solver()->getNumRows()) {
    58     matrix_ = *model->solver()->getMatrixByCol();
    59   }
    60   originalNumberRows_=model->solver()->getNumRows();
     54    model_ = model;
     55    // Get a copy of original matrix
     56    assert(model->solver());
     57    if (model->solver()->getNumRows()) {
     58        matrix_ = *model->solver()->getMatrixByCol();
     59    }
     60    originalNumberRows_ = model->solver()->getNumRows();
    6161}
    6262// Create C++ lines to get to current state
    63 void 
    64 CbcHeuristicGreedyCover::generateCpp( FILE * fp) 
    65 {
    66   CbcHeuristicGreedyCover other;
    67   fprintf(fp,"0#include \"CbcHeuristicGreedy.hpp\"\n");
    68   fprintf(fp,"3  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n");
    69   CbcHeuristic::generateCpp(fp,"heuristicGreedyCover");
    70   if (algorithm_!=other.algorithm_)
    71     fprintf(fp,"3  heuristicGreedyCover.setAlgorithm(%d);\n",algorithm_);
    72   else
    73     fprintf(fp,"4  heuristicGreedyCover.setAlgorithm(%d);\n",algorithm_);
    74   if (numberTimes_!=other.numberTimes_)
    75     fprintf(fp,"3  heuristicGreedyCover.setNumberTimes(%d);\n",numberTimes_);
    76   else
    77     fprintf(fp,"4  heuristicGreedyCover.setNumberTimes(%d);\n",numberTimes_);
    78   fprintf(fp,"3  cbcModel->addHeuristic(&heuristicGreedyCover);\n");
    79 }
    80 
    81 // Copy constructor 
     63void
     64CbcHeuristicGreedyCover::generateCpp( FILE * fp)
     65{
     66    CbcHeuristicGreedyCover other;
     67    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
     68    fprintf(fp, "3  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n");
     69    CbcHeuristic::generateCpp(fp, "heuristicGreedyCover");
     70    if (algorithm_ != other.algorithm_)
     71        fprintf(fp, "3  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
     72    else
     73        fprintf(fp, "4  heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_);
     74    if (numberTimes_ != other.numberTimes_)
     75        fprintf(fp, "3  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
     76    else
     77        fprintf(fp, "4  heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_);
     78    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyCover);\n");
     79}
     80
     81// Copy constructor
    8282CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(const CbcHeuristicGreedyCover & rhs)
    83 :
    84   CbcHeuristic(rhs),
    85   matrix_(rhs.matrix_),
    86   originalNumberRows_(rhs.originalNumberRows_),
    87   algorithm_(rhs.algorithm_),
    88   numberTimes_(rhs.numberTimes_)
    89 {
    90 }
    91 
    92 // Assignment operator 
    93 CbcHeuristicGreedyCover & 
    94 CbcHeuristicGreedyCover::operator=( const CbcHeuristicGreedyCover& rhs)
    95 {
    96   if (this != &rhs) {
    97     CbcHeuristic::operator=(rhs);
    98     matrix_=rhs.matrix_;
    99     originalNumberRows_=rhs.originalNumberRows_;
    100     algorithm_=rhs.algorithm_;
    101     numberTimes_=rhs.numberTimes_;
    102   }
    103   return *this;
     83        :
     84        CbcHeuristic(rhs),
     85        matrix_(rhs.matrix_),
     86        originalNumberRows_(rhs.originalNumberRows_),
     87        algorithm_(rhs.algorithm_),
     88        numberTimes_(rhs.numberTimes_)
     89{
     90}
     91
     92// Assignment operator
     93CbcHeuristicGreedyCover &
     94CbcHeuristicGreedyCover::operator=( const CbcHeuristicGreedyCover & rhs)
     95{
     96    if (this != &rhs) {
     97        CbcHeuristic::operator=(rhs);
     98        matrix_ = rhs.matrix_;
     99        originalNumberRows_ = rhs.originalNumberRows_;
     100        algorithm_ = rhs.algorithm_;
     101        numberTimes_ = rhs.numberTimes_;
     102    }
     103    return *this;
    104104}
    105105// Returns 1 if solution, 0 if not
    106106int
    107107CbcHeuristicGreedyCover::solution(double & solutionValue,
    108                          double * betterSolution)
    109 {
    110   numCouldRun_++;
    111   if (!model_)
    112     return 0;
    113   // See if to do
    114   if (!when()||(when()==1&&model_->phase()!=1))
    115     return 0; // switched off
    116   if (model_->getNodeCount()>numberTimes_)
    117     return 0;
    118   // See if at root node
    119   bool atRoot = model_->getNodeCount()==0;
    120   int passNumber = model_->getCurrentPassNumber();
    121   if (atRoot&&passNumber!=1)
    122     return 0;
    123   OsiSolverInterface * solver = model_->solver();
    124   const double * columnLower = solver->getColLower();
    125   const double * columnUpper = solver->getColUpper();
    126   // And original upper bounds in case we want to use them
    127   const double * originalUpper = model_->continuousSolver()->getColUpper();
    128   // But not if algorithm says so
    129   if ((algorithm_%10)==0)
    130     originalUpper = columnUpper;
    131   const double * rowLower = solver->getRowLower();
    132   const double * solution = solver->getColSolution();
    133   const double * objective = solver->getObjCoefficients();
    134   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    135   double primalTolerance;
    136   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    137 
    138   // This is number of rows when matrix was passed in
    139   int numberRows = originalNumberRows_;
    140   if (!numberRows)
    141     return 0; // switched off
    142 
    143   numRuns_++;
    144   assert (numberRows==matrix_.getNumRows());
    145   int iRow, iColumn;
    146   double direction = solver->getObjSense();
    147   double offset;
    148   solver->getDblParam(OsiObjOffset,offset);
    149   double newSolutionValue = -offset;
    150   int returnCode = 0;
    151 
    152   // Column copy
    153   const double * element = matrix_.getElements();
    154   const int * row = matrix_.getIndices();
    155   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    156   const int * columnLength = matrix_.getVectorLengths();
    157 
    158   // Get solution array for heuristic solution
    159   int numberColumns = solver->getNumCols();
    160   double * newSolution = new double [numberColumns];
    161   double * rowActivity = new double[numberRows];
    162   memset(rowActivity,0,numberRows*sizeof(double));
    163   bool allOnes=true;
    164   // Get rounded down solution
    165   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    166     CoinBigIndex j;
    167     double value = solution[iColumn];
    168     if (solver->isInteger(iColumn)) {
    169       // Round down integer
    170       if (fabs(floor(value+0.5)-value)<integerTolerance) {
    171         value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
    172       } else {
    173         value=CoinMax(floor(value),columnLower[iColumn]);
    174       }
    175     }
    176     // make sure clean
    177     value = CoinMin(value,columnUpper[iColumn]);
    178     value = CoinMax(value,columnLower[iColumn]);
    179     newSolution[iColumn]=value;
    180     double cost = direction * objective[iColumn];
    181     newSolutionValue += value*cost;
    182     for (j=columnStart[iColumn];
    183          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    184       int iRow=row[j];
    185       rowActivity[iRow] += value*element[j];
    186       if (element[j]!=1.0)
    187         allOnes=false;
    188     }
    189   }
    190   // See if we round up
    191   bool roundup = ((algorithm_%100)!=0);
    192   if (roundup&&allOnes) {
    193     // Get rounded up solution
    194     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    195       CoinBigIndex j;
    196       double value = solution[iColumn];
    197       if (solver->isInteger(iColumn)) {
    198         // but round up if no activity
    199         if (roundup&&value>=0.499999&&!newSolution[iColumn]) {
    200           bool choose=true;
    201           for (j=columnStart[iColumn];
    202                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    203             int iRow=row[j];
    204             if (rowActivity[iRow]) {
    205               choose=false;
    206               break;
    207             }
    208           }
    209           if (choose) {
    210             newSolution[iColumn]=1.0;
     108                                  double * betterSolution)
     109{
     110    numCouldRun_++;
     111    if (!model_)
     112        return 0;
     113    // See if to do
     114    if (!when() || (when() == 1 && model_->phase() != 1))
     115        return 0; // switched off
     116    if (model_->getNodeCount() > numberTimes_)
     117        return 0;
     118    // See if at root node
     119    bool atRoot = model_->getNodeCount() == 0;
     120    int passNumber = model_->getCurrentPassNumber();
     121    if (atRoot && passNumber != 1)
     122        return 0;
     123    OsiSolverInterface * solver = model_->solver();
     124    const double * columnLower = solver->getColLower();
     125    const double * columnUpper = solver->getColUpper();
     126    // And original upper bounds in case we want to use them
     127    const double * originalUpper = model_->continuousSolver()->getColUpper();
     128    // But not if algorithm says so
     129    if ((algorithm_ % 10) == 0)
     130        originalUpper = columnUpper;
     131    const double * rowLower = solver->getRowLower();
     132    const double * solution = solver->getColSolution();
     133    const double * objective = solver->getObjCoefficients();
     134    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     135    double primalTolerance;
     136    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
     137
     138    // This is number of rows when matrix was passed in
     139    int numberRows = originalNumberRows_;
     140    if (!numberRows)
     141        return 0; // switched off
     142
     143    numRuns_++;
     144    assert (numberRows == matrix_.getNumRows());
     145    int iRow, iColumn;
     146    double direction = solver->getObjSense();
     147    double offset;
     148    solver->getDblParam(OsiObjOffset, offset);
     149    double newSolutionValue = -offset;
     150    int returnCode = 0;
     151
     152    // Column copy
     153    const double * element = matrix_.getElements();
     154    const int * row = matrix_.getIndices();
     155    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
     156    const int * columnLength = matrix_.getVectorLengths();
     157
     158    // Get solution array for heuristic solution
     159    int numberColumns = solver->getNumCols();
     160    double * newSolution = new double [numberColumns];
     161    double * rowActivity = new double[numberRows];
     162    memset(rowActivity, 0, numberRows*sizeof(double));
     163    bool allOnes = true;
     164    // Get rounded down solution
     165    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     166        CoinBigIndex j;
     167        double value = solution[iColumn];
     168        if (solver->isInteger(iColumn)) {
     169            // Round down integer
     170            if (fabs(floor(value + 0.5) - value) < integerTolerance) {
     171                value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
     172            } else {
     173                value = CoinMax(floor(value), columnLower[iColumn]);
     174            }
     175        }
     176        // make sure clean
     177        value = CoinMin(value, columnUpper[iColumn]);
     178        value = CoinMax(value, columnLower[iColumn]);
     179        newSolution[iColumn] = value;
     180        double cost = direction * objective[iColumn];
     181        newSolutionValue += value * cost;
     182        for (j = columnStart[iColumn];
     183                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     184            int iRow = row[j];
     185            rowActivity[iRow] += value * element[j];
     186            if (element[j] != 1.0)
     187                allOnes = false;
     188        }
     189    }
     190    // See if we round up
     191    bool roundup = ((algorithm_ % 100) != 0);
     192    if (roundup && allOnes) {
     193        // Get rounded up solution
     194        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     195            CoinBigIndex j;
     196            double value = solution[iColumn];
     197            if (solver->isInteger(iColumn)) {
     198                // but round up if no activity
     199                if (roundup && value >= 0.499999 && !newSolution[iColumn]) {
     200                    bool choose = true;
     201                    for (j = columnStart[iColumn];
     202                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     203                        int iRow = row[j];
     204                        if (rowActivity[iRow]) {
     205                            choose = false;
     206                            break;
     207                        }
     208                    }
     209                    if (choose) {
     210                        newSolution[iColumn] = 1.0;
     211                        double cost = direction * objective[iColumn];
     212                        newSolutionValue += cost;
     213                        for (j = columnStart[iColumn];
     214                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     215                            int iRow = row[j];
     216                            rowActivity[iRow] += 1.0;
     217                        }
     218                    }
     219                }
     220            }
     221        }
     222    }
     223    // Get initial list
     224    int * which = new int [numberColumns];
     225    for (iColumn = 0; iColumn < numberColumns; iColumn++)
     226        which[iColumn] = iColumn;
     227    int numberLook = numberColumns;
     228    // See if we want to perturb more
     229    double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
     230    // Keep going round until a solution
     231    while (true) {
     232        // Get column with best ratio
     233        int bestColumn = -1;
     234        double bestRatio = COIN_DBL_MAX;
     235        double bestStepSize = 0.0;
     236        int newNumber = 0;
     237        for (int jColumn = 0; jColumn < numberLook; jColumn++) {
     238            int iColumn = which[jColumn];
     239            CoinBigIndex j;
     240            double value = newSolution[iColumn];
    211241            double cost = direction * objective[iColumn];
    212             newSolutionValue += cost;
    213             for (j=columnStart[iColumn];
    214                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    215               int iRow=row[j];
    216               rowActivity[iRow] += 1.0;
    217             }
    218           }
    219         }
    220       }
    221     }
    222   }
    223   // Get initial list
    224   int * which = new int [numberColumns];
    225   for (iColumn=0;iColumn<numberColumns;iColumn++)
    226     which[iColumn]=iColumn;
    227   int numberLook=numberColumns;
    228   // See if we want to perturb more
    229   double perturb = ((algorithm_%10)==0) ? 0.1 : 0.25;
    230   // Keep going round until a solution
    231   while (true) {
    232     // Get column with best ratio
    233     int bestColumn=-1;
    234     double bestRatio=COIN_DBL_MAX;
    235     double bestStepSize = 0.0;
    236     int newNumber=0;
    237     for (int jColumn=0;jColumn<numberLook;jColumn++) {
    238       int iColumn = which[jColumn];
    239       CoinBigIndex j;
    240       double value = newSolution[iColumn];
    241       double cost = direction * objective[iColumn];
    242       if (solver->isInteger(iColumn)) {
    243         // use current upper or original upper
    244         if (value+0.99<originalUpper[iColumn]) {
    245           double sum=0.0;
    246           int numberExact=0;
    247           for (j=columnStart[iColumn];
    248                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    249             int iRow=row[j];
    250             double gap = rowLower[iRow]-rowActivity[iRow];
    251             double elementValue = allOnes ? 1.0 : element[j];
    252             if (gap>1.0e-7) {
    253               sum += CoinMin(elementValue,gap);
    254               if (fabs(elementValue-gap)<1.0e-7)
    255                 numberExact++;
    256             }
    257           }
    258           // could bias if exact
    259           if (sum>0.0) {
    260             // add to next time
    261             which[newNumber++]=iColumn;
    262             double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
    263             // If at root choose first
    264             if (atRoot)
    265               ratio = iColumn;
    266             if (ratio<bestRatio) {
    267               bestRatio=ratio;
    268               bestColumn=iColumn;
    269               bestStepSize=1.0;
    270             }
    271           }
    272         }
    273       } else {
    274         // continuous
    275         if (value<columnUpper[iColumn]) {
    276           // Go through twice - first to get step length
    277           double step=1.0e50;
    278           for (j=columnStart[iColumn];
    279                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    280             int iRow=row[j];
    281             if (rowActivity[iRow]<rowLower[iRow]-1.0e-10&&
    282                 element[j]*step+rowActivity[iRow]>=rowLower[iRow]) {
    283               step = (rowLower[iRow]-rowActivity[iRow])/element[j];;
    284             }
    285           }
    286           // now ratio
    287           if (step<1.0e50) {
    288             // add to next time
    289             which[newNumber++]=iColumn;
    290             assert (step>0.0);
    291             double sum=0.0;
    292             for (j=columnStart[iColumn];
    293                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    294               int iRow=row[j];
    295               double newActivity = element[j]*step+rowActivity[iRow];
    296               if (rowActivity[iRow]<rowLower[iRow]-1.0e-10&&
    297                 newActivity>=rowLower[iRow]-1.0e-12) {
    298                 sum += element[j];
    299               }
    300             }
    301             assert (sum>0.0);
    302             double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
    303             if (ratio<bestRatio) {
    304               bestRatio=ratio;
    305               bestColumn=iColumn;
    306               bestStepSize=step;
    307             }
    308           }
    309         }
    310       }
    311     }
    312     if (bestColumn<0)
    313       break; // we have finished
    314     // Increase chosen column
    315     newSolution[bestColumn] += bestStepSize;
    316     double cost = direction * objective[bestColumn];
    317     newSolutionValue += bestStepSize*cost;
    318     for (CoinBigIndex j=columnStart[bestColumn];
    319          j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
    320       int iRow = row[j];
    321       rowActivity[iRow] += bestStepSize*element[j];
    322     }
    323   }
    324   delete [] which;
    325   if (newSolutionValue<solutionValue) {
    326     // check feasible
    327     memset(rowActivity,0,numberRows*sizeof(double));
    328     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    329       CoinBigIndex j;
    330       double value = newSolution[iColumn];
    331       if (value) {
    332         for (j=columnStart[iColumn];
    333              j<columnStart[iColumn]+columnLength[iColumn];j++) {
    334           int iRow=row[j];
    335           rowActivity[iRow] += value*element[j];
    336         }
    337       }
    338     }
    339     // check was approximately feasible
    340     bool feasible=true;
    341     for (iRow=0;iRow<numberRows;iRow++) {
    342       if(rowActivity[iRow]<rowLower[iRow]) {
    343         if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
    344           feasible = false;
    345       }
    346     }
    347     if (feasible) {
    348       // new solution
    349       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    350       solutionValue = newSolutionValue;
    351       //printf("** Solution of %g found by rounding\n",newSolutionValue);
    352       returnCode=1;
    353     } else {
    354       // Can easily happen
    355       //printf("Debug CbcHeuristicGreedyCover giving bad solution\n");
    356     }
    357   }
    358   delete [] newSolution;
    359   delete [] rowActivity;
    360   return returnCode;
     242            if (solver->isInteger(iColumn)) {
     243                // use current upper or original upper
     244                if (value + 0.99 < originalUpper[iColumn]) {
     245                    double sum = 0.0;
     246                    int numberExact = 0;
     247                    for (j = columnStart[iColumn];
     248                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     249                        int iRow = row[j];
     250                        double gap = rowLower[iRow] - rowActivity[iRow];
     251                        double elementValue = allOnes ? 1.0 : element[j];
     252                        if (gap > 1.0e-7) {
     253                            sum += CoinMin(elementValue, gap);
     254                            if (fabs(elementValue - gap) < 1.0e-7)
     255                                numberExact++;
     256                        }
     257                    }
     258                    // could bias if exact
     259                    if (sum > 0.0) {
     260                        // add to next time
     261                        which[newNumber++] = iColumn;
     262                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
     263                        // If at root choose first
     264                        if (atRoot)
     265                            ratio = iColumn;
     266                        if (ratio < bestRatio) {
     267                            bestRatio = ratio;
     268                            bestColumn = iColumn;
     269                            bestStepSize = 1.0;
     270                        }
     271                    }
     272                }
     273            } else {
     274                // continuous
     275                if (value < columnUpper[iColumn]) {
     276                    // Go through twice - first to get step length
     277                    double step = 1.0e50;
     278                    for (j = columnStart[iColumn];
     279                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     280                        int iRow = row[j];
     281                        if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 &&
     282                                element[j]*step + rowActivity[iRow] >= rowLower[iRow]) {
     283                            step = (rowLower[iRow] - rowActivity[iRow]) / element[j];;
     284                        }
     285                    }
     286                    // now ratio
     287                    if (step < 1.0e50) {
     288                        // add to next time
     289                        which[newNumber++] = iColumn;
     290                        assert (step > 0.0);
     291                        double sum = 0.0;
     292                        for (j = columnStart[iColumn];
     293                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     294                            int iRow = row[j];
     295                            double newActivity = element[j] * step + rowActivity[iRow];
     296                            if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 &&
     297                                    newActivity >= rowLower[iRow] - 1.0e-12) {
     298                                sum += element[j];
     299                            }
     300                        }
     301                        assert (sum > 0.0);
     302                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
     303                        if (ratio < bestRatio) {
     304                            bestRatio = ratio;
     305                            bestColumn = iColumn;
     306                            bestStepSize = step;
     307                        }
     308                    }
     309                }
     310            }
     311        }
     312        if (bestColumn < 0)
     313            break; // we have finished
     314        // Increase chosen column
     315        newSolution[bestColumn] += bestStepSize;
     316        double cost = direction * objective[bestColumn];
     317        newSolutionValue += bestStepSize * cost;
     318        for (CoinBigIndex j = columnStart[bestColumn];
     319                j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
     320            int iRow = row[j];
     321            rowActivity[iRow] += bestStepSize * element[j];
     322        }
     323    }
     324    delete [] which;
     325    if (newSolutionValue < solutionValue) {
     326        // check feasible
     327        memset(rowActivity, 0, numberRows*sizeof(double));
     328        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     329            CoinBigIndex j;
     330            double value = newSolution[iColumn];
     331            if (value) {
     332                for (j = columnStart[iColumn];
     333                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     334                    int iRow = row[j];
     335                    rowActivity[iRow] += value * element[j];
     336                }
     337            }
     338        }
     339        // check was approximately feasible
     340        bool feasible = true;
     341        for (iRow = 0; iRow < numberRows; iRow++) {
     342            if (rowActivity[iRow] < rowLower[iRow]) {
     343                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
     344                    feasible = false;
     345            }
     346        }
     347        if (feasible) {
     348            // new solution
     349            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
     350            solutionValue = newSolutionValue;
     351            //printf("** Solution of %g found by rounding\n",newSolutionValue);
     352            returnCode = 1;
     353        } else {
     354            // Can easily happen
     355            //printf("Debug CbcHeuristicGreedyCover giving bad solution\n");
     356        }
     357    }
     358    delete [] newSolution;
     359    delete [] rowActivity;
     360    return returnCode;
    361361}
    362362// update model
    363363void CbcHeuristicGreedyCover::setModel(CbcModel * model)
    364364{
    365   gutsOfConstructor(model);
    366   validate();
     365    gutsOfConstructor(model);
     366    validate();
    367367}
    368368// Resets stuff if model changes
    369 void 
     369void
    370370CbcHeuristicGreedyCover::resetModel(CbcModel * model)
    371371{
    372   gutsOfConstructor(model);
     372    gutsOfConstructor(model);
    373373}
    374374// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    375 void 
    376 CbcHeuristicGreedyCover::validate() 
    377 {
    378   if (model_&&when()<10) {
    379     if (model_->numberIntegers()!=
    380         model_->numberObjects()&&(model_->numberObjects()||
    381                                   (model_->specialOptions()&1024)==0)) {
    382       int numberOdd=0;
    383       for (int i=0;i<model_->numberObjects();i++) {
    384         if (!model_->object(i)->canDoHeuristics())
    385           numberOdd++;
    386       }
    387       if (numberOdd)
    388         setWhen(0);
    389     }
    390     // Only works if costs positive, coefficients positive and all rows G
    391     OsiSolverInterface * solver = model_->solver();
    392     const double * columnLower = solver->getColLower();
    393     const double * rowUpper = solver->getRowUpper();
    394     const double * objective = solver->getObjCoefficients();
    395     double direction = solver->getObjSense();
    396 
    397     int numberRows = solver->getNumRows();
    398     // Column copy
    399     const double * element = matrix_.getElements();
    400     const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    401     const int * columnLength = matrix_.getVectorLengths();
    402     bool good = true;
    403     for (int iRow=0;iRow<numberRows;iRow++) {
    404       if (rowUpper[iRow]<1.0e30)
    405         good = false;
    406     }
    407     int numberColumns = solver->getNumCols();
    408     for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    409       if (objective[iColumn]*direction<0.0)
    410         good=false;
    411       if (columnLower[iColumn]<0.0)
    412         good=false;
    413       CoinBigIndex j;
    414       for (j=columnStart[iColumn];
    415            j<columnStart[iColumn]+columnLength[iColumn];j++) {
    416         if (element[j]<0.0)
    417           good=false;
    418       }
    419     }
    420     if (!good)
    421       setWhen(0); // switch off
    422   }
     375void
     376CbcHeuristicGreedyCover::validate()
     377{
     378    if (model_ && when() < 10) {
     379        if (model_->numberIntegers() !=
     380                model_->numberObjects() && (model_->numberObjects() ||
     381                                            (model_->specialOptions()&1024) == 0)) {
     382            int numberOdd = 0;
     383            for (int i = 0; i < model_->numberObjects(); i++) {
     384                if (!model_->object(i)->canDoHeuristics())
     385                    numberOdd++;
     386            }
     387            if (numberOdd)
     388                setWhen(0);
     389        }
     390        // Only works if costs positive, coefficients positive and all rows G
     391        OsiSolverInterface * solver = model_->solver();
     392        const double * columnLower = solver->getColLower();
     393        const double * rowUpper = solver->getRowUpper();
     394        const double * objective = solver->getObjCoefficients();
     395        double direction = solver->getObjSense();
     396
     397        int numberRows = solver->getNumRows();
     398        // Column copy
     399        const double * element = matrix_.getElements();
     400        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
     401        const int * columnLength = matrix_.getVectorLengths();
     402        bool good = true;
     403        for (int iRow = 0; iRow < numberRows; iRow++) {
     404            if (rowUpper[iRow] < 1.0e30)
     405                good = false;
     406        }
     407        int numberColumns = solver->getNumCols();
     408        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     409            if (objective[iColumn]*direction < 0.0)
     410                good = false;
     411            if (columnLower[iColumn] < 0.0)
     412                good = false;
     413            CoinBigIndex j;
     414            for (j = columnStart[iColumn];
     415                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     416                if (element[j] < 0.0)
     417                    good = false;
     418            }
     419        }
     420        if (!good)
     421            setWhen(0); // switch off
     422    }
    423423}
    424424// Default Constructor
    425 CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality() 
    426   :CbcHeuristic()
    427 {
    428   // matrix  will automatically be empty
    429   fraction_=1.0; // no branch and bound
    430   originalNumberRows_=0;
    431   algorithm_=0;
    432   numberTimes_=100;
    433   whereFrom_=1;
     425CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality()
     426        : CbcHeuristic()
     427{
     428    // matrix  will automatically be empty
     429    fraction_ = 1.0; // no branch and bound
     430    originalNumberRows_ = 0;
     431    algorithm_ = 0;
     432    numberTimes_ = 100;
     433    whereFrom_ = 1;
    434434}
    435435
    436436// Constructor from model
    437437CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel & model)
    438   :CbcHeuristic(model)
    439 {
    440   // Get a copy of original matrix
    441   gutsOfConstructor(&model);
    442   fraction_=1.0; // no branch and bound
    443   algorithm_=0;
    444   numberTimes_=100;
    445   whereFrom_=1;
    446 }
    447 
    448 // Destructor 
     438        : CbcHeuristic(model)
     439{
     440    // Get a copy of original matrix
     441    gutsOfConstructor(&model);
     442    fraction_ = 1.0; // no branch and bound
     443    algorithm_ = 0;
     444    numberTimes_ = 100;
     445    whereFrom_ = 1;
     446}
     447
     448// Destructor
    449449CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality ()
    450450{
     
    455455CbcHeuristicGreedyEquality::clone() const
    456456{
    457   return new CbcHeuristicGreedyEquality(*this);
     457    return new CbcHeuristicGreedyEquality(*this);
    458458}
    459459// Guts of constructor from a CbcModel
    460 void 
     460void
    461461CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel * model)
    462462{
    463   model_=model;
    464   // Get a copy of original matrix
    465   assert(model->solver());
    466   if (model->solver()->getNumRows()) {
    467     matrix_ = *model->solver()->getMatrixByCol();
    468   }
    469   originalNumberRows_=model->solver()->getNumRows();
     463    model_ = model;
     464    // Get a copy of original matrix
     465    assert(model->solver());
     466    if (model->solver()->getNumRows()) {
     467        matrix_ = *model->solver()->getMatrixByCol();
     468    }
     469    originalNumberRows_ = model->solver()->getNumRows();
    470470}
    471471// Create C++ lines to get to current state
    472 void 
    473 CbcHeuristicGreedyEquality::generateCpp( FILE * fp) 
    474 {
    475   CbcHeuristicGreedyEquality other;
    476   fprintf(fp,"0#include \"CbcHeuristicGreedy.hpp\"\n");
    477   fprintf(fp,"3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
    478   CbcHeuristic::generateCpp(fp,"heuristicGreedyEquality");
    479   if (algorithm_!=other.algorithm_)
    480     fprintf(fp,"3  heuristicGreedyEquality.setAlgorithm(%d);\n",algorithm_);
    481   else
    482     fprintf(fp,"4  heuristicGreedyEquality.setAlgorithm(%d);\n",algorithm_);
    483   if (fraction_!=other.fraction_)
    484     fprintf(fp,"3  heuristicGreedyEquality.setFraction(%g);\n",fraction_);
    485   else
    486     fprintf(fp,"4  heuristicGreedyEquality.setFraction(%g);\n",fraction_);
    487   if (numberTimes_!=other.numberTimes_)
    488     fprintf(fp,"3  heuristicGreedyEquality.setNumberTimes(%d);\n",numberTimes_);
    489   else
    490     fprintf(fp,"4  heuristicGreedyEquality.setNumberTimes(%d);\n",numberTimes_);
    491   fprintf(fp,"3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
    492 }
    493 
    494 // Copy constructor 
     472void
     473CbcHeuristicGreedyEquality::generateCpp( FILE * fp)
     474{
     475    CbcHeuristicGreedyEquality other;
     476    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
     477    fprintf(fp, "3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
     478    CbcHeuristic::generateCpp(fp, "heuristicGreedyEquality");
     479    if (algorithm_ != other.algorithm_)
     480        fprintf(fp, "3  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
     481    else
     482        fprintf(fp, "4  heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_);
     483    if (fraction_ != other.fraction_)
     484        fprintf(fp, "3  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
     485    else
     486        fprintf(fp, "4  heuristicGreedyEquality.setFraction(%g);\n", fraction_);
     487    if (numberTimes_ != other.numberTimes_)
     488        fprintf(fp, "3  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
     489    else
     490        fprintf(fp, "4  heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_);
     491    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
     492}
     493
     494// Copy constructor
    495495CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality & rhs)
    496 :
    497   CbcHeuristic(rhs),
    498   matrix_(rhs.matrix_),
    499   fraction_(rhs.fraction_),
    500   originalNumberRows_(rhs.originalNumberRows_),
    501   algorithm_(rhs.algorithm_),
    502   numberTimes_(rhs.numberTimes_)
    503 {
    504 }
    505 
    506 // Assignment operator 
    507 CbcHeuristicGreedyEquality & 
    508 CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality& rhs)
    509 {
    510   if (this != &rhs) {
    511     CbcHeuristic::operator=(rhs);
    512     matrix_=rhs.matrix_;
    513     fraction_ = rhs.fraction_;
    514     originalNumberRows_=rhs.originalNumberRows_;
    515     algorithm_=rhs.algorithm_;
    516     numberTimes_=rhs.numberTimes_;
    517   }
    518   return *this;
     496        :
     497        CbcHeuristic(rhs),
     498        matrix_(rhs.matrix_),
     499        fraction_(rhs.fraction_),
     500        originalNumberRows_(rhs.originalNumberRows_),
     501        algorithm_(rhs.algorithm_),
     502        numberTimes_(rhs.numberTimes_)
     503{
     504}
     505
     506// Assignment operator
     507CbcHeuristicGreedyEquality &
     508CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality & rhs)
     509{
     510    if (this != &rhs) {
     511        CbcHeuristic::operator=(rhs);
     512        matrix_ = rhs.matrix_;
     513        fraction_ = rhs.fraction_;
     514        originalNumberRows_ = rhs.originalNumberRows_;
     515        algorithm_ = rhs.algorithm_;
     516        numberTimes_ = rhs.numberTimes_;
     517    }
     518    return *this;
    519519}
    520520// Returns 1 if solution, 0 if not
    521521int
    522522CbcHeuristicGreedyEquality::solution(double & solutionValue,
    523                          double * betterSolution)
    524 {
    525   numCouldRun_++;
    526   if (!model_)
    527     return 0;
    528   // See if to do
    529   if (!when()||(when()==1&&model_->phase()!=1))
    530     return 0; // switched off
    531   if (model_->getNodeCount()>numberTimes_)
    532     return 0;
    533   // See if at root node
    534   bool atRoot = model_->getNodeCount()==0;
    535   int passNumber = model_->getCurrentPassNumber();
    536   if (atRoot&&passNumber!=1)
    537     return 0;
    538   OsiSolverInterface * solver = model_->solver();
    539   const double * columnLower = solver->getColLower();
    540   const double * columnUpper = solver->getColUpper();
    541   // And original upper bounds in case we want to use them
    542   const double * originalUpper = model_->continuousSolver()->getColUpper();
    543   // But not if algorithm says so
    544   if ((algorithm_%10)==0)
    545     originalUpper = columnUpper;
    546   const double * rowLower = solver->getRowLower();
    547   const double * rowUpper = solver->getRowUpper();
    548   const double * solution = solver->getColSolution();
    549   const double * objective = solver->getObjCoefficients();
    550   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    551   double primalTolerance;
    552   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    553 
    554   // This is number of rows when matrix was passed in
    555   int numberRows = originalNumberRows_;
    556   if (!numberRows)
    557     return 0; // switched off
    558   numRuns_++;
    559 
    560   assert (numberRows==matrix_.getNumRows());
    561   int iRow, iColumn;
    562   double direction = solver->getObjSense();
    563   double offset;
    564   solver->getDblParam(OsiObjOffset,offset);
    565   double newSolutionValue = -offset;
    566   int returnCode = 0;
    567 
    568   // Column copy
    569   const double * element = matrix_.getElements();
    570   const int * row = matrix_.getIndices();
    571   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    572   const int * columnLength = matrix_.getVectorLengths();
    573 
    574   // Get solution array for heuristic solution
    575   int numberColumns = solver->getNumCols();
    576   double * newSolution = new double [numberColumns];
    577   double * rowActivity = new double[numberRows];
    578   memset(rowActivity,0,numberRows*sizeof(double));
    579   double rhsNeeded=0;
    580   for (iRow=0;iRow<numberRows;iRow++)
    581     rhsNeeded += rowUpper[iRow];
    582   rhsNeeded *= fraction_;
    583   bool allOnes=true;
    584   // Get rounded down solution
    585   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    586     CoinBigIndex j;
    587     double value = solution[iColumn];
    588     if (solver->isInteger(iColumn)) {
    589       // Round down integer
    590       if (fabs(floor(value+0.5)-value)<integerTolerance) {
    591         value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
    592       } else {
    593         value=CoinMax(floor(value),columnLower[iColumn]);
    594       }
    595     }
    596     // make sure clean
    597     value = CoinMin(value,columnUpper[iColumn]);
    598     value = CoinMax(value,columnLower[iColumn]);
    599     newSolution[iColumn]=value;
    600     double cost = direction * objective[iColumn];
    601     newSolutionValue += value*cost;
    602     for (j=columnStart[iColumn];
    603          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    604       int iRow=row[j];
    605       rowActivity[iRow] += value*element[j];
    606       rhsNeeded -= value*element[j];
    607       if (element[j]!=1.0)
    608         allOnes=false;
    609     }
    610   }
    611   // See if we round up
    612   bool roundup = ((algorithm_%100)!=0);
    613   if (roundup&&allOnes) {
    614     // Get rounded up solution
    615     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    616       CoinBigIndex j;
    617       double value = solution[iColumn];
    618       if (solver->isInteger(iColumn)) {
    619         // but round up if no activity
    620         if (roundup&&value>=0.6&&!newSolution[iColumn]) {
    621           bool choose=true;
    622           for (j=columnStart[iColumn];
    623                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    624             int iRow=row[j];
    625             if (rowActivity[iRow]) {
    626               choose=false;
    627               break;
    628             }
    629           }
    630           if (choose) {
    631             newSolution[iColumn]=1.0;
     523                                     double * betterSolution)
     524{
     525    numCouldRun_++;
     526    if (!model_)
     527        return 0;
     528    // See if to do
     529    if (!when() || (when() == 1 && model_->phase() != 1))
     530        return 0; // switched off
     531    if (model_->getNodeCount() > numberTimes_)
     532        return 0;
     533    // See if at root node
     534    bool atRoot = model_->getNodeCount() == 0;
     535    int passNumber = model_->getCurrentPassNumber();
     536    if (atRoot && passNumber != 1)
     537        return 0;
     538    OsiSolverInterface * solver = model_->solver();
     539    const double * columnLower = solver->getColLower();
     540    const double * columnUpper = solver->getColUpper();
     541    // And original upper bounds in case we want to use them
     542    const double * originalUpper = model_->continuousSolver()->getColUpper();
     543    // But not if algorithm says so
     544    if ((algorithm_ % 10) == 0)
     545        originalUpper = columnUpper;
     546    const double * rowLower = solver->getRowLower();
     547    const double * rowUpper = solver->getRowUpper();
     548    const double * solution = solver->getColSolution();
     549    const double * objective = solver->getObjCoefficients();
     550    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     551    double primalTolerance;
     552    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
     553
     554    // This is number of rows when matrix was passed in
     555    int numberRows = originalNumberRows_;
     556    if (!numberRows)
     557        return 0; // switched off
     558    numRuns_++;
     559
     560    assert (numberRows == matrix_.getNumRows());
     561    int iRow, iColumn;
     562    double direction = solver->getObjSense();
     563    double offset;
     564    solver->getDblParam(OsiObjOffset, offset);
     565    double newSolutionValue = -offset;
     566    int returnCode = 0;
     567
     568    // Column copy
     569    const double * element = matrix_.getElements();
     570    const int * row = matrix_.getIndices();
     571    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
     572    const int * columnLength = matrix_.getVectorLengths();
     573
     574    // Get solution array for heuristic solution
     575    int numberColumns = solver->getNumCols();
     576    double * newSolution = new double [numberColumns];
     577    double * rowActivity = new double[numberRows];
     578    memset(rowActivity, 0, numberRows*sizeof(double));
     579    double rhsNeeded = 0;
     580    for (iRow = 0; iRow < numberRows; iRow++)
     581        rhsNeeded += rowUpper[iRow];
     582    rhsNeeded *= fraction_;
     583    bool allOnes = true;
     584    // Get rounded down solution
     585    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     586        CoinBigIndex j;
     587        double value = solution[iColumn];
     588        if (solver->isInteger(iColumn)) {
     589            // Round down integer
     590            if (fabs(floor(value + 0.5) - value) < integerTolerance) {
     591                value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
     592            } else {
     593                value = CoinMax(floor(value), columnLower[iColumn]);
     594            }
     595        }
     596        // make sure clean
     597        value = CoinMin(value, columnUpper[iColumn]);
     598        value = CoinMax(value, columnLower[iColumn]);
     599        newSolution[iColumn] = value;
     600        double cost = direction * objective[iColumn];
     601        newSolutionValue += value * cost;
     602        for (j = columnStart[iColumn];
     603                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     604            int iRow = row[j];
     605            rowActivity[iRow] += value * element[j];
     606            rhsNeeded -= value * element[j];
     607            if (element[j] != 1.0)
     608                allOnes = false;
     609        }
     610    }
     611    // See if we round up
     612    bool roundup = ((algorithm_ % 100) != 0);
     613    if (roundup && allOnes) {
     614        // Get rounded up solution
     615        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     616            CoinBigIndex j;
     617            double value = solution[iColumn];
     618            if (solver->isInteger(iColumn)) {
     619                // but round up if no activity
     620                if (roundup && value >= 0.6 && !newSolution[iColumn]) {
     621                    bool choose = true;
     622                    for (j = columnStart[iColumn];
     623                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     624                        int iRow = row[j];
     625                        if (rowActivity[iRow]) {
     626                            choose = false;
     627                            break;
     628                        }
     629                    }
     630                    if (choose) {
     631                        newSolution[iColumn] = 1.0;
     632                        double cost = direction * objective[iColumn];
     633                        newSolutionValue += cost;
     634                        for (j = columnStart[iColumn];
     635                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     636                            int iRow = row[j];
     637                            rowActivity[iRow] += 1.0;
     638                            rhsNeeded -= 1.0;
     639                        }
     640                    }
     641                }
     642            }
     643        }
     644    }
     645    // Get initial list
     646    int * which = new int [numberColumns];
     647    for (iColumn = 0; iColumn < numberColumns; iColumn++)
     648        which[iColumn] = iColumn;
     649    int numberLook = numberColumns;
     650    // See if we want to perturb more
     651    double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25;
     652    // Keep going round until a solution
     653    while (true) {
     654        // Get column with best ratio
     655        int bestColumn = -1;
     656        double bestRatio = COIN_DBL_MAX;
     657        double bestStepSize = 0.0;
     658        int newNumber = 0;
     659        for (int jColumn = 0; jColumn < numberLook; jColumn++) {
     660            int iColumn = which[jColumn];
     661            CoinBigIndex j;
     662            double value = newSolution[iColumn];
    632663            double cost = direction * objective[iColumn];
    633             newSolutionValue += cost;
    634             for (j=columnStart[iColumn];
    635                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    636               int iRow=row[j];
    637               rowActivity[iRow] += 1.0;
    638               rhsNeeded -= 1.0;
    639             }
    640           }
    641         }
    642       }
    643     }
    644   }
    645   // Get initial list
    646   int * which = new int [numberColumns];
    647   for (iColumn=0;iColumn<numberColumns;iColumn++)
    648     which[iColumn]=iColumn;
    649   int numberLook=numberColumns;
    650   // See if we want to perturb more
    651   double perturb = ((algorithm_%10)==0) ? 0.1 : 0.25;
    652   // Keep going round until a solution
    653   while (true) {
    654     // Get column with best ratio
    655     int bestColumn=-1;
    656     double bestRatio=COIN_DBL_MAX;
    657     double bestStepSize = 0.0;
    658     int newNumber=0;
    659     for (int jColumn=0;jColumn<numberLook;jColumn++) {
    660       int iColumn = which[jColumn];
    661       CoinBigIndex j;
    662       double value = newSolution[iColumn];
    663       double cost = direction * objective[iColumn];
    664       if (solver->isInteger(iColumn)) {
    665         // use current upper or original upper
    666         if (value+0.9999<originalUpper[iColumn]) {
    667           double movement = 1.0;
    668           double sum=0.0;
    669           for (j=columnStart[iColumn];
    670                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    671             int iRow=row[j];
    672             double gap = rowUpper[iRow]-rowActivity[iRow];
    673             double elementValue = allOnes ? 1.0 : element[j];
    674             sum += elementValue;
    675             if (movement*elementValue>gap) {
    676               movement = gap/elementValue;
    677             }
    678           }
    679           if (movement>0.999999) {
    680             // add to next time
    681             which[newNumber++]=iColumn;
    682             double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
    683             // If at root
    684             if (atRoot) {
    685               if (fraction_==1.0)
    686                 ratio = iColumn; // choose first
    687               else
    688                 ratio = - solution[iColumn]; // choose largest
    689             }
    690             if (ratio<bestRatio) {
    691               bestRatio=ratio;
    692               bestColumn=iColumn;
    693               bestStepSize=1.0;
    694             }
    695           }
    696         }
    697       } else {
    698         // continuous
    699         if (value<columnUpper[iColumn]) {
    700           double movement=1.0e50;
    701           double sum=0.0;
    702           for (j=columnStart[iColumn];
    703                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    704             int iRow=row[j];
    705             if (element[j]*movement+rowActivity[iRow]>rowUpper[iRow]) {
    706               movement = (rowUpper[iRow]-rowActivity[iRow])/element[j];;
    707             }
    708             sum += element[j];
    709           }
    710           // now ratio
    711           if (movement>1.0e-7) {
    712             // add to next time
    713             which[newNumber++]=iColumn;
    714             double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
    715             if (ratio<bestRatio) {
    716               bestRatio=ratio;
    717               bestColumn=iColumn;
    718               bestStepSize=movement;
    719             }
    720           }
    721         }
    722       }
    723     }
    724     if (bestColumn<0)
    725       break; // we have finished
    726     // Increase chosen column
    727     newSolution[bestColumn] += bestStepSize;
    728     double cost = direction * objective[bestColumn];
    729     newSolutionValue += bestStepSize*cost;
    730     for (CoinBigIndex j=columnStart[bestColumn];
    731          j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
    732       int iRow = row[j];
    733       rowActivity[iRow] += bestStepSize*element[j];
    734       rhsNeeded -= bestStepSize*element[j];
    735     }
    736     if (rhsNeeded<1.0e-8)
    737       break;
    738   }
    739   delete [] which;
    740   if (fraction_<1.0&&rhsNeeded<1.0e-8&&newSolutionValue<solutionValue) {
    741     // do branch and cut
    742     // fix all nonzero
    743     OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
    744     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    745       if (newSolver->isInteger(iColumn))
    746         newSolver->setColLower(iColumn,newSolution[iColumn]);
    747     }
    748     int returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
    749                                          solutionValue,"CbcHeuristicGreedy");
    750     if (returnCode<0)
    751       returnCode=0; // returned on size
    752     if ((returnCode&2)!=0) {
    753       // could add cut
    754       returnCode &= ~2;
    755     }
    756     rhsNeeded = 1.0-returnCode;
    757     delete newSolver;
    758   }
    759   if (newSolutionValue<solutionValue&&rhsNeeded<1.0e-8) {
    760     // check feasible
    761     memset(rowActivity,0,numberRows*sizeof(double));
    762     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    763       CoinBigIndex j;
    764       double value = newSolution[iColumn];
    765       if (value) {
    766         for (j=columnStart[iColumn];
    767              j<columnStart[iColumn]+columnLength[iColumn];j++) {
    768           int iRow=row[j];
    769           rowActivity[iRow] += value*element[j];
    770         }
    771       }
    772     }
    773     // check was approximately feasible
    774     bool feasible=true;
    775     for (iRow=0;iRow<numberRows;iRow++) {
    776       if(rowActivity[iRow]<rowLower[iRow]) {
    777         if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
    778           feasible = false;
    779       }
    780     }
    781     if (feasible) {
    782       // new solution
    783       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    784       solutionValue = newSolutionValue;
    785       returnCode=1;
    786     }
    787   }
    788   delete [] newSolution;
    789   delete [] rowActivity;
    790   if (atRoot&&fraction_==1.0) {
    791     // try quick search
    792     fraction_=0.4;
    793     int newCode=this->solution(solutionValue,betterSolution);
    794     if (newCode)
    795       returnCode=1;
    796     fraction_=1.0;
    797   }
    798   return returnCode;
     664            if (solver->isInteger(iColumn)) {
     665                // use current upper or original upper
     666                if (value + 0.9999 < originalUpper[iColumn]) {
     667                    double movement = 1.0;
     668                    double sum = 0.0;
     669                    for (j = columnStart[iColumn];
     670                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     671                        int iRow = row[j];
     672                        double gap = rowUpper[iRow] - rowActivity[iRow];
     673                        double elementValue = allOnes ? 1.0 : element[j];
     674                        sum += elementValue;
     675                        if (movement*elementValue > gap) {
     676                            movement = gap / elementValue;
     677                        }
     678                    }
     679                    if (movement > 0.999999) {
     680                        // add to next time
     681                        which[newNumber++] = iColumn;
     682                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
     683                        // If at root
     684                        if (atRoot) {
     685                            if (fraction_ == 1.0)
     686                                ratio = iColumn; // choose first
     687                            else
     688                                ratio = - solution[iColumn]; // choose largest
     689                        }
     690                        if (ratio < bestRatio) {
     691                            bestRatio = ratio;
     692                            bestColumn = iColumn;
     693                            bestStepSize = 1.0;
     694                        }
     695                    }
     696                }
     697            } else {
     698                // continuous
     699                if (value < columnUpper[iColumn]) {
     700                    double movement = 1.0e50;
     701                    double sum = 0.0;
     702                    for (j = columnStart[iColumn];
     703                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     704                        int iRow = row[j];
     705                        if (element[j]*movement + rowActivity[iRow] > rowUpper[iRow]) {
     706                            movement = (rowUpper[iRow] - rowActivity[iRow]) / element[j];;
     707                        }
     708                        sum += element[j];
     709                    }
     710                    // now ratio
     711                    if (movement > 1.0e-7) {
     712                        // add to next time
     713                        which[newNumber++] = iColumn;
     714                        double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble());
     715                        if (ratio < bestRatio) {
     716                            bestRatio = ratio;
     717                            bestColumn = iColumn;
     718                            bestStepSize = movement;
     719                        }
     720                    }
     721                }
     722            }
     723        }
     724        if (bestColumn < 0)
     725            break; // we have finished
     726        // Increase chosen column
     727        newSolution[bestColumn] += bestStepSize;
     728        double cost = direction * objective[bestColumn];
     729        newSolutionValue += bestStepSize * cost;
     730        for (CoinBigIndex j = columnStart[bestColumn];
     731                j < columnStart[bestColumn] + columnLength[bestColumn]; j++) {
     732            int iRow = row[j];
     733            rowActivity[iRow] += bestStepSize * element[j];
     734            rhsNeeded -= bestStepSize * element[j];
     735        }
     736        if (rhsNeeded < 1.0e-8)
     737            break;
     738    }
     739    delete [] which;
     740    if (fraction_ < 1.0 && rhsNeeded < 1.0e-8 && newSolutionValue < solutionValue) {
     741        // do branch and cut
     742        // fix all nonzero
     743        OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
     744        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     745            if (newSolver->isInteger(iColumn))
     746                newSolver->setColLower(iColumn, newSolution[iColumn]);
     747        }
     748        int returnCode = smallBranchAndBound(newSolver, 200, newSolution, newSolutionValue,
     749                                             solutionValue, "CbcHeuristicGreedy");
     750        if (returnCode < 0)
     751            returnCode = 0; // returned on size
     752        if ((returnCode&2) != 0) {
     753            // could add cut
     754            returnCode &= ~2;
     755        }
     756        rhsNeeded = 1.0 - returnCode;
     757        delete newSolver;
     758    }
     759    if (newSolutionValue < solutionValue && rhsNeeded < 1.0e-8) {
     760        // check feasible
     761        memset(rowActivity, 0, numberRows*sizeof(double));
     762        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     763            CoinBigIndex j;
     764            double value = newSolution[iColumn];
     765            if (value) {
     766                for (j = columnStart[iColumn];
     767                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     768                    int iRow = row[j];
     769                    rowActivity[iRow] += value * element[j];
     770                }
     771            }
     772        }
     773        // check was approximately feasible
     774        bool feasible = true;
     775        for (iRow = 0; iRow < numberRows; iRow++) {
     776            if (rowActivity[iRow] < rowLower[iRow]) {
     777                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
     778                    feasible = false;
     779            }
     780        }
     781        if (feasible) {
     782            // new solution
     783            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
     784            solutionValue = newSolutionValue;
     785            returnCode = 1;
     786        }
     787    }
     788    delete [] newSolution;
     789    delete [] rowActivity;
     790    if (atRoot && fraction_ == 1.0) {
     791        // try quick search
     792        fraction_ = 0.4;
     793        int newCode = this->solution(solutionValue, betterSolution);
     794        if (newCode)
     795            returnCode = 1;
     796        fraction_ = 1.0;
     797    }
     798    return returnCode;
    799799}
    800800// update model
    801801void CbcHeuristicGreedyEquality::setModel(CbcModel * model)
    802802{
    803   gutsOfConstructor(model);
    804   validate();
     803    gutsOfConstructor(model);
     804    validate();
    805805}
    806806// Resets stuff if model changes
    807 void 
     807void
    808808CbcHeuristicGreedyEquality::resetModel(CbcModel * model)
    809809{
    810   gutsOfConstructor(model);
     810    gutsOfConstructor(model);
    811811}
    812812// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    813 void 
    814 CbcHeuristicGreedyEquality::validate() 
    815 {
    816   if (model_&&when()<10) {
    817     if (model_->numberIntegers()!=
    818         model_->numberObjects())
    819       setWhen(0);
    820     // Only works if costs positive, coefficients positive and all rows E or L
    821     // And if values are integer
    822     OsiSolverInterface * solver = model_->solver();
    823     const double * columnLower = solver->getColLower();
    824     const double * rowUpper = solver->getRowUpper();
    825     const double * rowLower = solver->getRowLower();
    826     const double * objective = solver->getObjCoefficients();
    827     double direction = solver->getObjSense();
    828 
    829     int numberRows = solver->getNumRows();
    830     // Column copy
    831     const double * element = matrix_.getElements();
    832     const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    833     const int * columnLength = matrix_.getVectorLengths();
    834     bool good = true;
    835     for (int iRow=0;iRow<numberRows;iRow++) {
    836       if (rowUpper[iRow]>1.0e30)
    837         good = false;
    838       if (rowLower[iRow]>0.0&&rowLower[iRow]!=rowUpper[iRow])
    839         good = false;
    840       if (floor(rowUpper[iRow]+0.5)!=rowUpper[iRow])
    841         good = false;
    842     }
    843     int numberColumns = solver->getNumCols();
    844     for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    845       if (objective[iColumn]*direction<0.0)
    846         good=false;
    847       if (columnLower[iColumn]<0.0)
    848         good=false;
    849       CoinBigIndex j;
    850       for (j=columnStart[iColumn];
    851            j<columnStart[iColumn]+columnLength[iColumn];j++) {
    852         if (element[j]<0.0)
    853           good=false;
    854         if (floor(element[j]+0.5)!=element[j])
    855           good = false;
    856       }
    857     }
    858     if (!good)
    859       setWhen(0); // switch off
    860   }
    861 }
    862 
    863  
     813void
     814CbcHeuristicGreedyEquality::validate()
     815{
     816    if (model_ && when() < 10) {
     817        if (model_->numberIntegers() !=
     818                model_->numberObjects())
     819            setWhen(0);
     820        // Only works if costs positive, coefficients positive and all rows E or L
     821        // And if values are integer
     822        OsiSolverInterface * solver = model_->solver();
     823        const double * columnLower = solver->getColLower();
     824        const double * rowUpper = solver->getRowUpper();
     825        const double * rowLower = solver->getRowLower();
     826        const double * objective = solver->getObjCoefficients();
     827        double direction = solver->getObjSense();
     828
     829        int numberRows = solver->getNumRows();
     830        // Column copy
     831        const double * element = matrix_.getElements();
     832        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
     833        const int * columnLength = matrix_.getVectorLengths();
     834        bool good = true;
     835        for (int iRow = 0; iRow < numberRows; iRow++) {
     836            if (rowUpper[iRow] > 1.0e30)
     837                good = false;
     838            if (rowLower[iRow] > 0.0 && rowLower[iRow] != rowUpper[iRow])
     839                good = false;
     840            if (floor(rowUpper[iRow] + 0.5) != rowUpper[iRow])
     841                good = false;
     842        }
     843        int numberColumns = solver->getNumCols();
     844        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     845            if (objective[iColumn]*direction < 0.0)
     846                good = false;
     847            if (columnLower[iColumn] < 0.0)
     848                good = false;
     849            CoinBigIndex j;
     850            for (j = columnStart[iColumn];
     851                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     852                if (element[j] < 0.0)
     853                    good = false;
     854                if (floor(element[j] + 0.5) != element[j])
     855                    good = false;
     856            }
     857        }
     858        if (!good)
     859            setWhen(0); // switch off
     860    }
     861}
     862
     863
Note: See TracChangeset for help on using the changeset viewer.