Ignore:
Timestamp:
Dec 30, 2010 12:45:15 PM (10 years ago)
Author:
forrest
Message:

add some heuristic variants

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcHeuristicRENS.cpp

    r1499 r1564  
    1414#include "CbcHeuristicRENS.hpp"
    1515#include "CoinWarmStartBasis.hpp"
     16#include "CoinSort.hpp"
    1617#include "CbcBranchActual.hpp"
    1718#include "CbcStrategy.hpp"
     
    8990
    9091    const double * currentSolution = solver->getColSolution();
    91     const double * dj = solver->getReducedCost();
     92    int numberColumns = solver->getNumCols();
     93    double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
    9294    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
    93     int numberColumns = newSolver->getNumCols();
    9495    double direction = newSolver->getObjSense();
    9596    int type = rensType_&15;
     
    9798    const double * colLower = newSolver->getColLower();
    9899    const double * colUpper = newSolver->getColUpper();
    99     if ((type&3)==3) {
    100         double total=0.0;
    101         int n=0;
    102         CoinWarmStartBasis * basis =
    103           dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
    104         if (basis) {
    105             for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    106               if (colUpper[iColumn]>colLower[iColumn]&&
    107                   basis->getStructStatus(iColumn) !=
    108                   CoinWarmStartBasis::basic) {
    109                   n++;
    110                   total += fabs(dj[iColumn]);
     100    int numberFixed = 0;
     101    if ((type&7)==3) {
     102      double total=0.0;
     103      int n=0;
     104      CoinWarmStartBasis * basis =
     105        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     106      if (basis&&basis->getNumArtificial()) {
     107        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     108          if (colUpper[iColumn]>colLower[iColumn]&&
     109              basis->getStructStatus(iColumn) !=
     110              CoinWarmStartBasis::basic) {
     111            n++;
     112            total += fabs(dj[iColumn]);
     113          }
     114        }
     115        if (n)
     116          djTolerance = (0.01*total)/static_cast<double>(n);
     117        delete basis;
     118      }
     119    } else if ((type&7)==4) {
     120      double * sort = new double [numberColumns];
     121      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     122        sort[iColumn]=1.0e30;
     123        if (colUpper[iColumn]>colLower[iColumn]) {
     124          sort[iColumn] = fabs(dj[iColumn]);
     125        }
     126      }
     127      std::sort(sort,sort+numberColumns);
     128      int last = static_cast<int>(numberColumns*fractionSmall_);
     129      djTolerance = sort[last];
     130      delete [] sort;
     131    } else if ((type&7)==5||(type&7)==6) {
     132      // SOS type fixing
     133      bool fixSets = (type&7)==5;
     134      CoinWarmStartBasis * basis =
     135        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     136      if (basis&&basis->getNumArtificial()) {
     137        //const double * rowLower = solver->getRowLower();
     138        const double * rowUpper = solver->getRowUpper();
     139       
     140        int numberRows = solver->getNumRows();
     141        // Column copy
     142        const CoinPackedMatrix * matrix = solver->getMatrixByCol();
     143        const double * element = matrix->getElements();
     144        const int * row = matrix->getIndices();
     145        const CoinBigIndex * columnStart = matrix->getVectorStarts();
     146        const int * columnLength = matrix->getVectorLengths();
     147        double * bestDj = new double [numberRows];
     148        for (int i=0;i<numberRows;i++) {
     149          if (rowUpper[i]==1.0)
     150            bestDj[i]=1.0e20;
     151          else
     152            bestDj[i]=1.0e30;
     153        }
     154        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     155          if (colUpper[iColumn]>colLower[iColumn]) {
     156            CoinBigIndex j;
     157            if (currentSolution[iColumn]>1.0e-6&&
     158                currentSolution[iColumn]<0.999999) {
     159              for (j = columnStart[iColumn];
     160                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     161                int iRow = row[j];
     162                bestDj[iRow]=1.0e30;
     163              }
     164            } else if ( basis->getStructStatus(iColumn) !=
     165              CoinWarmStartBasis::basic) {
     166              for (j = columnStart[iColumn];
     167                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     168                int iRow = row[j];
     169                if (bestDj[iRow]<1.0e30) {
     170                  if (element[j] != 1.0)
     171                    bestDj[iRow]=1.0e30;
     172                  else
     173                    bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
     174                }
    111175              }
    112176            }
    113             if (n)
    114                 djTolerance = (0.01*total)/static_cast<double>(n);
    115             delete basis;
    116         }
     177          }
     178        }
     179        int nSOS=0;
     180        double * sort = new double [numberRows];
     181        for (int i=0;i<numberRows;i++) {
     182          if (bestDj[i]<1.0e30) {
     183            sort[nSOS++]=bestDj[i];
     184          }
     185        }
     186        if (10*nSOS>8*numberRows) {
     187          std::sort(sort,sort+nSOS);
     188          int last = static_cast<int>(nSOS*0.9*fractionSmall_);
     189          double tolerance = sort[last];
     190          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     191            if (colUpper[iColumn]>colLower[iColumn]) {
     192              CoinBigIndex j;
     193              if (currentSolution[iColumn]<=1.0e-6||
     194                  currentSolution[iColumn]>=0.999999) {
     195                if (fixSets) {
     196                  for (j = columnStart[iColumn];
     197                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     198                    int iRow = row[j];
     199                    if (bestDj[iRow]<1.0e30&&bestDj[iRow]>=tolerance) {
     200                      numberFixed++;
     201                      if (currentSolution[iColumn]<=1.0e-6)
     202                        newSolver->setColUpper(iColumn,0.0);
     203                      else if (currentSolution[iColumn]>=0.999999)
     204                        newSolver->setColLower(iColumn,1.0);
     205                    }
     206                  }
     207                } else if (columnLength[iColumn]==1) {
     208                  // leave more slacks
     209                  int iRow = row[columnStart[iColumn]];
     210                  if (bestDj[iRow]<1.0e30) {
     211                    // fake dj
     212                    dj[iColumn] *= 0.000001;
     213                  }
     214                }
     215              }
     216            }
     217          }
     218          if (fixSets)
     219            djTolerance = 1.0e30;
     220          else
     221            djTolerance = tolerance;
     222        }
     223        delete basis;
     224        delete [] sort;
     225        delete [] bestDj;
     226      }
    117227    }
    118 
     228   
    119229    double primalTolerance;
    120230    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
    121231
    122232    int i;
    123     int numberFixed = 0;
    124233    int numberTightened = 0;
    125234    int numberAtBound = 0;
     
    183292#endif
    184293    }
     294    delete [] dj;
    185295    if (numberFixed > numberIntegers / 5) {
    186296        if ( numberFixed < numberColumns / 5) {
Note: See TracChangeset for help on using the changeset viewer.