/* $Id: CbcCutGenerator.cpp 1883 2013-04-06 13:33:15Z stefan $ */
// Copyright (C) 2003, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif
#include "CbcConfig.h"
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <cfloat>

#ifdef COIN_HAS_CLP
#include "OsiClpSolverInterface.hpp"
#else
#include "OsiSolverInterface.hpp"
#endif
//#define CGL_DEBUG 1
#ifdef CGL_DEBUG
#include "OsiRowCutDebugger.hpp"
#endif
#include "CbcModel.hpp"
#include "CbcMessage.hpp"
#include "CbcCutGenerator.hpp"
#include "CbcBranchDynamic.hpp"
#include "CglProbing.hpp"
#include "CoinTime.hpp"

// Default Constructor
CbcCutGenerator::CbcCutGenerator ()
        : timeInCutGenerator_(0.0),
        model_(NULL),
        generator_(NULL),
        generatorName_(NULL),
        whenCutGenerator_(-1),
        whenCutGeneratorInSub_(-100),
        switchOffIfLessThan_(0),
        depthCutGenerator_(-1),
        depthCutGeneratorInSub_(-1),
        inaccuracy_(0),
        numberTimes_(0),
        numberCuts_(0),
        numberElements_(0),
        numberColumnCuts_(0),
        numberCutsActive_(0),
        numberCutsAtRoot_(0),
        numberActiveCutsAtRoot_(0),
        numberShortCutsAtRoot_(0),
	switches_(1),
	maximumTries_(-1)
{
}
// Normal constructor
CbcCutGenerator::CbcCutGenerator(CbcModel * model, CglCutGenerator * generator,
                                 int howOften, const char * name,
                                 bool normal, bool atSolution,
                                 bool infeasible, int howOftenInSub,
                                 int whatDepth, int whatDepthInSub,
                                 int switchOffIfLessThan)
        :
        timeInCutGenerator_(0.0),
        depthCutGenerator_(whatDepth),
        depthCutGeneratorInSub_(whatDepthInSub),
        inaccuracy_(0),
        numberTimes_(0),
        numberCuts_(0),
        numberElements_(0),
        numberColumnCuts_(0),
        numberCutsActive_(0),
        numberCutsAtRoot_(0),
        numberActiveCutsAtRoot_(0),
        numberShortCutsAtRoot_(0),
        switches_(1),
	maximumTries_(-1)
{
    if (howOften < -1900) {
        setGlobalCuts(true);
        howOften += 2000;
    } else if (howOften < -900) {
        setGlobalCutsAtRoot(true);
        howOften += 1000;
    }
    model_ = model;
    generator_ = generator->clone();
    generator_->refreshSolver(model_->solver());
    setNeedsOptimalBasis(generator_->needsOptimalBasis());
    whenCutGenerator_ = howOften;
    whenCutGeneratorInSub_ = howOftenInSub;
    switchOffIfLessThan_ = switchOffIfLessThan;
    if (name)
        generatorName_ = CoinStrdup(name);
    else
        generatorName_ = CoinStrdup("Unknown");
    setNormal(normal);
    setAtSolution(atSolution);
    setWhenInfeasible(infeasible);
}

// Copy constructor
CbcCutGenerator::CbcCutGenerator ( const CbcCutGenerator & rhs)
{
    model_ = rhs.model_;
    generator_ = rhs.generator_->clone();
    //generator_->refreshSolver(model_->solver());
    whenCutGenerator_ = rhs.whenCutGenerator_;
    whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
    switchOffIfLessThan_ = rhs.switchOffIfLessThan_;
    depthCutGenerator_ = rhs.depthCutGenerator_;
    depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_;
    generatorName_ = CoinStrdup(rhs.generatorName_);
    switches_ = rhs.switches_;
    maximumTries_ = rhs.maximumTries_;
    timeInCutGenerator_ = rhs.timeInCutGenerator_;
    savedCuts_ = rhs.savedCuts_;
    inaccuracy_ = rhs.inaccuracy_;
    numberTimes_ = rhs.numberTimes_;
    numberCuts_ = rhs.numberCuts_;
    numberElements_ = rhs.numberElements_;
    numberColumnCuts_ = rhs.numberColumnCuts_;
    numberCutsActive_ = rhs.numberCutsActive_;
    numberCutsAtRoot_  = rhs.numberCutsAtRoot_;
    numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_;
    numberShortCutsAtRoot_ = rhs.numberShortCutsAtRoot_;
}

// Assignment operator
CbcCutGenerator &
CbcCutGenerator::operator=( const CbcCutGenerator & rhs)
{
    if (this != &rhs) {
        delete generator_;
        free(generatorName_);
        model_ = rhs.model_;
        generator_ = rhs.generator_->clone();
        generator_->refreshSolver(model_->solver());
        whenCutGenerator_ = rhs.whenCutGenerator_;
        whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
        switchOffIfLessThan_ = rhs.switchOffIfLessThan_;
        depthCutGenerator_ = rhs.depthCutGenerator_;
        depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_;
        generatorName_ = CoinStrdup(rhs.generatorName_);
        switches_ = rhs.switches_;
	maximumTries_ = rhs.maximumTries_;
        timeInCutGenerator_ = rhs.timeInCutGenerator_;
        savedCuts_ = rhs.savedCuts_;
        inaccuracy_ = rhs.inaccuracy_;
        numberTimes_ = rhs.numberTimes_;
        numberCuts_ = rhs.numberCuts_;
        numberElements_ = rhs.numberElements_;
        numberColumnCuts_ = rhs.numberColumnCuts_;
        numberCutsActive_ = rhs.numberCutsActive_;
        numberCutsAtRoot_  = rhs.numberCutsAtRoot_;
        numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_;
        numberShortCutsAtRoot_ = rhs.numberShortCutsAtRoot_;
    }
    return *this;
}

// Destructor
CbcCutGenerator::~CbcCutGenerator ()
{
    free(generatorName_);
    delete generator_;
}

/* This is used to refresh any inforamtion.
   It also refreshes the solver in the cut generator
   in case generator wants to do some work
*/
void
CbcCutGenerator::refreshModel(CbcModel * model)
{
    model_ = model;
    generator_->refreshSolver(model_->solver());
}
/* Generate cuts for the model data contained in si.
   The generated cuts are inserted into and returned in the
   collection of cuts cs.
*/
bool
CbcCutGenerator::generateCuts( OsiCuts & cs , int fullScan, OsiSolverInterface * solver, CbcNode * node)
{
    /*
	  Make some decisions about whether we'll generate cuts. First convert
	  whenCutGenerator_ to a set of canonical values for comparison to the node
	  count.

		 0 <	mod 1000000, with a result of 0 forced to 1
	   -99 <= <= 0	convert to 1
	  -100 =	Off, period
	*/
	int depth;
    if (node)
        depth = node->depth();
    else
        depth = 0;
    int howOften = whenCutGenerator_;
    if (dynamic_cast<CglProbing*>(generator_)) {
        if (howOften == -100 && model_->doCutsNow(3)) {
            howOften = 1; // do anyway
        }
    }
    if (howOften == -100)
        return false;
    int pass = model_->getCurrentPassNumber() - 1;
    if (maximumTries_>0) {
      // howOften means what it says
      if ((pass%howOften)!=0||depth)
	return false;
      else
	howOften=1;
    }
    if (howOften > 0)
        howOften = howOften % 1000000;
    else
        howOften = 1;
    if (!howOften)
        howOften = 1;
    bool returnCode = false;
    //OsiSolverInterface * solver = model_->solver();
    // Reset cuts on first pass
    if (!pass)
        savedCuts_ = OsiCuts();
    /*
	  Determine if we should generate cuts based on node count.
	*/
	bool doThis = (model_->getNodeCount() % howOften) == 0;
    /*
	  If the user has provided a depth specification, it will override the node
	  count specification.
	*/
	if (depthCutGenerator_ > 0) {
        doThis = (depth % depthCutGenerator_) == 0;
        if (depth < depthCutGenerator_)
            doThis = true; // and also at top of tree
    }
	/*
	  A few magic numbers ...

	  The distinction between -100 and 100 for howOften is that we can override 100
	  with fullScan. -100 means no cuts, period. As does the magic number -200 for
	  whenCutGeneratorInSub_.
	*/

    // But turn off if 100
    if (howOften == 100)
        doThis = false;
    // Switch off if special setting
    if (whenCutGeneratorInSub_ == -200 && model_->parentModel()) {
        fullScan = 0;
        doThis = false;
    }
    if (fullScan || doThis) {
        CoinThreadRandom * randomNumberGenerator = NULL;
#ifdef COIN_HAS_CLP
        {
            OsiClpSolverInterface * clpSolver
            = dynamic_cast<OsiClpSolverInterface *> (solver);
            if (clpSolver)
                randomNumberGenerator =
                    clpSolver->getModelPtr()->randomNumberGenerator();
        }
#endif
        double time1 = 0.0;
        if (timing())
            time1 = CoinCpuTime();
        //#define CBC_DEBUG
        int numberRowCutsBefore = cs.sizeRowCuts() ;
        int numberColumnCutsBefore = cs.sizeColCuts() ;
#ifdef JJF_ZERO
        int cutsBefore = cs.sizeCuts();
#endif
        CglTreeInfo info;
        info.level = depth;
        info.pass = pass;
        info.formulation_rows = model_->numberRowsAtContinuous();
        info.inTree = node != NULL;
        info.randomNumberGenerator = randomNumberGenerator;
        info.options = (globalCutsAtRoot()) ? 8 : 0;
        if (ineffectualCuts())
            info.options |= 32;
        if (globalCuts())
            info.options |= 16;
        if (fullScan < 0)
            info.options |= 128;
	if (whetherInMustCallAgainMode())
	  info.options |= 1024;
        // See if we want alternate set of cuts
        if ((model_->moreSpecialOptions()&16384) != 0)
            info.options |= 256;
        if (model_->parentModel())
            info.options |= 512;
        // above had &&!model_->parentModel()&&depth<2)
        incrementNumberTimesEntered();
        CglProbing* generator =
            dynamic_cast<CglProbing*>(generator_);
	//if (!depth&&!pass)
	//printf("Cut generator %s when %d\n",generatorName_,whenCutGenerator_);
        if (!generator) {
            // Pass across model information in case it could be useful
            //void * saveData = solver->getApplicationData();
            //solver->setApplicationData(model_);
            generator_->generateCuts(*solver, cs, info);
            //solver->setApplicationData(saveData);
        } else {
            // Probing - return tight column bounds
            CglTreeProbingInfo * info2 = model_->probingInfo();
            bool doCuts = false;
            if (info2 && !depth) {
                info2->options = (globalCutsAtRoot()) ? 8 : 0;
                info2->level = depth;
                info2->pass = pass;
                info2->formulation_rows = model_->numberRowsAtContinuous();
                info2->inTree = node != NULL;
                info2->randomNumberGenerator = randomNumberGenerator;
                generator->generateCutsAndModify(*solver, cs, info2);
                doCuts = true;
            } else if (depth) {
                /* The idea behind this is that probing may work in a different
                   way deep in tree.  So every now and then try various
                   combinations to see what works.
                */
#define TRY_NOW_AND_THEN
#ifdef TRY_NOW_AND_THEN
                if ((numberTimes_ == 200 || (numberTimes_ > 200 && (numberTimes_ % 2000) == 0))
                        && !model_->parentModel() && info.formulation_rows > 200) {
                    /* In tree, every now and then try various combinations
                       maxStack, maxProbe (last 5 digits)
                       123 is special and means CglProbing will try and
                       be intelligent.
                    */
                    int test[] = {
                        100123,
                        199999,
                        200123,
                        299999,
                        500123,
                        599999,
                        1000123,
                        1099999,
                        2000123,
                        2099999
                    };
                    int n = static_cast<int> (sizeof(test) / sizeof(int));
                    int saveStack = generator->getMaxLook();
                    int saveNumber = generator->getMaxProbe();
                    int kr1 = 0;
                    int kc1 = 0;
                    int bestStackTree = -1;
                    int bestNumberTree = -1;
                    for (int i = 0; i < n; i++) {
                        //OsiCuts cs2 = cs;
                        int stack = test[i] / 100000;
                        int number = test[i] - 100000 * stack;
                        generator->setMaxLook(stack);
                        generator->setMaxProbe(number);
                        int numberRowCutsBefore = cs.sizeRowCuts() ;
                        int numberColumnCutsBefore = cs.sizeColCuts() ;
                        generator_->generateCuts(*solver, cs, info);
                        int numberRowCuts = cs.sizeRowCuts() - numberRowCutsBefore ;
                        int numberColumnCuts = cs.sizeColCuts() - numberColumnCutsBefore ;
#ifdef CLP_INVESTIGATE
                        if (numberRowCuts < kr1 || numberColumnCuts < kc1)
                            printf("Odd ");
#endif
                        if (numberRowCuts > kr1 || numberColumnCuts > kc1) {
#ifdef CLP_INVESTIGATE
                            printf("*** ");
#endif
                            kr1 = numberRowCuts;
                            kc1 = numberColumnCuts;
                            bestStackTree = stack;
                            bestNumberTree = number;
                            doCuts = true;
                        }
#ifdef CLP_INVESTIGATE
                        printf("maxStack %d number %d gives %d row cuts and %d column cuts\n",
                               stack, number, numberRowCuts, numberColumnCuts);
#endif
                    }
                    generator->setMaxLook(saveStack);
                    generator->setMaxProbe(saveNumber);
                    if (bestStackTree > 0) {
                        generator->setMaxLook(bestStackTree);
                        generator->setMaxProbe(bestNumberTree);
#ifdef CLP_INVESTIGATE
                        printf("RRNumber %d -> %d, stack %d -> %d\n",
                               saveNumber, bestNumberTree, saveStack, bestStackTree);
#endif
                    } else {
                        // no good
                        generator->setMaxLook(0);
#ifdef CLP_INVESTIGATE
                        printf("RRSwitching off number %d -> %d, stack %d -> %d\n",
                               saveNumber, saveNumber, saveStack, 1);
#endif
                    }
                }
#endif
                if (generator->getMaxLook() > 0 && !doCuts) {
                    generator->generateCutsAndModify(*solver, cs, &info);
                    doCuts = true;
                }
            } else {
                // at root - don't always do
                if (pass < 15 || (pass&1) == 0) {
                    generator->generateCutsAndModify(*solver, cs, &info);
                    doCuts = true;
                }
            }
            if (doCuts && generator->tightLower()) {
                // probing may have tightened bounds - check
                const double * tightLower = generator->tightLower();
                const double * lower = solver->getColLower();
                const double * tightUpper = generator->tightUpper();
                const double * upper = solver->getColUpper();
                const double * solution = solver->getColSolution();
                int j;
                int numberColumns = solver->getNumCols();
                double primalTolerance = 1.0e-8;
                const char * tightenBounds = generator->tightenBounds();
#ifdef CGL_DEBUG
                const OsiRowCutDebugger * debugger = solver->getRowCutDebugger();
                if (debugger && debugger->onOptimalPath(*solver)) {
                    printf("On optimal path CbcCut\n");
                    int nCols = solver->getNumCols();
                    int i;
                    const double * optimal = debugger->optimalSolution();
                    const double * objective = solver->getObjCoefficients();
                    double objval1 = 0.0, objval2 = 0.0;
                    for (i = 0; i < nCols; i++) {
#if CGL_DEBUG>1
                        printf("%d %g %g %g %g\n", i, lower[i], solution[i], upper[i], optimal[i]);
#endif
                        objval1 += solution[i] * objective[i];
                        objval2 += optimal[i] * objective[i];
                        assert(optimal[i] >= lower[i] - 1.0e-5 && optimal[i] <= upper[i] + 1.0e-5);
                        assert(optimal[i] >= tightLower[i] - 1.0e-5 && optimal[i] <= tightUpper[i] + 1.0e-5);
                    }
                    printf("current obj %g, integer %g\n", objval1, objval2);
                }
#endif
                bool feasible = true;
                if ((model_->getThreadMode()&2) == 0) {
                    for (j = 0; j < numberColumns; j++) {
                        if (solver->isInteger(j)) {
                            if (tightUpper[j] < upper[j]) {
                                double nearest = floor(tightUpper[j] + 0.5);
                                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
                                solver->setColUpper(j, nearest);
                                if (nearest < solution[j] - primalTolerance)
                                    returnCode = true;
                            }
                            if (tightLower[j] > lower[j]) {
                                double nearest = floor(tightLower[j] + 0.5);
                                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
                                solver->setColLower(j, nearest);
                                if (nearest > solution[j] + primalTolerance)
                                    returnCode = true;
                            }
                        } else {
                            if (upper[j] > lower[j]) {
                                if (tightUpper[j] == tightLower[j]) {
                                    // fix
                                    //if (tightLower[j]!=lower[j])
                                    solver->setColLower(j, tightLower[j]);
                                    //if (tightUpper[j]!=upper[j])
                                    solver->setColUpper(j, tightUpper[j]);
                                    if (tightLower[j] > solution[j] + primalTolerance ||
                                            tightUpper[j] < solution[j] - primalTolerance)
                                        returnCode = true;
                                } else if (tightenBounds && tightenBounds[j]) {
                                    solver->setColLower(j, CoinMax(tightLower[j], lower[j]));
                                    solver->setColUpper(j, CoinMin(tightUpper[j], upper[j]));
                                    if (tightLower[j] > solution[j] + primalTolerance ||
                                            tightUpper[j] < solution[j] - primalTolerance)
                                        returnCode = true;
                                }
                            }
                        }
                        if (upper[j] < lower[j] - 1.0e-3) {
                            feasible = false;
                            break;
                        }
                    }
                } else {
                    CoinPackedVector lbs;
                    CoinPackedVector ubs;
                    int numberChanged = 0;
                    bool ifCut = false;
                    for (j = 0; j < numberColumns; j++) {
                        if (solver->isInteger(j)) {
                            if (tightUpper[j] < upper[j]) {
                                double nearest = floor(tightUpper[j] + 0.5);
                                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
                                ubs.insert(j, nearest);
                                numberChanged++;
                                if (nearest < solution[j] - primalTolerance)
                                    ifCut = true;
                            }
                            if (tightLower[j] > lower[j]) {
                                double nearest = floor(tightLower[j] + 0.5);
                                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
                                lbs.insert(j, nearest);
                                numberChanged++;
                                if (nearest > solution[j] + primalTolerance)
                                    ifCut = true;
                            }
                        } else {
                            if (upper[j] > lower[j]) {
                                if (tightUpper[j] == tightLower[j]) {
                                    // fix
                                    lbs.insert(j, tightLower[j]);
                                    ubs.insert(j, tightUpper[j]);
                                    if (tightLower[j] > solution[j] + primalTolerance ||
                                            tightUpper[j] < solution[j] - primalTolerance)
                                        ifCut = true;
                                } else if (tightenBounds && tightenBounds[j]) {
                                    lbs.insert(j, CoinMax(tightLower[j], lower[j]));
                                    ubs.insert(j, CoinMin(tightUpper[j], upper[j]));
                                    if (tightLower[j] > solution[j] + primalTolerance ||
                                            tightUpper[j] < solution[j] - primalTolerance)
                                        ifCut = true;
                                }
                            }
                        }
                        if (upper[j] < lower[j] - 1.0e-3) {
                            feasible = false;
                            break;
                        }
                    }
                    if (numberChanged) {
                        OsiColCut cc;
                        cc.setUbs(ubs);
                        cc.setLbs(lbs);
                        if (ifCut) {
                            cc.setEffectiveness(100.0);
                        } else {
                            cc.setEffectiveness(1.0e-5);
                        }
                        cs.insert(cc);
                    }
                }
                if (!feasible) {
                    // not feasible -add infeasible cut
                    OsiRowCut rc;
                    rc.setLb(COIN_DBL_MAX);
                    rc.setUb(0.0);
                    cs.insert(rc);
                }
            }
            //if (!solver->basisIsAvailable())
	    //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
            char * mark = new char[numberColumns];
            memset(mark, 0, numberColumns);
            int nLook = generator->numberThisTime();
            const int * lookedAt = generator->lookedAt();
            const int * fixedDown = generator->fixedDown();
            const int * fixedUp = generator->fixedUp();
            for (j = 0; j < nLook; j++)
                mark[lookedAt[j]] = 1;
            int numberObjects = model_->numberObjects();
            for (int i = 0; i < numberObjects; i++) {
                CbcSimpleIntegerDynamicPseudoCost * obj1 =
                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ;
                if (obj1) {
                    int iColumn = obj1->columnNumber();
                    if (mark[iColumn])
                        obj1->setProbingInformation(fixedDown[iColumn], fixedUp[iColumn]);
                }
            }
            delete [] mark;
#endif
        }
        CbcCutModifier * modifier = model_->cutModifier();
        if (modifier) {
            int numberRowCutsAfter = cs.sizeRowCuts() ;
            int k ;
            int nOdd = 0;
            //const OsiSolverInterface * solver = model_->solver();
            for (k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) {
                OsiRowCut & thisCut = cs.rowCut(k) ;
                int returnCode = modifier->modify(solver, thisCut);
                if (returnCode) {
                    nOdd++;
                    if (returnCode == 3)
                        cs.eraseRowCut(k);
                }
            }
            if (nOdd)
                COIN_DETAIL_PRINT(printf("Cut generator %s produced %d cuts of which %d were modified\n",
					 generatorName_, numberRowCutsAfter - numberRowCutsBefore, nOdd));
        }
        {
            // make all row cuts without test for duplicate
            int numberRowCutsAfter = cs.sizeRowCuts() ;
            int k ;
#ifdef CGL_DEBUG
            const OsiRowCutDebugger * debugger = solver->getRowCutDebugger();
#endif
            for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) {
                OsiRowCut * thisCut = cs.rowCutPtr(k) ;
#ifdef CGL_DEBUG
                if (debugger && debugger->onOptimalPath(*solver))
                    assert(!debugger->invalidCut(*thisCut));
#endif
                thisCut->mutableRow().setTestForDuplicateIndex(false);
            }
        }
        // Add in saved cuts if violated
        if (false && !depth) {
            const double * solution = solver->getColSolution();
            double primalTolerance = 1.0e-7;
            int numberCuts = savedCuts_.sizeRowCuts() ;
            for (int k = numberCuts - 1; k >= 0; k--) {
                const OsiRowCut * thisCut = savedCuts_.rowCutPtr(k) ;
                double sum = 0.0;
                int n = thisCut->row().getNumElements();
                const int * column = thisCut->row().getIndices();
                const double * element = thisCut->row().getElements();
                assert (n);
                for (int i = 0; i < n; i++) {
                    double value = element[i];
                    sum += value * solution[column[i]];
                }
                if (sum > thisCut->ub() + primalTolerance) {
                    sum = sum - thisCut->ub();
                } else if (sum < thisCut->lb() - primalTolerance) {
                    sum = thisCut->lb() - sum;
                } else {
                    sum = 0.0;
                }
                if (sum) {
                    // add to candidates and take out here
                    cs.insert(*thisCut);
                    savedCuts_.eraseRowCut(k);
                }
            }
        }
        if (!atSolution()) {
            int numberRowCutsAfter = cs.sizeRowCuts() ;
            int k ;
            int nEls = 0;
            int nCuts = numberRowCutsAfter - numberRowCutsBefore;
            // Remove NULL cuts!
            int nNull = 0;
            const double * solution = solver->getColSolution();
            bool feasible = true;
            double primalTolerance = 1.0e-7;
            int shortCut = (depth) ? -1 : generator_->maximumLengthOfCutInTree();
            for (k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) {
                const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
                double sum = 0.0;
                if (thisCut->lb() <= thisCut->ub()) {
                    int n = thisCut->row().getNumElements();
                    if (n <= shortCut)
                        numberShortCutsAtRoot_++;
                    const int * column = thisCut->row().getIndices();
                    const double * element = thisCut->row().getElements();
                    if (n <= 0) {
                        // infeasible cut - give up
                        feasible = false;
                        break;
                    }
                    nEls += n;
                    for (int i = 0; i < n; i++) {
                        double value = element[i];
                        sum += value * solution[column[i]];
                    }
                    if (sum > thisCut->ub() + primalTolerance) {
                        sum = sum - thisCut->ub();
                    } else if (sum < thisCut->lb() - primalTolerance) {
                        sum = thisCut->lb() - sum;
                    } else {
                        sum = 0.0;
                        cs.eraseRowCut(k);
                        nNull++;
                    }
                }
            }
            //if (nNull)
            //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_,
            //       nCuts,nEls,nNull);
            numberRowCutsAfter = cs.sizeRowCuts() ;
            nCuts = numberRowCutsAfter - numberRowCutsBefore;
            nEls = 0;
            for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) {
                const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
                int n = thisCut->row().getNumElements();
                nEls += n;
            }
            //printf("%s has %d cuts and %d elements\n",generatorName_,
            //     nCuts,nEls);
            int nElsNow = solver->getMatrixByCol()->getNumElements();
            int numberColumns = solver->getNumCols();
            int numberRows = solver->getNumRows();
            //double averagePerRow = static_cast<double>(nElsNow)/
            //static_cast<double>(numberRows);
            int nAdd;
            int nAdd2;
            int nReasonable;
            if (!model_->parentModel() && depth < 2) {
                if (inaccuracy_ < 3) {
                    nAdd = 10000;
                    if (pass > 0 && numberColumns > -500)
                        nAdd = CoinMin(nAdd, nElsNow + 2 * numberRows);
                } else {
                    nAdd = 10000;
                    if (pass > 0)
                        nAdd = CoinMin(nAdd, nElsNow + 2 * numberRows);
                }
                nAdd2 = 5 * numberColumns;
                nReasonable = CoinMax(nAdd2, nElsNow / 8 + nAdd);
                if (!depth && !pass) {
                    // allow more
                    nAdd += nElsNow / 2;
                    nAdd2 += nElsNow / 2;
                    nReasonable += nElsNow / 2;
                }
                //if (!depth&&ineffectualCuts())
                //nReasonable *= 2;
            } else {
                nAdd = 200;
                nAdd2 = 2 * numberColumns;
                nReasonable = CoinMax(nAdd2, nElsNow / 8 + nAdd);
            }
            //#define UNS_WEIGHT 0.1
#ifdef UNS_WEIGHT
            const double * colLower = solver->getColLower();
            const double * colUpper = solver->getColUpper();
#endif
            if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts && feasible) {
                //printf("need to remove cuts\n");
                // just add most effective
#ifndef JJF_ONE
                int nDelete = nEls - nReasonable;

                nElsNow = nEls;
                double * sort = new double [nCuts];
                int * which = new int [nCuts];
                // For parallel cuts
                double * element2 = new double [numberColumns];
                //#define USE_OBJECTIVE 2
#ifdef USE_OBJECTIVE
                const double *objective = solver->getObjCoefficients() ;
#if USE_OBJECTIVE>1
                double objNorm = 0.0;
                for (int i = 0; i < numberColumns; i++)
                    objNorm += objective[i] * objective[i];
                if (objNorm)
                    objNorm = 1.0 / sqrt(objNorm);
                else
                    objNorm = 1.0;
                objNorm *= 0.01; // downgrade
#endif
#endif
                CoinZeroN(element2, numberColumns);
                for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) {
                    const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
                    double sum = 0.0;
                    if (thisCut->lb() <= thisCut->ub()) {
                        int n = thisCut->row().getNumElements();
                        const int * column = thisCut->row().getIndices();
                        const double * element = thisCut->row().getElements();
                        assert (n);
#ifdef UNS_WEIGHT
                        double normU = 0.0;
                        double norm = 1.0e-3;
                        int nU = 0;
                        for (int i = 0; i < n; i++) {
                            double value = element[i];
                            int iColumn = column[i];
                            double solValue = solution[iColumn];
                            sum += value * solValue;
                            value *= value;
                            norm += value;
                            if (solValue > colLower[iColumn] + 1.0e-6 &&
                                    solValue < colUpper[iColumn] - 1.0e-6) {
                                normU += value;
                                nU++;
                            }
                        }
#ifdef JJF_ZERO
                        int nS = n - nU;
                        if (numberColumns > 20000) {
                            if (nS > 50) {
                                double ratio = 50.0 / nS;
                                normU /= ratio;
                            }
                        }
#endif
                        norm += UNS_WEIGHT * (normU - norm);
#else
                        double norm = 1.0e-3;
#ifdef USE_OBJECTIVE
                        double obj = 0.0;
#endif
                        for (int i = 0; i < n; i++) {
                            int iColumn = column[i];
                            double value = element[i];
                            sum += value * solution[iColumn];
                            norm += value * value;
#ifdef USE_OBJECTIVE
                            obj += value * objective[iColumn];
#endif
                        }
#endif
                        if (sum > thisCut->ub()) {
                            sum = sum - thisCut->ub();
                        } else if (sum < thisCut->lb()) {
                            sum = thisCut->lb() - sum;
                        } else {
                            sum = 0.0;
                        }
#ifdef USE_OBJECTIVE
                        if (sum) {
#if USE_OBJECTIVE==1
                            obj = CoinMax(1.0e-6, fabs(obj));
                            norm = sqrt(obj * norm);
                            //sum += fabs(obj)*invObjNorm;
                            //printf("sum %g norm %g normobj %g invNorm %g mod %g\n",
                            //     sum,norm,obj,invObjNorm,obj*invObjNorm);
                            // normalize
                            sum /= sqrt(norm);
#else
                            // normalize
                            norm = 1.0 / sqrt(norm);
                            sum = (sum + objNorm * obj) * norm;
#endif
                        }
#else
                        // normalize
                        sum /= sqrt(norm);
#endif
                        //sum /= pow(norm,0.3);
                        // adjust for length
                        //sum /= pow(reinterpret_cast<double>(n),0.2);
                        //sum /= sqrt((double) n);
                        // randomize
                        //double randomNumber =
                        //model_->randomNumberGenerator()->randomDouble();
                        //sum *= (0.5+randomNumber);
                    } else {
                        // keep
                        sum = COIN_DBL_MAX;
                    }
                    sort[k-numberRowCutsBefore] = sum;
                    which[k-numberRowCutsBefore] = k;
                }
                CoinSort_2(sort, sort + nCuts, which);
                // Now see which ones are too similar
                int nParallel = 0;
                double testValue = (depth > 1) ? 0.99 : 0.999999;
                for (k = 0; k < nCuts; k++) {
                    int j = which[k];
                    const OsiRowCut * thisCut = cs.rowCutPtr(j) ;
                    if (thisCut->lb() > thisCut->ub())
                        break; // cut is infeasible
                    int n = thisCut->row().getNumElements();
                    const int * column = thisCut->row().getIndices();
                    const double * element = thisCut->row().getElements();
                    assert (n);
                    double norm = 0.0;
                    double lb = thisCut->lb();
                    double ub = thisCut->ub();
                    for (int i = 0; i < n; i++) {
                        double value = element[i];
                        element2[column[i]] = value;
                        norm += value * value;
                    }
                    int kkk = CoinMin(nCuts, k + 5);
                    for (int kk = k + 1; kk < kkk; kk++) {
                        int jj = which[kk];
                        const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ;
                        if (thisCut2->lb() > thisCut2->ub())
                            break; // cut is infeasible
                        int nB = thisCut2->row().getNumElements();
                        const int * columnB = thisCut2->row().getIndices();
                        const double * elementB = thisCut2->row().getElements();
                        assert (nB);
                        double normB = 0.0;
                        double product = 0.0;
                        for (int i = 0; i < nB; i++) {
                            double value = elementB[i];
                            normB += value * value;
                            product += value * element2[columnB[i]];
                        }
                        if (product > 0.0 && product*product > testValue*norm*normB) {
                            bool parallel = true;
                            double lbB = thisCut2->lb();
                            double ubB = thisCut2->ub();
                            if ((lb < -1.0e20 && lbB > -1.0e20) ||
                                    (lbB < -1.0e20 && lb > -1.0e20))
                                parallel = false;
                            double tolerance;
                            tolerance = CoinMax(fabs(lb), fabs(lbB)) + 1.0e-6;
                            if (fabs(lb - lbB) > tolerance)
                                parallel = false;
                            if ((ub > 1.0e20 && ubB < 1.0e20) ||
                                    (ubB > 1.0e20 && ub < 1.0e20))
                                parallel = false;
                            tolerance = CoinMax(fabs(ub), fabs(ubB)) + 1.0e-6;
                            if (fabs(ub - ubB) > tolerance)
                                parallel = false;
                            if (parallel) {
                                nParallel++;
                                sort[k] = 0.0;
                                break;
                            }
                        }
                    }
                    for (int i = 0; i < n; i++) {
                        element2[column[i]] = 0.0;
                    }
                }
                delete [] element2;
                CoinSort_2(sort, sort + nCuts, which);
                k = 0;
                while (nDelete > 0 || !sort[k]) {
                    int iCut = which[k];
                    const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ;
                    int n = thisCut->row().getNumElements();
                    // may be best, just to save if short
                    if (false && n && sort[k]) {
                        // add to saved cuts
                        savedCuts_.insert(*thisCut);
                    }
                    nDelete -= n;
                    k++;
                    if (k >= nCuts)
                        break;
                }
                std::sort(which, which + k);
                k--;
                for (; k >= 0; k--) {
                    cs.eraseRowCut(which[k]);
                }
                delete [] sort;
                delete [] which;
                numberRowCutsAfter = cs.sizeRowCuts() ;
#else
                double * norm = new double [nCuts];
                int * which = new int [2*nCuts];
                double * score = new double [nCuts];
                double * ortho = new double [nCuts];
                int nIn = 0;
                int nOut = nCuts;
                // For parallel cuts
                double * element2 = new double [numberColumns];
                const double *objective = solver->getObjCoefficients() ;
                double objNorm = 0.0;
                for (int i = 0; i < numberColumns; i++)
                    objNorm += objective[i] * objective[i];
                if (objNorm)
                    objNorm = 1.0 / sqrt(objNorm);
                else
                    objNorm = 1.0;
                objNorm *= 0.1; // weight of 0.1
                CoinZeroN(element2, numberColumns);
                int numberRowCuts = numberRowCutsAfter - numberRowCutsBefore;
                int iBest = -1;
                double best = 0.0;
                int nPossible = 0;
                double testValue = (depth > 1) ? 0.7 : 0.5;
                for (k = 0; k < numberRowCuts; k++) {
                    const OsiRowCut * thisCut = cs.rowCutPtr(k + numberRowCutsBefore) ;
                    double sum = 0.0;
                    if (thisCut->lb() <= thisCut->ub()) {
                        int n = thisCut->row().getNumElements();
                        const int * column = thisCut->row().getIndices();
                        const double * element = thisCut->row().getElements();
                        assert (n);
                        double normThis = 1.0e-6;
                        double obj = 0.0;
                        for (int i = 0; i < n; i++) {
                            int iColumn = column[i];
                            double value = element[i];
                            sum += value * solution[iColumn];
                            normThis += value * value;
                            obj += value * objective[iColumn];
                        }
                        if (sum > thisCut->ub()) {
                            sum = sum - thisCut->ub();
                        } else if (sum < thisCut->lb()) {
                            sum = thisCut->lb() - sum;
                        } else {
                            sum = 0.0;
                        }
                        if (sum) {
                            normThis = 1.0 / sqrt(normThis);
                            norm[k] = normThis;
                            sum *= normThis;
                            obj *= normThis;
                            score[k] = sum + obj * objNorm;
                            ortho[k] = 1.0;
                        }
                    } else {
                        // keep and discard others
                        nIn = 1;
                        which[0] = k;
                        for (int j = 0; j < numberRowCuts; j++) {
                            if (j != k)
                                which[nOut++] = j;
                        }
                        iBest = -1;
                        break;
                    }
                    if (sum) {
                        if (score[k] > best) {
                            best = score[k];
                            iBest = nPossible;
                        }
                        which[nPossible++] = k;
                    } else {
                        which[nOut++] = k;
                    }
                }
                while (iBest >= 0) {
                    int kBest = which[iBest];
                    int j = which[nIn];
                    which[iBest] = j;
                    which[nIn++] = kBest;
                    const OsiRowCut * thisCut = cs.rowCutPtr(kBest + numberRowCutsBefore) ;
                    int n = thisCut->row().getNumElements();
                    nReasonable -= n;
                    if (nReasonable <= 0) {
                        for (k = nIn; k < nPossible; k++)
                            which[nOut++] = which[k];
                        break;
                    }
                    // Now see which ones are too similar and choose next
                    iBest = -1;
                    best = 0.0;
                    int nOld = nPossible;
                    nPossible = nIn;
                    const int * column = thisCut->row().getIndices();
                    const double * element = thisCut->row().getElements();
                    assert (n);
                    double normNew = norm[kBest];
                    for (int i = 0; i < n; i++) {
                        double value = element[i];
                        element2[column[i]] = value;
                    }
                    for (int j = nIn; j < nOld; j++) {
                        k = which[j];
                        const OsiRowCut * thisCut2 = cs.rowCutPtr(k + numberRowCutsBefore) ;
                        int nB = thisCut2->row().getNumElements();
                        const int * columnB = thisCut2->row().getIndices();
                        const double * elementB = thisCut2->row().getElements();
                        assert (nB);
                        double normB = norm[k];
                        double product = 0.0;
                        for (int i = 0; i < nB; i++) {
                            double value = elementB[i];
                            product += value * element2[columnB[i]];
                        }
                        double orthoScore = 1.0 - product * normNew * normB;
                        if (orthoScore >= testValue) {
                            ortho[k] = CoinMin(orthoScore, ortho[k]);
                            double test = score[k] + ortho[k];
                            if (test > best) {
                                best = score[k];
                                iBest = nPossible;
                            }
                            which[nPossible++] = k;
                        } else {
                            which[nOut++] = k;
                        }
                    }
                    for (int i = 0; i < n; i++) {
                        element2[column[i]] = 0.0;
                    }
                }
                delete [] score;
                delete [] ortho;
                std::sort(which + nCuts, which + nOut);
                k = nOut - 1;
                for (; k >= nCuts; k--) {
                    cs.eraseRowCut(which[k] + numberRowCutsBefore);
                }
                delete [] norm;
                delete [] which;
                numberRowCutsAfter = cs.sizeRowCuts() ;
#endif
            }
        }
#ifdef CBC_DEBUG
        {
            int numberRowCutsAfter = cs.sizeRowCuts() ;
            int k ;
            int nBad = 0;
            for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) {
                OsiRowCut thisCut = cs.rowCut(k) ;
                if (thisCut.lb() > thisCut.ub() ||
                        thisCut.lb() > 1.0e8 ||
                        thisCut.ub() < -1.0e8)
                    printf("cut from %s has bounds %g and %g!\n",
                           generatorName_, thisCut.lb(), thisCut.ub());
                if (thisCut.lb() <= thisCut.ub()) {
                    /* check size of elements.
                       We can allow smaller but this helps debug generators as it
                       is unsafe to have small elements */
                    int n = thisCut.row().getNumElements();
                    const int * column = thisCut.row().getIndices();
                    const double * element = thisCut.row().getElements();
                    assert (n);
                    for (int i = 0; i < n; i++) {
                        double value = element[i];
                        if (fabs(value) <= 1.0e-12 || fabs(value) >= 1.0e20)
                            nBad++;
                    }
                }
                if (nBad)
                    printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
                           generatorName_, numberRowCutsAfter - numberRowCutsBefore, nBad);
            }
        }
#endif
	int numberRowCutsAfter = cs.sizeRowCuts() ;
	int numberColumnCutsAfter = cs.sizeColCuts() ;
        if (numberRowCutsBefore < numberRowCutsAfter) {
            for (int k = numberRowCutsBefore; k < numberRowCutsAfter; k++) {
                OsiRowCut thisCut = cs.rowCut(k) ;
                int n = thisCut.row().getNumElements();
                numberElements_ += n;
            }
#ifdef JJF_ZERO
            printf("generator %s generated %d row cuts\n",
                   generatorName_, numberRowCutsAfter - numberRowCutsBefore);
#endif
            numberCuts_ += numberRowCutsAfter - numberRowCutsBefore;
        }
        if (numberColumnCutsBefore < numberColumnCutsAfter) {
#ifdef JJF_ZERO
            printf("generator %s generated %d column cuts\n",
                   generatorName_, numberColumnCutsAfter - numberColumnCutsBefore);
#endif
            numberColumnCuts_ += numberColumnCutsAfter - numberColumnCutsBefore;
        }
        if (timing())
            timeInCutGenerator_ += CoinCpuTime() - time1;
        // switch off if first time and no good
        if (node == NULL && !pass ) {
            if (numberRowCutsAfter - numberRowCutsBefore
		< switchOffIfLessThan_ /*&& numberCuts_ < switchOffIfLessThan_*/) {
	      // switch off
	      maximumTries_ = 0;
	      whenCutGenerator_=-100;
	      //whenCutGenerator_ = -100;
	      //whenCutGeneratorInSub_ = -200;
            }
        }
	if (maximumTries_>0) {
	  maximumTries_--;
	  if (!maximumTries_) 
	    whenCutGenerator_=-100;
	}
    }
    return returnCode;
}
void
CbcCutGenerator::setHowOften(int howOften)
{

    if (howOften >= 1000000) {
        // leave Probing every SCANCUTS_PROBING
        howOften = howOften % 1000000;
        CglProbing* generator =
            dynamic_cast<CglProbing*>(generator_);

        if (generator && howOften > SCANCUTS_PROBING)
            howOften = SCANCUTS_PROBING + 1000000;
        else
            howOften += 1000000;
    }
    whenCutGenerator_ = howOften;
}
void
CbcCutGenerator::setWhatDepth(int value)
{
    depthCutGenerator_ = value;
}
void
CbcCutGenerator::setWhatDepthInSub(int value)
{
    depthCutGeneratorInSub_ = value;
}
// Add in statistics from other
void 
CbcCutGenerator::addStatistics(const CbcCutGenerator * other)
{
  // Time in cut generator
  timeInCutGenerator_ += other->timeInCutGenerator_;
  // Number times cut generator entered
  numberTimes_ += other->numberTimes_;
  // Total number of cuts added
  numberCuts_ += other->numberCuts_;
  // Total number of elements added
  numberElements_ += other->numberElements_;
  // Total number of column cuts added
  numberColumnCuts_ += other->numberColumnCuts_;
  // Total number of cuts active after (at end of n cut passes at each node)
  numberCutsActive_ += other->numberCutsActive_;
  // Number of cuts generated at root
  numberCutsAtRoot_ += other->numberCutsAtRoot_;
  // Number of cuts active at root
  numberActiveCutsAtRoot_ += other->numberActiveCutsAtRoot_;
  // Number of short cuts at root
  numberShortCutsAtRoot_ += other->numberShortCutsAtRoot_;
}
// Scale back statistics by factor
void 
CbcCutGenerator::scaleBackStatistics(int factor)
{
  // leave time
  // Number times cut generator entered
  numberTimes_ = (numberTimes_+factor-1)/factor;
  // Total number of cuts added
  numberCuts_ = (numberCuts_+factor-1)/factor;
  // Total number of elements added
  numberElements_ = (numberElements_+factor-1)/factor;
  // Total number of column cuts added
  numberColumnCuts_ = (numberColumnCuts_+factor-1)/factor;
  // Total number of cuts active after (at end of n cut passes at each node)
  numberCutsActive_ = (numberCutsActive_+factor-1)/factor;
  // Number of cuts generated at root
  numberCutsAtRoot_ = (numberCutsAtRoot_+factor-1)/factor;
  // Number of cuts active at root
  numberActiveCutsAtRoot_ = (numberActiveCutsAtRoot_+factor-1)/factor;
  // Number of short cuts at root
  numberShortCutsAtRoot_ = (numberShortCutsAtRoot_+factor-1)/factor;
}
// Create C++ lines to get to current state
void
CbcCutGenerator::generateTuning( FILE * fp)
{
    fprintf(fp, "// Cbc tuning for generator %s\n", generatorName_);
    fprintf(fp, "   generator->setHowOften(%d);\n", whenCutGenerator_);
    fprintf(fp, "   generator->setSwitchOffIfLessThan(%d);\n", switchOffIfLessThan_);
    fprintf(fp, "   generator->setWhatDepth(%d);\n", depthCutGenerator_);
    fprintf(fp, "   generator->setInaccuracy(%d);\n", inaccuracy_);
    if (timing())
        fprintf(fp, "   generator->setTiming(true);\n");
    if (normal())
        fprintf(fp, "   generator->setNormal(true);\n");
    if (atSolution())
        fprintf(fp, "   generator->setAtSolution(true);\n");
    if (whenInfeasible())
        fprintf(fp, "   generator->setWhenInfeasible(true);\n");
    if (needsOptimalBasis())
        fprintf(fp, "   generator->setNeedsOptimalBasis(true);\n");
    if (mustCallAgain())
        fprintf(fp, "   generator->setMustCallAgain(true);\n");
    if (whetherToUse())
        fprintf(fp, "   generator->setWhetherToUse(true);\n");
}