Context Navigation

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

Last change on this file since 1587 was 1587, checked in by forrest, 10 years ago
minor changes to heuristics and fix a memory leak
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 111.8 KB

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

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

Download in other formats: