Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

CbcHeuristicFPump.cpp

Last change on this file was 1883, checked in by stefan, 7 years ago
sync with trunk rev 1882
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 133.5 KB

Rev	Line
[1271]	1	/* $Id: CbcHeuristicFPump.cpp 1883 2013-04-06 13:33:15Z unxusr $ */
[175]	2	// Copyright (C) 2004, International Business Machines
	3	// Corporation and others. All Rights Reserved.
[1573]	4	// This code is licensed under the terms of the Eclipse Public License (EPL).
	5
[175]	6	#if defined(_MSC_VER)
	7	// Turn off compiler warning about long names
	8	# pragma warning(disable:4786)
	9	#endif
	10	#include <cassert>
[904]	11	#include <cstdlib>
[175]	12	#include <cmath>
	13	#include <cfloat>
	14
	15	#include "OsiSolverInterface.hpp"
	16	#include "CbcModel.hpp"
[1088]	17	#ifdef COIN_HAS_CLP
	18	#include "OsiClpSolverInterface.hpp"
	19	#endif
[175]	20	#include "CbcMessage.hpp"
	21	#include "CbcHeuristicFPump.hpp"
	22	#include "CbcBranchActual.hpp"
[1271]	23	#include "CbcBranchDynamic.hpp"
[175]	24	#include "CoinHelperFunctions.hpp"
[738]	25	#include "CoinWarmStartBasis.hpp"
[175]	26	#include "CoinTime.hpp"
[747]	27	#include "CbcEventHandler.hpp"
[175]	28
	29
	30	// Default Constructor
[1286]	31	CbcHeuristicFPump::CbcHeuristicFPump()
	32	: CbcHeuristic(),
	33	startTime_(0.0),
	34	maximumTime_(0.0),
	35	fakeCutoff_(COIN_DBL_MAX),
	36	absoluteIncrement_(0.0),
	37	relativeIncrement_(0.0),
	38	defaultRounding_(0.49999),
	39	initialWeight_(0.0),
	40	weightFactor_(0.1),
	41	artificialCost_(COIN_DBL_MAX),
	42	iterationRatio_(0.0),
	43	reducedCostMultiplier_(1.0),
	44	maximumPasses_(100),
	45	maximumRetries_(1),
	46	accumulate_(0),
	47	fixOnReducedCosts_(1),
	48	roundExpensive_(false)
[175]	49	{
[1286]	50	setWhen(1);
[175]	51	}
	52
	53	// Constructor from model
	54	CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
[1286]	55	double downValue, bool roundExpensive)
	56	: CbcHeuristic(model),
	57	startTime_(0.0),
	58	maximumTime_(0.0),
	59	fakeCutoff_(COIN_DBL_MAX),
	60	absoluteIncrement_(0.0),
	61	relativeIncrement_(0.0),
	62	defaultRounding_(downValue),
	63	initialWeight_(0.0),
	64	weightFactor_(0.1),
	65	artificialCost_(COIN_DBL_MAX),
	66	iterationRatio_(0.0),
	67	reducedCostMultiplier_(1.0),
	68	maximumPasses_(100),
	69	maximumRetries_(1),
	70	accumulate_(0),
	71	fixOnReducedCosts_(1),
	72	roundExpensive_(roundExpensive)
[175]	73	{
[1286]	74	setWhen(1);
[175]	75	}
	76
[1286]	77	// Destructor
[175]	78	CbcHeuristicFPump::~CbcHeuristicFPump ()
	79	{
	80	}
	81
	82	// Clone
	83	CbcHeuristic *
	84	CbcHeuristicFPump::clone() const
	85	{
[1286]	86	return new CbcHeuristicFPump(*this);
[175]	87	}
[356]	88	// Create C++ lines to get to current state
[1286]	89	void
	90	CbcHeuristicFPump::generateCpp( FILE * fp)
[356]	91	{
[1286]	92	CbcHeuristicFPump other;
	93	fprintf(fp, "0#include \"CbcHeuristicFPump.hpp\"\n");
	94	fprintf(fp, "3 CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
	95	CbcHeuristic::generateCpp(fp, "heuristicFPump");
	96	if (maximumPasses_ != other.maximumPasses_)
	97	fprintf(fp, "3 heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
	98	else
	99	fprintf(fp, "4 heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
	100	if (maximumRetries_ != other.maximumRetries_)
	101	fprintf(fp, "3 heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
	102	else
	103	fprintf(fp, "4 heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
	104	if (accumulate_ != other.accumulate_)
	105	fprintf(fp, "3 heuristicFPump.setAccumulate(%d);\n", accumulate_);
	106	else
	107	fprintf(fp, "4 heuristicFPump.setAccumulate(%d);\n", accumulate_);
	108	if (fixOnReducedCosts_ != other.fixOnReducedCosts_)
	109	fprintf(fp, "3 heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
	110	else
	111	fprintf(fp, "4 heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
	112	if (maximumTime_ != other.maximumTime_)
	113	fprintf(fp, "3 heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
	114	else
	115	fprintf(fp, "4 heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
	116	if (fakeCutoff_ != other.fakeCutoff_)
	117	fprintf(fp, "3 heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
	118	else
	119	fprintf(fp, "4 heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
	120	if (absoluteIncrement_ != other.absoluteIncrement_)
	121	fprintf(fp, "3 heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
	122	else
	123	fprintf(fp, "4 heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
	124	if (relativeIncrement_ != other.relativeIncrement_)
	125	fprintf(fp, "3 heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
	126	else
	127	fprintf(fp, "4 heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
	128	if (defaultRounding_ != other.defaultRounding_)
	129	fprintf(fp, "3 heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
	130	else
	131	fprintf(fp, "4 heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
	132	if (initialWeight_ != other.initialWeight_)
	133	fprintf(fp, "3 heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
	134	else
	135	fprintf(fp, "4 heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
	136	if (weightFactor_ != other.weightFactor_)
	137	fprintf(fp, "3 heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
	138	else
	139	fprintf(fp, "4 heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
	140	if (artificialCost_ != other.artificialCost_)
	141	fprintf(fp, "3 heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
	142	else
	143	fprintf(fp, "4 heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
	144	if (iterationRatio_ != other.iterationRatio_)
	145	fprintf(fp, "3 heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
	146	else
	147	fprintf(fp, "4 heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
	148	if (reducedCostMultiplier_ != other.reducedCostMultiplier_)
	149	fprintf(fp, "3 heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
	150	else
	151	fprintf(fp, "4 heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
	152	fprintf(fp, "3 cbcModel->addHeuristic(&heuristicFPump);\n");
[356]	153	}
[175]	154
[1286]	155	// Copy constructor
[175]	156	CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
[1286]	157	:
	158	CbcHeuristic(rhs),
	159	startTime_(rhs.startTime_),
	160	maximumTime_(rhs.maximumTime_),
	161	fakeCutoff_(rhs.fakeCutoff_),
	162	absoluteIncrement_(rhs.absoluteIncrement_),
	163	relativeIncrement_(rhs.relativeIncrement_),
	164	defaultRounding_(rhs.defaultRounding_),
	165	initialWeight_(rhs.initialWeight_),
	166	weightFactor_(rhs.weightFactor_),
	167	artificialCost_(rhs.artificialCost_),
	168	iterationRatio_(rhs.iterationRatio_),
	169	reducedCostMultiplier_(rhs.reducedCostMultiplier_),
	170	maximumPasses_(rhs.maximumPasses_),
	171	maximumRetries_(rhs.maximumRetries_),
	172	accumulate_(rhs.accumulate_),
	173	fixOnReducedCosts_(rhs.fixOnReducedCosts_),
	174	roundExpensive_(rhs.roundExpensive_)
[175]	175	{
	176	}
[640]	177
[1286]	178	// Assignment operator
	179	CbcHeuristicFPump &
	180	CbcHeuristicFPump::operator=( const CbcHeuristicFPump & rhs)
[640]	181	{
[1286]	182	if (this != &rhs) {
	183	CbcHeuristic::operator=(rhs);
	184	startTime_ = rhs.startTime_;
	185	maximumTime_ = rhs.maximumTime_;
	186	fakeCutoff_ = rhs.fakeCutoff_;
	187	absoluteIncrement_ = rhs.absoluteIncrement_;
	188	relativeIncrement_ = rhs.relativeIncrement_;
	189	defaultRounding_ = rhs.defaultRounding_;
	190	initialWeight_ = rhs.initialWeight_;
	191	weightFactor_ = rhs.weightFactor_;
	192	artificialCost_ = rhs.artificialCost_;
	193	iterationRatio_ = rhs.iterationRatio_;
	194	reducedCostMultiplier_ = rhs.reducedCostMultiplier_;
	195	maximumPasses_ = rhs.maximumPasses_;
	196	maximumRetries_ = rhs.maximumRetries_;
	197	accumulate_ = rhs.accumulate_;
	198	fixOnReducedCosts_ = rhs.fixOnReducedCosts_;
	199	roundExpensive_ = rhs.roundExpensive_;
	200	}
	201	return *this;
[640]	202	}
	203
[175]	204	// Resets stuff if model changes
[1286]	205	void
[1271]	206	CbcHeuristicFPump::resetModel(CbcModel * )
[175]	207	{
	208	}
	209
	210	/************************BEGIN MAIN PROCEDURE *********************************/
	211
	212	// See if feasibility pump will give better solution
	213	// Sets value of solution
	214	// Returns 1 if solution, 0 if not
	215	int
	216	CbcHeuristicFPump::solution(double & solutionValue,
[1286]	217	double * betterSolution)
[175]	218	{
[1286]	219	startTime_ = CoinCpuTime();
	220	numCouldRun_++;
	221	double incomingObjective = solutionValue;
[698]	222	#define LEN_PRINT 200
[1286]	223	char pumpPrint[LEN_PRINT];
	224	pumpPrint[0] = '\0';
[1364]	225	/*
	226	Decide if we want to run. Standard values for when are described in
	227	CbcHeuristic.hpp. If we're off, or running only at root and this isn't the
	228	root, bail out.
	229
	230	The double test (against phase, then atRoot and passNumber) has a fair bit
	231	of redundancy, but the results will differ depending on whether we're
	232	actually at the root of the main search tree or at the root of a small tree
	233	(recursive call to branchAndBound).
	234
	235	FPump also supports some exotic values (11 -- 15) for when, described in
	236	CbcHeuristicFPump.hpp.
	237	*/
[1286]	238	if (!when() \|\| (when() == 1 && model_->phase() != 1))
	239	return 0; // switched off
	240	// See if at root node
	241	bool atRoot = model_->getNodeCount() == 0;
	242	int passNumber = model_->getCurrentPassNumber();
	243	// just do once
	244	if (!atRoot)
	245	return 0;
	246	int options = feasibilityPumpOptions_;
	247	if ((options % 1000000) > 0) {
	248	int kOption = options / 1000000;
	249	options = options % 1000000;
	250	/*
	251	Add 10 to do even if solution
	252	1 - do after cuts
	253	2 - do after cuts (not before)
	254	3 - not used do after every cut round (and after cuts)
	255	k not used do after every (k-2)th round
	256	*/
	257	if (kOption < 10 && model_->getSolutionCount())
	258	return 0;
	259	if (model_->getSolutionCount())
	260	kOption = kOption % 10;
	261	bool good;
	262	if (kOption == 1) {
	263	good = (passNumber == 999999);
	264	} else if (kOption == 2) {
	265	good = (passNumber == 999999);
	266	passNumber = 2; // so won't run before
	267	//} else if (kOption==3) {
	268	//good = true;
	269	} else {
	270	//good = (((passNumber-1)%(kOption-2))==0);
	271	good = false;
	272	}
	273	if (passNumber != 1 && !good)
	274	return 0;
[1271]	275	} else {
[1286]	276	if (passNumber != 1)
	277	return 0;
[1271]	278	}
[1286]	279	// loop round doing repeated pumps
	280	double cutoff;
	281	model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
	282	double direction = model_->solver()->getObjSense();
	283	cutoff *= direction;
	284	int numberBandBsolutions = 0;
	285	double firstCutoff = fabs(cutoff);
	286	cutoff = CoinMin(cutoff, solutionValue);
	287	// check plausible and space for rounded solution
	288	int numberColumns = model_->getNumCols();
	289	int numberIntegers = model_->numberIntegers();
	290	const int * integerVariableOrig = model_->integerVariable();
	291	double iterationLimit = -1.0;
	292	//iterationRatio_=1.0;
	293	if (iterationRatio_ > 0.0)
	294	iterationLimit = (2 * model_->solver()->getNumRows() + 2 * numberColumns) *
	295	iterationRatio_;
	296	int totalNumberIterations = 0;
	297	int averageIterationsPerTry = -1;
	298	int numberIterationsLastPass = 0;
	299	// 1. initially check 0-1
[1364]	300	/*
	301	I'm skeptical of the above comment, but it's likely accurate as the default.
	302	Bit 4 or bit 8 needs to be set in order to consider working with general
	303	integers.
	304	*/
[1286]	305	int i, j;
	306	int general = 0;
	307	int * integerVariable = new int[numberIntegers];
	308	const double * lower = model_->solver()->getColLower();
	309	const double * upper = model_->solver()->getColUpper();
	310	bool doGeneral = (accumulate_ & 4) != 0;
	311	j = 0;
[1364]	312	/*
	313	Scan the objects, recording the columns and counting general integers.
	314
	315	Seems like the NDEBUG tests could be made into an applicability test. If
	316	a scan of the objects reveals complex objects, just clean up and return
	317	failure.
	318	*/
[1286]	319	for (i = 0; i < numberIntegers; i++) {
	320	int iColumn = integerVariableOrig[i];
[356]	321	#ifndef NDEBUG
[1286]	322	const OsiObject * object = model_->object(i);
	323	const CbcSimpleInteger * integerObject =
	324	dynamic_cast<const CbcSimpleInteger *> (object);
	325	const OsiSimpleInteger * integerObject2 =
	326	dynamic_cast<const OsiSimpleInteger *> (object);
	327	assert(integerObject \|\| integerObject2);
[201]	328	#endif
[1286]	329	if (upper[iColumn] - lower[iColumn] > 1.000001) {
	330	general++;
	331	if (doGeneral)
	332	integerVariable[j++] = iColumn;
	333	} else {
	334	integerVariable[j++] = iColumn;
	335	}
[175]	336	}
[1364]	337	/*
	338	If 2/3 of integers are general integers, and we're not going to work with
	339	them, might as well go home.
	340
	341	The else case is unclear to me. We reach it if general integers are less than
	342	2/3 of the total, or if either of bit 4 or 8 is set. But only bit 8 is used
	343	in the decision. (Let manyGen = 1 if more than 2/3 of integers are general
	344	integers. Then a k-map on manyGen, bit4, and bit8 shows it clearly.)
	345
	346	So there's something odd here. In the case where bit4 = 1 and bit8 = 0,
	347	we've included general integers in integerVariable, but we're not going to
	348	process them.
	349	*/
[1286]	350	if (general3 > 2numberIntegers && !doGeneral) {
	351	delete [] integerVariable;
	352	return 0;
	353	} else if ((accumulate_&4) == 0) {
	354	doGeneral = false;
	355	j = 0;
	356	for (i = 0; i < numberIntegers; i++) {
	357	int iColumn = integerVariableOrig[i];
	358	if (upper[iColumn] - lower[iColumn] < 1.000001)
	359	integerVariable[j++] = iColumn;
	360	}
[1053]	361	}
[1286]	362	if (!general)
	363	doGeneral = false;
[1053]	364	#ifdef CLP_INVESTIGATE
[1286]	365	if (doGeneral)
	366	printf("DOing general with %d out of %d\n", general, numberIntegers);
[1053]	367	#endif
[1364]	368	/*
	369	This `closest solution' will satisfy integrality, but violate some other
	370	constraints?
	371	*/
[1286]	372	// For solution closest to feasible if none found
	373	int * closestSolution = general ? NULL : new int[numberIntegers];
	374	double closestObjectiveValue = COIN_DBL_MAX;
[1271]	375
[1286]	376	int numberIntegersOrig = numberIntegers;
	377	numberIntegers = j;
	378	double * newSolution = new double [numberColumns];
	379	double newSolutionValue = COIN_DBL_MAX;
[1802]	380	int maxSolutions = model_->getMaximumSolutions();
	381	int numberSolutions=0;
[1286]	382	bool solutionFound = false;
	383	int * usedColumn = NULL;
	384	double * lastSolution = NULL;
	385	int fixContinuous = 0;
	386	bool fixInternal = false;
	387	bool allSlack = false;
	388	if (when_ >= 21 && when_ <= 25) {
	389	when_ -= 10;
	390	allSlack = true;
[753]	391	}
[1286]	392	double time1 = CoinCpuTime();
[1364]	393	/*
	394	Obtain a relaxed lp solution.
	395	*/
[1286]	396	model_->solver()->resolve();
	397	if (!model_->solver()->isProvenOptimal()) {
	398	// presumably max time or some such
	399	return 0;
[1115]	400	}
[1286]	401	numRuns_++;
	402	if (cutoff < 1.0e50 && false) {
	403	// Fix on djs
	404	double direction = model_->solver()->getObjSense() ;
	405	double gap = cutoff - model_->solver()->getObjValue() * direction ;
	406	double tolerance;
	407	model_->solver()->getDblParam(OsiDualTolerance, tolerance) ;
	408	if (gap > 0.0) {
	409	gap += 100.0 * tolerance;
	410	int nFix = model_->solver()->reducedCostFix(gap);
	411	printf("dj fixing fixed %d variables\n", nFix);
	412	}
	413	}
[1364]	414	/*
	415	I have no idea why we're doing this, except perhaps that saveBasis will be
	416	automagically deleted on exit from the routine.
	417	*/
[1286]	418	CoinWarmStartBasis saveBasis;
	419	CoinWarmStartBasis * basis =
	420	dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ;
	421	if (basis) {
	422	saveBasis = * basis;
	423	delete basis;
	424	}
	425	double continuousObjectiveValue = model_->solver()->getObjValue() * model_->solver()->getObjSense();
	426	double * firstPerturbedObjective = NULL;
	427	double * firstPerturbedSolution = NULL;
	428	double firstPerturbedValue = COIN_DBL_MAX;
	429	if (when_ >= 11 && when_ <= 15) {
	430	fixInternal = when_ > 11 && when_ < 15;
	431	if (when_ < 13)
	432	fixContinuous = 0;
	433	else if (when_ != 14)
	434	fixContinuous = 1;
	435	else
	436	fixContinuous = 2;
	437	when_ = 1;
	438	if ((accumulate_&1) != 0) {
	439	usedColumn = new int [numberColumns];
	440	for (int i = 0; i < numberColumns; i++)
	441	usedColumn[i] = -1;
	442	}
	443	lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(), numberColumns);
	444	}
	445	int finalReturnCode = 0;
	446	int totalNumberPasses = 0;
	447	int numberTries = 0;
	448	CoinWarmStartBasis bestBasis;
	449	bool exitAll = false;
	450	//double saveBestObjective = model_->getMinimizationObjValue();
	451	OsiSolverInterface * solver = NULL;
	452	double artificialFactor = 0.00001;
	453	// also try rounding!
	454	double * roundingSolution = new double[numberColumns];
	455	double roundingObjective = solutionValue;
	456	CbcRounding roundingHeuristic(*model_);
	457	int dualPass = 0;
	458	int secondPassOpt = 0;
[1271]	459	#define RAND_RAND
	460	#ifdef RAND_RAND
[1286]	461	int offRandom = 0;
[1271]	462	#endif
[1286]	463	int maximumAllowed = -1;
	464	bool moreIterations = false;
	465	if (options > 0) {
	466	if (options >= 1000)
	467	maximumAllowed = options / 1000;
	468	int options2 = (options % 1000) / 100;
[1271]	469	#ifdef RAND_RAND
[1286]	470	offRandom = options2 & 1;
[1271]	471	#endif
[1286]	472	moreIterations = (options2 & 2) != 0;
	473	secondPassOpt = (options / 10) % 10;
	474	/* 1 to 7 - re-use solution
	475	8 use dual and current solution(ish)
	476	9 use dual and allslack
	477	1 - primal and mod obj
	478	2 - dual and mod obj
	479	3 - primal and no mod obj
	480	add 3 to redo current solution
	481	*/
	482	if (secondPassOpt >= 8) {
	483	dualPass = secondPassOpt - 7;
	484	secondPassOpt = 0;
	485	}
[1053]	486	}
[1286]	487	// Number of passes to do
	488	int maximumPasses = maximumPasses_;
[1271]	489	#ifdef COIN_HAS_CLP
[1813]	490	{
	491	OsiClpSolverInterface * clpSolver
[1286]	492	= dynamic_cast<OsiClpSolverInterface *> (model_->solver());
[1813]	493	if (clpSolver ) {
	494	if (maximumPasses == 30) {
	495	if (clpSolver->fakeObjective())
[1286]	496	maximumPasses = 100; // feasibility problem?
[1813]	497	}
[1883]	498	randomNumberGenerator_.randomize();
[1813]	499	if (model_->getRandomSeed()!=-1)
	500	clpSolver->getModelPtr()->setRandomSeed(randomNumberGenerator_.getSeed());
[1883]	501	clpSolver->getModelPtr()->randomNumberGenerator()->randomize();
[1813]	502	}
[1286]	503	}
[1271]	504	#endif
	505	#ifdef RAND_RAND
[1286]	506	double * randomFactor = new double [numberColumns];
	507	for (int i = 0; i < numberColumns; i++) {
	508	double value = floor(1.0e3 * randomNumberGenerator_.randomDouble());
	509	randomFactor[i] = 1.0 + value * 1.0e-4;
	510	}
[1271]	511	#endif
[1286]	512	// guess exact multiple of objective
	513	double exactMultiple = model_->getCutoffIncrement();
	514	exactMultiple *= 2520;
	515	if (fabs(exactMultiple / 0.999 - floor(exactMultiple / 0.999 + 0.5)) < 1.0e-9)
	516	exactMultiple /= 2520.0 * 0.999;
	517	else if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-9)
	518	exactMultiple /= 2520.0;
	519	else
	520	exactMultiple = 0.0;
	521	// check for rounding errors (only for integral case)
	522	if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-8)
	523	exactMultiple = floor(exactMultiple + 0.5);
	524	//printf("exact multiple %g\n",exactMultiple);
	525	// Clone solver for rounding
	526	OsiSolverInterface * clonedSolver = cloneBut(2); // wasmodel_->solver()->clone();
	527	while (!exitAll) {
	528	// Cutoff rhs
	529	double useRhs = COIN_DBL_MAX;
	530	double useOffset = 0.0;
	531	int numberPasses = 0;
	532	artificialFactor *= 10.0;
	533	int lastMove = (!numberTries) ? -10 : 1000000;
	534	double lastSumInfeas = COIN_DBL_MAX;
	535	numberTries++;
	536	// Clone solver - otherwise annoys root node computations
	537	solver = cloneBut(2); // was model_->solver()->clone();
[1053]	538	#ifdef COIN_HAS_CLP
[1286]	539	{
	540	OsiClpSolverInterface * clpSolver
	541	= dynamic_cast<OsiClpSolverInterface *> (solver);
	542	if (clpSolver) {
	543	// better to clean up using primal?
	544	ClpSimplex * lp = clpSolver->getModelPtr();
	545	int options = lp->specialOptions();
	546	lp->setSpecialOptions(options \| 8192);
	547	//lp->setSpecialOptions(options\|0x01000000);
[1271]	548	#ifdef CLP_INVESTIGATE
[1286]	549	clpSolver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
	550	lp->setLogLevel(CoinMax(1, lp->logLevel()));
[1271]	551	#endif
[1286]	552	}
	553	}
[1053]	554	#endif
[1286]	555	if (CoinMin(fakeCutoff_, cutoff) < 1.0e50) {
	556	// Fix on djs
	557	double direction = solver->getObjSense() ;
	558	double gap = CoinMin(fakeCutoff_, cutoff) - solver->getObjValue() * direction ;
	559	double tolerance;
	560	solver->getDblParam(OsiDualTolerance, tolerance) ;
	561	if (gap > 0.0 && (fixOnReducedCosts_ == 1 \|\| (numberTries == 1 && fixOnReducedCosts_ == 2))) {
	562	gap += 100.0 * tolerance;
	563	gap *= reducedCostMultiplier_;
	564	int nFix = solver->reducedCostFix(gap);
	565	if (nFix) {
	566	sprintf(pumpPrint, "Reduced cost fixing fixed %d variables on major pass %d", nFix, numberTries);
	567	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	568	<< pumpPrint
	569	<< CoinMessageEol;
	570	//pumpPrint[0]='\0';
	571	}
	572	}
	573	}
	574	// if cutoff exists then add constraint
	575	bool useCutoff = (fabs(cutoff) < 1.0e20 && (fakeCutoff_ != COIN_DBL_MAX \|\| numberTries > 1));
	576	// but there may be a close one
	577	if (firstCutoff < 2.0*solutionValue && numberTries == 1 && CoinMin(cutoff, fakeCutoff_) < 1.0e20)
	578	useCutoff = true;
	579	if (useCutoff) {
	580	double rhs = CoinMin(cutoff, fakeCutoff_);
	581	const double * objective = solver->getObjCoefficients();
	582	int numberColumns = solver->getNumCols();
	583	int * which = new int[numberColumns];
	584	double * els = new double[numberColumns];
	585	int nel = 0;
	586	for (int i = 0; i < numberColumns; i++) {
	587	double value = objective[i];
	588	if (value) {
	589	which[nel] = i;
	590	els[nel++] = direction * value;
	591	}
	592	}
	593	solver->getDblParam(OsiObjOffset, useOffset);
[753]	594	#ifdef COIN_DEVELOP
[1286]	595	if (useOffset)
	596	printf("CbcHeuristicFPump obj offset %g\n", useOffset);
[753]	597	#endif
[1286]	598	useOffset *= direction;
	599	// Tweak rhs and save
	600	useRhs = rhs;
[1393]	601	#ifdef JJF_ZERO
[1286]	602	double tempValue = 60.0 * useRhs;
	603	if (fabs(tempValue - floor(tempValue + 0.5)) < 1.0e-7 && rhs != fakeCutoff_) {
	604	// add a little
	605	useRhs += 1.0e-5;
	606	}
[1117]	607	#endif
[1286]	608	solver->addRow(nel, which, els, -COIN_DBL_MAX, useRhs + useOffset);
	609	delete [] which;
	610	delete [] els;
	611	bool takeHint;
	612	OsiHintStrength strength;
	613	solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
	614	solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
	615	solver->resolve();
	616	solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
	617	if (!solver->isProvenOptimal()) {
	618	// presumably max time or some such
	619	exitAll = true;
	620	break;
	621	}
	622	}
	623	solver->setDblParam(OsiDualObjectiveLimit, 1.0e50);
	624	solver->resolve();
	625	// Solver may not be feasible
	626	if (!solver->isProvenOptimal()) {
	627	exitAll = true;
	628	break;
	629	}
	630	const double * lower = solver->getColLower();
	631	const double * upper = solver->getColUpper();
	632	const double * solution = solver->getColSolution();
	633	if (lastSolution)
	634	memcpy(lastSolution, solution, numberColumns*sizeof(double));
	635	double primalTolerance;
	636	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
	637
	638	// 2 space for last rounded solutions
[640]	639	#define NUMBER_OLD 4
[1286]	640	double ** oldSolution = new double * [NUMBER_OLD];
	641	for (j = 0; j < NUMBER_OLD; j++) {
	642	oldSolution[j] = new double[numberColumns];
	643	for (i = 0; i < numberColumns; i++) oldSolution[j][i] = -COIN_DBL_MAX;
	644	}
	645
	646	// 3. Replace objective with an initial 0-valued objective
	647	double * saveObjective = new double [numberColumns];
	648	memcpy(saveObjective, solver->getObjCoefficients(), numberColumns*sizeof(double));
	649	for (i = 0; i < numberColumns; i++) {
	650	solver->setObjCoeff(i, 0.0);
	651	}
	652	bool finished = false;
	653	double direction = solver->getObjSense();
	654	int returnCode = 0;
	655	bool takeHint;
	656	OsiHintStrength strength;
	657	solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
	658	solver->setHintParam(OsiDoDualInResolve, false);
	659	//solver->messageHandler()->setLogLevel(0);
	660
	661	// 4. Save objective offset so we can see progress
	662	double saveOffset;
	663	solver->getDblParam(OsiObjOffset, saveOffset);
	664	// Get amount for original objective
	665	double scaleFactor = 0.0;
[833]	666	#ifdef COIN_DEVELOP
[1286]	667	double largestCost = 0.0;
	668	int nArtificial = 0;
[833]	669	#endif
[1286]	670	for (i = 0; i < numberColumns; i++) {
	671	double value = saveObjective[i];
	672	scaleFactor += value * value;
[833]	673	#ifdef COIN_DEVELOP
[1286]	674	largestCost = CoinMax(largestCost, fabs(value));
	675	if (value*direction >= artificialCost_)
	676	nArtificial++;
[833]	677	#endif
[1286]	678	}
	679	if (scaleFactor)
	680	scaleFactor = (initialWeight_ * sqrt(static_cast<double> (numberIntegers))) / sqrt(scaleFactor);
[1271]	681	#ifdef CLP_INVESTIGATE
[833]	682	#ifdef COIN_DEVELOP
[1286]	683	if (scaleFactor \|\| nArtificial)
	684	printf("Using %g fraction of original objective (decay %g) - largest %g - %d artificials\n", scaleFactor, weightFactor_,
	685	largestCost, nArtificial);
[1271]	686	#else
[1286]	687	if (scaleFactor)
	688	printf("Using %g fraction of original objective (decay %g)\n",
	689	scaleFactor, weightFactor_);
[833]	690	#endif
[1271]	691	#endif
[1286]	692	// This is an array of sums of infeasibilities so can see if "bobbling"
[1127]	693	#define SIZE_BOBBLE 20
[1286]	694	double saveSumInf[SIZE_BOBBLE];
	695	CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
	696	// 0 before doing anything
	697	int bobbleMode = 0;
	698	// 5. MAIN WHILE LOOP
	699	//bool newLineNeeded=false;
[1364]	700	/*
	701	finished occurs exactly twice in this routine: immediately above, where it's
	702	set to false, and here in the loop condition.
	703	*/
[1286]	704	while (!finished) {
	705	double newTrueSolutionValue = 0.0;
	706	double newSumInfeas = 0.0;
	707	int newNumberInfeas = 0;
	708	returnCode = 0;
	709	if (model_->maximumSecondsReached()) {
	710	exitAll = true;
	711	break;
	712	}
	713	// see what changed
	714	if (usedColumn) {
	715	for (i = 0; i < numberColumns; i++) {
	716	if (fabs(solution[i] - lastSolution[i]) > 1.0e-8)
	717	usedColumn[i] = numberPasses;
	718	lastSolution[i] = solution[i];
	719	}
	720	}
	721	if (averageIterationsPerTry >= 0) {
	722	int n = totalNumberIterations - numberIterationsLastPass;
	723	double perPass = totalNumberIterations /
	724	(totalNumberPasses + numberPasses + 1.0e-5);
	725	perPass /= (solver->getNumRows() + numberColumns);
	726	double test = moreIterations ? 0.3 : 0.05;
	727	if (n > CoinMax(20000, 3*averageIterationsPerTry)
	728	&& (switches_&2) == 0 && maximumPasses<200 && perPass>test) {
	729	exitAll = true;
	730	}
	731	}
	732	// Exit on exact total number if maximumPasses large
	733	if ((maximumPasses >= 200 \|\| (switches_&2) != 0)
	734	&& numberPasses + totalNumberPasses >=
	735	maximumPasses)
	736	exitAll = true;
	737	bool exitThis = false;
	738	if (iterationLimit < 0.0) {
	739	if (numberPasses >= maximumPasses) {
	740	// If going well then keep going if maximumPasses small
	741	if (lastMove < numberPasses - 4 \|\| lastMove == 1000000)
	742	exitThis = true;
	743	if (maximumPasses > 20 \|\| numberPasses >= 40)
	744	exitThis = true;
	745	}
	746	}
	747	if (iterationLimit > 0.0 && totalNumberIterations > iterationLimit
	748	&& numberPasses > 15) {
	749	// exiting on iteration count
	750	exitAll = true;
	751	} else if (maximumPasses<30 && numberPasses>100) {
	752	// too many passes anyway
	753	exitAll = true;
	754	}
	755	if (maximumTime_ > 0.0 && CoinCpuTime() >= startTime_ + maximumTime_) {
	756	exitAll = true;
	757	// force exit
	758	switches_ \|= 2048;
	759	}
	760	if (exitAll \|\| exitThis)
	761	break;
	762	memcpy(newSolution, solution, numberColumns*sizeof(double));
	763	int flip;
	764	if (numberPasses == 0 && false) {
	765	// always use same seed
	766	randomNumberGenerator_.setSeed(987654321);
	767	}
	768	returnCode = rounds(solver, newSolution,/saveObjective,/
	769	numberIntegers, integerVariable,
	770	/pumpPrint,/numberPasses,
	771	/roundExpensive_,/defaultRounding_, &flip);
	772	if (numberPasses == 0 && false) {
	773	// Make sure random will be different
	774	for (i = 1; i < numberTries; i++)
	775	randomNumberGenerator_.randomDouble();
	776	}
	777	numberPasses++;
	778	if (returnCode) {
	779	// SOLUTION IS INTEGER
	780	// Put back correct objective
	781	for (i = 0; i < numberColumns; i++)
	782	solver->setObjCoeff(i, saveObjective[i]);
[640]	783
[1286]	784	// solution - but may not be better
	785	// Compute using dot product
	786	solver->setDblParam(OsiObjOffset, saveOffset);
	787	newSolutionValue = -saveOffset;
	788	for ( i = 0 ; i < numberColumns ; i++ )
	789	newSolutionValue += saveObjective[i] * newSolution[i];
	790	newSolutionValue *= direction;
	791	sprintf(pumpPrint, "Solution found of %g", newSolutionValue);
	792	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	793	<< pumpPrint
	794	<< CoinMessageEol;
	795	//newLineNeeded=false;
	796	if (newSolutionValue < solutionValue) {
	797	double saveValue = solutionValue;
	798	if (!doGeneral) {
	799	int numberLeft = 0;
	800	for (i = 0; i < numberIntegersOrig; i++) {
	801	int iColumn = integerVariableOrig[i];
	802	double value = floor(newSolution[iColumn] + 0.5);
	803	if (solver->isBinary(iColumn)) {
	804	solver->setColLower(iColumn, value);
	805	solver->setColUpper(iColumn, value);
	806	} else {
	807	if (fabs(value - newSolution[iColumn]) > 1.0e-7)
	808	numberLeft++;
	809	}
	810	}
	811	if (numberLeft) {
	812	sprintf(pumpPrint, "Branch and bound needed to clear up %d general integers", numberLeft);
	813	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	814	<< pumpPrint
	815	<< CoinMessageEol;
	816	returnCode = smallBranchAndBound(solver, numberNodes_, newSolution, newSolutionValue,
	817	solutionValue, "CbcHeuristicFpump");
	818	if (returnCode < 0) {
	819	if (returnCode == -2)
	820	exitAll = true;
	821	returnCode = 0; // returned on size or event
	822	}
	823	if ((returnCode&2) != 0) {
	824	// could add cut
	825	returnCode &= ~2;
	826	}
	827	if (returnCode != 1)
	828	newSolutionValue = saveValue;
	829	if (returnCode && newSolutionValue < saveValue)
	830	numberBandBsolutions++;
	831	}
	832	}
	833	if (returnCode && newSolutionValue < saveValue) {
	834	memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
	835	solutionFound = true;
	836	if (exitNow(newSolutionValue))
	837	exitAll = true;
	838	CoinWarmStartBasis * basis =
	839	dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
	840	if (basis) {
	841	bestBasis = * basis;
	842	delete basis;
	843	int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
	844	if (action == 0) {
	845	double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
	846	double saveObjectiveValue = model_->getMinimizationObjValue();
	847	model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
	848	if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
	849	model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
	850	delete [] saveOldSolution;
	851	}
	852	if (action == 0 \|\| model_->maximumSecondsReached()) {
	853	exitAll = true; // exit
	854	break;
	855	}
	856	}
	857	if ((accumulate_&1) != 0) {
	858	model_->incrementUsed(betterSolution); // for local search
	859	}
	860	solutionValue = newSolutionValue;
	861	solutionFound = true;
[1802]	862	numberSolutions++;
	863	if (numberSolutions>=maxSolutions)
	864	exitAll = true;
[1286]	865	if (general && saveValue != newSolutionValue) {
	866	sprintf(pumpPrint, "Cleaned solution of %g", solutionValue);
	867	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	868	<< pumpPrint
	869	<< CoinMessageEol;
	870	}
	871	if (exitNow(newSolutionValue))
	872	exitAll = true;
	873	} else {
	874	sprintf(pumpPrint, "Mini branch and bound could not fix general integers");
	875	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	876	<< pumpPrint
	877	<< CoinMessageEol;
	878	}
	879	} else {
	880	sprintf(pumpPrint, "After further testing solution no better than previous of %g", solutionValue);
	881	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	882	<< pumpPrint
	883	<< CoinMessageEol;
	884	//newLineNeeded=false;
	885	returnCode = 0;
	886	}
	887	break;
	888	} else {
	889	// SOLUTION IS not INTEGER
	890	// 1. check for loop
	891	bool matched;
	892	for (int k = NUMBER_OLD - 1; k > 0; k--) {
	893	double * b = oldSolution[k];
	894	matched = true;
	895	for (i = 0; i < numberIntegers; i++) {
	896	int iColumn = integerVariable[i];
	897	if (newSolution[iColumn] != b[iColumn]) {
	898	matched = false;
	899	break;
	900	}
	901	}
	902	if (matched) break;
	903	}
	904	int numberPerturbed = 0;
	905	if (matched \|\| numberPasses % 100 == 0) {
	906	// perturbation
	907	//sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
	908	//newLineNeeded=true;
	909	double factorX[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
	910	double factor = 1.0;
	911	double target = -1.0;
	912	double * randomX = new double [numberIntegers];
	913	for (i = 0; i < numberIntegers; i++)
	914	randomX[i] = CoinMax(0.0, randomNumberGenerator_.randomDouble() - 0.3);
	915	for (int k = 0; k < 10; k++) {
[753]	916	#ifdef COIN_DEVELOP_x
[1286]	917	printf("kpass %d\n", k);
[750]	918	#endif
[1286]	919	int numberX[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
	920	for (i = 0; i < numberIntegers; i++) {
	921	int iColumn = integerVariable[i];
	922	double value = randomX[i];
	923	double difference = fabs(solution[iColumn] - newSolution[iColumn]);
	924	for (int j = 0; j < 10; j++) {
	925	if (difference + value*factorX[j] > 0.5)
	926	numberX[j]++;
	927	}
	928	}
	929	if (target < 0.0) {
	930	if (numberX[9] <= 200)
	931	break; // not very many changes
	932	target = CoinMax(200.0, CoinMin(0.05 * numberX[9], 1000.0));
	933	}
	934	int iX = -1;
	935	int iBand = -1;
	936	for (i = 0; i < 10; i++) {
[753]	937	#ifdef COIN_DEVELOP_x
[1286]	938	printf("** %d changed at %g\n", numberX[i], factorX[i]);
[750]	939	#endif
[1286]	940	if (numberX[i] >= target && numberX[i] < 2.0*target && iX < 0)
	941	iX = i;
	942	if (iBand<0 && numberX[i]>target) {
	943	iBand = i;
	944	factor = factorX[i];
	945	}
	946	}
	947	if (iX >= 0) {
	948	factor = factorX[iX];
	949	break;
	950	} else {
	951	assert (iBand >= 0);
	952	double hi = factor;
	953	double lo = (iBand > 0) ? factorX[iBand-1] : 0.0;
	954	double diff = (hi - lo) / 9.0;
	955	for (i = 0; i < 10; i++) {
	956	factorX[i] = lo;
	957	lo += diff;
	958	}
	959	}
	960	}
	961	for (i = 0; i < numberIntegers; i++) {
	962	int iColumn = integerVariable[i];
	963	double value = randomX[i];
	964	double difference = fabs(solution[iColumn] - newSolution[iColumn]);
	965	if (difference + value*factor > 0.5) {
	966	numberPerturbed++;
	967	if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
	968	newSolution[iColumn] += 1.0;
	969	} else if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
	970	newSolution[iColumn] -= 1.0;
	971	} else {
	972	// general integer
	973	if (difference + value > 0.75)
	974	newSolution[iColumn] += 1.0;
	975	else
	976	newSolution[iColumn] -= 1.0;
	977	}
	978	}
	979	}
	980	delete [] randomX;
	981	} else {
	982	for (j = NUMBER_OLD - 1; j > 0; j--) {
	983	for (i = 0; i < numberColumns; i++) oldSolution[j][i] = oldSolution[j-1][i];
	984	}
	985	for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
	986	}
	987
	988	// 2. update the objective function based on the new rounded solution
	989	double offset = 0.0;
	990	double costValue = (1.0 - scaleFactor) * solver->getObjSense();
	991	int numberChanged = 0;
	992	const double * oldObjective = solver->getObjCoefficients();
	993	for (i = 0; i < numberColumns; i++) {
	994	// below so we can keep original code and allow for objective
	995	int iColumn = i;
	996	// Special code for "artificials"
	997	if (direction*saveObjective[iColumn] >= artificialCost_) {
	998	//solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
	999	solver->setObjCoeff(iColumn, (artificialFactor*saveObjective[iColumn]) / artificialCost_);
	1000	}
	1001	if (!solver->isBinary(iColumn) && !doGeneral)
	1002	continue;
	1003	// deal with fixed variables (i.e., upper=lower)
	1004	if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance \|\| !solver->isInteger(iColumn)) {
	1005	//if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
	1006	//else solver->setObjCoeff(iColumn,costValue);
	1007	continue;
	1008	}
	1009	double newValue = 0.0;
	1010	if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
	1011	newValue = costValue + scaleFactor * saveObjective[iColumn];
	1012	} else {
	1013	if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
	1014	newValue = -costValue + scaleFactor * saveObjective[iColumn];
	1015	}
	1016	}
[1271]	1017	#ifdef RAND_RAND
[1286]	1018	if (!offRandom)
	1019	newValue *= randomFactor[iColumn];
[1271]	1020	#endif
[1286]	1021	if (newValue != oldObjective[iColumn]) {
	1022	numberChanged++;
	1023	}
	1024	solver->setObjCoeff(iColumn, newValue);
	1025	offset += costValue * newSolution[iColumn];
	1026	}
[1883]	1027	if (numberPasses==1 && !totalNumberPasses && (model_->specialOptions()&8388608)!=0) {
	1028	// doing multiple solvers - make a real difference - flip 5%
	1029	for (i = 0; i < numberIntegers; i++) {
	1030	int iColumn = integerVariable[i];
	1031	double value = floor(newSolution[iColumn]+0.5);
	1032	if (fabs(value-solution[iColumn])>primalTolerance) {
	1033	value = randomNumberGenerator_.randomDouble();
	1034	if(value<0.05) {
	1035	//printf("Flipping %d - random %g\n",iColumn,value);
	1036	solver->setObjCoeff(iColumn,-solver->getObjCoefficients()[iColumn]);
	1037	}
	1038	}
	1039	}
	1040	}
[1286]	1041	solver->setDblParam(OsiObjOffset, -offset);
	1042	if (!general && false) {
	1043	// Solve in two goes - first keep satisfied ones fixed
	1044	double * saveLower = new double [numberIntegers];
	1045	double * saveUpper = new double [numberIntegers];
	1046	for (i = 0; i < numberIntegers; i++) {
	1047	int iColumn = integerVariable[i];
	1048	saveLower[i] = COIN_DBL_MAX;
	1049	saveUpper[i] = -COIN_DBL_MAX;
	1050	if (solution[iColumn] < lower[iColumn] + primalTolerance) {
	1051	saveUpper[i] = upper[iColumn];
	1052	solver->setColUpper(iColumn, lower[iColumn]);
	1053	} else if (solution[iColumn] > upper[iColumn] - primalTolerance) {
	1054	saveLower[i] = lower[iColumn];
	1055	solver->setColLower(iColumn, upper[iColumn]);
	1056	}
	1057	}
	1058	solver->resolve();
	1059	if (!solver->isProvenOptimal()) {
	1060	// presumably max time or some such
	1061	exitAll = true;
	1062	break;
	1063	}
	1064	for (i = 0; i < numberIntegers; i++) {
	1065	int iColumn = integerVariable[i];
	1066	if (saveLower[i] != COIN_DBL_MAX)
	1067	solver->setColLower(iColumn, saveLower[i]);
	1068	if (saveUpper[i] != -COIN_DBL_MAX)
	1069	solver->setColUpper(iColumn, saveUpper[i]);
	1070	saveUpper[i] = -COIN_DBL_MAX;
	1071	}
	1072	memcpy(newSolution, solution, numberColumns*sizeof(double));
	1073	int flip;
	1074	returnCode = rounds(solver, newSolution,/saveObjective,/
	1075	numberIntegers, integerVariable,
	1076	/pumpPrint,/numberPasses,
	1077	/roundExpensive_,/defaultRounding_, &flip);
	1078	numberPasses++;
	1079	if (returnCode) {
	1080	// solution - but may not be better
	1081	// Compute using dot product
	1082	double newSolutionValue = -saveOffset;
	1083	for ( i = 0 ; i < numberColumns ; i++ )
	1084	newSolutionValue += saveObjective[i] * newSolution[i];
	1085	newSolutionValue *= direction;
	1086	sprintf(pumpPrint, "Intermediate solution found of %g", newSolutionValue);
	1087	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1088	<< pumpPrint
	1089	<< CoinMessageEol;
	1090	if (newSolutionValue < solutionValue) {
	1091	memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
	1092	CoinWarmStartBasis * basis =
	1093	dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
	1094	solutionFound = true;
[1802]	1095	numberSolutions++;
	1096	if (numberSolutions>=maxSolutions)
	1097	exitAll = true;
[1286]	1098	if (exitNow(newSolutionValue))
	1099	exitAll = true;
	1100	if (basis) {
	1101	bestBasis = * basis;
	1102	delete basis;
	1103	int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
	1104	if (!action) {
	1105	double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
	1106	double saveObjectiveValue = model_->getMinimizationObjValue();
	1107	model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
	1108	if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
	1109	model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
	1110	delete [] saveOldSolution;
	1111	}
	1112	if (!action \|\| model_->maximumSecondsReached()) {
	1113	exitAll = true; // exit
	1114	break;
	1115	}
	1116	}
	1117	if ((accumulate_&1) != 0) {
	1118	model_->incrementUsed(betterSolution); // for local search
	1119	}
	1120	solutionValue = newSolutionValue;
	1121	solutionFound = true;
[1802]	1122	numberSolutions++;
	1123	if (numberSolutions>=maxSolutions)
	1124	exitAll = true;
[1286]	1125	if (exitNow(newSolutionValue))
	1126	exitAll = true;
	1127	} else {
	1128	returnCode = 0;
	1129	}
	1130	}
	1131	}
	1132	int numberIterations = 0;
	1133	if (!doGeneral) {
	1134	// faster to do from all slack!!!!
	1135	if (allSlack) {
	1136	CoinWarmStartBasis dummy;
	1137	solver->setWarmStart(&dummy);
	1138	}
[750]	1139	#ifdef COIN_DEVELOP
[1286]	1140	printf("%d perturbed out of %d columns (%d changed)\n", numberPerturbed, numberColumns, numberChanged);
[750]	1141	#endif
[1286]	1142	bool takeHint;
	1143	OsiHintStrength strength;
	1144	solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
	1145	if (dualPass && numberChanged > 2) {
	1146	solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
	1147	if (dualPass == 1 && 2*numberChanged < numberColumns &&
	1148	(numberChanged < 5000 \|\| 6*numberChanged < numberColumns)) {
	1149	// but we need to make infeasible
	1150	CoinWarmStartBasis * basis =
	1151	dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
	1152	if (basis) {
	1153	// modify
	1154	const double * lower = solver->getColLower();
	1155	const double * upper = solver->getColUpper();
	1156	double * solution = CoinCopyOfArray(solver->getColSolution(),
	1157	numberColumns);
	1158	const double * objective = solver->getObjCoefficients();
	1159	int nChanged = 0;
	1160	for (i = 0; i < numberIntegersOrig; i++) {
	1161	int iColumn = integerVariableOrig[i];
[1271]	1162	#ifdef RAND_RAND
[1286]	1163	if (nChanged > numberChanged)
	1164	break;
[1271]	1165	#endif
[1286]	1166	if (objective[iColumn] > 0.0) {
	1167	if (basis->getStructStatus(iColumn) ==
	1168	CoinWarmStartBasis::atUpperBound) {
	1169	solution[iColumn] = lower[iColumn];
	1170	basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
	1171	nChanged++;
	1172	}
	1173	} else if (objective[iColumn] < 0.0) {
	1174	if (basis->getStructStatus(iColumn) ==
	1175	CoinWarmStartBasis::atLowerBound) {
	1176	solution[iColumn] = upper[iColumn];
	1177	basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
	1178	nChanged++;
	1179	}
	1180	}
	1181	}
	1182	if (!nChanged) {
	1183	for (i = 0; i < numberIntegersOrig; i++) {
	1184	int iColumn = integerVariableOrig[i];
	1185	if (objective[iColumn] > 0.0) {
	1186	if (basis->getStructStatus(iColumn) ==
	1187	CoinWarmStartBasis::basic) {
	1188	solution[iColumn] = lower[iColumn];
	1189	basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
	1190	break;
	1191	}
	1192	} else if (objective[iColumn] < 0.0) {
	1193	if (basis->getStructStatus(iColumn) ==
	1194	CoinWarmStartBasis::basic) {
	1195	solution[iColumn] = upper[iColumn];
	1196	basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
	1197	break;
	1198	}
	1199	}
	1200	}
	1201	}
	1202	solver->setColSolution(solution);
	1203	delete [] solution;
	1204	solver->setWarmStart(basis);
	1205	delete basis;
	1206	}
	1207	} else {
	1208	// faster to do from all slack!!!! ???
	1209	CoinWarmStartBasis dummy;
	1210	solver->setWarmStart(&dummy);
	1211	}
	1212	}
	1213	if (numberTries > 1 && numberPasses == 1 && firstPerturbedObjective) {
	1214	// Modify to use convex combination
	1215	// use basis from first time
	1216	solver->setWarmStart(&saveBasis);
	1217	// and objective
	1218	if (secondPassOpt < 3 \|\| (secondPassOpt >= 4 && secondPassOpt < 6))
	1219	solver->setObjective(firstPerturbedObjective);
	1220	// and solution
	1221	solver->setColSolution(firstPerturbedSolution);
	1222	//if (secondPassOpt==2\|\|secondPassOpt==5\|\|
	1223	if (firstPerturbedValue > cutoff)
	1224	solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
	1225	}
	1226	solver->resolve();
	1227	if (!solver->isProvenOptimal()) {
	1228	// presumably max time or some such
	1229	exitAll = true;
	1230	break;
	1231	}
	1232	solver->setHintParam(OsiDoDualInResolve, takeHint);
	1233	newTrueSolutionValue = -saveOffset;
	1234	newSumInfeas = 0.0;
	1235	newNumberInfeas = 0;
	1236	{
	1237	const double * newSolution = solver->getColSolution();
	1238	for ( i = 0 ; i < numberColumns ; i++ ) {
	1239	if (solver->isInteger(i)) {
	1240	double value = newSolution[i];
	1241	double nearest = floor(value + 0.5);
	1242	newSumInfeas += fabs(value - nearest);
	1243	if (fabs(value - nearest) > 1.0e-6)
	1244	newNumberInfeas++;
	1245	}
	1246	newTrueSolutionValue += saveObjective[i] * newSolution[i];
	1247	}
	1248	newTrueSolutionValue *= direction;
	1249	if (numberPasses == 1 && secondPassOpt) {
	1250	if (numberTries == 1 \|\| secondPassOpt > 3) {
	1251	// save basis
	1252	CoinWarmStartBasis * basis =
	1253	dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
	1254	if (basis) {
	1255	saveBasis = * basis;
	1256	delete basis;
	1257	}
	1258	delete [] firstPerturbedObjective;
	1259	delete [] firstPerturbedSolution;
	1260	firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
	1261	firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(), numberColumns);
	1262	firstPerturbedValue = newTrueSolutionValue;
	1263	}
	1264	}
	1265	if (newNumberInfeas && newNumberInfeas < 15) {
[1393]	1266	#ifdef JJF_ZERO
[1286]	1267	roundingObjective = solutionValue;
	1268	OsiSolverInterface * saveSolver = model_->swapSolver(solver);
	1269	double * currentObjective =
	1270	CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
	1271	solver->setObjective(saveObjective);
	1272	double saveOffset2;
	1273	solver->getDblParam(OsiObjOffset, saveOffset2);
	1274	solver->setDblParam(OsiObjOffset, saveOffset);
	1275	int ifSol = roundingHeuristic.solution(roundingObjective, roundingSolution);
	1276	solver->setObjective(currentObjective);
	1277	solver->setDblParam(OsiObjOffset, saveOffset2);
	1278	delete [] currentObjective;
	1279	model_->swapSolver(saveSolver);
	1280	if (ifSol > 0)
	1281	abort();
[1013]	1282	#endif
[1286]	1283	int numberRows = solver->getNumRows();
	1284	double * rowActivity = new double[numberRows];
	1285	memset(rowActivity, 0, numberRows*sizeof(double));
	1286	int * which = new int[newNumberInfeas];
	1287	int * stack = new int[newNumberInfeas+1];
	1288	double * baseValue = new double[newNumberInfeas];
	1289	int * whichRow = new int[numberRows];
	1290	double * rowValue = new double[numberRows];
	1291	memset(rowValue, 0, numberRows*sizeof(double));
	1292	int nRow = 0;
	1293	// Column copy
	1294	const double * element = solver->getMatrixByCol()->getElements();
	1295	const int * row = solver->getMatrixByCol()->getIndices();
	1296	const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
	1297	const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
	1298	int n = 0;
	1299	double contrib = 0.0;
	1300	for ( i = 0 ; i < numberColumns ; i++ ) {
	1301	double value = newSolution[i];
	1302	if (solver->isInteger(i)) {
	1303	double nearest = floor(value + 0.5);
	1304	if (fabs(value - nearest) > 1.0e-6) {
	1305	//printf("Column %d value %g\n",i,value);
	1306	for (CoinBigIndex j = columnStart[i];
	1307	j < columnStart[i] + columnLength[i]; j++) {
	1308	int iRow = row[j];
	1309	//printf("row %d element %g\n",iRow,element[j]);
	1310	if (!rowValue[iRow]) {
	1311	rowValue[iRow] = 1.0;
	1312	whichRow[nRow++] = iRow;
	1313	}
	1314	}
	1315	baseValue[n] = floor(value);
	1316	contrib += saveObjective[i] * value;
	1317	value = 0.0;
	1318	stack[n] = 0;
	1319	which[n++] = i;
	1320	}
	1321	}
	1322	for (CoinBigIndex j = columnStart[i];
	1323	j < columnStart[i] + columnLength[i]; j++) {
	1324	int iRow = row[j];
	1325	rowActivity[iRow] += value * element[j];
	1326	}
	1327	}
	1328	if (newNumberInfeas < 15) {
	1329	stack[n] = newNumberInfeas + 100;
	1330	int iStack = n;
	1331	memset(rowValue, 0, numberRows*sizeof(double));
	1332	const double * rowLower = solver->getRowLower();
	1333	const double * rowUpper = solver->getRowUpper();
	1334	while (iStack >= 0) {
	1335	double contrib2 = 0.0;
	1336	// Could do faster
	1337	for (int k = 0 ; k < n ; k++ ) {
	1338	i = which[k];
	1339	double value = baseValue[k] + stack[k];
	1340	contrib2 += saveObjective[i] * value;
	1341	for (CoinBigIndex j = columnStart[i];
	1342	j < columnStart[i] + columnLength[i]; j++) {
	1343	int iRow = row[j];
	1344	rowValue[iRow] += value * element[j];
	1345	}
	1346	}
	1347	// check if feasible
	1348	bool feasible = true;
	1349	for (int k = 0; k < nRow; k++) {
	1350	i = whichRow[k];
	1351	double value = rowValue[i] + rowActivity[i];
	1352	rowValue[i] = 0.0;
	1353	if (value < rowLower[i] - 1.0e-7 \|\|
	1354	value > rowUpper[i] + 1.0e-7)
	1355	feasible = false;
	1356	}
	1357	if (feasible) {
	1358	double newObj = newTrueSolutionValue * direction;
	1359	newObj += contrib2 - contrib;
	1360	newObj *= direction;
[1127]	1361	#ifdef COIN_DEVELOP
[1286]	1362	printf("FFFeasible! - obj %g\n", newObj);
[1127]	1363	#endif
[1286]	1364	if (newObj < roundingObjective - 1.0e-6) {
[1127]	1365	#ifdef COIN_DEVELOP
[1286]	1366	printf("FBetter\n");
[1127]	1367	#endif
[1286]	1368	roundingObjective = newObj;
	1369	memcpy(roundingSolution, newSolution, numberColumns*sizeof(double));
	1370	for (int k = 0 ; k < n ; k++ ) {
	1371	i = which[k];
	1372	double value = baseValue[k] + stack[k];
	1373	roundingSolution[i] = value;
	1374	}
	1375	}
	1376	}
	1377	while (iStack >= 0 && stack[iStack]) {
	1378	stack[iStack]--;
	1379	iStack--;
	1380	}
	1381	if (iStack >= 0) {
	1382	stack[iStack] = 1;
	1383	iStack = n;
	1384	stack[n] = 1;
	1385	}
	1386	}
	1387	}
	1388	delete [] rowActivity;
	1389	delete [] which;
	1390	delete [] stack;
	1391	delete [] baseValue;
	1392	delete [] whichRow;
	1393	delete [] rowValue;
	1394	}
	1395	}
	1396	if (true) {
	1397	OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
	1398	clonedSolver->setColSolution(solver->getColSolution());
	1399	CbcRounding heuristic1(*model_);
	1400	heuristic1.setHeuristicName("rounding in feaspump!");
	1401	heuristic1.setWhen(1);
	1402	roundingObjective = CoinMin(roundingObjective, solutionValue);
	1403	double testSolutionValue = newTrueSolutionValue;
	1404	int returnCode = heuristic1.solution(roundingObjective,
	1405	roundingSolution,
	1406	testSolutionValue) ;
	1407	if (returnCode == 1) {
[1127]	1408	#ifdef COIN_DEVELOP
[1286]	1409	printf("rounding obj of %g?\n", roundingObjective);
[1127]	1410	#endif
[1286]	1411	//roundingObjective = newSolutionValue;
	1412	//} else {
	1413	//roundingObjective = COIN_DBL_MAX;
	1414	}
	1415	model_->swapSolver(saveSolver);
	1416	}
	1417	if (!solver->isProvenOptimal()) {
	1418	// presumably max time or some such
	1419	exitAll = true;
	1420	break;
	1421	}
	1422	// in case very dubious solver
	1423	lower = solver->getColLower();
	1424	upper = solver->getColUpper();
	1425	solution = solver->getColSolution();
	1426	numberIterations = solver->getIterationCount();
	1427	} else {
	1428	int * addStart = new int[2*general+1];
	1429	int * addIndex = new int[4*general];
	1430	double * addElement = new double[4*general];
	1431	double * addLower = new double[2*general];
	1432	double * addUpper = new double[2*general];
	1433	double * obj = new double[general];
	1434	int nAdd = 0;
	1435	for (i = 0; i < numberIntegers; i++) {
	1436	int iColumn = integerVariable[i];
	1437	if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
	1438	newSolution[iColumn] < upper[iColumn] - primalTolerance) {
	1439	assert (upper[iColumn] - lower[iColumn] > 1.00001);
	1440	obj[nAdd] = 1.0;
	1441	addLower[nAdd] = 0.0;
	1442	addUpper[nAdd] = COIN_DBL_MAX;
	1443	nAdd++;
	1444	}
	1445	}
	1446	OsiSolverInterface * solver2 = solver;
	1447	if (nAdd) {
	1448	CoinZeroN(addStart, nAdd + 1);
	1449	solver2 = solver->clone();
	1450	solver2->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, obj);
	1451	// feasible solution
	1452	double * sol = new double[nAdd+numberColumns];
	1453	memcpy(sol, solution, numberColumns*sizeof(double));
	1454	// now rows
	1455	int nAdd = 0;
	1456	int nEl = 0;
	1457	int nAddRow = 0;
	1458	for (i = 0; i < numberIntegers; i++) {
	1459	int iColumn = integerVariable[i];
	1460	if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
	1461	newSolution[iColumn] < upper[iColumn] - primalTolerance) {
	1462	addLower[nAddRow] = -newSolution[iColumn];;
	1463	addUpper[nAddRow] = COIN_DBL_MAX;
	1464	addIndex[nEl] = iColumn;
	1465	addElement[nEl++] = -1.0;
	1466	addIndex[nEl] = numberColumns + nAdd;
	1467	addElement[nEl++] = 1.0;
	1468	nAddRow++;
	1469	addStart[nAddRow] = nEl;
	1470	addLower[nAddRow] = newSolution[iColumn];;
	1471	addUpper[nAddRow] = COIN_DBL_MAX;
	1472	addIndex[nEl] = iColumn;
	1473	addElement[nEl++] = 1.0;
	1474	addIndex[nEl] = numberColumns + nAdd;
	1475	addElement[nEl++] = 1.0;
	1476	nAddRow++;
	1477	addStart[nAddRow] = nEl;
	1478	sol[nAdd+numberColumns] = fabs(sol[iColumn] - newSolution[iColumn]);
	1479	nAdd++;
	1480	}
	1481	}
	1482	solver2->setColSolution(sol);
	1483	delete [] sol;
	1484	solver2->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
	1485	}
	1486	delete [] addStart;
	1487	delete [] addIndex;
	1488	delete [] addElement;
	1489	delete [] addLower;
	1490	delete [] addUpper;
	1491	delete [] obj;
	1492	solver2->resolve();
	1493	if (!solver2->isProvenOptimal()) {
	1494	// presumably max time or some such
	1495	exitAll = true;
	1496	break;
	1497	}
	1498	//assert (solver2->isProvenOptimal());
	1499	if (nAdd) {
	1500	solver->setColSolution(solver2->getColSolution());
	1501	numberIterations = solver2->getIterationCount();
	1502	delete solver2;
	1503	} else {
	1504	numberIterations = solver->getIterationCount();
	1505	}
	1506	lower = solver->getColLower();
	1507	upper = solver->getColUpper();
	1508	solution = solver->getColSolution();
	1509	newTrueSolutionValue = -saveOffset;
	1510	newSumInfeas = 0.0;
	1511	newNumberInfeas = 0;
	1512	{
	1513	const double * newSolution = solver->getColSolution();
	1514	for ( i = 0 ; i < numberColumns ; i++ ) {
	1515	if (solver->isInteger(i)) {
	1516	double value = newSolution[i];
	1517	double nearest = floor(value + 0.5);
	1518	newSumInfeas += fabs(value - nearest);
	1519	if (fabs(value - nearest) > 1.0e-6)
	1520	newNumberInfeas++;
	1521	}
	1522	newTrueSolutionValue += saveObjective[i] * newSolution[i];
	1523	}
	1524	newTrueSolutionValue *= direction;
	1525	}
	1526	}
	1527	if (lastMove != 1000000) {
	1528	if (newSumInfeas < lastSumInfeas) {
	1529	lastMove = numberPasses;
	1530	lastSumInfeas = newSumInfeas;
	1531	} else if (newSumInfeas > lastSumInfeas + 1.0e-5) {
	1532	lastMove = 1000000; // going up
	1533	}
	1534	}
	1535	totalNumberIterations += numberIterations;
	1536	if (solver->getNumRows() < 3000)
	1537	sprintf(pumpPrint, "Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d",
	1538	numberPasses + totalNumberPasses,
	1539	newSumInfeas, newNumberInfeas,
	1540	newTrueSolutionValue, numberIterations);
	1541	else
	1542	sprintf(pumpPrint, "Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses + totalNumberPasses,
	1543	model_->getCurrentSeconds(), newSumInfeas, newNumberInfeas,
	1544	newTrueSolutionValue, numberIterations);
	1545	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1546	<< pumpPrint
	1547	<< CoinMessageEol;
	1548	if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
	1549	int i;
	1550	const double * objective = solver->getObjCoefficients();
	1551	for (i = 0; i < numberIntegersOrig; i++) {
	1552	int iColumn = integerVariableOrig[i];
	1553	if (objective[iColumn] > 0.0)
	1554	closestSolution[i] = 0;
	1555	else
	1556	closestSolution[i] = 1;
	1557	}
	1558	closestObjectiveValue = solver->getObjValue();
	1559	}
	1560	// See if we need to think about changing rhs
	1561	if ((switches_&12) != 0 && useRhs < 1.0e50) {
	1562	double oldRhs = useRhs;
	1563	bool trying = false;
	1564	if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) {
	1565	if (solutionValue > 1.0e20) {
	1566	// only if no genuine solution
	1567	double gap = useRhs - continuousObjectiveValue;
	1568	useRhs += 0.1 * gap;
	1569	if (exactMultiple) {
	1570	useRhs = exactMultiple * ceil(useRhs / exactMultiple);
	1571	useRhs = CoinMax(useRhs, oldRhs + exactMultiple);
	1572	}
	1573	trying = true;
	1574	}
	1575	}
	1576	if ((switches_&8) != 0) {
	1577	// Put in new suminf and check
	1578	double largest = newSumInfeas;
	1579	double smallest = newSumInfeas;
	1580	for (int i = 0; i < SIZE_BOBBLE - 1; i++) {
	1581	double value = saveSumInf[i+1];
	1582	saveSumInf[i] = value;
	1583	largest = CoinMax(largest, value);
	1584	smallest = CoinMin(smallest, value);
	1585	}
	1586	saveSumInf[SIZE_BOBBLE-1] = newSumInfeas;
	1587	if (smallest*1.5 > largest && smallest > 2.0) {
	1588	if (bobbleMode == 0) {
	1589	// go closer
	1590	double gap = oldRhs - continuousObjectiveValue;
	1591	useRhs -= 0.4 * gap;
	1592	if (exactMultiple) {
	1593	double value = floor(useRhs / exactMultiple);
	1594	useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple);
	1595	}
	1596	if (useRhs < continuousObjectiveValue) {
	1597	// skip decrease
	1598	bobbleMode = 1;
	1599	useRhs = oldRhs;
	1600	}
	1601	}
	1602	if (bobbleMode) {
	1603	trying = true;
	1604	// weaken
	1605	if (solutionValue < 1.0e20) {
	1606	double gap = solutionValue - oldRhs;
	1607	useRhs += 0.3 * gap;
	1608	} else {
	1609	double gap = oldRhs - continuousObjectiveValue;
	1610	useRhs += 0.05 * gap;
	1611	}
	1612	if (exactMultiple) {
	1613	double value = ceil(useRhs / exactMultiple);
	1614	useRhs = CoinMin(value * exactMultiple,
	1615	solutionValue - exactMultiple);
	1616	}
	1617	}
	1618	bobbleMode++;
	1619	// reset
	1620	CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
	1621	}
	1622	}
	1623	if (useRhs != oldRhs) {
	1624	// tidy up
	1625	if (exactMultiple) {
	1626	double value = floor(useRhs / exactMultiple);
	1627	double bestPossible = ceil(continuousObjectiveValue / exactMultiple);
	1628	useRhs = CoinMax(value, bestPossible) * exactMultiple;
	1629	} else {
	1630	useRhs = CoinMax(useRhs, continuousObjectiveValue);
	1631	}
	1632	int k = solver->getNumRows() - 1;
	1633	solver->setRowUpper(k, useRhs + useOffset);
	1634	bool takeHint;
	1635	OsiHintStrength strength;
	1636	solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
	1637	if (useRhs < oldRhs) {
	1638	solver->setHintParam(OsiDoDualInResolve, true);
	1639	solver->resolve();
	1640	} else if (useRhs > oldRhs) {
	1641	solver->setHintParam(OsiDoDualInResolve, false);
	1642	solver->resolve();
	1643	}
	1644	solver->setHintParam(OsiDoDualInResolve, takeHint);
	1645	if (!solver->isProvenOptimal()) {
	1646	// presumably max time or some such
	1647	exitAll = true;
	1648	break;
	1649	}
	1650	} else if (trying) {
	1651	// doesn't look good
	1652	break;
	1653	}
	1654	}
	1655	}
	1656	// reduce scale factor
	1657	scaleFactor *= weightFactor_;
	1658	} // END WHILE
	1659	// see if rounding worked!
	1660	if (roundingObjective < solutionValue) {
	1661	if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) {
	1662	sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n",
	1663	roundingObjective, solutionValue);
	1664	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1665	<< pumpPrint
	1666	<< CoinMessageEol;
	1667	}
	1668	solutionValue = roundingObjective;
	1669	newSolutionValue = solutionValue;
	1670	memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
	1671	solutionFound = true;
[1802]	1672	numberSolutions++;
	1673	if (numberSolutions>=maxSolutions)
	1674	exitAll = true;
[1286]	1675	if (exitNow(roundingObjective))
	1676	exitAll = true;
	1677	}
	1678	if (!solutionFound) {
	1679	sprintf(pumpPrint, "No solution found this major pass");
	1680	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1681	<< pumpPrint
	1682	<< CoinMessageEol;
	1683	}
	1684	//}
	1685	delete solver;
	1686	solver = NULL;
	1687	for ( j = 0; j < NUMBER_OLD; j++)
	1688	delete [] oldSolution[j];
	1689	delete [] oldSolution;
	1690	delete [] saveObjective;
	1691	if (usedColumn && !exitAll) {
	1692	OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
[1839]	1693	#if 0 //def COIN_HAS_CLP
	1694	OsiClpSolverInterface * clpSolver
	1695	= dynamic_cast<OsiClpSolverInterface *> (newSolver);
	1696	if (clpSolver) {
	1697	ClpSimplex * simplex = clpSolver->getModelPtr();
	1698	simplex->writeMps("start.mps",2,1);
	1699	}
	1700	#endif
[1286]	1701	const double * colLower = newSolver->getColLower();
	1702	const double * colUpper = newSolver->getColUpper();
	1703	bool stopBAB = false;
	1704	int allowedPass = -1;
	1705	if (maximumAllowed > 0)
	1706	allowedPass = CoinMax(numberPasses - maximumAllowed, -1);
	1707	while (!stopBAB) {
	1708	stopBAB = true;
	1709	int i;
	1710	int nFix = 0;
	1711	int nFixI = 0;
	1712	int nFixC = 0;
	1713	int nFixC2 = 0;
	1714	for (i = 0; i < numberIntegersOrig; i++) {
	1715	int iColumn = integerVariableOrig[i];
	1716	//const OsiObject * object = model_->object(i);
	1717	//double originalLower;
	1718	//double originalUpper;
	1719	//getIntegerInformation( object,originalLower, originalUpper);
	1720	//assert(colLower[iColumn]==originalLower);
	1721	//newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
	1722	newSolver->setColLower(iColumn, colLower[iColumn]);
	1723	//assert(colUpper[iColumn]==originalUpper);
	1724	//newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
	1725	newSolver->setColUpper(iColumn, colUpper[iColumn]);
	1726	if (usedColumn[iColumn] <= allowedPass) {
	1727	double value = lastSolution[iColumn];
	1728	double nearest = floor(value + 0.5);
	1729	if (fabs(value - nearest) < 1.0e-7) {
	1730	if (nearest == colLower[iColumn]) {
	1731	newSolver->setColUpper(iColumn, colLower[iColumn]);
	1732	nFix++;
	1733	} else if (nearest == colUpper[iColumn]) {
	1734	newSolver->setColLower(iColumn, colUpper[iColumn]);
	1735	nFix++;
	1736	} else if (fixInternal) {
	1737	newSolver->setColLower(iColumn, nearest);
	1738	newSolver->setColUpper(iColumn, nearest);
	1739	nFix++;
	1740	nFixI++;
	1741	}
	1742	}
	1743	}
	1744	}
	1745	if (fixContinuous) {
	1746	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
	1747	if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
	1748	double value = lastSolution[iColumn];
	1749	if (value < colLower[iColumn] + 1.0e-8) {
	1750	newSolver->setColUpper(iColumn, colLower[iColumn]);
	1751	nFixC++;
	1752	} else if (value > colUpper[iColumn] - 1.0e-8) {
	1753	newSolver->setColLower(iColumn, colUpper[iColumn]);
	1754	nFixC++;
	1755	} else if (fixContinuous == 2) {
	1756	newSolver->setColLower(iColumn, value);
	1757	newSolver->setColUpper(iColumn, value);
	1758	nFixC++;
	1759	nFixC2++;
	1760	}
	1761	}
	1762	}
	1763	}
	1764	newSolver->initialSolve();
	1765	if (!newSolver->isProvenOptimal()) {
	1766	//newSolver->writeMps("bad.mps");
	1767	//assert (newSolver->isProvenOptimal());
	1768	exitAll = true;
	1769	break;
	1770	}
	1771	sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous",
	1772	nFix, nFixC);
	1773	if (nFixC2 + nFixI != 0)
	1774	sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous",
	1775	nFixI, nFixC2);
	1776	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1777	<< pumpPrint
	1778	<< CoinMessageEol;
	1779	double saveValue = newSolutionValue;
	1780	if (newSolutionValue - model_->getCutoffIncrement()
	1781	> continuousObjectiveValue - 1.0e-7) {
	1782	double saveFraction = fractionSmall_;
	1783	if (numberTries > 1 && !numberBandBsolutions)
	1784	fractionSmall_ *= 0.5;
	1785	// Give branch and bound a bit more freedom
	1786	double cutoff2 = newSolutionValue +
	1787	CoinMax(model_->getCutoffIncrement(), 1.0e-3);
[1839]	1788	#if 0
	1789	{
	1790	OsiClpSolverInterface * clpSolver
	1791	= dynamic_cast<OsiClpSolverInterface *> (newSolver);
	1792	if (clpSolver) {
	1793	ClpSimplex * simplex = clpSolver->getModelPtr();
	1794	simplex->writeMps("testA.mps",2,1);
	1795	}
	1796	}
	1797	#endif
[1286]	1798	int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
	1799	cutoff2, "CbcHeuristicLocalAfterFPump");
	1800	fractionSmall_ = saveFraction;
	1801	if (returnCode2 < 0) {
	1802	if (returnCode2 == -2) {
	1803	exitAll = true;
	1804	returnCode = 0;
	1805	} else {
	1806	returnCode2 = 0; // returned on size - try changing
	1807	//#define ROUND_AGAIN
[1053]	1808	#ifdef ROUND_AGAIN
[1286]	1809	if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
	1810	allowedPass = (numberPasses + allowedPass) >> 1;
	1811	sprintf(pumpPrint,
	1812	"Fixing all variables which were last changed on pass %d and trying again",
	1813	allowedPass);
	1814	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1815	<< pumpPrint
	1816	<< CoinMessageEol;
	1817	stopBAB = false;
	1818	continue;
	1819	}
[1053]	1820	#endif
[1286]	1821	}
	1822	}
	1823	if ((returnCode2&2) != 0) {
	1824	// could add cut
	1825	returnCode2 &= ~2;
	1826	}
[1870]	1827	if (returnCode2) {
[1286]	1828	numberBandBsolutions++;
[1870]	1829	// may not have got solution earlier
	1830	returnCode \|= 1;
	1831	}
[1286]	1832	} else {
	1833	// no need
	1834	exitAll = true;
	1835	//returnCode=0;
	1836	}
	1837	// recompute solution value
	1838	if (returnCode && true) {
[1839]	1839	#if 0
	1840	{
	1841	OsiClpSolverInterface * clpSolver
	1842	= dynamic_cast<OsiClpSolverInterface *> (newSolver);
	1843	if (clpSolver) {
	1844	ClpSimplex * simplex = clpSolver->getModelPtr();
	1845	simplex->writeMps("testB.mps",2,1);
	1846	}
	1847	}
	1848	#endif
[1286]	1849	delete newSolver;
	1850	newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
	1851	newSolutionValue = -saveOffset;
	1852	double newSumInfeas = 0.0;
	1853	const double * obj = newSolver->getObjCoefficients();
	1854	for (int i = 0 ; i < numberColumns ; i++ ) {
	1855	if (newSolver->isInteger(i)) {
	1856	double value = newSolution[i];
	1857	double nearest = floor(value + 0.5);
	1858	newSumInfeas += fabs(value - nearest);
	1859	}
	1860	newSolutionValue += obj[i] * newSolution[i];
	1861	}
	1862	newSolutionValue *= direction;
	1863	}
	1864	bool gotSolution = false;
	1865	if (returnCode && newSolutionValue < saveValue) {
	1866	sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
	1867	saveValue, newSolutionValue, model_->getCurrentSeconds());
	1868	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1869	<< pumpPrint
	1870	<< CoinMessageEol;
	1871	memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
	1872	gotSolution = true;
	1873	if (fixContinuous && nFixC + nFixC2 > 0) {
	1874	// may be able to do even better
	1875	int nFixed = 0;
	1876	const double * lower = model_->solver()->getColLower();
	1877	const double * upper = model_->solver()->getColUpper();
	1878	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
	1879	double value = newSolution[iColumn];
	1880	if (newSolver->isInteger(iColumn)) {
	1881	value = floor(newSolution[iColumn] + 0.5);
	1882	newSolver->setColLower(iColumn, value);
	1883	newSolver->setColUpper(iColumn, value);
	1884	nFixed++;
	1885	} else {
	1886	newSolver->setColLower(iColumn, lower[iColumn]);
	1887	newSolver->setColUpper(iColumn, upper[iColumn]);
	1888	if (value < lower[iColumn])
	1889	value = lower[iColumn];
	1890	else if (value > upper[iColumn])
	1891	value = upper[iColumn];
	1892
	1893	}
	1894	newSolution[iColumn] = value;
	1895	}
	1896	newSolver->setColSolution(newSolution);
	1897	//#define CLP_INVESTIGATE2
[1271]	1898	#ifdef CLP_INVESTIGATE2
[1286]	1899	{
	1900	// check
	1901	// get row activities
	1902	int numberRows = newSolver->getNumRows();
	1903	double * rowActivity = new double[numberRows];
	1904	memset(rowActivity, 0, numberRows*sizeof(double));
	1905	newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
	1906	double largestInfeasibility = primalTolerance;
	1907	double sumInfeasibility = 0.0;
	1908	int numberBadRows = 0;
	1909	const double * rowLower = newSolver->getRowLower();
	1910	const double * rowUpper = newSolver->getRowUpper();
	1911	for (i = 0 ; i < numberRows ; i++) {
	1912	double value;
	1913	value = rowLower[i] - rowActivity[i];
	1914	if (value > primalTolerance) {
	1915	numberBadRows++;
	1916	largestInfeasibility = CoinMax(largestInfeasibility, value);
	1917	sumInfeasibility += value;
	1918	}
	1919	value = rowActivity[i] - rowUpper[i];
	1920	if (value > primalTolerance) {
	1921	numberBadRows++;
	1922	largestInfeasibility = CoinMax(largestInfeasibility, value);
	1923	sumInfeasibility += value;
	1924	}
	1925	}
	1926	printf("%d bad rows, largest inf %g sum %g\n",
	1927	numberBadRows, largestInfeasibility, sumInfeasibility);
	1928	delete [] rowActivity;
	1929	}
[1271]	1930	#endif
	1931	#ifdef COIN_HAS_CLP
[1286]	1932	OsiClpSolverInterface * clpSolver
	1933	= dynamic_cast<OsiClpSolverInterface *> (newSolver);
	1934	if (clpSolver) {
	1935	ClpSimplex * simplex = clpSolver->getModelPtr();
	1936	//simplex->writeBasis("test.bas",true,2);
	1937	//simplex->writeMps("test.mps",2,1);
	1938	if (nFixed3 > numberColumns2)
	1939	simplex->allSlackBasis(); // may as well go from all slack
	1940	simplex->primal(1);
	1941	clpSolver->setWarmStart(NULL);
	1942	}
[1271]	1943	#endif
[1286]	1944	newSolver->initialSolve();
	1945	if (newSolver->isProvenOptimal()) {
	1946	double value = newSolver->getObjValue() * newSolver->getObjSense();
	1947	if (value < newSolutionValue) {
	1948	//newSolver->writeMps("query","mps");
[1393]	1949	#ifdef JJF_ZERO
[1286]	1950	{
	1951	double saveOffset;
	1952	newSolver->getDblParam(OsiObjOffset, saveOffset);
	1953	const double * obj = newSolver->getObjCoefficients();
	1954	double newTrueSolutionValue = -saveOffset;
	1955	double newSumInfeas = 0.0;
	1956	int numberColumns = newSolver->getNumCols();
	1957	const double * solution = newSolver->getColSolution();
	1958	for (int i = 0 ; i < numberColumns ; i++ ) {
	1959	if (newSolver->isInteger(i)) {
	1960	double value = solution[i];
	1961	double nearest = floor(value + 0.5);
	1962	newSumInfeas += fabs(value - nearest);
	1963	}
	1964	if (solution[i])
	1965	printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
	1966	newTrueSolutionValue);
	1967	newTrueSolutionValue += obj[i] * solution[i];
	1968	}
	1969	printf("obj %g - inf %g\n", newTrueSolutionValue,
	1970	newSumInfeas);
	1971	}
[1053]	1972	#endif
[1286]	1973	sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
	1974	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1975	<< pumpPrint
	1976	<< CoinMessageEol;
	1977	newSolutionValue = value;
	1978	memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
	1979	}
	1980	} else {
	1981	//newSolver->writeMps("bad3.mps");
[1839]	1982	sprintf(pumpPrint, "On closer inspection solution is not valid");
	1983	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1984	<< pumpPrint
	1985	<< CoinMessageEol;
[1286]	1986	exitAll = true;
	1987	break;
	1988	}
	1989	}
	1990	} else {
	1991	sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
	1992	model_->getCurrentSeconds());
	1993	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	1994	<< pumpPrint
	1995	<< CoinMessageEol;
	1996	if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
	1997	// may be able to do better
	1998	const double * lower = model_->solver()->getColLower();
	1999	const double * upper = model_->solver()->getColUpper();
	2000	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
	2001	if (newSolver->isInteger(iColumn)) {
	2002	double value = floor(newSolution[iColumn] + 0.5);
	2003	newSolver->setColLower(iColumn, value);
	2004	newSolver->setColUpper(iColumn, value);
	2005	} else {
	2006	newSolver->setColLower(iColumn, lower[iColumn]);
	2007	newSolver->setColUpper(iColumn, upper[iColumn]);
	2008	}
	2009	}
	2010	newSolver->initialSolve();
	2011	if (newSolver->isProvenOptimal()) {
	2012	double value = newSolver->getObjValue() * newSolver->getObjSense();
	2013	if (value < saveValue) {
	2014	sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
	2015	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	2016	<< pumpPrint
	2017	<< CoinMessageEol;
	2018	newSolutionValue = value;
	2019	memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
	2020	gotSolution = true;
	2021	}
	2022	}
	2023	}
	2024	}
	2025	if (gotSolution) {
	2026	if ((accumulate_&1) != 0) {
	2027	model_->incrementUsed(betterSolution); // for local search
	2028	}
	2029	solutionValue = newSolutionValue;
	2030	solutionFound = true;
[1802]	2031	numberSolutions++;
	2032	if (numberSolutions>=maxSolutions)
	2033	exitAll = true;
[1286]	2034	if (exitNow(newSolutionValue))
	2035	exitAll = true;
	2036	CoinWarmStartBasis * basis =
	2037	dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
	2038	if (basis) {
	2039	bestBasis = * basis;
	2040	delete basis;
	2041	int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
	2042	if (action == 0) {
	2043	double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
	2044	double saveObjectiveValue = model_->getMinimizationObjValue();
	2045	model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
	2046	if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
	2047	model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
	2048	delete [] saveOldSolution;
	2049	}
	2050	if (!action \|\| model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
	2051	exitAll = true; // exit
	2052	break;
	2053	}
	2054	}
	2055	}
	2056	} // end stopBAB while
	2057	delete newSolver;
	2058	}
	2059	if (solutionFound) finalReturnCode = 1;
	2060	cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
	2061	if (numberTries >= maximumRetries_ \|\| !solutionFound \|\| exitAll \|\| cutoff < continuousObjectiveValue + 1.0e-7) {
	2062	break;
	2063	} else {
	2064	solutionFound = false;
	2065	if (absoluteIncrement_ > 0.0 \|\| relativeIncrement_ > 0.0) {
	2066	double gap = relativeIncrement_ * fabs(solutionValue);
	2067	double change = CoinMax(gap, absoluteIncrement_);
	2068	cutoff = CoinMin(cutoff, solutionValue - change);
	2069	} else {
	2070	//double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
	2071	double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
	2072	cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
	2073	}
	2074	// But round down
	2075	if (exactMultiple)
	2076	cutoff = exactMultiple * floor(cutoff / exactMultiple);
	2077	if (cutoff < continuousObjectiveValue)
	2078	break;
	2079	sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
	2080	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	2081	<< pumpPrint
	2082	<< CoinMessageEol;
	2083	if ((accumulate_&3) < 2 && usedColumn) {
	2084	for (int i = 0; i < numberColumns; i++)
	2085	usedColumn[i] = -1;
	2086	}
	2087	averageIterationsPerTry = totalNumberIterations / numberTries;
	2088	numberIterationsLastPass = totalNumberIterations;
	2089	totalNumberPasses += numberPasses - 1;
	2090	}
[175]	2091	}
[1364]	2092	/*
	2093	End of the `exitAll' loop.
	2094	*/
[1271]	2095	#ifdef RAND_RAND
[1286]	2096	delete [] randomFactor;
[1271]	2097	#endif
[1286]	2098	delete solver; // probably NULL but do anyway
	2099	if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
	2100	usedColumn && !model_->maximumSecondsReached()) {
	2101	// try a bit of branch and bound
	2102	OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
	2103	const double * colLower = newSolver->getColLower();
	2104	const double * colUpper = newSolver->getColUpper();
	2105	int i;
	2106	double rhs = 0.0;
	2107	for (i = 0; i < numberIntegersOrig; i++) {
	2108	int iColumn = integerVariableOrig[i];
	2109	int direction = closestSolution[i];
	2110	closestSolution[i] = iColumn;
	2111	if (direction == 0) {
	2112	// keep close to LB
	2113	rhs += colLower[iColumn];
	2114	lastSolution[i] = 1.0;
	2115	} else {
	2116	// keep close to UB
	2117	rhs -= colUpper[iColumn];
	2118	lastSolution[i] = -1.0;
	2119	}
	2120	}
	2121	newSolver->addRow(numberIntegersOrig, closestSolution,
	2122	lastSolution, -COIN_DBL_MAX, rhs + 10.0);
	2123	//double saveValue = newSolutionValue;
	2124	//newSolver->writeMps("sub");
	2125	int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
	2126	newSolutionValue, "CbcHeuristicLocalAfterFPump");
	2127	if (returnCode < 0)
	2128	returnCode = 0; // returned on size
	2129	if ((returnCode&2) != 0) {
	2130	// could add cut
	2131	returnCode &= ~2;
	2132	}
	2133	if (returnCode) {
	2134	//printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
	2135	memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
	2136	//abort();
	2137	solutionValue = newSolutionValue;
	2138	solutionFound = true;
[1802]	2139	numberSolutions++;
	2140	if (numberSolutions>=maxSolutions)
	2141	exitAll = true;
[1286]	2142	if (exitNow(newSolutionValue))
	2143	exitAll = true;
	2144	}
	2145	delete newSolver;
[640]	2146	}
[1286]	2147	delete clonedSolver;
	2148	delete [] roundingSolution;
	2149	delete [] usedColumn;
	2150	delete [] lastSolution;
	2151	delete [] newSolution;
	2152	delete [] closestSolution;
	2153	delete [] integerVariable;
	2154	delete [] firstPerturbedObjective;
	2155	delete [] firstPerturbedSolution;
	2156	if (solutionValue == incomingObjective)
	2157	sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
	2158	model_->getCurrentSeconds(), CoinCpuTime() - time1);
[1802]	2159	else if (numberSolutions < maxSolutions)
	2160	sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
	2161	model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
	2162	else
	2163	sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
	2164	model_->getCurrentSeconds(), solutionValue,
	2165	numberSolutions,CoinCpuTime() - time1);
[1286]	2166	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
[1802]	2167	<< pumpPrint
	2168	<< CoinMessageEol;
[1286]	2169	if (bestBasis.getNumStructural())
	2170	model_->setBestSolutionBasis(bestBasis);
	2171	//model_->setMinimizationObjValue(saveBestObjective);
	2172	if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
	2173	model_->setSolutionCount(1); // for local search
	2174	model_->setNumberHeuristicSolutions(1);
[640]	2175	}
[1286]	2176	#ifdef COIN_DEVELOP
	2177	{
	2178	double ncol = model_->solver()->getNumCols();
	2179	double nrow = model_->solver()->getNumRows();
	2180	printf("XXX total iterations %d ratios - %g %g %g\n",
	2181	totalNumberIterations,
	2182	static_cast<double> (totalNumberIterations) / nrow,
	2183	static_cast<double> (totalNumberIterations) / ncol,
	2184	static_cast<double> (totalNumberIterations) / (2nrow + 2ncol));
[640]	2185	}
[931]	2186	#endif
[1286]	2187	return finalReturnCode;
[175]	2188	}
	2189
	2190	/************************END MAIN PROCEDURE *********************************/
	2191
	2192	// update model
	2193	void CbcHeuristicFPump::setModel(CbcModel * model)
	2194	{
[1286]	2195	model_ = model;
[175]	2196	}
	2197
	2198	/* Rounds solution - down if < downValue
	2199	returns 1 if current is a feasible solution
	2200	*/
[1286]	2201	int
	2202	CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
	2203	//const double * objective,
	2204	int numberIntegers, const int * integerVariable,
	2205	/char pumpPrint,*/ int iter,
	2206	/bool roundExpensive,/ double downValue, int *flip)
[175]	2207	{
[1286]	2208	double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
	2209	double primalTolerance ;
	2210	solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
[175]	2211
[1286]	2212	int i;
[175]	2213
[1286]	2214	const double * cost = solver->getObjCoefficients();
	2215	int flip_up = 0;
	2216	int flip_down = 0;
	2217	double v = randomNumberGenerator_.randomDouble() * 20.0;
	2218	int nn = 10 + static_cast<int> (v);
	2219	int nnv = 0;
	2220	int * list = new int [nn];
	2221	double * val = new double [nn];
	2222	for (i = 0; i < nn; i++) val[i] = .001;
[175]	2223
[1286]	2224	const double * rowLower = solver->getRowLower();
	2225	const double * rowUpper = solver->getRowUpper();
	2226	int numberRows = solver->getNumRows();
	2227	if (false && (iter&1) != 0) {
	2228	// Do set covering variables
	2229	const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
	2230	const double * elementByRow = matrixByRow->getElements();
	2231	const int * column = matrixByRow->getIndices();
	2232	const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
	2233	const int * rowLength = matrixByRow->getVectorLengths();
	2234	for (i = 0; i < numberRows; i++) {
	2235	if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
	2236	bool cover = true;
	2237	double largest = 0.0;
	2238	int jColumn = -1;
	2239	for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
	2240	int iColumn = column[k];
	2241	if (elementByRow[k] != 1.0 \|\| !solver->isInteger(iColumn)) {
	2242	cover = false;
	2243	break;
	2244	} else {
	2245	if (solution[iColumn]) {
	2246	double value = solution[iColumn] *
	2247	(randomNumberGenerator_.randomDouble() + 5.0);
	2248	if (value > largest) {
	2249	largest = value;
	2250	jColumn = iColumn;
	2251	}
	2252	}
	2253	}
	2254	}
	2255	if (cover) {
	2256	for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
	2257	int iColumn = column[k];
	2258	if (iColumn == jColumn)
	2259	solution[iColumn] = 1.0;
	2260	else
	2261	solution[iColumn] = 0.0;
	2262	}
	2263	}
	2264	}
	2265	}
	2266	}
	2267	int numberColumns = solver->getNumCols();
[1393]	2268	#ifdef JJF_ZERO
[1053]	2269	// Do set covering variables
	2270	const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
	2271	const double * elementByRow = matrixByRow->getElements();
	2272	const int * column = matrixByRow->getIndices();
	2273	const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
	2274	const int * rowLength = matrixByRow->getVectorLengths();
[1286]	2275	double * sortTemp = new double[numberColumns];
	2276	int * whichTemp = new int [numberColumns];
	2277	char * rowTemp = new char [numberRows];
	2278	memset(rowTemp, 0, numberRows);
	2279	for (i = 0; i < numberColumns; i++)
	2280	whichTemp[i] = -1;
	2281	int nSOS = 0;
	2282	for (i = 0; i < numberRows; i++) {
	2283	if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
	2284	bool cover = true;
	2285	for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
	2286	int iColumn = column[k];
	2287	if (elementByRow[k] != 1.0 \|\| !solver->isInteger(iColumn)) {
	2288	cover = false;
	2289	break;
	2290	}
	2291	}
	2292	if (cover) {
	2293	rowTemp[i] = 1;
	2294	nSOS++;
	2295	for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
	2296	int iColumn = column[k];
	2297	double value = solution[iColumn];
	2298	whichTemp[iColumn] = iColumn;
	2299	}
	2300	}
	2301	}
[1013]	2302	}
[1286]	2303	if (nSOS) {
	2304	// Column copy
	2305	const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
	2306	//const double * element = matrixByColumn->getElements();
	2307	const int * row = matrixByColumn->getIndices();
	2308	const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
	2309	const int * columnLength = matrixByColumn->getVectorLengths();
	2310	int nLook = 0;
	2311	for (i = 0; i < numberColumns; i++) {
	2312	if (whichTemp[i] >= 0) {
	2313	whichTemp[nLook] = i;
	2314	double value = solution[i];
	2315	if (value < 0.5)
	2316	value = (0.1 randomNumberGenerator_.randomDouble() + 0.3);
	2317	sortTemp[nLook++] = -value;
	2318	}
	2319	}
	2320	CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
	2321	double smallest = 1.0;
	2322	int nFix = 0;
	2323	int nOne = 0;
	2324	for (int j = 0; j < nLook; j++) {
	2325	int jColumn = whichTemp[j];
	2326	double thisValue = solution[jColumn];
	2327	if (!thisValue)
	2328	continue;
	2329	if (thisValue == 1.0)
	2330	nOne++;
	2331	smallest = CoinMin(smallest, thisValue);
	2332	solution[jColumn] = 1.0;
	2333	double largest = 0.0;
	2334	for (CoinBigIndex jEl = columnStart[jColumn];
	2335	jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
	2336	int jRow = row[jEl];
	2337	if (rowTemp[jRow]) {
	2338	for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
	2339	int iColumn = column[k];
	2340	if (solution[iColumn]) {
	2341	if (iColumn != jColumn) {
	2342	double value = solution[iColumn];
	2343	if (value > largest)
	2344	largest = value;
	2345	solution[iColumn] = 0.0;
	2346	}
	2347	}
	2348	}
	2349	}
	2350	}
	2351	if (largest > thisValue)
	2352	printf("%d was at %g - chosen over a value of %g\n",
	2353	jColumn, thisValue, largest);
	2354	nFix++;
	2355	}
	2356	printf("%d fixed out of %d (%d at one already)\n",
	2357	nFix, nLook, nOne);
[1013]	2358	}
[1286]	2359	delete [] sortTemp;
	2360	delete [] whichTemp;
	2361	delete [] rowTemp;
[1013]	2362	#endif
[1286]	2363	const double * columnLower = solver->getColLower();
	2364	const double * columnUpper = solver->getColUpper();
	2365	// Check if valid with current solution (allow for 0.99999999s)
	2366	for (i = 0; i < numberIntegers; i++) {
	2367	int iColumn = integerVariable[i];
	2368	double value = solution[iColumn];
	2369	double round = floor(value + 0.5);
	2370	if (fabs(value - round) > primalTolerance)
	2371	break;
[1053]	2372	}
[1286]	2373	if (i == numberIntegers) {
	2374	// may be able to use solution even if 0.99999's
	2375	double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
	2376	double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
	2377	double * saveSolution = CoinCopyOfArray(solution, numberColumns);
	2378	double * tempSolution = CoinCopyOfArray(solution, numberColumns);
	2379	CoinWarmStartBasis * saveBasis =
	2380	dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
	2381	for (i = 0; i < numberIntegers; i++) {
	2382	int iColumn = integerVariable[i];
	2383	double value = solution[iColumn];
	2384	double round = floor(value + 0.5);
	2385	solver->setColLower(iColumn, round);
	2386	solver->setColUpper(iColumn, round);
	2387	tempSolution[iColumn] = round;
	2388	}
	2389	solver->setColSolution(tempSolution);
	2390	delete [] tempSolution;
	2391	solver->resolve();
	2392	solver->setColLower(saveLower);
	2393	solver->setColUpper(saveUpper);
	2394	solver->setWarmStart(saveBasis);
	2395	delete [] saveLower;
	2396	delete [] saveUpper;
	2397	delete saveBasis;
	2398	if (!solver->isProvenOptimal()) {
	2399	solver->setColSolution(saveSolution);
	2400	}
	2401	delete [] saveSolution;
	2402	if (solver->isProvenOptimal()) {
	2403	// feasible
	2404	delete [] list;
	2405	delete [] val;
	2406	return 1;
	2407	}
[1127]	2408	}
[1286]	2409	//double * saveSolution = CoinCopyOfArray(solution,numberColumns);
	2410	// return rounded solution
	2411	for (i = 0; i < numberIntegers; i++) {
	2412	int iColumn = integerVariable[i];
	2413	double value = solution[iColumn];
	2414	double round = floor(value + primalTolerance);
	2415	if (value - round > downValue) round += 1.;
[1393]	2416	#ifndef JJF_ONE
[1286]	2417	if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
	2418	if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
[1013]	2419	#else
[1286]	2420	if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
	2421	if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
[1013]	2422	#endif
[1286]	2423	if (flip_up + flip_down == 0) {
	2424	for (int k = 0; k < nn; k++) {
	2425	if (fabs(value - round) > val[k]) {
	2426	nnv++;
	2427	for (int j = nn - 2; j >= k; j--) {
	2428	val[j+1] = val[j];
	2429	list[j+1] = list[j];
	2430	}
	2431	val[k] = fabs(value - round);
	2432	list[k] = iColumn;
	2433	break;
	2434	}
	2435	}
	2436	}
	2437	solution[iColumn] = round;
[175]	2438	}
	2439
[1286]	2440	if (nnv > nn) nnv = nn;
	2441	//if (iter != 0)
	2442	//sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
	2443	*flip = flip_up + flip_down;
[175]	2444
[1286]	2445	if (*flip == 0 && iter != 0) {
	2446	//sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
	2447	for (i = 0; i < nnv; i++) {
	2448	// was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
	2449	int index = list[i];
	2450	double value = solution[index];
	2451	if (value <= 1.0) {
	2452	solution[index] = 1.0 - value;
	2453	} else if (value < columnLower[index] + integerTolerance) {
	2454	solution[index] = value + 1.0;
	2455	} else if (value > columnUpper[index] - integerTolerance) {
	2456	solution[index] = value - 1.0;
	2457	} else {
	2458	solution[index] = value - 1.0;
	2459	}
	2460	}
	2461	*flip = nnv;
	2462	} else {
	2463	//sprintf(pumpPrint+strlen(pumpPrint)," ");
[1100]	2464	}
[1286]	2465	delete [] list;
	2466	delete [] val;
	2467	//iter++;
	2468
	2469	// get row activities
	2470	double * rowActivity = new double[numberRows];
	2471	memset(rowActivity, 0, numberRows*sizeof(double));
	2472	solver->getMatrixByCol()->times(solution, rowActivity) ;
	2473	double largestInfeasibility = primalTolerance;
	2474	double sumInfeasibility = 0.0;
	2475	int numberBadRows = 0;
	2476	for (i = 0 ; i < numberRows ; i++) {
	2477	double value;
	2478	value = rowLower[i] - rowActivity[i];
	2479	if (value > primalTolerance) {
	2480	numberBadRows++;
	2481	largestInfeasibility = CoinMax(largestInfeasibility, value);
	2482	sumInfeasibility += value;
	2483	}
	2484	value = rowActivity[i] - rowUpper[i];
	2485	if (value > primalTolerance) {
	2486	numberBadRows++;
	2487	largestInfeasibility = CoinMax(largestInfeasibility, value);
	2488	sumInfeasibility += value;
	2489	}
[1100]	2490	}
[1393]	2491	#ifdef JJF_ZERO
[1286]	2492	if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
	2493	// Can we improve by flipping
	2494	for (int iPass = 0; iPass < 10; iPass++) {
	2495	int numberColumns = solver->getNumCols();
	2496	const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
	2497	const double * element = matrixByCol->getElements();
	2498	const int * row = matrixByCol->getIndices();
	2499	const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
	2500	const int * columnLength = matrixByCol->getVectorLengths();
	2501	double oldSum = sumInfeasibility;
	2502	// First improve by moving continuous ones
	2503	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
	2504	if (!solver->isInteger(iColumn)) {
	2505	double solValue = solution[iColumn];
	2506	double thetaUp = columnUpper[iColumn] - solValue;
	2507	double improvementUp = 0.0;
	2508	if (thetaUp > primalTolerance) {
	2509	// can go up
	2510	for (CoinBigIndex j = columnStart[iColumn];
	2511	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2512	int iRow = row[j];
	2513	double distanceUp = rowUpper[iRow] - rowActivity[iRow];
	2514	double distanceDown = rowLower[iRow] - rowActivity[iRow];
	2515	double el = element[j];
	2516	if (el > 0.0) {
	2517	// positive element
	2518	if (distanceUp > 0.0) {
	2519	if (thetaUp*el > distanceUp)
	2520	thetaUp = distanceUp / el;
	2521	} else {
	2522	improvementUp -= el;
	2523	}
	2524	if (distanceDown > 0.0) {
	2525	if (thetaUp*el > distanceDown)
	2526	thetaUp = distanceDown / el;
	2527	improvementUp += el;
	2528	}
	2529	} else {
	2530	// negative element
	2531	if (distanceDown < 0.0) {
	2532	if (thetaUp*el < distanceDown)
	2533	thetaUp = distanceDown / el;
	2534	} else {
	2535	improvementUp += el;
	2536	}
	2537	if (distanceUp < 0.0) {
	2538	if (thetaUp*el < distanceUp)
	2539	thetaUp = distanceUp / el;
	2540	improvementUp -= el;
	2541	}
	2542	}
	2543	}
	2544	}
	2545	double thetaDown = solValue - columnLower[iColumn];
	2546	double improvementDown = 0.0;
	2547	if (thetaDown > primalTolerance) {
	2548	// can go down
	2549	for (CoinBigIndex j = columnStart[iColumn];
	2550	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2551	int iRow = row[j];
	2552	double distanceUp = rowUpper[iRow] - rowActivity[iRow];
	2553	double distanceDown = rowLower[iRow] - rowActivity[iRow];
	2554	double el = -element[j]; // not change in sign form up
	2555	if (el > 0.0) {
	2556	// positive element
	2557	if (distanceUp > 0.0) {
	2558	if (thetaDown*el > distanceUp)
	2559	thetaDown = distanceUp / el;
	2560	} else {
	2561	improvementDown -= el;
	2562	}
	2563	if (distanceDown > 0.0) {
	2564	if (thetaDown*el > distanceDown)
	2565	thetaDown = distanceDown / el;
	2566	improvementDown += el;
	2567	}
	2568	} else {
	2569	// negative element
	2570	if (distanceDown < 0.0) {
	2571	if (thetaDown*el < distanceDown)
	2572	thetaDown = distanceDown / el;
	2573	} else {
	2574	improvementDown += el;
	2575	}
	2576	if (distanceUp < 0.0) {
	2577	if (thetaDown*el < distanceUp)
	2578	thetaDown = distanceUp / el;
	2579	improvementDown -= el;
	2580	}
	2581	}
	2582	}
	2583	if (thetaUp < 1.0e-8)
	2584	improvementUp = 0.0;
	2585	if (thetaDown < 1.0e-8)
	2586	improvementDown = 0.0;
	2587	double theta;
	2588	if (improvementUp >= improvementDown) {
	2589	theta = thetaUp;
	2590	} else {
	2591	improvementUp = improvementDown;
	2592	theta = -thetaDown;
	2593	}
	2594	if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
	2595	// Could move
	2596	double oldSum = 0.0;
	2597	double newSum = 0.0;
	2598	solution[iColumn] += theta;
	2599	for (CoinBigIndex j = columnStart[iColumn];
	2600	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2601	int iRow = row[j];
	2602	double lower = rowLower[iRow];
	2603	double upper = rowUpper[iRow];
	2604	double value = rowActivity[iRow];
	2605	if (value > upper)
	2606	oldSum += value - upper;
	2607	else if (value < lower)
	2608	oldSum += lower - value;
	2609	value += theta * element[j];
	2610	rowActivity[iRow] = value;
	2611	if (value > upper)
	2612	newSum += value - upper;
	2613	else if (value < lower)
	2614	newSum += lower - value;
	2615	}
	2616	assert (newSum <= oldSum);
	2617	sumInfeasibility += newSum - oldSum;
	2618	}
	2619	}
	2620	}
	2621	}
	2622	// Now flip some integers?
[1393]	2623	#ifdef JJF_ZERO
[1286]	2624	for (i = 0; i < numberIntegers; i++) {
	2625	int iColumn = integerVariable[i];
	2626	double solValue = solution[iColumn];
	2627	assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
	2628	double improvementUp = 0.0;
	2629	if (columnUpper[iColumn] >= solValue + 1.0) {
	2630	// can go up
	2631	double oldSum = 0.0;
	2632	double newSum = 0.0;
	2633	for (CoinBigIndex j = columnStart[iColumn];
	2634	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2635	int iRow = row[j];
	2636	double lower = rowLower[iRow];
	2637	double upper = rowUpper[iRow];
	2638	double value = rowActivity[iRow];
	2639	if (value > upper)
	2640	oldSum += value - upper;
	2641	else if (value < lower)
	2642	oldSum += lower - value;
	2643	value += element[j];
	2644	if (value > upper)
	2645	newSum += value - upper;
	2646	else if (value < lower)
	2647	newSum += lower - value;
	2648	}
	2649	improvementUp = oldSum - newSum;
	2650	}
	2651	double improvementDown = 0.0;
	2652	if (columnLower[iColumn] <= solValue - 1.0) {
	2653	// can go down
	2654	double oldSum = 0.0;
	2655	double newSum = 0.0;
	2656	for (CoinBigIndex j = columnStart[iColumn];
	2657	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2658	int iRow = row[j];
	2659	double lower = rowLower[iRow];
	2660	double upper = rowUpper[iRow];
	2661	double value = rowActivity[iRow];
	2662	if (value > upper)
	2663	oldSum += value - upper;
	2664	else if (value < lower)
	2665	oldSum += lower - value;
	2666	value -= element[j];
	2667	if (value > upper)
	2668	newSum += value - upper;
	2669	else if (value < lower)
	2670	newSum += lower - value;
	2671	}
	2672	improvementDown = oldSum - newSum;
	2673	}
	2674	double theta;
	2675	if (improvementUp >= improvementDown) {
	2676	theta = 1.0;
	2677	} else {
	2678	improvementUp = improvementDown;
	2679	theta = -1.0;
	2680	}
	2681	if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
	2682	// Could move
	2683	double oldSum = 0.0;
	2684	double newSum = 0.0;
	2685	solution[iColumn] += theta;
	2686	for (CoinBigIndex j = columnStart[iColumn];
	2687	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2688	int iRow = row[j];
	2689	double lower = rowLower[iRow];
	2690	double upper = rowUpper[iRow];
	2691	double value = rowActivity[iRow];
	2692	if (value > upper)
	2693	oldSum += value - upper;
	2694	else if (value < lower)
	2695	oldSum += lower - value;
	2696	value += theta * element[j];
	2697	rowActivity[iRow] = value;
	2698	if (value > upper)
	2699	newSum += value - upper;
	2700	else if (value < lower)
	2701	newSum += lower - value;
	2702	}
	2703	assert (newSum <= oldSum);
	2704	sumInfeasibility += newSum - oldSum;
	2705	}
	2706	}
[1100]	2707	#else
[1286]	2708	int bestColumn = -1;
	2709	double bestImprovement = primalTolerance;
	2710	double theta = 0.0;
	2711	for (i = 0; i < numberIntegers; i++) {
	2712	int iColumn = integerVariable[i];
	2713	double solValue = solution[iColumn];
	2714	assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
	2715	double improvementUp = 0.0;
	2716	if (columnUpper[iColumn] >= solValue + 1.0) {
	2717	// can go up
	2718	double oldSum = 0.0;
	2719	double newSum = 0.0;
	2720	for (CoinBigIndex j = columnStart[iColumn];
	2721	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2722	int iRow = row[j];
	2723	double lower = rowLower[iRow];
	2724	double upper = rowUpper[iRow];
	2725	double value = rowActivity[iRow];
	2726	if (value > upper)
	2727	oldSum += value - upper;
	2728	else if (value < lower)
	2729	oldSum += lower - value;
	2730	value += element[j];
	2731	if (value > upper)
	2732	newSum += value - upper;
	2733	else if (value < lower)
	2734	newSum += lower - value;
	2735	}
	2736	improvementUp = oldSum - newSum;
	2737	}
	2738	double improvementDown = 0.0;
	2739	if (columnLower[iColumn] <= solValue - 1.0) {
	2740	// can go down
	2741	double oldSum = 0.0;
	2742	double newSum = 0.0;
	2743	for (CoinBigIndex j = columnStart[iColumn];
	2744	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2745	int iRow = row[j];
	2746	double lower = rowLower[iRow];
	2747	double upper = rowUpper[iRow];
	2748	double value = rowActivity[iRow];
	2749	if (value > upper)
	2750	oldSum += value - upper;
	2751	else if (value < lower)
	2752	oldSum += lower - value;
	2753	value -= element[j];
	2754	if (value > upper)
	2755	newSum += value - upper;
	2756	else if (value < lower)
	2757	newSum += lower - value;
	2758	}
	2759	improvementDown = oldSum - newSum;
	2760	}
	2761	double improvement = CoinMax(improvementUp, improvementDown);
	2762	if (improvement > bestImprovement) {
	2763	bestImprovement = improvement;
	2764	bestColumn = iColumn;
	2765	if (improvementUp > improvementDown)
	2766	theta = 1.0;
	2767	else
	2768	theta = -1.0;
	2769	}
	2770	}
	2771	if (bestColumn >= 0) {
	2772	// Could move
	2773	int iColumn = bestColumn;
	2774	double oldSum = 0.0;
	2775	double newSum = 0.0;
	2776	solution[iColumn] += theta;
	2777	for (CoinBigIndex j = columnStart[iColumn];
	2778	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
	2779	int iRow = row[j];
	2780	double lower = rowLower[iRow];
	2781	double upper = rowUpper[iRow];
	2782	double value = rowActivity[iRow];
	2783	if (value > upper)
	2784	oldSum += value - upper;
	2785	else if (value < lower)
	2786	oldSum += lower - value;
	2787	value += theta * element[j];
	2788	rowActivity[iRow] = value;
	2789	if (value > upper)
	2790	newSum += value - upper;
	2791	else if (value < lower)
	2792	newSum += lower - value;
	2793	}
	2794	assert (newSum <= oldSum);
	2795	sumInfeasibility += newSum - oldSum;
	2796	}
[1100]	2797	#endif
[1286]	2798	if (oldSum <= sumInfeasibility + primalTolerance)
	2799	break; // no good
	2800	}
[1100]	2801	}
[1286]	2802	//delete [] saveSolution;
[1100]	2803	#endif
[1286]	2804	delete [] rowActivity;
	2805	return (largestInfeasibility > primalTolerance) ? 0 : 1;
[175]	2806	}
	2807	// Set maximum Time (default off) - also sets starttime to current
[1286]	2808	void
[175]	2809	CbcHeuristicFPump::setMaximumTime(double value)
	2810	{
[1286]	2811	startTime_ = CoinCpuTime();
	2812	maximumTime_ = value;
[175]	2813	}
	2814
[1088]	2815	# ifdef COIN_HAS_CLP
[1286]	2816
[1088]	2817	//#############################################################################
	2818	// Constructors / Destructor / Assignment
	2819	//#############################################################################
	2820
	2821	//-------------------------------------------------------------------
[1286]	2822	// Default Constructor
[1088]	2823	//-------------------------------------------------------------------
[1286]	2824	CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
	2825	: OsiClpDisasterHandler(),
	2826	cbcModel_(model)
[1088]	2827	{
[1286]	2828	if (model) {
	2829	osiModel_
	2830	= dynamic_cast<OsiClpSolverInterface *> (model->solver());
	2831	if (osiModel_)
	2832	setSimplex(osiModel_->getModelPtr());
	2833	}
[1088]	2834	}
	2835
	2836	//-------------------------------------------------------------------
[1286]	2837	// Copy constructor
[1088]	2838	//-------------------------------------------------------------------
[1286]	2839	CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
	2840	: OsiClpDisasterHandler(rhs),
	2841	cbcModel_(rhs.cbcModel_)
	2842	{
[1088]	2843	}
	2844
	2845
	2846	//-------------------------------------------------------------------
[1286]	2847	// Destructor
[1088]	2848	//-------------------------------------------------------------------
	2849	CbcDisasterHandler::~CbcDisasterHandler ()
	2850	{
	2851	}
	2852
	2853	//----------------------------------------------------------------
[1286]	2854	// Assignment operator
[1088]	2855	//-------------------------------------------------------------------
	2856	CbcDisasterHandler &
[1286]	2857	CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
[1088]	2858	{
[1286]	2859	if (this != &rhs) {
	2860	OsiClpDisasterHandler::operator=(rhs);
	2861	cbcModel_ = rhs.cbcModel_;
	2862	}
	2863	return *this;
[1088]	2864	}
	2865	//-------------------------------------------------------------------
	2866	// Clone
	2867	//-------------------------------------------------------------------
	2868	ClpDisasterHandler * CbcDisasterHandler::clone() const
	2869	{
[1286]	2870	return new CbcDisasterHandler(*this);
[1088]	2871	}
	2872	// Type of disaster 0 can fix, 1 abort
[1286]	2873	int
[1088]	2874	CbcDisasterHandler::typeOfDisaster()
	2875	{
[1286]	2876	if (!cbcModel_->parentModel() &&
	2877	(cbcModel_->specialOptions()&2048) == 0) {
	2878	return 0;
	2879	} else {
	2880	if (cbcModel_->parentModel())
	2881	cbcModel_->setMaximumNodes(0);
	2882	return 1;
	2883	}
[1088]	2884	}
	2885	/* set model. */
[1286]	2886	void
[1088]	2887	CbcDisasterHandler::setCbcModel(CbcModel * model)
	2888	{
[1286]	2889	cbcModel_ = model;
	2890	if (model) {
	2891	osiModel_
	2892	= dynamic_cast<OsiClpSolverInterface *> (model->solver());
	2893	if (osiModel_)
	2894	setSimplex(osiModel_->getModelPtr());
	2895	else
	2896	setSimplex(NULL);
	2897	}
[1088]	2898	}
	2899	#endif
[1286]	2900

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: stable/2.8/Cbc/src/CbcHeuristicFPump.cpp

Download in other formats: