Changeset 1656 for trunk

Timestamp:

May 31, 2011 4:10:30 AM (10 years ago)

Author:

forrest

Message:

allow end cuts and lagomory

Location:

trunk/Cbc/src

Files:

• CbcCutGenerator.cpp (modified) (2 diffs)
• CbcCutGenerator.hpp (modified) (2 diffs)
• CbcHeuristicRENS.cpp (modified) (2 diffs)
• CbcModel.cpp (modified) (10 diffs)
• CbcModel.hpp (modified) (6 diffs)
• CbcSolver.cpp (modified) (7 diffs)

trunk/Cbc/src/CbcCutGenerator.cpp

@@ -276 +276 @@
         if (fullScan < 0)
             info.options |= 128;
+        if (whetherInMustCallAgainMode())
+          info.options |= 1024;
         // See if we want alternate set of cuts
         if ((model_->moreSpecialOptions()&16384) != 0)
@@ -539 +541 @@
             }
             //if (!solver->basisIsAvailable())
-            //returnCode=true;
+            //returnCode=true;
+            if (!returnCode) {
+              // bounds changed but still optimal
+#ifdef COIN_HAS_CLP
+              OsiClpSolverInterface * clpSolver
+                = dynamic_cast<OsiClpSolverInterface *> (solver);
+              if (clpSolver) {
+                clpSolver->setLastAlgorithm(2);
+              }
+#endif
+            }
 #ifdef JJF_ZERO
             // Pass across info to pseudocosts

trunk/Cbc/src/CbcCutGenerator.hpp

@@ -167 +167 @@
     }
 
+    /// Get switches (for debug)
+    inline int switches() const {
+        return switches_;
+    }
     /// Get whether the cut generator should be called in the normal place
     inline bool normal() const {
@@ -319 +323 @@
         switches_ &= ~1024;
         switches_ |= yesNo ? 1024 : 0;
+    }
+    /// Whether in must call again mode (or after others)
+    inline bool whetherInMustCallAgainMode() const {
+        return (switches_&2048) != 0;
+    }
+    /// Set whether in must call again mode (or after others)
+    inline void setWhetherInMustCallAgainMode(bool yesNo) {
+        switches_ &= ~2048;
+        switches_ |= yesNo ? 2048 : 0;
+    }
+    /// Whether to call at end
+    inline bool whetherCallAtEnd() const {
+        return (switches_&4096) != 0;
+    }
+    /// Set whether to call at end
+    inline void setWhetherCallAtEnd(bool yesNo) {
+        switches_ &= ~4096;
+        switches_ |= yesNo ? 4096 : 0;
     }
     /// Number of cuts generated at root

trunk/Cbc/src/CbcHeuristicRENS.cpp

@@ -852 +852 @@
         returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
                                          model_->getCutoff(), "CbcHeuristicRENS");
-        if (returnCode < 0) {
-            returnCode = 0; // returned on size
+        if (returnCode < 0 || returnCode == 0) {
 #ifdef RENS_FIX_CONTINUOUS
             if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
@@ -905 +904 @@
                 returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
                                                  model_->getCutoff(), "CbcHeuristicRENS");
-#endif
-            }
+            }
+#endif
+            if (returnCode < 0 || returnCode == 0) {
+                // Do passes fixing up those >0.9 and
+                // down those < 0.05
+#define RENS_PASS 3
+              //#define KEEP_GOING
+#ifdef KEEP_GOING
+                double * saveLower = CoinCopyOfArray(colLower,numberColumns);
+                double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
+                bool badPass=false;
+                int nSolved=0;
+#endif
+                for (int iPass=0;iPass<RENS_PASS;iPass++) {
+                  int nFixed=0;
+                  int nFixedAlready=0;
+                  int nFixedContinuous=0;
+                  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+                    if (colUpper[iColumn]>colLower[iColumn]) {
+                      if (newSolver->isInteger(iColumn)) {
+                        double value = currentSolution[iColumn];
+                        double fixTo = floor(value+0.1);
+                        if (fixTo>value || value-fixTo < 0.05) {
+                          // above 0.9 or below 0.05
+                          nFixed++;
+                          newSolver->setColLower(iColumn, fixTo);
+                          newSolver->setColUpper(iColumn, fixTo);
+                        }
+                      }
+                    } else if (newSolver->isInteger(iColumn)) {
+                      nFixedAlready++;
+                    } else {
+                      nFixedContinuous++;
+                    }
+                  }
+#ifdef CLP_INVESTIGATE2
+                  printf("%d more integers fixed (total %d) plus %d continuous\n",
+                         nFixed,nFixed+nFixedAlready,nFixedContinuous);
+#endif
+#ifdef KEEP_GOING
+                  if (nFixed) {
+                    newSolver->resolve();
+                    if (!newSolver->isProvenOptimal()) {
+                      badPass=true;
+                      break;
+                    } else {
+                      nSolved++;
+                      memcpy(saveLower,colLower,numberColumns*sizeof(double));
+                      memcpy(saveUpper,colUpper,numberColumns*sizeof(double));
+                    }
+                  } else {
+                    break;
+                  }
+#else
+                  if (nFixed) {
+                    newSolver->resolve();
+                    if (!newSolver->isProvenOptimal()) {
+                      returnCode=0;
+                      break;
+                    }
+                    returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
+                                                     model_->getCutoff(), "CbcHeuristicRENS");
+                  } else {
+                    returnCode=0;
+                  }
+                  if (returnCode>=0)
+                    break;
+                }
+                if (returnCode < 0) 
+                returnCode = 0; // returned on size
+#endif
+            }
+#ifdef KEEP_GOING
+            if (badPass) {
+              newSolver->setColLower(saveLower);
+              newSolver->setColUpper(saveUpper);
+              newSolver->resolve();
+            }
+            delete [] saveLower;
+            delete [] saveUpper;
+            if (nSolved)
+              returnCode = 
+                smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
+                                    model_->getCutoff(), "CbcHeuristicRENS");
+            else
+              returnCode=0;
+        }
+#endif
         }
         //printf("return code %d",returnCode);

trunk/Cbc/src/CbcModel.cpp

@@ -1981 +1981 @@
 #endif
     bool feasible;
+    numberSolves_ = 0 ;
     // If NLP then we assume already solved outside branchAndbound
     if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
@@ -2223 +2224 @@
         if (clpSolver)
             clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
+#ifdef CLP_RESOLVE
+        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
+          resolveClp(clpSolver,0);
+        }
+#endif
     }
 #endif
@@ -2332 +2338 @@
     */
     numberIterations_ = 0 ;
-    numberSolves_ = 0 ;
     numberNodes_ = 0 ;
     numberNodes2_ = 0 ;
@@ -2561 +2566 @@
         feasible = false;
 
+#ifdef CLP_RESOLVE
+    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
+      OsiClpSolverInterface * clpSolver
+        = dynamic_cast<OsiClpSolverInterface *> (solver_);
+      if (clpSolver)
+        resolveClp(clpSolver,0);
+    }
+#endif
     // make cut generators less aggressive
     for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
@@ -6920 +6933 @@
     // Whether to increase minimum drop
     bool increaseDrop = (moreSpecialOptions_ & 8192) != 0;
+    for (int i = 0; i < numberCutGenerators_; i++) 
+      generator_[i]->setWhetherInMustCallAgainMode(false);
     /*
       Begin cut generation loop. Cuts generated during each iteration are
@@ -6943 +6958 @@
                 if (!generator_[i]->mustCallAgain())
                     generator_[i]->setSwitchedOff(true);
+                else
+                  generator_[i]->setWhetherInMustCallAgainMode(true);
             }
         }
@@ -7530 +7547 @@
             keepGoing=false;
         }
+        if (numberTries <=0 && feasible && !keepGoing && !parentModel_ && !numberNodes_) {
+          for (int i = 0; i < numberCutGenerators_; i++) {
+            if (generator_[i]->whetherCallAtEnd()
+                &&!generator_[i]->whetherInMustCallAgainMode()) {
+              keepGoing= true;
+              break;
+            }
+          }
+        }
     } while (numberTries > 0 || keepGoing) ;
     /*
@@ -8345 +8371 @@
         int numberRowCutsAfter = numberRowCutsBefore;
         int numberColumnCutsAfter = numberColumnCutsBefore;
+        /*printf("GEN %d %s switches %d\n",
+               i,generator_[i]->cutGeneratorName(),
+               generator_[i]->switches());*/
         bool generate = generator_[i]->normal();
         // skip if not optimal and should be (maybe a cut generator has fixed variables)
@@ -8375 +8404 @@
             }
             if (numberRowCutsBefore < numberRowCutsAfter &&
-                    generator_[i]->mustCallAgain() && status >= 0)
+                generator_[i]->mustCallAgain() && status >= 0)
+              /*printf("%s before %d after %d must %c atend %c off %c endmode %c\n",
+                   generator_[i]->cutGeneratorName(),
+                   numberRowCutsBefore,numberRowCutsAfter,
+                   generator_[i]->mustCallAgain() ? 'Y': 'N',
+                   generator_[i]->whetherCallAtEnd() ? 'Y': 'N',
+                   generator_[i]->switchedOff() ? 'Y': 'N',
+                   generator_[i]->whetherInMustCallAgainMode() ? 'Y': 'N');*/
+            if (numberRowCutsBefore < numberRowCutsAfter &&
+                generator_[i]->mustCallAgain() && status >= 0)
                 status = 1 ; // say must go round
             // Check last cut to see if infeasible
@@ -12107 +12145 @@
     return solver->isProvenOptimal() ? 1 : 0;
 }
+#ifdef CLP_RESOLVE
+// Special purpose resolve
+int 
+CbcModel::resolveClp(OsiClpSolverInterface * clpSolver, int type)
+{
+    numberSolves_++;
+    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
+    int save = clpSimplex->specialOptions();
+    clpSimplex->setSpecialOptions(save | 0x11000000); // say is Cbc (and in branch and bound)
+    int save2 = clpSolver->specialOptions();
+    clpSolver->resolve();
+    if (!numberNodes_) {
+      double error = CoinMax(clpSimplex->largestDualError(),
+                             clpSimplex->largestPrimalError());
+      if (error > 1.0e-2 || !clpSolver->isProvenOptimal()) {
+#ifdef CLP_INVESTIGATE
+        printf("Problem was %s largest dual error %g largest primal %g - safer cuts\n",
+               clpSolver->isProvenOptimal() ? "optimal" : "!infeasible",
+               clpSimplex->largestDualError(),
+               clpSimplex->largestPrimalError());
+#endif
+        if (!clpSolver->isProvenOptimal()) {
+          clpSolver->setSpecialOptions(save2 | 2048);
+          clpSimplex->allSlackBasis(true);
+          clpSolver->resolve();
+          if (!clpSolver->isProvenOptimal()) {
+            bool takeHint;
+            OsiHintStrength strength;
+            clpSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
+            clpSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
+            clpSolver->resolve();
+            clpSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
+          }
+        }
+        // make cuts safer
+        for (int iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
+          CglCutGenerator * generator = generator_[iCutGenerator]->generator();
+          CglGomory * cgl1 = dynamic_cast<CglGomory *>(generator);
+          if (cgl1) {
+            cgl1->setLimitAtRoot(cgl1->getLimit());
+          }
+          CglTwomir * cgl2 = dynamic_cast<CglTwomir *>(generator);
+          if (cgl2) {
+            generator_[iCutGenerator]->setHowOften(-100);
+          }
+        }
+      }
+    }
+    clpSolver->setSpecialOptions(save2);
+#ifdef CLP_INVESTIGATE
+    if (clpSimplex->numberIterations() > 1000)
+      printf("node %d took %d iterations\n", numberNodes_, clpSimplex->numberIterations());
+#endif
+    if (type==0 && clpSolver->isProvenOptimal()) {
+      ClpSimplex newModel(*clpSimplex);
+      newModel.primal();
+      int numberColumns = newModel.numberColumns();
+      int numberRows = newModel.numberRows();
+      double * obj = new double [numberColumns];
+      int * which = new int [numberColumns];
+      const double * solution = clpSimplex->primalColumnSolution();
+      double rhs=1.0e-8;
+      int numberObj=0;
+      double integerTolerance = getDblParam(CbcIntegerTolerance) ;
+      double * objective = newModel.objective();
+      for (int i=0;i<numberColumns;i++) {
+        if (objective[i]) {
+          rhs += objective[i]*solution[i];
+          obj[numberObj]=objective[i];
+          which[numberObj++]=i;
+          objective[i]=0.0;
+        }
+      }
+      if (numberObj) {
+        newModel.addRow(numberObj,which,obj,-COIN_DBL_MAX,rhs);
+      }
+      delete [] obj;
+      delete [] which;
+      double * lower = newModel.columnLower();
+      double * upper = newModel.columnUpper();
+      int numberInf = 0;
+      int numberLb = 0;
+      int numberUb = 0;
+      int numberInt = 0;
+      double sumInf=0.0;
+      for (int i=0;i<numberIntegers_;i++) {
+        int iSequence = integerVariable_[i];
+        double value = solution[iSequence];
+        value = CoinMax(value, lower[iSequence]);
+        value = CoinMin(value, upper[iSequence]);
+        double nearest = floor(value + 0.5);
+        if (value<lower[iSequence]+integerTolerance) {
+          objective[iSequence]=1.0;
+          numberLb++;
+        } else if (value>upper[iSequence]-integerTolerance) {
+          objective[iSequence]=-1.0;
+          numberUb++;
+        } else if (fabs(value - nearest) <= integerTolerance) {
+          // fix??
+          lower[iSequence]=nearest;
+          upper[iSequence]=nearest;
+          numberInt++;
+        } else {
+          lower[iSequence]=floor(value);
+          upper[iSequence]=ceil(value);
+          if (value>nearest) {
+            objective[iSequence]=1.0;
+            sumInf += value-nearest;
+          } else {
+            objective[iSequence]=-1.0;
+            sumInf -= value-nearest;
+          }
+          numberInf++;
+        }
+      }
+      printf("XX %d inf (sum %g), %d at lb %d at ub %d other integer\n",
+             numberInf,sumInf,numberLb,numberUb,numberInt);
+      if (numberInf) {
+        newModel.primal(1);
+        if (!newModel.isProvenOptimal()) {
+          printf("not optimal - scaling issue - switch off\n");
+          clpSimplex->setSpecialOptions(save);
+          if (clpSimplex->status()==4)
+            clpSimplex->setProblemStatus(1);
+          return clpSolver->isProvenOptimal() ? 1 : 0;
+        }
+        //newModel.writeMps("bad.mps");
+        //assert (newModel.isProvenOptimal());
+        printf("%d iterations\n",newModel.numberIterations());
+        int numberInf2 = 0;
+        int numberLb2 = 0;
+        int numberUb2 = 0;
+        int numberInt2 = 0;
+        double sumInf2=0.0;
+        const double * solution = newModel.primalColumnSolution();
+        const double * lower = clpSimplex->columnLower();
+        const double * upper = clpSimplex->columnUpper();
+        for (int i=0;i<numberIntegers_;i++) {
+          int iSequence = integerVariable_[i];
+          double value = solution[iSequence];
+          value = CoinMax(value, lower[iSequence]);
+          value = CoinMin(value, upper[iSequence]);
+          double nearest = floor(value + 0.5);
+          if (value<lower[iSequence]+integerTolerance) {
+            numberLb2++;
+          } else if (value>upper[iSequence]-integerTolerance) {
+            numberUb2++;
+          } else if (fabs(value - nearest) <= integerTolerance) {
+            numberInt2++;
+          } else {
+            if (value>nearest) {
+              sumInf2 += value-nearest;
+            } else {
+              sumInf2 -= value-nearest;
+            }
+            numberInf2++;
+          }
+        }
+        printf("XXX %d inf (sum %g), %d at lb %d at ub %d other integer\n",
+               numberInf2,sumInf2,numberLb2,numberUb2,numberInt2);
+        if (sumInf2<sumInf*0.95) {
+          printf("XXXX suminf reduced from %g (%d) to %g (%d)\n",
+                 sumInf,numberInf,sumInf2,numberInf2);
+          if (numberObj) {
+            newModel.deleteRows(1,&numberRows);
+          }
+          memcpy(newModel.objective(),
+                 clpSimplex->objective(),
+                 numberColumns*sizeof(double));
+          memcpy(newModel.columnLower(),
+                 clpSimplex->columnLower(),
+                 numberColumns*sizeof(double));
+          memcpy(newModel.columnUpper(),
+                 clpSimplex->columnUpper(),
+                 numberColumns*sizeof(double));
+          newModel.setClpScaledMatrix(NULL);
+          newModel.primal(1);
+          printf("%d iterations\n",newModel.numberIterations());
+          int numberInf3 = 0;
+          int numberLb3 = 0;
+          int numberUb3 = 0;
+          int numberInt3 = 0;
+          double sumInf3=0.0;
+          const double * solution = newModel.primalColumnSolution();
+          const double * lower = clpSimplex->columnLower();
+          const double * upper = clpSimplex->columnUpper();
+          for (int i=0;i<numberIntegers_;i++) {
+            int iSequence = integerVariable_[i];
+            double value = solution[iSequence];
+            value = CoinMax(value, lower[iSequence]);
+            value = CoinMin(value, upper[iSequence]);
+            double nearest = floor(value + 0.5);
+            if (value<lower[iSequence]+integerTolerance) {
+              numberLb3++;
+            } else if (value>upper[iSequence]-integerTolerance) {
+              numberUb3++;
+            } else if (fabs(value - nearest) <= integerTolerance) {
+              numberInt3++;
+            } else {
+              if (value>nearest) {
+                sumInf3 += value-nearest;
+              } else {
+                sumInf3 -= value-nearest;
+              }
+              numberInf3++;
+            }
+          }
+          printf("XXXXX %d inf (sum %g), %d at lb %d at ub %d other integer\n",
+                 numberInf3,sumInf3,numberLb3,numberUb3,numberInt3);
+          if (sumInf3<sumInf*0.95) {
+            memcpy(clpSimplex->primalColumnSolution(),
+                   newModel.primalColumnSolution(),
+                   numberColumns*sizeof(double));
+            memcpy(clpSimplex->dualColumnSolution(),
+                   newModel.dualColumnSolution(),
+                   numberColumns*sizeof(double));
+            memcpy(clpSimplex->primalRowSolution(),
+                   newModel.primalRowSolution(),
+                   numberRows*sizeof(double));
+            memcpy(clpSimplex->dualRowSolution(),
+                   newModel.dualRowSolution(),
+                   numberRows*sizeof(double));
+            memcpy(clpSimplex->statusArray(),
+                   newModel.statusArray(),
+                   (numberColumns+numberRows)*sizeof(unsigned char));
+            clpSolver->setWarmStart(NULL);
+          }
+        }
+      }
+    }
+    clpSimplex->setSpecialOptions(save);
+    if (clpSimplex->status()==4)
+      clpSimplex->setProblemStatus(1);
+    return clpSolver->isProvenOptimal() ? 1 : 0;
+}
+#endif
 /*!
     \todo It'd be really nice if there were an overload for this method that

trunk/Cbc/src/CbcModel.hpp

@@ -18 +18 @@
 #include "ClpDualRowPivot.hpp"
 
-//class OsiSolverInterface;
 
 class CbcCutGenerator;
@@ -39 +38 @@
 class CglPreProcess;
 # ifdef COIN_HAS_CLP
+class OsiClpSolverInterface;
 class ClpNodeStuff;
 #endif
@@ -1716 +1716 @@
         16 bit (65536) - Original model had integer bounds
         17 bit (131072) - Perturbation switched off
+        18 bit (262144) - donor CbcModel
+        19 bit (524288) - recipient CbcModel
     */
     inline void setSpecialOptions(int value) {
@@ -1737 +1739 @@
         65536 no cuts in preprocessing
         131072 Time limits elapsed
+        18 bit (262144) - Perturb fathom nodes
+        19 bit (524288) - No limit on fathom nodes
+        20 bit (1048576) - Reduce sum of infeasibilities before cuts
+        21 bit (2097152) - Reduce sum of infeasibilities after cuts
     */
     inline void setMoreSpecialOptions(int value) {
@@ -2004 +2010 @@
     /// Encapsulates solver resolve
     int resolve(OsiSolverInterface * solver);
+#ifdef CLP_RESOLVE
+    /// Special purpose resolve
+    int resolveClp(OsiClpSolverInterface * solver, int type);
+#endif
 
     /** Encapsulates choosing a variable -
@@ -2457 +2467 @@
         18 bit (262144) - Perturb fathom nodes
         19 bit (524288) - No limit on fathom nodes
+        20 bit (1048576) - Reduce sum of infeasibilities before cuts
+        21 bit (2097152) - Reduce sum of infeasibilities after cuts
     */
     int moreSpecialOptions_;

trunk/Cbc/src/CbcSolver.cpp

@@ -1130 +1130 @@
     int numberParameters_ = numberParameters;
     CbcModel & model_ = model;
+#ifdef CBC_USE_INITIAL_TIME
+    if (model_.useElapsedTime())
+      model_.setDblParam(CbcModel::CbcStartSeconds, CoinGetTimeOfDay());
+    else
+      model_.setDblParam(CbcModel::CbcStartSeconds, CoinCpuTime());
+#endif
     CbcModel * babModel_ = NULL;
     int returnMode = 1;
@@ -1753 +1759 @@
                 numberGoodCommands++;
                 if (type == CBC_PARAM_ACTION_BAB && goodModel) {
+#ifndef CBC_USE_INITIAL_TIME
                   if (model_.useElapsedTime())
                     model_.setDblParam(CbcModel::CbcStartSeconds, CoinGetTimeOfDay());
                   else
                     model_.setDblParam(CbcModel::CbcStartSeconds, CoinCpuTime());
+#endif
                     // check if any integers
 #ifndef CBC_OTHER_SOLVER
@@ -2044 +2052 @@
                                 pumpChanged = true;
                             else if (parameters_[iParam].type() == CBC_PARAM_INT_EXPERIMENT) {
+                                int addFlags=0;
+                                if (value>=10) {
+                                  addFlags = 1048576*(value/10);
+                                  value = value % 10;
+                                  parameters[whichParam(CBC_PARAM_INT_EXPERIMENT, numberParameters, parameters)].setIntValue(value);
+                                }
                                 if (value >= 1) {
                                     int values[]={24003,280003,792003,24003,24003};
@@ -2053 +2067 @@
                                       parameters_[iParam].setCurrentOption("on");
                                     }
-                                    int extra4 = values[value-1];
+                                    int extra4 = values[value-1]+addFlags;
                                     parameters[whichParam(CBC_PARAM_INT_EXTRA4, numberParameters, parameters)].setIntValue(extra4);
                                     if (!noPrinting_) {
@@ -4002 +4016 @@
                             int bothFlags = CoinMax(CoinMin(experimentFlag, 1), strategyFlag);
                             // add cut generators if wanted
-                            int switches[20];
-                            int accuracyFlag[20];
+                            int switches[30];
+                            int accuracyFlag[30];
+                            char doAtEnd[30];
+                            memset(doAtEnd,0,30);
                             int numberGenerators = 0;
                             int translate[] = { -100, -1, -99, -98, 1, -1098, -999, 1, 1, 1, -1};
@@ -4078 +4094 @@
                                     }
                                 }
-                                int extra3 = parameters_[whichParam(CBC_PARAM_INT_EXTRA3, numberParameters_, parameters_)].intValue();
-                                if (extra3>=100) {
-                                  // replace
-                                  gomoryGen.passInOriginalSolver(babModel_->solver());
-                                  gomoryGen.setGomoryType(2);
-                                  extra3 = -1;
-                                  parameters_[whichParam(CBC_PARAM_INT_EXTRA3, numberParameters_, parameters_)].setIntValue(extra3);
-                                  babModel_->addCutGenerator(&gomoryGen, translate[gomoryAction], "GomoryL");
-                                } else {
+                                int laGomory = parameters_[whichParam(CBC_PARAM_STR_LAGOMORYCUTS, numberParameters_, parameters_)].currentOptionAsInteger();
+                                int gType = translate[gomoryAction];
+                                if (!laGomory) {
+                                  // Normal
                                   babModel_->addCutGenerator(&gomoryGen, translate[gomoryAction], "Gomory");
-                                }
-                                accuracyFlag[numberGenerators] = 3;
-                                switches[numberGenerators++] = 0;
-                                if (extra3>=10) {
-                                  // just root if 10
-                                  int itype=-99;
-                                  if (extra3>=20) {
-                                    extra3-=10;
-                                    itype = translate[gomoryAction];
-                                  }
-                                  gomoryGen.passInOriginalSolver(babModel_->solver());
-                                  babModel_->addCutGenerator(&gomoryGen, itype, "GomoryL2");
                                   accuracyFlag[numberGenerators] = 3;
                                   switches[numberGenerators++] = 0;
-                                  extra3 -= 10;
-                                  parameters_[whichParam(CBC_PARAM_INT_EXTRA3, numberParameters_, parameters_)].setIntValue(extra3);
+                                } else {
+                                  laGomory--;
+                                  int type = (laGomory % 3)+1;
+                                  int when = laGomory/3;
+                                  char atEnd = (when<2) ? 1 : 0;
+                                  int gomoryTypeMajor = 0;
+                                  if (when<3) {
+                                    // normal as well
+                                    babModel_->addCutGenerator(&gomoryGen, gType, "Gomory");
+                                    accuracyFlag[numberGenerators] = 3;
+                                    switches[numberGenerators++] = 0;
+                                    if (when==2)
+                                      gomoryTypeMajor=10;
+                                  } else {
+                                    when--; // so on
+                                    gomoryTypeMajor=20;
+                                  }
+                                  if (!when)
+                                    gType=-99; // root
+                                  gomoryGen.passInOriginalSolver(babModel_->solver());
+                                  if ((type&1) !=0) {
+                                    // clean
+                                    gomoryGen.setGomoryType(gomoryTypeMajor+1);
+                                    babModel_->addCutGenerator(&gomoryGen, gType, "GomoryL1");
+                                    accuracyFlag[numberGenerators] = 3;
+                                    doAtEnd[numberGenerators]=atEnd;
+                                    switches[numberGenerators++] = 0;
+                                  }
+                                  if ((type&2) !=0) {
+                                    // simple
+                                    gomoryGen.setGomoryType(gomoryTypeMajor+2);
+                                    babModel_->addCutGenerator(&gomoryGen, gType, "GomoryL2");
+                                    accuracyFlag[numberGenerators] = 3;
+                                    doAtEnd[numberGenerators]=atEnd;
+                                    switches[numberGenerators++] = 0;
+                                  }
                                 }
                             }
@@ -4190 +4223 @@
                                 //generator->setWhetherToUse(true);
                                 generator->setInaccuracy(accuracyFlag[iGenerator]);
+                                if (doAtEnd[iGenerator]) {
+                                  generator->setWhetherCallAtEnd(true);
+                                  generator->setMustCallAgain(true);
+                                }
                                 generator->setTiming(true);
                                 if (cutDepth >= 0)