Changeset 2404 for trunk/Cbc


Ignore:
Timestamp:
Oct 28, 2018 11:57:03 AM (12 months ago)
Author:
forrest
Message:

try and improve nauty

Location:
trunk/Cbc/src
Files:
3 edited

Legend:

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

    r2382 r2404  
    1963019630        mipStart_.push_back( std::pair<std::string, double>( std::string(colNames[i]), colValues[i]) );
    1963119631}
     19632#ifdef COIN_HAS_NTY
     19633// get rid of all
     19634void
     19635CbcModel::zapSymmetry() {
     19636  delete symmetryInfo_;
     19637  symmetryInfo_=NULL;
     19638}
     19639#endif
    1963219640/* Add SOS info to solver -
    1963319641   Overwrites SOS information in solver with information
  • trunk/Cbc/src/CbcModel.hpp

    r2382 r2404  
    23542354        maximumNumberIterations_ = value;
    23552355    }
     2356#ifdef COIN_HAS_NTY
    23562357    /// Symmetry information
    23572358    inline CbcSymmetry * symmetryInfo() const
    23582359    { return symmetryInfo_;} 
     2360    /// get rid of all
     2361    void zapSymmetry();
     2362#endif
    23592363    /// Set depth for fast nodes
    23602364    inline void setFastNodeDepth(int value) {
  • trunk/Cbc/src/CbcSolver.cpp

    r2402 r2404  
    824824                        std::string & check);
    825825static void generateCode(CbcModel * model, const char * fileName, int type, int preProcess);
     826#ifdef COIN_HAS_NTY
     827// returns number of constraints added
     828static int nautiedConstraints (CbcModel & model, int maxPass);
     829#endif
    826830
    827831// dummy fake main programs for UserClp and UserCbc
     
    72427246                                babModel_->setMoreSpecialOptions2(parameters_[whichParam(CBC_PARAM_INT_MOREMOREMIPOPTIONS, numberParameters_, parameters_)].intValue());
    72437247#ifdef COIN_HAS_NTY
     7248                                int nautyAdded=0;
    72447249                                {
    72457250                                  int jParam = whichParam(CBC_PARAM_STR_ORBITAL,
     
    72477252                                  if(parameters_[jParam].currentOptionAsInteger()) {
    72487253                                    int k = parameters_[jParam].currentOptionAsInteger();
    7249                                     babModel_->setMoreSpecialOptions2(babModel_->moreSpecialOptions2() | (k*128));
     7254                                    if (k<4) {
     7255                                      babModel_->setMoreSpecialOptions2(babModel_->moreSpecialOptions2() | (k*128));
     7256                                    } else {
     7257#define MAX_NAUTY_PASS 2000
     7258                                      nautyAdded=nautiedConstraints(*babModel_,
     7259                                                                    MAX_NAUTY_PASS);
     7260                                    }
    72507261                                  }
    72517262                                }
     
    72567267                                }
    72577268                                babModel_->branchAndBound(statistics);
     7269#ifdef COIN_HAS_NTY
     7270                                if (nautyAdded) {
     7271                                  int * which = new int [nautyAdded];
     7272                                  int numberOldRows=
     7273                                    babModel_->solver()->getNumRows()-nautyAdded;
     7274                                  for (int i=0;i<nautyAdded;i++)
     7275                                    which[i]=i+numberOldRows;
     7276                                  babModel_->solver()->deleteRows(nautyAdded,which);
     7277                                  delete [] which;
     7278                                  babModel_->solver()->resolve();
     7279                                }
     7280#endif
    72587281                                if (truncateColumns<babModel_->solver()->getNumCols()) {
    72597282                                  OsiSolverInterface * solverX = babModel_->solver();
     
    96629685                                    statistics_nprocessedrows, statistics_nprocessedcols);
    96639686                            for (int i = 0; i < statistics_number_generators; i++)
    9664                                 fprintf(fp, ",%d", statistics_number_cuts[i]!=NULL ? statistics_number_cuts[i] : 0);
     9687                                fprintf(fp, ",%d", statistics_number_cuts!=NULL ? statistics_number_cuts[i] : 0);
    96659688                            fprintf(fp, ",");
    96669689                            for (int i = 1; i < argc; i++) {
     
    1241212435#endif
    1241312436}
     12437#ifdef COIN_HAS_NTY
     12438#include "CbcSymmetry.hpp"
     12439// returns number of constraints added
     12440static int nautiedConstraints (CbcModel & model,int maxPass)
     12441{
     12442  bool changed =true;
     12443  int numberAdded=0;
     12444  int numberPasses=0;
     12445  int changeType = 0; //(more2&(128|256))>>7;
     12446  OsiSolverInterface * solverOriginal = model.solver();
     12447#define REALLY_CHANGE
     12448#ifdef REALLY_CHANGE
     12449  OsiSolverInterface * solver = solverOriginal;
     12450#else
     12451  int numberOriginalRows=solverOriginal->getNumRows();
     12452  OsiSolverInterface * solver = solverOriginal->clone();
     12453#endif
     12454  while (changed) {
     12455    changed=false;
     12456    CbcSymmetry symmetryInfo;
     12457    //symmetryInfo.setModel(&model);
     12458    // for now strong is just on counts - use user option
     12459    //int maxN=5000000;
     12460    //OsiSolverInterface * solver = model.solver();
     12461    symmetryInfo.setupSymmetry(*solver);
     12462    int numberGenerators = symmetryInfo.statsOrbits(&model,0);
     12463    if (numberGenerators) {
     12464      //symmetryInfo.Print_Orbits();
     12465      int numberUsefulOrbits = symmetryInfo.numberUsefulOrbits();
     12466      if (numberUsefulOrbits) {
     12467        symmetryInfo.Compute_Symmetry();
     12468        symmetryInfo.fillOrbits(/*true*/);
     12469        const int * orbits = symmetryInfo.whichOrbit();
     12470        int numberUsefulOrbits = symmetryInfo.numberUsefulOrbits();
     12471        int * counts = new int[numberUsefulOrbits];
     12472        memset(counts,0,numberUsefulOrbits*sizeof(int));
     12473        int numberColumns = solver->getNumCols();
     12474        int numberUseful=0;
     12475        if (changeType==1) {
     12476          // just 0-1
     12477          for (int i=0;i<numberColumns;i++) {
     12478            int iOrbit=orbits[i];
     12479            if (iOrbit>=0) {
     12480              if (solver->isBinary(i)) {
     12481                counts[iOrbit]++;
     12482                numberUseful++;
     12483              }
     12484            }
     12485          }
     12486        } else if (changeType==2) {
     12487          // just integer
     12488          for (int i=0;i<numberColumns;i++) {
     12489            int iOrbit=orbits[i];
     12490            if (iOrbit>=0) {
     12491              if (solver->isInteger(i)) {
     12492                counts[iOrbit]++;
     12493                numberUseful++;
     12494              }
     12495            }
     12496          }
     12497        } else {
     12498          // all
     12499          for (int i=0;i<numberColumns;i++) {
     12500            int iOrbit=orbits[i];
     12501            if (iOrbit>=0) {
     12502              counts[iOrbit]++;
     12503              numberUseful++;
     12504            }
     12505          }
     12506        }
     12507        int iOrbit=-1;
     12508#define LONGEST 1
     12509#if LONGEST
     12510        // choose longest
     12511        int maxOrbit=0;
     12512        for (int i=0;i<numberUsefulOrbits;i++) {
     12513          if (counts[i]>maxOrbit) {
     12514            maxOrbit=counts[i];
     12515            iOrbit=i;
     12516          }
     12517        }
     12518#else
     12519        // choose closest to 2
     12520        int minOrbit=numberColumns+1;
     12521        for (int i=0;i<numberUsefulOrbits;i++) {
     12522          if (counts[i]>1&&counts[i]<minOrbit) {
     12523            minOrbit=counts[i];
     12524            iOrbit=i;
     12525          }
     12526        }
     12527#endif
     12528        delete [] counts;
     12529        if (!numberUseful)
     12530          break;
     12531        // take largest
     12532        const double * solution = solver->getColSolution();
     12533        double * size = new double [numberColumns];
     12534        int * which = new int [numberColumns];
     12535        int nIn=0;
     12536        for (int i=0;i<numberColumns;i++) {
     12537          if (orbits[i]==iOrbit) {
     12538            size[nIn]=-solution[i];
     12539            which[nIn++]=i;
     12540          }
     12541        }
     12542        if (nIn>1) {
     12543          //printf("Using orbit length %d\n",nIn);
     12544          CoinSort_2(size,size+nIn,which);
     12545          size[0]=1.0;
     12546          size[1]=-1.0;
     12547#if LONGEST == 0
     12548          solver->addRow(2,which,size,0.0,COIN_DBL_MAX);
     12549          numberAdded++;
     12550#elif LONGEST == 1
     12551          for (int i=0;i<nIn-1;i++) {
     12552            solver->addRow(2,which+i,size,0.0,COIN_DBL_MAX);
     12553            numberAdded++;
     12554          }
     12555#else
     12556          for (int i=0;i<nIn-1;i++) {
     12557            solver->addRow(2,which,size,0.0,COIN_DBL_MAX);
     12558            which[1]=which[2+i];
     12559            numberAdded++;
     12560          }
     12561#endif
     12562          numberPasses++;
     12563          if (numberPasses<3/*maxPass*/)
     12564            changed=true;
     12565        }
     12566        delete [] size;
     12567        delete [] which;
     12568      }
     12569    }
     12570  }
     12571  if (numberAdded) {
     12572    char general[100];
     12573    if (numberPasses<maxPass)
     12574      sprintf(general,"%d constraints added in %d passes",numberAdded,
     12575             numberPasses);
     12576    else
     12577      sprintf(general,"%d constraints added in %d passes (maximum) - must be better way",numberAdded,
     12578             numberPasses);
     12579    model.messageHandler()->message(CBC_GENERAL,
     12580                                    model.messages())
     12581      << general << CoinMessageEol ;
     12582    //printf("saving model on nauty\n");
     12583    //solver->writeMps("nauty");
     12584  }
     12585#ifndef REALLY_CHANGE
     12586  CbcRowCuts * globalCuts = model.globalCuts();
     12587  int numberRows = solver->getNumRows();
     12588  if (numberRows>numberOriginalRows) {
     12589    const CoinPackedMatrix * rowCopy = solver->getMatrixByRow();
     12590    const int * column = rowCopy->getIndices();
     12591    const int * rowLength = rowCopy->getVectorLengths();
     12592    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     12593    const double * elements = rowCopy->getElements();
     12594    const double * rowLower = solver->getRowLower();
     12595    const double * rowUpper = solver->getRowUpper();
     12596    for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
     12597      OsiRowCut rc;
     12598      rc.setLb(rowLower[iRow]);
     12599      rc.setUb(rowUpper[iRow]);
     12600      CoinBigIndex start = rowStart[iRow];
     12601      rc.setRow(rowLength[iRow],column+start,elements+start,false);
     12602      globalCuts->addCutIfNotDuplicate(rc);
     12603    }
     12604  }
     12605  delete solver;
     12606#endif
     12607  return numberAdded;
     12608}
     12609#endif
    1241412610/*
    1241512611  Version 1.00.00 November 16 2005.
Note: See TracChangeset for help on using the changeset viewer.