Changeset 1879 for trunk


Ignore:
Timestamp:
Sep 13, 2012 12:31:33 PM (7 years ago)
Author:
forrest
Message:

fix a few problems with Aboca

Location:
trunk/Clp/src
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/AbcPrimalColumnSteepest.cpp

    r1878 r1879  
    585585     }
    586586#endif
    587 #if 1
     587#if 0
    588588     for (int iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    589589       if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
  • trunk/Clp/src/AbcSimplex.cpp

    r1878 r1879  
    40504050  // add values pass options
    40514051  minimumThetaMovement_=1.0e-12;
     4052  if ((specialOptions_&8192)!=0)
     4053    minimumThetaMovement_=1.0e-15;
    40524054  int returnCode=static_cast<AbcSimplexPrimal *> (this)->primal(ifValuesPass);
    40534055  stateOfProblem_ &= ~VALUES_PASS;
     
    40644066  minimumThetaMovement_=1.0e-12;
    40654067  int iPass=0;
    4066   while (problemStatus_==10) {
     4068  while (problemStatus_==10&&minimumThetaMovement_>1.0e-15) {
    40674069    iPass++;
    40684070    if (minimumThetaMovement_==1.0e-12)
     
    40794081      perturbation_=101;
    40804082    baseIteration_=numberIterations_;
    4081     static_cast<AbcSimplexDual *> (this)->dual();
     4083    if ((specialOptions_&8192)==0)
     4084      static_cast<AbcSimplexDual *> (this)->dual();
    40824085    baseIteration_=numberIterations_;
    40834086    if (problemStatus_==10) {
    4084       perturbation_=100;
     4087      if (initialSumInfeasibilities_>=0.0) {
     4088        if (savePerturbation>=100) {
     4089          perturbation_=100;
     4090        } else {
     4091          if (savePerturbation==50)
     4092            perturbation_=CoinMax(51,HEAVY_PERTURBATION-4-iPass); //smaller
     4093          else
     4094            perturbation_=CoinMax(51,savePerturbation-1-iPass); //smaller
     4095        }
     4096      } else {
     4097        specialOptions_ |= 8192; // stop going to dual
     4098        perturbation_=savePerturbation;
     4099        abcFactorization_->setPivots(100000); // force factorization
     4100        initialSumInfeasibilities_=1.23456789e-5;
     4101        // clean pivots
     4102        int numberBasic=0;
     4103        int * pivot=pivotVariable();
     4104        for (int i=0;i<numberTotal_;i++) {
     4105          if (getInternalStatus(i)==basic)
     4106            pivot[numberBasic++]=i;
     4107        }
     4108        assert (numberBasic==numberRows_);
     4109      }
     4110      if (iPass>2)
     4111        perturbation_=100;
    40854112      copyFromSaved(14);
    40864113      baseIteration_=numberIterations_;
     
    59245951      << numberTimes_
    59255952      << CoinMessageEol;
     5953    printf("loop detected %d times out of %d\n",numberBadTimes_,numberTimes_);
    59265954    numberBadTimes_++;
    59275955    if (numberBadTimes_ < 10) {
     
    59736001          model->setSequenceIn(iSequence);
    59746002          model->setFlagged(iSequence);
     6003          model->setLastBadIteration(model->numberIterations());
    59756004          model->setSequenceIn(save);
    59766005          //printf("flagging %d from loop\n",iSequence);
  • trunk/Clp/src/AbcSimplexDual.cpp

    r1878 r1879  
    208208#define CLEAN_FIXED 0
    209209// Startup part of dual (may be extended to other algorithms)
     210// To force to follow another run put logfile name here and define
     211#define FORCE_FOLLOW
     212//#ifdef FORCE_FOLLOW
     213static FILE * fpFollow = NULL;
     214static const char * forceFile = NULL;
     215static int force_in = -1;
     216static int force_out = -1;
     217static int force_way_in = -1;
     218static int force_way_out = -1;
     219static int force_iterations = 0;
     220int force_argc=0;
     221const char ** force_argv=NULL;
     222#endif
    210223void
    211224AbcSimplexDual::startupSolve()
    212225{
     226#ifdef FORCE_FOLLOW
     227  int ifld;
     228  for (ifld=1;ifld<force_argc;ifld++) {
     229    if (strstr(argv[ifld],".log")) {
     230      forceFile=argv[ifld];
     231      break;
     232    }
     233  }
     234  if (!fpFollow && strlen(forceFile)) {
     235    fpFollow = fopen(forceFile, "r");
     236    assert (fpFollow);
     237  }
     238  if (fpFollow) {
     239    int numberRead= atoi(argv[ifld+1]);
     240    force_iterations=atoi(argv[ifld+1]);
     241    printf("skipping %d iterations and then doing %d\n",numberRead,force_iterations);
     242    for (int iteration=0;iteration<numberRead;iteration++) {
     243      // read to next Clp0102
     244      char temp[300];
     245      while (fgets(temp, 250, fpFollow)) {
     246        if (!strncmp(temp, "dirOut", 6))
     247          break;
     248      }
     249      sscanf(temp+7 , "%d%*s%d%*s%*s%d",
     250             &force_way_out, &force_way_in);
     251      fgets(temp, 250, fpFollow);
     252      char cin, cout;
     253      int whichIteration;
     254      sscanf(temp , "%d%*f%*s%*c%c%d%*s%*c%c%d",
     255             &whichIteration, &cin, &force_in, &cout, &force_out);
     256      assert (whichIterations==iteration+1);
     257      if (cin == 'C')
     258        force_in += numberColumns_;
     259      if (cout == 'C')
     260        force_out += numberColumns_;
     261      printf("Iteration %d will force %d out (way %d) and %d in (way %d)\n",
     262             whichIteration, force_out, force_way_out,force_in,force_way_in);
     263      assert (getInternalStatus(force_out)==basic);
     264      assert (getInternalStatus(force_in)!=basic);
     265      setInternalStatus(force_in,basic);
     266      if (force_way_out==-1)
     267        setInternalStatus(force_out,atUpperBound);
     268      else
     269        setInternalStatus(force_out,atLowerBound);
     270    }
     271    // clean up
     272    int numberBasic=0;
     273    for (int i=0;i<numberTotal_;i++) {
     274      if (getInternalStatus(i)==basic) {
     275        abcPivotVariable_[numberBasic++]=i;
     276      } else if (getInternalStatus(i)==atLowerBound) {
     277        if (abcLower_[i]<-1.0e30) {
     278          abcLower_[i]=abcUpper_[i]-dualBound_;
     279          setFakeLower(i);
     280        }
     281      } else if (getInternalStatus(i)==atUpperBound) {
     282        if (abcUpper_[i]>1.0e30) {
     283          abcUpper_[i]=abcLower_[i]+dualBound_;
     284          setFakeUpper(i);
     285        }
     286      }
     287    }
     288    assert (numberBasic==numberRows_);
     289  }
     290#endif
    213291  // initialize - no values pass and algorithm_ is -1
    214292  // put in standard form (and make row copy)
     
    27072785#endif
    27082786          setFlagged(sequenceOut_);
     2787          lastBadIteration_ = numberIterations_; // say be more cautious
    27092788          // so can't happen immediately again
    27102789          sequenceOut_=-1;
     
    27442823  double savePrimalInfeasibilities=sumPrimalInfeasibilities_;
    27452824  if ((moreSpecialOptions_&16384)!=0&&(moreSpecialOptions_&8192)==0) {
    2746     if (!numberIterations_) {
     2825    if (numberIterations_==baseIteration_) {
    27472826      // If primal feasible and looks suited to primal go to primal
    27482827      if (!sumPrimalInfeasibilities_&&sumDualInfeasibilities_>1.0&&
     
    28472926#endif
    28482927          setFlagged(sequenceOut_);
     2928          lastBadIteration_ = numberIterations_; // say be more cautious
    28492929#if ABC_NORMAL_DEBUG>0
    28502930        } else {
     
    54605540    int saveStatus=problemStatus_;
    54615541    // check update
     5542#if MOVE_REPLACE_PART1A < 0
     5543    checkReplacePart1();
     5544#endif
    54625545    int updateStatus=abcFactorization_->checkReplacePart2(pivotRow_,btranAlpha_,alpha_,ftAlpha_);
    54635546#if ABC_DEBUG
     
    56265709#endif
    56275710        setFlagged(sequenceOut_);
     5711        lastBadIteration_ = numberIterations_; // say be more cautious
    56285712        // so can't happen immediately again
    56295713        sequenceOut_=-1;
  • trunk/Clp/src/AbcSimplexParallel.cpp

    r1878 r1879  
    11321132  // get sign for finding row of tableau
    11331133  // create as packed
    1134 #ifdef MOVE_REPLACE_PART1A
     1134#if MOVE_REPLACE_PART1A > 0
    11351135  if (!abcFactorization_->usingFT()) {
    11361136#endif
    11371137    usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
    11381138    abcFactorization_->updateColumnTranspose(usefulArray_[arrayForBtran_]);
    1139 #ifdef MOVE_REPLACE_PART1A
     1139#if MOVE_REPLACE_PART1A > 0
    11401140  } else {
    11411141    cilk_spawn abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
     
    11781178  //cilk
    11791179  getTableauColumnPart1Cilk();
    1180 #ifndef MOVE_REPLACE_PART1A
     1180#if MOVE_REPLACE_PART1A <= 0
    11811181  cilk_spawn getTableauColumnPart2();
     1182#if MOVE_REPLACE_PART1A == 0
    11821183  cilk_spawn checkReplacePart1();
     1184#endif
    11831185  numberFlipped=flipBounds();
    11841186  cilk_sync;
  • trunk/Clp/src/AbcSimplexPrimal.cpp

    r1878 r1879  
    737737  double lastSumInfeasibility = COIN_DBL_MAX;
    738738  int lastNumberInfeasibility = 1;
     739#ifndef CLP_CAUTION
     740#define CLP_CAUTION 1
     741#endif
     742#if CLP_CAUTION
     743  double lastAverageInfeasibility = sumDualInfeasibilities_ /
     744    static_cast<double>(numberDualInfeasibilities_ + 1);
     745#endif
    739746  if (numberIterations_&&type) {
    740747    lastSumInfeasibility = abcNonLinearCost_->sumInfeasibilities();
    741748    lastNumberInfeasibility = abcNonLinearCost_->numberInfeasibilities();
     749  } else {
     750    lastAverageInfeasibility=1.0e10;
    742751  }
    743752  bool ifValuesPass=(stateOfProblem_&VALUES_PASS)!=0;
     
    810819    // put back original costs and then check
    811820    // createRim(4); // costs do not change
    812 #ifndef CLP_CAUTION
    813 #define CLP_CAUTION 1
    814 #endif
    815 #if CLP_CAUTION
    816     double lastAverageInfeasibility = sumDualInfeasibilities_ /
    817       static_cast<double>(numberDualInfeasibilities_ + 10);
    818 #endif
    819821    if (ifValuesPass&&numberIterations_==baseIteration_) {
    820822      abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
     
    12201222  // we may wish to say it is optimal even if infeasible
    12211223  bool alwaysOptimal = (specialOptions_ & 1) != 0;
     1224#if CLP_CAUTION
     1225  // If twice nearly there ...
     1226  if (lastAverageInfeasibility<2.0*dualTolerance_) {
     1227    double averageInfeasibility = sumDualInfeasibilities_ /
     1228      static_cast<double>(numberDualInfeasibilities_ + 1);
     1229    printf("last av %g now %g\n",lastAverageInfeasibility,
     1230           averageInfeasibility);
     1231    if (averageInfeasibility<2.0*dualTolerance_)
     1232      sumOfRelaxedDualInfeasibilities_ = 0.0;
     1233  }
     1234#endif
    12221235  // give code benefit of doubt
    12231236  if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
     
    12701283      //may need infeasiblity cost changed
    12711284      // we can see if we can construct a ray
    1272       // make up a new objective
    1273       double saveWeight = infeasibilityCost_;
    1274       // save nonlinear cost as we are going to switch off costs
    1275       AbcNonLinearCost * nonLinear = abcNonLinearCost_;
    12761285      // do twice to make sure Primal solution has settled
    12771286      // put non-basics to bounds in case tolerance moved
     
    12801289      abcNonLinearCost_->checkInfeasibilities(0.0);
    12811290      gutsOfPrimalSolution(3);
    1282      
    1283       infeasibilityCost_ = 1.0e100;
    12841291      // put back original costs
    12851292      CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);;
    12861293      abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
     1294      gutsOfPrimalSolution(3);
    12871295      // may have fixed infeasibilities - double check
    12881296      if (abcNonLinearCost_->numberInfeasibilities() == 0) {
    12891297        // carry on
    12901298        problemStatus_ = -1;
    1291         infeasibilityCost_ = saveWeight;
    12921299        abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
    12931300      } else {
    1294         abcNonLinearCost_ = NULL;
    1295         // scale
    1296         int i;
    1297         for (i = 0; i < numberRows_ + numberColumns_; i++)
    1298           abcCost_[i] *= 1.0e-95;
    1299         gutsOfPrimalSolution(3);
    1300         abcNonLinearCost_ = nonLinear;
    1301         infeasibilityCost_ = saveWeight;
    13021301        if ((infeasibilityCost_ >= 1.0e18 ||
    1303              numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
     1302             numberDualInfeasibilities_ == 0) && perturbation_ == 101
     1303            && (specialOptions_&8192)==0) {
    13041304          goToDual = unPerturb(); // stop any further perturbation
    13051305#ifndef TRY_ABC_GUS
     
    13071307            goToDual = false;
    13081308#endif
    1309           abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
    13101309          numberDualInfeasibilities_ = 1; // carry on
    13111310          problemStatus_ = -1;
    13121311        } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e-2
    13131312#ifndef TRY_ABC_GUS
    1314                    &&(moreSpecialOptions_ & 256) == 0
     1313                   &&((moreSpecialOptions_ & 256) == 0&&(specialOptions_ & 8192) == 0)
     1314#else
     1315                   &&(specialOptions_ & 8192) == 0
    13151316#endif
    13161317                   ) {
    13171318          goToDual = true;
    1318           abcFactorization_->pivotTolerance(CoinMax(0.9, abcFactorization_->pivotTolerance()));
     1319          abcFactorization_->pivotTolerance(CoinMax(0.5, abcFactorization_->pivotTolerance()));
    13191320        }
    13201321        if (!goToDual) {
     
    13471348            // put back original costs and then check
    13481349            CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);;
    1349 #ifndef TRY_ABC_GUS
    13501350            abcNonLinearCost_->checkInfeasibilities(0.0);
    13511351            gutsOfPrimalSolution(3);
    13521352            problemStatus_ = -1; //continue
     1353#ifndef TRY_ABC_GUS
    13531354            goToDual = false;
    13541355#else
    1355             goToDual=true;
     1356            if((specialOptions_&8192)==0&&!sumOfRelaxedDualInfeasibilities_)
     1357              goToDual=true;
    13561358#endif
    13571359          } else {
     
    13671369#ifndef TRY_ABC_GUS
    13681370        if ((numberRows_ > 20000 || numberDualInfeasibilities_) && !numberTimesOptimal_)
     1371#else
     1372          if ((specialOptions_&8192)!=0)
     1373#endif
    13691374          goToDual = false; // Better to carry on a bit longer
    1370 #endif
    13711375        lastCleaned_ = -1; // carry on
    13721376      }
     
    15181522    moreSpecialOptions_ &= ~16; // clear check accuracy flag
    15191523#ifndef TRY_ABC_GUS
    1520   if ((moreSpecialOptions_ & 256) != 0||saveSum)
     1524  if ((moreSpecialOptions_ & 256) != 0||saveSum||(specialOptions_ & 8192) != 0)
     1525    goToDual=false;
     1526#else
     1527  if ((specialOptions_ & 8192) != 0)
    15211528    goToDual=false;
    15221529#endif
     
    33293336          sequenceOut_ = -1;
    33303337          sequenceIn_=-1;
     3338          if (abcFactorization_->pivots()<10&&abcFactorization_->pivotTolerance()<0.25)
     3339            abcFactorization_->saferTolerances(1.0,-1.03);
    33313340          break;
    33323341        } else {
     
    33443353            // Make safer?
    33453354            double tolerance=abcFactorization_->pivotTolerance();
    3346             abcFactorization_->saferTolerances (-0.99, -1.03);
     3355            abcFactorization_->saferTolerances (1.0, -1.03);
    33473356#endif
    33483357            abcProgress_.clearBadTimes();
  • trunk/Clp/src/ClpPresolve.cpp

    r1878 r1879  
    566566        // skip if no cost (should be able to get rid of)
    567567        if (!cost[icol]) {
    568           printf("should be able to get rid of %d with no cost\n",icol);
     568          PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n",icol));
    569569          continue;
    570570        }
    571571        // skip if negative cost for now
    572572        if (cost[icol]<0.0) {
    573           printf("code for negative cost\n");
     573          PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
    574574          continue;
    575575        }
     
    655655        //#define PRINT_VALUES
    656656        if (!binding0||!binding1) {
    657           printf("Row redundant for column %d\n",icol);
     657          PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n",icol));
    658658        } else {
    659659#ifdef PRINT_VALUES
  • trunk/Clp/src/ClpSimplex.hpp

    r1878 r1879  
    12711271     inline int lastBadIteration() const {
    12721272          return lastBadIteration_;
     1273     }
     1274     /// Set so we know when to be cautious
     1275     inline void setLastBadIteration(int value) {
     1276          lastBadIteration_=value;
    12731277     }
    12741278     /// Progress flag - at present 0 bit says artificials out
  • trunk/Clp/src/ClpSolve.cpp

    r1878 r1879  
    548548      this->dual(0);
    549549    else
    550       this->primal(startUp);
     550      this->primal(startUp ? 1 : 0);
    551551  } else {
    552552    AbcSimplex * abcModel2=new AbcSimplex(*this);
     
    560560    int crashState=abcModel2->abcState()&(256+512+1024);
    561561    abcModel2->setAbcState(CLP_ABC_WANTED|crashState);
    562     int ifValuesPass=startUp;
     562    int ifValuesPass=startUp ? 1 : 0;
    563563    // temp
    564564    if (fabs(this->primalTolerance()-1.001e-6)<0.999e-9) {
     
    672672#endif
    673673#endif
    674     if (!solveType)
     674    if (!solveType) {
    675675      abcModel2->ClpSimplex::doAbcDual();
    676     else
     676    } else {
     677      int saveOptions=abcModel2->specialOptions();
     678      if (startUp==2)
     679        abcModel2->setSpecialOptions(8192|saveOptions);
    677680      abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
     681      abcModel2->setSpecialOptions(saveOptions);
     682    }
    678683#if ABC_INSTRUMENT
    679684    double timeCpu=CoinCpuTime()-startTimeCpu;
     
    31163121                        primal(1);
    31173122#else
    3118                         dealWithAbc(1,1,interrupt);
     3123                        dealWithAbc(1,2,interrupt);
    31193124#endif
    31203125                    }
     
    31343139          } else if (rcode >= 0) {
    31353140#ifdef ABC_INHERIT
    3136             dealWithAbc(1,1,true);
     3141            dealWithAbc(1,2,true);
    31373142#else
    31383143            primal(1);
Note: See TracChangeset for help on using the changeset viewer.