Context Navigation

source: trunk/Cbc/src/CbcHeuristicFPump.cpp @ 2094

Last change on this file since 2094 was 2094, checked in by forrest, 3 years ago
for memory leaks and heuristics and some experimental stuff
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 140.7 KB

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

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

Download in other formats: