Context Navigation

source: trunk/Cbc/src/CbcHeuristic.cpp @ 2280

Last change on this file since 2280 was 2280, checked in by forrest, 5 years ago
allow heuristics to see if integers are 'optional'
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 119.9 KB

Line
1	/* $Id: CbcHeuristic.cpp 2280 2016-06-14 14:39:54Z forrest $ */
2	// Copyright (C) 2002, 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
11	#include "CbcConfig.h"
12
13	#include <cassert>
14	#include <cstdlib>
15	#include <cmath>
16	#include <cfloat>
17
18	//#define PRINT_DEBUG
19	#ifdef COIN_HAS_CLP
20	#include "OsiClpSolverInterface.hpp"
21	#endif
22	#include "CbcModel.hpp"
23	#include "CbcMessage.hpp"
24	#include "CbcHeuristic.hpp"
25	#include "CbcHeuristicFPump.hpp"
26	#include "CbcHeuristicRINS.hpp"
27	#include "CbcEventHandler.hpp"
28	#include "CbcStrategy.hpp"
29	#include "CglPreProcess.hpp"
30	#include "CglGomory.hpp"
31	#include "CglProbing.hpp"
32	#include "OsiAuxInfo.hpp"
33	#include "OsiRowCutDebugger.hpp"
34	#include "OsiPresolve.hpp"
35	#include "CbcBranchActual.hpp"
36	#include "CbcCutGenerator.hpp"
37	//==============================================================================
38
39	CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
40	{
41	numObjects_ = rhs.numObjects_;
42	brObj_ = new CbcBranchingObject*[numObjects_];
43	for (int i = 0; i < numObjects_; ++i) {
44	brObj_[i] = rhs.brObj_[i]->clone();
45	}
46	}
47
48	void
49	CbcHeuristicNodeList::gutsOfDelete()
50	{
51	for (int i = (static_cast<int>(nodes_.size())) - 1; i >= 0; --i) {
52	delete nodes_[i];
53	}
54	}
55
56	void
57	CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
58	{
59	append(rhs);
60	}
61
62	CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
63	{
64	gutsOfCopy(rhs);
65	}
66
67	CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
68	(const CbcHeuristicNodeList & rhs)
69	{
70	if (this != &rhs) {
71	gutsOfDelete();
72	gutsOfCopy(rhs);
73	}
74	return *this;
75	}
76
77	CbcHeuristicNodeList::~CbcHeuristicNodeList()
78	{
79	gutsOfDelete();
80	}
81
82	void
83	CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
84	{
85	nodes_.push_back(node);
86	node = NULL;
87	}
88
89	void
90	CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
91	{
92	nodes_.reserve(nodes_.size() + nodes.size());
93	for (int i = 0; i < nodes.size(); ++i) {
94	CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
95	append(node);
96	}
97	}
98
99	//==============================================================================
100	#define DEFAULT_WHERE ((255-2-16)*(1+256))
101	// Default Constructor
102	CbcHeuristic::CbcHeuristic() :
103	model_(NULL),
104	when_(2),
105	numberNodes_(200),
106	feasibilityPumpOptions_(-1),
107	fractionSmall_(1.0),
108	heuristicName_("Unknown"),
109	howOften_(1),
110	decayFactor_(0.0),
111	switches_(0),
112	whereFrom_(DEFAULT_WHERE),
113	shallowDepth_(1),
114	howOftenShallow_(1),
115	numInvocationsInShallow_(0),
116	numInvocationsInDeep_(0),
117	lastRunDeep_(0),
118	numRuns_(0),
119	minDistanceToRun_(1),
120	runNodes_(),
121	numCouldRun_(0),
122	numberSolutionsFound_(0),
123	numberNodesDone_(0),
124	inputSolution_(NULL)
125	{
126	// As CbcHeuristic virtual need to modify .cpp if above change
127	}
128
129	// Constructor from model
130	CbcHeuristic::CbcHeuristic(CbcModel & model) :
131	model_(&model),
132	when_(2),
133	numberNodes_(200),
134	feasibilityPumpOptions_(-1),
135	fractionSmall_(1.0),
136	heuristicName_("Unknown"),
137	howOften_(1),
138	decayFactor_(0.0),
139	switches_(0),
140	whereFrom_(DEFAULT_WHERE),
141	shallowDepth_(1),
142	howOftenShallow_(1),
143	numInvocationsInShallow_(0),
144	numInvocationsInDeep_(0),
145	lastRunDeep_(0),
146	numRuns_(0),
147	minDistanceToRun_(1),
148	runNodes_(),
149	numCouldRun_(0),
150	numberSolutionsFound_(0),
151	numberNodesDone_(0),
152	inputSolution_(NULL)
153	{
154	}
155
156	void
157	CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
158	{
159	model_ = rhs.model_;
160	when_ = rhs.when_;
161	numberNodes_ = rhs.numberNodes_;
162	feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
163	fractionSmall_ = rhs.fractionSmall_;
164	randomNumberGenerator_ = rhs.randomNumberGenerator_;
165	heuristicName_ = rhs.heuristicName_;
166	howOften_ = rhs.howOften_;
167	decayFactor_ = rhs.decayFactor_;
168	switches_ = rhs.switches_;
169	whereFrom_ = rhs.whereFrom_;
170	shallowDepth_ = rhs.shallowDepth_;
171	howOftenShallow_ = rhs.howOftenShallow_;
172	numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
173	numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
174	lastRunDeep_ = rhs.lastRunDeep_;
175	numRuns_ = rhs.numRuns_;
176	numCouldRun_ = rhs.numCouldRun_;
177	minDistanceToRun_ = rhs.minDistanceToRun_;
178	runNodes_ = rhs.runNodes_;
179	numberSolutionsFound_ = rhs.numberSolutionsFound_;
180	numberNodesDone_ = rhs.numberNodesDone_;
181	if (rhs.inputSolution_) {
182	int numberColumns = model_->getNumCols();
183	setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
184	}
185	}
186	// Copy constructor
187	CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
188	{
189	inputSolution_ = NULL;
190	gutsOfCopy(rhs);
191	}
192
193	// Assignment operator
194	CbcHeuristic &
195	CbcHeuristic::operator=( const CbcHeuristic & rhs)
196	{
197	if (this != &rhs) {
198	gutsOfDelete();
199	gutsOfCopy(rhs);
200	}
201	return *this;
202	}
203
204	void CbcHeurDebugNodes(CbcModel* model_)
205	{
206	CbcNode* node = model_->currentNode();
207	CbcNodeInfo* nodeInfo = node->nodeInfo();
208	std::cout << "===============================================================\n";
209	while (nodeInfo) {
210	const CbcNode* node = nodeInfo->owner();
211	printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
212	{
213	const CbcIntegerBranchingObject* brPrint =
214	dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch());
215	if (!brPrint) {
216	printf(" parentBranch: NULL\n");
217	} else {
218	const double* downBounds = brPrint->downBounds();
219	const double* upBounds = brPrint->upBounds();
220	int variable = brPrint->variable();
221	int way = brPrint->way();
222	printf(" parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
223	variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
224	static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
225	}
226	}
227	if (! node) {
228	printf(" owner: NULL\n");
229	} else {
230	printf(" owner: node %i depth %i onTree %i active %i",
231	node->nodeNumber(), node->depth(), node->onTree(), node->active());
232	const OsiBranchingObject* osibr =
233	nodeInfo->owner()->branchingObject();
234	const CbcBranchingObject* cbcbr =
235	dynamic_cast<const CbcBranchingObject*>(osibr);
236	const CbcIntegerBranchingObject* brPrint =
237	dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr);
238	if (!brPrint) {
239	printf(" ownerBranch: NULL\n");
240	} else {
241	const double* downBounds = brPrint->downBounds();
242	const double* upBounds = brPrint->upBounds();
243	int variable = brPrint->variable();
244	int way = brPrint->way();
245	printf(" ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
246	variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
247	static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
248	}
249	}
250	nodeInfo = nodeInfo->parent();
251	}
252	}
253
254	void
255	CbcHeuristic::debugNodes()
256	{
257	CbcHeurDebugNodes(model_);
258	}
259
260	void
261	CbcHeuristic::printDistanceToNodes()
262	{
263	const CbcNode* currentNode = model_->currentNode();
264	if (currentNode != NULL) {
265	CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
266	for (int i = runNodes_.size() - 1; i >= 0; --i) {
267	nodeDesc->distance(runNodes_.node(i));
268	}
269	runNodes_.append(nodeDesc);
270	}
271	}
272
273	bool
274	CbcHeuristic::shouldHeurRun(int whereFrom)
275	{
276	assert (whereFrom >= 0 && whereFrom < 16);
277	// take off 8 (code - likes new solution)
278	whereFrom &= 7;
279	if ((whereFrom_&(1 << whereFrom)) == 0)
280	return false;
281	// No longer used for original purpose - so use for ever run at all JJF
282	#ifndef JJF_ONE
283	// Don't run if hot start or no rows!
284	if (model_ && (model_->hotstartSolution()\|\|!model_->getNumRows()))
285	return false;
286	else
287	return true;
288	#else
289	#ifdef JJF_ZERO
290	const CbcNode* currentNode = model_->currentNode();
291	if (currentNode == NULL) {
292	return false;
293	}
294
295	debugNodes();
296	// return false;
297
298	const int depth = currentNode->depth();
299	#else
300	int depth = model_->currentDepth();
301	#endif
302
303	const int nodeCount = model_->getNodeCount(); // FIXME: check that this is
304	// correct in parallel
305
306	if (nodeCount == 0 \|\| depth <= shallowDepth_) {
307	// what to do when we are in the shallow part of the tree
308	if (model_->getCurrentPassNumber() == 1) {
309	// first time in the node...
310	numInvocationsInShallow_ = 0;
311	}
312	++numInvocationsInShallow_;
313	// Very large howOftenShallow_ will give the original test:
314	// (model_->getCurrentPassNumber() != 1)
315	// if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
316	if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
317	return false;
318	}
319	// LL: should we save these nodes in the list of nodes where the heur was
320	// LL: run?
321	#ifndef JJF_ONE
322	if (currentNode != NULL) {
323	// Get where we are and create the appropriate CbcHeuristicNode object
324	CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
325	runNodes_.append(nodeDesc);
326	}
327	#endif
328	} else {
329	// deeper in the tree
330	if (model_->getCurrentPassNumber() == 1) {
331	// first time in the node...
332	++numInvocationsInDeep_;
333	}
334	if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
335	return false;
336	}
337	if (model_->getCurrentPassNumber() > 1) {
338	// Run the heuristic only when first entering the node.
339	// LL: I don't think this is right. It should run just before strong
340	// LL: branching, I believe.
341	return false;
342	}
343	// Get where we are and create the appropriate CbcHeuristicNode object
344	CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
345	//#ifdef PRINT_DEBUG
346	#ifndef JJF_ONE
347	const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
348	#else
349	const double minDistanceToRun = minDistanceToRun_;
350	#endif
351	#ifdef PRINT_DEBUG
352	double minDistance = nodeDesc->minDistance(runNodes_);
353	std::cout << "minDistance = " << minDistance
354	<< ", minDistanceToRun = " << minDistanceToRun << std::endl;
355	#endif
356	if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
357	delete nodeDesc;
358	return false;
359	}
360	runNodes_.append(nodeDesc);
361	lastRunDeep_ = numInvocationsInDeep_;
362	// ++lastRunDeep_;
363	}
364	++numRuns_;
365	return true;
366	#endif
367	}
368
369	bool
370	CbcHeuristic::shouldHeurRun_randomChoice()
371	{
372	if (!when_)
373	return false;
374	int depth = model_->currentDepth();
375	// when_ -999 is special marker to force to run
376	if (depth != 0 && when_ != -999) {
377	const double numerator = depth * depth;
378	const double denominator = exp(depth * log(2.0));
379	double probability = numerator / denominator;
380	double randomNumber = randomNumberGenerator_.randomDouble();
381	int when = when_ % 100;
382	if (when > 2 && when < 8) {
383	/* JJF adjustments
384	3 only at root and if no solution
385	4 only at root and if this heuristic has not got solution
386	5 decay (but only if no solution)
387	6 if depth <3 or decay
388	7 run up to 2 times if solution found 4 otherwise
389	*/
390	switch (when) {
391	case 3:
392	default:
393	if (model_->bestSolution())
394	probability = -1.0;
395	break;
396	case 4:
397	if (numberSolutionsFound_)
398	probability = -1.0;
399	break;
400	case 5:
401	assert (decayFactor_);
402	if (model_->bestSolution()) {
403	probability = -1.0;
404	} else if (numCouldRun_ > 1000) {
405	decayFactor_ *= 0.99;
406	probability *= decayFactor_;
407	}
408	break;
409	case 6:
410	if (depth >= 3) {
411	if ((numCouldRun_ % howOften_) == 0 &&
412	numberSolutionsFound_*howOften_ < numCouldRun_) {
413	//#define COIN_DEVELOP
414	#ifdef COIN_DEVELOP
415	int old = howOften_;
416	#endif
417	howOften_ = CoinMin(CoinMax(static_cast<int> (howOften_ * 1.1), howOften_ + 1), 1000000);
418	#ifdef COIN_DEVELOP
419	printf("Howoften changed from %d to %d for %s\n",
420	old, howOften_, heuristicName_.c_str());
421	#endif
422	}
423	probability = 1.0 / howOften_;
424	if (model_->bestSolution())
425	probability *= 0.5;
426	} else {
427	probability = 1.1;
428	}
429	break;
430	case 7:
431	if ((model_->bestSolution() && numRuns_ >= 2) \|\| numRuns_ >= 4)
432	probability = -1.0;
433	break;
434	}
435	}
436	if (randomNumber > probability)
437	return false;
438
439	if (model_->getCurrentPassNumber() > 1)
440	return false;
441	#ifdef COIN_DEVELOP
442	printf("Running %s, random %g probability %g\n",
443	heuristicName_.c_str(), randomNumber, probability);
444	#endif
445	} else {
446	#ifdef COIN_DEVELOP
447	printf("Running %s, depth %d when %d\n",
448	heuristicName_.c_str(), depth, when_);
449	#endif
450	}
451	++numRuns_;
452	return true;
453	}
454
455	// Resets stuff if model changes
456	void
457	CbcHeuristic::resetModel(CbcModel * model)
458	{
459	model_ = model;
460	}
461	// Set seed
462	void
463	CbcHeuristic::setSeed(int value)
464	{
465	if (value==0) {
466	double time = fabs(CoinGetTimeOfDay());
467	while (time>=COIN_INT_MAX)
468	time *= 0.5;
469	value = static_cast<int>(time);
470	char printArray[100];
471	sprintf(printArray, "using time of day seed was changed from %d to %d",
472	randomNumberGenerator_.getSeed(), value);
473	if (model_)
474	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
475	<< printArray
476	<< CoinMessageEol;
477	}
478	randomNumberGenerator_.setSeed(value);
479	}
480	// Get seed
481	int
482	CbcHeuristic::getSeed() const
483	{
484	return randomNumberGenerator_.getSeed();
485	}
486
487	// Create C++ lines to get to current state
488	void
489	CbcHeuristic::generateCpp( FILE * fp, const char * heuristic)
490	{
491	// hard coded as CbcHeuristic virtual
492	if (when_ != 2)
493	fprintf(fp, "3 %s.setWhen(%d);\n", heuristic, when_);
494	else
495	fprintf(fp, "4 %s.setWhen(%d);\n", heuristic, when_);
496	if (numberNodes_ != 200)
497	fprintf(fp, "3 %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
498	else
499	fprintf(fp, "4 %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
500	if (feasibilityPumpOptions_ != -1)
501	fprintf(fp, "3 %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
502	else
503	fprintf(fp, "4 %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
504	if (fractionSmall_ != 1.0)
505	fprintf(fp, "3 %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
506	else
507	fprintf(fp, "4 %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
508	if (heuristicName_ != "Unknown")
509	fprintf(fp, "3 %s.setHeuristicName(\"%s\");\n",
510	heuristic, heuristicName_.c_str()) ;
511	else
512	fprintf(fp, "4 %s.setHeuristicName(\"%s\");\n",
513	heuristic, heuristicName_.c_str()) ;
514	if (decayFactor_ != 0.0)
515	fprintf(fp, "3 %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
516	else
517	fprintf(fp, "4 %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
518	if (switches_ != 0)
519	fprintf(fp, "3 %s.setSwitches(%d);\n", heuristic, switches_);
520	else
521	fprintf(fp, "4 %s.setSwitches(%d);\n", heuristic, switches_);
522	if (whereFrom_ != DEFAULT_WHERE)
523	fprintf(fp, "3 %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
524	else
525	fprintf(fp, "4 %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
526	if (shallowDepth_ != 1)
527	fprintf(fp, "3 %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
528	else
529	fprintf(fp, "4 %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
530	if (howOftenShallow_ != 1)
531	fprintf(fp, "3 %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
532	else
533	fprintf(fp, "4 %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
534	if (minDistanceToRun_ != 1)
535	fprintf(fp, "3 %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
536	else
537	fprintf(fp, "4 %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
538	}
539	// Destructor
540	CbcHeuristic::~CbcHeuristic ()
541	{
542	delete [] inputSolution_;
543	}
544
545	// update model
546	void CbcHeuristic::setModel(CbcModel * model)
547	{
548	model_ = model;
549	}
550	/* Clone but ..
551	type 0 clone solver, 1 clone continuous solver
552	Add 2 to say without integer variables which are at low priority
553	Add 4 to say quite likely infeasible so give up easily.*/
554	OsiSolverInterface *
555	CbcHeuristic::cloneBut(int type)
556	{
557	OsiSolverInterface * solver;
558	if ((type&1) == 0 \|\| !model_->continuousSolver())
559	solver = model_->solver()->clone();
560	else
561	solver = model_->continuousSolver()->clone();
562	#ifdef COIN_HAS_CLP
563	OsiClpSolverInterface * clpSolver
564	= dynamic_cast<OsiClpSolverInterface *> (solver);
565	#endif
566	if ((type&2) != 0) {
567	int n = model_->numberObjects();
568	int priority = model_->continuousPriority();
569	if (priority < COIN_INT_MAX) {
570	for (int i = 0; i < n; i++) {
571	const OsiObject * obj = model_->object(i);
572	const CbcSimpleInteger * thisOne =
573	dynamic_cast <const CbcSimpleInteger *> (obj);
574	if (thisOne) {
575	int iColumn = thisOne->columnNumber();
576	if (thisOne->priority() >= priority)
577	solver->setContinuous(iColumn);
578	}
579	}
580	}
581	#ifdef COIN_HAS_CLP
582	if (clpSolver) {
583	for (int i = 0; i < n; i++) {
584	const OsiObject * obj = model_->object(i);
585	const CbcSimpleInteger * thisOne =
586	dynamic_cast <const CbcSimpleInteger *> (obj);
587	if (thisOne) {
588	int iColumn = thisOne->columnNumber();
589	if (clpSolver->isOptionalInteger(iColumn))
590	clpSolver->setContinuous(iColumn);
591	}
592	}
593	}
594	#endif
595	}
596	#ifdef COIN_HAS_CLP
597	if ((type&4) != 0 && clpSolver) {
598	int options = clpSolver->getModelPtr()->moreSpecialOptions();
599	clpSolver->getModelPtr()->setMoreSpecialOptions(options \| 64);
600	}
601	#endif
602	return solver;
603	}
604	// Whether to exit at once on gap
605	bool
606	CbcHeuristic::exitNow(double bestObjective) const
607	{
608	if ((switches_&2048) != 0) {
609	// exit may be forced - but unset for next time
610	switches_ &= ~2048;
611	if ((switches_&1024) != 0)
612	return true;
613	} else if ((switches_&1) == 0) {
614	return false;
615	}
616	// See if can stop on gap
617	OsiSolverInterface * solver = model_->solver();
618	double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
619	double absGap = CoinMax(model_->getAllowableGap(),
620	model_->getHeuristicGap());
621	double fracGap = CoinMax(model_->getAllowableFractionGap(),
622	model_->getHeuristicFractionGap());
623	double testGap = CoinMax(absGap, fracGap *
624	CoinMax(fabs(bestObjective),
625	fabs(bestPossibleObjective)));
626
627	if (bestObjective - bestPossibleObjective < testGap
628	&& model_->getCutoffIncrement() >= 0.0) {
629	return true;
630	} else {
631	return false;
632	}
633	}
634	#ifdef HISTORY_STATISTICS
635	extern bool getHistoryStatistics_;
636	#endif
637	static double sizeRatio(int numberRowsNow, int numberColumnsNow,
638	int numberRowsStart, int numberColumnsStart)
639	{
640	double valueNow;
641	if (numberRowsNow*10 > numberColumnsNow \|\| numberColumnsNow < 200) {
642	valueNow = 2 * numberRowsNow + numberColumnsNow;
643	} else {
644	// long and thin - rows are more important
645	if (numberRowsNow*40 > numberColumnsNow)
646	valueNow = 10 * numberRowsNow + numberColumnsNow;
647	else
648	valueNow = 200 * numberRowsNow + numberColumnsNow;
649	}
650	double valueStart;
651	if (numberRowsStart*10 > numberColumnsStart \|\| numberColumnsStart < 200) {
652	valueStart = 2 * numberRowsStart + numberColumnsStart;
653	} else {
654	// long and thin - rows are more important
655	if (numberRowsStart*40 > numberColumnsStart)
656	valueStart = 10 * numberRowsStart + numberColumnsStart;
657	else
658	valueStart = 200 * numberRowsStart + numberColumnsStart;
659	}
660	//printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
661	// valueNow,numberRowsNow,numberColumnsNow,
662	// valueStart,numberRowsStart,numberColumnsStart);
663	if (10numberRowsNow < 8numberRowsStart \|\| 10numberColumnsNow < 7numberColumnsStart)
664	return valueNow / valueStart;
665	else if (10numberRowsNow < 9numberRowsStart)
666	return 1.1*(valueNow / valueStart);
667	else if (numberRowsNow < numberRowsStart)
668	return 1.5*(valueNow / valueStart);
669	else
670	return 2.0*(valueNow / valueStart);
671	}
672
673	//static int saveModel=0;
674	// Do mini branch and bound (return 1 if solution)
675	int
676	CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver, int numberNodes,
677	double * newSolution, double & newSolutionValue,
678	double cutoff, std::string name) const
679	{
680	CbcEventHandler *eventHandler = model_->getEventHandler() ;
681	// Use this fraction
682	double fractionSmall = fractionSmall_;
683	int maximumSolutions = model_->getMaximumSolutions();
684	int iterationMultiplier = 100;
685	if (eventHandler) {
686	typedef struct {
687	double fractionSmall;
688	double spareDouble[3];
689	OsiSolverInterface * solver;
690	void * sparePointer[2];
691	int numberNodes;
692	int maximumSolutions;
693	int iterationMultiplier;
694	int howOften;
695	int spareInt[3];
696	} SmallMod;
697	SmallMod temp;
698	temp.solver=solver;
699	temp.fractionSmall=fractionSmall;
700	temp.numberNodes=numberNodes;
701	temp.iterationMultiplier=iterationMultiplier;
702	temp.howOften=howOften_;
703	temp.maximumSolutions=maximumSolutions;
704	CbcEventHandler::CbcAction status =
705	eventHandler->event(CbcEventHandler::smallBranchAndBound,
706	&temp);
707	if (status==CbcEventHandler::killSolution)
708	return -1;
709	if (status==CbcEventHandler::takeAction) {
710	fractionSmall=temp.fractionSmall;
711	numberNodes=temp.numberNodes;
712	iterationMultiplier=temp.iterationMultiplier;
713	howOften_=temp.howOften;
714	maximumSolutions=temp.maximumSolutions;
715	}
716	}
717	#if 0
718	if (saveModel \|\| model_->getMaximumSolutions()==100) {
719	printf("writing model\n");
720	solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
721	}
722	#endif
723	// size before
724	int shiftRows = 0;
725	if (numberNodes < 0)
726	shiftRows = solver->getNumRows() - numberNodes_;
727	int numberRowsStart = solver->getNumRows() - shiftRows;
728	int numberColumnsStart = solver->getNumCols();
729	#ifdef CLP_INVESTIGATE
730	printf("%s has %d rows, %d columns\n",
731	name.c_str(), solver->getNumRows(), solver->getNumCols());
732	#endif
733	double before = 2 * numberRowsStart + numberColumnsStart;
734	if (before > 40000.0) {
735	// fairly large - be more conservative
736	double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
737	if (multiplier < 1.0) {
738	fractionSmall *= multiplier;
739	#ifdef CLP_INVESTIGATE
740	printf("changing fractionSmall from %g to %g for %s\n",
741	fractionSmall_, fractionSmall, name.c_str());
742	#endif
743	}
744	}
745	#ifdef COIN_HAS_CLP
746	OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
747	if (osiclp && (osiclp->specialOptions()&65536) == 0) {
748	// go faster stripes
749	if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500) {
750	osiclp->setupForRepeatedUse(2, 0);
751	} else {
752	osiclp->setupForRepeatedUse(0, 0);
753	}
754	// Turn this off if you get problems
755	// Used to be automatically set
756	osiclp->setSpecialOptions(osiclp->specialOptions() \| (128 + 64 - 128));
757	ClpSimplex * lpSolver = osiclp->getModelPtr();
758	lpSolver->setSpecialOptions(lpSolver->specialOptions() \| 0x01000000); // say is Cbc (and in branch and bound)
759	lpSolver->setSpecialOptions(lpSolver->specialOptions() \|
760	(/16384+/4096 + 512 + 128));
761	}
762	#endif
763	#ifdef HISTORY_STATISTICS
764	getHistoryStatistics_ = false;
765	#endif
766	#ifdef COIN_DEVELOP
767	int status = 0;
768	#endif
769	int logLevel = model_->logLevel();
770	#define LEN_PRINT 250
771	char generalPrint[LEN_PRINT];
772	// Do presolve to see if possible
773	int numberColumns = solver->getNumCols();
774	char * reset = NULL;
775	int returnCode = 1;
776	int saveModelOptions = model_->specialOptions();
777	//assert ((saveModelOptions&2048) == 0);
778	model_->setSpecialOptions(saveModelOptions \| 2048);
779	if (fractionSmall<1.0) {
780	int saveLogLevel = solver->messageHandler()->logLevel();
781	if (saveLogLevel == 1)
782	solver->messageHandler()->setLogLevel(0);
783	OsiPresolve * pinfo = new OsiPresolve();
784	int presolveActions = 0;
785	// Allow dual stuff on integers
786	presolveActions = 1;
787	// Do not allow all +1 to be tampered with
788	//if (allPlusOnes)
789	//presolveActions \|= 2;
790	// allow transfer of costs
791	// presolveActions \|= 4;
792	pinfo->setPresolveActions(presolveActions);
793	OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
794	delete pinfo;
795	// see if too big
796
797	if (presolvedModel) {
798	int afterRows = presolvedModel->getNumRows();
799	int afterCols = presolvedModel->getNumCols();
800	//#define COIN_DEVELOP
801	#ifdef COIN_DEVELOP_z
802	if (numberNodes < 0) {
803	solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
804	presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
805	}
806	#endif
807	delete presolvedModel;
808	double ratio = sizeRatio(afterRows - shiftRows, afterCols,
809	numberRowsStart, numberColumnsStart);
810	double after = 2 * afterRows + afterCols;
811	if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
812	// Need code to try again to compress further using used
813	const int * used = model_->usedInSolution();
814	int maxUsed = 0;
815	int iColumn;
816	const double * lower = solver->getColLower();
817	const double * upper = solver->getColUpper();
818	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
819	if (upper[iColumn] > lower[iColumn]) {
820	if (solver->isBinary(iColumn))
821	maxUsed = CoinMax(maxUsed, used[iColumn]);
822	}
823	}
824	if (maxUsed) {
825	reset = new char [numberColumns];
826	int nFix = 0;
827	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
828	reset[iColumn] = 0;
829	if (upper[iColumn] > lower[iColumn]) {
830	if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
831	bool setValue = true;
832	if (maxUsed == 1) {
833	double randomNumber = randomNumberGenerator_.randomDouble();
834	if (randomNumber > 0.3)
835	setValue = false;
836	}
837	if (setValue) {
838	reset[iColumn] = 1;
839	solver->setColLower(iColumn, 1.0);
840	nFix++;
841	}
842	}
843	}
844	}
845	pinfo = new OsiPresolve();
846	presolveActions = 0;
847	// Allow dual stuff on integers
848	presolveActions = 1;
849	// Do not allow all +1 to be tampered with
850	//if (allPlusOnes)
851	//presolveActions \|= 2;
852	// allow transfer of costs
853	// presolveActions \|= 4;
854	pinfo->setPresolveActions(presolveActions);
855	presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
856	delete pinfo;
857	if (presolvedModel) {
858	// see if too big
859	int afterRows2 = presolvedModel->getNumRows();
860	int afterCols2 = presolvedModel->getNumCols();
861	delete presolvedModel;
862	double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
863	numberRowsStart, numberColumnsStart);
864	double after = 2 * afterRows2 + afterCols2;
865	if (ratio > fractionSmall && (after > 300 \|\| numberNodes < 0)) {
866	sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
867	solver->getNumRows(), solver->getNumCols(),
868	afterRows, afterCols, nFix, afterRows2, afterCols2);
869	// If much too big - give up
870	if (ratio > 0.75)
871	returnCode = -1;
872	} else {
873	sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
874	solver->getNumRows(), solver->getNumCols(),
875	afterRows, afterCols, nFix, afterRows2, afterCols2);
876	}
877	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
878	<< generalPrint
879	<< CoinMessageEol;
880	} else {
881	returnCode = 2; // infeasible
882	}
883	}
884	} else if (ratio > fractionSmall && after > 300 && numberNodes >=0) {
885	returnCode = -1;
886	}
887	} else {
888	returnCode = 2; // infeasible
889	}
890	solver->messageHandler()->setLogLevel(saveLogLevel);
891	}
892	if (returnCode == 2 \|\| returnCode == -1) {
893	model_->setSpecialOptions(saveModelOptions);
894	delete [] reset;
895	#ifdef HISTORY_STATISTICS
896	getHistoryStatistics_ = true;
897	#endif
898	//printf("small no good\n");
899	return returnCode;
900	}
901	// Reduce printout
902	bool takeHint;
903	OsiHintStrength strength;
904	solver->getHintParam(OsiDoReducePrint, takeHint, strength);
905	solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
906	solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
907	double signedCutoff = cutoff*solver->getObjSense();
908	solver->setDblParam(OsiDualObjectiveLimit, signedCutoff);
909	solver->initialSolve();
910	if (solver->isProvenOptimal()) {
911	CglPreProcess process;
912	OsiSolverInterface * solver2 = NULL;
913	if ((model_->moreSpecialOptions()&65536)!=0)
914	process.setOptions(2+4+8+16); // no cuts
915	else
916	process.setOptions(16); // no complicated dupcol stuff
917	/* Do not try and produce equality cliques and
918	do up to 2 passes (normally) 5 if restart */
919	int numberPasses = 2;
920	if ((model_->moreSpecialOptions2()&16)!=0) {
921	// quick
922	process.setOptions(2+4+8+16); // no cuts
923	numberPasses = 1;
924	}
925	if (numberNodes < 0) {
926	numberPasses = 5;
927	// Say some rows cuts
928	int numberRows = solver->getNumRows();
929	if (numberNodes_ < numberRows && true /* think */) {
930	char * type = new char[numberRows];
931	memset(type, 0, numberNodes_);
932	memset(type + numberNodes_, 1, numberRows - numberNodes_);
933	process.passInRowTypes(type, numberRows);
934	delete [] type;
935	}
936	}
937	if (logLevel <= 1)
938	process.messageHandler()->setLogLevel(0);
939	if (!solver->defaultHandler()&&
940	solver->messageHandler()->logLevel(0)!=-1000)
941	process.passInMessageHandler(solver->messageHandler());
942	#ifdef CGL_DEBUG
943	/*
944	We're debugging. (specialOptions 1)
945	*/
946	if ((model_->specialOptions()&1) != 0) {
947	const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
948	if (debugger) {
949	process.setApplicationData(const_cast<double *>(debugger->optimalSolution()));
950	}
951	}
952	#endif
953	solver2 = process.preProcessNonDefault(*solver, false,
954	numberPasses);
955	if (!solver2) {
956	if (logLevel > 1)
957	printf("Pre-processing says infeasible\n");
958	returnCode = 2; // so will be infeasible
959	} else {
960	#ifdef COIN_DEVELOP_z
961	if (numberNodes < 0) {
962	solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
963	}
964	#endif
965	// see if too big
966	double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
967	numberRowsStart, numberColumnsStart);
968	double after = 2 * solver2->getNumRows() + solver2->getNumCols();
969	if (ratio > fractionSmall && (after > 300 \|\| numberNodes < 0)) {
970	sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
971	solver->getNumRows(), solver->getNumCols(),
972	solver2->getNumRows(), solver2->getNumCols());
973	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
974	<< generalPrint
975	<< CoinMessageEol;
976	returnCode = -1;
977	//printf("small no good2\n");
978	} else {
979	sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
980	solver->getNumRows(), solver->getNumCols(),
981	solver2->getNumRows(), solver2->getNumCols());
982	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
983	<< generalPrint
984	<< CoinMessageEol;
985	}
986	#ifdef CGL_DEBUG
987	if ((model_->specialOptions()&1) != 0) {
988	const OsiRowCutDebugger *debugger = solver2->getRowCutDebugger() ;
989	if (debugger) {
990	printf("On optimal path after preprocessing\n");
991	}
992	}
993	#endif
994	if (returnCode == 1) {
995	solver2->resolve();
996	CbcModel model(*solver2);
997	double startTime=model_->getDblParam(CbcModel::CbcStartSeconds);
998	model.setDblParam(CbcModel::CbcStartSeconds,startTime);
999	// move seed across
1000	model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
1001	if (numberNodes >= 0) {
1002	// normal
1003	model.setSpecialOptions(saveModelOptions \| 2048);
1004	if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
1005	model.setLogLevel(0);
1006	else
1007	model.setLogLevel(logLevel);
1008	// No small fathoming
1009	model.setFastNodeDepth(-1);
1010	model.setCutoff(signedCutoff);
1011	model.setStrongStrategy(0);
1012	// Don't do if original fraction > 1.0 and too large
1013	if (fractionSmall_>1.0 && fractionSmall_ < 1000000.0) {
1014	/* 1.4 means -1 nodes if >.4
1015	2.4 means -1 nodes if >.5 and 0 otherwise
1016	3.4 means -1 nodes if >.6 and 0 or 5
1017	4.4 means -1 nodes if >.7 and 0, 5 or 10
1018	*/
1019	double fraction = fractionSmall_-floor(fractionSmall_);
1020	if (ratio>fraction) {
1021	int type = static_cast<int>(floor(fractionSmall_*0.1));
1022	int over = static_cast<int>(ceil(ratio-fraction));
1023	int maxNodes[]={-1,0,5,10};
1024	if (type>over)
1025	numberNodes=maxNodes[type-over];
1026	else
1027	numberNodes=-1;
1028	}
1029	}
1030	model.setMaximumNodes(numberNodes);
1031	model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1032	if ((saveModelOptions&2048) == 0)
1033	model.setMoreSpecialOptions(model_->moreSpecialOptions());
1034	model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
1035	// off conflict analysis
1036	model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304);
1037
1038	// Lightweight
1039	CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
1040	model.setStrategy(strategy);
1041	model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
1042	model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
1043	model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
1044	// Set best solution (even if bad for this submodel)
1045	if (model_->bestSolution()) {
1046	const double * bestSolution = model_->bestSolution();
1047	int numberColumns2 = model.solver()->getNumCols();
1048	double * bestSolution2 = new double [numberColumns2];
1049	const int * originalColumns = process.originalColumns();
1050	for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
1051	int jColumn = originalColumns[iColumn];
1052	bestSolution2[iColumn] = bestSolution[jColumn];
1053	}
1054	model.setBestSolution(bestSolution2,numberColumns2,
1055	1.0e50,
1056	false);
1057	model.setSolutionCount(1);
1058	maximumSolutions++;
1059	delete [] bestSolution2;
1060	}
1061	} else {
1062	model.setSpecialOptions(saveModelOptions);
1063	model_->messageHandler()->message(CBC_RESTART, model_->messages())
1064	<< solver2->getNumRows() << solver2->getNumCols()
1065	<< CoinMessageEol;
1066	// going for full search and copy across more stuff
1067	model.gutsOfCopy(*model_, 2);
1068	#ifdef CGL_DEBUG
1069	if ((model_->specialOptions()&1) != 0) {
1070	const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ;
1071	if (debugger) {
1072	printf("On optimal path BB\n");
1073	}
1074	}
1075	#endif
1076	assert (!model_->heuristicModel());
1077	model_->setHeuristicModel(&model);
1078	for (int i = 0; i < model.numberCutGenerators(); i++) {
1079	CbcCutGenerator * generator = model.cutGenerator(i);
1080	CglGomory * gomory = dynamic_cast<CglGomory *>
1081	(generator->generator());
1082	if (gomory&&gomory->originalSolver())
1083	gomory->passInOriginalSolver(model.solver());
1084	generator->setTiming(true);
1085	// Turn on if was turned on
1086	int iOften = model_->cutGenerator(i)->howOften();
1087	#ifdef CLP_INVESTIGATE
1088	printf("Gen %d often %d %d\n",
1089	i, generator->howOften(),
1090	iOften);
1091	#endif
1092	if (iOften > 0)
1093	generator->setHowOften(iOften % 1000000);
1094	if (model_->cutGenerator(i)->howOftenInSub() == -200)
1095	generator->setHowOften(-100);
1096	}
1097	model.setCutoff(signedCutoff);
1098	// make sure can't do nested search! but allow heuristics
1099	model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) \| 1024);
1100	bool takeHint;
1101	OsiHintStrength strength;
1102	// Switch off printing if asked to
1103	model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
1104	model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
1105	// no cut generators if none in parent
1106	CbcStrategyDefault
1107	strategy(model_->numberCutGenerators() ? 1 : -1,
1108	model_->numberStrong(),
1109	model_->numberBeforeTrust());
1110	// Set up pre-processing - no
1111	strategy.setupPreProcessing(0); // was (4);
1112	model.setStrategy(strategy);
1113	//model.solver()->writeMps("crunched");
1114	int numberCuts = process.cuts().sizeRowCuts();
1115	if (numberCuts) {
1116	// add in cuts
1117	CglStored cuts = process.cuts();
1118	model.addCutGenerator(&cuts, 1, "Stored from first");
1119	model.cutGenerator(model.numberCutGenerators()-1)->setGlobalCuts(true);
1120	}
1121	}
1122	// Do search
1123	if (logLevel > 1)
1124	model_->messageHandler()->message(CBC_START_SUB, model_->messages())
1125	<< name
1126	<< model.getMaximumNodes()
1127	<< CoinMessageEol;
1128	// probably faster to use a basis to get integer solutions
1129	model.setSpecialOptions(model.specialOptions() \| 2);
1130	#ifdef CBC_THREAD
1131	if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) {
1132	// See if at root node
1133	bool atRoot = model_->getNodeCount() == 0;
1134	int passNumber = model_->getCurrentPassNumber();
1135	if (atRoot && passNumber == 1)
1136	model.setNumberThreads(model_->getNumberThreads());
1137	}
1138	#endif
1139	model.setParentModel(*model_);
1140	model.setMaximumSolutions(maximumSolutions);
1141	model.setOriginalColumns(process.originalColumns());
1142	model.setSearchStrategy(-1);
1143	// If no feasibility pump then insert a lightweight one
1144	if (feasibilityPumpOptions_ >= 0 \|\| feasibilityPumpOptions_ == -2) {
1145	CbcHeuristicFPump * fpump = NULL;
1146	for (int i = 0; i < model.numberHeuristics(); i++) {
1147	CbcHeuristicFPump* pump =
1148	dynamic_cast<CbcHeuristicFPump*>(model.heuristic(i));
1149	if (pump) {
1150	fpump = pump;
1151	break;
1152	}
1153	}
1154	if (!fpump) {
1155	CbcHeuristicFPump heuristic4;
1156	// use any cutoff
1157	heuristic4.setFakeCutoff(0.5*COIN_DBL_MAX);
1158	if (fractionSmall_<=1.0)
1159	heuristic4.setMaximumPasses(10);
1160	int pumpTune = feasibilityPumpOptions_;
1161	if (pumpTune==-2)
1162	pumpTune = 4; // proximity
1163	if (pumpTune > 0) {
1164	/*
1165	>=10000000 for using obj
1166	>=1000000 use as accumulate switch
1167	>=1000 use index+1 as number of large loops
1168	>=100 use 0.05 objvalue as increment
1169	%100 == 10,20 etc for experimentation
1170	1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
1171	4 and static continuous, 5 as 3 but no internal integers
1172	6 as 3 but all slack basis!
1173	*/
1174	double value = solver2->getObjSense() * solver2->getObjValue();
1175	int w = pumpTune / 10;
1176	int ix = w % 10;
1177	w /= 10;
1178	int c = w % 10;
1179	w /= 10;
1180	int r = w;
1181	int accumulate = r / 1000;
1182	r -= 1000 * accumulate;
1183	if (accumulate >= 10) {
1184	int which = accumulate / 10;
1185	accumulate -= 10 * which;
1186	which--;
1187	// weights and factors
1188	double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
1189	double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
1190	heuristic4.setInitialWeight(weight[which]);
1191	heuristic4.setWeightFactor(factor[which]);
1192	}
1193	// fake cutoff
1194	if (c) {
1195	double cutoff;
1196	solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
1197	cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
1198	heuristic4.setFakeCutoff(cutoff);
1199	}
1200	if (r) {
1201	// also set increment
1202	//double increment = (0.01i+0.005)(fabs(value)+1.0e-12);
1203	double increment = 0.0;
1204	heuristic4.setAbsoluteIncrement(increment);
1205	heuristic4.setAccumulate(accumulate);
1206	heuristic4.setMaximumRetries(r + 1);
1207	}
1208	pumpTune = pumpTune % 100;
1209	if (pumpTune == 6)
1210	pumpTune = 13;
1211	if (pumpTune != 13)
1212	pumpTune = pumpTune % 10;
1213	heuristic4.setWhen(pumpTune);
1214	if (ix) {
1215	heuristic4.setFeasibilityPumpOptions(ix*10);
1216	}
1217	}
1218	model.addHeuristic(&heuristic4, "feasibility pump", 0);
1219
1220	}
1221	} else if (feasibilityPumpOptions_==-3) {
1222	// add all (except this)
1223	for (int i = 0; i < model_->numberHeuristics(); i++) {
1224	if (strcmp(heuristicName(),model_->heuristic(i)->heuristicName()))
1225	model.addHeuristic(model_->heuristic(i));
1226	}
1227	}
1228	// modify heuristics
1229	for (int i = 0; i < model.numberHeuristics(); i++) {
1230	// reset lastNode
1231	CbcHeuristicRINS * rins =
1232	dynamic_cast<CbcHeuristicRINS*>(model.heuristic(i));
1233	if (rins) {
1234	rins->setLastNode(-1000);
1235	rins->setSolutionCount(0);
1236	}
1237	}
1238	//printf("sol %x\n",inputSolution_);
1239	#ifdef CGL_DEBUG
1240	if ((model_->specialOptions()&1) != 0) {
1241	const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ;
1242	if (debugger) {
1243	printf("On optimal path CC\n");
1244	}
1245	}
1246	#endif
1247	if (inputSolution_) {
1248	// translate and add a serendipity heuristic
1249	int numberColumns = solver2->getNumCols();
1250	const int * which = process.originalColumns();
1251	OsiSolverInterface * solver3 = solver2->clone();
1252	for (int i = 0; i < numberColumns; i++) {
1253	if (isHeuristicInteger(solver3,i)) {
1254	int k = which[i];
1255	double value = inputSolution_[k];
1256	//if (value)
1257	//printf("orig col %d now %d val %g\n",
1258	// k,i,value);
1259	solver3->setColLower(i, value);
1260	solver3->setColUpper(i, value);
1261	}
1262	}
1263	solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1264	solver3->resolve();
1265	if (!solver3->isProvenOptimal()) {
1266	// Try just setting nonzeros
1267	OsiSolverInterface * solver4 = solver2->clone();
1268	for (int i = 0; i < numberColumns; i++) {
1269	if (isHeuristicInteger(solver4,i)) {
1270	int k = which[i];
1271	double value = floor(inputSolution_[k] + 0.5);
1272	if (value) {
1273	solver3->setColLower(i, value);
1274	solver3->setColUpper(i, value);
1275	}
1276	}
1277	}
1278	solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1279	solver4->resolve();
1280	int nBad = -1;
1281	if (solver4->isProvenOptimal()) {
1282	nBad = 0;
1283	const double * solution = solver4->getColSolution();
1284	for (int i = 0; i < numberColumns; i++) {
1285	if (isHeuristicInteger(solver4,i)) {
1286	double value = floor(solution[i] + 0.5);
1287	if (fabs(value - solution[i]) > 1.0e-6)
1288	nBad++;
1289	}
1290	}
1291	}
1292	if (nBad) {
1293	delete solver4;
1294	} else {
1295	delete solver3;
1296	solver3 = solver4;
1297	}
1298	}
1299	if (solver3->isProvenOptimal()) {
1300	// good
1301	CbcSerendipity heuristic(model);
1302	double value = solver3->getObjSense() * solver3->getObjValue();
1303	heuristic.setInputSolution(solver3->getColSolution(), value);
1304	value = value + 1.0e-7*(1.0 + fabs(value));
1305	value *= solver3->getObjSense();
1306	model.setCutoff(value);
1307	model.addHeuristic(&heuristic, "Previous solution", 0);
1308	//printf("added seren\n");
1309	} else {
1310	double value = model_->getMinimizationObjValue();
1311	value = value + 1.0e-7*(1.0 + fabs(value));
1312	value *= solver3->getObjSense();
1313	model.setCutoff(value);
1314	sprintf(generalPrint, "Unable to insert previous solution - using cutoff of %g",
1315	value);
1316	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1317	<< generalPrint
1318	<< CoinMessageEol;
1319	#ifdef CLP_INVESTIGATE
1320	printf("NOT added seren\n");
1321	solver3->writeMps("bad_seren");
1322	solver->writeMps("orig_seren");
1323	#endif
1324	}
1325	delete solver3;
1326	}
1327	if (model_->searchStrategy() == 2) {
1328	model.setNumberStrong(5);
1329	model.setNumberBeforeTrust(5);
1330	}
1331	if (model.getNumCols()) {
1332	if (numberNodes >= 0) {
1333	setCutAndHeuristicOptions(model);
1334	// not too many iterations
1335	model.setMaximumNumberIterations(iterationMultiplier*(numberNodes + 10));
1336	// Not fast stuff
1337	model.setFastNodeDepth(-1);
1338	//model.solver()->writeMps("before");
1339	} else if (model.fastNodeDepth() >= 1000000) {
1340	// already set
1341	model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
1342	}
1343	model.setWhenCuts(999998);
1344	#define ALWAYS_DUAL
1345	#ifdef ALWAYS_DUAL
1346	OsiSolverInterface * solverD = model.solver();
1347	bool takeHint;
1348	OsiHintStrength strength;
1349	solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
1350	solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
1351	#endif
1352	model.passInEventHandler(model_->getEventHandler());
1353	// say model_ is sitting there
1354	int saveOptions = model_->specialOptions();
1355	model_->setSpecialOptions(saveOptions\|1048576);
1356	// and switch off debugger
1357	model.setSpecialOptions(model.specialOptions()&(~1));
1358	#if 0 //def COIN_HAS_CLP
1359	OsiClpSolverInterface * clpSolver
1360	= dynamic_cast<OsiClpSolverInterface *> (model.solver());
1361	if (clpSolver)
1362	clpSolver->zapDebugger();
1363	#endif
1364	#ifdef CONFLICT_CUTS
1365	if ((model_->moreSpecialOptions()&4194304)!=0)
1366	model.zapGlobalCuts();
1367	#endif
1368	#ifdef CGL_DEBUG
1369	if ((model_->specialOptions()&1) != 0) {
1370	const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ;
1371	if (debugger) {
1372	printf("On optimal path DD\n");
1373	}
1374	}
1375	#endif
1376	model.branchAndBound();
1377	model_->setHeuristicModel(NULL);
1378	model_->setSpecialOptions(saveOptions);
1379	#ifdef ALWAYS_DUAL
1380	solverD = model.solver();
1381	solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
1382	#endif
1383	numberNodesDone_ = model.getNodeCount();
1384	#ifdef COIN_DEVELOP
1385	printf("sub branch %d nodes, %d iterations - max %d\n",
1386	model.getNodeCount(), model.getIterationCount(),
1387	100*(numberNodes + 10));
1388	#endif
1389	if (numberNodes < 0) {
1390	model_->incrementIterationCount(model.getIterationCount());
1391	model_->incrementNodeCount(model.getNodeCount());
1392	// update best solution (in case ctrl-c)
1393	// !!! not a good idea - think a bit harder
1394	//model_->setMinimizationObjValue(model.getMinimizationObjValue());
1395	for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
1396	CbcCutGenerator * generator = model.cutGenerator(iGenerator);
1397	sprintf(generalPrint,
1398	"%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
1399	generator->cutGeneratorName(),
1400	generator->numberTimesEntered(),
1401	generator->numberCutsInTotal() +
1402	generator->numberColumnCuts(),
1403	generator->numberCutsActive(),
1404	generator->timeInCutGenerator());
1405	CglStored * stored = dynamic_cast<CglStored*>(generator->generator());
1406	if (stored && !generator->numberCutsInTotal())
1407	continue;
1408	#ifndef CLP_INVESTIGATE
1409	CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
1410	if (implication)
1411	continue;
1412	#endif
1413	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1414	<< generalPrint
1415	<< CoinMessageEol;
1416	}
1417	}
1418	} else {
1419	// empty model
1420	model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
1421	}
1422	if (logLevel > 1)
1423	model_->messageHandler()->message(CBC_END_SUB, model_->messages())
1424	<< name
1425	<< CoinMessageEol;
1426	if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
1427	// solution
1428	if (model.getNumCols())
1429	returnCode = model.isProvenOptimal() ? 3 : 1;
1430	else
1431	returnCode = 3;
1432	// post process
1433	#ifdef COIN_HAS_CLP
1434	OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
1435	if (clpSolver) {
1436	ClpSimplex * lpSolver = clpSolver->getModelPtr();
1437	lpSolver->setSpecialOptions(lpSolver->specialOptions() \| 0x01000000); // say is Cbc (and in branch and bound)
1438	}
1439	#endif
1440	//if (fractionSmall_ < 1000000.0)
1441	process.postProcess(*model.solver());
1442	if (solver->isProvenOptimal() && solver->getObjValue()*solver->getObjSense() < cutoff) {
1443	// Solution now back in solver
1444	int numberColumns = solver->getNumCols();
1445	memcpy(newSolution, solver->getColSolution(),
1446	numberColumns*sizeof(double));
1447	newSolutionValue = model.getMinimizationObjValue();
1448	} else {
1449	// odd - but no good
1450	returnCode = 0; // so will be infeasible
1451	}
1452	} else {
1453	// no good
1454	returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
1455	}
1456	int totalNumberIterations = model.getIterationCount() +
1457	process.numberIterationsPre() +
1458	process.numberIterationsPost();
1459	if (totalNumberIterations > 100*(numberNodes + 10)
1460	&& fractionSmall_ < 1000000.0) {
1461	// only allow smaller problems
1462	fractionSmall = fractionSmall_;
1463	fractionSmall_ *= 0.9;
1464	#ifdef CLP_INVESTIGATE
1465	printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
1466	fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
1467	#endif
1468	}
1469	if (model.status() == 5)
1470	model_->sayEventHappened();
1471	#ifdef COIN_DEVELOP
1472	if (model.isProvenInfeasible())
1473	status = 1;
1474	else if (model.isProvenOptimal())
1475	status = 2;
1476	#endif
1477	}
1478	}
1479	} else {
1480	returnCode = 2; // infeasible finished
1481	if (logLevel > 1){
1482	printf("Infeasible on initial solve\n");
1483	}
1484	}
1485	model_->setSpecialOptions(saveModelOptions);
1486	model_->setLogLevel(logLevel);
1487	if (returnCode == 1 \|\| returnCode == 2) {
1488	OsiSolverInterface * solverC = model_->continuousSolver();
1489	if (false && solverC) {
1490	const double * lower = solver->getColLower();
1491	const double * upper = solver->getColUpper();
1492	const double * lowerC = solverC->getColLower();
1493	const double * upperC = solverC->getColUpper();
1494	bool good = true;
1495	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1496	if (isHeuristicInteger(solverC,iColumn)) {
1497	if (lower[iColumn] > lowerC[iColumn] &&
1498	upper[iColumn] < upperC[iColumn]) {
1499	good = false;
1500	printf("CUT - can't add\n");
1501	break;
1502	}
1503	}
1504	}
1505	if (good) {
1506	double * cut = new double [numberColumns];
1507	int * which = new int [numberColumns];
1508	double rhs = -1.0;
1509	int n = 0;
1510	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1511	if (isHeuristicInteger(solverC,iColumn)) {
1512	if (lower[iColumn] == upperC[iColumn]) {
1513	rhs += lower[iColumn];
1514	cut[n] = 1.0;
1515	which[n++] = iColumn;
1516	} else if (upper[iColumn] == lowerC[iColumn]) {
1517	rhs -= upper[iColumn];
1518	cut[n] = -1.0;
1519	which[n++] = iColumn;
1520	}
1521	}
1522	}
1523	printf("CUT has %d entries\n", n);
1524	OsiRowCut newCut;
1525	newCut.setLb(-COIN_DBL_MAX);
1526	newCut.setUb(rhs);
1527	newCut.setRow(n, which, cut, false);
1528	model_->makeGlobalCut(newCut);
1529	delete [] cut;
1530	delete [] which;
1531	}
1532	}
1533	#ifdef COIN_DEVELOP
1534	if (status == 1)
1535	printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
1536	else if (status == 2)
1537	printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
1538	#endif
1539	}
1540	if (reset) {
1541	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1542	if (reset[iColumn])
1543	solver->setColLower(iColumn, 0.0);
1544	}
1545	delete [] reset;
1546	}
1547	#ifdef HISTORY_STATISTICS
1548	getHistoryStatistics_ = true;
1549	#endif
1550	solver->setHintParam(OsiDoReducePrint, takeHint, strength);
1551	return returnCode;
1552	}
1553	// Set input solution
1554	void
1555	CbcHeuristic::setInputSolution(const double * solution, double objValue)
1556	{
1557	delete [] inputSolution_;
1558	inputSolution_ = NULL;
1559	if (model_ && solution) {
1560	int numberColumns = model_->getNumCols();
1561	inputSolution_ = new double [numberColumns+1];
1562	memcpy(inputSolution_, solution, numberColumns*sizeof(double));
1563	inputSolution_[numberColumns] = objValue;
1564	}
1565	}
1566
1567	//##############################################################################
1568
1569	inline int compare3BranchingObjects(const CbcBranchingObject* br0,
1570	const CbcBranchingObject* br1)
1571	{
1572	const int t0 = br0->type();
1573	const int t1 = br1->type();
1574	if (t0 < t1) {
1575	return -1;
1576	}
1577	if (t0 > t1) {
1578	return 1;
1579	}
1580	return br0->compareOriginalObject(br1);
1581	}
1582
1583	//==============================================================================
1584
1585	inline bool compareBranchingObjects(const CbcBranchingObject* br0,
1586	const CbcBranchingObject* br1)
1587	{
1588	return compare3BranchingObjects(br0, br1) < 0;
1589	}
1590
1591	//==============================================================================
1592
1593	void
1594	CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
1595	{
1596	// CbcHeurDebugNodes(&model);
1597	CbcNode* node = model.currentNode();
1598	brObj_ = new CbcBranchingObject*[node->depth()];
1599	CbcNodeInfo* nodeInfo = node->nodeInfo();
1600	int cnt = 0;
1601	while (nodeInfo->parentBranch() != NULL) {
1602	const OsiBranchingObject* br = nodeInfo->parentBranch();
1603	const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
1604	if (! cbcbr) {
1605	throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
1606	"gutsOfConstructor",
1607	"CbcHeuristicNode",
1608	__FILE__, __LINE__);
1609	}
1610	brObj_[cnt] = cbcbr->clone();
1611	brObj_[cnt]->previousBranch();
1612	++cnt;
1613	nodeInfo = nodeInfo->parent();
1614	}
1615	std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
1616	if (cnt <= 1) {
1617	numObjects_ = cnt;
1618	} else {
1619	numObjects_ = 0;
1620	CbcBranchingObject* br = NULL; // What should this be?
1621	for (int i = 1; i < cnt; ++i) {
1622	if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
1623	int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
1624	switch (comp) {
1625	case CbcRangeSame: // the same range
1626	case CbcRangeDisjoint: // disjoint decisions
1627	// should not happen! we are on a chain!
1628	abort();
1629	case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
1630	delete brObj_[i];
1631	break;
1632	case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
1633	delete brObj_[numObjects_];
1634	brObj_[numObjects_] = brObj_[i];
1635	break;
1636	case CbcRangeOverlap: // overlap
1637	delete brObj_[i];
1638	delete brObj_[numObjects_];
1639	brObj_[numObjects_] = br;
1640	break;
1641	}
1642	continue;
1643	} else {
1644	brObj_[++numObjects_] = brObj_[i];
1645	}
1646	}
1647	++numObjects_;
1648	}
1649	}
1650
1651	//==============================================================================
1652
1653	CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
1654	{
1655	gutsOfConstructor(model);
1656	}
1657
1658	//==============================================================================
1659
1660	double
1661	CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
1662	{
1663
1664	const double disjointWeight = 1;
1665	const double overlapWeight = 0.4;
1666	const double subsetWeight = 0.2;
1667	int countDisjointWeight = 0;
1668	int countOverlapWeight = 0;
1669	int countSubsetWeight = 0;
1670	int i = 0;
1671	int j = 0;
1672	double dist = 0.0;
1673	#ifdef PRINT_DEBUG
1674	printf(" numObjects_ = %i, node->numObjects_ = %i\n",
1675	numObjects_, node->numObjects_);
1676	#endif
1677	while ( i < numObjects_ && j < node->numObjects_) {
1678	CbcBranchingObject* br0 = brObj_[i];
1679	const CbcBranchingObject* br1 = node->brObj_[j];
1680	#ifdef PRINT_DEBUG
1681	const CbcIntegerBranchingObject* brPrint0 =
1682	dynamic_cast<const CbcIntegerBranchingObject*>(br0);
1683	const double* downBounds = brPrint0->downBounds();
1684	const double* upBounds = brPrint0->upBounds();
1685	int variable = brPrint0->variable();
1686	int way = brPrint0->way();
1687	printf(" br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1688	variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1689	static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1690	const CbcIntegerBranchingObject* brPrint1 =
1691	dynamic_cast<const CbcIntegerBranchingObject*>(br1);
1692	downBounds = brPrint1->downBounds();
1693	upBounds = brPrint1->upBounds();
1694	variable = brPrint1->variable();
1695	way = brPrint1->way();
1696	printf(" br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1697	variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1698	static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1699	#endif
1700	const int brComp = compare3BranchingObjects(br0, br1);
1701	if (brComp < 0) {
1702	dist += subsetWeight;
1703	countSubsetWeight++;
1704	++i;
1705	} else if (brComp > 0) {
1706	dist += subsetWeight;
1707	countSubsetWeight++;
1708	++j;
1709	} else {
1710	const int comp = br0->compareBranchingObject(br1, false);
1711	switch (comp) {
1712	case CbcRangeSame:
1713	// do nothing
1714	break;
1715	case CbcRangeDisjoint: // disjoint decisions
1716	dist += disjointWeight;
1717	countDisjointWeight++;
1718	break;
1719	case CbcRangeSubset: // subset one way or another
1720	case CbcRangeSuperset:
1721	dist += subsetWeight;
1722	countSubsetWeight++;
1723	break;
1724	case CbcRangeOverlap: // overlap
1725	dist += overlapWeight;
1726	countOverlapWeight++;
1727	break;
1728	}
1729	++i;
1730	++j;
1731	}
1732	}
1733	dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
1734	countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
1735	COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
1736	countOverlapWeight, countDisjointWeight));
1737	return dist;
1738	}
1739
1740	//==============================================================================
1741
1742	CbcHeuristicNode::~CbcHeuristicNode()
1743	{
1744	for (int i = 0; i < numObjects_; ++i) {
1745	delete brObj_[i];
1746	}
1747	delete [] brObj_;
1748	}
1749
1750	//==============================================================================
1751
1752	double
1753	CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
1754	{
1755	double minDist = COIN_DBL_MAX;
1756	for (int i = nodeList.size() - 1; i >= 0; --i) {
1757	minDist = CoinMin(minDist, distance(nodeList.node(i)));
1758	}
1759	return minDist;
1760	}
1761
1762	//==============================================================================
1763
1764	bool
1765	CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
1766	const double threshold) const
1767	{
1768	for (int i = nodeList.size() - 1; i >= 0; --i) {
1769	if (distance(nodeList.node(i)) >= threshold) {
1770	continue;
1771	} else {
1772	return true;
1773	}
1774	}
1775	return false;
1776	}
1777
1778	//==============================================================================
1779
1780	double
1781	CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
1782	{
1783	if (nodeList.size() == 0) {
1784	return COIN_DBL_MAX;
1785	}
1786	double sumDist = 0;
1787	for (int i = nodeList.size() - 1; i >= 0; --i) {
1788	sumDist += distance(nodeList.node(i));
1789	}
1790	return sumDist / nodeList.size();
1791	}
1792
1793	//##############################################################################
1794
1795	// Default Constructor
1796	CbcRounding::CbcRounding()
1797	: CbcHeuristic()
1798	{
1799	// matrix and row copy will automatically be empty
1800	seed_ = 7654321;
1801	down_ = NULL;
1802	up_ = NULL;
1803	equal_ = NULL;
1804	//whereFrom_ \|= 16*(1+256); // allow more often
1805	}
1806
1807	// Constructor from model
1808	CbcRounding::CbcRounding(CbcModel & model)
1809	: CbcHeuristic(model)
1810	{
1811	// Get a copy of original matrix (and by row for rounding);
1812	assert(model.solver());
1813	if (model.solver()->getNumRows()) {
1814	matrix_ = *model.solver()->getMatrixByCol();
1815	matrixByRow_ = *model.solver()->getMatrixByRow();
1816	validate();
1817	}
1818	down_ = NULL;
1819	up_ = NULL;
1820	equal_ = NULL;
1821	seed_ = 7654321;
1822	//whereFrom_ \|= 16*(1+256); // allow more often
1823	}
1824
1825	// Destructor
1826	CbcRounding::~CbcRounding ()
1827	{
1828	delete [] down_;
1829	delete [] up_;
1830	delete [] equal_;
1831	}
1832
1833	// Clone
1834	CbcHeuristic *
1835	CbcRounding::clone() const
1836	{
1837	return new CbcRounding(*this);
1838	}
1839	// Create C++ lines to get to current state
1840	void
1841	CbcRounding::generateCpp( FILE * fp)
1842	{
1843	CbcRounding other;
1844	fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
1845	fprintf(fp, "3 CbcRounding rounding(*cbcModel);\n");
1846	CbcHeuristic::generateCpp(fp, "rounding");
1847	if (seed_ != other.seed_)
1848	fprintf(fp, "3 rounding.setSeed(%d);\n", seed_);
1849	else
1850	fprintf(fp, "4 rounding.setSeed(%d);\n", seed_);
1851	fprintf(fp, "3 cbcModel->addHeuristic(&rounding);\n");
1852	}
1853	//#define NEW_ROUNDING
1854	// Copy constructor
1855	CbcRounding::CbcRounding(const CbcRounding & rhs)
1856	:
1857	CbcHeuristic(rhs),
1858	matrix_(rhs.matrix_),
1859	matrixByRow_(rhs.matrixByRow_),
1860	seed_(rhs.seed_)
1861	{
1862	#ifdef NEW_ROUNDING
1863	int numberColumns = matrix_.getNumCols();
1864	down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1865	up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1866	equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1867	#else
1868	down_ = NULL;
1869	up_ = NULL;
1870	equal_ = NULL;
1871	#endif
1872	}
1873
1874	// Assignment operator
1875	CbcRounding &
1876	CbcRounding::operator=( const CbcRounding & rhs)
1877	{
1878	if (this != &rhs) {
1879	CbcHeuristic::operator=(rhs);
1880	matrix_ = rhs.matrix_;
1881	matrixByRow_ = rhs.matrixByRow_;
1882	#ifdef NEW_ROUNDING
1883	delete [] down_;
1884	delete [] up_;
1885	delete [] equal_;
1886	int numberColumns = matrix_.getNumCols();
1887	down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1888	up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1889	equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1890	#else
1891	down_ = NULL;
1892	up_ = NULL;
1893	equal_ = NULL;
1894	#endif
1895	seed_ = rhs.seed_;
1896	}
1897	return *this;
1898	}
1899
1900	// Resets stuff if model changes
1901	void
1902	CbcRounding::resetModel(CbcModel * model)
1903	{
1904	model_ = model;
1905	// Get a copy of original matrix (and by row for rounding);
1906	assert(model_->solver());
1907	matrix_ = *model_->solver()->getMatrixByCol();
1908	matrixByRow_ = *model_->solver()->getMatrixByRow();
1909	validate();
1910	}
1911	/* Check whether the heuristic should run at all
1912	0 - before cuts at root node (or from doHeuristics)
1913	1 - during cuts at root
1914	2 - after root node cuts
1915	3 - after cuts at other nodes
1916	4 - during cuts at other nodes
1917	8 added if previous heuristic in loop found solution
1918	*/
1919	bool
1920	CbcRounding::shouldHeurRun(int whereFrom)
1921	{
1922	if (whereFrom!=4) {
1923	return CbcHeuristic::shouldHeurRun(whereFrom);
1924	} else {
1925	numCouldRun_++;
1926	return shouldHeurRun_randomChoice();
1927	}
1928	}
1929	// See if rounding will give solution
1930	// Sets value of solution
1931	// Assumes rhs for original matrix still okay
1932	// At present only works with integers
1933	// Fix values if asked for
1934	// Returns 1 if solution, 0 if not
1935	int
1936	CbcRounding::solution(double & solutionValue,
1937	double * betterSolution)
1938	{
1939
1940	numCouldRun_++;
1941	// See if to do
1942	if (!when() \|\| (when() % 10 == 1 && model_->phase() != 1) \|\|
1943	(when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1944	return 0; // switched off
1945	numRuns_++;
1946	#ifdef HEURISTIC_INFORM
1947	printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
1948	heuristicName(),numRuns_,numCouldRun_,when_);
1949	#endif
1950	OsiSolverInterface * solver = model_->solver();
1951	double direction = solver->getObjSense();
1952	double newSolutionValue = direction * solver->getObjValue();
1953	return solution(solutionValue, betterSolution, newSolutionValue);
1954	}
1955	// See if rounding will give solution
1956	// Sets value of solution
1957	// Assumes rhs for original matrix still okay
1958	// At present only works with integers
1959	// Fix values if asked for
1960	// Returns 1 if solution, 0 if not
1961	int
1962	CbcRounding::solution(double & solutionValue,
1963	double * betterSolution,
1964	double newSolutionValue)
1965	{
1966
1967	// See if to do
1968	if (!when() \|\| (when() % 10 == 1 && model_->phase() != 1) \|\|
1969	(when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1970	return 0; // switched off
1971	OsiSolverInterface * solver = model_->solver();
1972	const double * lower = solver->getColLower();
1973	const double * upper = solver->getColUpper();
1974	const double * rowLower = solver->getRowLower();
1975	const double * rowUpper = solver->getRowUpper();
1976	const double * solution = solver->getColSolution();
1977	const double * objective = solver->getObjCoefficients();
1978	double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1979	double primalTolerance;
1980	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1981	double useTolerance = primalTolerance;
1982
1983	int numberRows = matrix_.getNumRows();
1984	assert (numberRows <= solver->getNumRows());
1985	if (numberRows == 0){
1986	return 0;
1987	}
1988	int numberIntegers = model_->numberIntegers();
1989	const int * integerVariable = model_->integerVariable();
1990	int i;
1991	double direction = solver->getObjSense();
1992	//double newSolutionValue = direction*solver->getObjValue();
1993	int returnCode = 0;
1994	// Column copy
1995	const double * element = matrix_.getElements();
1996	const int * row = matrix_.getIndices();
1997	const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1998	const int * columnLength = matrix_.getVectorLengths();
1999	// Row copy
2000	const double * elementByRow = matrixByRow_.getElements();
2001	const int * column = matrixByRow_.getIndices();
2002	const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2003	const int * rowLength = matrixByRow_.getVectorLengths();
2004
2005	// Get solution array for heuristic solution
2006	int numberColumns = solver->getNumCols();
2007	double * newSolution = new double [numberColumns];
2008	memcpy(newSolution, solution, numberColumns*sizeof(double));
2009
2010	double * rowActivity = new double[numberRows];
2011	memset(rowActivity, 0, numberRows*sizeof(double));
2012	for (i = 0; i < numberColumns; i++) {
2013	int j;
2014	double value = newSolution[i];
2015	if (value < lower[i]) {
2016	value = lower[i];
2017	newSolution[i] = value;
2018	} else if (value > upper[i]) {
2019	value = upper[i];
2020	newSolution[i] = value;
2021	}
2022	if (value) {
2023	for (j = columnStart[i];
2024	j < columnStart[i] + columnLength[i]; j++) {
2025	int iRow = row[j];
2026	rowActivity[iRow] += value * element[j];
2027	}
2028	}
2029	}
2030	// check was feasible - if not adjust (cleaning may move)
2031	for (i = 0; i < numberRows; i++) {
2032	if (rowActivity[i] < rowLower[i]) {
2033	//assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2034	rowActivity[i] = rowLower[i];
2035	} else if (rowActivity[i] > rowUpper[i]) {
2036	//assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2037	rowActivity[i] = rowUpper[i];
2038	}
2039	}
2040	for (i = 0; i < numberIntegers; i++) {
2041	int iColumn = integerVariable[i];
2042	double value = newSolution[iColumn];
2043	//double thisTolerance = integerTolerance;
2044	if (fabs(floor(value + 0.5) - value) > integerTolerance) {
2045	double below = floor(value);
2046	double newValue = newSolution[iColumn];
2047	double cost = direction * objective[iColumn];
2048	double move;
2049	if (cost > 0.0) {
2050	// try up
2051	move = 1.0 - (value - below);
2052	} else if (cost < 0.0) {
2053	// try down
2054	move = below - value;
2055	} else {
2056	// won't be able to move unless we can grab another variable
2057	double randomNumber = randomNumberGenerator_.randomDouble();
2058	// which way?
2059	if (randomNumber < 0.5)
2060	move = below - value;
2061	else
2062	move = 1.0 - (value - below);
2063	}
2064	newValue += move;
2065	newSolution[iColumn] = newValue;
2066	newSolutionValue += move * cost;
2067	int j;
2068	for (j = columnStart[iColumn];
2069	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2070	int iRow = row[j];
2071	rowActivity[iRow] += move * element[j];
2072	}
2073	}
2074	}
2075
2076	double penalty = 0.0;
2077	// see if feasible - just using singletons
2078	for (i = 0; i < numberRows; i++) {
2079	double value = rowActivity[i];
2080	double thisInfeasibility = 0.0;
2081	if (value < rowLower[i] - primalTolerance)
2082	thisInfeasibility = value - rowLower[i];
2083	else if (value > rowUpper[i] + primalTolerance)
2084	thisInfeasibility = value - rowUpper[i];
2085	if (thisInfeasibility) {
2086	// See if there are any slacks I can use to fix up
2087	// maybe put in coding for multiple slacks?
2088	double bestCost = 1.0e50;
2089	int k;
2090	int iBest = -1;
2091	double addCost = 0.0;
2092	double newValue = 0.0;
2093	double changeRowActivity = 0.0;
2094	double absInfeasibility = fabs(thisInfeasibility);
2095	for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2096	int iColumn = column[k];
2097	// See if all elements help
2098	if (columnLength[iColumn] == 1) {
2099	double currentValue = newSolution[iColumn];
2100	double elementValue = elementByRow[k];
2101	double lowerValue = lower[iColumn];
2102	double upperValue = upper[iColumn];
2103	double gap = rowUpper[i] - rowLower[i];
2104	double absElement = fabs(elementValue);
2105	if (thisInfeasibility*elementValue > 0.0) {
2106	// we want to reduce
2107	if ((currentValue - lowerValue)*absElement >= absInfeasibility) {
2108	// possible - check if integer
2109	double distance = absInfeasibility / absElement;
2110	double thisCost = -direction * objective[iColumn] * distance;
2111	if (isHeuristicInteger(solver,iColumn)) {
2112	distance = ceil(distance - useTolerance);
2113	if (currentValue - distance >= lowerValue - useTolerance) {
2114	if (absInfeasibility - distance*absElement < -gap - useTolerance)
2115	thisCost = 1.0e100; // no good
2116	else
2117	thisCost = -direction * objective[iColumn] * distance;
2118	} else {
2119	thisCost = 1.0e100; // no good
2120	}
2121	}
2122	if (thisCost < bestCost) {
2123	bestCost = thisCost;
2124	iBest = iColumn;
2125	addCost = thisCost;
2126	newValue = currentValue - distance;
2127	changeRowActivity = -distance * elementValue;
2128	}
2129	}
2130	} else {
2131	// we want to increase
2132	if ((upperValue - currentValue)*absElement >= absInfeasibility) {
2133	// possible - check if integer
2134	double distance = absInfeasibility / absElement;
2135	double thisCost = direction * objective[iColumn] * distance;
2136	if (isHeuristicInteger(solver,iColumn)) {
2137	distance = ceil(distance - 1.0e-7);
2138	assert (currentValue - distance <= upperValue + useTolerance);
2139	if (absInfeasibility - distance*absElement < -gap - useTolerance)
2140	thisCost = 1.0e100; // no good
2141	else
2142	thisCost = direction * objective[iColumn] * distance;
2143	}
2144	if (thisCost < bestCost) {
2145	bestCost = thisCost;
2146	iBest = iColumn;
2147	addCost = thisCost;
2148	newValue = currentValue + distance;
2149	changeRowActivity = distance * elementValue;
2150	}
2151	}
2152	}
2153	}
2154	}
2155	if (iBest >= 0) {
2156	/*printf("Infeasibility of %g on row %d cost %g\n",
2157	thisInfeasibility,i,addCost);*/
2158	newSolution[iBest] = newValue;
2159	thisInfeasibility = 0.0;
2160	newSolutionValue += addCost;
2161	rowActivity[i] += changeRowActivity;
2162	}
2163	penalty += fabs(thisInfeasibility);
2164	}
2165	}
2166	if (penalty) {
2167	// see if feasible using any
2168	// first continuous
2169	double penaltyChange = 0.0;
2170	int iColumn;
2171	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2172	if (isHeuristicInteger(solver,iColumn))
2173	continue;
2174	double currentValue = newSolution[iColumn];
2175	double lowerValue = lower[iColumn];
2176	double upperValue = upper[iColumn];
2177	int j;
2178	int anyBadDown = 0;
2179	int anyBadUp = 0;
2180	double upImprovement = 0.0;
2181	double downImprovement = 0.0;
2182	for (j = columnStart[iColumn];
2183	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2184	int iRow = row[j];
2185	if (rowUpper[iRow] > rowLower[iRow]) {
2186	double value = element[j];
2187	if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2188	// infeasible above
2189	downImprovement += value;
2190	upImprovement -= value;
2191	if (value > 0.0)
2192	anyBadUp++;
2193	else
2194	anyBadDown++;
2195	} else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2196	// feasible at ub
2197	if (value > 0.0) {
2198	upImprovement -= value;
2199	anyBadUp++;
2200	} else {
2201	downImprovement += value;
2202	anyBadDown++;
2203	}
2204	} else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2205	// feasible in interior
2206	} else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2207	// feasible at lb
2208	if (value < 0.0) {
2209	upImprovement += value;
2210	anyBadUp++;
2211	} else {
2212	downImprovement -= value;
2213	anyBadDown++;
2214	}
2215	} else {
2216	// infeasible below
2217	downImprovement -= value;
2218	upImprovement += value;
2219	if (value < 0.0)
2220	anyBadUp++;
2221	else
2222	anyBadDown++;
2223	}
2224	} else {
2225	// equality row
2226	double value = element[j];
2227	if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2228	// infeasible above
2229	downImprovement += value;
2230	upImprovement -= value;
2231	if (value > 0.0)
2232	anyBadUp++;
2233	else
2234	anyBadDown++;
2235	} else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2236	// infeasible below
2237	downImprovement -= value;
2238	upImprovement += value;
2239	if (value < 0.0)
2240	anyBadUp++;
2241	else
2242	anyBadDown++;
2243	} else {
2244	// feasible - no good
2245	anyBadUp = -1;
2246	anyBadDown = -1;
2247	break;
2248	}
2249	}
2250	}
2251	// could change tests for anyBad
2252	if (anyBadUp)
2253	upImprovement = 0.0;
2254	if (anyBadDown)
2255	downImprovement = 0.0;
2256	double way = 0.0;
2257	double improvement = 0.0;
2258	if (downImprovement > 0.0 && currentValue > lowerValue) {
2259	way = -1.0;
2260	improvement = downImprovement;
2261	} else if (upImprovement > 0.0 && currentValue < upperValue) {
2262	way = 1.0;
2263	improvement = upImprovement;
2264	}
2265	if (way) {
2266	// can improve
2267	double distance;
2268	if (way > 0.0)
2269	distance = upperValue - currentValue;
2270	else
2271	distance = currentValue - lowerValue;
2272	for (j = columnStart[iColumn];
2273	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2274	int iRow = row[j];
2275	double value = element[j] * way;
2276	if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2277	// infeasible above
2278	assert (value < 0.0);
2279	double gap = rowActivity[iRow] - rowUpper[iRow];
2280	if (gap + value*distance < 0.0)
2281	distance = -gap / value;
2282	} else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2283	// infeasible below
2284	assert (value > 0.0);
2285	double gap = rowActivity[iRow] - rowLower[iRow];
2286	if (gap + value*distance > 0.0)
2287	distance = -gap / value;
2288	} else {
2289	// feasible
2290	if (value > 0) {
2291	double gap = rowActivity[iRow] - rowUpper[iRow];
2292	if (gap + value*distance > 0.0)
2293	distance = -gap / value;
2294	} else {
2295	double gap = rowActivity[iRow] - rowLower[iRow];
2296	if (gap + value*distance < 0.0)
2297	distance = -gap / value;
2298	}
2299	}
2300	}
2301	//move
2302	penaltyChange += improvement * distance;
2303	distance *= way;
2304	newSolution[iColumn] += distance;
2305	newSolutionValue += direction * objective[iColumn] * distance;
2306	for (j = columnStart[iColumn];
2307	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2308	int iRow = row[j];
2309	double value = element[j];
2310	rowActivity[iRow] += distance * value;
2311	}
2312	}
2313	}
2314	// and now all if improving
2315	double lastChange = penaltyChange ? 1.0 : 0.0;
2316	int numberPasses=0;
2317	while (lastChange > 1.0e-2 && numberPasses < 1000) {
2318	lastChange = 0;
2319	numberPasses++;
2320	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2321	bool isInteger = isHeuristicInteger(solver,iColumn);
2322	double currentValue = newSolution[iColumn];
2323	double lowerValue = lower[iColumn];
2324	double upperValue = upper[iColumn];
2325	int j;
2326	int anyBadDown = 0;
2327	int anyBadUp = 0;
2328	double upImprovement = 0.0;
2329	double downImprovement = 0.0;
2330	for (j = columnStart[iColumn];
2331	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2332	int iRow = row[j];
2333	double value = element[j];
2334	if (isInteger) {
2335	if (value > 0.0) {
2336	if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
2337	anyBadUp++;
2338	if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
2339	anyBadDown++;
2340	} else {
2341	if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
2342	anyBadDown++;
2343	if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
2344	anyBadUp++;
2345	}
2346	}
2347	if (rowUpper[iRow] > rowLower[iRow]) {
2348	if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2349	// infeasible above
2350	downImprovement += value;
2351	upImprovement -= value;
2352	if (value > 0.0)
2353	anyBadUp++;
2354	else
2355	anyBadDown++;
2356	} else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2357	// feasible at ub
2358	if (value > 0.0) {
2359	upImprovement -= value;
2360	anyBadUp++;
2361	} else {
2362	downImprovement += value;
2363	anyBadDown++;
2364	}
2365	} else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2366	// feasible in interior
2367	} else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2368	// feasible at lb
2369	if (value < 0.0) {
2370	upImprovement += value;
2371	anyBadUp++;
2372	} else {
2373	downImprovement -= value;
2374	anyBadDown++;
2375	}
2376	} else {
2377	// infeasible below
2378	downImprovement -= value;
2379	upImprovement += value;
2380	if (value < 0.0)
2381	anyBadUp++;
2382	else
2383	anyBadDown++;
2384	}
2385	} else {
2386	// equality row
2387	if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2388	// infeasible above
2389	downImprovement += value;
2390	upImprovement -= value;
2391	if (value > 0.0)
2392	anyBadUp++;
2393	else
2394	anyBadDown++;
2395	} else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2396	// infeasible below
2397	downImprovement -= value;
2398	upImprovement += value;
2399	if (value < 0.0)
2400	anyBadUp++;
2401	else
2402	anyBadDown++;
2403	} else {
2404	// feasible - no good
2405	anyBadUp = -1;
2406	anyBadDown = -1;
2407	break;
2408	}
2409	}
2410	}
2411	// could change tests for anyBad
2412	if (anyBadUp)
2413	upImprovement = 0.0;
2414	if (anyBadDown)
2415	downImprovement = 0.0;
2416	double way = 0.0;
2417	double improvement = 0.0;
2418	if (downImprovement > 0.0 && currentValue > lowerValue) {
2419	way = -1.0;
2420	improvement = downImprovement;
2421	if (isInteger&&currentValue<lowerValue+0.99)
2422	continue; // no good
2423	} else if (upImprovement > 0.0 && currentValue < upperValue) {
2424	way = 1.0;
2425	improvement = upImprovement;
2426	if (isInteger&&currentValue>upperValue-0.99)
2427	continue; // no good
2428	}
2429	if (way) {
2430	// can improve
2431	double distance = COIN_DBL_MAX;
2432	for (j = columnStart[iColumn];
2433	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2434	int iRow = row[j];
2435	double value = element[j] * way;
2436	if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2437	// infeasible above
2438	assert (value < 0.0);
2439	double gap = rowActivity[iRow] - rowUpper[iRow];
2440	if (gap + value*distance < 0.0) {
2441	// If integer then has to move by 1
2442	if (!isInteger)
2443	distance = -gap / value;
2444	else
2445	distance = CoinMax(-gap / value, 1.0);
2446	}
2447	} else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2448	// infeasible below
2449	assert (value > 0.0);
2450	double gap = rowActivity[iRow] - rowLower[iRow];
2451	if (gap + value*distance > 0.0) {
2452	// If integer then has to move by 1
2453	if (!isInteger)
2454	distance = -gap / value;
2455	else
2456	distance = CoinMax(-gap / value, 1.0);
2457	}
2458	} else {
2459	// feasible
2460	if (value > 0) {
2461	double gap = rowActivity[iRow] - rowUpper[iRow];
2462	if (gap + value*distance > 0.0)
2463	distance = -gap / value;
2464	} else {
2465	double gap = rowActivity[iRow] - rowLower[iRow];
2466	if (gap + value*distance < 0.0)
2467	distance = -gap / value;
2468	}
2469	}
2470	}
2471	if (isInteger)
2472	distance = floor(distance + 1.05e-8);
2473	if (!distance) {
2474	// should never happen
2475	//printf("zero distance in CbcRounding - debug\n");
2476	}
2477	//move
2478	lastChange += improvement * distance;
2479	distance *= way;
2480	newSolution[iColumn] += distance;
2481	newSolutionValue += direction * objective[iColumn] * distance;
2482	for (j = columnStart[iColumn];
2483	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2484	int iRow = row[j];
2485	double value = element[j];
2486	rowActivity[iRow] += distance * value;
2487	}
2488	}
2489	}
2490	penaltyChange += lastChange;
2491	}
2492	penalty -= penaltyChange;
2493	if (penalty < 1.0e-5*fabs(penaltyChange)) {
2494	// recompute
2495	penalty = 0.0;
2496	for (i = 0; i < numberRows; i++) {
2497	double value = rowActivity[i];
2498	if (value < rowLower[i] - primalTolerance)
2499	penalty += rowLower[i] - value;
2500	else if (value > rowUpper[i] + primalTolerance)
2501	penalty += value - rowUpper[i];
2502	}
2503	}
2504	}
2505
2506	// Could also set SOS (using random) and repeat
2507	if (!penalty) {
2508	// See if we can do better
2509	//seed_++;
2510	//CoinSeedRandom(seed_);
2511	// Random number between 0 and 1.
2512	double randomNumber = randomNumberGenerator_.randomDouble();
2513	int iPass;
2514	int start[2];
2515	int end[2];
2516	int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers)));
2517	start[0] = iRandom;
2518	end[0] = numberIntegers;
2519	start[1] = 0;
2520	end[1] = iRandom;
2521	for (iPass = 0; iPass < 2; iPass++) {
2522	int i;
2523	for (i = start[iPass]; i < end[iPass]; i++) {
2524	int iColumn = integerVariable[i];
2525	#ifndef NDEBUG
2526	double value = newSolution[iColumn];
2527	assert (fabs(floor(value + 0.5) - value) < integerTolerance);
2528	#endif
2529	double cost = direction * objective[iColumn];
2530	double move = 0.0;
2531	if (cost > 0.0)
2532	move = -1.0;
2533	else if (cost < 0.0)
2534	move = 1.0;
2535	while (move) {
2536	bool good = true;
2537	double newValue = newSolution[iColumn] + move;
2538	if (newValue < lower[iColumn] - useTolerance \|\|
2539	newValue > upper[iColumn] + useTolerance) {
2540	move = 0.0;
2541	} else {
2542	// see if we can move
2543	int j;
2544	for (j = columnStart[iColumn];
2545	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2546	int iRow = row[j];
2547	double newActivity = rowActivity[iRow] + move * element[j];
2548	if (newActivity < rowLower[iRow] - primalTolerance \|\|
2549	newActivity > rowUpper[iRow] + primalTolerance) {
2550	good = false;
2551	break;
2552	}
2553	}
2554	if (good) {
2555	newSolution[iColumn] = newValue;
2556	newSolutionValue += move * cost;
2557	int j;
2558	for (j = columnStart[iColumn];
2559	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2560	int iRow = row[j];
2561	rowActivity[iRow] += move * element[j];
2562	}
2563	} else {
2564	move = 0.0;
2565	}
2566	}
2567	}
2568	}
2569	}
2570	// Just in case of some stupidity
2571	double objOffset = 0.0;
2572	solver->getDblParam(OsiObjOffset, objOffset);
2573	newSolutionValue = -objOffset;
2574	for ( i = 0 ; i < numberColumns ; i++ )
2575	newSolutionValue += objective[i] * newSolution[i];
2576	newSolutionValue *= direction;
2577	//printf("new solution value %g %g\n",newSolutionValue,solutionValue);
2578	if (newSolutionValue < solutionValue) {
2579	// paranoid check
2580	memset(rowActivity, 0, numberRows*sizeof(double));
2581	for (i = 0; i < numberColumns; i++) {
2582	int j;
2583	double value = newSolution[i];
2584	if (value) {
2585	for (j = columnStart[i];
2586	j < columnStart[i] + columnLength[i]; j++) {
2587	int iRow = row[j];
2588	rowActivity[iRow] += value * element[j];
2589	}
2590	}
2591	}
2592	// check was approximately feasible
2593	bool feasible = true;
2594	for (i = 0; i < numberRows; i++) {
2595	if (rowActivity[i] < rowLower[i]) {
2596	if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
2597	feasible = false;
2598	} else if (rowActivity[i] > rowUpper[i]) {
2599	if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
2600	feasible = false;
2601	}
2602	}
2603	if (feasible) {
2604	// new solution
2605	memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2606	solutionValue = newSolutionValue;
2607	//printf("** Solution of %g found by rounding\n",newSolutionValue);
2608	returnCode = 1;
2609	} else {
2610	// Can easily happen
2611	//printf("Debug CbcRounding giving bad solution\n");
2612	}
2613	}
2614	}
2615	#ifdef NEW_ROUNDING
2616	if (!returnCode) {
2617	#ifdef JJF_ZERO
2618	// back to starting point
2619	memcpy(newSolution, solution, numberColumns*sizeof(double));
2620	memset(rowActivity, 0, numberRows*sizeof(double));
2621	for (i = 0; i < numberColumns; i++) {
2622	int j;
2623	double value = newSolution[i];
2624	if (value < lower[i]) {
2625	value = lower[i];
2626	newSolution[i] = value;
2627	} else if (value > upper[i]) {
2628	value = upper[i];
2629	newSolution[i] = value;
2630	}
2631	if (value) {
2632	for (j = columnStart[i];
2633	j < columnStart[i] + columnLength[i]; j++) {
2634	int iRow = row[j];
2635	rowActivity[iRow] += value * element[j];
2636	}
2637	}
2638	}
2639	// check was feasible - if not adjust (cleaning may move)
2640	for (i = 0; i < numberRows; i++) {
2641	if (rowActivity[i] < rowLower[i]) {
2642	//assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2643	rowActivity[i] = rowLower[i];
2644	} else if (rowActivity[i] > rowUpper[i]) {
2645	//assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2646	rowActivity[i] = rowUpper[i];
2647	}
2648	}
2649	#endif
2650	int * candidate = new int [numberColumns];
2651	int nCandidate = 0;
2652	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2653	bool isInteger = isHeuristicInteger(solver,iColumn);
2654	if (isInteger) {
2655	double currentValue = newSolution[iColumn];
2656	if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
2657	candidate[nCandidate++] = iColumn;
2658	}
2659	}
2660	if (true) {
2661	// Rounding as in Berthold
2662	while (nCandidate) {
2663	double infeasibility = 1.0e-7;
2664	int iRow = -1;
2665	for (i = 0; i < numberRows; i++) {
2666	double value = 0.0;
2667	if (rowActivity[i] < rowLower[i]) {
2668	value = rowLower[i] - rowActivity[i];
2669	} else if (rowActivity[i] > rowUpper[i]) {
2670	value = rowActivity[i] - rowUpper[i];
2671	}
2672	if (value > infeasibility) {
2673	infeasibility = value;
2674	iRow = i;
2675	}
2676	}
2677	if (iRow >= 0) {
2678	// infeasible
2679	} else {
2680	// feasible
2681	}
2682	}
2683	} else {
2684	// Shifting as in Berthold
2685	}
2686	delete [] candidate;
2687	}
2688	#endif
2689	delete [] newSolution;
2690	delete [] rowActivity;
2691	return returnCode;
2692	}
2693	// update model
2694	void CbcRounding::setModel(CbcModel * model)
2695	{
2696	model_ = model;
2697	// Get a copy of original matrix (and by row for rounding);
2698	assert(model_->solver());
2699	if (model_->solver()->getNumRows()) {
2700	matrix_ = *model_->solver()->getMatrixByCol();
2701	matrixByRow_ = *model_->solver()->getMatrixByRow();
2702	// make sure model okay for heuristic
2703	validate();
2704	}
2705	}
2706	// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2707	void
2708	CbcRounding::validate()
2709	{
2710	if (model_ && (when() % 100) < 10) {
2711	if (model_->numberIntegers() !=
2712	model_->numberObjects() && (model_->numberObjects() \|\|
2713	(model_->specialOptions()&1024) == 0)) {
2714	int numberOdd = 0;
2715	for (int i = 0; i < model_->numberObjects(); i++) {
2716	if (!model_->object(i)->canDoHeuristics())
2717	numberOdd++;
2718	}
2719	if (numberOdd)
2720	setWhen(0);
2721	}
2722	}
2723	#ifdef NEW_ROUNDING
2724	int numberColumns = matrix_.getNumCols();
2725	down_ = new unsigned short [numberColumns];
2726	up_ = new unsigned short [numberColumns];
2727	equal_ = new unsigned short [numberColumns];
2728	// Column copy
2729	const double * element = matrix_.getElements();
2730	const int * row = matrix_.getIndices();
2731	const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2732	const int * columnLength = matrix_.getVectorLengths();
2733	const double * rowLower = model.solver()->getRowLower();
2734	const double * rowUpper = model.solver()->getRowUpper();
2735	for (int i = 0; i < numberColumns; i++) {
2736	int down = 0;
2737	int up = 0;
2738	int equal = 0;
2739	if (columnLength[i] > 65535) {
2740	equal[0] = 65535;
2741	break; // unlikely to work
2742	}
2743	for (CoinBigIndex j = columnStart[i];
2744	j < columnStart[i] + columnLength[i]; j++) {
2745	int iRow = row[j];
2746	if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
2747	equal++;
2748	} else if (element[j] > 0.0) {
2749	if (rowUpper[iRow] < 1.0e20)
2750	up++;
2751	else
2752	down--;
2753	} else {
2754	if (rowLower[iRow] > -1.0e20)
2755	up++;
2756	else
2757	down--;
2758	}
2759	}
2760	down_[i] = (unsigned short) down;
2761	up_[i] = (unsigned short) up;
2762	equal_[i] = (unsigned short) equal;
2763	}
2764	#else
2765	down_ = NULL;
2766	up_ = NULL;
2767	equal_ = NULL;
2768	#endif
2769	}
2770
2771	// Default Constructor
2772	CbcHeuristicPartial::CbcHeuristicPartial()
2773	: CbcHeuristic()
2774	{
2775	fixPriority_ = 10000;
2776	}
2777
2778	// Constructor from model
2779	CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
2780	: CbcHeuristic(model)
2781	{
2782	fixPriority_ = fixPriority;
2783	setNumberNodes(numberNodes);
2784	validate();
2785	}
2786
2787	// Destructor
2788	CbcHeuristicPartial::~CbcHeuristicPartial ()
2789	{
2790	}
2791
2792	// Clone
2793	CbcHeuristic *
2794	CbcHeuristicPartial::clone() const
2795	{
2796	return new CbcHeuristicPartial(*this);
2797	}
2798	// Create C++ lines to get to current state
2799	void
2800	CbcHeuristicPartial::generateCpp( FILE * fp)
2801	{
2802	CbcHeuristicPartial other;
2803	fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2804	fprintf(fp, "3 CbcHeuristicPartial partial(*cbcModel);\n");
2805	CbcHeuristic::generateCpp(fp, "partial");
2806	if (fixPriority_ != other.fixPriority_)
2807	fprintf(fp, "3 partial.setFixPriority(%d);\n", fixPriority_);
2808	else
2809	fprintf(fp, "4 partial.setFixPriority(%d);\n", fixPriority_);
2810	fprintf(fp, "3 cbcModel->addHeuristic(&partial);\n");
2811	}
2812	//#define NEW_PARTIAL
2813	// Copy constructor
2814	CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
2815	:
2816	CbcHeuristic(rhs),
2817	fixPriority_(rhs.fixPriority_)
2818	{
2819	}
2820
2821	// Assignment operator
2822	CbcHeuristicPartial &
2823	CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs)
2824	{
2825	if (this != &rhs) {
2826	CbcHeuristic::operator=(rhs);
2827	fixPriority_ = rhs.fixPriority_;
2828	}
2829	return *this;
2830	}
2831
2832	// Resets stuff if model changes
2833	void
2834	CbcHeuristicPartial::resetModel(CbcModel * model)
2835	{
2836	model_ = model;
2837	// Get a copy of original matrix (and by row for partial);
2838	assert(model_->solver());
2839	validate();
2840	}
2841	// See if partial will give solution
2842	// Sets value of solution
2843	// Assumes rhs for original matrix still okay
2844	// At present only works with integers
2845	// Fix values if asked for
2846	// Returns 1 if solution, 0 if not
2847	int
2848	CbcHeuristicPartial::solution(double & solutionValue,
2849	double * betterSolution)
2850	{
2851	// Return if already done
2852	if (fixPriority_ < 0)
2853	return 0; // switched off
2854	#ifdef HEURISTIC_INFORM
2855	printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
2856	heuristicName(),numRuns_,numCouldRun_,when_);
2857	#endif
2858	const double * hotstartSolution = model_->hotstartSolution();
2859	const int * hotstartPriorities = model_->hotstartPriorities();
2860	if (!hotstartSolution)
2861	return 0;
2862	OsiSolverInterface * solver = model_->solver();
2863
2864	int numberIntegers = model_->numberIntegers();
2865	const int * integerVariable = model_->integerVariable();
2866
2867	OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
2868	const double * colLower = newSolver->getColLower();
2869	const double * colUpper = newSolver->getColUpper();
2870
2871	double primalTolerance;
2872	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
2873
2874	int i;
2875	int numberFixed = 0;
2876	int returnCode = 0;
2877
2878	for (i = 0; i < numberIntegers; i++) {
2879	int iColumn = integerVariable[i];
2880	if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
2881	double value = hotstartSolution[iColumn];
2882	double lower = colLower[iColumn];
2883	double upper = colUpper[iColumn];
2884	value = CoinMax(value, lower);
2885	value = CoinMin(value, upper);
2886	if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
2887	value = floor(value + 0.5);
2888	newSolver->setColLower(iColumn, value);
2889	newSolver->setColUpper(iColumn, value);
2890	numberFixed++;
2891	}
2892	}
2893	}
2894	if (numberFixed > numberIntegers / 5 - 100000000) {
2895	#ifdef COIN_DEVELOP
2896	printf("%d integers fixed\n", numberFixed);
2897	#endif
2898	returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
2899	model_->getCutoff(), "CbcHeuristicPartial");
2900	if (returnCode < 0)
2901	returnCode = 0; // returned on size
2902	//printf("return code %d",returnCode);
2903	if ((returnCode&2) != 0) {
2904	// could add cut
2905	returnCode &= ~2;
2906	//printf("could add cut with %d elements (if all 0-1)\n",nFix);
2907	} else {
2908	//printf("\n");
2909	}
2910	}
2911	fixPriority_ = -1; // switch off
2912
2913	delete newSolver;
2914	return returnCode;
2915	}
2916	// update model
2917	void CbcHeuristicPartial::setModel(CbcModel * model)
2918	{
2919	model_ = model;
2920	assert(model_->solver());
2921	// make sure model okay for heuristic
2922	validate();
2923	}
2924	// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2925	void
2926	CbcHeuristicPartial::validate()
2927	{
2928	if (model_ && (when() % 100) < 10) {
2929	if (model_->numberIntegers() !=
2930	model_->numberObjects())
2931	setWhen(0);
2932	}
2933	}
2934	bool
2935	CbcHeuristicPartial::shouldHeurRun(int /whereFrom/)
2936	{
2937	return true;
2938	}
2939
2940	// Default Constructor
2941	CbcSerendipity::CbcSerendipity()
2942	: CbcHeuristic()
2943	{
2944	}
2945
2946	// Constructor from model
2947	CbcSerendipity::CbcSerendipity(CbcModel & model)
2948	: CbcHeuristic(model)
2949	{
2950	}
2951
2952	// Destructor
2953	CbcSerendipity::~CbcSerendipity ()
2954	{
2955	}
2956
2957	// Clone
2958	CbcHeuristic *
2959	CbcSerendipity::clone() const
2960	{
2961	return new CbcSerendipity(*this);
2962	}
2963	// Create C++ lines to get to current state
2964	void
2965	CbcSerendipity::generateCpp( FILE * fp)
2966	{
2967	fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2968	fprintf(fp, "3 CbcSerendipity serendipity(*cbcModel);\n");
2969	CbcHeuristic::generateCpp(fp, "serendipity");
2970	fprintf(fp, "3 cbcModel->addHeuristic(&serendipity);\n");
2971	}
2972
2973	// Copy constructor
2974	CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
2975	:
2976	CbcHeuristic(rhs)
2977	{
2978	}
2979
2980	// Assignment operator
2981	CbcSerendipity &
2982	CbcSerendipity::operator=( const CbcSerendipity & rhs)
2983	{
2984	if (this != &rhs) {
2985	CbcHeuristic::operator=(rhs);
2986	}
2987	return *this;
2988	}
2989
2990	// Returns 1 if solution, 0 if not
2991	int
2992	CbcSerendipity::solution(double & solutionValue,
2993	double * betterSolution)
2994	{
2995	if (!model_)
2996	return 0;
2997	#ifdef HEURISTIC_INFORM
2998	printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
2999	heuristicName(),numRuns_,numCouldRun_,when_);
3000	#endif
3001	if (!inputSolution_) {
3002	// get information on solver type
3003	OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
3004	OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
3005	if (auxiliaryInfo) {
3006	return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols());
3007	} else {
3008	return 0;
3009	}
3010	} else {
3011	int numberColumns = model_->getNumCols();
3012	double value = inputSolution_[numberColumns];
3013	int returnCode = 0;
3014	if (value < solutionValue) {
3015	solutionValue = value;
3016	memcpy(betterSolution, inputSolution_, numberColumns*sizeof(double));
3017	returnCode = 1;
3018	}
3019	delete [] inputSolution_;
3020	inputSolution_ = NULL;
3021	model_ = NULL; // switch off
3022	return returnCode;
3023	}
3024	}
3025	// update model
3026	void CbcSerendipity::setModel(CbcModel * model)
3027	{
3028	model_ = model;
3029	}
3030	// Resets stuff if model changes
3031	void
3032	CbcSerendipity::resetModel(CbcModel * model)
3033	{
3034	model_ = model;
3035	}
3036
3037
3038	// Default Constructor
3039	CbcHeuristicJustOne::CbcHeuristicJustOne()
3040	: CbcHeuristic(),
3041	probabilities_(NULL),
3042	heuristic_(NULL),
3043	numberHeuristics_(0)
3044	{
3045	}
3046
3047	// Constructor from model
3048	CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel & model)
3049	: CbcHeuristic(model),
3050	probabilities_(NULL),
3051	heuristic_(NULL),
3052	numberHeuristics_(0)
3053	{
3054	}
3055
3056	// Destructor
3057	CbcHeuristicJustOne::~CbcHeuristicJustOne ()
3058	{
3059	for (int i = 0; i < numberHeuristics_; i++)
3060	delete heuristic_[i];
3061	delete [] heuristic_;
3062	delete [] probabilities_;
3063	}
3064
3065	// Clone
3066	CbcHeuristicJustOne *
3067	CbcHeuristicJustOne::clone() const
3068	{
3069	return new CbcHeuristicJustOne(*this);
3070	}
3071
3072	// Create C++ lines to get to current state
3073	void
3074	CbcHeuristicJustOne::generateCpp( FILE * fp)
3075	{
3076	CbcHeuristicJustOne other;
3077	fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n");
3078	fprintf(fp, "3 CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n");
3079	CbcHeuristic::generateCpp(fp, "heuristicJustOne");
3080	fprintf(fp, "3 cbcModel->addHeuristic(&heuristicJustOne);\n");
3081	}
3082
3083	// Copy constructor
3084	CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne & rhs)
3085	:
3086	CbcHeuristic(rhs),
3087	probabilities_(NULL),
3088	heuristic_(NULL),
3089	numberHeuristics_(rhs.numberHeuristics_)
3090	{
3091	if (numberHeuristics_) {
3092	probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
3093	heuristic_ = new CbcHeuristic * [numberHeuristics_];
3094	for (int i = 0; i < numberHeuristics_; i++)
3095	heuristic_[i] = rhs.heuristic_[i]->clone();
3096	}
3097	}
3098
3099	// Assignment operator
3100	CbcHeuristicJustOne &
3101	CbcHeuristicJustOne::operator=( const CbcHeuristicJustOne & rhs)
3102	{
3103	if (this != &rhs) {
3104	CbcHeuristic::operator=(rhs);
3105	for (int i = 0; i < numberHeuristics_; i++)
3106	delete heuristic_[i];
3107	delete [] heuristic_;
3108	delete [] probabilities_;
3109	probabilities_ = NULL;
3110	heuristic_ = NULL;
3111	numberHeuristics_ = rhs.numberHeuristics_;
3112	if (numberHeuristics_) {
3113	probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
3114	heuristic_ = new CbcHeuristic * [numberHeuristics_];
3115	for (int i = 0; i < numberHeuristics_; i++)
3116	heuristic_[i] = rhs.heuristic_[i]->clone();
3117	}
3118	}
3119	return *this;
3120	}
3121	// Sets value of solution
3122	// Returns 1 if solution, 0 if not
3123	int
3124	CbcHeuristicJustOne::solution(double & solutionValue,
3125	double * betterSolution)
3126	{
3127	#ifdef DIVE_DEBUG
3128	std::cout << "solutionValue = " << solutionValue << std::endl;
3129	#endif
3130	++numCouldRun_;
3131
3132	// test if the heuristic can run
3133	if (!shouldHeurRun_randomChoice() \|\| !numberHeuristics_)
3134	return 0;
3135	double randomNumber = randomNumberGenerator_.randomDouble();
3136	int i;
3137	for (i = 0; i < numberHeuristics_; i++) {
3138	if (randomNumber < probabilities_[i])
3139	break;
3140	}
3141	assert (i < numberHeuristics_);
3142	int returnCode;
3143	//model_->unsetDivingHasRun();
3144	#ifdef COIN_DEVELOP
3145	printf("JustOne running %s\n",
3146	heuristic_[i]->heuristicName());
3147	#endif
3148	returnCode = heuristic_[i]->solution(solutionValue, betterSolution);
3149	#ifdef COIN_DEVELOP
3150	if (returnCode)
3151	printf("JustOne running %s found solution\n",
3152	heuristic_[i]->heuristicName());
3153	#endif
3154	return returnCode;
3155	}
3156	// Resets stuff if model changes
3157	void
3158	CbcHeuristicJustOne::resetModel(CbcModel * model)
3159	{
3160	CbcHeuristic::resetModel(model);
3161	for (int i = 0; i < numberHeuristics_; i++)
3162	heuristic_[i]->resetModel(model);
3163	}
3164	// update model (This is needed if cliques update matrix etc)
3165	void
3166	CbcHeuristicJustOne::setModel(CbcModel * model)
3167	{
3168	CbcHeuristic::setModel(model);
3169	for (int i = 0; i < numberHeuristics_; i++)
3170	heuristic_[i]->setModel(model);
3171	}
3172	// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
3173	void
3174	CbcHeuristicJustOne::validate()
3175	{
3176	CbcHeuristic::validate();
3177	for (int i = 0; i < numberHeuristics_; i++)
3178	heuristic_[i]->validate();
3179	}
3180	// Adds an heuristic with probability
3181	void
3182	CbcHeuristicJustOne::addHeuristic(const CbcHeuristic * heuristic, double probability)
3183	{
3184	CbcHeuristic * thisOne = heuristic->clone();
3185	thisOne->setWhen(-999);
3186	CbcHeuristic ** tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1,
3187	numberHeuristics_);
3188	delete [] heuristic_;
3189	heuristic_ = tempH;
3190	heuristic_[numberHeuristics_] = thisOne;
3191	double * tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1,
3192	numberHeuristics_);
3193	delete [] probabilities_;
3194	probabilities_ = tempP;
3195	probabilities_[numberHeuristics_] = probability;
3196	numberHeuristics_++;
3197	}
3198	// Normalize probabilities
3199	void
3200	CbcHeuristicJustOne::normalizeProbabilities()
3201	{
3202	double sum = 0.0;
3203	for (int i = 0; i < numberHeuristics_; i++)
3204	sum += probabilities_[i];
3205	double multiplier = 1.0 / sum;
3206	sum = 0.0;
3207	for (int i = 0; i < numberHeuristics_; i++) {
3208	sum += probabilities_[i];
3209	probabilities_[i] = sum * multiplier;
3210	}
3211	assert (fabs(probabilities_[numberHeuristics_-1] - 1.0) < 1.0e-5);
3212	probabilities_[numberHeuristics_-1] = 1.000001;
3213	}
3214

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

Download in other formats: