Context Navigation

CbcModel.cpp @ 2351

Last change on this file since 2351 was 2351, checked in by forrest, 3 years ago
fix when obviously infeasible at root node
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 758.1 KB

Line
1	/* $Id: CbcModel.cpp 2351 2017-12-19 10:08:16Z 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 <string>
14	//#define CBC_DEBUG 1
15	//#define CHECK_CUT_COUNTS
16	//#define CHECK_NODE
17	//#define CHECK_NODE_FULL
18	//#define NODE_LOG
19	//#define GLOBAL_CUTS_JUST_POINTERS
20	#ifdef CGL_DEBUG_GOMORY
21	extern int gomory_try;
22	#endif
23	#include <cassert>
24	#include <cmath>
25	#include <cfloat>
26	#ifdef COIN_HAS_CLP
27	// include Presolve from Clp
28	#include "ClpPresolve.hpp"
29	#include "OsiClpSolverInterface.hpp"
30	#include "ClpNode.hpp"
31	#include "ClpDualRowDantzig.hpp"
32	#include "ClpSimplexPrimal.hpp"
33	#endif
34
35	#include "CbcEventHandler.hpp"
36
37	#include "OsiSolverInterface.hpp"
38	#include "OsiAuxInfo.hpp"
39	#include "OsiSolverBranch.hpp"
40	#include "OsiChooseVariable.hpp"
41	#include "CoinWarmStartBasis.hpp"
42	#include "CoinPackedMatrix.hpp"
43	#include "CoinHelperFunctions.hpp"
44	#include "CbcBranchActual.hpp"
45	#include "CbcBranchDynamic.hpp"
46	#include "CbcHeuristic.hpp"
47	#include "CbcHeuristicFPump.hpp"
48	#include "CbcHeuristicRINS.hpp"
49	#include "CbcHeuristicDive.hpp"
50	#include "CbcModel.hpp"
51	#include "CbcTreeLocal.hpp"
52	#include "CbcStatistics.hpp"
53	#include "CbcStrategy.hpp"
54	#include "CbcMessage.hpp"
55	#include "OsiRowCut.hpp"
56	#include "OsiColCut.hpp"
57	#include "OsiRowCutDebugger.hpp"
58	#include "OsiCuts.hpp"
59	#include "CbcCountRowCut.hpp"
60	#include "CbcCutGenerator.hpp"
61	#include "CbcFeasibilityBase.hpp"
62	#include "CbcFathom.hpp"
63	#include "CbcFullNodeInfo.hpp"
64	#ifdef COIN_HAS_NTY
65	#include "CbcSymmetry.hpp"
66	#endif
67	// include Probing
68	#include "CglProbing.hpp"
69	#include "CglGomory.hpp"
70	#include "CglTwomir.hpp"
71	// include preprocessing
72	#include "CglPreProcess.hpp"
73	#include "CglDuplicateRow.hpp"
74	#include "CglStored.hpp"
75	#include "CglClique.hpp"
76	#include "CglKnapsackCover.hpp"
77
78	#include "CoinTime.hpp"
79	#include "CoinMpsIO.hpp"
80
81	#include "CbcCompareActual.hpp"
82	#include "CbcTree.hpp"
83	// This may be dummy
84	#include "CbcThread.hpp"
85	/* Various functions local to CbcModel.cpp */
86
87	static void * doRootCbcThread(void * voidInfo);
88
89	namespace {
90
91	//-------------------------------------------------------------------
92	// Returns the greatest common denominator of two
93	// positive integers, a and b, found using Euclid's algorithm
94	//-------------------------------------------------------------------
95	static int gcd(int a, int b)
96	{
97	int remainder = -1;
98	// make sure a<=b (will always remain so)
99	if (a > b) {
100	// Swap a and b
101	int temp = a;
102	a = b;
103	b = temp;
104	}
105	// if zero then gcd is nonzero (zero may occur in rhs of packed)
106	if (!a) {
107	if (b) {
108	return b;
109	} else {
110	printf("**** gcd given two zeros!!\n");
111	abort();
112	}
113	}
114	while (remainder) {
115	remainder = b % a;
116	b = a;
117	a = remainder;
118	}
119	return b;
120	}
121
122
123
124	#ifdef CHECK_NODE_FULL
125
126	/*
127	Routine to verify that tree linkage is correct. The invariant that is tested
128	is
129
130	reference count = (number of actual references) + (number of branches left)
131
132	The routine builds a set of paired arrays, info and count, by traversing the
133	tree. Each CbcNodeInfo is recorded in info, and the number of times it is
134	referenced (via the parent field) is recorded in count. Then a final check is
135	made to see if the numberPointingToThis_ field agrees.
136	*/
137
138	void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
139
140	{
141	if (model.getNodeCount() == 661) return;
142	printf("*** CHECKING tree after %d nodes\n", model.getNodeCount()) ;
143
144	int j ;
145	int nNodes = branchingTree->size() ;
146	# define MAXINFO 1000
147	int *count = new int [MAXINFO] ;
148	CbcNodeInfo *info = new CbcNodeInfo[MAXINFO] ;
149	int nInfo = 0 ;
150	/*
151	Collect all CbcNodeInfo objects in info, by starting from each live node and
152	traversing back to the root. Nodes in the live set should have unexplored
153	branches remaining.
154
155	TODO: The `while (nodeInfo)' loop could be made to break on reaching a
156	common ancester (nodeInfo is found in info[k]). Alternatively, the
157	check could change to signal an error if nodeInfo is not found above a
158	common ancestor.
159	*/
160	for (j = 0 ; j < nNodes ; j++) {
161	CbcNode *node = branchingTree->nodePointer(j) ;
162	if (!node)
163	continue;
164	CbcNodeInfo *nodeInfo = node->nodeInfo() ;
165	int change = node->nodeInfo()->numberBranchesLeft() ;
166	assert(change) ;
167	while (nodeInfo) {
168	int k ;
169	for (k = 0 ; k < nInfo ; k++) {
170	if (nodeInfo == info[k]) break ;
171	}
172	if (k == nInfo) {
173	assert(nInfo < MAXINFO) ;
174	nInfo++ ;
175	info[k] = nodeInfo ;
176	count[k] = 0 ;
177	}
178	nodeInfo = nodeInfo->parent() ;
179	}
180	}
181	/*
182	Walk the info array. For each nodeInfo, look up its parent in info and
183	increment the corresponding count.
184	*/
185	for (j = 0 ; j < nInfo ; j++) {
186	CbcNodeInfo *nodeInfo = info[j] ;
187	nodeInfo = nodeInfo->parent() ;
188	if (nodeInfo) {
189	int k ;
190	for (k = 0 ; k < nInfo ; k++) {
191	if (nodeInfo == info[k]) break ;
192	}
193	assert (k < nInfo) ;
194	count[k]++ ;
195	}
196	}
197	/*
198	Walk the info array one more time and check that the invariant holds. The
199	number of references (numberPointingToThis()) should equal the sum of the
200	number of actual references (held in count[]) plus the number of potential
201	references (unexplored branches, numberBranchesLeft()).
202	*/
203	for (j = 0; j < nInfo; j++) {
204	CbcNodeInfo * nodeInfo = info[j] ;
205	if (nodeInfo) {
206	int k ;
207	for (k = 0; k < nInfo; k++)
208	if (nodeInfo == info[k])
209	break ;
210	printf("Nodeinfo %x - %d left, %d count\n",
211	nodeInfo,
212	nodeInfo->numberBranchesLeft(),
213	nodeInfo->numberPointingToThis()) ;
214	assert(nodeInfo->numberPointingToThis() ==
215	count[k] + nodeInfo->numberBranchesLeft()) ;
216	}
217	}
218
219	delete [] count ;
220	delete [] info ;
221
222	return ;
223	}
224
225	#endif /* CHECK_NODE_FULL */
226
227
228
229	#ifdef CHECK_CUT_COUNTS
230
231	/*
232	Routine to verify that cut reference counts are correct.
233	*/
234	void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
235
236	{
237	printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount()) ;
238
239	int j ;
240	int nNodes = branchingTree->size() ;
241
242	/*
243	cut.tempNumber_ exists for the purpose of doing this verification. Clear it
244	in all cuts. We traverse the tree by starting from each live node and working
245	back to the root. At each CbcNodeInfo, check for cuts.
246	*/
247	for (j = 0 ; j < nNodes ; j++) {
248	CbcNode *node = branchingTree->nodePointer(j) ;
249	CbcNodeInfo * nodeInfo = node->nodeInfo() ;
250	assert (node->nodeInfo()->numberBranchesLeft()) ;
251	while (nodeInfo) {
252	int k ;
253	for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
254	CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
255	if (cut) cut->tempNumber_ = 0;
256	}
257	nodeInfo = nodeInfo->parent() ;
258	}
259	}
260	/*
261	Walk the live set again, this time collecting the list of cuts in use at each
262	node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
263	that when we recreate the basis for a node, we compress out the slack cuts.
264	*/
265	for (j = 0 ; j < nNodes ; j++) {
266	CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
267	CbcNode *node = branchingTree->nodePointer(j) ;
268	CbcNodeInfo *nodeInfo = node->nodeInfo();
269	int change = node->nodeInfo()->numberBranchesLeft() ;
270	printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
271	node->nodeInfo(), node->columnNumber(), node->way(),
272	node->objectiveValue()) ;
273
274	model.addCuts1(node, debugws) ;
275
276	int i ;
277	int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
278	CbcCountRowCut **addedCuts = model.addedCuts() ;
279	for (i = 0 ; i < model.currentNumberCuts() ; i++) {
280	CoinWarmStartBasis::Status status =
281	debugws->getArtifStatus(i + numberRowsAtContinuous) ;
282	if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
283	addedCuts[i]->tempNumber_ += change ;
284	}
285	}
286
287	while (nodeInfo) {
288	nodeInfo = nodeInfo->parent() ;
289	if (nodeInfo) printf(" -> %x", nodeInfo);
290	}
291	printf("\n") ;
292	delete debugws ;
293	}
294	/*
295	The moment of truth: We've tallied up the references by direct scan of the search tree. Check for agreement with the count in the cut.
296
297	TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
298	*/
299	for (j = 0 ; j < nNodes ; j++) {
300	CbcNode *node = branchingTree->nodePointer(j) ;
301	CbcNodeInfo *nodeInfo = node->nodeInfo();
302	while (nodeInfo) {
303	int k ;
304	for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
305	CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
306	if (cut && cut->tempNumber_ >= 0) {
307	if (cut->tempNumber_ != cut->numberPointingToThis())
308	printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
309	cut, cut->tempNumber_, cut->numberPointingToThis()) ;
310	else
311	printf(" match %x %d %x %d %d\n", nodeInfo, k,
312	cut, cut->tempNumber_, cut->numberPointingToThis()) ;
313	cut->tempNumber_ = -1 ;
314	}
315	}
316	nodeInfo = nodeInfo->parent() ;
317	}
318	}
319
320	return ;
321	}
322
323	#endif /* CHECK_CUT_COUNTS */
324
325
326	#ifdef CHECK_CUT_SIZE
327
328	/*
329	Routine to verify that cut reference counts are correct.
330	*/
331	void verifyCutSize (const CbcTree * branchingTree, CbcModel &model)
332	{
333
334	int j ;
335	int nNodes = branchingTree->size() ;
336	int totalCuts = 0;
337
338	/*
339	cut.tempNumber_ exists for the purpose of doing this verification. Clear it
340	in all cuts. We traverse the tree by starting from each live node and working
341	back to the root. At each CbcNodeInfo, check for cuts.
342	*/
343	for (j = 0 ; j < nNodes ; j++) {
344	CbcNode *node = branchingTree->nodePointer(j) ;
345	CbcNodeInfo * nodeInfo = node->nodeInfo() ;
346	assert (node->nodeInfo()->numberBranchesLeft()) ;
347	while (nodeInfo) {
348	totalCuts += nodeInfo->numberCuts();
349	nodeInfo = nodeInfo->parent() ;
350	}
351	}
352	printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts) ;
353	return ;
354	}
355
356	#endif /* CHECK_CUT_SIZE */
357
358	}
359
360	/* End unnamed namespace for CbcModel.cpp */
361
362
363	void
364	CbcModel::analyzeObjective ()
365	/*
366	Try to find a minimum change in the objective function. The first scan
367	checks that there are no continuous variables with non-zero coefficients,
368	and grabs the largest objective coefficient associated with an unfixed
369	integer variable. The second scan attempts to scale up the objective
370	coefficients to a point where they are sufficiently close to integer that
371	we can pretend they are integer, and calculate a gcd over the coefficients
372	of interest. This will be the minimum increment for the scaled coefficients.
373	The final action is to scale the increment back for the original coefficients
374	and install it, if it's better than the existing value.
375
376	John's note: We could do better than this.
377
378	John's second note - apologies for changing s to z
379	*/
380	{
381	const double *objective = getObjCoefficients() ;
382	const double *lower = getColLower() ;
383	const double *upper = getColUpper() ;
384	/*
385	Scan continuous and integer variables to see if continuous
386	are cover or network with integral rhs.
387	*/
388	double continuousMultiplier = 1.0;
389	double * coeffMultiplier = NULL;
390	double largestObj = 0.0;
391	double smallestObj = COIN_DBL_MAX;
392	{
393	const double *rowLower = getRowLower() ;
394	const double *rowUpper = getRowUpper() ;
395	int numberRows = solver_->getNumRows() ;
396	double * rhs = new double [numberRows];
397	memset(rhs, 0, numberRows*sizeof(double));
398	int iColumn;
399	int numberColumns = solver_->getNumCols() ;
400	// Column copy of matrix
401	int problemType = -1;
402	const double * element = solver_->getMatrixByCol()->getElements();
403	const int * row = solver_->getMatrixByCol()->getIndices();
404	const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
405	const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
406	int numberInteger = 0;
407	int numberIntegerObj = 0;
408	int numberGeneralIntegerObj = 0;
409	int numberIntegerWeight = 0;
410	int numberContinuousObj = 0;
411	double cost = COIN_DBL_MAX;
412	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
413	if (upper[iColumn] == lower[iColumn]) {
414	CoinBigIndex start = columnStart[iColumn];
415	CoinBigIndex end = start + columnLength[iColumn];
416	for (CoinBigIndex j = start; j < end; j++) {
417	int iRow = row[j];
418	rhs[iRow] += lower[iColumn] * element[j];
419	}
420	} else {
421	double objValue = objective[iColumn];
422	if (solver_->isInteger(iColumn))
423	numberInteger++;
424	if (objValue) {
425	if (!solver_->isInteger(iColumn)) {
426	numberContinuousObj++;
427	} else {
428	largestObj = CoinMax(largestObj, fabs(objValue));
429	smallestObj = CoinMin(smallestObj, fabs(objValue));
430	numberIntegerObj++;
431	if (cost == COIN_DBL_MAX)
432	cost = objValue;
433	else if (cost != objValue)
434	cost = -COIN_DBL_MAX;
435	int gap = static_cast<int> (upper[iColumn] - lower[iColumn]);
436	if (gap > 1) {
437	numberGeneralIntegerObj++;
438	numberIntegerWeight += gap;
439	}
440	}
441	}
442	}
443	}
444	int iType = 0;
445	if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 &&
446	numberIntegerObj*3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
447	iType = 1 + 4 + (((moreSpecialOptions_&536870912)==0) ? 2 : 0);
448	else if (!numberContinuousObj && numberIntegerObj <= 100 &&
449	numberIntegerObj*5 < numberObjects_ && numberIntegerWeight <= 100 &&
450	!parentModel_ &&
451	solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
452	iType = 4 + (((moreSpecialOptions_&536870912)==0) ? 2 : 0);
453	else if (!numberContinuousObj && numberIntegerObj <= 100 &&
454	numberIntegerObj*5 < numberObjects_ &&
455	!parentModel_ &&
456	solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
457	iType = 8;
458	int iTest = getMaximumNodes();
459	if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) {
460	iType = iTest - 987654320;
461	printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g - %d continuous\n",
462	numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj);
463	if (iType == 9)
464	exit(77);
465	if (numberContinuousObj)
466	iType = 0;
467	}
468
469	//if (!numberContinuousObj&&(numberIntegerObj<=5\|\|cost!=-COIN_DBL_MAX)&&
470	//numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_->getNumRows()>100) {
471	if (iType) {
472	/*
473	A) put high priority on (if none)
474	B) create artificial objective (if clp)
475	*/
476	int iPriority = -1;
477	for (int i = 0; i < numberObjects_; i++) {
478	int k = object_[i]->priority();
479	if (iPriority == -1)
480	iPriority = k;
481	else if (iPriority != k)
482	iPriority = -2;
483	}
484	bool branchOnSatisfied = ((iType & 1) != 0);
485	bool createFake = ((iType & 2) != 0);
486	bool randomCost = ((iType & 4) != 0);
487	if (iPriority >= 0) {
488	char general[200];
489	if (cost == -COIN_DBL_MAX) {
490	sprintf(general, "%d integer variables out of %d objects (%d integer) have costs - high priority",
491	numberIntegerObj, numberObjects_, numberInteger);
492	} else if (cost == COIN_DBL_MAX) {
493	sprintf(general, "No integer variables out of %d objects (%d integer) have costs",
494	numberObjects_, numberInteger);
495	branchOnSatisfied = false;
496	} else {
497	sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g - high priority",
498	numberIntegerObj, numberObjects_, numberInteger, cost);
499	}
500	messageHandler()->message(CBC_GENERAL,
501	messages())
502	<< general << CoinMessageEol ;
503	sprintf(general, "branch on satisfied %c create fake objective %c random cost %c",
504	branchOnSatisfied ? 'Y' : 'N',
505	createFake ? 'Y' : 'N',
506	randomCost ? 'Y' : 'N');
507	messageHandler()->message(CBC_GENERAL,
508	messages())
509	<< general << CoinMessageEol ;
510	// switch off clp type branching
511	// no ? fastNodeDepth_ = -1;
512	int highPriority = (branchOnSatisfied) ? -999 : 100;
513	for (int i = 0; i < numberObjects_; i++) {
514	CbcSimpleInteger * thisOne = dynamic_cast <CbcSimpleInteger *> (object_[i]);
515	object_[i]->setPriority(1000);
516	if (thisOne) {
517	int iColumn = thisOne->columnNumber();
518	if (objective[iColumn])
519	thisOne->setPriority(highPriority);
520	}
521	}
522	}
523	#ifdef COIN_HAS_CLP
524	OsiClpSolverInterface * clpSolver
525	= dynamic_cast<OsiClpSolverInterface *> (solver_);
526	if (clpSolver && createFake) {
527	// Create artificial objective to be used when all else fixed
528	int numberColumns = clpSolver->getNumCols();
529	double * fakeObj = new double [numberColumns];
530	// Column copy
531	const CoinPackedMatrix * matrixByCol = clpSolver->getMatrixByCol();
532	//const double * element = matrixByCol.getElements();
533	//const int * row = matrixByCol.getIndices();
534	//const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
535	const int * columnLength = matrixByCol->getVectorLengths();
536	const double * solution = clpSolver->getColSolution();
537	#ifdef JJF_ZERO
538	int nAtBound = 0;
539	for (int i = 0; i < numberColumns; i++) {
540	double lowerValue = lower[i];
541	double upperValue = upper[i];
542	if (clpSolver->isInteger(i)) {
543	double lowerValue = lower[i];
544	double upperValue = upper[i];
545	double value = solution[i];
546	if (value < lowerValue + 1.0e-6 \|\|
547	value > upperValue - 1.0e-6)
548	nAtBound++;
549	}
550	}
551	#endif
552	/*
553	Generate a random objective function for problems where the given objective
554	function is not terribly useful. (Nearly feasible, single integer variable,
555	that sort of thing.
556	*/
557	CoinDrand48(true, 1234567);
558	for (int i = 0; i < numberColumns; i++) {
559	double lowerValue = lower[i];
560	double upperValue = upper[i];
561	double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
562	: i + 1 + columnLength[i] * 1000;
563	value *= 0.001;
564	//value += columnLength[i];
565	if (lowerValue > -1.0e5 \|\| upperValue < 1.0e5) {
566	if (fabs(lowerValue) > fabs(upperValue))
567	value = - value;
568	if (clpSolver->isInteger(i)) {
569	double solValue = solution[i];
570	// Better to add in 0.5 or 1.0??
571	if (solValue < lowerValue + 1.0e-6)
572	value = fabs(value) + 0.5; //fabs(value*1.5);
573	else if (solValue > upperValue - 1.0e-6)
574	value = -fabs(value) - 0.5; //-fabs(value*1.5);
575	}
576	} else {
577	value = 0.0;
578	}
579	fakeObj[i] = value;
580	}
581	// pass to solver
582	clpSolver->setFakeObjective(fakeObj);
583	delete [] fakeObj;
584	}
585	#endif
586	} else if (largestObj < smallestObj*5.0 && !parentModel_ &&
587	!numberContinuousObj &&
588	!numberGeneralIntegerObj &&
589	numberIntegerObj*2 < CoinMin(numberColumns,20)) {
590	// up priorities on costed
591	int iPriority = -1;
592	for (int i = 0; i < numberObjects_; i++) {
593	int k = object_[i]->priority();
594	if (iPriority == -1)
595	iPriority = k;
596	else if (iPriority != k)
597	iPriority = -2;
598	}
599	if (iPriority >= 100) {
600	#if CBC_USEFUL_PRINTING>1
601	printf("Setting variables with obj to high priority\n");
602	#endif
603	for (int i = 0; i < numberObjects_; i++) {
604	CbcSimpleInteger * obj =
605	dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
606	if (obj) {
607	int iColumn = obj->columnNumber();
608	if (objective[iColumn])
609	object_[i]->setPriority(iPriority - 1);
610	}
611	}
612	}
613	}
614	int iRow;
615	for (iRow = 0; iRow < numberRows; iRow++) {
616	if (rowLower[iRow] > -1.0e20 &&
617	fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
618	continuousMultiplier = 0.0;
619	break;
620	}
621	if (rowUpper[iRow] < 1.0e20 &&
622	fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
623	continuousMultiplier = 0.0;
624	break;
625	}
626	// set rhs to limiting value
627	if (rowLower[iRow] != rowUpper[iRow]) {
628	if (rowLower[iRow] > -1.0e20) {
629	if (rowUpper[iRow] < 1.0e20) {
630	// no good
631	continuousMultiplier = 0.0;
632	break;
633	} else {
634	rhs[iRow] = rowLower[iRow] - rhs[iRow];
635	if (problemType < 0)
636	problemType = 3; // set cover
637	else if (problemType != 3)
638	problemType = 4;
639	}
640	} else {
641	rhs[iRow] = rowUpper[iRow] - rhs[iRow];
642	if (problemType < 0)
643	problemType = 1; // set partitioning <=
644	else if (problemType != 1)
645	problemType = 4;
646	}
647	} else {
648	rhs[iRow] = rowUpper[iRow] - rhs[iRow];
649	if (problemType < 0)
650	problemType = 3; // set partitioning ==
651	else if (problemType != 2)
652	problemType = 2;
653	}
654	if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
655	problemType = 4;
656	}
657	if (continuousMultiplier) {
658	// 1 network, 2 cover, 4 negative cover
659	int possible = 7;
660	bool unitRhs = true;
661	// See which rows could be set cover
662	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
663	if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
664	CoinBigIndex start = columnStart[iColumn];
665	CoinBigIndex end = start + columnLength[iColumn];
666	for (CoinBigIndex j = start; j < end; j++) {
667	double value = element[j];
668	if (value == 1.0) {
669	} else if (value == -1.0) {
670	rhs[row[j]] = -0.5;
671	} else {
672	rhs[row[j]] = -COIN_DBL_MAX;
673	}
674	}
675	}
676	}
677	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
678	if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
679	if (!isInteger(iColumn)) {
680	CoinBigIndex start = columnStart[iColumn];
681	CoinBigIndex end = start + columnLength[iColumn];
682	double rhsValue = 0.0;
683	// 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
684	int type = 0;
685	for (CoinBigIndex j = start; j < end; j++) {
686	double value = element[j];
687	if (fabs(value) != 1.0) {
688	type = 3;
689	break;
690	} else if (value == 1.0) {
691	if (!type)
692	type = 1;
693	else if (type != 1)
694	type = 2;
695	} else {
696	if (!type)
697	type = -1;
698	else if (type != -1)
699	type = 2;
700	}
701	int iRow = row[j];
702	if (rhs[iRow] == -COIN_DBL_MAX) {
703	type = 3;
704	break;
705	} else if (rhs[iRow] == -0.5) {
706	// different values
707	unitRhs = false;
708	} else if (rhsValue) {
709	if (rhsValue != rhs[iRow])
710	unitRhs = false;
711	} else {
712	rhsValue = rhs[iRow];
713	}
714	}
715	// if no elements OK
716	if (type == 3) {
717	// no good
718	possible = 0;
719	break;
720	} else if (type == 2) {
721	if (end - start > 2) {
722	// no good
723	possible = 0;
724	break;
725	} else {
726	// only network
727	possible &= 1;
728	if (!possible)
729	break;
730	}
731	} else if (type == 1) {
732	// only cover
733	possible &= 2;
734	if (!possible)
735	break;
736	} else if (type == -1) {
737	// only negative cover
738	possible &= 4;
739	if (!possible)
740	break;
741	}
742	}
743	}
744	}
745	if ((possible == 2 \|\| possible == 4) && !unitRhs) {
746	#if COIN_DEVELOP>1
747	printf("XXXXXX Continuous all +1 but different rhs\n");
748	#endif
749	possible = 0;
750	}
751	// may be all integer
752	if (possible != 7) {
753	if (!possible)
754	continuousMultiplier = 0.0;
755	else if (possible == 1)
756	continuousMultiplier = 1.0;
757	else
758	continuousMultiplier = 0.0; // 0.5 was incorrect;
759	#if COIN_DEVELOP>1
760	if (continuousMultiplier)
761	printf("XXXXXX multiplier of %g\n", continuousMultiplier);
762	#endif
763	if (continuousMultiplier == 0.5) {
764	coeffMultiplier = new double [numberColumns];
765	bool allOne = true;
766	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
767	coeffMultiplier[iColumn] = 1.0;
768	if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
769	if (!isInteger(iColumn)) {
770	CoinBigIndex start = columnStart[iColumn];
771	int iRow = row[start];
772	double value = rhs[iRow];
773	assert (value >= 0.0);
774	if (value != 0.0 && value != 1.0)
775	allOne = false;
776	coeffMultiplier[iColumn] = 0.5 * value;
777	}
778	}
779	}
780	if (allOne) {
781	// back to old way
782	delete [] coeffMultiplier;
783	coeffMultiplier = NULL;
784	}
785	}
786	} else {
787	// all integer
788	problemType_ = problemType;
789	#if COIN_DEVELOP>1
790	printf("Problem type is %d\n", problemType_);
791	#endif
792	}
793	}
794
795	// But try again
796	if (continuousMultiplier < 1.0) {
797	memset(rhs, 0, numberRows*sizeof(double));
798	int * count = new int [numberRows];
799	memset(count, 0, numberRows*sizeof(int));
800	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
801	CoinBigIndex start = columnStart[iColumn];
802	CoinBigIndex end = start + columnLength[iColumn];
803	if (upper[iColumn] == lower[iColumn]) {
804	for (CoinBigIndex j = start; j < end; j++) {
805	int iRow = row[j];
806	rhs[iRow] += lower[iColumn] * element[j];
807	}
808	} else if (solver_->isInteger(iColumn)) {
809	for (CoinBigIndex j = start; j < end; j++) {
810	int iRow = row[j];
811	if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
812	rhs[iRow] = COIN_DBL_MAX;
813	}
814	} else {
815	for (CoinBigIndex j = start; j < end; j++) {
816	int iRow = row[j];
817	count[iRow]++;
818	if (fabs(element[j]) != 1.0)
819	rhs[iRow] = COIN_DBL_MAX;
820	}
821	}
822	}
823	// now look at continuous
824	bool allGood = true;
825	double direction = solver_->getObjSense() ;
826	int numberObj = 0;
827	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
828	if (upper[iColumn] > lower[iColumn]) {
829	double objValue = objective[iColumn] * direction;
830	if (objValue && !solver_->isInteger(iColumn)) {
831	numberObj++;
832	CoinBigIndex start = columnStart[iColumn];
833	CoinBigIndex end = start + columnLength[iColumn];
834	if (objValue > 0.0) {
835	// wants to be as low as possible
836	if (lower[iColumn] < -1.0e10 \|\| fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
837	allGood = false;
838	break;
839	} else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
840	allGood = false;
841	break;
842	}
843	bool singletonRow = true;
844	bool equality = false;
845	for (CoinBigIndex j = start; j < end; j++) {
846	int iRow = row[j];
847	if (count[iRow] > 1)
848	singletonRow = false;
849	else if (rowLower[iRow] == rowUpper[iRow])
850	equality = true;
851	double rhsValue = rhs[iRow];
852	double lowerValue = rowLower[iRow];
853	double upperValue = rowUpper[iRow];
854	if (rhsValue < 1.0e20) {
855	if (lowerValue > -1.0e20)
856	lowerValue -= rhsValue;
857	if (upperValue < 1.0e20)
858	upperValue -= rhsValue;
859	}
860	if (fabs(rhsValue) > 1.0e20 \|\| fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
861	\|\| fabs(element[j]) != 1.0) {
862	// no good
863	allGood = false;
864	break;
865	}
866	if (element[j] > 0.0) {
867	if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
868	// no good
869	allGood = false;
870	break;
871	}
872	} else {
873	if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
874	// no good
875	allGood = false;
876	break;
877	}
878	}
879	}
880	if (!singletonRow && end > start + 1 && !equality)
881	allGood = false;
882	if (!allGood)
883	break;
884	} else {
885	// wants to be as high as possible
886	if (upper[iColumn] > 1.0e10 \|\| fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
887	allGood = false;
888	break;
889	} else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
890	allGood = false;
891	break;
892	}
893	bool singletonRow = true;
894	bool equality = false;
895	for (CoinBigIndex j = start; j < end; j++) {
896	int iRow = row[j];
897	if (count[iRow] > 1)
898	singletonRow = false;
899	else if (rowLower[iRow] == rowUpper[iRow])
900	equality = true;
901	double rhsValue = rhs[iRow];
902	double lowerValue = rowLower[iRow];
903	double upperValue = rowUpper[iRow];
904	if (rhsValue < 1.0e20) {
905	if (lowerValue > -1.0e20)
906	lowerValue -= rhsValue;
907	if (upperValue < 1.0e20)
908	upperValue -= rhsValue;
909	}
910	if (fabs(rhsValue) > 1.0e20 \|\| fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
911	\|\| fabs(element[j]) != 1.0) {
912	// no good
913	allGood = false;
914	break;
915	}
916	if (element[j] < 0.0) {
917	if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
918	// no good
919	allGood = false;
920	break;
921	}
922	} else {
923	if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
924	// no good
925	allGood = false;
926	break;
927	}
928	}
929	}
930	if (!singletonRow && end > start + 1 && !equality)
931	allGood = false;
932	if (!allGood)
933	break;
934	}
935	}
936	}
937	}
938	delete [] count;
939	if (allGood) {
940	#if COIN_DEVELOP>1
941	if (numberObj)
942	printf("YYYY analysis says all continuous with costs will be integer\n");
943	#endif
944	continuousMultiplier = 1.0;
945	}
946	}
947	delete [] rhs;
948	}
949	/*
950	Take a first scan to see if there are unfixed continuous variables in the
951	objective. If so, the minimum objective change could be arbitrarily small.
952	Also pick off the maximum coefficient of an unfixed integer variable.
953
954	If the objective is found to contain only integer variables, set the
955	fathoming discipline to strict.
956	*/
957	double maximumCost = 0.0 ;
958	//double trueIncrement=0.0;
959	int iColumn ;
960	int numberColumns = getNumCols() ;
961	double scaleFactor = 1.0; // due to rhs etc
962	/*
963	Original model did not have integer bounds.
964	*/
965	if ((specialOptions_&65536) == 0) {
966	/* be on safe side (later look carefully as may be able to
967	to get 0.5 say if bounds are multiples of 0.5 */
968	for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
969	if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
970	double value;
971	value = fabs(lower[iColumn]);
972	if (floor(value + 0.5) != value) {
973	scaleFactor = CoinMin(scaleFactor, 0.5);
974	if (floor(2.0value + 0.5) != 2.0value) {
975	scaleFactor = CoinMin(scaleFactor, 0.25);
976	if (floor(4.0value + 0.5) != 4.0value) {
977	scaleFactor = 0.0;
978	}
979	}
980	}
981	value = fabs(upper[iColumn]);
982	if (floor(value + 0.5) != value) {
983	scaleFactor = CoinMin(scaleFactor, 0.5);
984	if (floor(2.0value + 0.5) != 2.0value) {
985	scaleFactor = CoinMin(scaleFactor, 0.25);
986	if (floor(4.0value + 0.5) != 4.0value) {
987	scaleFactor = 0.0;
988	}
989	}
990	}
991	}
992	}
993	}
994	bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0 ;
995	if (possibleMultiple) {
996	for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
997	if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
998	maximumCost = CoinMax(maximumCost, fabs(objective[iColumn])) ;
999	}
1000	}
1001	}
1002	setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple) ;
1003	/*
1004	If a nontrivial increment is possible, try and figure it out. We're looking
1005	for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
1006	variables. Since the c<j> might not be integers, try and inflate them
1007	sufficiently that they look like integers (and we'll deflate the gcd
1008	later).
1009
1010	2520.0 is used as it is a nice multiple of 2,3,5,7
1011	*/
1012	if (possibleMultiple && maximumCost) {
1013	int increment = 0 ;
1014	double multiplier = 2520.0 ;
1015	while (10.0multipliermaximumCost < 1.0e8)
1016	multiplier *= 10.0 ;
1017	int bigIntegers = 0; // Count of large costs which are integer
1018	for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
1019	if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
1020	double objValue = fabs(objective[iColumn]);
1021	if (!isInteger(iColumn)) {
1022	if (!coeffMultiplier)
1023	objValue *= continuousMultiplier;
1024	else
1025	objValue *= coeffMultiplier[iColumn];
1026	}
1027	if (objValue) {
1028	double value = objValue * multiplier ;
1029	if (value < 2.1e9) {
1030	int nearest = static_cast<int> (floor(value + 0.5)) ;
1031	if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
1032	increment = 0 ;
1033	break ;
1034	} else if (!increment) {
1035	increment = nearest ;
1036	} else {
1037	increment = gcd(increment, nearest) ;
1038	}
1039	} else {
1040	// large value - may still be multiple of 1.0
1041	if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
1042	increment = 0;
1043	break;
1044	} else {
1045	bigIntegers++;
1046	}
1047	}
1048	}
1049	}
1050	}
1051	delete [] coeffMultiplier;
1052	/*
1053	If the increment beats the current value for objective change, install it.
1054	*/
1055	if (increment) {
1056	double value = increment ;
1057	double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
1058	if (bigIntegers) {
1059	// allow for 1.0
1060	increment = gcd(increment, static_cast<int> (multiplier));
1061	value = increment;
1062	}
1063	value /= multiplier ;
1064	value *= scaleFactor;
1065	//trueIncrement=CoinMax(cutoff,value);;
1066	if (value*0.999 > cutoff) {
1067	messageHandler()->message(CBC_INTEGERINCREMENT,
1068	messages())
1069	<< value << CoinMessageEol ;
1070	setDblParam(CbcModel::CbcCutoffIncrement, CoinMax(value*0.999,value-1.0e-4)) ;
1071	}
1072	}
1073	}
1074
1075	return ;
1076	}
1077
1078	/*
1079	saveModel called (carved out of) BranchandBound
1080	*/
1081	void CbcModel::saveModel(OsiSolverInterface * saveSolver, double * checkCutoffForRestart, bool * feasible)
1082	{
1083	if (saveSolver && (specialOptions_&32768) != 0) {
1084	// See if worth trying reduction
1085	*checkCutoffForRestart = getCutoff();
1086	bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
1087	(*checkCutoffForRestart < 1.0e20);
1088	int numberColumns = getNumCols();
1089	if (tryNewSearch) {
1090	#if CBC_USEFUL_PRINTING>1
1091	printf("after %d nodes, cutoff %g - looking\n",
1092	numberNodes_, getCutoff());
1093	#endif
1094	saveSolver->resolve();
1095	double direction = saveSolver->getObjSense() ;
1096	double gap = checkCutoffForRestart - saveSolver->getObjValue() direction ;
1097	double tolerance;
1098	saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
1099	if (gap <= 0.0)
1100	gap = tolerance;
1101	gap += 100.0 * tolerance;
1102	double integerTolerance = getDblParam(CbcIntegerTolerance) ;
1103
1104	const double *lower = saveSolver->getColLower() ;
1105	const double *upper = saveSolver->getColUpper() ;
1106	const double *solution = saveSolver->getColSolution() ;
1107	const double *reducedCost = saveSolver->getReducedCost() ;
1108
1109	int numberFixed = 0 ;
1110	int numberFixed2 = 0;
1111	for (int i = 0 ; i < numberIntegers_ ; i++) {
1112	int iColumn = integerVariable_[i] ;
1113	double djValue = direction * reducedCost[iColumn] ;
1114	if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1115	if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1116	saveSolver->setColUpper(iColumn, lower[iColumn]) ;
1117	numberFixed++ ;
1118	} else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1119	saveSolver->setColLower(iColumn, upper[iColumn]) ;
1120	numberFixed++ ;
1121	}
1122	} else {
1123	numberFixed2++;
1124	}
1125	}
1126	#ifdef COIN_DEVELOP
1127	/*
1128	We're debugging. (specialOptions 1)
1129	*/
1130	if ((specialOptions_&1) != 0) {
1131	const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
1132	if (debugger) {
1133	printf("Contains optimal\n") ;
1134	OsiSolverInterface * temp = saveSolver->clone();
1135	const double * solution = debugger->optimalSolution();
1136	const double *lower = temp->getColLower() ;
1137	const double *upper = temp->getColUpper() ;
1138	int n = temp->getNumCols();
1139	for (int i = 0; i < n; i++) {
1140	if (temp->isInteger(i)) {
1141	double value = floor(solution[i] + 0.5);
1142	assert (value >= lower[i] && value <= upper[i]);
1143	temp->setColLower(i, value);
1144	temp->setColUpper(i, value);
1145	}
1146	}
1147	temp->writeMps("reduced_fix");
1148	delete temp;
1149	saveSolver->writeMps("reduced");
1150	} else {
1151	abort();
1152	}
1153	}
1154	printf("Restart could fix %d integers (%d already fixed)\n",
1155	numberFixed + numberFixed2, numberFixed2);
1156	#endif
1157	numberFixed += numberFixed2;
1158	if (numberFixed*20 < numberColumns)
1159	tryNewSearch = false;
1160	}
1161	if (tryNewSearch) {
1162	// back to solver without cuts?
1163	OsiSolverInterface * solver2 = continuousSolver_->clone();
1164	const double *lower = saveSolver->getColLower() ;
1165	const double *upper = saveSolver->getColUpper() ;
1166	for (int i = 0 ; i < numberIntegers_ ; i++) {
1167	int iColumn = integerVariable_[i] ;
1168	solver2->setColLower(iColumn, lower[iColumn]);
1169	solver2->setColUpper(iColumn, upper[iColumn]);
1170	}
1171	// swap
1172	delete saveSolver;
1173	saveSolver = solver2;
1174	double * newSolution = new double[numberColumns];
1175	double objectiveValue = *checkCutoffForRestart;
1176	CbcSerendipity heuristic(*this);
1177	if (bestSolution_)
1178	heuristic.setInputSolution(bestSolution_, bestObjective_);
1179	heuristic.setFractionSmall(0.9);
1180	heuristic.setFeasibilityPumpOptions(1008013);
1181	// Use numberNodes to say how many are original rows
1182	heuristic.setNumberNodes(continuousSolver_->getNumRows());
1183	#ifdef COIN_DEVELOP
1184	if (continuousSolver_->getNumRows() <
1185	saveSolver->getNumRows())
1186	printf("%d rows added ZZZZZ\n",
1187	solver_->getNumRows() - continuousSolver_->getNumRows());
1188	#endif
1189	int returnCode = heuristic.smallBranchAndBound(saveSolver,
1190	-1, newSolution,
1191	objectiveValue,
1192	*checkCutoffForRestart, "Reduce");
1193	if (returnCode < 0) {
1194	#ifdef COIN_DEVELOP
1195	printf("Restart - not small enough to do search after fixing\n");
1196	#endif
1197	delete [] newSolution;
1198	} else {
1199	if ((returnCode&1) != 0) {
1200	// increment number of solutions so other heuristics can test
1201	numberSolutions_++;
1202	numberHeuristicSolutions_++;
1203	lastHeuristic_ = NULL;
1204	setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
1205	}
1206	delete [] newSolution;
1207	*feasible = false; // stop search
1208	}
1209	#if 0 // probably not needed def CBC_THREAD
1210	if (master_) {
1211	lockThread();
1212	if (parallelMode() > 0) {
1213	while (master_->waitForThreadsInTree(0)) {
1214	lockThread();
1215	double dummyBest;
1216	tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
1217	//unlockThread();
1218	}
1219	}
1220	master_->waitForThreadsInTree(2);
1221	delete master_;
1222	master_ = NULL;
1223	masterThread_ = NULL;
1224	}
1225	#endif
1226	}
1227	}
1228	}
1229	/*
1230	Adds integers, called from BranchandBound()
1231	*/
1232	void CbcModel::AddIntegers()
1233	{
1234	int numberColumns = continuousSolver_->getNumCols();
1235	int numberRows = continuousSolver_->getNumRows();
1236	int numberOriginalIntegers = numberIntegers_;
1237	int * del = new int [CoinMax(numberColumns, numberRows)];
1238	int * original = new int [numberColumns];
1239	char * possibleRow = new char [numberRows];
1240	{
1241	const CoinPackedMatrix * rowCopy = continuousSolver_->getMatrixByRow();
1242	const int * column = rowCopy->getIndices();
1243	const int * rowLength = rowCopy->getVectorLengths();
1244	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1245	const double * rowLower = continuousSolver_->getRowLower();
1246	const double * rowUpper = continuousSolver_->getRowUpper();
1247	const double * element = rowCopy->getElements();
1248	for (int i = 0; i < numberRows; i++) {
1249	int nLeft = 0;
1250	bool possible = false;
1251	if (rowLower[i] < -1.0e20) {
1252	double value = rowUpper[i];
1253	if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1254	possible = true;
1255	} else if (rowUpper[i] > 1.0e20) {
1256	double value = rowLower[i];
1257	if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1258	possible = true;
1259	} else {
1260	double value = rowUpper[i];
1261	if (rowLower[i] == rowUpper[i] &&
1262	fabs(value - floor(value + 0.5)) < 1.0e-8)
1263	possible = true;
1264	}
1265	double allSame = (possible) ? 0.0 : -1.0;
1266	for (CoinBigIndex j = rowStart[i];
1267	j < rowStart[i] + rowLength[i]; j++) {
1268	int iColumn = column[j];
1269	if (continuousSolver_->isInteger(iColumn)) {
1270	if (fabs(element[j]) != 1.0)
1271	possible = false;
1272	} else {
1273	nLeft++;
1274	if (!allSame) {
1275	allSame = fabs(element[j]);
1276	} else if (allSame>0.0) {
1277	if (allSame!=fabs(element[j]))
1278	allSame = -1.0;
1279	}
1280	}
1281	}
1282	if (nLeft == rowLength[i] && allSame > 0.0)
1283	possibleRow[i] = 2;
1284	else if (possible \|\| !nLeft)
1285	possibleRow[i] = 1;
1286	else
1287	possibleRow[i] = 0;
1288	}
1289	}
1290	int nDel = 0;
1291	for (int i = 0; i < numberColumns; i++) {
1292	original[i] = i;
1293	if (continuousSolver_->isInteger(i))
1294	del[nDel++] = i;
1295	}
1296	{
1297	// we must not exclude current best solution (rounding errors)
1298	// also not if large values
1299	const int * row = continuousSolver_->getMatrixByCol()->getIndices();
1300	const CoinBigIndex * columnStart = continuousSolver_->getMatrixByCol()->getVectorStarts();
1301	const int * columnLength = continuousSolver_->getMatrixByCol()->getVectorLengths();
1302	const double * solution = continuousSolver_->getColSolution();
1303	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1304	if (!continuousSolver_->isInteger(iColumn)) {
1305	double value = bestSolution_ ? bestSolution_[iColumn] : 0.0;
1306	double value2 = solution[iColumn];
1307	if (fabs(value-floor(value+0.5))>1.0e-8 \|\|
1308	fabs(value2)>1.0e3) {
1309	CoinBigIndex start = columnStart[iColumn];
1310	CoinBigIndex end = start + columnLength[iColumn];
1311	for (CoinBigIndex j = start; j < end; j++) {
1312	int iRow = row[j];
1313	possibleRow[iRow]=0;
1314	}
1315	}
1316	}
1317	}
1318	}
1319	int nExtra = 0;
1320	OsiSolverInterface * copy1 = continuousSolver_->clone();
1321	int nPass = 0;
1322	while (nDel && nPass < 10) {
1323	nPass++;
1324	OsiSolverInterface * copy2 = copy1->clone();
1325	int nLeft = 0;
1326	for (int i = 0; i < nDel; i++)
1327	original[del[i]] = -1;
1328	for (int i = 0; i < numberColumns; i++) {
1329	int kOrig = original[i];
1330	if (kOrig >= 0)
1331	original[nLeft++] = kOrig;
1332	}
1333	assert (nLeft == numberColumns - nDel);
1334	copy2->deleteCols(nDel, del);
1335	numberColumns = copy2->getNumCols();
1336	const CoinPackedMatrix * rowCopy = copy2->getMatrixByRow();
1337	numberRows = rowCopy->getNumRows();
1338	const int * column = rowCopy->getIndices();
1339	const int * rowLength = rowCopy->getVectorLengths();
1340	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1341	const double * rowLower = copy2->getRowLower();
1342	const double * rowUpper = copy2->getRowUpper();
1343	const double * element = rowCopy->getElements();
1344	const CoinPackedMatrix * columnCopy = copy2->getMatrixByCol();
1345	const int * columnLength = columnCopy->getVectorLengths();
1346	nDel = 0;
1347	// Could do gcd stuff on ones with costs
1348	for (int i = 0; i < numberRows; i++) {
1349	if (!rowLength[i]) {
1350	del[nDel++] = i;
1351	} else if (possibleRow[i]) {
1352	if (rowLength[i] == 1) {
1353	int k = rowStart[i];
1354	int iColumn = column[k];
1355	if (!copy2->isInteger(iColumn)) {
1356	double mult = 1.0 / fabs(element[k]);
1357	if (rowLower[i] < -1.0e20) {
1358	// treat rhs as multiple of 1 unless elements all same
1359	double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
1360	if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1361	del[nDel++] = i;
1362	if (columnLength[iColumn] == 1) {
1363	copy2->setInteger(iColumn);
1364	int kOrig = original[iColumn];
1365	setOptionalInteger(kOrig);
1366	}
1367	}
1368	} else if (rowUpper[i] > 1.0e20) {
1369	// treat rhs as multiple of 1 unless elements all same
1370	double value = ((possibleRow[i]==2) ? rowLower[i] : 1.0) * mult;
1371	if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1372	del[nDel++] = i;
1373	if (columnLength[iColumn] == 1) {
1374	copy2->setInteger(iColumn);
1375	int kOrig = original[iColumn];
1376	setOptionalInteger(kOrig);
1377	}
1378	}
1379	} else {
1380	// treat rhs as multiple of 1 unless elements all same
1381	double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
1382	if (rowLower[i] == rowUpper[i] &&
1383	fabs(value - floor(value + 0.5)) < 1.0e-8) {
1384	del[nDel++] = i;
1385	copy2->setInteger(iColumn);
1386	int kOrig = original[iColumn];
1387	setOptionalInteger(kOrig);
1388	}
1389	}
1390	}
1391	} else {
1392	// only if all singletons
1393	bool possible = false;
1394	if (rowLower[i] < -1.0e20) {
1395	double value = rowUpper[i];
1396	if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1397	possible = true;
1398	} else if (rowUpper[i] > 1.0e20) {
1399	double value = rowLower[i];
1400	if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1401	possible = true;
1402	} else {
1403	double value = rowUpper[i];
1404	if (rowLower[i] == rowUpper[i] &&
1405	fabs(value - floor(value + 0.5)) < 1.0e-8)
1406	possible = true;
1407	}
1408	if (possible) {
1409	for (CoinBigIndex j = rowStart[i];
1410	j < rowStart[i] + rowLength[i]; j++) {
1411	int iColumn = column[j];
1412	if (columnLength[iColumn] != 1 \|\| fabs(element[j]) != 1.0) {
1413	possible = false;
1414	break;
1415	}
1416	}
1417	if (possible) {
1418	for (CoinBigIndex j = rowStart[i];
1419	j < rowStart[i] + rowLength[i]; j++) {
1420	int iColumn = column[j];
1421	if (!copy2->isInteger(iColumn)) {
1422	copy2->setInteger(iColumn);
1423	int kOrig = original[iColumn];
1424	setOptionalInteger(kOrig);
1425	}
1426	}
1427	del[nDel++] = i;
1428	}
1429	}
1430	}
1431	}
1432	}
1433	if (nDel) {
1434	copy2->deleteRows(nDel, del);
1435	// pack down possible
1436	int n=0;
1437	for (int i=0;i<nDel;i++)
1438	possibleRow[del[i]]=-1;
1439	for (int i=0;i<numberRows;i++) {
1440	if (possibleRow[i]>=0)
1441	possibleRow[n++]=possibleRow[i];
1442	}
1443	}
1444	if (nDel != numberRows) {
1445	nDel = 0;
1446	for (int i = 0; i < numberColumns; i++) {
1447	if (copy2->isInteger(i)) {
1448	del[nDel++] = i;
1449	nExtra++;
1450	}
1451	}
1452	} else {
1453	nDel = 0;
1454	}
1455	delete copy1;
1456	copy1 = copy2->clone();
1457	delete copy2;
1458	}
1459	// See if what's left is a network
1460	bool couldBeNetwork = false;
1461	if (copy1->getNumRows() && copy1->getNumCols()) {
1462	#ifdef COIN_HAS_CLP
1463	OsiClpSolverInterface * clpSolver
1464	= dynamic_cast<OsiClpSolverInterface *> (copy1);
1465	if (false && clpSolver) {
1466	numberRows = clpSolver->getNumRows();
1467	char * rotate = new char[numberRows];
1468	int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
1469	delete [] rotate;
1470	#if CBC_USEFUL_PRINTING>1
1471	printf("INTA network %d rows out of %d\n", n, numberRows);
1472	#endif
1473	if (CoinAbs(n) == numberRows) {
1474	couldBeNetwork = true;
1475	for (int i = 0; i < numberRows; i++) {
1476	if (!possibleRow[i]) {
1477	couldBeNetwork = false;
1478	#if CBC_USEFUL_PRINTING>1
1479	printf("but row %d is bad\n", i);
1480	#endif
1481	break;
1482	}
1483	}
1484	}
1485	} else
1486	#endif
1487	{
1488	numberColumns = copy1->getNumCols();
1489	numberRows = copy1->getNumRows();
1490	const double * rowLower = copy1->getRowLower();
1491	const double * rowUpper = copy1->getRowUpper();
1492	couldBeNetwork = true;
1493	for (int i = 0; i < numberRows; i++) {
1494	if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
1495	couldBeNetwork = false;
1496	break;
1497	}
1498	if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
1499	couldBeNetwork = false;
1500	break;
1501	}
1502	if (possibleRow[i]==0) {
1503	couldBeNetwork = false;
1504	break;
1505	}
1506	}
1507	if (couldBeNetwork) {
1508	const CoinPackedMatrix * matrixByCol = copy1->getMatrixByCol();
1509	const double * element = matrixByCol->getElements();
1510	//const int * row = matrixByCol->getIndices();
1511	const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
1512	const int * columnLength = matrixByCol->getVectorLengths();
1513	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1514	CoinBigIndex start = columnStart[iColumn];
1515	CoinBigIndex end = start + columnLength[iColumn];
1516	if (end > start + 2) {
1517	couldBeNetwork = false;
1518	break;
1519	}
1520	int type = 0;
1521	for (CoinBigIndex j = start; j < end; j++) {
1522	double value = element[j];
1523	if (fabs(value) != 1.0) {
1524	couldBeNetwork = false;
1525	break;
1526	} else if (value == 1.0) {
1527	if ((type&1) == 0)
1528	type \|= 1;
1529	else
1530	type = 7;
1531	} else if (value == -1.0) {
1532	if ((type&2) == 0)
1533	type \|= 2;
1534	else
1535	type = 7;
1536	}
1537	}
1538	if (type > 3) {
1539	couldBeNetwork = false;
1540	break;
1541	}
1542	}
1543	}
1544	}
1545	}
1546	if (couldBeNetwork) {
1547	for (int i = 0; i < numberColumns; i++)
1548	setOptionalInteger(original[i]);
1549	}
1550	if (nExtra \|\| couldBeNetwork) {
1551	numberColumns = copy1->getNumCols();
1552	numberRows = copy1->getNumRows();
1553	if (!numberColumns \|\| !numberRows) {
1554	int numberColumns = solver_->getNumCols();
1555	for (int i = 0; i < numberColumns; i++)
1556	assert(solver_->isInteger(i));
1557	}
1558	#if CBC_USEFUL_PRINTING>1
1559	if (couldBeNetwork \|\| nExtra)
1560	printf("INTA %d extra integers, %d left%s\n", nExtra,
1561	numberColumns,
1562	couldBeNetwork ? ", all network" : "");
1563	#endif
1564	findIntegers(true, 2);
1565	convertToDynamic();
1566	}
1567	#if CBC_USEFUL_PRINTING>1
1568	if (!couldBeNetwork && copy1->getNumCols() &&
1569	copy1->getNumRows()) {
1570	printf("INTA %d rows and %d columns remain\n",
1571	copy1->getNumRows(), copy1->getNumCols());
1572	if (copy1->getNumCols() < 200) {
1573	copy1->writeMps("moreint");
1574	printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
1575	copy1->getNumRows(), copy1->getNumCols());
1576	}
1577	}
1578	#endif
1579	delete copy1;
1580	delete [] del;
1581	delete [] original;
1582	delete [] possibleRow;
1583	// double check increment
1584	analyzeObjective();
1585	// If any changes - tell code
1586	if(numberOriginalIntegers<numberIntegers_)
1587	synchronizeModel();
1588	}
1589	/**
1590	\todo
1591	Normally, it looks like we enter here from command dispatch in the main
1592	routine, after calling the solver for an initial solution
1593	(CbcModel::initialSolve, which simply calls the solver's initialSolve
1594	routine.) The first thing we do is call resolve. Presumably there are
1595	circumstances where this is nontrivial? There's also a call from
1596	CbcModel::originalModel (tied up with integer presolve), which should be
1597	checked.
1598
1599	*/
1600
1601	#ifdef CONFLICT_CUTS
1602	#if PRINT_CONFLICT==1
1603	static int numberConflictCuts=0;
1604	static int lastNumberConflictCuts=0;
1605	static double lengthConflictCuts=0.0;
1606	#endif
1607	#endif
1608	/*
1609	The overall flow can be divided into three stages:
1610	* Prep: Check that the lp relaxation remains feasible at the root. If so,
1611	do all the setup for B&C.
1612	* Process the root node: Generate cuts, apply heuristics, and in general do
1613	the best we can to resolve the problem without B&C.
1614	* Do B&C search until we hit a limit or exhaust the search tree.
1615
1616	Keep in mind that in general there is no node in the search tree that
1617	corresponds to the active subproblem. The active subproblem is represented
1618	by the current state of the model, of the solver, and of the constraint
1619	system held by the solver.
1620	*/
1621	#ifdef COIN_HAS_CPX
1622	#include "OsiCpxSolverInterface.hpp"
1623	#include "cplex.h"
1624	#endif
1625	void CbcModel::branchAndBound(int doStatistics)
1626
1627	{
1628	if (!parentModel_) {
1629	/*
1630	Capture a time stamp before we start (unless set).
1631	*/
1632	if (!dblParam_[CbcStartSeconds]) {
1633	if (!useElapsedTime())
1634	dblParam_[CbcStartSeconds] = CoinCpuTime();
1635	else
1636	dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
1637	}
1638	}
1639	dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
1640	dblParam_[CbcSumChange] = 0.0;
1641	dblParam_[CbcLargestChange] = 0.0;
1642	intParam_[CbcNumberBranches] = 0;
1643	double lastBestPossibleObjective=-COIN_DBL_MAX;
1644	// when to check for restart
1645	int nextCheckRestart=50;
1646	// Force minimization !!!!
1647	bool flipObjective = (solver_->getObjSense()<0.0);
1648	if (flipObjective)
1649	flipModel();
1650	dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
1651	strongInfo_[0] = 0;
1652	strongInfo_[1] = 0;
1653	strongInfo_[2] = 0;
1654	strongInfo_[3] = 0;
1655	strongInfo_[4] = 0;
1656	strongInfo_[5] = 0;
1657	strongInfo_[6] = 0;
1658	numberStrongIterations_ = 0;
1659	currentNode_ = NULL;
1660	// See if should do cuts old way
1661	if (parallelMode() < 0) {
1662	specialOptions_ \|= 4096 + 8192;
1663	} else if (parallelMode() > 0) {
1664	specialOptions_ \|= 4096;
1665	}
1666	int saveMoreSpecialOptions = moreSpecialOptions_;
1667	if (dynamic_cast<CbcTreeLocal *> (tree_))
1668	specialOptions_ \|= 4096 + 8192;
1669	#ifdef COIN_HAS_CLP
1670	{
1671	OsiClpSolverInterface * clpSolver
1672	= dynamic_cast<OsiClpSolverInterface *> (solver_);
1673	if (clpSolver) {
1674	// pass in disaster handler
1675	CbcDisasterHandler handler(this);
1676	clpSolver->passInDisasterHandler(&handler);
1677	// Initialise solvers seed (unless users says not)
1678	if ((specialOptions_&4194304)==0)
1679	clpSolver->getModelPtr()->setRandomSeed(1234567);
1680	#ifdef JJF_ZERO
1681	// reduce factorization frequency
1682	int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1683	clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
1684	#endif
1685	}
1686	}
1687	#endif
1688	// original solver (only set if pre-processing)
1689	OsiSolverInterface * originalSolver = NULL;
1690	int numberOriginalObjects = numberObjects_;
1691	OsiObject ** originalObject = NULL;
1692	// Save whether there were any objects
1693	bool noObjects = (numberObjects_ == 0);
1694	// Set up strategies
1695	/*
1696	See if the user has supplied a strategy object and deal with it if present.
1697	The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1698	perform integer preprocessing, if requested.
1699
1700	We need to hang on to a pointer to solver_. setupOther will assign a
1701	preprocessed solver to model, but will instruct assignSolver not to trash the
1702	existing one.
1703	*/
1704	if (strategy_) {
1705	// May do preprocessing
1706	originalSolver = solver_;
1707	strategy_->setupOther(*this);
1708	if (strategy_->preProcessState()) {
1709	// pre-processing done
1710	if (strategy_->preProcessState() < 0) {
1711	// infeasible (or unbounded)
1712	status_ = 0 ;
1713	if (!solver_->isProvenDualInfeasible()) {
1714	handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1715	secondaryStatus_ = 1;
1716	} else {
1717	handler_->message(CBC_UNBOUNDED,
1718	messages_) << CoinMessageEol ;
1719	secondaryStatus_ = 7;
1720	}
1721	originalContinuousObjective_ = COIN_DBL_MAX;
1722	if (flipObjective)
1723	flipModel();
1724	return ;
1725	} else if (numberObjects_ && object_) {
1726	numberOriginalObjects = numberObjects_;
1727	// redo sequence
1728	numberIntegers_ = 0;
1729	int numberColumns = getNumCols();
1730	int nOrig = originalSolver->getNumCols();
1731	CglPreProcess * process = strategy_->process();
1732	assert (process);
1733	const int * originalColumns = process->originalColumns();
1734	// allow for cliques etc
1735	nOrig = CoinMax(nOrig, originalColumns[numberColumns-1] + 1);
1736	// try and redo debugger
1737	OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1738	if (debugger) {
1739	if (numberColumns<=debugger->numberColumns())
1740	debugger->redoSolution(numberColumns, originalColumns);
1741	else
1742	debugger=NULL; // no idea how to handle (SOS?)
1743	}
1744	// User-provided solution might have been best. Synchronise.
1745	if (bestSolution_) {
1746	// need to redo - in case no better found in BAB
1747	// just get integer part right
1748	for (int i = 0; i < numberColumns; i++) {
1749	int jColumn = originalColumns[i];
1750	bestSolution_[i] = bestSolution_[jColumn];
1751	}
1752	}
1753	originalObject = object_;
1754	// object number or -1
1755	int * temp = new int[nOrig];
1756	int iColumn;
1757	for (iColumn = 0; iColumn < nOrig; iColumn++)
1758	temp[iColumn] = -1;
1759	int iObject;
1760	int nNonInt = 0;
1761	for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1762	iColumn = originalObject[iObject]->columnNumber();
1763	if (iColumn < 0) {
1764	nNonInt++;
1765	} else {
1766	temp[iColumn] = iObject;
1767	}
1768	}
1769	int numberNewIntegers = 0;
1770	int numberOldIntegers = 0;
1771	int numberOldOther = 0;
1772	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1773	int jColumn = originalColumns[iColumn];
1774	if (temp[jColumn] >= 0) {
1775	int iObject = temp[jColumn];
1776	CbcSimpleInteger * obj =
1777	dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1778	if (obj)
1779	numberOldIntegers++;
1780	else
1781	numberOldOther++;
1782	} else if (isInteger(iColumn)) {
1783	numberNewIntegers++;
1784	}
1785	}
1786	/*
1787	Allocate an array to hold the indices of the integer variables.
1788	Make a large enough array for all objects
1789	*/
1790	numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1791	object_ = new OsiObject * [numberObjects_];
1792	delete [] integerVariable_;
1793	integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
1794	/*
1795	Walk the variables again, filling in the indices and creating objects for
1796	the integer variables. Initially, the objects hold the index and upper &
1797	lower bounds.
1798	*/
1799	numberIntegers_ = 0;
1800	int n = originalColumns[numberColumns-1] + 1;
1801	int * backward = new int[n];
1802	int i;
1803	for ( i = 0; i < n; i++)
1804	backward[i] = -1;
1805	for (i = 0; i < numberColumns; i++)
1806	backward[originalColumns[i]] = i;
1807	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1808	int jColumn = originalColumns[iColumn];
1809	if (temp[jColumn] >= 0) {
1810	int iObject = temp[jColumn];
1811	CbcSimpleInteger * obj =
1812	dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1813	if (obj) {
1814	object_[numberIntegers_] = originalObject[iObject]->clone();
1815	// redo ids etc
1816	//object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1817	object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1818	integerVariable_[numberIntegers_++] = iColumn;
1819	}
1820	} else if (isInteger(iColumn)) {
1821	object_[numberIntegers_] =
1822	new CbcSimpleInteger(this, iColumn);
1823	integerVariable_[numberIntegers_++] = iColumn;
1824	}
1825	}
1826	delete [] backward;
1827	numberObjects_ = numberIntegers_;
1828	// Now append other column stuff
1829	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1830	int jColumn = originalColumns[iColumn];
1831	if (temp[jColumn] >= 0) {
1832	int iObject = temp[jColumn];
1833	CbcSimpleInteger * obj =
1834	dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1835	if (!obj) {
1836	object_[numberObjects_] = originalObject[iObject]->clone();
1837	// redo ids etc
1838	CbcObject * obj =
1839	dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1840	assert (obj);
1841	obj->redoSequenceEtc(this, numberColumns, originalColumns);
1842	numberObjects_++;
1843	}
1844	}
1845	}
1846	// now append non column stuff
1847	for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1848	iColumn = originalObject[iObject]->columnNumber();
1849	if (iColumn < 0) {
1850	// already has column numbers changed
1851	object_[numberObjects_] = originalObject[iObject]->clone();
1852	#ifdef JJF_ZERO
1853	// redo ids etc
1854	CbcObject * obj =
1855	dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1856	assert (obj);
1857	obj->redoSequenceEtc(this, numberColumns, originalColumns);
1858	#endif
1859	numberObjects_++;
1860	}
1861	}
1862	delete [] temp;
1863	if (!numberObjects_)
1864	handler_->message(CBC_NOINT, messages_) << CoinMessageEol ;
1865	} else {
1866	int numberColumns = getNumCols();
1867	CglPreProcess * process = strategy_->process();
1868	assert (process);
1869	const int * originalColumns = process->originalColumns();
1870	// try and redo debugger
1871	OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1872	if (debugger)
1873	debugger->redoSolution(numberColumns, originalColumns);
1874	}
1875	} else {
1876	//no preprocessing
1877	originalSolver = NULL;
1878	}
1879	strategy_->setupCutGenerators(*this);
1880	strategy_->setupHeuristics(*this);
1881	// Set strategy print level to models
1882	strategy_->setupPrinting(*this, handler_->logLevel());
1883	}
1884	eventHappened_ = false;
1885	CbcEventHandler *eventHandler = getEventHandler() ;
1886	if (eventHandler)
1887	eventHandler->setModel(this);
1888	#define CLIQUE_ANALYSIS
1889	#ifdef CLIQUE_ANALYSIS
1890	// set up for probing
1891	// If we're doing clever stuff with cliques, additional info here.
1892	if (!parentModel_)
1893	probingInfo_ = new CglTreeProbingInfo(solver_);
1894	else
1895	probingInfo_ = NULL;
1896	#else
1897	probingInfo_ = NULL;
1898	#endif
1899
1900	// Try for dominated columns
1901	if ((specialOptions_&64) != 0) {
1902	CglDuplicateRow dupcuts(solver_);
1903	dupcuts.setMode(2);
1904	CglStored * storedCuts = dupcuts.outDuplicates(solver_);
1905	if (storedCuts) {
1906	COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
1907	addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1908	true, false, false, -200);
1909	}
1910	}
1911	if (!nodeCompare_)
1912	nodeCompare_ = new CbcCompareDefault();;
1913	// See if hot start wanted
1914	CbcCompareBase * saveCompare = NULL;
1915	// User supplied hotstart. Adapt for preprocessing.
1916	if (hotstartSolution_) {
1917	if (strategy_ && strategy_->preProcessState() > 0) {
1918	CglPreProcess * process = strategy_->process();
1919	assert (process);
1920	int n = solver_->getNumCols();
1921	const int * originalColumns = process->originalColumns();
1922	// columns should be in order ... but
1923	double * tempS = new double[n];
1924	for (int i = 0; i < n; i++) {
1925	int iColumn = originalColumns[i];
1926	tempS[i] = hotstartSolution_[iColumn];
1927	}
1928	delete [] hotstartSolution_;
1929	hotstartSolution_ = tempS;
1930	if (hotstartPriorities_) {
1931	int * tempP = new int [n];
1932	for (int i = 0; i < n; i++) {
1933	int iColumn = originalColumns[i];
1934	tempP[i] = hotstartPriorities_[iColumn];
1935	}
1936	delete [] hotstartPriorities_;
1937	hotstartPriorities_ = tempP;
1938	}
1939	}
1940	saveCompare = nodeCompare_;
1941	// depth first
1942	nodeCompare_ = new CbcCompareDepth();
1943	}
1944	if (!problemFeasibility_)
1945	problemFeasibility_ = new CbcFeasibilityBase();
1946	# ifdef CBC_DEBUG
1947	std::string problemName ;
1948	solver_->getStrParam(OsiProbName, problemName) ;
1949	printf("Problem name - %s\n", problemName.c_str()) ;
1950	solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
1951	# endif
1952	/*
1953	Assume we're done, and see if we're proven wrong.
1954	*/
1955	status_ = 0 ;
1956	secondaryStatus_ = 0;
1957	phase_ = 0;
1958	/*
1959	Scan the variables, noting the integer variables. Create an
1960	CbcSimpleInteger object for each integer variable.
1961	*/
1962	findIntegers(false) ;
1963	// Say not dynamic pseudo costs
1964	ownership_ &= ~0x40000000;
1965	// If dynamic pseudo costs then do
1966	if (numberBeforeTrust_)
1967	convertToDynamic();
1968	// Set up char array to say if integer (speed)
1969	delete [] integerInfo_;
1970	{
1971	int n = solver_->getNumCols();
1972	integerInfo_ = new char [n];
1973	for (int i = 0; i < n; i++) {
1974	if (solver_->isInteger(i))
1975	integerInfo_[i] = 1;
1976	else
1977	integerInfo_[i] = 0;
1978	}
1979	}
1980	if (preferredWay_) {
1981	// set all unset ones
1982	for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1983	CbcObject * obj =
1984	dynamic_cast <CbcObject *>(object_[iObject]) ;
1985	if (obj && !obj->preferredWay())
1986	obj->setPreferredWay(preferredWay_);
1987	}
1988	}
1989	/*
1990	Ensure that objects on the lists of OsiObjects, heuristics, and cut
1991	generators attached to this model all refer to this model.
1992	*/
1993	synchronizeModel() ;
1994	if (!solverCharacteristics_) {
1995	OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1996	if (solverCharacteristics) {
1997	solverCharacteristics_ = solverCharacteristics;
1998	} else {
1999	// replace in solver
2000	OsiBabSolver defaultC;
2001	solver_->setAuxiliaryInfo(&defaultC);
2002	solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
2003	}
2004	}
2005
2006	solverCharacteristics_->setSolver(solver_);
2007	// Set so we can tell we are in initial phase in resolve
2008	continuousObjective_ = -COIN_DBL_MAX ;
2009	/*
2010	Solve the relaxation.
2011
2012	Apparently there are circumstances where this will be non-trivial --- i.e.,
2013	we've done something since initialSolve that's trashed the solution to the
2014	continuous relaxation.
2015	*/
2016	/* Tell solver we are in Branch and Cut
2017	Could use last parameter for subtle differences */
2018	solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
2019	#ifdef COIN_HAS_CLP
2020	{
2021	OsiClpSolverInterface * clpSolver
2022	= dynamic_cast<OsiClpSolverInterface *> (solver_);
2023	if (clpSolver) {
2024	ClpSimplex * clpSimplex = clpSolver->getModelPtr();
2025	if ((specialOptions_&32) == 0) {
2026	// take off names (unless going to be saving)
2027	if (numberAnalyzeIterations_>=0\|\|(-numberAnalyzeIterations_&64)==0)
2028	clpSimplex->dropNames();
2029	}
2030	// no crunch if mostly continuous
2031	if ((clpSolver->specialOptions()&(1 + 8)) != (1 + 8)) {
2032	int numberColumns = solver_->getNumCols();
2033	if (numberColumns > 1000 && numberIntegers_*4 < numberColumns)
2034	clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
2035	}
2036	//#define NO_CRUNCH
2037	#ifdef NO_CRUNCH
2038	printf("TEMP switching off crunch\n");
2039	int iOpt = clpSolver->specialOptions();
2040	iOpt &= ~1;
2041	iOpt \|= 65536;
2042	clpSolver->setSpecialOptions(iOpt);
2043	#endif
2044	}
2045	}
2046	#endif
2047	bool feasible;
2048	numberSolves_ = 0 ;
2049	{
2050	// check
2051	int numberOdd = 0;
2052	for (int i = 0; i < numberObjects_; i++) {
2053	CbcSimpleInteger * obj =
2054	dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
2055	if (!obj)
2056	numberOdd++;
2057	}
2058	if (numberOdd)
2059	moreSpecialOptions_ \|= 1073741824;
2060	}
2061	// If NLP then we assume already solved outside branchAndbound
2062	if (!solverCharacteristics_->solverType() \|\| solverCharacteristics_->solverType() == 4) {
2063	feasible = resolve(NULL, 0) != 0 ;
2064	} else {
2065	// pick up given status
2066	feasible = (solver_->isProvenOptimal() &&
2067	!solver_->isDualObjectiveLimitReached()) ;
2068	}
2069	if (problemFeasibility_->feasible(this, 0) < 0) {
2070	feasible = false; // pretend infeasible
2071	}
2072	numberSavedSolutions_ = 0;
2073	int saveNumberStrong = numberStrong_;
2074	int saveNumberBeforeTrust = numberBeforeTrust_;
2075	/*
2076	If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2077	continue with processing the root node.
2078	*/
2079	if (!feasible) {
2080	status_ = 0 ;
2081	if (!solver_->isProvenDualInfeasible()) {
2082	handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
2083	secondaryStatus_ = 1;
2084	} else {
2085	handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
2086	secondaryStatus_ = 7;
2087	}
2088	originalContinuousObjective_ = COIN_DBL_MAX;
2089	if (bestSolution_ &&
2090	((specialOptions_&8388608)==0\|\|(specialOptions_&2048)!=0)) {
2091	// best solution found by various heuristics - set solution
2092	char general[200];
2093	sprintf(general,"Solution of %g already found by heuristic",
2094	bestObjective_);
2095	messageHandler()->message(CBC_GENERAL,
2096	messages())
2097	<< general << CoinMessageEol ;
2098	setCutoff(1.0e50) ; // As best solution should be worse than cutoff
2099	// change cutoff as constraint if wanted
2100	if (cutoffRowNumber_>=0) {
2101	if (solver_->getNumRows()>cutoffRowNumber_)
2102	solver_->setRowUpper(cutoffRowNumber_,1.0e50);
2103	}
2104	// also in continuousSolver_
2105	if (continuousSolver_) {
2106	// Solvers know about direction
2107	double direction = solver_->getObjSense();
2108	continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
2109	} else {
2110	continuousSolver_ = solver_->clone();
2111	}
2112	phase_ = 5;
2113	double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2114	if ((specialOptions_&4) == 0)
2115	bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
2116	setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
2117	continuousSolver_->resolve() ;
2118	if (!continuousSolver_->isProvenOptimal()) {
2119	continuousSolver_->messageHandler()->setLogLevel(1) ;
2120	continuousSolver_->initialSolve() ;
2121	}
2122	delete solver_ ;
2123	solverCharacteristics_=NULL;
2124	solver_ = continuousSolver_ ;
2125	setPointers(solver_);
2126	continuousSolver_ = NULL ;
2127	}
2128	solverCharacteristics_ = NULL;
2129	if (flipObjective)
2130	flipModel();
2131	return ;
2132	} else if (!numberObjects_) {
2133	// nothing to do
2134	// Undo preprocessing performed during BaB.
2135	if (strategy_ && strategy_->preProcessState() > 0) {
2136	// undo preprocessing
2137	CglPreProcess * process = strategy_->process();
2138	assert (process);
2139	int n = originalSolver->getNumCols();
2140	if (bestSolution_) {
2141	delete [] bestSolution_;
2142	bestSolution_ = new double [n];
2143	process->postProcess(*solver_);
2144	}
2145	strategy_->deletePreProcess();
2146	// Solution now back in originalSolver
2147	delete solver_;
2148	solver_ = originalSolver;
2149	if (bestSolution_) {
2150	bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2151	memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
2152	}
2153	// put back original objects if there were any
2154	if (originalObject) {
2155	int iColumn;
2156	assert (ownObjects_);
2157	for (iColumn = 0; iColumn < numberObjects_; iColumn++)
2158	delete object_[iColumn];
2159	delete [] object_;
2160	numberObjects_ = numberOriginalObjects;
2161	object_ = originalObject;
2162	delete [] integerVariable_;
2163	numberIntegers_ = 0;
2164	for (iColumn = 0; iColumn < n; iColumn++) {
2165	if (solver_->isInteger(iColumn))
2166	numberIntegers_++;
2167	}
2168	integerVariable_ = new int[numberIntegers_];
2169	numberIntegers_ = 0;
2170	for (iColumn = 0; iColumn < n; iColumn++) {
2171	if (solver_->isInteger(iColumn))
2172	integerVariable_[numberIntegers_++] = iColumn;
2173	}
2174	}
2175	}
2176	if (flipObjective)
2177	flipModel();
2178	solverCharacteristics_ = NULL;
2179	bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2180	int numberColumns = solver_->getNumCols();
2181	delete [] bestSolution_;
2182	bestSolution_ = new double[numberColumns];
2183	CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2184	return ;
2185	}
2186	/*
2187	See if we're using the Osi side of the branching hierarchy. If so, either
2188	convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2189	first case, CbcModel owns the objects on the object_ list. In the second
2190	case, the solver holds the objects and object_ simply points to the
2191	solver's list.
2192
2193	080417 The conversion code here (the block protected by `if (obj)') cannot
2194	possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2195	all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2196	-> all other Cbc object classes. It's structurally impossible for any Osi
2197	object to descend from CbcObject. The only thing I can see is that this is
2198	really dead code, and object detection is now handled from the Osi side.
2199	*/
2200	// Convert to Osi if wanted
2201	//OsiBranchingInformation * persistentInfo = NULL;
2202	if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2203	//persistentInfo = new OsiBranchingInformation(solver_);
2204	if (numberOriginalObjects) {
2205	for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2206	CbcObject * obj =
2207	dynamic_cast <CbcObject *>(object_[iObject]) ;
2208	if (obj) {
2209	CbcSimpleInteger * obj2 =
2210	dynamic_cast <CbcSimpleInteger *>(obj) ;
2211	if (obj2) {
2212	// back to Osi land
2213	object_[iObject] = obj2->osiObject();
2214	delete obj;
2215	} else {
2216	OsiSimpleInteger * obj3 =
2217	dynamic_cast <OsiSimpleInteger *>(obj) ;
2218	if (!obj3) {
2219	OsiSOS * obj4 =
2220	dynamic_cast <OsiSOS *>(obj) ;
2221	if (!obj4) {
2222	CbcSOS * obj5 =
2223	dynamic_cast <CbcSOS *>(obj) ;
2224	if (obj5) {
2225	// back to Osi land
2226	object_[iObject] = obj5->osiObject(solver_);
2227	} else {
2228	printf("Code up CbcObject type in Osi land\n");
2229	abort();
2230	}
2231	}
2232	}
2233	}
2234	}
2235	}
2236	// and add to solver
2237	//if (!solver_->numberObjects()) {
2238	solver_->addObjects(numberObjects_, object_);
2239	//} else {
2240	//if (solver_->numberObjects()!=numberOriginalObjects) {
2241	//printf("should have trapped that solver has objects before\n");
2242	//abort();
2243	//}
2244	//}
2245	} else {
2246	/*
2247	As of 080104, findIntegersAndSOS is misleading --- the default OSI
2248	implementation finds only integers.
2249	*/
2250	// do from solver
2251	deleteObjects(false);
2252	solver_->findIntegersAndSOS(false);
2253	numberObjects_ = solver_->numberObjects();
2254	object_ = solver_->objects();
2255	ownObjects_ = false;
2256	}
2257	branchingMethod_->chooseMethod()->setSolver(solver_);
2258	}
2259	// take off heuristics if have to (some do not work with SOS, for example)
2260	// object should know what's safe.
2261	{
2262	int numberOdd = 0;
2263	int numberSOS = 0;
2264	for (int i = 0; i < numberObjects_; i++) {
2265	if (!object_[i]->canDoHeuristics())
2266	numberOdd++;
2267	CbcSOS * obj =
2268	dynamic_cast <CbcSOS *>(object_[i]) ;
2269	if (obj)
2270	numberSOS++;
2271	}
2272	if (numberOdd) {
2273	if (numberHeuristics_ && (specialOptions_&1024)==0 ) {
2274	int k = 0;
2275	for (int i = 0; i < numberHeuristics_; i++) {
2276	if (!heuristic_[i]->canDealWithOdd())
2277	delete heuristic_[i];
2278	else
2279	heuristic_[k++] = heuristic_[i];
2280	}
2281	if (!k) {
2282	delete [] heuristic_;
2283	heuristic_ = NULL;
2284	}
2285	numberHeuristics_ = k;
2286	handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ;
2287	}
2288	// If odd switch off AddIntegers
2289	specialOptions_ &= ~65536;
2290	// switch off fast nodes for now
2291	fastNodeDepth_ = -1;
2292	moreSpecialOptions_ &= ~33554432; // no diving
2293	} else if (numberSOS) {
2294	specialOptions_ \|= 128; // say can do SOS in dynamic mode
2295	// switch off fast nodes for now
2296	fastNodeDepth_ = -1;
2297	moreSpecialOptions_ &= ~33554432; // no diving
2298	}
2299	if (numberThreads_ > 0) {
2300	// switch off fast nodes for now
2301	//fastNodeDepth_ = -1;
2302	}
2303	}
2304	// Save objective (just so user can access it)
2305	originalContinuousObjective_ = solver_->getObjValue()* solver_->getObjSense();
2306	bestPossibleObjective_ = originalContinuousObjective_;
2307	sumChangeObjective1_ = 0.0;
2308	sumChangeObjective2_ = 0.0;
2309	/*
2310	OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2311	Assuming it recognises the problem, when called upon it will check a cut to
2312	see if it cuts off the optimal answer.
2313	*/
2314	// If debugger exists set specialOptions_ bit
2315	if (solver_->getRowCutDebuggerAlways()) {
2316	specialOptions_ \|= 1;
2317	}
2318
2319	# ifdef CBC_DEBUG
2320	if ((specialOptions_&1) == 0)
2321	solver_->activateRowCutDebugger(problemName.c_str()) ;
2322	if (solver_->getRowCutDebuggerAlways())
2323	specialOptions_ \|= 1;
2324	# endif
2325
2326	/*
2327	Begin setup to process a feasible root node.
2328	*/
2329	bestObjective_ = CoinMin(bestObjective_, 1.0e50) ;
2330	if (!bestSolution_) {
2331	numberSolutions_ = 0 ;
2332	numberHeuristicSolutions_ = 0 ;
2333	}
2334	stateOfSearch_ = 0;
2335	// Everything is minimization
2336	{
2337	// needed to sync cutoffs
2338	double value ;
2339	solver_->getDblParam(OsiDualObjectiveLimit, value) ;
2340	dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2341	}
2342	double cutoff = getCutoff() ;
2343	double direction = solver_->getObjSense() ;
2344	dblParam_[CbcOptimizationDirection] = direction;
2345	if (cutoff < 1.0e20 && direction < 0.0)
2346	messageHandler()->message(CBC_CUTOFF_WARNING1,
2347	messages())
2348	<< cutoff << -cutoff << CoinMessageEol ;
2349	if (cutoff > bestObjective_)
2350	cutoff = bestObjective_ ;
2351	setCutoff(cutoff) ;
2352	/*
2353	We probably already have a current solution, but just in case ...
2354	*/
2355	int numberColumns = getNumCols() ;
2356	if (!currentSolution_)
2357	currentSolution_ = new double[numberColumns] ;
2358	testSolution_ = currentSolution_;
2359	/*
2360	Create a copy of the solver, thus capturing the original (root node)
2361	constraint system (aka the continuous system).
2362	*/
2363	delete continuousSolver_;
2364	continuousSolver_ = solver_->clone() ;
2365	#ifdef COIN_HAS_NTY
2366	// maybe allow on fix and restart later
2367	if ((moreSpecialOptions2_&(128\|256))!=0&&!parentModel_) {
2368	symmetryInfo_ = new CbcSymmetry();
2369	symmetryInfo_->setupSymmetry(*continuousSolver_);
2370	int numberGenerators = symmetryInfo_->statsOrbits(this,0);
2371	if (!symmetryInfo_->numberUsefulOrbits()&&(moreSpecialOptions2_&(128\|256))!=(128\|256)) {
2372	delete symmetryInfo_;
2373	symmetryInfo_=NULL;
2374	moreSpecialOptions2_ &= ~(128\|256);
2375	}
2376	if ((moreSpecialOptions2_&(128\|256))==(128\|256)) {
2377	//moreSpecialOptions2_ &= ~256;
2378	}
2379	}
2380	#endif
2381
2382	// add cutoff as constraint if wanted
2383	if (cutoffRowNumber_==-2) {
2384	if (!parentModel_) {
2385	int numberColumns=solver_->getNumCols();
2386	double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
2387	int * indices = new int [numberColumns];
2388	int n=0;
2389	for (int i=0;i<numberColumns;i++) {
2390	if (obj[i]) {
2391	indices[n]=i;
2392	obj[n++]=obj[i];
2393	}
2394	}
2395	if (n) {
2396	double cutoff=getCutoff();
2397	// relax a little bit
2398	cutoff += 1.0e-4;
2399	double offset;
2400	solver_->getDblParam(OsiObjOffset, offset);
2401	cutoffRowNumber_ = solver_->getNumRows();
2402	solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
2403	} else {
2404	// no objective!
2405	cutoffRowNumber_ = -1;
2406	}
2407	delete [] indices;
2408	delete [] obj;
2409	} else {
2410	// switch off
2411	cutoffRowNumber_ = -1;
2412	}
2413	}
2414	numberRowsAtContinuous_ = getNumRows() ;
2415	solver_->saveBaseModel();
2416	/*
2417	Check the objective to see if we can deduce a nontrivial increment. If
2418	it's better than the current value for CbcCutoffIncrement, it'll be
2419	installed.
2420	*/
2421	if (solverCharacteristics_->reducedCostsAccurate())
2422	analyzeObjective() ;
2423	{
2424	// may be able to change cutoff now
2425	double cutoff = getCutoff();
2426	double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2427	if (cutoff > bestObjective_ - increment) {
2428	cutoff = bestObjective_ - increment ;
2429	setCutoff(cutoff) ;
2430	}
2431	}
2432	#ifdef COIN_HAS_CLP
2433	// Possible save of pivot method
2434	ClpDualRowPivot * savePivotMethod = NULL;
2435	{
2436	// pass tolerance and increment to solver
2437	OsiClpSolverInterface * clpSolver
2438	= dynamic_cast<OsiClpSolverInterface *> (solver_);
2439	if (clpSolver)
2440	clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2441	#ifdef CLP_RESOLVE
2442	if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
2443	resolveClp(clpSolver,0);
2444	}
2445	#endif
2446	}
2447	#endif
2448	/*
2449	Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2450	the active subproblem. whichGenerator will be used to record the generator
2451	that produced a given cut.
2452	*/
2453	#define INITIAL_MAXIMUM_WHICH 1000
2454	maximumWhich_ = INITIAL_MAXIMUM_WHICH ;
2455	delete [] whichGenerator_;
2456	whichGenerator_ = new int[maximumWhich_] ;
2457	memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2458	maximumNumberCuts_ = 0 ;
2459	currentNumberCuts_ = 0 ;
2460	delete [] addedCuts_ ;
2461	addedCuts_ = NULL ;
2462	OsiObject ** saveObjects = NULL;
2463	maximumRows_ = numberRowsAtContinuous_;
2464	currentDepth_ = 0;
2465	workingBasis_.resize(maximumRows_, numberColumns);
2466	/*
2467	Set up an empty heap and associated data structures to hold the live set
2468	(problems which require further exploration).
2469	*/
2470	CbcCompareDefault * compareActual
2471	= dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2472	if (compareActual) {
2473	compareActual->setBestPossible(direction*solver_->getObjValue());
2474	compareActual->setCutoff(getCutoff());
2475	#ifdef JJF_ZERO
2476	if (false && !numberThreads_ && !parentModel_) {
2477	printf("CbcTreeArray ? threads ? parentArray\n");
2478	// Setup new style tree
2479	delete tree_;
2480	tree_ = new CbcTreeArray();
2481	}
2482	#endif
2483	}
2484	tree_->setComparison(*nodeCompare_) ;
2485	/*
2486	Used to record the path from a node to the root of the search tree, so that
2487	we can then traverse from the root to the node when restoring a subproblem.
2488	*/
2489	maximumDepth_ = 10 ;
2490	delete [] walkback_ ;
2491	walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2492	lastDepth_ = 0;
2493	delete [] lastNodeInfo_ ;
2494	lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2495	delete [] lastNumberCuts_ ;
2496	lastNumberCuts_ = new int [maximumDepth_] ;
2497	maximumCuts_ = 100;
2498	lastNumberCuts2_ = 0;
2499	delete [] lastCut_;
2500	lastCut_ = new const OsiRowCut * [maximumCuts_];
2501	/*
2502	Used to generate bound edits for CbcPartialNodeInfo.
2503	*/
2504	double * lowerBefore = new double [numberColumns] ;
2505	double * upperBefore = new double [numberColumns] ;
2506	/*
2507	Set up to run heuristics and generate cuts at the root node. The heavy
2508	lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2509
2510	To start, tell cut generators they can be a bit more aggressive at the
2511	root node.
2512
2513	QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2514	with cuts at root'. Is phase_ = 1 the correct indication when
2515	doHeurisiticsAtRoot is called to run heuristics outside of the main
2516	cut / heurisitc / reoptimise loop in solveWithCuts?
2517
2518	Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2519	lifting. It will iterate a generate/reoptimise loop (including reduced cost
2520	fixing) until no cuts are generated, the change in objective falls off, or
2521	the limit on the number of rounds of cut generation is exceeded.
2522
2523	At the end of all this, any cuts will be recorded in cuts and also
2524	installed in the solver's constraint system. We'll have reoptimised, and
2525	removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2526	adjusted accordingly).
2527
2528	Tell cut generators they can be a bit more aggressive at root node
2529
2530	TODO: Why don't we make a copy of the solution after solveWithCuts?
2531	TODO: If numberUnsatisfied == 0, don't we have a solution?
2532	*/
2533	phase_ = 1;
2534	int iCutGenerator;
2535	for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2536	// If parallel switch off global cuts
2537	if (numberThreads_) {
2538	generator_[iCutGenerator]->setGlobalCuts(false);
2539	generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2540	}
2541	CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2542	generator->setAggressiveness(generator->getAggressiveness() + 100);
2543	if (!generator->canDoGlobalCuts())
2544	generator->setGlobalCuts(false);
2545	}
2546	OsiCuts cuts ;
2547	int anyAction = -1 ;
2548	numberOldActiveCuts_ = 0 ;
2549	numberNewCuts_ = 0 ;
2550	// Array to mark solution
2551	delete [] usedInSolution_;
2552	usedInSolution_ = new int[numberColumns];
2553	CoinZeroN(usedInSolution_, numberColumns);
2554	/*
2555	For printing totals and for CbcNode (numberNodes_)
2556	*/
2557	numberIterations_ = 0 ;
2558	numberNodes_ = 0 ;
2559	numberNodes2_ = 0 ;
2560	maximumStatistics_ = 0;
2561	maximumDepthActual_ = 0;
2562	numberDJFixed_ = 0.0;
2563	if (!parentModel_) {
2564	if ((specialOptions_&262144) != 0) {
2565	// create empty stored cuts
2566	//storedRowCuts_ = new CglStored(solver_->getNumCols());
2567	} else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2568	// tighten and set best solution
2569	// A) tight bounds on integer variables
2570	/*
2571	storedRowCuts_ are coming in from outside, probably for nonlinear.
2572	John was unsure about origin.
2573	*/
2574	const double * lower = solver_->getColLower();
2575	const double * upper = solver_->getColUpper();
2576	const double * tightLower = storedRowCuts_->tightLower();
2577	const double * tightUpper = storedRowCuts_->tightUpper();
2578	int nTightened = 0;
2579	for (int i = 0; i < numberIntegers_; i++) {
2580	int iColumn = integerVariable_[i];
2581	if (tightLower[iColumn] > lower[iColumn]) {
2582	nTightened++;
2583	solver_->setColLower(iColumn, tightLower[iColumn]);
2584	}
2585	if (tightUpper[iColumn] < upper[iColumn]) {
2586	nTightened++;
2587	solver_->setColUpper(iColumn, tightUpper[iColumn]);
2588	}
2589	}
2590	if (nTightened)
2591	COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2592	if (storedRowCuts_->bestObjective() < bestObjective_) {
2593	// B) best solution
2594	double objValue = storedRowCuts_->bestObjective();
2595	setBestSolution(CBC_SOLUTION, objValue,
2596	storedRowCuts_->bestSolution()) ;
2597	// Do heuristics
2598	// Allow RINS
2599	for (int i = 0; i < numberHeuristics_; i++) {
2600	CbcHeuristicRINS * rins
2601	= dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2602	if (rins) {
2603	rins->setLastNode(-100);
2604	}
2605	}
2606	}
2607	}
2608	}
2609	#ifdef SWITCH_VARIABLES
2610	// see if any switching variables
2611	if (numberIntegers_<solver_->getNumCols())
2612	findSwitching();
2613	#endif
2614	/*
2615	Run heuristics at the root. This is the only opportunity to run FPump; it
2616	will be removed from the heuristics list by doHeuristicsAtRoot.
2617	*/
2618	// See if multiple runs wanted
2619	CbcModel ** rootModels=NULL;
2620	if (!parentModel_&&multipleRootTries_%100) {
2621	double rootTimeCpu=CoinCpuTime();
2622	double startTimeRoot=CoinGetTimeOfDay();
2623	int numberRootThreads=1;
2624	/* undocumented fine tuning
2625	aabbcc where cc is number of tries
2626	bb if nonzero is number of threads
2627	aa if nonzero just do heuristics
2628	*/
2629	int numberModels = multipleRootTries_%100;
2630	#ifdef CBC_THREAD
2631	numberRootThreads = (multipleRootTries_/100)%100;
2632	if (!numberRootThreads)
2633	numberRootThreads=numberModels;
2634	#endif
2635	int otherOptions = (multipleRootTries_/10000)%100;
2636	rootModels = new CbcModel * [numberModels];
2637	unsigned int newSeed = randomSeed_;
2638	if (newSeed==0) {
2639	double time = fabs(CoinGetTimeOfDay());
2640	while (time>=COIN_INT_MAX)
2641	time *= 0.5;
2642	newSeed = static_cast<unsigned int>(time);
2643	} else if (newSeed<0) {
2644	newSeed = 123456789;
2645	#ifdef COIN_HAS_CLP
2646	OsiClpSolverInterface * clpSolver
2647	= dynamic_cast<OsiClpSolverInterface *> (solver_);
2648	if (clpSolver) {
2649	newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2650	}
2651	#endif
2652	}
2653	CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
2654	for (int i=0;i<numberModels;i++) {
2655	rootModels[i]=new CbcModel(*this);
2656	rootModels[i]->setNumberThreads(0);
2657	rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2658	rootModels[i]->setRandomSeed(newSeed+10000000*i);
2659	rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
2660	rootModels[i]->setMultipleRootTries(0);
2661	// use seed
2662	rootModels[i]->setSpecialOptions(specialOptions_ \|(4194304\|8388608));
2663	rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ &
2664	(~(134217728\|4194304)));
2665	rootModels[i]->setMoreSpecialOptions2(moreSpecialOptions2_ &
2666	(~(128\|256)));
2667	rootModels[i]->solver_->setWarmStart(basis);
2668	#ifdef COIN_HAS_CLP
2669	OsiClpSolverInterface * clpSolver
2670	= dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
2671	#define NEW_RANDOM_BASIS
2672	#ifdef NEW_RANDOM_BASIS
2673	if (i==0)
2674	continue;
2675	#endif
2676	if (clpSolver) {
2677	ClpSimplex * simplex = clpSolver->getModelPtr();
2678	if (defaultHandler_)
2679	simplex->setDefaultMessageHandler();
2680	simplex->setRandomSeed(newSeed+20000000*i);
2681	simplex->allSlackBasis();
2682	int logLevel=simplex->logLevel();
2683	if (logLevel==1)
2684	simplex->setLogLevel(0);
2685	if (i!=0) {
2686	#ifdef NEW_RANDOM_BASIS
2687	int numberRows = simplex->numberRows();
2688	int throwOut=20;//2+numberRows/100;
2689	for (int iThrow=0;iThrow<throwOut;iThrow++) {
2690	double random = simplex->randomNumberGenerator()->randomDouble();
2691	int iStart=static_cast<int>(random*numberRows);
2692	for (int j=iStart;j<numberRows;j++) {
2693	if (simplex->getRowStatus(j)!=ClpSimplex::basic) {
2694	simplex->setRowStatus(j,ClpSimplex::basic);
2695	break;
2696	}
2697	}
2698	}
2699	clpSolver->setWarmStart(NULL);
2700	#else
2701	double random = simplex->randomNumberGenerator()->randomDouble();
2702	int bias = static_cast<int>(random*(numberIterations/4));
2703	simplex->setMaximumIterations(numberIterations/2+bias);
2704	simplex->primal();
2705	simplex->setMaximumIterations(COIN_INT_MAX);
2706	simplex->dual();
2707	#endif
2708	} else {
2709	#ifndef NEW_RANDOM_BASIS
2710	simplex->primal();
2711	#endif
2712	#endif
2713	}
2714	#ifdef NEW_RANDOM_BASIS
2715	simplex->setLogLevel(logLevel);
2716	clpSolver->setWarmStart(NULL);
2717	#endif
2718	}
2719	for (int j=0;j<numberHeuristics_;j++)
2720	rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed()+100000000*i);
2721	for (int j=0;j<numberCutGenerators_;j++)
2722	rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2723	}
2724	delete basis;
2725	#ifdef CBC_THREAD
2726	if (numberRootThreads==1) {
2727	#endif
2728	for (int iModel=0;iModel<numberModels;iModel++) {
2729	doRootCbcThread(rootModels[iModel]);
2730	// see if solved at root node
2731	if (rootModels[iModel]->getMaximumNodes()) {
2732	feasible=false;
2733	break;
2734	}
2735	}
2736	#ifdef CBC_THREAD
2737	} else {
2738	Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2739	for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
2740	bool finished=false;
2741	for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2742	pthread_create(&(threadId[iModel-kModel].thr), NULL,
2743	doRootCbcThread,
2744	rootModels[iModel]);
2745	}
2746	// wait
2747	for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2748	pthread_join(threadId[iModel-kModel].thr, NULL);
2749	}
2750	// see if solved at root node
2751	for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2752	if (rootModels[iModel]->getMaximumNodes())
2753	finished=true;
2754	}
2755	if (finished) {
2756	feasible=false;
2757	break;
2758	}
2759	}
2760	delete [] threadId;
2761	}
2762	#endif
2763	// sort solutions
2764	int * which = new int [numberModels];
2765	double * value = new double [numberModels];
2766	int numberSolutions=0;
2767	for (int iModel=0;iModel<numberModels;iModel++) {
2768	if (rootModels[iModel]->bestSolution()) {
2769	which[numberSolutions]=iModel;
2770	value[numberSolutions++]=
2771	-rootModels[iModel]->getMinimizationObjValue();
2772	}
2773	}
2774	char general[100];
2775	rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2776	if (numberRootThreads==1)
2777	sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2778	rootTimeCpu);
2779	else
2780	sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2781	rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2782	messageHandler()->message(CBC_GENERAL,
2783	messages())
2784	<< general << CoinMessageEol ;
2785	CoinSort_2(value,value+numberSolutions,which);
2786	// to get name
2787	CbcHeuristicRINS dummyHeuristic;
2788	dummyHeuristic.setHeuristicName("Multiple root solvers");
2789	lastHeuristic_=&dummyHeuristic;
2790	for (int i=0;i<numberSolutions;i++) {
2791	double objValue = - value[i];
2792	if (objValue<getCutoff()) {
2793	int iModel=which[i];
2794	setBestSolution(CBC_ROUNDING,objValue,
2795	rootModels[iModel]->bestSolution());
2796	}
2797	}
2798	lastHeuristic_=NULL;
2799	delete [] which;
2800	delete [] value;
2801	}
2802	// Do heuristics
2803	if (numberObjects_&&!rootModels)
2804	doHeuristicsAtRoot();
2805	if (solverCharacteristics_->solutionAddsCuts()) {
2806	// With some heuristics solver needs a resolve here
2807	solver_->resolve();
2808	if(!isProvenOptimal()){
2809	solver_->initialSolve();
2810	}
2811	}
2812	/*
2813	Grepping through the code, it would appear that this is a command line
2814	debugging hook. There's no obvious place in the code where this is set to
2815	a negative value.
2816
2817	User hook, says John.
2818	*/
2819	if ( intParam_[CbcMaxNumNode] < 0
2820	\|\|numberSolutions_>=getMaximumSolutions())
2821	eventHappened_ = true; // stop as fast as possible
2822	stoppedOnGap_ = false ;
2823	// See if can stop on gap
2824	bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2825	if(canStopOnGap()) {
2826	if (bestPossibleObjective_ < getCutoff())
2827	stoppedOnGap_ = true ;
2828	feasible = false;
2829	//eventHappened_=true; // stop as fast as possible
2830	}
2831	/*
2832	Set up for statistics collection, if requested. Standard values are
2833	documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2834	CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2835	command line debugging hook.
2836	*/
2837	statistics_ = NULL;
2838	// Do on switch
2839	if (doStatistics > 0 && doStatistics <= 100) {
2840	maximumStatistics_ = 10000;
2841	statistics_ = new CbcStatistics * [maximumStatistics_];
2842	memset(statistics_, 0, maximumStatistics_sizeof(CbcStatistics ));
2843	}
2844	// See if we can add integers
2845	// taken out in stable as causes problems
2846	if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_ && false)
2847	AddIntegers();
2848	/*
2849	Do an initial round of cut generation for the root node. Depending on the
2850	type of underlying solver, we may want to do this even if the initial query
2851	to the objects indicates they're satisfied.
2852
2853	solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2854	loop (including reduced cost fixing) until no cuts are generated, the
2855	change in objective falls off, or the limit on the number of rounds of cut
2856	generation is exceeded.
2857
2858	At the end of all this, any cuts will be recorded in cuts and also
2859	installed in the solver's constraint system. We'll have reoptimised, and
2860	removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2861	adjusted accordingly).
2862	*/
2863	int iObject ;
2864	int numberUnsatisfied = 0 ;
2865	delete [] currentSolution_;
2866	currentSolution_ = new double [numberColumns];
2867	testSolution_ = currentSolution_;
2868	memcpy(currentSolution_, solver_->getColSolution(),
2869	numberColumns*sizeof(double)) ;
2870	// point to useful information
2871	OsiBranchingInformation usefulInfo = usefulInformation();
2872
2873	for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2874	double infeasibility =
2875	object_[iObject]->checkInfeasibility(&usefulInfo) ;
2876	if (infeasibility ) numberUnsatisfied++ ;
2877	}
2878	// replace solverType
2879	double * tightBounds = NULL;
2880	if (solverCharacteristics_->tryCuts()) {
2881
2882	if (numberUnsatisfied) {
2883	// User event
2884	if (!eventHappened_ && feasible) {
2885	if (rootModels) {
2886	// for fixings
2887	int numberColumns=solver_->getNumCols();
2888	tightBounds = new double [2*numberColumns];
2889	{
2890	const double * lower = solver_->getColLower();
2891	const double * upper = solver_->getColUpper();
2892	for (int i=0;i<numberColumns;i++) {
2893	tightBounds[2*i+0]=lower[i];
2894	tightBounds[2*i+1]=upper[i];
2895	}
2896	}
2897	int numberModels = multipleRootTries_%100;
2898	const OsiSolverInterface ** solvers = new
2899	const OsiSolverInterface * [numberModels];
2900	int numberRows=continuousSolver_->getNumRows();
2901	int maxCuts=0;
2902	for (int i=0;i<numberModels;i++) {
2903	solvers[i]=rootModels[i]->solver();
2904	const double * lower = solvers[i]->getColLower();
2905	const double * upper = solvers[i]->getColUpper();
2906	for (int j=0;j<numberColumns;j++) {
2907	tightBounds[2j+0]=CoinMax(lower[j],tightBounds[2j+0]);
2908	tightBounds[2j+1]=CoinMin(upper[j],tightBounds[2j+1]);
2909	}
2910	int numberRows2=solvers[i]->getNumRows();
2911	assert (numberRows2>=numberRows);
2912	maxCuts += numberRows2-numberRows;
2913	// accumulate statistics
2914	for (int j=0;j<numberCutGenerators_;j++) {
2915	generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2916	}
2917	}
2918	for (int j=0;j<numberCutGenerators_;j++) {
2919	generator_[j]->scaleBackStatistics(numberModels);
2920	}
2921	//CbcRowCuts rowCut(maxCuts);
2922	const OsiRowCutDebugger *debugger = NULL;
2923	if ((specialOptions_&1) != 0)
2924	debugger = solver_->getRowCutDebugger() ;
2925	for (int iModel=0;iModel<numberModels;iModel++) {
2926	int numberRows2=solvers[iModel]->getNumRows();
2927	const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
2928	const int * rowLength = rowCopy->getVectorLengths();
2929	const double * elements = rowCopy->getElements();
2930	const int * column = rowCopy->getIndices();
2931	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2932	const double * rowLower = solvers[iModel]->getRowLower();
2933	const double * rowUpper = solvers[iModel]->getRowUpper();
2934	for (int iRow=numberRows;iRow<numberRows2;iRow++) {
2935	OsiRowCut rc;
2936	rc.setLb(rowLower[iRow]);
2937	rc.setUb(rowUpper[iRow]);
2938	CoinBigIndex start = rowStart[iRow];
2939	rc.setRow(rowLength[iRow],column+start,elements+start,false);
2940	if (debugger)
2941	CoinAssert (!debugger->invalidCut(rc));
2942	globalCuts_.addCutIfNotDuplicate(rc);
2943	}
2944	//int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2945	//numberCuts += cutsAdded;
2946	//printf("Model %d gave %d cuts (out of %d possible)\n",
2947	// iModel,cutsAdded,numberRows2-numberRows);
2948	}
2949	// normally replace global cuts
2950	//if (!globalCuts_.())
2951	//globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2952	//rowCut.addCuts(globalCuts_);
2953	int nTightened=0;
2954	assert(feasible);
2955	{
2956	double tolerance=1.0e-5;
2957	const double * lower = solver_->getColLower();
2958	const double * upper = solver_->getColUpper();
2959	for (int i=0;i<numberColumns;i++) {
2960	if (tightBounds[2i+0]>tightBounds[2i+1]+1.0e-9) {
2961	feasible=false;
2962	char general[200];
2963	sprintf(general,"Solvers give infeasible bounds on %d %g,%g was %g,%g - search finished\n",
2964	i,tightBounds[2i+0],tightBounds[2i+1],lower[i],upper[i]);
2965	messageHandler()->message(CBC_GENERAL,messages())
2966	<< general << CoinMessageEol ;
2967	break;
2968	}
2969	double oldLower=lower[i];
2970	double oldUpper=upper[i];
2971	if (tightBounds[2*i+0]>oldLower+tolerance) {
2972	nTightened++;
2973	solver_->setColLower(i,tightBounds[2*i+0]);
2974	}
2975	if (tightBounds[2*i+1]<oldUpper-tolerance) {
2976	nTightened++;
2977	solver_->setColUpper(i,tightBounds[2*i+1]);
2978	}
2979	}
2980	}
2981	delete [] tightBounds;
2982	tightBounds=NULL;
2983	char printBuffer[200];
2984	sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
2985	numberModels,globalCuts_.sizeRowCuts(),maxCuts);
2986	messageHandler()->message(CBC_GENERAL, messages())
2987	<< printBuffer << CoinMessageEol ;
2988	if (nTightened) {
2989	sprintf(printBuffer,"%d bounds were tightened",
2990	nTightened);
2991	messageHandler()->message(CBC_GENERAL, messages())
2992	<< printBuffer << CoinMessageEol ;
2993	}
2994	delete [] solvers;
2995	}
2996	if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
2997	// load cuts from file
2998	FILE * fp = fopen("global.cuts","rb");
2999	if (fp) {
3000	size_t nRead;
3001	int numberColumns=solver_->getNumCols();
3002	int numCols;
3003	nRead = fread(&numCols, sizeof(int), 1, fp);
3004	if (nRead != 1)
3005	throw("Error in fread");
3006	if (numberColumns!=numCols) {
3007	printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3008	fclose(fp);
3009	} else {
3010	// If rootModel just do some
3011	double threshold=-1.0;
3012	if (!multipleRootTries_)
3013	threshold=0.5;
3014	int initialCuts=0;
3015	int initialGlobal = globalCuts_.sizeRowCuts();
3016	double * elements = new double [numberColumns+2];
3017	int * indices = new int [numberColumns];
3018	int numberEntries=1;
3019	while (numberEntries>0) {
3020	nRead = fread(&numberEntries, sizeof(int), 1, fp);
3021	if (nRead != 1)
3022	throw("Error in fread");
3023	double randomNumber=randomNumberGenerator_.randomDouble();
3024	if (numberEntries>0) {
3025	initialCuts++;
3026	nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3027	if (nRead != static_cast<size_t>(numberEntries+2))
3028	throw("Error in fread");
3029	nRead = fread(indices, sizeof(int), numberEntries, fp);
3030	if (nRead != static_cast<size_t>(numberEntries))
3031	throw("Error in fread");
3032	if (randomNumber>threshold) {
3033	OsiRowCut rc;
3034	rc.setLb(elements[numberEntries]);
3035	rc.setUb(elements[numberEntries+1]);
3036	rc.setRow(numberEntries,indices,elements,
3037	false);
3038	rc.setGloballyValidAsInteger(2);
3039	globalCuts_.addCutIfNotDuplicate(rc) ;
3040	}
3041	}
3042	}
3043	fclose(fp);
3044	// fixes
3045	int nTightened=0;
3046	fp = fopen("global.fix","rb");
3047	if (fp) {
3048	nRead = fread(indices, sizeof(int), 2, fp);
3049	if (nRead != 2)
3050	throw("Error in fread");
3051	if (numberColumns!=indices[0]) {
3052	printf("Mismatch on columns %d %d\n",numberColumns,
3053	indices[0]);
3054	} else {
3055	indices[0]=1;
3056	while (indices[0]>=0) {
3057	nRead = fread(indices, sizeof(int), 2, fp);
3058	if (nRead != 2)
3059	throw("Error in fread");
3060	int iColumn=indices[0];
3061	if (iColumn>=0) {
3062	nTightened++;
3063	nRead = fread(elements, sizeof(double), 4, fp);
3064	if (nRead != 4)
3065	throw("Error in fread");
3066	solver_->setColLower(iColumn,elements[0]);
3067	solver_->setColUpper(iColumn,elements[1]);
3068	}
3069	}
3070	}
3071	}
3072	if (fp)
3073	fclose(fp);
3074	char printBuffer[200];
3075	sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
3076	initialCuts,
3077	globalCuts_.sizeRowCuts()-initialGlobal,nTightened);
3078	messageHandler()->message(CBC_GENERAL, messages())
3079	<< printBuffer << CoinMessageEol ;
3080	delete [] elements;
3081	delete [] indices;
3082	}
3083	}
3084	}
3085	if (feasible)
3086	feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3087	NULL);
3088	if (multipleRootTries_&&
3089	(moreSpecialOptions_&134217728)!=0) {
3090	FILE * fp=NULL;
3091	size_t nRead;
3092	int numberColumns=solver_->getNumCols();
3093	int initialCuts=0;
3094	if ((moreSpecialOptions_&134217728)!=0) {
3095	// append so go down to end
3096	fp = fopen("global.cuts","r+b");
3097	if (fp) {
3098	int numCols;
3099	nRead = fread(&numCols, sizeof(int), 1, fp);
3100	if (nRead != 1)
3101	throw("Error in fread");
3102	if (numberColumns!=numCols) {
3103	printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3104	fclose(fp);
3105	fp=NULL;
3106	}
3107	}
3108	}
3109	double * elements = new double [numberColumns+2];
3110	int * indices = new int [numberColumns];
3111	if (fp) {
3112	int numberEntries=1;
3113	while (numberEntries>0) {
3114	fpos_t position;
3115	fgetpos(fp, &position);
3116	nRead = fread(&numberEntries, sizeof(int), 1, fp);
3117	if (nRead != 1)
3118	throw("Error in fread");
3119	if (numberEntries>0) {
3120	initialCuts++;
3121	nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3122	if (nRead != static_cast<size_t>(numberEntries+2))
3123	throw("Error in fread");
3124	nRead = fread(indices, sizeof(int), numberEntries, fp);
3125	if (nRead != static_cast<size_t>(numberEntries))
3126	throw("Error in fread");
3127	} else {
3128	// end
3129	fsetpos(fp, &position);
3130	}
3131	}
3132	} else {
3133	fp = fopen("global.cuts","wb");
3134	size_t nWrite;
3135	nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
3136	if (nWrite != 1)
3137	throw("Error in fwrite");
3138	}
3139	size_t nWrite;
3140	// now append binding cuts
3141	int numberC=continuousSolver_->getNumRows();
3142	int numberRows=solver_->getNumRows();
3143	printf("Saving %d cuts (up from %d)\n",
3144	initialCuts+numberRows-numberC,initialCuts);
3145	const double * rowLower = solver_->getRowLower();
3146	const double * rowUpper = solver_->getRowUpper();
3147	// Row copy
3148	CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3149	const double * elementByRow = matrixByRow.getElements();
3150	const int * column = matrixByRow.getIndices();
3151	const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3152	const int * rowLength = matrixByRow.getVectorLengths();
3153	for (int iRow=numberC;iRow<numberRows;iRow++) {
3154	int n=rowLength[iRow];
3155	assert (n);
3156	CoinBigIndex start=rowStart[iRow];
3157	memcpy(elements,elementByRow+start,n*sizeof(double));
3158	memcpy(indices,column+start,n*sizeof(int));
3159	elements[n]=rowLower[iRow];
3160	elements[n+1]=rowUpper[iRow];
3161	nWrite=fwrite(&n,sizeof(int),1,fp);
3162	if (nWrite != 1)
3163	throw("Error in fwrite");
3164	nWrite=fwrite(elements,sizeof(double),n+2,fp);
3165	if (nWrite != static_cast<size_t>(n+2))
3166	throw("Error in fwrite");
3167	nWrite=fwrite(indices,sizeof(int),n,fp);
3168	if (nWrite != static_cast<size_t>(n))
3169	throw("Error in fwrite");
3170	}
3171	// eof marker
3172	int eofMarker=-1;
3173	nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
3174	if (nWrite != 1)
3175	throw("Error in fwrite");
3176	fclose(fp);
3177	// do tighter bounds (? later extra to original columns)
3178	int nTightened=0;
3179	const double * lower = solver_->getColLower();
3180	const double * upper = solver_->getColUpper();
3181	const double * originalLower = continuousSolver_->getColLower();
3182	const double * originalUpper = continuousSolver_->getColUpper();
3183	double tolerance=1.0e-5;
3184	for (int i=0;i<numberColumns;i++) {
3185	if (lower[i]>originalLower[i]+tolerance) {
3186	nTightened++;
3187	}
3188	if (upper[i]<originalUpper[i]-tolerance) {
3189	nTightened++;
3190	}
3191	}
3192	if (nTightened) {
3193	fp = fopen("global.fix","wb");
3194	size_t nWrite;
3195	indices[0]=numberColumns;
3196	if (originalColumns_)
3197	indices[1]=COIN_INT_MAX;
3198	else
3199	indices[1]=-1;
3200	nWrite=fwrite(indices,sizeof(int),2,fp);
3201	if (nWrite != 2)
3202	throw("Error in fwrite");
3203	for (int i=0;i<numberColumns;i++) {
3204	int nTightened=0;
3205	if (lower[i]>originalLower[i]+tolerance) {
3206	nTightened++;
3207	}
3208	if (upper[i]<originalUpper[i]-tolerance) {
3209	nTightened++;
3210	}
3211	if (nTightened) {
3212	indices[0]=i;
3213	if (originalColumns_)
3214	indices[1]=originalColumns_[i];
3215	elements[0]=lower[i];
3216	elements[1]=upper[i];
3217	elements[2]=originalLower[i];
3218	elements[3]=originalUpper[i];
3219	nWrite=fwrite(indices,sizeof(int),2,fp);
3220	if (nWrite != 2)
3221	throw("Error in fwrite");
3222	nWrite=fwrite(elements,sizeof(double),4,fp);
3223	if (nWrite != 4)
3224	throw("Error in fwrite");
3225	}
3226	}
3227	// eof marker
3228	indices[0]=-1;
3229	nWrite=fwrite(indices,sizeof(int),2,fp);
3230	if (nWrite != 2)
3231	throw("Error in fwrite");
3232	fclose(fp);
3233	}
3234	delete [] elements;
3235	delete [] indices;
3236	}
3237	if ((specialOptions_&524288) != 0 && !parentModel_
3238	&& storedRowCuts_) {
3239	if (feasible) {
3240	/* pick up stuff and try again
3241	add cuts, maybe keep around
3242	do best solution and if so new heuristics
3243	obviously tighten bounds
3244	*/
3245	// A and B probably done on entry
3246	// A) tight bounds on integer variables
3247	const double * lower = solver_->getColLower();
3248	const double * upper = solver_->getColUpper();
3249	const double * tightLower = storedRowCuts_->tightLower();
3250	const double * tightUpper = storedRowCuts_->tightUpper();
3251	int nTightened = 0;
3252	for (int i = 0; i < numberIntegers_; i++) {
3253	int iColumn = integerVariable_[i];
3254	if (tightLower[iColumn] > lower[iColumn]) {
3255	nTightened++;
3256	solver_->setColLower(iColumn, tightLower[iColumn]);
3257	}
3258	if (tightUpper[iColumn] < upper[iColumn]) {
3259	nTightened++;
3260	solver_->setColUpper(iColumn, tightUpper[iColumn]);
3261	}
3262	}
3263	if (nTightened)
3264	COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3265	if (storedRowCuts_->bestObjective() < bestObjective_) {
3266	// B) best solution
3267	double objValue = storedRowCuts_->bestObjective();
3268	setBestSolution(CBC_SOLUTION, objValue,
3269	storedRowCuts_->bestSolution()) ;
3270	// Do heuristics
3271	// Allow RINS
3272	for (int i = 0; i < numberHeuristics_; i++) {
3273	CbcHeuristicRINS * rins
3274	= dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3275	if (rins) {
3276	rins->setLastNode(-100);
3277	}
3278	}
3279	doHeuristicsAtRoot();
3280	}
3281	#ifdef JJF_ZERO
3282	int nCuts = storedRowCuts_->sizeRowCuts();
3283	// add to global list
3284	for (int i = 0; i < nCuts; i++) {
3285	OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3286	newCut.setGloballyValidAsInteger(2);
3287	newCut.mutableRow().setTestForDuplicateIndex(false);
3288	globalCuts_.insert(newCut) ;
3289	}
3290	#else
3291	addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3292	true, false, false, -200);
3293	#endif
3294	// Set cuts as active
3295	delete [] addedCuts_ ;
3296	maximumNumberCuts_ = cuts.sizeRowCuts();
3297	if (maximumNumberCuts_) {
3298	addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3299	} else {
3300	addedCuts_ = NULL;
3301	}
3302	for (int i = 0; i < maximumNumberCuts_; i++)
3303	addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3304	NULL, -1, -1, 2);
3305	COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3306	cuts = OsiCuts();
3307	currentNumberCuts_ = maximumNumberCuts_;
3308	feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3309	NULL);
3310	for (int i = 0; i < maximumNumberCuts_; i++)
3311	delete addedCuts_[i];
3312	}
3313	delete storedRowCuts_;
3314	storedRowCuts_ = NULL;
3315	}
3316	} else {
3317	feasible = false;
3318	}
3319	} else if (solverCharacteristics_->solutionAddsCuts() \|\|
3320	solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3321	// may generate cuts and turn the solution
3322	//to an infeasible one
3323	feasible = solveWithCuts(cuts, 2,
3324	NULL);
3325	}
3326	}
3327	if (rootModels) {
3328	int numberModels = multipleRootTries_%100;
3329	for (int i=0;i<numberModels;i++)
3330	delete rootModels[i];
3331	delete [] rootModels;
3332	}
3333	// check extra info on feasibility
3334	if (!solverCharacteristics_->mipFeasible())
3335	feasible = false;
3336	// If max nodes==0 - don't do strong branching
3337	if (!getMaximumNodes()) {
3338	if (feasible)
3339	feasible=false;
3340	else
3341	setMaximumNodes(1); //allow to stop on success
3342	}
3343	topOfTree_=NULL;
3344	#ifdef CLP_RESOLVE
3345	if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3346	OsiClpSolverInterface * clpSolver
3347	= dynamic_cast<OsiClpSolverInterface *> (solver_);
3348	if (clpSolver)
3349	resolveClp(clpSolver,0);
3350	}
3351	#endif
3352	// make cut generators less aggressive
3353	for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3354	CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3355	generator->setAggressiveness(generator->getAggressiveness() - 100);
3356	}
3357	currentNumberCuts_ = numberNewCuts_ ;
3358	if (solverCharacteristics_->solutionAddsCuts()) {
3359	// With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3360	solver_->resolve();
3361	if(!isProvenOptimal()){
3362	solver_->initialSolve();
3363	}
3364	}
3365	// See if can stop on gap
3366	bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3367	if(canStopOnGap()) {
3368	if (bestPossibleObjective_ < getCutoff())
3369	stoppedOnGap_ = true ;
3370	feasible = false;
3371	}
3372	// User event
3373	if (eventHappened_)
3374	feasible = false;
3375	#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
3376	/*
3377	This is the notion of using Cbc stuff to get going, then calling cplex to
3378	finish off.
3379	*/
3380	if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3381	// Use Cplex to do search!
3382	double time1 = CoinCpuTime();
3383	OsiClpSolverInterface * clpSolver
3384	= dynamic_cast<OsiClpSolverInterface *> (solver_);
3385	OsiCpxSolverInterface cpxSolver;
3386	double direction = clpSolver->getObjSense();
3387	cpxSolver.setObjSense(direction);
3388	// load up cplex
3389	const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3390	const double * rowLower = continuousSolver_->getRowLower();
3391	const double * rowUpper = continuousSolver_->getRowUpper();
3392	const double * columnLower = continuousSolver_->getColLower();
3393	const double * columnUpper = continuousSolver_->getColUpper();
3394	const double * objective = continuousSolver_->getObjCoefficients();
3395	cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3396	objective, rowLower, rowUpper);
3397	double * setSol = new double [numberIntegers_];
3398	int * setVar = new int [numberIntegers_];
3399	// cplex doesn't know about objective offset
3400	double offset = clpSolver->getModelPtr()->objectiveOffset();
3401	for (int i = 0; i < numberIntegers_; i++) {
3402	int iColumn = integerVariable_[i];
3403	cpxSolver.setInteger(iColumn);
3404	if (bestSolution_) {
3405	setSol[i] = bestSolution_[iColumn];
3406	setVar[i] = iColumn;
3407	}
3408	}
3409	CPXENVptr env = cpxSolver.getEnvironmentPtr();
3410	CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3411	cpxSolver.switchToMIP();
3412	if (bestSolution_) {
3413	#if 0
3414	CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3415	#else
3416	int zero = 0;
3417	CPXaddmipstarts(env, lpPtr, 1, numberIntegers_, &zero, setVar, setSol, NULL, NULL);
3418	#endif
3419	}
3420	if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3421	// add cuts
3422	const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3423	const double * rhs = clpSolver->getRightHandSide();
3424	const char * rowSense = clpSolver->getRowSense();
3425	const double * elementByRow = matrix->getElements();
3426	const int * column = matrix->getIndices();
3427	const CoinBigIndex * rowStart = matrix->getVectorStarts();
3428	const int * rowLength = matrix->getVectorLengths();
3429	int nStart = continuousSolver_->getNumRows();
3430	int nRows = clpSolver->getNumRows();
3431	int size = rowStart[nRows-1] + rowLength[nRows-1] -
3432	rowStart[nStart];
3433	int nAdd = 0;
3434	double * rmatval = new double [size];
3435	int * rmatind = new int [size];
3436	int * rmatbeg = new int [nRows-nStart+1];
3437	size = 0;
3438	rmatbeg[0] = 0;
3439	for (int i = nStart; i < nRows; i++) {
3440	for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3441	rmatind[size] = column[k];
3442	rmatval[size++] = elementByRow[k];
3443	}
3444	nAdd++;
3445	rmatbeg[nAdd] = size;
3446	}
3447	CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3448	rhs, rowSense, rmatbeg,
3449	rmatind, rmatval, NULL);
3450	CPXsetintparam( env, CPX_PARAM_REDUCE,
3451	// CPX_PREREDUCE_NOPRIMALORDUAL (0)
3452	CPX_PREREDUCE_PRIMALONLY);
3453	}
3454	if (getCutoff() < 1.0e50) {
3455	double useCutoff = getCutoff() + offset;
3456	if (bestObjective_ < 1.0e50)
3457	useCutoff = bestObjective_ + offset + 1.0e-7;
3458	cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3459	direction);
3460	if ( direction > 0.0 )
3461	CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3462	else
3463	CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3464	}
3465	CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3466	delete [] setSol;
3467	delete [] setVar;
3468	char printBuffer[200];
3469	if (offset) {
3470	sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3471	-offset);
3472	messageHandler()->message(CBC_GENERAL, messages())
3473	<< printBuffer << CoinMessageEol ;
3474	cpxSolver.setDblParam(OsiObjOffset, offset);
3475	}
3476	cpxSolver.branchAndBound();
3477	double timeTaken = CoinCpuTime() - time1;
3478	sprintf(printBuffer, "Cplex took %g seconds",
3479	timeTaken);
3480	messageHandler()->message(CBC_GENERAL, messages())
3481	<< printBuffer << CoinMessageEol ;
3482	numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3483	numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3484	double value = cpxSolver.getObjValue() * direction;
3485	if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3486	feasible = true;
3487	clpSolver->setWarmStart(NULL);
3488	// try and do solution
3489	double * newSolution =
3490	CoinCopyOfArray(cpxSolver.getColSolution(),
3491	getNumCols());
3492	setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3493	delete [] newSolution;
3494	}
3495	feasible = false;
3496	}
3497	#endif
3498	if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3499	// try idiotic idea
3500	CbcObject * obj = new CbcIdiotBranch(this);
3501	obj->setPriority(1); // temp
3502	addObjects(1, &obj);
3503	delete obj;
3504	}
3505
3506	/*
3507	A hook to use clp to quickly explore some part of the tree.
3508	*/
3509	if (fastNodeDepth_ == 1000 &&/!parentModel_/(specialOptions_&2048) == 0) {
3510	fastNodeDepth_ = -1;
3511	CbcObject * obj =
3512	new CbcFollowOn(this);
3513	obj->setPriority(1);
3514	addObjects(1, &obj);
3515	delete obj;
3516	}
3517	int saveNumberSolves = numberSolves_;
3518	int saveNumberIterations = numberIterations_;
3519	if ((fastNodeDepth_ >= 0\|\|(moreSpecialOptions_&33554432)!=0)
3520	&&/!parentModel_/(specialOptions_&2048) == 0) {
3521	// add in a general depth object doClp
3522	int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3523	if ((moreSpecialOptions_&33554432)!=0)
3524	type=12;
3525	else
3526	fastNodeDepth_ += 1000000; // mark as done
3527	CbcObject * obj =
3528	new CbcGeneralDepth(this, type);
3529	addObjects(1, &obj);
3530	delete obj;
3531	// fake number of objects
3532	numberObjects_--;
3533	if (parallelMode() < -1) {
3534	// But make sure position is correct
3535	OsiObject * obj2 = object_[numberObjects_];
3536	obj = dynamic_cast<CbcObject *> (obj2);
3537	assert (obj);
3538	obj->setPosition(numberObjects_);
3539	}
3540	}
3541	#ifdef COIN_HAS_CLP
3542	#ifdef NO_CRUNCH
3543	if (true) {
3544	OsiClpSolverInterface * clpSolver
3545	= dynamic_cast<OsiClpSolverInterface *> (solver_);
3546	if (clpSolver && !parentModel_) {
3547	ClpSimplex * clpSimplex = clpSolver->getModelPtr();
3548	clpSimplex->setSpecialOptions(clpSimplex->specialOptions() \| 131072);
3549	//clpSimplex->startPermanentArrays();
3550	clpSimplex->setPersistenceFlag(2);
3551	}
3552	}
3553	#endif
3554	#endif
3555	// Save copy of solver
3556	OsiSolverInterface * saveSolver = NULL;
3557	if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
3558	saveSolver = solver_->clone();
3559	double checkCutoffForRestart = 1.0e100;
3560	saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3561	if ((specialOptions_&262144) != 0 && !parentModel_) {
3562	// Save stuff and return!
3563	storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3564	solver_->getColLower(),
3565	solver_->getColUpper());
3566	delete [] lowerBefore;
3567	delete [] upperBefore;
3568	delete saveSolver;
3569	if (flipObjective)
3570	flipModel();
3571	return;
3572	}
3573	/*
3574	We've taken the continuous relaxation as far as we can. Time to branch.
3575	The first order of business is to actually create a node. chooseBranch
3576	currently uses strong branching to evaluate branch object candidates,
3577	unless forced back to simple branching. If chooseBranch concludes that a
3578	branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3579	== -2) when forced to integer values, it returns here immediately.
3580
3581	Monotone variables trigger a call to resolve(). If the problem remains
3582	feasible, try again to choose a branching variable. At the end of the loop,
3583	resolved == true indicates that some variables were fixed.
3584
3585	Loss of feasibility will result in the deletion of newNode.
3586	*/
3587
3588	bool resolved = false ;
3589	CbcNode *newNode = NULL ;
3590	numberFixedAtRoot_ = 0;
3591	numberFixedNow_ = 0;
3592	if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
3593	#ifdef COIN_HAS_CLP
3594	OsiClpSolverInterface * clpSolver
3595	= dynamic_cast<OsiClpSolverInterface *> (solver_);
3596	if (clpSolver) {
3597	if (getCutoff()>1.0e20) {
3598	printf("Zapping costs\n");
3599	int numberColumns=solver_->getNumCols();
3600	double * zeroCost = new double [numberColumns];
3601	// could make small random
3602	memset(zeroCost,0,numberColumns*sizeof(double));
3603	solver_->setObjective(zeroCost);
3604	double objValue = solver_->getObjValue();
3605	solver_->setDblParam(OsiObjOffset,-objValue);
3606	clpSolver->getModelPtr()->setObjectiveValue(objValue);
3607	delete [] zeroCost;
3608	} else {
3609	moreSpecialOptions2_ &= ~2;
3610	}
3611	} else {
3612	#endif
3613	moreSpecialOptions2_ &= ~2;
3614	#ifdef COIN_HAS_CLP
3615	}
3616	#endif
3617	}
3618	int numberIterationsAtContinuous = numberIterations_;
3619	//solverCharacteristics_->setSolver(solver_);
3620	if (feasible) {
3621	// mark all cuts as globally valid
3622	int numberCuts=cuts.sizeRowCuts();
3623	resizeWhichGenerator(0,numberCuts);
3624	for (int i=0;i<numberCuts;i++) {
3625	cuts.rowCutPtr(i)->setGloballyValid();
3626	whichGenerator_[i]=20000+(whichGenerator_[i]%10000);
3627	}
3628	#define HOTSTART -1
3629	#if HOTSTART<0
3630	if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
3631	(moreSpecialOptions_&1024) != 0) {
3632	// Set priorities so only branch on ones we need to
3633	// use djs and see if only few branches needed
3634	#ifndef NDEBUG
3635	double integerTolerance = getIntegerTolerance() ;
3636	#endif
3637	bool possible = true;
3638	const double * saveLower = continuousSolver_->getColLower();
3639	const double * saveUpper = continuousSolver_->getColUpper();
3640	for (int i = 0; i < numberObjects_; i++) {
3641	const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
3642	if (thisOne) {
3643	int iColumn = thisOne->columnNumber();
3644	if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3645	possible = false;
3646	break;
3647	}
3648	} else {
3649	possible = false;
3650	break;
3651	}
3652	}
3653	if (possible) {
3654	OsiSolverInterface * solver = continuousSolver_->clone();
3655	int numberColumns = solver->getNumCols();
3656	for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3657	double value = bestSolution_[iColumn] ;
3658	value = CoinMax(value, saveLower[iColumn]) ;
3659	value = CoinMin(value, saveUpper[iColumn]) ;
3660	value = floor(value + 0.5);
3661	if (solver->isInteger(iColumn)) {
3662	solver->setColLower(iColumn, value);
3663	solver->setColUpper(iColumn, value);
3664	}
3665	}
3666	solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3667	// objlim and all slack
3668	double direction = solver->getObjSense();
3669	solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
3670	CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
3671	solver->setWarmStart(basis);
3672	delete basis;
3673	bool changed = true;
3674	hotstartPriorities_ = new int [numberColumns];
3675	for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3676	hotstartPriorities_[iColumn] = 1;
3677	while (changed) {
3678	changed = false;
3679	solver->resolve();
3680	if (!solver->isProvenOptimal()) {
3681	possible = false;
3682	break;
3683	}
3684	const double * dj = solver->getReducedCost();
3685	const double * colLower = solver->getColLower();
3686	const double * colUpper = solver->getColUpper();
3687	const double * solution = solver->getColSolution();
3688	int nAtLbNatural = 0;
3689	int nAtUbNatural = 0;
3690	int nZeroDj = 0;
3691	int nForced = 0;
3692	for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3693	double value = solution[iColumn] ;
3694	value = CoinMax(value, saveLower[iColumn]) ;
3695	value = CoinMin(value, saveUpper[iColumn]) ;
3696	if (solver->isInteger(iColumn)) {
3697	assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
3698	if (hotstartPriorities_[iColumn] == 1) {
3699	if (dj[iColumn] < -1.0e-6) {
3700	// negative dj
3701	if (saveUpper[iColumn] == colUpper[iColumn]) {
3702	nAtUbNatural++;
3703	hotstartPriorities_[iColumn] = 2;
3704	solver->setColLower(iColumn, saveLower[iColumn]);
3705	solver->setColUpper(iColumn, saveUpper[iColumn]);
3706	} else {
3707	nForced++;
3708	}
3709	} else if (dj[iColumn] > 1.0e-6) {
3710	// positive dj
3711	if (saveLower[iColumn] == colLower[iColumn]) {
3712	nAtLbNatural++;
3713	hotstartPriorities_[iColumn] = 2;
3714	solver->setColLower(iColumn, saveLower[iColumn]);
3715	solver->setColUpper(iColumn, saveUpper[iColumn]);
3716	} else {
3717	nForced++;
3718	}
3719	} else {
3720	// zero dj
3721	nZeroDj++;
3722	}
3723	}
3724	}
3725	}
3726	#if CBC_USEFUL_PRINTING>1
3727	printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3728	nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3729	#endif
3730	if (nAtLbNatural \|\| nAtUbNatural) {
3731	changed = true;
3732	} else {
3733	if (nForced + nZeroDj > 5000 \|\|
3734	(nForced + nZeroDj)*2 > numberIntegers_)
3735	possible = false;
3736	}
3737	}
3738	delete solver;
3739	}
3740	if (possible) {
3741	setHotstartSolution(bestSolution_);
3742	if (!saveCompare) {
3743	// create depth first comparison
3744	saveCompare = nodeCompare_;
3745	// depth first
3746	nodeCompare_ = new CbcCompareDepth();
3747	tree_->setComparison(*nodeCompare_) ;
3748	}
3749	} else {
3750	delete [] hotstartPriorities_;
3751	hotstartPriorities_ = NULL;
3752	}
3753	}
3754	#endif
3755	#if HOTSTART>0
3756	if (hotstartSolution_ && !hotstartPriorities_) {
3757	// Set up hot start
3758	OsiSolverInterface * solver = solver_->clone();
3759	double direction = solver_->getObjSense() ;
3760	int numberColumns = solver->getNumCols();
3761	double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3762	double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3763	// move solution
3764	solver->setColSolution(hotstartSolution_);
3765	// point to useful information
3766	const double * saveSolution = testSolution_;
3767	testSolution_ = solver->getColSolution();
3768	OsiBranchingInformation usefulInfo = usefulInformation();
3769	testSolution_ = saveSolution;
3770	/*
3771	Run through the objects and use feasibleRegion() to set variable bounds
3772	so as to fix the variables specified in the objects at their value in this
3773	solution. Since the object list contains (at least) one object for every
3774	integer variable, this has the effect of fixing all integer variables.
3775	*/
3776	for (int i = 0; i < numberObjects_; i++)
3777	object_[i]->feasibleRegion(solver, &usefulInfo);
3778	solver->resolve();
3779	assert (solver->isProvenOptimal());
3780	double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
3781	const double * dj = solver->getReducedCost();
3782	const double * colLower = solver->getColLower();
3783	const double * colUpper = solver->getColUpper();
3784	const double * solution = solver->getColSolution();
3785	int nAtLbNatural = 0;
3786	int nAtUbNatural = 0;
3787	int nAtLbNaturalZero = 0;
3788	int nAtUbNaturalZero = 0;
3789	int nAtLbFixed = 0;
3790	int nAtUbFixed = 0;
3791	int nAtOther = 0;
3792	int nAtOtherNatural = 0;
3793	int nNotNeeded = 0;
3794	delete [] hotstartSolution_;
3795	hotstartSolution_ = new double [numberColumns];
3796	delete [] hotstartPriorities_;
3797	hotstartPriorities_ = new int [numberColumns];
3798	int * order = (int *) saveUpper;
3799	int nFix = 0;
3800	double bestRatio = COIN_DBL_MAX;
3801	for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3802	double value = solution[iColumn] ;
3803	value = CoinMax(value, saveLower[iColumn]) ;
3804	value = CoinMin(value, saveUpper[iColumn]) ;
3805	double sortValue = COIN_DBL_MAX;
3806	if (solver->isInteger(iColumn)) {
3807	assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
3808	double value2 = floor(value + 0.5);
3809	if (dj[iColumn] < -1.0e-6) {
3810	// negative dj
3811	//assert (value2==colUpper[iColumn]);
3812	if (saveUpper[iColumn] == colUpper[iColumn]) {
3813	nAtUbNatural++;
3814	sortValue = 0.0;
3815	double value = -dj[iColumn];
3816	if (value > gap)
3817	nFix++;
3818	else if (gap < value*bestRatio)
3819	bestRatio = gap / value;
3820	if (saveLower[iColumn] != colLower[iColumn]) {
3821	nNotNeeded++;
3822	sortValue = 1.0e20;
3823	}
3824	} else if (saveLower[iColumn] == colUpper[iColumn]) {
3825	nAtLbFixed++;
3826	sortValue = dj[iColumn];
3827	} else {
3828	nAtOther++;
3829	sortValue = 0.0;
3830	if (saveLower[iColumn] != colLower[iColumn] &&
3831	saveUpper[iColumn] != colUpper[iColumn]) {
3832	nNotNeeded++;
3833	sortValue = 1.0e20;
3834	}
3835	}
3836	} else if (dj[iColumn] > 1.0e-6) {
3837	// positive dj
3838	//assert (value2==colLower[iColumn]);
3839	if (saveLower[iColumn] == colLower[iColumn]) {
3840	nAtLbNatural++;
3841	sortValue = 0.0;
3842	double value = dj[iColumn];
3843	if (value > gap)
3844	nFix++;
3845	else if (gap < value*bestRatio)
3846	bestRatio = gap / value;
3847	if (saveUpper[iColumn] != colUpper[iColumn]) {
3848	nNotNeeded++;
3849	sortValue = 1.0e20;
3850	}
3851	} else if (saveUpper[iColumn] == colLower[iColumn]) {
3852	nAtUbFixed++;
3853	sortValue = -dj[iColumn];
3854	} else {
3855	nAtOther++;
3856	sortValue = 0.0;
3857	if (saveLower[iColumn] != colLower[iColumn] &&
3858	saveUpper[iColumn] != colUpper[iColumn]) {
3859	nNotNeeded++;
3860	sortValue = 1.0e20;
3861	}
3862	}
3863	} else {
3864	// zero dj
3865	if (value2 == saveUpper[iColumn]) {
3866	nAtUbNaturalZero++;
3867	sortValue = 0.0;
3868	if (saveLower[iColumn] != colLower[iColumn]) {
3869	nNotNeeded++;
3870	sortValue = 1.0e20;
3871	}
3872	} else if (value2 == saveLower[iColumn]) {
3873	nAtLbNaturalZero++;
3874	sortValue = 0.0;
3875	} else {
3876	nAtOtherNatural++;
3877	sortValue = 0.0;
3878	if (saveLower[iColumn] != colLower[iColumn] &&
3879	saveUpper[iColumn] != colUpper[iColumn]) {
3880	nNotNeeded++;
3881	sortValue = 1.0e20;
3882	}
3883	}
3884	}
3885	#if HOTSTART==3
3886	sortValue = -fabs(dj[iColumn]);
3887	#endif
3888	}
3889	hotstartSolution_[iColumn] = value ;
3890	saveLower[iColumn] = sortValue;
3891	order[iColumn] = iColumn;
3892	}
3893	COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3894	nFix, bestRatio, gap));
3895	int nNeg = 0;
3896	CoinSort_2(saveLower, saveLower + numberColumns, order);
3897	for (int i = 0; i < numberColumns; i++) {
3898	if (saveLower[i] < 0.0) {
3899	nNeg++;
3900	#if HOTSTART==2\|\|HOTSTART==3
3901	// swap sign ?
3902	saveLower[i] = -saveLower[i];
3903	#endif
3904	}
3905	}
3906	CoinSort_2(saveLower, saveLower + nNeg, order);
3907	for (int i = 0; i < numberColumns; i++) {
3908	#if HOTSTART==1
3909	hotstartPriorities_[order[i]] = 100;
3910	#else
3911	hotstartPriorities_[order[i]] = -(i + 1);
3912	#endif
3913	}
3914	COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3915	nAtLbNatural,
3916	nAtUbNatural,
3917	nAtLbNaturalZero,
3918	nAtUbNaturalZero,
3919	nAtLbFixed,
3920	nAtUbFixed,
3921	nAtOther,
3922	nAtOtherNatural, nNotNeeded, nNeg));
3923	delete [] saveLower;
3924	delete [] saveUpper;
3925	if (!saveCompare) {
3926	// create depth first comparison
3927	saveCompare = nodeCompare_;
3928	// depth first
3929	nodeCompare_ = new CbcCompareDepth();
3930	tree_->setComparison(*nodeCompare_) ;
3931	}
3932	}
3933	#endif
3934	newNode = new CbcNode ;
3935	// Set objective value (not so obvious if NLP etc)
3936	setObjectiveValue(newNode, NULL);
3937	anyAction = -1 ;
3938	// To make depth available we may need a fake node
3939	CbcNode fakeNode;
3940	if (!currentNode_) {
3941	// Not true if sub trees assert (!numberNodes_);
3942	currentNode_ = &fakeNode;
3943	}
3944	phase_ = 3;
3945	// only allow 1000 passes
3946	int numberPassesLeft = 1000;
3947	// This is first crude step
3948	if (numberAnalyzeIterations_ && !parentModel_) {
3949	delete [] analyzeResults_;
3950	//int numberColumns = solver_->getNumCols();
3951	analyzeResults_ = new double [5*numberIntegers_];
3952	numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3953	if (numberFixedAtRoot_ > 0) {
3954	COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3955	setPointers(solver_);
3956	numberFixedNow_ = numberFixedAtRoot_;
3957	} else if (numberFixedAtRoot_ < 0) {
3958	COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3959	anyAction = -2;
3960	delete newNode ;
3961	newNode = NULL ;
3962	feasible = false ;
3963	}
3964	}
3965	OsiSolverBranch * branches = NULL;
3966	if (feasible)
3967	anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3968	NULL, NULL, NULL, branches);
3969	if (anyAction == -2 \|\| newNode->objectiveValue() >= cutoff) {
3970	if (anyAction != -2) {
3971	// zap parent nodeInfo
3972	#ifdef COIN_DEVELOP
3973	printf("zapping CbcNodeInfo %x\n", newNode->nodeInfo()->parent());
3974	#endif
3975	if (newNode->nodeInfo())
3976	newNode->nodeInfo()->nullParent();
3977	}
3978	delete newNode ;
3979	newNode = NULL ;
3980	feasible = false ;
3981	}
3982	}
3983	if (newNode && probingInfo_) {
3984	int number01 = probingInfo_->numberIntegers();
3985	//const fixEntry * entry = probingInfo_->fixEntries();
3986	const int * toZero = probingInfo_->toZero();
3987	//const int * toOne = probingInfo_->toOne();
3988	//const int * integerVariable = probingInfo_->integerVariable();
3989	if (toZero[number01]) {
3990	CglTreeProbingInfo info(*probingInfo_);
3991	if ((moreSpecialOptions2_&64)!=0&&!parentModel_) {
3992	/*
3993	Marginal idea. Further exploration probably good. Build some extra
3994	cliques from probing info. Not quite worth the effort?
3995	*/
3996	CglProbing generator1;
3997	generator1.setUsingObjective(false);
3998	generator1.setMaxPass(1);
3999	generator1.setMaxPassRoot(1);
4000	generator1.setMaxLook(100);
4001	generator1.setRowCuts(3);
4002	generator1.setMaxElements(300);
4003	generator1.setMaxProbeRoot(solver_->getNumCols());
4004	CoinThreadRandom randomGenerator;
4005	//CglTreeProbingInfo info(solver_);
4006	info.level = 0;
4007	info.formulation_rows = solver_->getNumRows();
4008	info.inTree = false;
4009	info.randomNumberGenerator=&randomGenerator;
4010	info.pass=4;
4011	generator1.setMode(8);
4012	OsiCuts cs;
4013	generator1.generateCutsAndModify(*solver_,cs,&info);
4014	// very clunky
4015	OsiSolverInterface * temp = generator1.cliqueModel(solver_,2);
4016	CglPreProcess dummy;
4017	OsiSolverInterface * newSolver=dummy.cliqueIt(*temp,0.0001);
4018	delete temp;
4019	OsiSolverInterface * fake = NULL;
4020	if (newSolver) {
4021	#if 0
4022	int numberCliques = generator1.numberCliques();
4023	cliqueEntry * entry = generator1.cliqueEntry();
4024	cliqueType * type = new cliqueType [numberCliques];
4025	int * start = new int [numberCliques+1];
4026	start[numberCliques]=2*numberCliques;
4027	int n=0;
4028	for (int i=0;i<numberCliques;i++) {
4029	start[i]=2*i;
4030	setOneFixesInCliqueEntry(entry[2*i],true);
4031	setOneFixesInCliqueEntry(entry[2*i+1],true);
4032	type[i]=0;
4033	}
4034	fake = info.analyze(*solver_, 1,numberCliques,start,
4035	entry,type);
4036	delete [] type;
4037	delete [] entry;
4038	#else
4039	fake = info.analyze(*newSolver, 1,-1);
4040	#endif
4041	delete newSolver;
4042	} else {
4043	fake = info.analyze(*solver_, 1);
4044	}
4045	if (fake) {
4046	//fake->writeMps("fake");
4047	CglFakeClique cliqueGen(fake);
4048	cliqueGen.setStarCliqueReport(false);
4049	cliqueGen.setRowCliqueReport(false);
4050	cliqueGen.setMinViolation(0.1);
4051	addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
4052	generator_[numberCutGenerators_-1]->setTiming(true);
4053	for (int i = 0; i < numberCutGenerators_; i++) {
4054	CglKnapsackCover * cutGen =
4055	dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
4056	if (cutGen) {
4057	cutGen->createCliques(*fake,2,200,false);
4058	}
4059	}
4060	}
4061	}
4062	if (probingInfo_->packDown()) {
4063	#if CBC_USEFUL_PRINTING>1
4064	printf("%d implications on %d 0-1\n", toZero[number01], number01);
4065	#endif
4066	// Create a cut generator that remembers implications discovered at root.
4067	CglImplication implication(probingInfo_);
4068	addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
4069	generator_[numberCutGenerators_-1]->setGlobalCuts(true);
4070	generator_[numberCutGenerators_-1]->setTiming(true);
4071	} else {
4072	delete probingInfo_;
4073	probingInfo_ = NULL;
4074	}
4075	} else {
4076	delete probingInfo_;
4077
4078	probingInfo_ = NULL;
4079	}
4080	}
4081	/*
4082	At this point, the root subproblem is infeasible or fathomed by bound
4083	(newNode == NULL), or we're live with an objective value that satisfies the
4084	current objective cutoff.
4085	*/
4086	assert (!newNode \|\| newNode->objectiveValue() <= cutoff) ;
4087	// Save address of root node as we don't want to delete it
4088	/*
4089	The common case is that the lp relaxation is feasible but doesn't satisfy
4090	integrality (i.e., newNode->branchingObject(), indicating we've been able to
4091	select a branching variable). Remove any cuts that have gone slack due to
4092	forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
4093	basis (via createInfo()) and stash the new cuts in the nodeInfo (via
4094	addCuts()). If, by some miracle, we have an integral solution at the root
4095	(newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
4096	a valid solution for use by setBestSolution().
4097	*/
4098	CoinWarmStartBasis *lastws = NULL ;
4099	if (feasible && newNode->branchingObject()) {
4100	if (resolved) {
4101	takeOffCuts(cuts, false, NULL) ;
4102	# ifdef CHECK_CUT_COUNTS
4103	{ printf("Number of rows after chooseBranch fix (root)"
4104	"(active only) %d\n",
4105	numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
4106	const CoinWarmStartBasis* debugws =
4107	dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4108	debugws->print() ;
4109	delete debugws ;
4110	}
4111	# endif
4112	}
4113	//newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
4114	//newNode->nodeInfo()->addCuts(cuts,
4115	// newNode->numberBranches(),whichGenerator_) ;
4116	if (lastws) delete lastws ;
4117	lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4118	}
4119	/*
4120	Continuous data to be used later
4121	*/
4122	continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4123	continuousInfeasibilities_ = 0 ;
4124	if (newNode) {
4125	continuousObjective_ = newNode->objectiveValue() ;
4126	delete [] continuousSolution_;
4127	continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4128	numberColumns);
4129	continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
4130	}
4131	/*
4132	Bound may have changed so reset in objects
4133	*/
4134	{
4135	int i ;
4136	for (i = 0; i < numberObjects_; i++)
4137	object_[i]->resetBounds(solver_) ;
4138	}
4139	/*
4140	Feasible? Then we should have either a live node prepped for future
4141	expansion (indicated by variable() >= 0), or (miracle of miracles) an
4142	integral solution at the root node.
4143
4144	initializeInfo sets the reference counts in the nodeInfo object. Since
4145	this node is still live, push it onto the heap that holds the live set.
4146	*/
4147	if (newNode) {
4148	if (newNode->branchingObject()) {
4149	newNode->initializeInfo() ;
4150	tree_->push(newNode) ;
4151	// save pointer to root node - so can pick up bounds
4152	if (!topOfTree_)
4153	topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
4154	if (statistics_) {
4155	if (numberNodes2_ == maximumStatistics_) {
4156	maximumStatistics_ = 2 * maximumStatistics_;
4157	CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
4158	memset(temp, 0, maximumStatistics_sizeof(CbcStatistics ));
4159	memcpy(temp, statistics_, numberNodes2_sizeof(CbcStatistics ));
4160	delete [] statistics_;
4161	statistics_ = temp;
4162	}
4163	assert (!statistics_[numberNodes2_]);
4164	statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4165	}
4166	numberNodes2_++;
4167	# ifdef CHECK_NODE
4168	printf("Node %x on tree\n", newNode) ;
4169	# endif
4170	} else {
4171	// continuous is integer
4172	double objectiveValue = newNode->objectiveValue();
4173	setBestSolution(CBC_SOLUTION, objectiveValue,
4174	solver_->getColSolution()) ;
4175	if (eventHandler) {
4176	// we are stopping anyway so no need to test return code
4177	eventHandler->event(CbcEventHandler::solution);
4178	}
4179	delete newNode ;
4180	newNode = NULL ;
4181	}
4182	}
4183
4184	if (printFrequency_ <= 0) {
4185	printFrequency_ = 1000 ;
4186	if (getNumCols() > 2000)
4187	printFrequency_ = 100 ;
4188	}
4189	/*
4190	It is possible that strong branching fixes one variable and then the code goes round
4191	again and again. This can take too long. So we need to warn user - just once.
4192	*/
4193	numberLongStrong_ = 0;
4194	CbcNode * createdNode = NULL;
4195	#ifdef CBC_THREAD
4196	if ((specialOptions_&2048) != 0)
4197	numberThreads_ = 0;
4198	if (numberThreads_ ) {
4199	nodeCompare_->sayThreaded(); // need to use addresses
4200	master_ = new CbcBaseModel(*this,
4201	(parallelMode() < -1) ? 1 : 0);
4202	masterThread_ = master_->masterThread();
4203	}
4204	#endif
4205	#ifdef COIN_HAS_CLP
4206	{
4207	OsiClpSolverInterface * clpSolver
4208	= dynamic_cast<OsiClpSolverInterface *> (solver_);
4209	if (clpSolver && !parentModel_) {
4210	clpSolver->computeLargestAway();
4211	}
4212	}
4213	#endif
4214	/*
4215	At last, the actual branch-and-cut search loop, which will iterate until
4216	the live set is empty or we hit some limit (integrality gap, time, node
4217	count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4218	solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4219	add the node to the live set.
4220
4221	The first action is to winnow the live set to remove nodes which are worse
4222	than the current objective cutoff.
4223	*/
4224	if (solver_->getRowCutDebuggerAlways()) {
4225	OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
4226	const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
4227	if (!debugger) {
4228	// infeasible!!
4229	printf("before search\n");
4230	const double * lower = solver_->getColLower();
4231	const double * upper = solver_->getColUpper();
4232	const double * solution = debuggerX->optimalSolution();
4233	int numberColumns = solver_->getNumCols();
4234	for (int i = 0; i < numberColumns; i++) {
4235	if (solver_->isInteger(i)) {
4236	if (solution[i] < lower[i] - 1.0e-6 \|\| solution[i] > upper[i] + 1.0e-6)
4237	printf("**** ");
4238	printf("%d %g <= %g <= %g\n",
4239	i, lower[i], solution[i], upper[i]);
4240	}
4241	}
4242	//abort();
4243	}
4244	}
4245	{
4246	// may be able to change cutoff now
4247	double cutoff = getCutoff();
4248	double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4249	if (cutoff > bestObjective_ - increment) {
4250	cutoff = bestObjective_ - increment ;
4251	setCutoff(cutoff) ;
4252	}
4253	}
4254	#ifdef CBC_THREAD
4255	bool goneParallel = false;
4256	#endif
4257	#define MAX_DEL_NODE 1
4258	CbcNode * delNode[MAX_DEL_NODE+1];
4259	int nDeleteNode = 0;
4260	// For Printing etc when parallel
4261	int lastEvery1000 = 0;
4262	int lastPrintEvery = 0;
4263	int numberConsecutiveInfeasible = 0;
4264	#define PERTURB_IN_FATHOM
4265	#ifdef PERTURB_IN_FATHOM
4266	// allow in fathom
4267	if ((moreSpecialOptions_& 262144) != 0)
4268	specialOptions_ \|= 131072;
4269	#endif
4270	while (true) {
4271	lockThread();
4272	#ifdef COIN_HAS_CLP
4273	// See if we want dantzig row choice
4274	goToDantzig(100, savePivotMethod);
4275	#endif
4276	if (tree_->empty()) {
4277	#ifdef CBC_THREAD
4278	if (parallelMode() > 0 && master_) {
4279	int anyLeft = master_->waitForThreadsInTree(0);
4280	if (!anyLeft) {
4281	master_->stopThreads(-1);
4282	break;
4283	}
4284	} else {
4285	break;
4286	}
4287	#else
4288	break;
4289	#endif
4290	} else {
4291	unlockThread();
4292	}
4293	// If done 50/100 nodes see if worth trying reduction
4294	if (numberNodes_ >= nextCheckRestart) {
4295	if (nextCheckRestart<100)
4296	nextCheckRestart=100;
4297	else
4298	nextCheckRestart=COIN_INT_MAX;
4299	#ifdef COIN_HAS_CLP
4300	OsiClpSolverInterface * clpSolver
4301	= dynamic_cast<OsiClpSolverInterface *> (solver_);
4302	if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4303	ClpSimplex * simplex = clpSolver->getModelPtr();
4304	int perturbation = simplex->perturbation();
4305	#if CBC_USEFUL_PRINTING>1
4306	printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4307	numberIterations_, saveNumberIterations,
4308	numberSolves_, saveNumberSolves, perturbation);
4309	#endif
4310	if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4311	8*(numberSolves_ - saveNumberSolves)) {
4312	// switch off perturbation
4313	simplex->setPerturbation(100);
4314	#if CBC_USEFUL_PRINTING>1
4315	printf("Perturbation switched off\n");
4316	#endif
4317	}
4318	}
4319	#endif
4320	/*
4321	Decide if we want to do a restart.
4322	*/
4323	if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4324	bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4325	(getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4326	int numberColumns = getNumCols();
4327	if (tryNewSearch) {
4328	// adding increment back allows current best - tiny bit weaker
4329	checkCutoffForRestart = getCutoff() + getCutoffIncrement() ;
4330	#if CBC_USEFUL_PRINTING>1
4331	printf("after %d nodes, cutoff %g - looking\n",
4332	numberNodes_, getCutoff());
4333	#endif
4334	saveSolver->resolve();
4335	double direction = saveSolver->getObjSense() ;
4336	double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4337	double tolerance;
4338	saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4339	if (gap <= 0.0)
4340	gap = tolerance;
4341	gap += 100.0 * tolerance;
4342	double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4343
4344	const double *lower = saveSolver->getColLower() ;
4345	const double *upper = saveSolver->getColUpper() ;
4346	const double *solution = saveSolver->getColSolution() ;
4347	const double *reducedCost = saveSolver->getReducedCost() ;
4348
4349	int numberFixed = 0 ;
4350	int numberFixed2 = 0;
4351	#ifdef COIN_DEVELOP
4352	printf("gap %g\n", gap);
4353	#endif
4354	for (int i = 0 ; i < numberIntegers_ ; i++) {
4355	int iColumn = integerVariable_[i] ;
4356	double djValue = direction * reducedCost[iColumn] ;
4357	if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4358	if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4359	//printf("%d to lb on dj of %g - bounds %g %g\n",
4360	// iColumn,djValue,lower[iColumn],upper[iColumn]);
4361	saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4362	numberFixed++ ;
4363	} else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4364	//printf("%d to ub on dj of %g - bounds %g %g\n",
4365	// iColumn,djValue,lower[iColumn],upper[iColumn]);
4366	saveSolver->setColLower(iColumn, upper[iColumn]) ;
4367	numberFixed++ ;
4368	}
4369	} else {
4370	//printf("%d has dj of %g - already fixed to %g\n",
4371	// iColumn,djValue,lower[iColumn]);
4372	numberFixed2++;
4373	}
4374	}
4375	#ifdef COIN_DEVELOP
4376	if ((specialOptions_&1) != 0) {
4377	const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4378	if (debugger) {
4379	printf("Contains optimal\n") ;
4380	OsiSolverInterface * temp = saveSolver->clone();
4381	const double * solution = debugger->optimalSolution();
4382	const double *lower = temp->getColLower() ;
4383	const double *upper = temp->getColUpper() ;
4384	int n = temp->getNumCols();
4385	for (int i = 0; i < n; i++) {
4386	if (temp->isInteger(i)) {
4387	double value = floor(solution[i] + 0.5);
4388	assert (value >= lower[i] && value <= upper[i]);
4389	temp->setColLower(i, value);
4390	temp->setColUpper(i, value);
4391	}
4392	}
4393	temp->writeMps("reduced_fix");
4394	delete temp;
4395	saveSolver->writeMps("reduced");
4396	} else {
4397	abort();
4398	}
4399	}
4400	printf("Restart could fix %d integers (%d already fixed)\n",
4401	numberFixed + numberFixed2, numberFixed2);
4402	#endif
4403	numberFixed += numberFixed2;
4404	if (numberFixed10 < numberColumns && numberFixed4 < numberIntegers_)
4405	tryNewSearch = false;
4406	}
4407	#ifdef CONFLICT_CUTS
4408	// temporary
4409	if ((moreSpecialOptions_&4194304)!=0)
4410	tryNewSearch=false;
4411	#endif
4412	if (tryNewSearch) {
4413	// back to solver without cuts?
4414	OsiSolverInterface * solver2 = saveSolver->clone();
4415	const double *lower = saveSolver->getColLower() ;
4416	const double *upper = saveSolver->getColUpper() ;
4417	for (int i = 0 ; i < numberIntegers_ ; i++) {
4418	int iColumn = integerVariable_[i] ;
4419	solver2->setColLower(iColumn, lower[iColumn]);
4420	solver2->setColUpper(iColumn, upper[iColumn]);
4421	}
4422	// swap
4423	delete saveSolver;
4424	saveSolver = solver2;
4425	double * newSolution = new double[numberColumns];
4426	double objectiveValue = checkCutoffForRestart;
4427	// Save the best solution so far.
4428	CbcSerendipity heuristic(*this);
4429	if (bestSolution_)
4430	heuristic.setInputSolution(bestSolution_, bestObjective_);
4431	// Magic number
4432	heuristic.setFractionSmall(0.8);
4433	// `pumpTune' to stand-alone solver for explanations.
4434	heuristic.setFeasibilityPumpOptions(1008013);
4435	// Use numberNodes to say how many are original rows
4436	heuristic.setNumberNodes(continuousSolver_->getNumRows());
4437	#ifdef COIN_DEVELOP
4438	if (continuousSolver_->getNumRows() <
4439	solver_->getNumRows())
4440	printf("%d rows added ZZZZZ\n",
4441	solver_->getNumRows() - continuousSolver_->getNumRows());
4442	#endif
4443	int returnCode = heuristic.smallBranchAndBound(saveSolver,
4444	-1, newSolution,
4445	objectiveValue,
4446	checkCutoffForRestart, "Reduce");
4447	if (returnCode < 0) {
4448	#ifdef COIN_DEVELOP
4449	printf("Restart - not small enough to do search after fixing\n");
4450	#endif
4451	delete [] newSolution;
4452	} else {
4453	// 1 for sol'n, 2 for finished, 3 for both
4454	if ((returnCode&1) != 0) {
4455	// increment number of solutions so other heuristics can test
4456	numberSolutions_++;
4457	numberHeuristicSolutions_++;
4458	lastHeuristic_ = NULL;
4459	setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4460	}
4461	delete [] newSolution;
4462	#ifdef CBC_THREAD
4463	if (master_) {
4464	lockThread();
4465	if (parallelMode() > 0) {
4466	while (master_->waitForThreadsInTree(0)) {
4467	lockThread();
4468	double dummyBest;
4469	tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4470	//unlockThread();
4471	}
4472	} else {
4473	double dummyBest;
4474	tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4475	}
4476	master_->waitForThreadsInTree(2);
4477	delete master_;
4478	master_ = NULL;
4479	masterThread_ = NULL;
4480	}
4481	#endif
4482	if (tree_->size()) {
4483	double dummyBest;
4484	tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4485	}
4486	break;
4487	}
4488	}
4489	delete saveSolver;
4490	saveSolver = NULL;
4491	}
4492	}
4493	/*
4494	Check for abort on limits: node count, solution count, time, integrality gap.
4495	*/
4496	if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4497	numberSolutions_ < intParam_[CbcMaxNumSol] &&
4498	!maximumSecondsReached() &&
4499	!stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 \|\|
4500	numberIterations_ < maximumNumberIterations_))) {
4501	// out of loop
4502	break;
4503	}
4504	#ifdef BONMIN
4505	assert(!solverCharacteristics_->solutionAddsCuts() \|\| solverCharacteristics_->mipFeasible());
4506	#endif
4507	// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4508	// we need to replace the comparison function to dive, and that requires reordering to retain the
4509	// heap property.
4510	#define DIVE_WHEN 1000
4511	#define DIVE_STOP 2000
4512	int kNode = numberNodes_ % 4000;
4513	if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4514	if (!parallelMode()) {
4515	if (kNode == DIVE_WHEN + 1 \|\| numberConsecutiveInfeasible > 1) {
4516	CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4517	(nodeCompare_);
4518	// Don't interfere if user has replaced the compare function.
4519	if (compare) {
4520	//printf("Redoing tree\n");
4521	compare->startDive(this);
4522	numberConsecutiveInfeasible = 0;
4523	}
4524	}
4525	}
4526	}
4527	// replace current cutoff?
4528	if (cutoff > getCutoff()) {
4529	double newCutoff = getCutoff();
4530	if (analyzeResults_) {
4531	// see if we could fix any (more)
4532	int n = 0;
4533	double * newLower = analyzeResults_;
4534	double * objLower = newLower + numberIntegers_;
4535	double * newUpper = objLower + numberIntegers_;
4536	double * objUpper = newUpper + numberIntegers_;
4537	for (int i = 0; i < numberIntegers_; i++) {
4538	if (objLower[i] > newCutoff) {
4539	n++;
4540	if (objUpper[i] > newCutoff) {
4541	newCutoff = -COIN_DBL_MAX;
4542	break;
4543	}
4544	// add as global cut
4545	objLower[i]=-COIN_DBL_MAX;
4546	OsiRowCut rc;
4547	rc.setLb(newLower[i]);
4548	rc.setUb(COIN_DBL_MAX);
4549	double one=1.0;
4550	rc.setRow(1,integerVariable_+i,&one,false);
4551	rc.setGloballyValidAsInteger(2);
4552	globalCuts_.addCutIfNotDuplicate(rc) ;
4553	} else if (objUpper[i] > newCutoff) {
4554	n++;
4555	// add as global cut
4556	objUpper[i]=-COIN_DBL_MAX;
4557	OsiRowCut rc;
4558	rc.setLb(-COIN_DBL_MAX);
4559	rc.setUb(newUpper[i]);
4560	double one=1.0;
4561	rc.setRow(1,integerVariable_+i,&one,false);
4562	rc.setGloballyValidAsInteger(2);
4563	globalCuts_.addCutIfNotDuplicate(rc) ;
4564	}
4565	}
4566	if (newCutoff == -COIN_DBL_MAX) {
4567	COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4568	} else if (n > numberFixedNow_) {
4569	COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4570	numberFixedNow_ = n;
4571	}
4572	}
4573	if (eventHandler) {
4574	if (!eventHandler->event(CbcEventHandler::solution)) {
4575	eventHappened_ = true; // exit
4576	}
4577	newCutoff = getCutoff();
4578	}
4579	lockThread();
4580	/*
4581	Clean the tree to reflect the new solution, then see if the
4582	node comparison predicate wants to make any changes. If so,
4583	call setComparison for the side effect of rebuilding the heap.
4584	*/
4585	tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4586	if (nodeCompare_->newSolution(this) \|\|
4587	nodeCompare_->newSolution(this,continuousObjective_,
4588	continuousInfeasibilities_)) {
4589	tree_->setComparison(*nodeCompare_) ;
4590	}
4591	if (tree_->empty()) {
4592	continue;
4593	}
4594	unlockThread();
4595	}
4596	cutoff = getCutoff() ;
4597	/*
4598	Periodic activities: Opportunities to
4599	+ tweak the nodeCompare criteria,
4600	+ check if we've closed the integrality gap enough to quit,
4601	+ print a summary line to let the user know we're working
4602	*/
4603	if (numberNodes_ >= lastEvery1000) {
4604	lockThread();
4605	#ifdef COIN_HAS_CLP
4606	// See if we want dantzig row choice
4607	goToDantzig(1000, savePivotMethod);
4608	#endif
4609	lastEvery1000 = numberNodes_ + 1000;
4610	bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4611	#ifdef CHECK_CUT_SIZE
4612	verifyCutSize (tree_, *this);
4613	#endif
4614	// redo tree if requested
4615	if (redoTree)
4616	tree_->setComparison(*nodeCompare_) ;
4617	unlockThread();
4618	}
4619	// Had hotstart before, now switched off
4620	if (saveCompare && !hotstartSolution_) {
4621	// hotstart switched off
4622	delete nodeCompare_; // off depth first
4623	nodeCompare_ = saveCompare;
4624	saveCompare = NULL;
4625	// redo tree
4626	lockThread();
4627	tree_->setComparison(*nodeCompare_) ;
4628	unlockThread();
4629	}
4630	if (numberNodes_ >= lastPrintEvery) {
4631	lastPrintEvery = numberNodes_ + printFrequency_;
4632	lockThread();
4633	int nNodes = tree_->size() ;
4634
4635	//MODIF PIERRE
4636	bestPossibleObjective_ = tree_->getBestPossibleObjective();
4637	#ifdef CBC_THREAD
4638	if (parallelMode() > 0 && master_) {
4639	// need to adjust for ones not on tree
4640	int numberThreads = master_->numberThreads();
4641	for (int i=0;i<numberThreads;i++) {
4642	CbcThread * child = master_->child(i);
4643	if (child->node()) {
4644	// adjust
4645	double value = child->node()->objectiveValue();
4646	bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4647	}
4648	}
4649	}
4650	#endif
4651	unlockThread();
4652	#if CBC_USEFUL_PRINTING>1
4653	if (getCutoff() < 1.0e20) {
4654	if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4655	!parentModel_)
4656	printf("model cutoff in status %g, best %g, increment %g\n",
4657	getCutoff(), bestObjective_, getCutoffIncrement());
4658	assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4659	1.0e-6 + 1.0e-10*fabs(bestObjective_));
4660	}
4661	#endif
4662	if (!intParam_[CbcPrinting]) {
4663	// Parallel may not have any nodes
4664	if (!nNodes)
4665	bestPossibleObjective_ = lastBestPossibleObjective;
4666	else
4667	lastBestPossibleObjective = bestPossibleObjective_;
4668	messageHandler()->message(CBC_STATUS, messages())
4669	<< numberNodes_ << CoinMax(nNodes,1) << bestObjective_ << bestPossibleObjective_
4670	<< getCurrentSeconds()
4671	<< CoinMessageEol ;
4672	} else if (intParam_[CbcPrinting] == 1) {
4673	messageHandler()->message(CBC_STATUS2, messages())
4674	<< numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4675	<< tree_->lastDepth() << tree_->lastUnsatisfied()
4676	<< tree_->lastObjective() << numberIterations_
4677	<< getCurrentSeconds()
4678	<< CoinMessageEol ;
4679	} else if (!numberExtraIterations_) {
4680	messageHandler()->message(CBC_STATUS2, messages())
4681	<< numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4682	<< tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4683	<< getCurrentSeconds()
4684	<< CoinMessageEol ;
4685	} else {
4686	messageHandler()->message(CBC_STATUS3, messages())
4687	<< numberNodes_ << numberFathoms_ << numberExtraNodes_ << nNodes
4688	<< bestObjective_ << bestPossibleObjective_
4689	<< tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4690	<< getCurrentSeconds()
4691	<< CoinMessageEol ;
4692	}
4693	#ifdef COIN_HAS_NTY
4694	if (symmetryInfo_)
4695	symmetryInfo_->statsOrbits(this,1);
4696	#endif
4697	#if PRINT_CONFLICT==1
4698	if (numberConflictCuts>lastNumberConflictCuts) {
4699	double length = lengthConflictCuts/numberConflictCuts;
4700	printf("%d new conflict cuts - total %d - average length %g\n",
4701	numberConflictCuts-lastNumberConflictCuts,
4702	numberConflictCuts,length);
4703	lastNumberConflictCuts = numberConflictCuts;
4704	}
4705	#endif
4706	if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4707	eventHappened_ = true; // exit
4708	}
4709	}
4710	// See if can stop on gap
4711	if(canStopOnGap()) {
4712	stoppedOnGap_ = true ;
4713	}
4714
4715	#ifdef CHECK_NODE_FULL
4716	verifyTreeNodes(tree_, *this) ;
4717	# endif
4718	# ifdef CHECK_CUT_COUNTS
4719	verifyCutCounts(tree_, *this) ;
4720	# endif
4721	/*
4722	Now we come to the meat of the loop. To create the active subproblem, we'll
4723	pop the most promising node in the live set, rebuild the subproblem it
4724	represents, and then execute the current arm of the branch to create the
4725	active subproblem.
4726	*/
4727	CbcNode * node = NULL;
4728	#ifdef CBC_THREAD
4729	if (!parallelMode() \|\| parallelMode() == -1) {
4730	#endif
4731	node = tree_->bestNode(cutoff) ;
4732	// Possible one on tree worse than cutoff
4733	// Weird comparison function can leave ineligible nodes on tree
4734	if (!node \|\| node->objectiveValue() > cutoff)
4735	continue;
4736	// Do main work of solving node here
4737	doOneNode(this, node, createdNode);
4738	#ifdef JJF_ZERO
4739	if (node) {
4740	if (createdNode) {
4741	printf("Node %d depth %d, created %d depth %d\n",
4742	node->nodeNumber(), node->depth(),
4743	createdNode->nodeNumber(), createdNode->depth());
4744	} else {
4745	printf("Node %d depth %d, no created node\n",
4746	node->nodeNumber(), node->depth());
4747	}
4748	} else if (createdNode) {
4749	printf("Node exhausted, created %d depth %d\n",
4750	createdNode->nodeNumber(), createdNode->depth());
4751	} else {
4752	printf("Node exhausted, no created node\n");
4753	numberConsecutiveInfeasible = 2;
4754	}
4755	#endif
4756	//if (createdNode)
4757	//numberConsecutiveInfeasible=0;
4758	//else
4759	//numberConsecutiveInfeasible++;
4760	#ifdef CBC_THREAD
4761	} else if (parallelMode() > 0) {
4762	//lockThread();
4763	//node = tree_->bestNode(cutoff) ;
4764	// Possible one on tree worse than cutoff
4765	if (true \|\| !node \|\| node->objectiveValue() > cutoff) {
4766	assert (master_);
4767	if (master_) {
4768	int anyLeft = master_->waitForThreadsInTree(1);
4769	// may need to go round again
4770	if (anyLeft) {
4771	continue;
4772	} else {
4773	master_->stopThreads(-1);
4774	}
4775	}
4776	}
4777	//unlockThread();
4778	} else {
4779	// Deterministic parallel
4780	if ((tree_->size() < CoinMax(numberThreads_, 8)\|\|
4781	hotstartSolution_) && !goneParallel) {
4782	node = tree_->bestNode(cutoff) ;
4783	// Possible one on tree worse than cutoff
4784	if (!node \|\| node->objectiveValue() > cutoff)
4785	continue;
4786	// Do main work of solving node here
4787	doOneNode(this, node, createdNode);
4788	assert (createdNode);
4789	if (!createdNode->active()) {
4790	delete createdNode;
4791	createdNode = NULL;
4792	} else {
4793	// Say one more pointing to this
4794	node->nodeInfo()->increment() ;
4795	tree_->push(createdNode) ;
4796	}
4797	if (node->active()) {
4798	assert (node->nodeInfo());
4799	if (node->nodeInfo()->numberBranchesLeft()) {
4800	tree_->push(node) ;
4801	} else {
4802	node->setActive(false);
4803	}
4804	} else {
4805	if (node->nodeInfo()) {
4806	if (!node->nodeInfo()->numberBranchesLeft())
4807	node->nodeInfo()->allBranchesGone(); // can clean up
4808	// So will delete underlying stuff
4809	node->setActive(true);
4810	}
4811	delNode[nDeleteNode++] = node;
4812	node = NULL;
4813	}
4814	if (nDeleteNode >= MAX_DEL_NODE) {
4815	for (int i = 0; i < nDeleteNode; i++) {
4816	//printf("trying to del %d %x\n",i,delNode[i]);
4817	delete delNode[i];
4818	//printf("done to del %d %x\n",i,delNode[i]);
4819	}
4820	nDeleteNode = 0;
4821	}
4822	} else {
4823	// Split and solve
4824	master_->deterministicParallel();
4825	goneParallel = true;
4826	}
4827	}
4828	#endif
4829	}
4830	if (nDeleteNode) {
4831	for (int i = 0; i < nDeleteNode; i++) {
4832	delete delNode[i];
4833	}
4834	nDeleteNode = 0;
4835	}
4836	#ifdef CBC_THREAD
4837	if (master_) {
4838	master_->stopThreads(-1);
4839	master_->waitForThreadsInTree(2);
4840	// adjust time to allow for children on some systems
4841	//dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4842	}
4843	#endif
4844	/*
4845	End of the non-abort actions. The next block of code is executed if we've
4846	aborted because we hit one of the limits. Clean up by deleting the live set
4847	and break out of the node processing loop. Note that on an abort, node may
4848	have been pushed back onto the tree for further processing, in which case
4849	it'll be deleted in cleanTree. We need to check.
4850	*/
4851	if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4852	numberSolutions_ < intParam_[CbcMaxNumSol] &&
4853	!maximumSecondsReached() &&
4854	!stoppedOnGap_ &&
4855	!eventHappened_ &&
4856	(maximumNumberIterations_ < 0 \|\| numberIterations_ < maximumNumberIterations_))
4857	) {
4858	if (tree_->size()) {
4859	double dummyBest;
4860	tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4861	#if 0 // Does not seem to be needed def CBC_THREAD
4862	if (parallelMode() > 0 && master_) {
4863	// see if any dangling nodes
4864	int numberThreads = master_->numberThreads();
4865	for (int i=0;i<numberThreads;i++) {
4866	CbcThread * child = master_->child(i);
4867	//if (child->createdNode())
4868	//printf("CHILD_NODE %p\n",child->createdNode());
4869	delete child->createdNode();
4870	}
4871	}
4872	#endif
4873	}
4874	delete nextRowCut_;
4875	/* order is important here:
4876	* maximumSecondsReached() should be checked before eventHappened_ and
4877	* isNodeLimitReached() should be checked after eventHappened_
4878	* reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4879	* and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4880	*/
4881	if (stoppedOnGap_) {
4882	messageHandler()->message(CBC_GAP, messages())
4883	<< bestObjective_ - bestPossibleObjective_
4884	<< dblParam_[CbcAllowableGap]
4885	<< dblParam_[CbcAllowableFractionGap]*100.0
4886	<< CoinMessageEol ;
4887	secondaryStatus_ = 2;
4888	status_ = 0 ;
4889	} else if (maximumSecondsReached()) {
4890	handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
4891	secondaryStatus_ = 4;
4892	status_ = 1 ;
4893	} else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4894	handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
4895	secondaryStatus_ = 6;
4896	status_ = 1 ;
4897	} else if (isNodeLimitReached()) {
4898	handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
4899	secondaryStatus_ = 3;
4900	status_ = 1 ;
4901	} else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4902	handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
4903	secondaryStatus_ = 8;
4904	status_ = 1 ;
4905	} else {
4906	handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
4907	secondaryStatus_ = 5;
4908	status_ = 5 ;
4909	}
4910	}
4911	#ifdef CBC_THREAD
4912	if (master_) {
4913	delete master_;
4914	master_ = NULL;
4915	masterThread_ = NULL;
4916	}
4917	#endif
4918	/*
4919	That's it, we've exhausted the search tree, or broken out of the loop because
4920	we hit some limit on evaluation.
4921
4922	We may have got an intelligent tree so give it one more chance
4923	*/
4924	// Tell solver we are not in Branch and Cut
4925	solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4926	tree_->endSearch();
4927	// If we did any sub trees - did we give up on any?
4928	if ( numberStoppedSubTrees_)
4929	status_ = 1;
4930	numberNodes_ += numberExtraNodes_;
4931	numberIterations_ += numberExtraIterations_;
4932	if (eventHandler) {
4933	eventHandler->event(CbcEventHandler::endSearch);
4934	}
4935	if (!status_) {
4936	// Set best possible unless stopped on gap
4937	if (secondaryStatus_ != 2)
4938	bestPossibleObjective_ = bestObjective_;
4939	handler_->message(CBC_END_GOOD, messages_)
4940	<< bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
4941	<< CoinMessageEol ;
4942	} else {
4943	handler_->message(CBC_END, messages_)
4944	<< bestObjective_ << bestPossibleObjective_
4945	<< numberIterations_ << numberNodes_ << getCurrentSeconds()
4946	<< CoinMessageEol ;
4947	}
4948	if ((moreSpecialOptions_&4194304)!=0) {
4949	// Conflict cuts
4950	int numberCuts = globalCuts_.sizeRowCuts();
4951	int nConflict=0;
4952	double sizeConflict = 0.0;
4953	for (int i=0;i<numberCuts;i++) {
4954	OsiRowCut2 * cut = globalCuts_.cut(i);
4955	if (cut->whichRow()==1) {
4956	nConflict++;
4957	sizeConflict += cut->row().getNumElements();
4958	}
4959	}
4960	if (nConflict) {
4961	sizeConflict /= nConflict;
4962	char general[200];
4963	sprintf(general, "%d conflict cuts generated - average length %g",
4964	nConflict,sizeConflict);
4965	messageHandler()->message(CBC_GENERAL,
4966	messages())
4967	<< general << CoinMessageEol ;
4968	}
4969	}
4970	if (numberStrongIterations_)
4971	handler_->message(CBC_STRONG_STATS, messages_)
4972	<< strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
4973	<< strongInfo_[1] << CoinMessageEol ;
4974	if (!numberExtraNodes_)
4975	handler_->message(CBC_OTHER_STATS, messages_)
4976	<< maximumDepthActual_
4977	<< numberDJFixed_ << CoinMessageEol ;
4978	else
4979	handler_->message(CBC_OTHER_STATS2, messages_)
4980	<< maximumDepthActual_
4981	<< numberDJFixed_ << numberFathoms_ << numberExtraNodes_ << numberExtraIterations_
4982	<< CoinMessageEol ;
4983	#ifdef COIN_HAS_NTY
4984	if (symmetryInfo_)
4985	symmetryInfo_->statsOrbits(this,1);
4986	#endif
4987	if (doStatistics == 100) {
4988	for (int i = 0; i < numberObjects_; i++) {
4989	CbcSimpleIntegerDynamicPseudoCost * obj =
4990	dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4991	if (obj)
4992	obj->print();
4993	}
4994	}
4995	if (statistics_) {
4996	// report in some way
4997	int * lookup = new int[numberObjects_];
4998	int i;
4999	for (i = 0; i < numberObjects_; i++)
5000	lookup[i] = -1;
5001	bool goodIds = false; //true;
5002	for (i = 0; i < numberObjects_; i++) {
5003	int iColumn = object_[i]->columnNumber();
5004	if (iColumn >= 0 && iColumn < numberColumns) {
5005	if (lookup[i] == -1) {
5006	lookup[i] = iColumn;
5007	} else {
5008	goodIds = false;
5009	break;
5010	}
5011	} else {
5012	goodIds = false;
5013	break;
5014	}
5015	}
5016	if (!goodIds) {
5017	delete [] lookup;
5018	lookup = NULL;
5019	}
5020	if (doStatistics >= 3) {
5021	printf(" node parent depth column value obj inf\n");
5022	for ( i = 0; i < numberNodes2_; i++) {
5023	statistics_[i]->print(lookup);
5024	}
5025	}
5026	if (doStatistics > 1) {
5027	// Find last solution
5028	int k;
5029	for (k = numberNodes2_ - 1; k >= 0; k--) {
5030	if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
5031	!statistics_[k]->endingInfeasibility())
5032	break;
5033	}
5034	if (k >= 0) {
5035	int depth = statistics_[k]->depth();
5036	int * which = new int[depth+1];
5037	for (i = depth; i >= 0; i--) {
5038	which[i] = k;
5039	k = statistics_[k]->parentNode();
5040	}
5041	printf(" node parent depth column value obj inf\n");
5042	for (i = 0; i <= depth; i++) {
5043	statistics_[which[i]]->print(lookup);
5044	}
5045	delete [] which;
5046	}
5047	}
5048	// now summary
5049	int maxDepth = 0;
5050	double averageSolutionDepth = 0.0;
5051	int numberSolutions = 0;
5052	double averageCutoffDepth = 0.0;
5053	double averageSolvedDepth = 0.0;
5054	int numberCutoff = 0;
5055	int numberDown = 0;
5056	int numberFirstDown = 0;
5057	double averageInfDown = 0.0;
5058	double averageObjDown = 0.0;
5059	int numberCutoffDown = 0;
5060	int numberUp = 0;
5061	int numberFirstUp = 0;
5062	double averageInfUp = 0.0;
5063	double averageObjUp = 0.0;
5064	int numberCutoffUp = 0;
5065	double averageNumberIterations1 = 0.0;
5066	double averageValue = 0.0;
5067	for ( i = 0; i < numberNodes2_; i++) {
5068	int depth = statistics_[i]->depth();
5069	int way = statistics_[i]->way();
5070	double value = statistics_[i]->value();
5071	double startingObjective = statistics_[i]->startingObjective();
5072	int startingInfeasibility = statistics_[i]->startingInfeasibility();
5073	double endingObjective = statistics_[i]->endingObjective();
5074	int endingInfeasibility = statistics_[i]->endingInfeasibility();
5075	maxDepth = CoinMax(depth, maxDepth);
5076	// Only for completed
5077	averageNumberIterations1 += statistics_[i]->numberIterations();
5078	averageValue += value;
5079	if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
5080	numberSolutions++;
5081	averageSolutionDepth += depth;
5082	}
5083	if (endingObjective == COIN_DBL_MAX) {
5084	numberCutoff++;
5085	averageCutoffDepth += depth;
5086	if (way < 0) {
5087	numberDown++;
5088	numberCutoffDown++;
5089	if (way == -1)
5090	numberFirstDown++;
5091	} else {
5092	numberUp++;
5093	numberCutoffUp++;
5094	if (way == 1)
5095	numberFirstUp++;
5096	}
5097	} else {
5098	averageSolvedDepth += depth;
5099	if (way < 0) {
5100	numberDown++;
5101	averageInfDown += startingInfeasibility - endingInfeasibility;
5102	averageObjDown += endingObjective - startingObjective;
5103	if (way == -1)
5104	numberFirstDown++;
5105	} else {
5106	numberUp++;
5107	averageInfUp += startingInfeasibility - endingInfeasibility;
5108	averageObjUp += endingObjective - startingObjective;
5109	if (way == 1)
5110	numberFirstUp++;
5111	}
5112	}
5113	}
5114	// Now print
5115	if (numberSolutions)
5116	averageSolutionDepth /= static_cast<double> (numberSolutions);
5117	int numberSolved = numberNodes2_ - numberCutoff;
5118	double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
5119	- numberIterationsAtContinuous;
5120	if (numberCutoff) {
5121	averageCutoffDepth /= static_cast<double> (numberCutoff);
5122	averageNumberIterations2 /= static_cast<double> (numberCutoff);
5123	}
5124	if (numberNodes2_)
5125	averageValue /= static_cast<double> (numberNodes2_);
5126	if (numberSolved) {
5127	averageNumberIterations1 /= static_cast<double> (numberSolved);
5128	averageSolvedDepth /= static_cast<double> (numberSolved);
5129	}
5130	printf("%d solution(s) were found (by branching) at an average depth of %g\n",
5131	numberSolutions, averageSolutionDepth);
5132	printf("average value of variable being branched on was %g\n",
5133	averageValue);
5134	printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
5135	numberCutoff, averageCutoffDepth, averageNumberIterations2);
5136	printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
5137	numberSolved, averageSolvedDepth, averageNumberIterations1);
5138	if (numberDown) {
5139	averageInfDown /= static_cast<double> (numberDown);
5140	averageObjDown /= static_cast<double> (numberDown);
5141	}
5142	printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5143	numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
5144	averageInfDown, averageObjDown);
5145	if (numberUp) {
5146	averageInfUp /= static_cast<double> (numberUp);
5147	averageObjUp /= static_cast<double> (numberUp);
5148	}
5149	printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5150	numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
5151	averageInfUp, averageObjUp);
5152	for ( i = 0; i < numberNodes2_; i++)
5153	delete statistics_[i];
5154	delete [] statistics_;
5155	statistics_ = NULL;
5156	maximumStatistics_ = 0;
5157	delete [] lookup;
5158	}
5159	/*
5160	If we think we have a solution, restore and confirm it with a call to
5161	setBestSolution(). We need to reset the cutoff value so as not to fathom
5162	the solution on bounds. Note that calling setBestSolution( ..., true)
5163	leaves the continuousSolver_ bounds vectors fixed at the solution value.
5164
5165	Running resolve() here is a failsafe --- setBestSolution has already
5166	reoptimised using the continuousSolver_. If for some reason we fail to
5167	prove optimality, run the problem again after instructing the solver to
5168	tell us more.
5169
5170	If all looks good, replace solver_ with continuousSolver_, so that the
5171	outside world will be able to obtain information about the solution using
5172	public methods.
5173
5174	Don't replace if we are trying to save cuts
5175	*/
5176	if (bestSolution_ && (solverCharacteristics_->solverType() < 2 \|\| solverCharacteristics_->solverType() == 4) &&
5177	((specialOptions_&8388608)==0\|\|(specialOptions_&2048)!=0)) {
5178	setCutoff(1.0e50) ; // As best solution should be worse than cutoff
5179	// change cutoff as constraint if wanted
5180	if (cutoffRowNumber_>=0) {
5181	if (solver_->getNumRows()>cutoffRowNumber_)
5182	solver_->setRowUpper(cutoffRowNumber_,1.0e50);
5183	}
5184	// also in continuousSolver_
5185	if (continuousSolver_) {
5186	// Solvers know about direction
5187	double direction = solver_->getObjSense();
5188	continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
5189	}
5190	phase_ = 5;
5191	double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
5192	if ((specialOptions_&4) == 0)
5193	bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
5194	setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
5195	continuousSolver_->resolve() ;
5196	if (!continuousSolver_->isProvenOptimal()) {
5197	continuousSolver_->messageHandler()->setLogLevel(2) ;
5198	continuousSolver_->initialSolve() ;
5199	}
5200	delete solver_ ;
5201	// above deletes solverCharacteristics_
5202	solverCharacteristics_ = NULL;
5203	solver_ = continuousSolver_ ;
5204	setPointers(solver_);
5205	continu