Context Navigation

CbcSolver.cpp @ 824

Last change on this file since 824 was 824, checked in by forrest, 13 years ago
for unit test return code
File size: 328.0 KB

Line
1	// Copyright (C) 2007, International Business Machines
2	// Corporation and others. All Rights Reserved.
3
4	#include "CbcConfig.h"
5	#include "CoinPragma.hpp"
6
7	#include <cassert>
8	#include <cstdio>
9	#include <cmath>
10	#include <cfloat>
11	#include <cstring>
12	#include <iostream>
13
14	#include "CoinPragma.hpp"
15	#include "CoinHelperFunctions.hpp"
16	// Version
17	#define CBCVERSION "2.00.00"
18
19	#include "CoinMpsIO.hpp"
20	#include "CoinModel.hpp"
21
22	#include "ClpFactorization.hpp"
23	#include "ClpQuadraticObjective.hpp"
24	#include "CoinTime.hpp"
25	#include "ClpSimplex.hpp"
26	#include "ClpSimplexOther.hpp"
27	#include "ClpSolve.hpp"
28	#include "ClpMessage.hpp"
29	#include "ClpPackedMatrix.hpp"
30	#include "ClpPlusMinusOneMatrix.hpp"
31	#include "ClpNetworkMatrix.hpp"
32	#include "ClpDualRowSteepest.hpp"
33	#include "ClpDualRowDantzig.hpp"
34	#include "ClpLinearObjective.hpp"
35	#include "ClpPrimalColumnSteepest.hpp"
36	#include "ClpPrimalColumnDantzig.hpp"
37	#include "ClpPresolve.hpp"
38	#include "CbcOrClpParam.hpp"
39	#include "OsiRowCutDebugger.hpp"
40	#include "OsiChooseVariable.hpp"
41	#include "OsiAuxInfo.hpp"
42	//#define USER_HAS_FAKE_CLP
43	//#define USER_HAS_FAKE_CBC
44	//#define CLP_MALLOC_STATISTICS
45	#ifdef CLP_MALLOC_STATISTICS
46	#include <malloc.h>
47	#include <exception>
48	#include <new>
49	static double malloc_times=0.0;
50	static double malloc_total=0.0;
51	static int malloc_amount[]={0,32,128,256,1024,4096,16384,65536,262144,INT_MAX};
52	static int malloc_n=10;
53	double malloc_counts[10]={0,0,0,0,0,0,0,0,0,0};
54	void * operator new (size_t size) throw (std::bad_alloc)
55	{
56	malloc_times ++;
57	malloc_total += size;
58	int i;
59	for (i=0;i<malloc_n;i++) {
60	if ((int) size<=malloc_amount[i]) {
61	malloc_counts[i]++;
62	break;
63	}
64	}
65	void * p =malloc(size);
66	//char * xx = (char *) p;
67	//memset(xx,0,size);
68	// Initialize random seed
69	//CoinSeedRandom(987654321);
70	return p;
71	}
72	void operator delete (void *p) throw()
73	{
74	free(p);
75	}
76	static void malloc_stats2()
77	{
78	double average = malloc_total/malloc_times;
79	printf("count %g bytes %g - average %g\n",malloc_times,malloc_total,average);
80	for (int i=0;i<malloc_n;i++)
81	printf("%g ",malloc_counts[i]);
82	printf("\n");
83	malloc_times=0.0;
84	malloc_total=0.0;
85	memset(malloc_counts,0,sizeof(malloc_counts));
86	}
87	#endif
88	//#define DMALLOC
89	#ifdef DMALLOC
90	#include "dmalloc.h"
91	#endif
92	#ifdef WSSMP_BARRIER
93	#define FOREIGN_BARRIER
94	#endif
95	#ifdef UFL_BARRIER
96	#define FOREIGN_BARRIER
97	#endif
98	#ifdef TAUCS_BARRIER
99	#define FOREIGN_BARRIER
100	#endif
101	#include "CoinWarmStartBasis.hpp"
102
103	#include "OsiSolverInterface.hpp"
104	#include "OsiCuts.hpp"
105	#include "OsiRowCut.hpp"
106	#include "OsiColCut.hpp"
107	#ifndef COIN_HAS_LINK
108	#define COIN_HAS_LINK
109	#endif
110	#ifdef COIN_HAS_LINK
111	#include "CbcLinked.hpp"
112	#endif
113	#include "CglPreProcess.hpp"
114	#include "CglCutGenerator.hpp"
115	#include "CglGomory.hpp"
116	#include "CglProbing.hpp"
117	#include "CglKnapsackCover.hpp"
118	#include "CglRedSplit.hpp"
119	#include "CglClique.hpp"
120	#include "CglFlowCover.hpp"
121	#include "CglMixedIntegerRounding2.hpp"
122	#include "CglTwomir.hpp"
123	#include "CglDuplicateRow.hpp"
124	#include "CglStored.hpp"
125	#include "CglLandP.hpp"
126	#include "CglResidualCapacity.hpp"
127
128	#include "CbcModel.hpp"
129	#include "CbcHeuristic.hpp"
130	#include "CbcHeuristicLocal.hpp"
131	#include "CbcHeuristicGreedy.hpp"
132	#include "CbcHeuristicFPump.hpp"
133	#include "CbcHeuristicRINS.hpp"
134	#include "CbcTreeLocal.hpp"
135	#include "CbcCompareActual.hpp"
136	#include "CbcBranchActual.hpp"
137	#include "CbcBranchLotsize.hpp"
138	#include "CbcOrClpParam.hpp"
139	#include "CbcCutGenerator.hpp"
140	#include "CbcStrategy.hpp"
141	#include "CbcBranchLotsize.hpp"
142
143	#include "OsiClpSolverInterface.hpp"
144	#include "CbcSolver.hpp"
145	CbcSolver::CbcSolver()
146	: babModel_(NULL),
147	userFunction_(NULL),
148	statusUserFunction_(NULL),
149	originalSolver_(NULL),
150	originalCoinModel_(NULL),
151	cutGenerator_(NULL),
152	numberUserFunctions_(0),
153	numberCutGenerators_(0),
154	startTime_(CoinCpuTime()),
155	parameters_(NULL),
156	numberParameters_(0),
157	doMiplib_(false),
158	noPrinting_(false),
159	readMode_(1)
160	{
161	callBack_ = new CbcStopNow();
162	fillParameters();
163	}
164	CbcSolver::CbcSolver(const OsiClpSolverInterface & solver)
165	: babModel_(NULL),
166	userFunction_(NULL),
167	statusUserFunction_(NULL),
168	originalSolver_(NULL),
169	originalCoinModel_(NULL),
170	cutGenerator_(NULL),
171	numberUserFunctions_(0),
172	numberCutGenerators_(0),
173	startTime_(CoinCpuTime()),
174	parameters_(NULL),
175	numberParameters_(0),
176	doMiplib_(false),
177	noPrinting_(false),
178	readMode_(1)
179	{
180	callBack_ = new CbcStopNow();
181	model_ = CbcModel(solver);
182	fillParameters();
183	}
184	CbcSolver::CbcSolver(const CbcModel & solver)
185	: babModel_(NULL),
186	userFunction_(NULL),
187	statusUserFunction_(NULL),
188	originalSolver_(NULL),
189	originalCoinModel_(NULL),
190	cutGenerator_(NULL),
191	numberUserFunctions_(0),
192	numberCutGenerators_(0),
193	startTime_(CoinCpuTime()),
194	parameters_(NULL),
195	numberParameters_(0),
196	doMiplib_(false),
197	noPrinting_(false),
198	readMode_(1)
199	{
200	callBack_ = new CbcStopNow();
201	model_ = solver;
202	fillParameters();
203	}
204	CbcSolver::~CbcSolver()
205	{
206	int i;
207	for (i=0;i<numberUserFunctions_;i++)
208	delete userFunction_[i];
209	delete [] userFunction_;
210	for (i=0;i<numberCutGenerators_;i++)
211	delete cutGenerator_[i];
212	delete [] cutGenerator_;
213	delete [] statusUserFunction_;
214	delete originalSolver_;
215	delete originalCoinModel_;
216	delete babModel_;
217	delete [] parameters_;
218	delete callBack_;
219	}
220	// Copy constructor
221	CbcSolver::CbcSolver ( const CbcSolver & rhs)
222	: model_(rhs.model_),
223	babModel_(NULL),
224	userFunction_(NULL),
225	statusUserFunction_(NULL),
226	numberUserFunctions_(rhs.numberUserFunctions_),
227	startTime_(CoinCpuTime()),
228	parameters_(NULL),
229	numberParameters_(rhs.numberParameters_),
230	doMiplib_(rhs.doMiplib_),
231	noPrinting_(rhs.noPrinting_),
232	readMode_(rhs.readMode_)
233	{
234	fillParameters();
235	if (rhs.babModel_)
236	babModel_ = new CbcModel(*rhs.babModel_);
237	userFunction_ = new CbcUser * [numberUserFunctions_];
238	int i;
239	for (i=0;i<numberUserFunctions_;i++)
240	userFunction_[i] = rhs.userFunction_[i]->clone();
241	for (i=0;i<numberParameters_;i++)
242	parameters_[i]=rhs.parameters_[i];
243	for (i=0;i<numberCutGenerators_;i++)
244	cutGenerator_[i] = rhs.cutGenerator_[i]->clone();
245	callBack_ = rhs.callBack_->clone();
246	originalSolver_ = NULL;
247	if (rhs.originalSolver_) {
248	OsiSolverInterface * temp = rhs.originalSolver_->clone();
249	originalSolver_ = dynamic_cast<OsiClpSolverInterface *> (temp);
250	assert (originalSolver_);
251	}
252	originalCoinModel_ = NULL;
253	if (rhs.originalCoinModel_)
254	originalCoinModel_ = new CoinModel(*rhs.originalCoinModel_);
255	}
256	// Assignment operator
257	CbcSolver &
258	CbcSolver::operator=(const CbcSolver & rhs)
259	{
260	if (this != &rhs) {
261	int i;
262	for (i=0;i<numberUserFunctions_;i++)
263	delete userFunction_[i];
264	delete [] userFunction_;
265	for (i=0;i<numberCutGenerators_;i++)
266	delete cutGenerator_[i];
267	delete [] cutGenerator_;
268	delete [] statusUserFunction_;
269	delete originalSolver_;
270	delete originalCoinModel_;
271	statusUserFunction_ = NULL;
272	delete babModel_;
273	delete [] parameters_;
274	delete callBack_;
275	numberUserFunctions_ = rhs.numberUserFunctions_;
276	startTime_ = rhs.startTime_;
277	numberParameters_= rhs.numberParameters_;
278	for (i=0;i<numberParameters_;i++)
279	parameters_[i]=rhs.parameters_[i];
280	for (i=0;i<numberCutGenerators_;i++)
281	cutGenerator_[i] = rhs.cutGenerator_[i]->clone();
282	noPrinting_ = rhs.noPrinting_;
283	readMode_ = rhs.readMode_;
284	doMiplib_ = rhs.doMiplib_;
285	model_ = rhs.model_;
286	if (rhs.babModel_)
287	babModel_ = new CbcModel(*rhs.babModel_);
288	else
289	babModel_ = NULL;
290	userFunction_ = new CbcUser * [numberUserFunctions_];
291	for (i=0;i<numberUserFunctions_;i++)
292	userFunction_[i] = rhs.userFunction_[i]->clone();
293	callBack_ = rhs.callBack_->clone();
294	originalSolver_ = NULL;
295	if (rhs.originalSolver_) {
296	OsiSolverInterface * temp = rhs.originalSolver_->clone();
297	originalSolver_ = dynamic_cast<OsiClpSolverInterface *> (temp);
298	assert (originalSolver_);
299	}
300	originalCoinModel_ = NULL;
301	if (rhs.originalCoinModel_)
302	originalCoinModel_ = new CoinModel(*rhs.originalCoinModel_);
303	}
304	return *this;
305	}
306	// Get int value
307	int CbcSolver::intValue(CbcOrClpParameterType type) const
308	{
309	return parameters_[whichParam(type,numberParameters_,parameters_)].intValue();
310	}
311	// Set int value
312	void CbcSolver::setIntValue(CbcOrClpParameterType type,int value)
313	{
314	parameters_[whichParam(type,numberParameters_,parameters_)].setIntValue(value);
315	}
316	// Get double value
317	double CbcSolver::doubleValue(CbcOrClpParameterType type) const
318	{
319	return parameters_[whichParam(type,numberParameters_,parameters_)].doubleValue();
320	}
321	// Set double value
322	void CbcSolver::setDoubleValue(CbcOrClpParameterType type,double value)
323	{
324	parameters_[whichParam(type,numberParameters_,parameters_)].setDoubleValue(value);
325	}
326	// User function (NULL if no match)
327	CbcUser * CbcSolver::userFunction(const char * name) const
328	{
329	int i;
330	for (i=0;i<numberUserFunctions_;i++) {
331	if (!strcmp(name,userFunction_[i]->name().c_str()))
332	break;
333	}
334	if (i<numberUserFunctions_)
335	return userFunction_[i];
336	else
337	return NULL;
338	}
339	void CbcSolver::fillParameters()
340	{
341	int maxParam = 200;
342	CbcOrClpParam * parameters = new CbcOrClpParam [maxParam];
343	numberParameters_=0 ;
344	establishParams(numberParameters_,parameters) ;
345	assert (numberParameters_<=maxParam);
346	parameters_ = new CbcOrClpParam [numberParameters_];
347	int i;
348	for (i=0;i<numberParameters_;i++)
349	parameters_[i]=parameters[i];
350	delete [] parameters;
351	const char dirsep = CoinFindDirSeparator();
352	std::string directory;
353	std::string dirSample;
354	std::string dirNetlib;
355	std::string dirMiplib;
356	if (dirsep == '/') {
357	directory = "./";
358	dirSample = "../../Data/Sample/";
359	dirNetlib = "../../Data/Netlib/";
360	dirMiplib = "../../Data/miplib3/";
361	} else {
362	directory = ".\\";
363	dirSample = "..\\..\\Data\\Sample\\";
364	dirNetlib = "..\\..\\Data\\Netlib\\";
365	dirMiplib = "..\\..\\Data\\miplib3\\";
366	}
367	std::string defaultDirectory = directory;
368	std::string importFile ="";
369	std::string exportFile ="default.mps";
370	std::string importBasisFile ="";
371	std::string importPriorityFile ="";
372	std::string debugFile="";
373	std::string printMask="";
374	std::string exportBasisFile ="default.bas";
375	std::string saveFile ="default.prob";
376	std::string restoreFile ="default.prob";
377	std::string solutionFile ="stdout";
378	std::string solutionSaveFile ="solution.file";
379	int doIdiot=-1;
380	int outputFormat=2;
381	int substitution=3;
382	int dualize=3;
383	int preSolve=5;
384	int doSprint=-1;
385	int testOsiParameters=-1;
386	int createSolver=0;
387	ClpSimplex * lpSolver;
388	OsiClpSolverInterface * clpSolver;
389	if (model_.solver()) {
390	clpSolver = dynamic_cast<OsiClpSolverInterface *> (model_.solver());
391	assert (clpSolver);
392	lpSolver = clpSolver->getModelPtr();
393	assert (lpSolver);
394	} else {
395	lpSolver = new ClpSimplex();
396	clpSolver = new OsiClpSolverInterface(lpSolver,true);
397	createSolver =1 ;
398	}
399	parameters_[whichParam(BASISIN,numberParameters_,parameters_)].setStringValue(importBasisFile);
400	parameters_[whichParam(PRIORITYIN,numberParameters_,parameters_)].setStringValue(importPriorityFile);
401	parameters_[whichParam(BASISOUT,numberParameters_,parameters_)].setStringValue(exportBasisFile);
402	parameters_[whichParam(DEBUG,numberParameters_,parameters_)].setStringValue(debugFile);
403	parameters_[whichParam(PRINTMASK,numberParameters_,parameters_)].setStringValue(printMask);
404	parameters_[whichParam(DIRECTORY,numberParameters_,parameters_)].setStringValue(directory);
405	parameters_[whichParam(DIRSAMPLE,numberParameters_,parameters_)].setStringValue(dirSample);
406	parameters_[whichParam(DIRNETLIB,numberParameters_,parameters_)].setStringValue(dirNetlib);
407	parameters_[whichParam(DIRMIPLIB,numberParameters_,parameters_)].setStringValue(dirMiplib);
408	parameters_[whichParam(DUALBOUND,numberParameters_,parameters_)].setDoubleValue(lpSolver->dualBound());
409	parameters_[whichParam(DUALTOLERANCE,numberParameters_,parameters_)].setDoubleValue(lpSolver->dualTolerance());
410	parameters_[whichParam(EXPORT,numberParameters_,parameters_)].setStringValue(exportFile);
411	parameters_[whichParam(IDIOT,numberParameters_,parameters_)].setIntValue(doIdiot);
412	parameters_[whichParam(IMPORT,numberParameters_,parameters_)].setStringValue(importFile);
413	parameters_[whichParam(PRESOLVETOLERANCE,numberParameters_,parameters_)].setDoubleValue(1.0e-8);
414	int iParam = whichParam(SOLVERLOGLEVEL,numberParameters_,parameters_);
415	int value=1;
416	clpSolver->messageHandler()->setLogLevel(1) ;
417	lpSolver->setLogLevel(1);
418	parameters_[iParam].setIntValue(value);
419	iParam = whichParam(LOGLEVEL,numberParameters_,parameters_);
420	model_.messageHandler()->setLogLevel(value);
421	parameters_[iParam].setIntValue(value);
422	parameters_[whichParam(MAXFACTOR,numberParameters_,parameters_)].setIntValue(lpSolver->factorizationFrequency());
423	parameters_[whichParam(MAXITERATION,numberParameters_,parameters_)].setIntValue(lpSolver->maximumIterations());
424	parameters_[whichParam(OUTPUTFORMAT,numberParameters_,parameters_)].setIntValue(outputFormat);
425	parameters_[whichParam(PRESOLVEPASS,numberParameters_,parameters_)].setIntValue(preSolve);
426	parameters_[whichParam(PERTVALUE,numberParameters_,parameters_)].setIntValue(lpSolver->perturbation());
427	parameters_[whichParam(PRIMALTOLERANCE,numberParameters_,parameters_)].setDoubleValue(lpSolver->primalTolerance());
428	parameters_[whichParam(PRIMALWEIGHT,numberParameters_,parameters_)].setDoubleValue(lpSolver->infeasibilityCost());
429	parameters_[whichParam(RESTORE,numberParameters_,parameters_)].setStringValue(restoreFile);
430	parameters_[whichParam(SAVE,numberParameters_,parameters_)].setStringValue(saveFile);
431	//parameters_[whichParam(TIMELIMIT,numberParameters_,parameters_)].setDoubleValue(1.0e8);
432	parameters_[whichParam(TIMELIMIT_BAB,numberParameters_,parameters_)].setDoubleValue(1.0e8);
433	parameters_[whichParam(SOLUTION,numberParameters_,parameters_)].setStringValue(solutionFile);
434	parameters_[whichParam(SAVESOL,numberParameters_,parameters_)].setStringValue(solutionSaveFile);
435	parameters_[whichParam(SPRINT,numberParameters_,parameters_)].setIntValue(doSprint);
436	parameters_[whichParam(SUBSTITUTION,numberParameters_,parameters_)].setIntValue(substitution);
437	parameters_[whichParam(DUALIZE,numberParameters_,parameters_)].setIntValue(dualize);
438	parameters_[whichParam(NUMBERBEFORE,numberParameters_,parameters_)].setIntValue(model_.numberBeforeTrust());
439	parameters_[whichParam(MAXNODES,numberParameters_,parameters_)].setIntValue(model_.getMaximumNodes());
440	parameters_[whichParam(STRONGBRANCHING,numberParameters_,parameters_)].setIntValue(model_.numberStrong());
441	parameters_[whichParam(INFEASIBILITYWEIGHT,numberParameters_,parameters_)].setDoubleValue(model_.getDblParam(CbcModel::CbcInfeasibilityWeight));
442	parameters_[whichParam(INTEGERTOLERANCE,numberParameters_,parameters_)].setDoubleValue(model_.getDblParam(CbcModel::CbcIntegerTolerance));
443	parameters_[whichParam(INCREMENT,numberParameters_,parameters_)].setDoubleValue(model_.getDblParam(CbcModel::CbcCutoffIncrement));
444	parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].setIntValue(testOsiParameters);
445	parameters_[whichParam(FPUMPTUNE,numberParameters_,parameters_)].setIntValue(1003);
446	#ifdef CBC_THREAD
447	parameters_[whichParam(THREADS,numberParameters_,parameters_)].setIntValue(0);
448	#endif
449	// Set up likely cut generators and defaults
450	parameters_[whichParam(PREPROCESS,numberParameters_,parameters_)].setCurrentOption("sos");
451	parameters_[whichParam(MIPOPTIONS,numberParameters_,parameters_)].setIntValue(128\|64\|1);
452	parameters_[whichParam(MIPOPTIONS,numberParameters_,parameters_)].setIntValue(1);
453	parameters_[whichParam(CUTPASSINTREE,numberParameters_,parameters_)].setIntValue(1);
454	parameters_[whichParam(MOREMIPOPTIONS,numberParameters_,parameters_)].setIntValue(-1);
455	parameters_[whichParam(MAXHOTITS,numberParameters_,parameters_)].setIntValue(100);
456	parameters_[whichParam(CUTSSTRATEGY,numberParameters_,parameters_)].setCurrentOption("on");
457	parameters_[whichParam(HEURISTICSTRATEGY,numberParameters_,parameters_)].setCurrentOption("on");
458	parameters_[whichParam(NODESTRATEGY,numberParameters_,parameters_)].setCurrentOption("fewest");
459	parameters_[whichParam(GOMORYCUTS,numberParameters_,parameters_)].setCurrentOption("ifmove");
460	parameters_[whichParam(PROBINGCUTS,numberParameters_,parameters_)].setCurrentOption("ifmove");
461	parameters_[whichParam(KNAPSACKCUTS,numberParameters_,parameters_)].setCurrentOption("ifmove");
462	parameters_[whichParam(REDSPLITCUTS,numberParameters_,parameters_)].setCurrentOption("off");
463	parameters_[whichParam(CLIQUECUTS,numberParameters_,parameters_)].setCurrentOption("ifmove");
464	parameters_[whichParam(MIXEDCUTS,numberParameters_,parameters_)].setCurrentOption("ifmove");
465	parameters_[whichParam(FLOWCUTS,numberParameters_,parameters_)].setCurrentOption("ifmove");
466	parameters_[whichParam(TWOMIRCUTS,numberParameters_,parameters_)].setCurrentOption("root");
467	parameters_[whichParam(LANDPCUTS,numberParameters_,parameters_)].setCurrentOption("off");
468	parameters_[whichParam(RESIDCUTS,numberParameters_,parameters_)].setCurrentOption("off");
469	parameters_[whichParam(ROUNDING,numberParameters_,parameters_)].setCurrentOption("on");
470	parameters_[whichParam(FPUMP,numberParameters_,parameters_)].setCurrentOption("on");
471	parameters_[whichParam(GREEDY,numberParameters_,parameters_)].setCurrentOption("on");
472	parameters_[whichParam(COMBINE,numberParameters_,parameters_)].setCurrentOption("on");
473	parameters_[whichParam(RINS,numberParameters_,parameters_)].setCurrentOption("off");
474	parameters_[whichParam(LOCALTREE,numberParameters_,parameters_)].setCurrentOption("off");
475	parameters_[whichParam(COSTSTRATEGY,numberParameters_,parameters_)].setCurrentOption("off");
476	if (createSolver)
477	delete clpSolver;
478	}
479	void CbcSolver::fillValuesInSolver()
480	{
481	OsiSolverInterface * solver = model_.solver();
482	OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
483	assert (clpSolver);
484	noPrinting_ = (clpSolver->getModelPtr()->logLevel()==0);
485	CoinMessageHandler * generalMessageHandler = clpSolver->messageHandler();
486	generalMessageHandler->setPrefix(true);
487	ClpSimplex * lpSolver = clpSolver->getModelPtr();
488	lpSolver->setPerturbation(50);
489	lpSolver->messageHandler()->setPrefix(false);
490	parameters_[whichParam(DUALBOUND,numberParameters_,parameters_)].setDoubleValue(lpSolver->dualBound());
491	parameters_[whichParam(DUALTOLERANCE,numberParameters_,parameters_)].setDoubleValue(lpSolver->dualTolerance());
492	int iParam = whichParam(SOLVERLOGLEVEL,numberParameters_,parameters_);
493	int value=parameters_[iParam].intValue();
494	clpSolver->messageHandler()->setLogLevel(value) ;
495	lpSolver->setLogLevel(value);
496	iParam = whichParam(LOGLEVEL,numberParameters_,parameters_);
497	value=parameters_[iParam].intValue();
498	model_.messageHandler()->setLogLevel(value);
499	parameters_[whichParam(LOGLEVEL,numberParameters_,parameters_)].setIntValue(model_.logLevel());
500	parameters_[whichParam(SOLVERLOGLEVEL,numberParameters_,parameters_)].setIntValue(lpSolver->logLevel());
501	parameters_[whichParam(MAXFACTOR,numberParameters_,parameters_)].setIntValue(lpSolver->factorizationFrequency());
502	parameters_[whichParam(MAXITERATION,numberParameters_,parameters_)].setIntValue(lpSolver->maximumIterations());
503	parameters_[whichParam(PERTVALUE,numberParameters_,parameters_)].setIntValue(lpSolver->perturbation());
504	parameters_[whichParam(PRIMALTOLERANCE,numberParameters_,parameters_)].setDoubleValue(lpSolver->primalTolerance());
505	parameters_[whichParam(PRIMALWEIGHT,numberParameters_,parameters_)].setDoubleValue(lpSolver->infeasibilityCost());
506	parameters_[whichParam(NUMBERBEFORE,numberParameters_,parameters_)].setIntValue(model_.numberBeforeTrust());
507	parameters_[whichParam(MAXNODES,numberParameters_,parameters_)].setIntValue(model_.getMaximumNodes());
508	parameters_[whichParam(STRONGBRANCHING,numberParameters_,parameters_)].setIntValue(model_.numberStrong());
509	parameters_[whichParam(INFEASIBILITYWEIGHT,numberParameters_,parameters_)].setDoubleValue(model_.getDblParam(CbcModel::CbcInfeasibilityWeight));
510	parameters_[whichParam(INTEGERTOLERANCE,numberParameters_,parameters_)].setDoubleValue(model_.getDblParam(CbcModel::CbcIntegerTolerance));
511	parameters_[whichParam(INCREMENT,numberParameters_,parameters_)].setDoubleValue(model_.getDblParam(CbcModel::CbcCutoffIncrement));
512	}
513	// Add user function
514	void
515	CbcSolver::addUserFunction(CbcUser * function)
516	{
517	CbcUser ** temp = new CbcUser * [numberUserFunctions_+1];
518	int i;
519	for (i=0;i<numberUserFunctions_;i++)
520	temp[i]=userFunction_[i];
521	delete [] userFunction_;
522	userFunction_ = temp;
523	userFunction_[numberUserFunctions_++] = function->clone();
524	delete [] statusUserFunction_;
525	statusUserFunction_ = NULL;
526	}
527	// Set user call back
528	void
529	CbcSolver::setUserCallBack(CbcStopNow * function)
530	{
531	delete callBack_;
532	callBack_ = function->clone();
533	}
534	// Copy of model on initial load (will contain output solutions)
535	void
536	CbcSolver::setOriginalSolver(OsiClpSolverInterface * originalSolver)
537	{
538	delete originalSolver_;
539	OsiSolverInterface * temp = originalSolver->clone();
540	originalSolver_ = dynamic_cast<OsiClpSolverInterface *> (temp);
541	assert (originalSolver_);
542
543	}
544	// Copy of model on initial load
545	void
546	CbcSolver::setOriginalCoinModel(CoinModel * originalCoinModel)
547	{
548	delete originalCoinModel_;
549	originalCoinModel_ = new CoinModel(*originalCoinModel);
550	}
551	// Add cut generator
552	void
553	CbcSolver::addCutGenerator(CglCutGenerator * generator)
554	{
555	CglCutGenerator ** temp = new CglCutGenerator * [numberCutGenerators_+1];
556	int i;
557	for (i=0;i<numberCutGenerators_;i++)
558	temp[i]=cutGenerator_[i];
559	delete [] cutGenerator_;
560	cutGenerator_ = temp;
561	cutGenerator_[numberCutGenerators_++] = generator->clone();
562	}
563	// User stuff (base class)
564	CbcUser::CbcUser()
565	: coinModel_(NULL),
566	userName_("null")
567	{
568	}
569	CbcUser::~CbcUser()
570	{
571	delete coinModel_;
572	}
573	// Copy constructor
574	CbcUser::CbcUser ( const CbcUser & rhs)
575	{
576	if (rhs.coinModel_)
577	coinModel_ = new CoinModel(*rhs.coinModel_);
578	else
579	coinModel_ = NULL;
580	userName_ = rhs.userName_;
581	}
582	// Assignment operator
583	CbcUser &
584	CbcUser::operator=(const CbcUser & rhs)
585	{
586	if (this != &rhs) {
587	if (rhs.coinModel_)
588	coinModel_ = new CoinModel(*rhs.coinModel_);
589	else
590	coinModel_ = NULL;
591	userName_ = rhs.userName_;
592	}
593	return *this;
594	}
595	/* Updates model_ from babModel_ according to returnMode
596	returnMode -
597	0 model and solver untouched - babModel updated
598	1 model updated - just with solution basis etc
599	2 model updated i.e. as babModel (babModel NULL) (only use without preprocessing!)
600	*/
601	void
602	CbcSolver::updateModel(ClpSimplex * model2, int returnMode)
603	{
604	if (!returnMode)
605	return;
606	if (returnMode==2&&babModel_)
607	model_ = *babModel_;
608	if (model2) {
609	// Only continuous valid
610	// update with basis etc
611	// map states
612	/* clp status
613	-1 - unknown e.g. before solve or if postSolve says not optimal
614	0 - optimal
615	1 - primal infeasible
616	2 - dual infeasible
617	3 - stopped on iterations or time
618	4 - stopped due to errors
619	5 - stopped by event handler (virtual int ClpEventHandler::event()) */
620	/* cbc status
621	-1 before branchAndBound
622	0 finished - check isProvenOptimal or isProvenInfeasible to see if solution found
623	(or check value of best solution)
624	1 stopped - on maxnodes, maxsols, maxtime
625	2 difficulties so run was abandoned
626	(5 event user programmed event occurred) */
627	/* clp secondary status of problem - may get extended
628	0 - none
629	1 - primal infeasible because dual limit reached OR probably primal
630	infeasible but can't prove it (main status 4)
631	2 - scaled problem optimal - unscaled problem has primal infeasibilities
632	3 - scaled problem optimal - unscaled problem has dual infeasibilities
633	4 - scaled problem optimal - unscaled problem has primal and dual infeasibilities
634	5 - giving up in primal with flagged variables
635	6 - failed due to empty problem check
636	7 - postSolve says not optimal
637	8 - failed due to bad element check
638	9 - status was 3 and stopped on time
639	100 up - translation of enum from ClpEventHandler
640	*/
641	/* cbc secondary status of problem
642	-1 unset (status_ will also be -1)
643	0 search completed with solution
644	1 linear relaxation not feasible (or worse than cutoff)
645	2 stopped on gap
646	3 stopped on nodes
647	4 stopped on time
648	5 stopped on user event
649	6 stopped on solutions
650	7 linear relaxation unbounded
651	*/
652	int iStatus = model2->status();
653	int iStatus2 = model2->secondaryStatus();
654	if (iStatus==0) {
655	iStatus2=0;
656	} else if (iStatus==1) {
657	iStatus=0;
658	iStatus2=1; // say infeasible
659	} else if (iStatus==2) {
660	iStatus=0;
661	iStatus2=7; // say unbounded
662	} else if (iStatus==3) {
663	iStatus=1;
664	if (iStatus2==9)
665	iStatus2=4;
666	else
667	iStatus2=3; // Use nodes - as closer than solutions
668	} else if (iStatus==4) {
669	iStatus=2; // difficulties
670	iStatus2=0;
671	}
672	model_.setProblemStatus(iStatus);
673	model_.setSecondaryStatus(iStatus2);
674	OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model_.solver());
675	ClpSimplex * lpSolver = clpSolver->getModelPtr();
676	if (model2!=lpSolver) {
677	lpSolver->moveInfo(*model2);
678	}
679	clpSolver->setWarmStart(NULL); // synchronize bases
680	if (originalSolver_) {
681	ClpSimplex * lpSolver2 = originalSolver_->getModelPtr();
682	assert (model2!=lpSolver2);
683	lpSolver2->moveInfo(*model2);
684	originalSolver_->setWarmStart(NULL); // synchronize bases
685	}
686	} else if (returnMode==1) {
687	OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model_.solver());
688	ClpSimplex * lpSolver = clpSolver->getModelPtr();
689	if (babModel_) {
690	model_.moveInfo(*babModel_);
691	int numberColumns = babModel_->getNumCols();
692	if (babModel_->bestSolution())
693	model_.setBestSolution(babModel_->bestSolution(),numberColumns,babModel_->getObjValue());
694	OsiClpSolverInterface * clpSolver1 = dynamic_cast< OsiClpSolverInterface*> (babModel_->solver());
695	ClpSimplex * lpSolver1 = clpSolver1->getModelPtr();
696	if (lpSolver1!=lpSolver&&model_.bestSolution()) {
697	lpSolver->moveInfo(*lpSolver1);
698	}
699	}
700	clpSolver->setWarmStart(NULL); // synchronize bases
701	}
702	if (returnMode==2) {
703	delete babModel_;
704	babModel_=NULL;
705	}
706	}
707	// Stop now stuff (base class)
708	CbcStopNow::CbcStopNow()
709	{
710	}
711	CbcStopNow::~CbcStopNow()
712	{
713	}
714	// Copy constructor
715	CbcStopNow::CbcStopNow ( const CbcStopNow & rhs)
716	{
717	}
718	// Assignment operator
719	CbcStopNow &
720	CbcStopNow::operator=(const CbcStopNow & rhs)
721	{
722	if (this != &rhs) {
723	}
724	return *this;
725	}
726	// Clone
727	CbcStopNow *
728	CbcStopNow::clone() const
729	{
730	return new CbcStopNow(*this);
731	}
732	//#undef NEW_STYLE_SOLVER
733	//#define NEW_STYLE_SOLVER
734	//#undef COIN_HAS_ASL
735	#ifdef COIN_HAS_ASL
736	#include "Cbc_ampl.h"
737	#endif
738	static double totalTime=0.0;
739	static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
740	static bool maskMatches(const int * starts, char ** masks,
741	std::string & check);
742	static void generateCode(CbcModel * model, const char * fileName,int type,int preProcess);
743	//#ifdef NDEBUG
744	//#undef NDEBUG
745	//#endif
746	// define (probably dummy fake main programs for UserClp and UserCbc
747	void fakeMain (ClpSimplex & model,OsiSolverInterface & osiSolver, CbcModel & babSolver);
748	void fakeMain2 (ClpSimplex & model,OsiClpSolverInterface & osiSolver,int options);
749
750	// Allow for interrupts
751	// But is this threadsafe ? (so switched off by option)
752
753	#include "CoinSignal.hpp"
754	static CbcModel * currentBranchModel = NULL;
755
756	extern "C" {
757	static void signal_handler(int whichSignal)
758	{
759	if (currentBranchModel!=NULL) {
760	currentBranchModel->setMaximumNodes(0); // stop at next node
761	currentBranchModel->setMaximumSeconds(0.0); // stop
762	}
763	return;
764	}
765	}
766	//#define CBC_SIG_TRAP
767	#ifdef CBC_SIG_TRAP
768	#include <setjmp.h>
769	static sigjmp_buf cbc_seg_buffer;
770	extern "C" {
771	static void signal_handler_error(int whichSignal)
772	{
773	siglongjmp(cbc_seg_buffer,1);
774	}
775	}
776	#endif
777	#if 0
778	/* Updates model_ from babModel_ according to returnMode
779	returnMode -
780	0 model and solver untouched
781	1 model updated - just with solution basis etc
782	*/
783	static void updateModel(CbcModel & model_,int returnMode)
784	{
785	if (!returnMode)
786	return;
787	assert (returnMode==1);
788	}
789	#endif
790	int CbcOrClpRead_mode=1;
791	FILE * CbcOrClpReadCommand=stdin;
792	static bool noPrinting=false;
793	#ifdef NEW_STYLE_SOLVER
794	int * CbcSolver::analyze(OsiClpSolverInterface * solverMod, int & numberChanged, double & increment,
795	bool changeInt, CoinMessageHandler * generalMessageHandler)
796	#else
797	static int * analyze(OsiClpSolverInterface * solverMod, int & numberChanged, double & increment,
798	bool changeInt, CoinMessageHandler * generalMessageHandler)
799	#endif
800	{
801	#ifndef NEW_STYLE_SOLVER
802	bool noPrinting_ = noPrinting;
803	#endif
804	OsiSolverInterface * solver = solverMod->clone();
805	char generalPrint[200];
806	if (0) {
807	// just get increment
808	CbcModel model(*solver);
809	model.analyzeObjective();
810	double increment2=model.getCutoffIncrement();
811	printf("initial cutoff increment %g\n",increment2);
812	}
813	const double *objective = solver->getObjCoefficients() ;
814	const double *lower = solver->getColLower() ;
815	const double *upper = solver->getColUpper() ;
816	int numberColumns = solver->getNumCols() ;
817	int numberRows = solver->getNumRows();
818	double direction = solver->getObjSense();
819	int iRow,iColumn;
820
821	// Row copy
822	CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
823	const double * elementByRow = matrixByRow.getElements();
824	const int * column = matrixByRow.getIndices();
825	const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
826	const int * rowLength = matrixByRow.getVectorLengths();
827
828	// Column copy
829	CoinPackedMatrix matrixByCol(*solver->getMatrixByCol());
830	const double * element = matrixByCol.getElements();
831	const int * row = matrixByCol.getIndices();
832	const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
833	const int * columnLength = matrixByCol.getVectorLengths();
834
835	const double * rowLower = solver->getRowLower();
836	const double * rowUpper = solver->getRowUpper();
837
838	char * ignore = new char [numberRows];
839	int * changed = new int[numberColumns];
840	int * which = new int[numberRows];
841	double * changeRhs = new double[numberRows];
842	memset(changeRhs,0,numberRows*sizeof(double));
843	memset(ignore,0,numberRows);
844	numberChanged=0;
845	int numberInteger=0;
846	for (iColumn=0;iColumn<numberColumns;iColumn++) {
847	if (upper[iColumn] > lower[iColumn]+1.0e-8&&solver->isInteger(iColumn))
848	numberInteger++;
849	}
850	bool finished=false;
851	while (!finished) {
852	int saveNumberChanged = numberChanged;
853	for (iRow=0;iRow<numberRows;iRow++) {
854	int numberContinuous=0;
855	double value1=0.0,value2=0.0;
856	bool allIntegerCoeff=true;
857	double sumFixed=0.0;
858	int jColumn1=-1,jColumn2=-1;
859	for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
860	int jColumn = column[j];
861	double value = elementByRow[j];
862	if (upper[jColumn] > lower[jColumn]+1.0e-8) {
863	if (!solver->isInteger(jColumn)) {
864	if (numberContinuous==0) {
865	jColumn1=jColumn;
866	value1=value;
867	} else {
868	jColumn2=jColumn;
869	value2=value;
870	}
871	numberContinuous++;
872	} else {
873	if (fabs(value-floor(value+0.5))>1.0e-12)
874	allIntegerCoeff=false;
875	}
876	} else {
877	sumFixed += lower[jColumn]*value;
878	}
879	}
880	double low = rowLower[iRow];
881	if (low>-1.0e20) {
882	low -= sumFixed;
883	if (fabs(low-floor(low+0.5))>1.0e-12)
884	allIntegerCoeff=false;
885	}
886	double up = rowUpper[iRow];
887	if (up<1.0e20) {
888	up -= sumFixed;
889	if (fabs(up-floor(up+0.5))>1.0e-12)
890	allIntegerCoeff=false;
891	}
892	if (!allIntegerCoeff)
893	continue; // can't do
894	if (numberContinuous==1) {
895	// see if really integer
896	// This does not allow for complicated cases
897	if (low==up) {
898	if (fabs(value1)>1.0e-3) {
899	value1 = 1.0/value1;
900	if (fabs(value1-floor(value1+0.5))<1.0e-12) {
901	// integer
902	changed[numberChanged++]=jColumn1;
903	solver->setInteger(jColumn1);
904	if (upper[jColumn1]>1.0e20)
905	solver->setColUpper(jColumn1,1.0e20);
906	if (lower[jColumn1]<-1.0e20)
907	solver->setColLower(jColumn1,-1.0e20);
908	}
909	}
910	} else {
911	if (fabs(value1)>1.0e-3) {
912	value1 = 1.0/value1;
913	if (fabs(value1-floor(value1+0.5))<1.0e-12) {
914	// This constraint will not stop it being integer
915	ignore[iRow]=1;
916	}
917	}
918	}
919	} else if (numberContinuous==2) {
920	if (low==up) {
921	/* need general theory - for now just look at 2 cases -
922	1 - +- 1 one in column and just costs i.e. matching objective
923	2 - +- 1 two in column but feeds into G/L row which will try and minimize
924	*/
925	if (fabs(value1)==1.0&&value1*value2==-1.0&&!lower[jColumn1]
926	&&!lower[jColumn2]) {
927	int n=0;
928	int i;
929	double objChange=direction*(objective[jColumn1]+objective[jColumn2]);
930	double bound = CoinMin(upper[jColumn1],upper[jColumn2]);
931	bound = CoinMin(bound,1.0e20);
932	for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
933	int jRow = row[i];
934	double value = element[i];
935	if (jRow!=iRow) {
936	which[n++]=jRow;
937	changeRhs[jRow]=value;
938	}
939	}
940	for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
941	int jRow = row[i];
942	double value = element[i];
943	if (jRow!=iRow) {
944	if (!changeRhs[jRow]) {
945	which[n++]=jRow;
946	changeRhs[jRow]=value;
947	} else {
948	changeRhs[jRow]+=value;
949	}
950	}
951	}
952	if (objChange>=0.0) {
953	// see if all rows OK
954	bool good=true;
955	for (i=0;i<n;i++) {
956	int jRow = which[i];
957	double value = changeRhs[jRow];
958	if (value) {
959	value *= bound;
960	if (rowLength[jRow]==1) {
961	if (value>0.0) {
962	double rhs = rowLower[jRow];
963	if (rhs>0.0) {
964	double ratio =rhs/value;
965	if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
966	good=false;
967	}
968	} else {
969	double rhs = rowUpper[jRow];
970	if (rhs<0.0) {
971	double ratio =rhs/value;
972	if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
973	good=false;
974	}
975	}
976	} else if (rowLength[jRow]==2) {
977	if (value>0.0) {
978	if (rowLower[jRow]>-1.0e20)
979	good=false;
980	} else {
981	if (rowUpper[jRow]<1.0e20)
982	good=false;
983	}
984	} else {
985	good=false;
986	}
987	}
988	}
989	if (good) {
990	// both can be integer
991	changed[numberChanged++]=jColumn1;
992	solver->setInteger(jColumn1);
993	if (upper[jColumn1]>1.0e20)
994	solver->setColUpper(jColumn1,1.0e20);
995	if (lower[jColumn1]<-1.0e20)
996	solver->setColLower(jColumn1,-1.0e20);
997	changed[numberChanged++]=jColumn2;
998	solver->setInteger(jColumn2);
999	if (upper[jColumn2]>1.0e20)
1000	solver->setColUpper(jColumn2,1.0e20);
1001	if (lower[jColumn2]<-1.0e20)
1002	solver->setColLower(jColumn2,-1.0e20);
1003	}
1004	}
1005	// clear
1006	for (i=0;i<n;i++) {
1007	changeRhs[which[i]]=0.0;
1008	}
1009	}
1010	}
1011	}
1012	}
1013	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1014	if (upper[iColumn] > lower[iColumn]+1.0e-8&&!solver->isInteger(iColumn)) {
1015	double value;
1016	value = upper[iColumn];
1017	if (value<1.0e20&&fabs(value-floor(value+0.5))>1.0e-12)
1018	continue;
1019	value = lower[iColumn];
1020	if (value>-1.0e20&&fabs(value-floor(value+0.5))>1.0e-12)
1021	continue;
1022	bool integer=true;
1023	for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
1024	int iRow = row[j];
1025	if (!ignore[iRow]) {
1026	integer=false;
1027	break;
1028	}
1029	}
1030	if (integer) {
1031	// integer
1032	changed[numberChanged++]=iColumn;
1033	solver->setInteger(iColumn);
1034	if (upper[iColumn]>1.0e20)
1035	solver->setColUpper(iColumn,1.0e20);
1036	if (lower[iColumn]<-1.0e20)
1037	solver->setColLower(iColumn,-1.0e20);
1038	}
1039	}
1040	}
1041	finished = numberChanged==saveNumberChanged;
1042	}
1043	delete [] which;
1044	delete [] changeRhs;
1045	delete [] ignore;
1046	//if (numberInteger&&!noPrinting_)
1047	//printf("%d integer variables",numberInteger);
1048	if (changeInt) {
1049	//if (!noPrinting_) {
1050	//if (numberChanged)
1051	// printf(" and %d variables made integer\n",numberChanged);
1052	//else
1053	// printf("\n");
1054	//}
1055	delete [] ignore;
1056	//increment=0.0;
1057	if (!numberChanged) {
1058	delete [] changed;
1059	delete solver;
1060	return NULL;
1061	} else {
1062	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1063	if (solver->isInteger(iColumn))
1064	solverMod->setInteger(iColumn);
1065	}
1066	delete solver;
1067	return changed;
1068	}
1069	} else {
1070	//if (!noPrinting_) {
1071	//if (numberChanged)
1072	// printf(" and %d variables could be made integer\n",numberChanged);
1073	//else
1074	// printf("\n");
1075	//}
1076	// just get increment
1077	int logLevel=generalMessageHandler->logLevel();
1078	CbcModel model(*solver);
1079	model.passInMessageHandler(generalMessageHandler);
1080	if (noPrinting_)
1081	model.setLogLevel(0);
1082	model.analyzeObjective();
1083	generalMessageHandler->setLogLevel(logLevel);
1084	double increment2=model.getCutoffIncrement();
1085	if (increment2>increment&&increment2>0.0) {
1086	if (!noPrinting_) {
1087	sprintf(generalPrint,"Cutoff increment increased from %g to %g",increment,increment2);
1088	CoinMessages generalMessages = solverMod->getModelPtr()->messages();
1089	generalMessageHandler->message(CLP_GENERAL,generalMessages)
1090	<< generalPrint
1091	<<CoinMessageEol;
1092	}
1093	increment=increment2;
1094	}
1095	delete solver;
1096	numberChanged=0;
1097	delete [] changed;
1098	return NULL;
1099	}
1100	}
1101	#if 1
1102	#include "ClpSimplexOther.hpp"
1103
1104	// Crunch down model
1105	static void
1106	crunchIt(ClpSimplex * model)
1107	{
1108	#if 0
1109	model->dual();
1110	#else
1111	int numberColumns = model->numberColumns();
1112	int numberRows = model->numberRows();
1113	// Use dual region
1114	double * rhs = model->dualRowSolution();
1115	int * whichRow = new int[3*numberRows];
1116	int * whichColumn = new int[2*numberColumns];
1117	int nBound;
1118	ClpSimplex * small = ((ClpSimplexOther *) model)->crunch(rhs,whichRow,whichColumn,
1119	nBound,false,false);
1120	if (small) {
1121	small->dual();
1122	if (small->problemStatus()==0) {
1123	model->setProblemStatus(0);
1124	((ClpSimplexOther ) model)->afterCrunch(small,whichRow,whichColumn,nBound);
1125	} else if (small->problemStatus()!=3) {
1126	model->setProblemStatus(1);
1127	} else {
1128	if (small->problemStatus()==3) {
1129	// may be problems
1130	small->computeObjectiveValue();
1131	model->setObjectiveValue(small->objectiveValue());
1132	model->setProblemStatus(3);
1133	} else {
1134	model->setProblemStatus(3);
1135	}
1136	}
1137	delete small;
1138	} else {
1139	model->setProblemStatus(1);
1140	}
1141	delete [] whichRow;
1142	delete [] whichColumn;
1143	#endif
1144	}
1145	/*
1146	On input
1147	doAction - 0 just fix in original and return NULL
1148	1 return fixed non-presolved solver
1149	2 as one but use presolve Inside this
1150	3 do heuristics and set best solution
1151	4 do BAB and just set best solution
1152	10+ then use lastSolution and relax a few
1153	-2 cleanup afterwards if using 2
1154	On output - number fixed
1155	*/
1156	static OsiClpSolverInterface * fixVubs(CbcModel & model, int skipZero2,
1157	int & doAction, CoinMessageHandler * generalMessageHandler,
1158	const double * lastSolution, double dextra[5])
1159	{
1160	if (doAction==11&&!lastSolution)
1161	lastSolution = model.bestSolution();
1162	assert (((doAction==1\|\|doAction==2)&&!lastSolution)\|\|(doAction==11&&lastSolution));
1163	double fractionIntFixed = dextra[3];
1164	double fractionFixed = dextra[4];
1165	double time1 = CoinCpuTime();
1166	OsiSolverInterface * originalSolver = model.solver();
1167	OsiClpSolverInterface * originalClpSolver = dynamic_cast< OsiClpSolverInterface*> (originalSolver);
1168	ClpSimplex * originalLpSolver = originalClpSolver->getModelPtr();
1169	int * originalColumns=NULL;
1170	OsiClpSolverInterface * clpSolver;
1171	ClpSimplex * lpSolver;
1172	ClpPresolve pinfo;
1173	if (doAction==2) {
1174	doAction=1;
1175	lpSolver = pinfo.presolvedModel(*originalLpSolver,1.0e-8,true,10);
1176	assert (lpSolver);
1177	clpSolver = new OsiClpSolverInterface(lpSolver,true);
1178	assert(lpSolver == clpSolver->getModelPtr());
1179	//int numberColumns = lpSolver->numberColumns();
1180	int numberColumns = lpSolver->numberColumns();
1181	originalColumns = CoinCopyOfArray(pinfo.originalColumns(),numberColumns);
1182	} else {
1183	OsiSolverInterface * solver = originalSolver->clone();
1184	clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
1185	lpSolver = clpSolver->getModelPtr();
1186	}
1187	// Tighten bounds
1188	lpSolver->tightenPrimalBounds(0.0,11,true);
1189	int numberColumns = clpSolver->getNumCols() ;
1190	double * saveColumnLower = CoinCopyOfArray(lpSolver->columnLower(),numberColumns);
1191	double * saveColumnUpper = CoinCopyOfArray(lpSolver->columnUpper(),numberColumns);
1192	//char generalPrint[200];
1193	const double *objective = lpSolver->getObjCoefficients() ;
1194	double *columnLower = lpSolver->columnLower() ;
1195	double *columnUpper = lpSolver->columnUpper() ;
1196	int numberRows = clpSolver->getNumRows();
1197	int iRow,iColumn;
1198
1199	// Row copy
1200	CoinPackedMatrix matrixByRow(*clpSolver->getMatrixByRow());
1201	const double * elementByRow = matrixByRow.getElements();
1202	const int * column = matrixByRow.getIndices();
1203	const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1204	const int * rowLength = matrixByRow.getVectorLengths();
1205
1206	// Column copy
1207	CoinPackedMatrix matrixByCol(*clpSolver->getMatrixByCol());
1208	//const double * element = matrixByCol.getElements();
1209	const int * row = matrixByCol.getIndices();
1210	const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
1211	const int * columnLength = matrixByCol.getVectorLengths();
1212
1213	const double * rowLower = clpSolver->getRowLower();
1214	const double * rowUpper = clpSolver->getRowUpper();
1215
1216	// Get maximum size of VUB tree
1217	// otherColumn is one fixed to 0 if this one zero
1218	int nEl = matrixByCol.getNumElements();
1219	CoinBigIndex * fixColumn = new CoinBigIndex [numberColumns+1];
1220	int * otherColumn = new int [nEl];
1221	int * fix = new int[numberColumns];
1222	char * mark = new char [numberColumns];
1223	memset(mark,0,numberColumns);
1224	int numberInteger=0;
1225	int numberOther=0;
1226	fixColumn[0]=0;
1227	double large = lpSolver->largeValue(); // treat bounds > this as infinite
1228	#ifndef NDEBUG
1229	double large2= 1.0e10*large;
1230	#endif
1231	double tolerance = lpSolver->primalTolerance();
1232	int * check = new int[numberRows];
1233	for (iRow=0;iRow<numberRows;iRow++) {
1234	check[iRow]=-2; // don't check
1235	if (rowLower[iRow]<-1.0e6&&rowUpper[iRow]>1.0e6)
1236	continue;// unlikely
1237	// possible row
1238	int numberPositive=0;
1239	int iPositive=-1;
1240	int numberNegative=0;
1241	int iNegative=-1;
1242	CoinBigIndex rStart = rowStart[iRow];
1243	CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
1244	CoinBigIndex j;
1245	int kColumn;
1246	for (j = rStart; j < rEnd; ++j) {
1247	double value=elementByRow[j];
1248	kColumn = column[j];
1249	if (columnUpper[kColumn]>columnLower[kColumn]) {
1250	if (value > 0.0) {
1251	numberPositive++;
1252	iPositive=kColumn;
1253	} else {
1254	numberNegative++;
1255	iNegative=kColumn;
1256	}
1257	}
1258	}
1259	if (numberPositive==1&&numberNegative==1)
1260	check[iRow]=-1; // try both
1261	if (numberPositive==1&&rowLower[iRow]>-1.0e20)
1262	check[iRow]=iPositive;
1263	else if (numberNegative==1&&rowUpper[iRow]<1.0e20)
1264	check[iRow]=iNegative;
1265	}
1266	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1267	fix[iColumn]=-1;
1268	if (columnUpper[iColumn] > columnLower[iColumn]+1.0e-8) {
1269	if (clpSolver->isInteger(iColumn))
1270	numberInteger++;
1271	if (columnLower[iColumn]==0.0) {
1272	bool infeasible=false;
1273	fix[iColumn]=0;
1274	// fake upper bound
1275	double saveUpper = columnUpper[iColumn];
1276	columnUpper[iColumn]=0.0;
1277	for (CoinBigIndex i=columnStart[iColumn];
1278	i<columnStart[iColumn]+columnLength[iColumn];i++) {
1279	iRow = row[i];
1280	if (check[iRow]!=-1&&check[iRow]!=iColumn)
1281	continue; // unlikely
1282	// possible row
1283	int infiniteUpper = 0;
1284	int infiniteLower = 0;
1285	double maximumUp = 0.0;
1286	double maximumDown = 0.0;
1287	double newBound;
1288	CoinBigIndex rStart = rowStart[iRow];
1289	CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
1290	CoinBigIndex j;
1291	int kColumn;
1292	// Compute possible lower and upper ranges
1293	for (j = rStart; j < rEnd; ++j) {
1294	double value=elementByRow[j];
1295	kColumn = column[j];
1296	if (value > 0.0) {
1297	if (columnUpper[kColumn] >= large) {
1298	++infiniteUpper;
1299	} else {
1300	maximumUp += columnUpper[kColumn] * value;
1301	}
1302	if (columnLower[kColumn] <= -large) {
1303	++infiniteLower;
1304	} else {
1305	maximumDown += columnLower[kColumn] * value;
1306	}
1307	} else if (value<0.0) {
1308	if (columnUpper[kColumn] >= large) {
1309	++infiniteLower;
1310	} else {
1311	maximumDown += columnUpper[kColumn] * value;
1312	}
1313	if (columnLower[kColumn] <= -large) {
1314	++infiniteUpper;
1315	} else {
1316	maximumUp += columnLower[kColumn] * value;
1317	}
1318	}
1319	}
1320	// Build in a margin of error
1321	maximumUp += 1.0e-8*fabs(maximumUp);
1322	maximumDown -= 1.0e-8*fabs(maximumDown);
1323	double maxUp = maximumUp+infiniteUpper*1.0e31;
1324	double maxDown = maximumDown-infiniteLower*1.0e31;
1325	if (maxUp <= rowUpper[iRow] + tolerance &&
1326	maxDown >= rowLower[iRow] - tolerance) {
1327	//printf("Redundant row in vubs %d\n",iRow);
1328	} else {
1329	if (maxUp < rowLower[iRow] -100.0*tolerance \|\|
1330	maxDown > rowUpper[iRow]+100.0*tolerance) {
1331	infeasible=true;
1332	break;
1333	}
1334	double lower = rowLower[iRow];
1335	double upper = rowUpper[iRow];
1336	for (j = rStart; j < rEnd; ++j) {
1337	double value=elementByRow[j];
1338	kColumn = column[j];
1339	double nowLower = columnLower[kColumn];
1340	double nowUpper = columnUpper[kColumn];
1341	if (value > 0.0) {
1342	// positive value
1343	if (lower>-large) {
1344	if (!infiniteUpper) {
1345	assert(nowUpper < large2);
1346	newBound = nowUpper +
1347	(lower - maximumUp) / value;
1348	// relax if original was large
1349	if (fabs(maximumUp)>1.0e8)
1350	newBound -= 1.0e-12*fabs(maximumUp);
1351	} else if (infiniteUpper==1&&nowUpper>large) {
1352	newBound = (lower -maximumUp) / value;
1353	// relax if original was large
1354	if (fabs(maximumUp)>1.0e8)
1355	newBound -= 1.0e-12*fabs(maximumUp);
1356	} else {
1357	newBound = -COIN_DBL_MAX;
1358	}
1359	if (newBound > nowLower + 1.0e-12&&newBound>-large) {
1360	// Tighten the lower bound
1361	// check infeasible (relaxed)
1362	if (nowUpper < newBound) {
1363	if (nowUpper - newBound <
1364	-100.0*tolerance) {
1365	infeasible=true;
1366	break;
1367	}
1368	}
1369	}
1370	}
1371	if (upper <large) {
1372	if (!infiniteLower) {
1373	assert(nowLower >- large2);
1374	newBound = nowLower +
1375	(upper - maximumDown) / value;
1376	// relax if original was large
1377	if (fabs(maximumDown)>1.0e8)
1378	newBound += 1.0e-12*fabs(maximumDown);
1379	} else if (infiniteLower==1&&nowLower<-large) {
1380	newBound = (upper - maximumDown) / value;
1381	// relax if original was large
1382	if (fabs(maximumDown)>1.0e8)
1383	newBound += 1.0e-12*fabs(maximumDown);
1384	} else {
1385	newBound = COIN_DBL_MAX;
1386	}
1387	if (newBound < nowUpper - 1.0e-12&&newBound<large) {
1388	// Tighten the upper bound
1389	// check infeasible (relaxed)
1390	if (nowLower > newBound) {
1391	if (newBound - nowLower <
1392	-100.0*tolerance) {
1393	infeasible=true;
1394	break;
1395	} else {
1396	newBound=nowLower;
1397	}
1398	}
1399	if (!newBound\|\|(clpSolver->isInteger(kColumn)&&newBound<0.999)) {
1400	// fix to zero
1401	if (!mark[kColumn]) {
1402	otherColumn[numberOther++]=kColumn;
1403	mark[kColumn]=1;
1404	if (check[iRow]==-1)
1405	check[iRow]=iColumn;
1406	else
1407	assert(check[iRow]==iColumn);
1408	}
1409	}
1410	}
1411	}
1412	} else {
1413	// negative value
1414	if (lower>-large) {
1415	if (!infiniteUpper) {
1416	assert(nowLower < large2);
1417	newBound = nowLower +
1418	(lower - maximumUp) / value;
1419	// relax if original was large
1420	if (fabs(maximumUp)>1.0e8)
1421	newBound += 1.0e-12*fabs(maximumUp);
1422	} else if (infiniteUpper==1&&nowLower<-large) {
1423	newBound = (lower -maximumUp) / value;
1424	// relax if original was large
1425	if (fabs(maximumUp)>1.0e8)
1426	newBound += 1.0e-12*fabs(maximumUp);
1427	} else {
1428	newBound = COIN_DBL_MAX;
1429	}
1430	if (newBound < nowUpper - 1.0e-12&&newBound<large) {
1431	// Tighten the upper bound
1432	// check infeasible (relaxed)
1433	if (nowLower > newBound) {
1434	if (newBound - nowLower <
1435	-100.0*tolerance) {
1436	infeasible=true;
1437	break;
1438	} else {
1439	newBound=nowLower;
1440	}
1441	}
1442	if (!newBound\|\|(clpSolver->isInteger(kColumn)&&newBound<0.999)) {
1443	// fix to zero
1444	if (!mark[kColumn]) {
1445	otherColumn[numberOther++]=kColumn;
1446	mark[kColumn]=1;
1447	if (check[iRow]==-1)
1448	check[iRow]=iColumn;
1449	else
1450	assert(check[iRow]==iColumn);
1451	}
1452	}
1453	}
1454	}
1455	if (upper <large) {
1456	if (!infiniteLower) {
1457	assert(nowUpper < large2);
1458	newBound = nowUpper +
1459	(upper - maximumDown) / value;
1460	// relax if original was large
1461	if (fabs(maximumDown)>1.0e8)
1462	newBound -= 1.0e-12*fabs(maximumDown);
1463	} else if (infiniteLower==1&&nowUpper>large) {
1464	newBound = (upper - maximumDown) / value;
1465	// relax if original was large
1466	if (fabs(maximumDown)>1.0e8)
1467	newBound -= 1.0e-12*fabs(maximumDown);
1468	} else {
1469	newBound = -COIN_DBL_MAX;
1470	}
1471	if (newBound > nowLower + 1.0e-12&&newBound>-large) {
1472	// Tighten the lower bound
1473	// check infeasible (relaxed)
1474	if (nowUpper < newBound) {
1475	if (nowUpper - newBound <
1476	-100.0*tolerance) {
1477	infeasible=true;
1478	break;
1479	}
1480	}
1481	}
1482	}
1483	}
1484	}
1485	}
1486	}
1487	for (int i=fixColumn[iColumn];i<numberOther;i++)
1488	mark[otherColumn[i]]=0;
1489	// reset bound unless infeasible
1490	if (!infeasible\|\|!clpSolver->isInteger(iColumn))
1491	columnUpper[iColumn]=saveUpper;
1492	else if (clpSolver->isInteger(iColumn))
1493	columnLower[iColumn]=1.0;
1494	}
1495	}
1496	fixColumn[iColumn+1]=numberOther;
1497	}
1498	delete [] check;
1499	delete [] mark;
1500	// Now do reverse way
1501	int * counts = new int [numberColumns];
1502	CoinZeroN(counts,numberColumns);
1503	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1504	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++)
1505	counts[otherColumn[i]]++;
1506	}
1507	numberOther=0;
1508	CoinBigIndex * fixColumn2 = new CoinBigIndex [numberColumns+1];
1509	int * otherColumn2 = new int [fixColumn[numberColumns]];
1510	fixColumn2[0]=0;
1511	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1512	numberOther += counts[iColumn];
1513	counts[iColumn]=0;
1514	fixColumn2[iColumn+1]=numberOther;
1515	}
1516	// Create other way
1517	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1518	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
1519	int jColumn=otherColumn[i];
1520	CoinBigIndex put = fixColumn2[jColumn] + counts[jColumn];
1521	counts[jColumn]++;
1522	otherColumn2[put]=iColumn;
1523	}
1524	}
1525	// get top layer i.e. those which are not fixed by any other
1526	int kLayer=0;
1527	while (true) {
1528	int numberLayered=0;
1529	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1530	if (fix[iColumn]==kLayer) {
1531	for (int i=fixColumn2[iColumn];i<fixColumn2[iColumn+1];i++) {
1532	int jColumn=otherColumn2[i];
1533	if (fix[jColumn]==kLayer) {
1534	fix[iColumn]=kLayer+100;
1535	}
1536	}
1537	}
1538	if (fix[iColumn]==kLayer) {
1539	numberLayered++;
1540	}
1541	}
1542	if (numberLayered) {
1543	kLayer+=100;
1544	} else {
1545	break;
1546	}
1547	}
1548	for (int iPass=0;iPass<2;iPass++) {
1549	for (int jLayer=0;jLayer<kLayer;jLayer++) {
1550	int check[]={-1,0,1,2,3,4,5,10,50,100,500,1000,5000,10000,COIN_INT_MAX};
1551	int nCheck = (int) (sizeof(check)/sizeof(int));
1552	int countsI[20];
1553	int countsC[20];
1554	assert (nCheck<=20);
1555	memset(countsI,0,nCheck*sizeof(int));
1556	memset(countsC,0,nCheck*sizeof(int));
1557	check[nCheck-1]=numberColumns;
1558	int numberLayered=0;
1559	int numberInteger=0;
1560	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1561	if (fix[iColumn]==jLayer) {
1562	numberLayered++;
1563	int nFix = fixColumn[iColumn+1]-fixColumn[iColumn];
1564	if (iPass) {
1565	// just integers
1566	nFix=0;
1567	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
1568	int jColumn=otherColumn[i];
1569	if (clpSolver->isInteger(jColumn))
1570	nFix++;
1571	}
1572	}
1573	int iFix;
1574	for (iFix=0;iFix<nCheck;iFix++) {
1575	if (nFix<=check[iFix])
1576	break;
1577	}
1578	assert (iFix<nCheck);
1579	if (clpSolver->isInteger(iColumn)) {
1580	numberInteger++;
1581	countsI[iFix]++;
1582	} else {
1583	countsC[iFix]++;
1584	}
1585	}
1586	}
1587	if (numberLayered) {
1588	printf("%d (%d integer) at priority %d\n",numberLayered,numberInteger,1+(jLayer/100));
1589	char buffer[50];
1590	for (int i=1;i<nCheck;i++) {
1591	if (countsI[i]\|\|countsC[i]) {
1592	if (i==1)
1593	sprintf(buffer," == zero ");
1594	else if (i<nCheck-1)
1595	sprintf(buffer,"> %6d and <= %6d ",check[i-1],check[i]);
1596	else
1597	sprintf(buffer,"> %6d ",check[i-1]);
1598	printf("%s %8d integers and %8d continuous\n",buffer,countsI[i],countsC[i]);
1599	}
1600	}
1601	}
1602	}
1603	}
1604	delete [] counts;
1605	// Now do fixing
1606	{
1607	// switch off presolve and up weight
1608	ClpSolve solveOptions;
1609	//solveOptions.setPresolveType(ClpSolve::presolveOff,0);
1610	solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
1611	//solveOptions.setSpecialOption(1,3,30); // sprint
1612	int numberColumns = lpSolver->numberColumns();
1613	int iColumn;
1614	bool allSlack=true;
1615	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1616	if (lpSolver->getColumnStatus(iColumn)==ClpSimplex::basic) {
1617	allSlack=false;
1618	break;
1619	}
1620	}
1621	if (allSlack)
1622	solveOptions.setSpecialOption(1,2,50); // idiot
1623	lpSolver->setInfeasibilityCost(1.0e11);
1624	lpSolver->defaultFactorizationFrequency();
1625	if (doAction!=11)
1626	lpSolver->initialSolve(solveOptions);
1627	double * columnLower = lpSolver->columnLower();
1628	double * columnUpper = lpSolver->columnUpper();
1629	double * fullSolution = lpSolver->primalColumnSolution();
1630	const double * dj = lpSolver->dualColumnSolution();
1631	int iPass=0;
1632	#define MAXPROB 2
1633	ClpSimplex models[MAXPROB];
1634	int pass[MAXPROB];
1635	int kPass=-1;
1636	int kLayer=0;
1637	int skipZero=0;
1638	if (skipZero2==-1)
1639	skipZero2=40; //-1;
1640	/* 0 fixed to 0 by choice
1641	1 lb of 1 by choice
1642	2 fixed to 0 by another
1643	3 as 2 but this go
1644	-1 free
1645	*/
1646	char * state = new char [numberColumns];
1647	for (iColumn=0;iColumn<numberColumns;iColumn++)
1648	state[iColumn]=-1;
1649	while (true) {
1650	double largest=-0.1;
1651	double smallest=1.1;
1652	int iLargest=-1;
1653	int iSmallest=-1;
1654	int atZero=0;
1655	int atOne=0;
1656	int toZero=0;
1657	int toOne=0;
1658	int numberFree=0;
1659	int numberGreater=0;
1660	columnLower = lpSolver->columnLower();
1661	columnUpper = lpSolver->columnUpper();
1662	fullSolution = lpSolver->primalColumnSolution();
1663	if (doAction==11) {
1664	{
1665	double * columnLower = lpSolver->columnLower();
1666	double * columnUpper = lpSolver->columnUpper();
1667	// lpSolver->dual();
1668	memcpy(columnLower,saveColumnLower,numberColumns*sizeof(double));
1669	memcpy(columnUpper,saveColumnUpper,numberColumns*sizeof(double));
1670	// lpSolver->dual();
1671	int iColumn;
1672	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1673	if (columnUpper[iColumn] > columnLower[iColumn]+1.0e-8) {
1674	if (clpSolver->isInteger(iColumn)) {
1675	double value = lastSolution[iColumn];
1676	int iValue = (int) (value+0.5);
1677	assert (fabs(value-((double) iValue))<1.0e-3);
1678	assert (iValue>=columnLower[iColumn]&&
1679	iValue<=columnUpper[iColumn]);
1680	columnLower[iColumn]=iValue;
1681	columnUpper[iColumn]=iValue;
1682	}
1683	}
1684	}
1685	lpSolver->initialSolve(solveOptions);
1686	memcpy(columnLower,saveColumnLower,numberColumns*sizeof(double));
1687	memcpy(columnUpper,saveColumnUpper,numberColumns*sizeof(double));
1688	}
1689	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1690	if (columnUpper[iColumn] > columnLower[iColumn]+1.0e-8) {
1691	if (clpSolver->isInteger(iColumn)) {
1692	double value = lastSolution[iColumn];
1693	int iValue = (int) (value+0.5);
1694	assert (fabs(value-((double) iValue))<1.0e-3);
1695	assert (iValue>=columnLower[iColumn]&&
1696	iValue<=columnUpper[iColumn]);
1697	if (!fix[iColumn]) {
1698	if (iValue==0) {
1699	state[iColumn]=0;
1700	assert (!columnLower[iColumn]);
1701	columnUpper[iColumn]=0.0;
1702	} else if (iValue==1) {
1703	state[iColumn]=1;
1704	columnLower[iColumn]=1.0;
1705	} else {
1706	// leave fixed
1707	columnLower[iColumn]=iValue;
1708	columnUpper[iColumn]=iValue;
1709	}
1710	} else if (iValue==0) {
1711	state[iColumn]=10;
1712	columnUpper[iColumn]=0.0;
1713	} else {
1714	// leave fixed
1715	columnLower[iColumn]=iValue;
1716	columnUpper[iColumn]=iValue;
1717	}
1718	}
1719	}
1720	}
1721	int jLayer=0;
1722	int nFixed=-1;
1723	int nTotalFixed=0;
1724	while (nFixed) {
1725	nFixed=0;
1726	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1727	if (columnUpper[iColumn]==0.0&&fix[iColumn]==jLayer) {
1728	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
1729	int jColumn=otherColumn[i];
1730	if (columnUpper[jColumn]) {
1731	bool canFix=true;
1732	for (int k=fixColumn2[jColumn];k<fixColumn2[jColumn+1];k++) {
1733	int kColumn=otherColumn2[k];
1734	if (state[kColumn]==1) {
1735	canFix=false;
1736	break;
1737	}
1738	}
1739	if (canFix) {
1740	columnUpper[jColumn]=0.0;
1741	nFixed++;
1742	}
1743	}
1744	}
1745	}
1746	}
1747	nTotalFixed += nFixed;
1748	jLayer += 100;
1749	}
1750	printf("This fixes %d variables in lower priorities\n",nTotalFixed);
1751	break;
1752	}
1753	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1754	if (!clpSolver->isInteger(iColumn)\|\|fix[iColumn]>kLayer)
1755	continue;
1756	// skip if fixes nothing
1757	if (fixColumn[iColumn+1]-fixColumn[iColumn]<=skipZero2)
1758	continue;
1759	double value = fullSolution[iColumn];
1760	if (value>1.00001) {
1761	numberGreater++;
1762	continue;
1763	}
1764	double lower = columnLower[iColumn];
1765	double upper = columnUpper[iColumn];
1766	if (lower==upper) {
1767	if (lower)
1768	atOne++;
1769	else
1770	atZero++;
1771	continue;
1772	}
1773	if (value<1.0e-7) {
1774	toZero++;
1775	columnUpper[iColumn]=0.0;
1776	state[iColumn]=10;
1777	continue;
1778	}
1779	if (value>1.0-1.0e-7) {
1780	toOne++;
1781	columnLower[iColumn]=1.0;
1782	state[iColumn]=1;
1783	continue;
1784	}
1785	numberFree++;
1786	// skip if fixes nothing
1787	if (fixColumn[iColumn+1]-fixColumn[iColumn]<=skipZero)
1788	continue;
1789	if (value<smallest) {
1790	smallest=value;
1791	iSmallest=iColumn;
1792	}
1793	if (value>largest) {
1794	largest=value;
1795	iLargest=iColumn;
1796	}
1797	}
1798	if (toZero\|\|toOne)
1799	printf("%d at 0 fixed and %d at one fixed\n",toZero,toOne);
1800	printf("%d variables free, %d fixed to 0, %d to 1 - smallest %g, largest %g\n",
1801	numberFree,atZero,atOne,smallest,largest);
1802	if (numberGreater&&!iPass)
1803	printf("%d variables have value > 1.0\n",numberGreater);
1804	//skipZero2=0; // leave 0 fixing
1805	int jLayer=0;
1806	int nFixed=-1;
1807	int nTotalFixed=0;
1808	while (nFixed) {
1809	nFixed=0;
1810	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1811	if (columnUpper[iColumn]==0.0&&fix[iColumn]==jLayer) {
1812	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
1813	int jColumn=otherColumn[i];
1814	if (columnUpper[jColumn]) {
1815	bool canFix=true;
1816	for (int k=fixColumn2[jColumn];k<fixColumn2[jColumn+1];k++) {
1817	int kColumn=otherColumn2[k];
1818	if (state[kColumn]==1) {
1819	canFix=false;
1820	break;
1821	}
1822	}
1823	if (canFix) {
1824	columnUpper[jColumn]=0.0;
1825	nFixed++;
1826	}
1827	}
1828	}
1829	}
1830	}
1831	nTotalFixed += nFixed;
1832	jLayer += 100;
1833	}
1834	printf("This fixes %d variables in lower priorities\n",nTotalFixed);
1835	if (iLargest<0)
1836	break;
1837	double movement;
1838	int way;
1839	if (smallest<=1.0-largest&&smallest<0.2) {
1840	columnUpper[iSmallest]=0.0;
1841	state[iSmallest]=0;
1842	movement=smallest;
1843	way=-1;
1844	} else {
1845	columnLower[iLargest]=1.0;
1846	state[iLargest]=1;
1847	movement=1.0-largest;
1848	way=1;
1849	}
1850	double saveObj = lpSolver->objectiveValue();
1851	iPass++;
1852	kPass = iPass%MAXPROB;
1853	models[kPass]=*lpSolver;
1854	if (way==-1) {
1855	// fix others
1856	for (int i=fixColumn[iSmallest];i<fixColumn[iSmallest+1];i++) {
1857	int jColumn=otherColumn[i];
1858	if (state[jColumn]==-1) {
1859	columnUpper[jColumn]=0.0;
1860	state[jColumn]=3;
1861	}
1862	}
1863	}
1864	pass[kPass]=iPass;
1865	double maxCostUp = COIN_DBL_MAX;
1866	objective = lpSolver->getObjCoefficients() ;
1867	if (way==-1)
1868	maxCostUp= (1.0-movement)*objective[iSmallest];
1869	lpSolver->setDualObjectiveLimit(saveObj+maxCostUp);
1870	crunchIt(lpSolver);
1871	double moveObj = lpSolver->objectiveValue()-saveObj;
1872	printf("movement %s was %g costing %g\n",
1873	(smallest<=1.0-largest) ? "down" : "up",movement,moveObj);
1874	if (way==-1&&(moveObj>=(1.0-movement)*objective[iSmallest]\|\|lpSolver->status())) {
1875	// go up
1876	columnLower = models[kPass].columnLower();
1877	columnUpper = models[kPass].columnUpper();
1878	columnLower[iSmallest]=1.0;
1879	columnUpper[iSmallest]=saveColumnUpper[iSmallest];
1880	*lpSolver=models[kPass];
1881	columnLower = lpSolver->columnLower();
1882	columnUpper = lpSolver->columnUpper();
1883	fullSolution = lpSolver->primalColumnSolution();
1884	dj = lpSolver->dualColumnSolution();
1885	columnLower[iSmallest]=1.0;
1886	columnUpper[iSmallest]=saveColumnUpper[iSmallest];
1887	state[iSmallest]=1;
1888	// unfix others
1889	for (int i=fixColumn[iSmallest];i<fixColumn[iSmallest+1];i++) {
1890	int jColumn=otherColumn[i];
1891	if (state[jColumn]==3) {
1892	columnUpper[jColumn]=saveColumnUpper[jColumn];
1893	state[jColumn]=-1;
1894	}
1895	}
1896	crunchIt(lpSolver);
1897	}
1898	models[kPass]=*lpSolver;
1899	}
1900	lpSolver->dual();
1901	printf("Fixing took %g seconds\n",CoinCpuTime()-time1);
1902	columnLower = lpSolver->columnLower();
1903	columnUpper = lpSolver->columnUpper();
1904	fullSolution = lpSolver->primalColumnSolution();
1905	dj = lpSolver->dualColumnSolution();
1906	int * sort = new int[numberColumns];
1907	double * dsort = new double[numberColumns];
1908	int chunk=20;
1909	int iRelax=0;
1910	//double fractionFixed=6.0/8.0;
1911	// relax while lots fixed
1912	while (true) {
1913	if (skipZero2>10&&doAction<10)
1914	break;
1915	iRelax++;
1916	int n=0;
1917	double sum0=0.0;
1918	double sum00=0.0;
1919	double sum1=0.0;
1920	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1921	if (!clpSolver->isInteger(iColumn)\|\|fix[iColumn]>kLayer)
1922	continue;
1923	// skip if fixes nothing
1924	if (fixColumn[iColumn+1]-fixColumn[iColumn]==0&&doAction<10)
1925	continue;
1926	double djValue = dj[iColumn];
1927	if (state[iColumn]==1) {
1928	assert (columnLower[iColumn]);
1929	assert (fullSolution[iColumn]>0.1);
1930	if (djValue>0.0) {
1931	//printf("YY dj of %d at %g is %g\n",iColumn,value,djValue);
1932	sum1 += djValue;
1933	sort[n]=iColumn;
1934	dsort[n++]=-djValue;
1935	} else {
1936	//printf("dj of %d at %g is %g\n",iColumn,value,djValue);
1937	}
1938	} else if (state[iColumn]==0\|\|state[iColumn]==10) {
1939	assert (fullSolution[iColumn]<0.1);
1940	assert (!columnUpper[iColumn]);
1941	double otherValue=0.0;
1942	int nn=0;
1943	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
1944	int jColumn=otherColumn[i];
1945	if (columnUpper[jColumn]==0.0) {
1946	if (dj[jColumn]<-1.0e-5) {
1947	nn++;
1948	otherValue += dj[jColumn]; // really need to look at rest
1949	}
1950	}
1951	}
1952	if (djValue<-1.0e-2\|\|otherValue<-1.0e-2) {
1953	//printf("XX dj of %d at %g is %g - %d out of %d contribute %g\n",iColumn,value,djValue,
1954	// nn,fixColumn[iColumn+1]-fixColumn[iColumn],otherValue);
1955	if (djValue<1.0e-8) {
1956	sum0 -= djValue;
1957	sum00 -= otherValue;
1958	sort[n]=iColumn;
1959	if (djValue<-1.0e-2)
1960	dsort[n++]=djValue+otherValue;
1961	else
1962	dsort[n++]=djValue+0.001*otherValue;
1963	}
1964	} else {
1965	//printf("dj of %d at %g is %g - no contribution from %d\n",iColumn,value,djValue,
1966	// fixColumn[iColumn+1]-fixColumn[iColumn]);
1967	}
1968	}
1969	}
1970	CoinSort_2(dsort,dsort+n,sort);
1971	double * originalColumnLower = saveColumnLower;
1972	double * originalColumnUpper = saveColumnUpper;
1973	double * lo = CoinCopyOfArray(columnLower,numberColumns);
1974	double * up = CoinCopyOfArray(columnUpper,numberColumns);
1975	for (int k=0;k<CoinMin(chunk,n);k++) {
1976	iColumn = sort[k];
1977	state[iColumn]=-2;
1978	}
1979	memcpy(columnLower,originalColumnLower,numberColumns*sizeof(double));
1980	memcpy(columnUpper,originalColumnUpper,numberColumns*sizeof(double));
1981	int nFixed=0;
1982	int nFixed0=0;
1983	int nFixed1=0;
1984	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1985	if (state[iColumn]==0\|\|state[iColumn]==10) {
1986	columnUpper[iColumn]=0.0;
1987	assert (lo[iColumn]==0.0);
1988	nFixed++;
1989	nFixed0++;
1990	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
1991	int jColumn=otherColumn[i];
1992	if (columnUpper[jColumn]) {
1993	bool canFix=true;
1994	for (int k=fixColumn2[jColumn];k<fixColumn2[jColumn+1];k++) {
1995	int kColumn=otherColumn2[k];
1996	if (state[kColumn]==1\|\|state[kColumn]==-2) {
1997	canFix=false;
1998	break;
1999	}
2000	}
2001	if (canFix) {
2002	columnUpper[jColumn]=0.0;
2003	assert (lo[jColumn]==0.0);
2004	nFixed++;
2005	}
2006	}
2007	}
2008	} else if (state[iColumn]==1) {
2009	columnLower[iColumn]=1.0;
2010	nFixed1++;
2011	}
2012	}
2013	printf("%d fixed %d orig 0 %d 1\n",nFixed,nFixed0,nFixed1);
2014	int jLayer=0;
2015	nFixed=-1;
2016	int nTotalFixed=0;
2017	while (nFixed) {
2018	nFixed=0;
2019	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
2020	if (columnUpper[iColumn]==0.0&&fix[iColumn]==jLayer) {
2021	for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
2022	int jColumn=otherColumn[i];
2023	if (columnUpper[jColumn]) {
2024	bool canFix=true;
2025	for (int k=fixColumn2[jColumn];k<fixColumn2[jColumn+1];k++) {
2026	int kColumn=otherColumn2[k];
2027	if (state[kColumn]==1\|\|state[kColumn]==-2) {
2028	canFix=false;
2029	break;
2030	}
2031	}
2032	if (canFix) {
2033	columnUpper[jColumn]=0.0;
2034	assert (lo[jColumn]==0.0);
2035	nFixed++;
2036	}
2037	}
2038	}
2039	}
2040	}
2041	nTotalFixed += nFixed;
2042	jLayer += 100;
2043	}
2044	nFixed=0;
2045	int nFixedI=0;
2046	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
2047	if (columnLower[iColumn]==columnUpper[iColumn]) {
2048	if (clpSolver->isInteger(iColumn))
2049	nFixedI++;
2050	nFixed++;
2051	}
2052	}
2053	printf("This fixes %d variables in lower priorities - total %d (%d integer) - all target %d, int target %d\n",
2054	nTotalFixed,nFixed,nFixedI,(int)(fractionFixednumberColumns),(int) (fractionIntFixednumberInteger));
2055	int nBad=0;
2056	int nRelax=0;
2057	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
2058	if (lo[iColumn]<columnLower[iColumn]\|\|
2059	up[iColumn]>columnUpper[iColumn]) {
2060	printf("bad %d old %g %g, new %g %g\n",iColumn,lo[iColumn],up[iColumn],
2061	columnLower[iColumn],columnUpper[iColumn]);
2062	nBad++;
2063	}
2064	if (lo[iColumn]>columnLower[iColumn]\|\|
2065	up[iColumn]<columnUpper[iColumn]) {
2066	nRelax++;
2067	}
2068	}
2069	printf("%d relaxed\n",nRelax);
2070	if (iRelax>20&&nRelax==chunk)
2071	nRelax=0;
2072	if (iRelax>50)
2073	nRelax=0;
2074	assert (!nBad);
2075	delete [] lo;
2076	delete [] up;
2077	lpSolver->primal(1);
2078	if (nFixed<fractionFixednumberColumns\|\|nFixedI<fractionIntFixednumberInteger\|\|!nRelax)
2079	break;
2080	}
2081	delete [] state;
2082	delete [] sort;
2083	delete [] dsort;
2084	}
2085	delete [] fix;
2086	delete [] fixColumn;
2087	delete [] otherColumn;
2088	delete [] otherColumn2;
2089	delete [] fixColumn2;
2090	// See if was presolved
2091	if (originalColumns) {
2092	columnLower = lpSolver->columnLower();
2093	columnUpper = lpSolver->columnUpper();
2094	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
2095	saveColumnLower[iColumn] = columnLower[iColumn];
2096	saveColumnUpper[iColumn] = columnUpper[iColumn];
2097	}
2098	pinfo.postsolve(true);
2099	columnLower = originalLpSolver->columnLower();
2100	columnUpper = originalLpSolver->columnUpper();
2101	double * newColumnLower = lpSolver->columnLower();
2102	double * newColumnUpper = lpSolver->columnUpper();
2103	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2104	int jColumn = originalColumns[iColumn];
2105	columnLower[jColumn] = CoinMax(columnLower[jColumn],newColumnLower[iColumn]);
2106	columnUpper[jColumn] = CoinMin(columnUpper[jColumn],newColumnUpper[iColumn]);
2107	}
2108	numberColumns = originalLpSolver->numberColumns();
2109	delete [] originalColumns;
2110	}
2111	delete [] saveColumnLower;
2112	delete [] saveColumnUpper;
2113	if (!originalColumns) {
2114	// Basis
2115	memcpy(originalLpSolver->statusArray(),lpSolver->statusArray(),numberRows+numberColumns);
2116	memcpy(originalLpSolver->primalColumnSolution(),lpSolver->primalColumnSolution(),numberColumns*sizeof(double));
2117	memcpy(originalLpSolver->primalRowSolution(),lpSolver->primalRowSolution(),numberRows*sizeof(double));
2118	// Fix in solver
2119	columnLower = lpSolver->columnLower();
2120	columnUpper = lpSolver->columnUpper();
2121	}
2122	double * originalColumnLower = originalLpSolver->columnLower();
2123	double * originalColumnUpper = originalLpSolver->columnUpper();
2124	// number fixed
2125	doAction=0;
2126	for ( iColumn=0;iColumn<numberColumns;iColumn++) {
2127	originalColumnLower[iColumn] = columnLower[iColumn];
2128	originalColumnUpper[iColumn] = columnUpper[iColumn];
2129	if (columnLower[iColumn]==columnUpper[iColumn])
2130	doAction++;
2131	}
2132	printf("%d fixed by vub preprocessing\n",doAction);
2133	if (originalColumns) {
2134	originalLpSolver->initialSolve();
2135	}
2136	delete clpSolver;
2137	return NULL;
2138	}
2139	#endif
2140	#ifdef COIN_HAS_LINK
2141	/* Returns OsiSolverInterface (User should delete)
2142	On entry numberKnapsack is maximum number of Total entries
2143	*/
2144	static OsiSolverInterface *
2145	expandKnapsack(CoinModel & model, int * whichColumn, int * knapsackStart,
2146	int * knapsackRow, int &numberKnapsack,
2147	CglStored & stored, int logLevel,
2148	int fixedPriority, int SOSPriority,CoinModel & tightenedModel)
2149	{
2150	int maxTotal = numberKnapsack;
2151	// load from coin model
2152	OsiSolverLink *si = new OsiSolverLink();
2153	OsiSolverInterface * finalModel=NULL;
2154	si->setDefaultMeshSize(0.001);
2155	// need some relative granularity
2156	si->setDefaultBound(100.0);
2157	si->setDefaultMeshSize(0.01);
2158	si->setDefaultBound(100000.0);
2159	si->setIntegerPriority(1000);
2160	si->setBiLinearPriority(10000);
2161	si->load(model,true,logLevel);
2162	// get priorities
2163	const int * priorities=model.priorities();
2164	int numberColumns = model.numberColumns();
2165	if (priorities) {
2166	OsiObject ** objects = si->objects();
2167	int numberObjects = si->numberObjects();
2168	for (int iObj = 0;iObj<numberObjects;iObj++) {
2169	int iColumn = objects[iObj]->columnNumber();
2170	if (iColumn>=0&&iColumn<numberColumns) {
2171	#ifndef NDEBUG
2172	OsiSimpleInteger * obj =
2173	dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
2174	#endif
2175	assert (obj);
2176	int iPriority = priorities[iColumn];
2177	if (iPriority>0)
2178	objects[iObj]->setPriority(iPriority);
2179	}
2180	}
2181	if (fixedPriority>0) {
2182	si->setFixedPriority(fixedPriority);
2183	}
2184	if (SOSPriority<0)
2185	SOSPriority=100000;
2186	}
2187	CoinModel coinModel=*si->coinModel();
2188	assert(coinModel.numberRows()>0);
2189	tightenedModel = coinModel;
2190	int numberRows = coinModel.numberRows();
2191	// Mark variables
2192	int * whichKnapsack = new int [numberColumns];
2193	int iRow,iColumn;
2194	for (iColumn=0;iColumn<numberColumns;iColumn++)
2195	whichKnapsack[iColumn]=-1;
2196	int kRow;
2197	bool badModel=false;
2198	// analyze
2199	if (logLevel>1) {
2200	for (iRow=0;iRow<numberRows;iRow++) {
2201	/* Just obvious one at first
2202	positive non unit coefficients
2203	all integer
2204	positive rowUpper
2205	for now - linear (but further down in code may use nonlinear)
2206	column bounds should be tight
2207	*/
2208	//double lower = coinModel.getRowLower(iRow);
2209	double upper = coinModel.getRowUpper(iRow);
2210	if (upper<1.0e10) {
2211	CoinModelLink triple=coinModel.firstInRow(iRow);
2212	bool possible=true;
2213	int n=0;
2214	int n1=0;
2215	while (triple.column()>=0) {
2216	int iColumn = triple.column();
2217	const char * el = coinModel.getElementAsString(iRow,iColumn);
2218	if (!strcmp("Numeric",el)) {
2219	if (coinModel.columnLower(iColumn)==coinModel.columnUpper(iColumn)) {
2220	triple=coinModel.next(triple);
2221	continue; // fixed
2222	}
2223	double value=coinModel.getElement(iRow,iColumn);
2224	if (value<0.0) {
2225	possible=false;
2226	} else {
2227	n++;
2228	if (value==1.0)
2229	n1++;
2230	if (coinModel.columnLower(iColumn)<0.0)
2231	possible=false;
2232	if (!coinModel.isInteger(iColumn))
2233	possible=false;
2234	if (whichKnapsack[iColumn]>=0)
2235	possible=false;
2236	}
2237	} else {
2238	possible=false; // non linear
2239	}
2240	triple=coinModel.next(triple);
2241	}
2242	if (n-n1>1&&possible) {
2243	double lower = coinModel.getRowLower(iRow);
2244	double upper = coinModel.getRowUpper(iRow);
2245	CoinModelLink triple=coinModel.firstInRow(iRow);
2246	while (triple.column()>=0) {
2247	int iColumn = triple.column();
2248	lower -= coinModel.columnLower(iColumn)*triple.value();
2249	upper -= coinModel.columnLower(iColumn)*triple.value();
2250	triple=coinModel.next(triple);
2251	}
2252	printf("%d is possible %g <=",iRow,lower);
2253	// print
2254	triple=coinModel.firstInRow(iRow);
2255	while (triple.column()>=0) {
2256	int iColumn = triple.column();
2257	if (coinModel.columnLower(iColumn)!=coinModel.columnUpper(iColumn))
2258	printf(" (%d,el %g up %g)",iColumn,triple.value(),
2259	coinModel.columnUpper(iColumn)-coinModel.columnLower(iColumn));
2260	triple=coinModel.next(triple);
2261	}
2262	printf(" <= %g\n",upper);
2263	}
2264	}
2265	}
2266	}
2267	numberKnapsack=0;
2268	for (kRow=0;kRow<numberRows;kRow++) {
2269	iRow=kRow;
2270	/* Just obvious one at first
2271	positive non unit coefficients
2272	all integer
2273	positive rowUpper
2274	for now - linear (but further down in code may use nonlinear)
2275	column bounds should be tight
2276	*/
2277	//double lower = coinModel.getRowLower(iRow);
2278	double upper = coinModel.getRowUpper(iRow);
2279	if (upper<1.0e10) {
2280	CoinModelLink triple=coinModel.firstInRow(iRow);
2281	bool possible=true;
2282	int n=0;
2283	int n1=0;
2284	while (triple.column()>=0) {
2285	int iColumn = triple.column();
2286	const char * el = coinModel.getElementAsString(iRow,iColumn);
2287	if (!strcmp("Numeric",el)) {
2288	if (coinModel.columnLower(iColumn)==coinModel.columnUpper(iColumn)) {
2289	triple=coinModel.next(triple);
2290	continue; // fixed
2291	}
2292	double value=coinModel.getElement(iRow,iColumn);
2293	if (value<0.0) {
2294	possible=false;
2295	} else {
2296	n++;
2297	if (value==1.0)
2298	n1++;
2299	if (coinModel.columnLower(iColumn)<0.0)
2300	possible=false;
2301	if (!coinModel.isInteger(iColumn))
2302	possible=false;
2303	if (whichKnapsack[iColumn]>=0)
2304	possible=false;
2305	}
2306	} else {
2307	possible=false; // non linear
2308	}
2309	triple=coinModel.next(triple);
2310	}
2311	if (n-n1>1&&possible) {
2312	// try
2313	CoinModelLink triple=coinModel.firstInRow(iRow);
2314	while (triple.column()>=0) {
2315	int iColumn = triple.column();
2316	if (coinModel.columnLower(iColumn)!=coinModel.columnUpper(iColumn))
2317	whichKnapsack[iColumn]=numberKnapsack;
2318	triple=coinModel.next(triple);
2319	}
2320	knapsackRow[numberKnapsack++]=iRow;
2321	}
2322	}
2323	}
2324	if (logLevel>0)
2325	printf("%d out of %d candidate rows are possible\n",numberKnapsack,numberRows);
2326	// Check whether we can get rid of nonlinearities
2327	/* mark rows
2328	-2 in knapsack and other variables
2329	-1 not involved
2330	n only in knapsack n
2331	*/
2332	int * markRow = new int [numberRows];
2333	for (iRow=0;iRow<numberRows;iRow++)
2334	markRow[iRow]=-1;
2335	int canDo=1; // OK and linear
2336	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2337	CoinModelLink triple=coinModel.firstInColumn(iColumn);
2338	int iKnapsack = whichKnapsack[iColumn];
2339	bool linear=true;
2340	// See if quadratic objective
2341	const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
2342	if (strcmp(expr,"Numeric")) {
2343	linear=false;
2344	}
2345	while (triple.row()>=0) {
2346	int iRow = triple.row();
2347	if (iKnapsack>=0) {
2348	if (markRow[iRow]==-1) {
2349	markRow[iRow]=iKnapsack;
2350	} else if (markRow[iRow]!=iKnapsack) {
2351	markRow[iRow]=-2;
2352	}
2353	}
2354	const char * expr = coinModel.getElementAsString(iRow,iColumn);
2355	if (strcmp(expr,"Numeric")) {
2356	linear=false;
2357	}
2358	triple=coinModel.next(triple);
2359	}
2360	if (!linear) {
2361	if (whichKnapsack[iColumn]<0) {
2362	canDo=0;
2363	break;
2364	} else {
2365	canDo=2;
2366	}
2367	}
2368	}
2369	int * markKnapsack = NULL;
2370	double * coefficient = NULL;
2371	double * linear = NULL;
2372	int * whichRow = NULL;
2373	int * lookupRow = NULL;
2374	badModel=(canDo==0);
2375	if (numberKnapsack&&canDo) {
2376	/* double check - OK if
2377	no nonlinear
2378	nonlinear only on columns in knapsack
2379	nonlinear only on columns in knapsack * ONE other - same for all in knapsack
2380	AND that is only row connected to knapsack
2381	(theoretically could split knapsack if two other and small numbers)
2382	also ONE could be ONE expression - not just a variable
2383	*/
2384	int iKnapsack;
2385	markKnapsack = new int [numberKnapsack];
2386	coefficient = new double [numberKnapsack];
2387	linear = new double [numberColumns];
2388	for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++)
2389	markKnapsack[iKnapsack]=-1;
2390	if (canDo==2) {
2391	for (iRow=-1;iRow<numberRows;iRow++) {
2392	int numberOdd;
2393	CoinPackedMatrix * row = coinModel.quadraticRow(iRow,linear,numberOdd);
2394	if (row) {
2395	// see if valid
2396	const double * element = row->getElements();
2397	const int * column = row->getIndices();
2398	const CoinBigIndex * columnStart = row->getVectorStarts();
2399	const int * columnLength = row->getVectorLengths();
2400	int numberLook = row->getNumCols();
2401	for (int i=0;i<numberLook;i++) {
2402	int iKnapsack=whichKnapsack[i];
2403	if (iKnapsack<0) {
2404	// might be able to swap - but for now can't have knapsack in
2405	for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
2406	int iColumn = column[j];
2407	if (whichKnapsack[iColumn]>=0) {
2408	canDo=0; // no good
2409	badModel=true;
2410	break;
2411	}
2412	}
2413	} else {
2414	// OK if in same knapsack - or maybe just one
2415	int marked=markKnapsack[iKnapsack];
2416	for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
2417	int iColumn = column[j];
2418	if (whichKnapsack[iColumn]!=iKnapsack&&whichKnapsack[iColumn]>=0) {
2419	canDo=0; // no good
2420	badModel=true;
2421	break;
2422	} else if (marked==-1) {
2423	markKnapsack[iKnapsack]=iColumn;
2424	marked=iColumn;
2425	coefficient[iKnapsack]=element[j];
2426	coinModel.associateElement(coinModel.columnName(iColumn),1.0);
2427	} else if (marked!=iColumn) {
2428	badModel=true;
2429	canDo=0; // no good
2430	break;
2431	} else {
2432	// could manage with different coefficients - but for now ...
2433	assert(coefficient[iKnapsack]==element[j]);
2434	}
2435	}
2436	}
2437	}
2438	delete row;
2439	}
2440	}
2441	}
2442	if (canDo) {
2443	// for any rows which are cuts
2444	whichRow = new int [numberRows];
2445	lookupRow = new int [numberRows];
2446	bool someNonlinear=false;
2447	double maxCoefficient=1.0;
2448	for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
2449	if (markKnapsack[iKnapsack]>=0) {
2450	someNonlinear=true;
2451	int iColumn = markKnapsack[iKnapsack];
2452	maxCoefficient = CoinMax(maxCoefficient,fabs(coefficient[iKnapsack]*coinModel.columnUpper(iColumn)));
2453	}
2454	}
2455	if (someNonlinear) {
2456	// associate all columns to stop possible error messages
2457	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2458	coinModel.associateElement(coinModel.columnName(iColumn),1.0);
2459	}
2460	}
2461	ClpSimplex tempModel;
2462	tempModel.loadProblem(coinModel);
2463	// Create final model - first without knapsacks
2464	int nCol=0;
2465	int nRow=0;
2466	for (iRow=0;iRow<numberRows;iRow++) {
2467	if (markRow[iRow]<0) {
2468	lookupRow[iRow]=nRow;
2469	whichRow[nRow++]=iRow;
2470	} else {
2471	lookupRow[iRow]=-1;
2472	}
2473	}
2474	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2475	if (whichKnapsack[iColumn]<0)
2476	whichColumn[nCol++]=iColumn;
2477	}
2478	ClpSimplex finalModelX(&tempModel,nRow,whichRow,nCol,whichColumn,false,false,false);
2479	OsiClpSolverInterface finalModelY(&finalModelX,true);
2480	finalModel = finalModelY.clone();
2481	finalModelY.releaseClp();
2482	// Put back priorities
2483	const int * priorities=model.priorities();
2484	if (priorities) {
2485	finalModel->findIntegers(false);
2486	OsiObject ** objects = finalModel->objects();
2487	int numberObjects = finalModel->numberObjects();
2488	for (int iObj = 0;iObj<numberObjects;iObj++) {
2489	int iColumn = objects[iObj]->columnNumber();
2490	if (iColumn>=0&&iColumn<nCol) {
2491	#ifndef NDEBUG
2492	OsiSimpleInteger * obj =
2493	dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
2494	#endif
2495	assert (obj);
2496	int iPriority = priorities[whichColumn[iColumn]];
2497	if (iPriority>0)
2498	objects[iObj]->setPriority(iPriority);
2499	}
2500	}
2501	}
2502	for (iRow=0;iRow<numberRows;iRow++) {
2503	whichRow[iRow]=iRow;
2504	}
2505	int numberOther=finalModel->getNumCols();
2506	int nLargest=0;
2507	int nelLargest=0;
2508	int nTotal=0;
2509	for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
2510	iRow = knapsackRow[iKnapsack];
2511	int nCreate = maxTotal;
2512	int nelCreate=coinModel.expandKnapsack(iRow,nCreate,NULL,NULL,NULL,NULL);
2513	if (nelCreate<0)
2514	badModel=true;
2515	nTotal+=nCreate;
2516	nLargest = CoinMax(nLargest,nCreate);
2517	nelLargest = CoinMax(nelLargest,nelCreate);
2518	}
2519	if (nTotal>maxTotal)
2520	badModel=true;
2521	if (!badModel) {
2522	// Now arrays for building
2523	nelLargest = CoinMax(nelLargest,nLargest)+1;
2524	double * buildObj = new double [nLargest];
2525	double * buildElement = new double [nelLargest];
2526	int * buildStart = new int[nLargest+1];
2527	int * buildRow = new int[nelLargest];
2528	// alow for integers in knapsacks
2529	OsiObject ** object = new OsiObject * [numberKnapsack+nTotal];
2530	int nSOS=0;
2531	int nObj=numberKnapsack;
2532	for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
2533	knapsackStart[iKnapsack]=finalModel->getNumCols();
2534	iRow = knapsackRow[iKnapsack];
2535	int nCreate = 10000;
2536	coinModel.expandKnapsack(iRow,nCreate,buildObj,buildStart,buildRow,buildElement);
2537	// Redo row numbers
2538	for (iColumn=0;iColumn<nCreate;iColumn++) {
2539	for (int j=buildStart[iColumn];j<buildStart[iColumn+1];j++) {
2540	int jRow=buildRow[j];
2541	jRow=lookupRow[jRow];
2542	assert (jRow>=0&&jRow<nRow);
2543	buildRow[j]=jRow;
2544	}
2545	}
2546	finalModel->addCols(nCreate,buildStart,buildRow,buildElement,NULL,NULL,buildObj);
2547	int numberFinal=finalModel->getNumCols();
2548	for (iColumn=numberOther;iColumn<numberFinal;iColumn++) {
2549	if (markKnapsack[iKnapsack]<0) {
2550	finalModel->setColUpper(iColumn,maxCoefficient);
2551	finalModel->setInteger(iColumn);
2552	} else {
2553	finalModel->setColUpper(iColumn,maxCoefficient+1.0);
2554	finalModel->setInteger(iColumn);
2555	}
2556	OsiSimpleInteger * sosObject = new OsiSimpleInteger(finalModel,iColumn);
2557	sosObject->setPriority(1000000);
2558	object[nObj++]=sosObject;
2559	buildRow[iColumn-numberOther]=iColumn;
2560	buildElement[iColumn-numberOther]=1.0;
2561	}
2562	if (markKnapsack[iKnapsack]<0) {
2563	// convexity row
2564	finalModel->addRow(numberFinal-numberOther,buildRow,buildElement,1.0,1.0);
2565	} else {
2566	int iColumn = markKnapsack[iKnapsack];
2567	int n=numberFinal-numberOther;
2568	buildRow[n]=iColumn;
2569	buildElement[n++]=-fabs(coefficient[iKnapsack]);
2570	// convexity row (sort of)
2571	finalModel->addRow(n,buildRow,buildElement,0.0,0.0);
2572	OsiSOS * sosObject = new OsiSOS(finalModel,n-1,buildRow,NULL,1);
2573	sosObject->setPriority(iKnapsack+SOSPriority);
2574	// Say not integral even if is (switch off heuristics)
2575	sosObject->setIntegerValued(false);
2576	object[nSOS++]=sosObject;
2577	}
2578	numberOther=numberFinal;
2579	}
2580	finalModel->addObjects(nObj,object);
2581	for (iKnapsack=0;iKnapsack<nObj;iKnapsack++)
2582	delete object[iKnapsack];
2583	delete [] object;
2584	// Can we move any rows to cuts
2585	const int * cutMarker = coinModel.cutMarker();
2586	if (cutMarker&&0) {
2587	printf("AMPL CUTS OFF until global cuts fixed\n");
2588	cutMarker=NULL;
2589	}
2590	if (cutMarker) {
2591	// Row copy
2592	const CoinPackedMatrix * matrixByRow = finalModel->getMatrixByRow();
2593	const double * elementByRow = matrixByRow->getElements();
2594	const int * column = matrixByRow->getIndices();
2595	const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2596	const int * rowLength = matrixByRow->getVectorLengths();
2597
2598	const double * rowLower = finalModel->getRowLower();
2599	const double * rowUpper = finalModel->getRowUpper();
2600	int nDelete=0;
2601	for (iRow=0;iRow<numberRows;iRow++) {
2602	if (cutMarker[iRow]&&lookupRow[iRow]>=0) {
2603	int jRow=lookupRow[iRow];
2604	whichRow[nDelete++]=jRow;
2605	int start = rowStart[jRow];
2606	stored.addCut(rowLower[jRow],rowUpper[jRow],
2607	rowLength[jRow],column+start,elementByRow+start);
2608	}
2609	}
2610	finalModel->deleteRows(nDelete,whichRow);
2611	}
2612	knapsackStart[numberKnapsack]=finalModel->getNumCols();
2613	delete [] buildObj;
2614	delete [] buildElement;
2615	delete [] buildStart;
2616	delete [] buildRow;
2617	finalModel->writeMps("full");
2618	}
2619	}
2620	}
2621	delete [] whichKnapsack;
2622	delete [] markRow;
2623	delete [] markKnapsack;
2624	delete [] coefficient;
2625	delete [] linear;
2626	delete [] whichRow;
2627	delete [] lookupRow;
2628	delete si;
2629	si=NULL;
2630	if (!badModel&&finalModel) {
2631	finalModel->setDblParam(OsiObjOffset,coinModel.objectiveOffset());
2632	return finalModel;
2633	} else {
2634	delete finalModel;
2635	printf("can't make knapsacks - did you set fixedPriority (extra1)\n");
2636	return NULL;
2637	}
2638	}
2639	// Fills in original solution (coinModel length)
2640	static void
2641	afterKnapsack(const CoinModel & coinModel2, const int * whichColumn, const int * knapsackStart,
2642	const int * knapsackRow, int numberKnapsack,
2643	const double * knapsackSolution, double * solution, int logLevel)
2644	{
2645	CoinModel coinModel = coinModel2;
2646	int numberColumns = coinModel.numberColumns();
2647	int iColumn;
2648	// associate all columns to stop possible error messages
2649	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2650	coinModel.associateElement(coinModel.columnName(iColumn),1.0);
2651	}
2652	CoinZeroN(solution,numberColumns);
2653	int nCol=knapsackStart[0];
2654	for (iColumn=0;iColumn<nCol;iColumn++) {
2655	int jColumn = whichColumn[iColumn];
2656	solution[jColumn]=knapsackSolution[iColumn];
2657	}
2658	int * buildRow = new int [numberColumns]; // wild overkill
2659	double * buildElement = new double [numberColumns];
2660	int iKnapsack;
2661	for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) {
2662	int k=-1;
2663	double value=0.0;
2664	for (iColumn=knapsackStart[iKnapsack];iColumn<knapsackStart[iKnapsack+1];iColumn++) {
2665	if (knapsackSolution[iColumn]>1.0e-5) {
2666	if (k>=0) {
2667	printf("Two nonzero values for knapsack %d at (%d,%g) and (%d,%g)\n",iKnapsack,
2668	k,knapsackSolution[k],iColumn,knapsackSolution[iColumn]);
2669	abort();
2670	}
2671	k=iColumn;
2672	value=floor(knapsackSolution[iColumn]+0.5);
2673	assert (fabs(value-knapsackSolution[iColumn])<1.0e-5);
2674	}
2675	}
2676	if (k>=0) {
2677	int iRow = knapsackRow[iKnapsack];
2678	int nCreate = 10000;
2679	int nel=coinModel.expandKnapsack(iRow,nCreate,NULL,NULL,buildRow,buildElement,k-knapsackStart[iKnapsack]);
2680	assert (nel);
2681	if (logLevel>0)
2682	printf("expanded column %d in knapsack %d has %d nonzero entries:\n",
2683	k-knapsackStart[iKnapsack],iKnapsack,nel);
2684	for (int i=0;i<nel;i++) {
2685	int jColumn = buildRow[i];
2686	double value = buildElement[i];
2687	if (logLevel>0)
2688	printf("%d - original %d has value %g\n",i,jColumn,value);
2689	solution[jColumn]=value;
2690	}
2691	}
2692	}
2693	delete [] buildRow;
2694	delete [] buildElement;
2695	#if 0
2696	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2697	if (solution[iColumn]>1.0e-5&&coinModel.isInteger(iColumn))
2698	printf("%d %g\n",iColumn,solution[iColumn]);
2699	}
2700	#endif
2701	}
2702	#endif
2703	#if 0
2704	static int outDupRow(OsiSolverInterface * solver)
2705	{
2706	CglDuplicateRow dupCuts(solver);
2707	CglTreeInfo info;
2708	info.level = 0;
2709	info.pass = 0;
2710	int numberRows = solver->getNumRows();
2711	info.formulation_rows = numberRows;
2712	info.inTree = false;
2713	info.strengthenRow= NULL;
2714	info.pass = 0;
2715	OsiCuts cs;
2716	dupCuts.generateCuts(*solver,cs,info);
2717	const int * duplicate = dupCuts.duplicate();
2718	// Get rid of duplicate rows
2719	int * which = new int[numberRows];
2720	int numberDrop=0;
2721	for (int iRow=0;iRow<numberRows;iRow++) {
2722	if (duplicate[iRow]==-2\|\|duplicate[iRow]>=0)
2723	which[numberDrop++]=iRow;
2724	}
2725	if (numberDrop) {
2726	solver->deleteRows(numberDrop,which);
2727	}
2728	delete [] which;
2729	// see if we have any column cuts
2730	int numberColumnCuts = cs.sizeColCuts() ;
2731	const double * columnLower = solver->getColLower();
2732	const double * columnUpper = solver->getColUpper();
2733	for (int k = 0;k<numberColumnCuts;k++) {
2734	OsiColCut * thisCut = cs.colCutPtr(k) ;
2735	const CoinPackedVector & lbs = thisCut->lbs() ;
2736	const CoinPackedVector & ubs = thisCut->ubs() ;
2737	int j ;
2738	int n ;
2739	const int * which ;
2740	const double * values ;
2741	n = lbs.getNumElements() ;
2742	which = lbs.getIndices() ;
2743	values = lbs.getElements() ;
2744	for (j = 0;j<n;j++) {
2745	int iColumn = which[j] ;
2746	if (values[j]>columnLower[iColumn])
2747	solver->setColLower(iColumn,values[j]) ;
2748	}
2749	n = ubs.getNumElements() ;
2750	which = ubs.getIndices() ;
2751	values = ubs.getElements() ;
2752	for (j = 0;j<n;j++) {
2753	int iColumn = which[j] ;
2754	if (values[j]<columnUpper[iColumn])
2755	solver->setColUpper(iColumn,values[j]) ;
2756	}
2757	}
2758	return numberDrop;
2759	}
2760	#endif
2761	void checkSOS(CbcModel * babModel, const OsiSolverInterface * solver)
2762	{
2763	#ifdef COIN_DEVELOP
2764	if (!babModel->ownObjects())
2765	return;
2766	//const double *objective = solver->getObjCoefficients() ;
2767	const double *columnLower = solver->getColLower() ;
2768	const double * columnUpper = solver->getColUpper() ;
2769	const double * solution = solver->getColSolution();
2770	//int numberColumns = solver->getNumCols() ;
2771	//int numberRows = solver->getNumRows();
2772	//double direction = solver->getObjSense();
2773	//int iRow,iColumn;
2774
2775	// Row copy
2776	CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
2777	//const double * elementByRow = matrixByRow.getElements();
2778	//const int * column = matrixByRow.getIndices();
2779	//const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2780	const int * rowLength = matrixByRow.getVectorLengths();
2781
2782	// Column copy
2783	CoinPackedMatrix matrixByCol(*solver->getMatrixByCol());
2784	const double * element = matrixByCol.getElements();
2785	const int * row = matrixByCol.getIndices();
2786	const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
2787	const int * columnLength = matrixByCol.getVectorLengths();
2788
2789	const double * rowLower = solver->getRowLower();
2790	const double * rowUpper = solver->getRowUpper();
2791	OsiObject ** objects = babModel->objects();
2792	int numberObjects = babModel->numberObjects();
2793	for (int iObj = 0;iObj<numberObjects;iObj++) {
2794	CbcSOS * objSOS =
2795	dynamic_cast <CbcSOS *>(objects[iObj]) ;
2796	if (objSOS) {
2797	int n=objSOS->numberMembers();
2798	const int * which = objSOS->members();
2799	const double * weight = objSOS->weights();
2800	int type = objSOS->sosType();
2801	// convexity row?
2802	int iColumn;
2803	iColumn=which[0];
2804	int j;
2805	int convex=-1;
2806	for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
2807	int iRow = row[j];
2808	double value = element[j];
2809	if (rowLower[iRow]==1.0&&rowUpper[iRow]==1.0&&
2810	value==1.0) {
2811	// possible
2812	if (rowLength[iRow]==n) {
2813	if (convex==-1)
2814	convex=iRow;
2815	else
2816	convex=-2;
2817	}
2818	}
2819	}
2820	printf ("set %d of type %d has %d members - possible convexity row %d\n",
2821	iObj,type,n,convex);
2822	for (int i=0;i<n;i++) {
2823	iColumn = which[i];
2824	int convex2=-1;
2825	for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
2826	int iRow = row[j];
2827	if (iRow==convex) {
2828	double value = element[j];
2829	if (value==1.0) {
2830	convex2=iRow;
2831	}
2832	}
2833	}
2834	if (convex2<0&&convex>=0) {
2835	printf("odd convexity row\n");
2836	convex=-2;
2837	}
2838	printf("col %d has weight %g and value %g, bounds %g %g\n",
2839	iColumn,weight[i],solution[iColumn],columnLower[iColumn],
2840	columnUpper[iColumn]);
2841	}
2842	}
2843	}
2844	#endif
2845	}
2846	#ifndef NEW_STYLE_SOLVER
2847	int callCbc1(const char * input2, CbcModel & model, int callBack(CbcModel * currentSolver, int whereFrom))
2848	{
2849	char * input = strdup(input2);
2850	int length = strlen(input);
2851	bool blank = input[0]=='0';
2852	int n=blank ? 0 : 1;
2853	for (int i=0;i<length;i++) {
2854	if (blank) {
2855	// look for next non blank
2856	if (input[i]==' ') {
2857	continue;
2858	} else {
2859	n++;
2860	blank=false;
2861	}
2862	} else {
2863	// look for next blank
2864	if (input[i]!=' ') {
2865	continue;
2866	} else {
2867	blank=true;
2868	}
2869	}
2870	}
2871	char ** argv = new char * [n+2];
2872	argv[0]=strdup("cbc");
2873	int i=0;
2874	while(input[i]==' ')
2875	i++;
2876	for (int j=0;j<n;j++) {
2877	int saveI=i;
2878	for (;i<length;i++) {
2879	// look for next blank
2880	if (input[i]!=' ') {
2881	continue;
2882	} else {
2883	break;
2884	}
2885	}
2886	input[i]='\0';
2887	argv[j+1]=strdup(input+saveI);
2888	while(input[i]==' ')
2889	i++;
2890	}
2891	argv[n+1]=strdup("-quit");
2892	free(input);
2893	totalTime=0.0;
2894	currentBranchModel = NULL;
2895	CbcOrClpRead_mode=1;
2896	CbcOrClpReadCommand=stdin;
2897	noPrinting=false;
2898	int returnCode = CbcMain1(n+2,const_cast<const char **>(argv),model,callBack);
2899	for (int k=0;k<n+2;k++)
2900	free(argv[k]);
2901	delete [] argv;
2902	return returnCode;
2903	}
2904	int callCbc1(const std::string input2, CbcModel & babSolver)
2905	{
2906	char * input3 = strdup(input2.c_str());
2907	int returnCode=callCbc1(input3,babSolver);
2908	free(input3);
2909	return returnCode;
2910	}
2911	int callCbc(const char * input2, CbcModel & babSolver)
2912	{
2913	CbcMain0(babSolver);
2914	return callCbc1(input2,babSolver);
2915	}
2916	int callCbc(const std::string input2, CbcModel & babSolver)
2917	{
2918	char * input3 = strdup(input2.c_str());
2919	CbcMain0(babSolver);
2920	int returnCode=callCbc1(input3,babSolver);
2921	free(input3);
2922	return returnCode;
2923	}
2924	int callCbc(const char * input2, OsiClpSolverInterface& solver1)
2925	{
2926	CbcModel model(solver1);
2927	return callCbc(input2,model);
2928	}
2929	int callCbc(const char * input2)
2930	{
2931	{
2932	OsiClpSolverInterface solver1;
2933	return callCbc(input2,solver1);
2934	}
2935	}
2936	int callCbc(const std::string input2, OsiClpSolverInterface& solver1)
2937	{
2938	char * input3 = strdup(input2.c_str());
2939	int returnCode=callCbc(input3,solver1);
2940	free(input3);
2941	return returnCode;
2942	}
2943	int callCbc(const std::string input2)
2944	{
2945	char * input3 = strdup(input2.c_str());
2946	OsiClpSolverInterface solver1;
2947	int returnCode=callCbc(input3,solver1);
2948	free(input3);
2949	return returnCode;
2950	}
2951	static int dummyCallBack(CbcModel * model, int whereFrom)
2952	{
2953	return 0;
2954	}
2955	int CbcMain1 (int argc, const char *argv[],
2956	CbcModel & model)
2957	{
2958	return CbcMain1(argc,argv,model,dummyCallBack);
2959	}
2960	int callCbc1(const std::string input2, CbcModel & babSolver, int callBack(CbcModel * currentSolver, int whereFrom))
2961	{
2962	char * input3 = strdup(input2.c_str());
2963	int returnCode=callCbc1(input3,babSolver,callBack);
2964	free(input3);
2965	return returnCode;
2966	}
2967	int callCbc1(const char * input2, CbcModel & model)
2968	{
2969	return callCbc1(input2,model,dummyCallBack);
2970	}
2971	int CbcMain (int argc, const char *argv[],
2972	CbcModel & model)
2973	{
2974	CbcMain0(model);
2975	return CbcMain1(argc,argv,model);
2976	}
2977	#define CBCMAXPARAMETERS 200
2978	static CbcOrClpParam parameters[CBCMAXPARAMETERS];
2979	static int numberParameters=0 ;
2980	void CbcMain0 (CbcModel & model)
2981	{
2982	OsiClpSolverInterface * originalSolver = dynamic_cast<OsiClpSolverInterface *> (model.solver());
2983	assert (originalSolver);
2984	CoinMessageHandler * generalMessageHandler = originalSolver->messageHandler();
2985	generalMessageHandler->setPrefix(true);
2986	OsiSolverInterface * solver = model.solver();
2987	OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
2988	ClpSimplex * lpSolver = clpSolver->getModelPtr();
2989	lpSolver->setPerturbation(50);
2990	lpSolver->messageHandler()->setPrefix(false);
2991	establishParams(numberParameters,parameters) ;
2992	const char dirsep = CoinFindDirSeparator();
2993	std::string directory;
2994	std::string dirSample;
2995	std::string dirNetlib;
2996	std::string dirMiplib;
2997	if (dirsep == '/') {
2998	directory = "./";
2999	dirSample = "../../Data/Sample/";
3000	dirNetlib = "../../Data/Netlib/";
3001	dirMiplib = "../../Data/miplib3/";
3002	} else {
3003	directory = ".\\";
3004	dirSample = "..\\..\\Data\\Sample\\";
3005	dirNetlib = "..\\..\\Data\\Netlib\\";
3006	dirMiplib = "..\\..\\Data\\miplib3\\";
3007	}
3008	std::string defaultDirectory = directory;
3009	std::string importFile ="";
3010	std::string exportFile ="default.mps";
3011	std::string importBasisFile ="";
3012	std::string importPriorityFile ="";
3013	std::string debugFile="";
3014	std::string printMask="";
3015	std::string exportBasisFile ="default.bas";
3016	std::string saveFile ="default.prob";
3017	std::string restoreFile ="default.prob";
3018	std::string solutionFile ="stdout";
3019	std::string solutionSaveFile ="solution.file";
3020	int doIdiot=-1;
3021	int outputFormat=2;
3022	int substitution=3;
3023	int dualize=3;
3024	int preSolve=5;
3025	int doSprint=-1;
3026	int testOsiParameters=-1;
3027	parameters[whichParam(BASISIN,numberParameters,parameters)].setStringValue(importBasisFile);
3028	parameters[whichParam(PRIORITYIN,numberParameters,parameters)].setStringValue(importPriorityFile);
3029	parameters[whichParam(BASISOUT,numberParameters,parameters)].setStringValue(exportBasisFile);
3030	parameters[whichParam(DEBUG,numberParameters,parameters)].setStringValue(debugFile);
3031	parameters[whichParam(PRINTMASK,numberParameters,parameters)].setStringValue(printMask);
3032	parameters[whichParam(DIRECTORY,numberParameters,parameters)].setStringValue(directory);
3033	parameters[whichParam(DIRSAMPLE,numberParameters,parameters)].setStringValue(dirSample);
3034	parameters[whichParam(DIRNETLIB,numberParameters,parameters)].setStringValue(dirNetlib);
3035	parameters[whichParam(DIRMIPLIB,numberParameters,parameters)].setStringValue(dirMiplib);
3036	parameters[whichParam(DUALBOUND,numberParameters,parameters)].setDoubleValue(lpSolver->dualBound());
3037	parameters[whichParam(DUALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->dualTolerance());
3038	parameters[whichParam(EXPORT,numberParameters,parameters)].setStringValue(exportFile);
3039	parameters[whichParam(IDIOT,numberParameters,parameters)].setIntValue(doIdiot);
3040	parameters[whichParam(IMPORT,numberParameters,parameters)].setStringValue(importFile);
3041	parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].setDoubleValue(1.0e-8);
3042	int slog = whichParam(SOLVERLOGLEVEL,numberParameters,parameters);
3043	int log = whichParam(LOGLEVEL,numberParameters,parameters);
3044	parameters[slog].setIntValue(1);
3045	clpSolver->messageHandler()->setLogLevel(1) ;
3046	model.messageHandler()->setLogLevel(1);
3047	lpSolver->setLogLevel(1);
3048	parameters[log].setIntValue(1);
3049	parameters[whichParam(MAXFACTOR,numberParameters,parameters)].setIntValue(lpSolver->factorizationFrequency());
3050	parameters[whichParam(MAXITERATION,numberParameters,parameters)].setIntValue(lpSolver->maximumIterations());
3051	parameters[whichParam(OUTPUTFORMAT,numberParameters,parameters)].setIntValue(outputFormat);
3052	parameters[whichParam(PRESOLVEPASS,numberParameters,parameters)].setIntValue(preSolve);
3053	parameters[whichParam(PERTVALUE,numberParameters,parameters)].setIntValue(lpSolver->perturbation());
3054	parameters[whichParam(PRIMALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->primalTolerance());
3055	parameters[whichParam(PRIMALWEIGHT,numberParameters,parameters)].setDoubleValue(lpSolver->infeasibilityCost());
3056	parameters[whichParam(RESTORE,numberParameters,parameters)].setStringValue(restoreFile);
3057	parameters[whichParam(SAVE,numberParameters,parameters)].setStringValue(saveFile);
3058	//parameters[whichParam(TIMELIMIT,numberParameters,parameters)].setDoubleValue(1.0e8);
3059	parameters[whichParam(TIMELIMIT_BAB,numberParameters,parameters)].setDoubleValue(1.0e8);
3060	parameters[whichParam(SOLUTION,numberParameters,parameters)].setStringValue(solutionFile);
3061	parameters[whichParam(SAVESOL,numberParameters,parameters)].setStringValue(solutionSaveFile);
3062	parameters[whichParam(SPRINT,numberParameters,parameters)].setIntValue(doSprint);
3063	parameters[whichParam(SUBSTITUTION,numberParameters,parameters)].setIntValue(substitution);
3064	parameters[whichParam(DUALIZE,numberParameters,parameters)].setIntValue(dualize);
3065	model.setNumberBeforeTrust(10);
3066	parameters[whichParam(NUMBERBEFORE,numberParameters,parameters)].setIntValue(5);
3067	parameters[whichParam(MAXNODES,numberParameters,parameters)].setIntValue(model.getMaximumNodes());
3068	model.setNumberStrong(5);
3069	parameters[whichParam(STRONGBRANCHING,numberParameters,parameters)].setIntValue(model.numberStrong());
3070	parameters[whichParam(INFEASIBILITYWEIGHT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcInfeasibilityWeight));
3071	parameters[whichParam(INTEGERTOLERANCE,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcIntegerTolerance));
3072	parameters[whichParam(INCREMENT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcCutoffIncrement));
3073	parameters[whichParam(TESTOSI,numberParameters,parameters)].setIntValue(testOsiParameters);
3074	parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].setIntValue(1003);
3075	#ifdef CBC_THREAD
3076	parameters[whichParam(THREADS,numberParameters,parameters)].setIntValue(0);
3077	#endif
3078	// Set up likely cut generators and defaults
3079	parameters[whichParam(PREPROCESS,numberParameters,parameters)].setCurrentOption("sos");
3080	parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].setIntValue(128\|64\|1);
3081	parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].setIntValue(1);
3082	parameters[whichParam(CUTPASSINTREE,numberParameters,parameters)].setIntValue(1);
3083	parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].setIntValue(-1);
3084	parameters[whichParam(MAXHOTITS,numberParameters,parameters)].setIntValue(100);
3085	parameters[whichParam(CUTSSTRATEGY,numberParameters,parameters)].setCurrentOption("on");
3086	parameters[whichParam(HEURISTICSTRATEGY,numberParameters,parameters)].setCurrentOption("on");
3087	parameters[whichParam(NODESTRATEGY,numberParameters,parameters)].setCurrentOption("fewest");
3088	parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
3089	parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
3090	parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
3091	parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption("off");
3092	parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption("ifmove");
3093	parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
3094	parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
3095	parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption("root");
3096	parameters[whichParam(LANDPCUTS,numberParameters,parameters)].setCurrentOption("off");
3097	parameters[whichParam(RESIDCUTS,numberParameters,parameters)].setCurrentOption("off");
3098	parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption("on");
3099	parameters[whichParam(FPUMP,numberParameters,parameters)].setCurrentOption("on");
3100	parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption("on");
3101	parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption("on");
3102	parameters[whichParam(RINS,numberParameters,parameters)].setCurrentOption("off");
3103	parameters[whichParam(LOCALTREE,numberParameters,parameters)].setCurrentOption("off");
3104	parameters[whichParam(COSTSTRATEGY,numberParameters,parameters)].setCurrentOption("off");
3105	}
3106	#endif
3107	/* 1 - add heuristics to model
3108	2 - do heuristics (and set cutoff and best solution)
3109	3 - for miplib test so skip some
3110	*/
3111	#ifndef NEW_STYLE_SOLVER
3112	static int doHeuristics(CbcModel * model,int type)
3113	#else
3114	int
3115	CbcSolver::doHeuristics(CbcModel * model,int type)
3116	#endif
3117	{
3118	#ifndef NEW_STYLE_SOLVER
3119	CbcOrClpParam * parameters_ = parameters;
3120	int numberParameters_ = numberParameters;
3121	#endif
3122	bool anyToDo=false;
3123	int logLevel = parameters_[whichParam(LOGLEVEL,numberParameters_,parameters_)].intValue();
3124	int useFpump = parameters_[whichParam(FPUMP,numberParameters_,parameters_)].currentOptionAsInteger();
3125	int useRounding = parameters_[whichParam(ROUNDING,numberParameters_,parameters_)].currentOptionAsInteger();
3126	int useGreedy = parameters_[whichParam(GREEDY,numberParameters_,parameters_)].currentOptionAsInteger();
3127	int useCombine = parameters_[whichParam(COMBINE,numberParameters_,parameters_)].currentOptionAsInteger();
3128	int useRINS = parameters_[whichParam(RINS,numberParameters_,parameters_)].currentOptionAsInteger();
3129	// FPump done first as it only works if no solution
3130	int kType = (type<3) ? type : 1;
3131	if (useFpump>=kType) {
3132	anyToDo=true;
3133	CbcHeuristicFPump heuristic4(*model);
3134	heuristic4.setFractionSmall(0.5);
3135	double dextra3 = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
3136	if (dextra3)
3137	heuristic4.setFractionSmall(dextra3);
3138	heuristic4.setMaximumPasses(parameters_[whichParam(FPUMPITS,numberParameters_,parameters_)].intValue());
3139	int pumpTune=parameters_[whichParam(FPUMPTUNE,numberParameters_,parameters_)].intValue();
3140	if (pumpTune>0) {
3141	/*
3142	>=10000000 for using obj
3143	>=1000000 use as accumulate switch
3144	>=1000 use index+1 as number of large loops
3145	>=100 use 0.05 objvalue as increment
3146	>=10 use +0.1 objvalue for cutoff (add)
3147	1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
3148	4 and static continuous, 5 as 3 but no internal integers
3149	6 as 3 but all slack basis!
3150	*/
3151	double value = model->solver()->getObjSense()*model->solver()->getObjValue();
3152	int w = pumpTune/10;
3153	int c = w % 10;
3154	w /= 10;
3155	int i = w % 10;
3156	w /= 10;
3157	int r = w;
3158	int accumulate = r/1000;
3159	r -= 1000*accumulate;
3160	if (accumulate>=10) {
3161	int which = accumulate/10;
3162	accumulate -= 10*which;
3163	which--;
3164	// weights and factors
3165	double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
3166	double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
3167	heuristic4.setInitialWeight(weight[which]);
3168	heuristic4.setWeightFactor(factor[which]);
3169	}
3170	// fake cutoff
3171	if (logLevel>1)
3172	printf("Setting ");
3173	if (c) {
3174	double cutoff;
3175	model->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
3176	cutoff = CoinMin(cutoff,value + 0.1fabs(value)c);
3177	double dextra1 = parameters_[whichParam(DEXTRA1,numberParameters_,parameters_)].doubleValue();
3178	if (dextra1)
3179	cutoff=dextra1;
3180	heuristic4.setFakeCutoff(cutoff);
3181	if (logLevel>1)
3182	printf("fake cutoff of %g ",cutoff);
3183	}
3184	if (i\|\|r) {
3185	// also set increment
3186	//double increment = (0.01i+0.005)(fabs(value)+1.0e-12);
3187	double increment = 0.0;
3188	double dextra2 = parameters_[whichParam(DEXTRA2,numberParameters_,parameters_)].doubleValue();
3189	if (dextra2)
3190	increment = dextra2;
3191	heuristic4.setAbsoluteIncrement(increment);
3192	heuristic4.setAccumulate(accumulate);
3193	heuristic4.setMaximumRetries(r+1);
3194	if (logLevel>1) {
3195	if (i)
3196	printf("increment of %g ",heuristic4.absoluteIncrement());
3197	if (accumulate)
3198	printf("accumulate of %d ",accumulate);
3199	printf("%d retries ",r+2);
3200	}
3201	}
3202	pumpTune = pumpTune%100;
3203	if (logLevel>1)
3204	printf("and setting when to %d\n",pumpTune+10);
3205	if (pumpTune==6)
3206	pumpTune =13;
3207	heuristic4.setWhen(pumpTune+10);
3208	}
3209	heuristic4.setHeuristicName("feasibility pump");
3210	//#define ROLF
3211	#ifdef ROLF
3212	CbcHeuristicFPump pump(*model);
3213	pump.setMaximumTime(60);
3214	pump.setMaximumPasses(100);
3215	pump.setMaximumRetries(1);
3216	pump.setFixOnReducedCosts(0);
3217	pump.setHeuristicName("Feasibility pump");
3218	pump.setFractionSmall(1.0);
3219	pump.setWhen(13);
3220	model->addHeuristic(&pump);
3221	#else
3222	model->addHeuristic(&heuristic4);
3223	#endif
3224	}
3225	if (useRounding>=type) {
3226	CbcRounding heuristic1(*model);
3227	heuristic1.setHeuristicName("rounding");
3228	model->addHeuristic(&heuristic1) ;
3229	anyToDo=true;
3230	}
3231	if (useCombine>=type) {
3232	CbcHeuristicLocal heuristic2(*model);
3233	heuristic2.setHeuristicName("combine solutions");
3234	heuristic2.setFractionSmall(0.6);
3235	heuristic2.setSearchType(1);
3236	model->addHeuristic(&heuristic2);
3237	anyToDo=true;
3238	}
3239	if (useGreedy>=type) {
3240	CbcHeuristicGreedyCover heuristic3(*model);
3241	heuristic3.setHeuristicName("greedy cover");
3242	CbcHeuristicGreedyEquality heuristic3a(*model);
3243	heuristic3a.setHeuristicName("greedy equality");
3244	model->addHeuristic(&heuristic3);
3245	model->addHeuristic(&heuristic3a);
3246	anyToDo=true;
3247	}
3248	if (useRINS>=kType) {
3249	CbcHeuristicRINS heuristic5(*model);
3250	heuristic5.setHeuristicName("RINS");
3251	heuristic5.setFractionSmall(0.6);
3252	model->addHeuristic(&heuristic5) ;
3253	anyToDo=true;
3254	}
3255	if (type==2&&anyToDo) {
3256	// Do heuristics
3257	#if 1
3258	// clean copy
3259	CbcModel model2(*model);
3260	// But get rid of heuristics in model
3261	model->doHeuristicsAtRoot(2);
3262	if (logLevel<=1)
3263	model2.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
3264	OsiBabSolver defaultC;
3265	//solver_->setAuxiliaryInfo(&defaultC);
3266	model2.passInSolverCharacteristics(&defaultC);
3267	// Save bounds
3268	int numberColumns = model2.solver()->getNumCols();
3269	model2.createContinuousSolver();
3270	bool cleanModel = !model2.numberIntegers()&&!model2.numberObjects();
3271	model2.findIntegers(false);
3272	model2.doHeuristicsAtRoot(1);
3273	if (cleanModel)
3274	model2.zapIntegerInformation(false);
3275	if (model2.bestSolution()) {
3276	double value = model2.getMinimizationObjValue();
3277	model->setCutoff(value);
3278	model->setBestSolution(model2.bestSolution(),numberColumns,value);
3279	model->setSolutionCount(1);
3280	model->setNumberHeuristicSolutions(1);
3281	}
3282	#else
3283	if (logLevel<=1)
3284	model->solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
3285	OsiBabSolver defaultC;
3286	//solver_->setAuxiliaryInfo(&defaultC);
3287	model->passInSolverCharacteristics(&defaultC);
3288	// Save bounds
3289	int numberColumns = model->solver()->getNumCols();
3290	model->createContinuousSolver();
3291	bool cleanModel = !model->numberIntegers()&&!model->numberObjects();
3292	model->findIntegers(false);
3293	model->doHeuristicsAtRoot(1);
3294	if (cleanModel)
3295	model->zapIntegerInformation(false);
3296	#endif
3297	return 0;
3298	} else {
3299	return 0;
3300	}
3301	}
3302
3303	int CbcClpUnitTest (const CbcModel & saveModel,
3304	std::string& dirMiplib, bool unitTestOnly);
3305	#ifdef NEW_STYLE_SOLVER
3306	/* This takes a list of commands, does "stuff" and returns
3307	returnMode -
3308	0 model and solver untouched - babModel updated
3309	1 model updated - just with solution basis etc
3310	2 model updated i.e. as babModel (babModel NULL) (only use without preprocessing)
3311	*/
3312	int
3313	CbcSolver::solve(const char * input2, int returnMode)
3314	{
3315	char * input = strdup(input2);
3316	int length = strlen(input);
3317	bool blank = input[0]=='0';
3318	int n=blank ? 0 : 1;
3319	for (int i=0;i<length;i++) {
3320	if (blank) {
3321	// look for next non blank
3322	if (input[i]==' ') {
3323	continue;
3324	} else {
3325	n++;
3326	blank=false;
3327	}
3328	} else {
3329	// look for next blank
3330	if (input[i]!=' ') {
3331	continue;
3332	} else {
3333	blank=true;
3334	}
3335	}
3336	}
3337	char ** argv = new char * [n+2];
3338	argv[0]=strdup("cbc");
3339	int i=0;
3340	while(input[i]==' ')
3341	i++;
3342	for (int j=0;j<n;j++) {
3343	int saveI=i;
3344	for (;i<length;i++) {
3345	// look for next blank
3346	if (input[i]!=' ') {
3347	continue;
3348	} else {
3349	break;
3350	}
3351	}
3352	char save = input[i];
3353	input[i]='\0';
3354	argv[j+1]=strdup(input+saveI);
3355	input[i]=save;
3356	while(input[i]==' ')
3357	i++;
3358	}
3359	argv[n+1]=strdup("-quit");
3360	free(input);
3361	int returnCode = solve(n+2,const_cast<const char **>(argv),returnMode);
3362	for (int k=0;k<n+2;k++)
3363	free(argv[k]);
3364	delete [] argv;
3365	return returnCode;
3366	}
3367	#endif
3368	/* Meaning of whereFrom:
3369	1 after initial solve by dualsimplex etc
3370	2 after preprocessing
3371	3 just before branchAndBound (so user can override)
3372	4 just after branchAndBound (before postprocessing)
3373	5 after postprocessing
3374	6 after a user called heuristic phase
3375	*/
3376
3377	#ifndef NEW_STYLE_SOLVER
3378	int CbcMain1 (int argc, const char *argv[],
3379	CbcModel & model,
3380	int callBack(CbcModel * currentSolver, int whereFrom))
3381	#else
3382	/* This takes a list of commands, does "stuff" and returns
3383	returnMode -
3384	0 model and solver untouched - babModel updated
3385	1 model updated - just with solution basis etc
3386	2 model updated i.e. as babModel (babModel NULL)
3387	*/
3388	int
3389	CbcSolver::solve (int argc, const char *argv[], int returnMode)
3390	#endif
3391	{
3392	#ifndef NEW_STYLE_SOLVER
3393	CbcOrClpParam * parameters_ = parameters;
3394	int numberParameters_ = numberParameters;
3395	CbcModel & model_ = model;
3396	CbcModel * babModel_ = NULL;
3397	int returnMode=1;
3398	int statusUserFunction_[1];
3399	int numberUserFunctions_=1; // to allow for ampl
3400	#else
3401	delete babModel_;
3402	babModel_ = NULL;
3403	CbcOrClpRead_mode=1;
3404	delete [] statusUserFunction_;
3405	statusUserFunction_ = new int [numberUserFunctions_];
3406	int iUser;
3407	#endif
3408	memset(statusUserFunction_,0,numberUserFunctions_*sizeof(int));
3409	/* Note
3410	This is meant as a stand-alone executable to do as much of coin as possible.
3411	It should only have one solver known to it.
3412	*/
3413	OsiClpSolverInterface * originalSolver = dynamic_cast<OsiClpSolverInterface *> (model_.solver());
3414	assert (originalSolver);
3415	CoinMessageHandler * generalMessageHandler = model_.messageHandler();
3416	generalMessageHandler->setPrefix(false);
3417	// Move handler across if not default
3418	if (!originalSolver->defaultHandler()&&originalSolver->getModelPtr()->defaultHandler())
3419	originalSolver->getModelPtr()->passInMessageHandler(originalSolver->messageHandler());
3420	CoinMessages generalMessages = originalSolver->getModelPtr()->messages();
3421	char generalPrint[10000];
3422	if (originalSolver->getModelPtr()->logLevel()==0)
3423	noPrinting=true;
3424	#ifndef NEW_STYLE_SOLVER
3425	bool noPrinting_=noPrinting;
3426	#endif
3427	// see if log in list
3428	for (int i=1;i<argc;i++) {
3429	if (!strncmp(argv[i],"log",3)) {
3430	const char * equals = strchr(argv[i],'=');
3431	if (equals&&atoi(equals+1)>0)
3432	noPrinting_=false;
3433	else
3434	noPrinting_=true;
3435	break;
3436	} else if (!strncmp(argv[i],"-log",4)&&i<argc-1) {
3437	if (atoi(argv[i+1])>0)
3438	noPrinting_=false;
3439	else
3440	noPrinting_=true;
3441	break;
3442	}
3443	}
3444	double time0;
3445	{
3446	double time1 = CoinCpuTime(),time2;
3447	time0=time1;
3448	bool goodModel=(originalSolver->getNumCols()) ? true : false;
3449
3450	CoinSighandler_t saveSignal=SIG_DFL;
3451	// register signal handler
3452	saveSignal = signal(SIGINT,signal_handler);
3453	// Set up all non-standard stuff
3454	int cutPass=-1234567;
3455	int cutPassInTree=-1234567;
3456	int tunePreProcess=5;
3457	int testOsiParameters=-1;
3458	// 0 normal, 1 from ampl or MIQP etc (2 allows cuts)
3459	int complicatedInteger=0;
3460	OsiSolverInterface * solver = model_.solver();
3461	OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3462	ClpSimplex * lpSolver = clpSolver->getModelPtr();
3463	if (noPrinting_) {
3464	setCbcOrClpPrinting(false);
3465	lpSolver->setLogLevel(0);
3466	}
3467	// For priorities etc
3468	int * priorities=NULL;
3469	int * branchDirection=NULL;
3470	double * pseudoDown=NULL;
3471	double * pseudoUp=NULL;
3472	double * solutionIn = NULL;
3473	int * prioritiesIn = NULL;
3474	int numberSOS = 0;
3475	int * sosStart = NULL;
3476	int * sosIndices = NULL;
3477	char * sosType = NULL;
3478	double * sosReference = NULL;
3479	int * cut=NULL;
3480	int * sosPriority=NULL;
3481	CglStored storedAmpl;
3482	CoinModel * coinModel = NULL;
3483	CoinModel saveCoinModel;
3484	CoinModel saveTightenedModel;
3485	int * whichColumn = NULL;
3486	int * knapsackStart=NULL;
3487	int * knapsackRow=NULL;
3488	int numberKnapsack=0;
3489	#ifdef NEW_STYLE_SOLVER
3490	int numberInputs=0;
3491	readMode_=CbcOrClpRead_mode;
3492	for (iUser=0;iUser<numberUserFunctions_;iUser++) {
3493	int status = userFunction_[iUser]->importData(this,argc,const_cast<char **>(argv));
3494	if (status>=0) {
3495	if (!status) {
3496	numberInputs++;
3497	statusUserFunction_[iUser]=1;
3498	goodModel=true;
3499	solver = model_.solver();
3500	clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3501	lpSolver = clpSolver->getModelPtr();
3502	} else {
3503	printf("Bad input from user function %s\n",userFunction_[iUser]->name().c_str());
3504	abort();
3505	}
3506	}
3507	}
3508	if (numberInputs>1) {
3509	printf("Two or more user inputs!\n");
3510	abort();
3511	}
3512	if (originalCoinModel_)
3513	complicatedInteger=1;
3514	testOsiParameters = intValue(TESTOSI);
3515	if (noPrinting_) {
3516	model_.messageHandler()->setLogLevel(0);
3517	setCbcOrClpPrinting(false);
3518	}
3519	CbcOrClpRead_mode=readMode_;
3520	#endif
3521	#ifdef COIN_HAS_ASL
3522	ampl_info info;
3523	{
3524	memset(&info,0,sizeof(info));
3525	if (argc>2&&!strcmp(argv[2],"-AMPL")) {
3526	statusUserFunction_[0]=1;
3527	// see if log in list
3528	noPrinting_=true;
3529	for (int i=1;i<argc;i++) {
3530	if (!strncmp(argv[i],"log",3)) {
3531	const char * equals = strchr(argv[i],'=');
3532	if (equals&&atoi(equals+1)>0) {
3533	noPrinting_=false;
3534	info.logLevel=atoi(equals+1);
3535	int log = whichParam(LOGLEVEL,numberParameters_,parameters_);
3536	parameters_[log].setIntValue(info.logLevel);
3537	// mark so won't be overWritten
3538	info.numberRows=-1234567;
3539	break;
3540	}
3541	}
3542	}
3543
3544	union { void * voidModel; CoinModel * model; } coinModelStart;
3545	coinModelStart.model=NULL;
3546	int returnCode = readAmpl(&info,argc, const_cast<char **>(argv),& coinModelStart.voidModel);
3547	coinModel=coinModelStart.model;
3548	if (returnCode)
3549	return returnCode;
3550	CbcOrClpRead_mode=2; // so will start with parameters
3551	// see if log in list (including environment)
3552	for (int i=1;i<info.numberArguments;i++) {
3553	if (!strcmp(info.arguments[i],"log")) {
3554	if (i<info.numberArguments-1&&atoi(info.arguments[i+1])>0)
3555	noPrinting_=false;
3556	break;
3557	}
3558	}
3559	if (noPrinting_) {
3560	model_.messageHandler()->setLogLevel(0);
3561	setCbcOrClpPrinting(false);
3562	}
3563	if (!noPrinting_)
3564	printf("%d rows, %d columns and %d elements\n",
3565	info.numberRows,info.numberColumns,info.numberElements);
3566	#ifdef COIN_HAS_LINK
3567	if (!coinModel) {
3568	#endif
3569	solver->loadProblem(info.numberColumns,info.numberRows,info.starts,
3570	info.rows,info.elements,
3571	info.columnLower,info.columnUpper,info.objective,
3572	info.rowLower,info.rowUpper);
3573	// take off cuts if ampl wants that
3574	if (info.cut&&0) {
3575	printf("AMPL CUTS OFF until global cuts fixed\n");
3576	info.cut=NULL;
3577	}
3578	if (info.cut) {
3579	int numberRows = info.numberRows;
3580	int * whichRow = new int [numberRows];
3581	// Row copy
3582	const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
3583	const double * elementByRow = matrixByRow->getElements();
3584	const int * column = matrixByRow->getIndices();
3585	const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
3586	const int * rowLength = matrixByRow->getVectorLengths();
3587
3588	const double * rowLower = solver->getRowLower();
3589	const double * rowUpper = solver->getRowUpper();
3590	int nDelete=0;
3591	for (int iRow=0;iRow<numberRows;iRow++) {
3592	if (info.cut[iRow]) {
3593	whichRow[nDelete++]=iRow;
3594	int start = rowStart[iRow];
3595	storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
3596	rowLength[iRow],column+start,elementByRow+start);
3597	}
3598	}
3599	solver->deleteRows(nDelete,whichRow);
3600	delete [] whichRow;
3601	}
3602	#ifdef COIN_HAS_LINK
3603	} else {
3604	// save
3605	saveCoinModel = *coinModel;
3606	// load from coin model
3607	OsiSolverLink solver1;
3608	OsiSolverInterface * solver2 = solver1.clone();
3609	model_.assignSolver(solver2,false);
3610	OsiSolverLink * si =
3611	dynamic_cast<OsiSolverLink *>(model_.solver()) ;
3612	assert (si != NULL);
3613	si->setDefaultMeshSize(0.001);
3614	// need some relative granularity
3615	si->setDefaultBound(100.0);
3616	double dextra3 = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
3617	if (dextra3)
3618	si->setDefaultMeshSize(dextra3);
3619	si->setDefaultBound(100000.0);
3620	si->setIntegerPriority(1000);
3621	si->setBiLinearPriority(10000);
3622	CoinModel * model2 = (CoinModel *) coinModel;
3623	int logLevel = parameters_[whichParam(LOGLEVEL,numberParameters_,parameters_)].intValue();
3624	si->load(*model2,true,logLevel);
3625	// redo
3626	solver = model_.solver();
3627	clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3628	lpSolver = clpSolver->getModelPtr();
3629	clpSolver->messageHandler()->setLogLevel(0) ;
3630	testOsiParameters=0;
3631	parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].setIntValue(0);
3632	complicatedInteger=1;
3633	if (info.cut) {
3634	printf("Sorry - can't do cuts with LOS as ruins delicate row order\n");
3635	abort();
3636	int numberRows = info.numberRows;
3637	int * whichRow = new int [numberRows];
3638	// Row copy
3639	const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
3640	const double * elementByRow = matrixByRow->getElements();
3641	const int * column = matrixByRow->getIndices();
3642	const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
3643	const int * rowLength = matrixByRow->getVectorLengths();
3644
3645	const double * rowLower = solver->getRowLower();
3646	const double * rowUpper = solver->getRowUpper();
3647	int nDelete=0;
3648	for (int iRow=0;iRow<numberRows;iRow++) {
3649	if (info.cut[iRow]) {
3650	whichRow[nDelete++]=iRow;
3651	int start = rowStart[iRow];
3652	storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
3653	rowLength[iRow],column+start,elementByRow+start);
3654	}
3655	}
3656	solver->deleteRows(nDelete,whichRow);
3657	// and special matrix
3658	si->cleanMatrix()->deleteRows(nDelete,whichRow);
3659	delete [] whichRow;
3660	}
3661	}
3662	#endif
3663	// If we had a solution use it
3664	if (info.primalSolution) {
3665	solver->setColSolution(info.primalSolution);
3666	}
3667	// status
3668	if (info.rowStatus) {
3669	unsigned char * statusArray = lpSolver->statusArray();
3670	int i;
3671	for (i=0;i<info.numberColumns;i++)
3672	statusArray[i]=(char)info.columnStatus[i];
3673	statusArray+=info.numberColumns;
3674	for (i=0;i<info.numberRows;i++)
3675	statusArray[i]=(char)info.rowStatus[i];
3676	CoinWarmStartBasis * basis = lpSolver->getBasis();
3677	solver->setWarmStart(basis);
3678	delete basis;
3679	}
3680	freeArrays1(&info);
3681	// modify objective if necessary
3682	solver->setObjSense(info.direction);
3683	solver->setDblParam(OsiObjOffset,info.offset);
3684	if (info.offset) {
3685	sprintf(generalPrint,"Ampl objective offset is %g",
3686	info.offset);
3687	generalMessageHandler->message(CLP_GENERAL,generalMessages)
3688	<< generalPrint
3689	<<CoinMessageEol;
3690	}
3691	// Set integer variables (unless nonlinear when set)
3692	if (!info.nonLinear) {
3693	for (int i=info.numberColumns-info.numberIntegers;
3694	i<info.numberColumns;i++)
3695	solver->setInteger(i);
3696	}
3697	goodModel=true;
3698	// change argc etc
3699	argc = info.numberArguments;
3700	argv = const_cast<const char **>(info.arguments);
3701	}
3702	}
3703	#endif
3704	// default action on import
3705	int allowImportErrors=0;
3706	int keepImportNames=1;
3707	int doIdiot=-1;
3708	int outputFormat=2;
3709	int slpValue=-1;
3710	int cppValue=-1;
3711	int printOptions=0;
3712	int printMode=0;
3713	int presolveOptions=0;
3714	int substitution=3;
3715	int dualize=3;
3716	int doCrash=0;
3717	int doVector=0;
3718	int doSprint=-1;
3719	int doScaling=4;
3720	// set reasonable defaults
3721	int preSolve=5;
3722	int preProcess=4;
3723	bool useStrategy=false;
3724	bool preSolveFile=false;
3725	bool strongChanged=false;
3726
3727	double djFix=1.0e100;
3728	double gapRatio=1.0e100;
3729	double tightenFactor=0.0;
3730	const char dirsep = CoinFindDirSeparator();
3731	std::string directory;
3732	std::string dirSample;
3733	std::string dirNetlib;
3734	std::string dirMiplib;
3735	if (dirsep == '/') {
3736	directory = "./";
3737	dirSample = "../../Data/Sample/";
3738	dirNetlib = "../../Data/Netlib/";
3739	dirMiplib = "../../Data/miplib3/";
3740	} else {
3741	directory = ".\\";
3742	dirSample = "..\\..\\Data\\Sample\\";
3743	dirNetlib = "..\\..\\Data\\Netlib\\";
3744	dirMiplib = "..\\..\\Data\\miplib3\\";
3745	}
3746	std::string defaultDirectory = directory;
3747	std::string importFile ="";
3748	std::string exportFile ="default.mps";
3749	std::string importBasisFile ="";
3750	std::string importPriorityFile ="";
3751	std::string debugFile="";
3752	std::string printMask="";
3753	double * debugValues = NULL;
3754	int numberDebugValues = -1;
3755	int basisHasValues=0;
3756	std::string exportBasisFile ="default.bas";
3757	std::string saveFile ="default.prob";
3758	std::string restoreFile ="default.prob";
3759	std::string solutionFile ="stdout";
3760	std::string solutionSaveFile ="solution.file";
3761	int slog = whichParam(SOLVERLOGLEVEL,numberParameters_,parameters_);
3762	int log = whichParam(LOGLEVEL,numberParameters_,parameters_);
3763	double normalIncrement=model_.getCutoffIncrement();;
3764	if (testOsiParameters>=0) {
3765	// trying nonlinear - switch off some stuff
3766	preProcess=0;
3767	}
3768	// Set up likely cut generators and defaults
3769	int nodeStrategy=0;
3770	int doSOS=1;
3771	int verbose=0;
3772	CglGomory gomoryGen;
3773	// try larger limit
3774	gomoryGen.setLimitAtRoot(512);
3775	gomoryGen.setLimit(50);
3776	// set default action (0=off,1=on,2=root)
3777	int gomoryAction=3;
3778
3779	CglProbing probingGen;
3780	probingGen.setUsingObjective(1);
3781	probingGen.setMaxPass(1);
3782	probingGen.setMaxPassRoot(1);
3783	// Number of unsatisfied variables to look at
3784	probingGen.setMaxProbe(10);
3785	probingGen.setMaxProbeRoot(50);
3786	// How far to follow the consequences
3787	probingGen.setMaxLook(10);
3788	probingGen.setMaxLookRoot(50);
3789	probingGen.setMaxLookRoot(10);
3790	// Only look at rows with fewer than this number of elements
3791	probingGen.setMaxElements(200);
3792	probingGen.setRowCuts(3);
3793	// set default action (0=off,1=on,2=root)
3794	int probingAction=1;
3795
3796	CglKnapsackCover knapsackGen;
3797	//knapsackGen.switchOnExpensive();
3798	// set default action (0=off,1=on,2=root)
3799	int knapsackAction=3;
3800
3801	CglRedSplit redsplitGen;
3802	//redsplitGen.setLimit(100);
3803	// set default action (0=off,1=on,2=root)
3804	// Off as seems to give some bad cuts
3805	int redsplitAction=0;
3806
3807	CglClique cliqueGen(false,true);
3808	cliqueGen.setStarCliqueReport(false);
3809	cliqueGen.setRowCliqueReport(false);
3810	cliqueGen.setMinViolation(0.1);
3811	// set default action (0=off,1=on,2=root)
3812	int cliqueAction=3;
3813
3814	CglMixedIntegerRounding2 mixedGen;
3815	// set default action (0=off,1=on,2=root)
3816	int mixedAction=3;
3817
3818	CglFlowCover flowGen;
3819	// set default action (0=off,1=on,2=root)
3820	int flowAction=3;
3821
3822	CglTwomir twomirGen;
3823	twomirGen.setMaxElements(250);
3824	// set default action (0=off,1=on,2=root)
3825	int twomirAction=2;
3826	CglLandP landpGen;
3827	// set default action (0=off,1=on,2=root)
3828	int landpAction=0;
3829	CglResidualCapacity residualCapacityGen;
3830	// set default action (0=off,1=on,2=root)
3831	int residualCapacityAction=0;
3832	// Stored cuts
3833	bool storedCuts = false;
3834
3835	int useCosts=0;
3836	// don't use input solution
3837	int useSolution=0;
3838
3839	// total number of commands read
3840	int numberGoodCommands=0;
3841	// Set false if user does anything advanced
3842	bool defaultSettings=true;
3843
3844	// Hidden stuff for barrier
3845	int choleskyType = 0;
3846	int gamma=0;
3847	int scaleBarrier=0;
3848	int doKKT=0;
3849	int crossover=2; // do crossover unless quadratic
3850	// For names
3851	int lengthName = 0;
3852	std::vector<std::string> rowNames;
3853	std::vector<std::string> columnNames;
3854
3855	std::string field;
3856	if (!noPrinting_) {
3857	sprintf(generalPrint,"Coin Cbc and Clp Solver version %s, build %s",
3858	CBCVERSION,__DATE__);
3859	generalMessageHandler->message(CLP_GENERAL,generalMessages)
3860	<< generalPrint
3861	<<CoinMessageEol;
3862	// Print command line
3863	if (argc>1) {
3864	sprintf(generalPrint,"command line - ");
3865	for (int i=0;i<argc;i++) {
3866	if (!argv[i])
3867	break;
3868	sprintf(generalPrint+strlen(generalPrint),"%s ",argv[i]);
3869	}
3870	generalMessageHandler->message(CLP_GENERAL,generalMessages)
3871	<< generalPrint
3872	<<CoinMessageEol;
3873	}
3874	}
3875	while (1) {
3876	// next command
3877	field=CoinReadGetCommand(argc,argv);
3878	// Reset time
3879	time1 = CoinCpuTime();
3880	// adjust field if has odd trailing characters
3881	char temp [200];
3882	strcpy(temp,field.c_str());
3883	int length = strlen(temp);
3884	for (int k=length-1;k>=0;k--) {
3885	if (temp[k]<' ')
3886	length--;
3887	else
3888	break;
3889	}
3890	temp[length]='\0';
3891	field=temp;
3892	// exit if null or similar
3893	if (!field.length()) {
3894	if (numberGoodCommands==1&&goodModel) {
3895	// we just had file name - do branch and bound
3896	field="branch";
3897	} else if (!numberGoodCommands) {
3898	// let's give the sucker a hint
3899	std::cout
3900	<<"CoinSolver takes input from arguments ( - switches to stdin)"
3901	<<std::endl
3902	<<"Enter ? for list of commands or help"<<std::endl;
3903	field="-";
3904	} else {
3905	break;
3906	}
3907	}
3908
3909	// see if ? at end
3910	int numberQuery=0;
3911	if (field!="?"&&field!="???") {
3912	int length = field.length();
3913	int i;
3914	for (i=length-1;i>0;i--) {
3915	if (field[i]=='?')
3916	numberQuery++;
3917	else
3918	break;
3919	}
3920	field=field.substr(0,length-numberQuery);
3921	}
3922	// find out if valid command
3923	int iParam;
3924	int numberMatches=0;
3925	int firstMatch=-1;
3926	for ( iParam=0; iParam<numberParameters_; iParam++ ) {
3927	int match = parameters_[iParam].matches(field);
3928	if (match==1) {
3929	numberMatches = 1;
3930	firstMatch=iParam;
3931	break;
3932	} else {
3933	if (match&&firstMatch<0)
3934	firstMatch=iParam;
3935	numberMatches += match>>1;
3936	}
3937	}
3938	if (iParam<numberParameters_&&!numberQuery) {
3939	// found
3940	CbcOrClpParam found = parameters_[iParam];
3941	CbcOrClpParameterType type = found.type();
3942	int valid;
3943	numberGoodCommands++;
3944	if (type==BAB&&goodModel) {
3945	// check if any integers
3946	#ifdef COIN_HAS_ASL
3947	if (info.numberSos&&doSOS&&statusUserFunction_[0]) {
3948	// SOS
3949	numberSOS = info.numberSos;
3950	}
3951	#endif
3952	lpSolver = clpSolver->getModelPtr();
3953	if (!lpSolver->integerInformation()&&!numberSOS&&
3954	!clpSolver->numberSOS()&&!model_.numberObjects()&&!clpSolver->numberObjects())
3955	type=DUALSIMPLEX;
3956	}
3957	if (type==GENERALQUERY) {
3958	bool evenHidden=false;
3959	if ((verbose&8)!=0) {
3960	// even hidden
3961	evenHidden = true;
3962	verbose &= ~8;
3963	}
3964	#ifdef COIN_HAS_ASL
3965	if (verbose<4&&statusUserFunction_[0])
3966	verbose +=4;
3967	#endif
3968	if (verbose<4) {
3969	std::cout<<"In argument list keywords have leading - "
3970	", -stdin or just - switches to stdin"<<std::endl;
3971	std::cout<<"One command per line (and no -)"<<std::endl;
3972	std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
3973	std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
3974	std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
3975	std::cout<<"abcd value sets value"<<std::endl;
3976	std::cout<<"Commands are:"<<std::endl;
3977	} else {
3978	std::cout<<"Cbc options are set within AMPL with commands like:"<<std::endl<<std::endl;
3979	std::cout<<" option cbc_options \"cuts=root log=2 feas=on slog=1\""<<std::endl<<std::endl;
3980	std::cout<<"only maximize, dual, primal, help and quit are recognized without ="<<std::endl;
3981	}
3982	int maxAcross=5;
3983	if ((verbose%4)!=0)
3984	maxAcross=1;
3985	int limits[]={1,51,101,151,201,251,301,351,401};
3986	std::vector<std::string> types;
3987	types.push_back("Double parameters:");
3988	types.push_back("Branch and Cut double parameters:");
3989	types.push_back("Integer parameters:");
3990	types.push_back("Branch and Cut integer parameters:");
3991	types.push_back("Keyword parameters:");
3992	types.push_back("Branch and Cut keyword parameters:");
3993	types.push_back("Actions or string parameters:");
3994	types.push_back("Branch and Cut actions:");
3995	int iType;
3996	for (iType=0;iType<8;iType++) {
3997	int across=0;
3998	if ((verbose%4)!=0)
3999	std::cout<<std::endl;
4000	std::cout<<types[iType]<<std::endl;
4001	if ((verbose&2)!=0)
4002	std::cout<<std::endl;
4003	for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4004	int type = parameters_[iParam].type();
4005	if ((parameters_[iParam].displayThis()\|\|evenHidden)&&
4006	type>=limits[iType]
4007	&&type<limits[iType+1]) {
4008	// but skip if not useful for ampl (and in ampl mode)
4009	if (verbose>=4&&(parameters_[iParam].whereUsed()&4)==0)
4010	continue;
4011	if (!across) {
4012	if ((verbose&2)==0)
4013	std::cout<<" ";
4014	else
4015	std::cout<<"Command ";
4016	}
4017	std::cout<<parameters_[iParam].matchName()<<" ";
4018	across++;
4019	if (across==maxAcross) {
4020	across=0;
4021	if ((verbose%4)!=0) {
4022	// put out description as well
4023	if ((verbose&1)!=0)
4024	std::cout<<parameters_[iParam].shortHelp();
4025	std::cout<<std::endl;
4026	if ((verbose&2)!=0) {
4027	std::cout<<"---- description"<<std::endl;
4028	parameters_[iParam].printLongHelp();
4029	std::cout<<"----"<<std::endl<<std::endl;
4030	}
4031	} else {
4032	std::cout<<std::endl;
4033	}
4034	}
4035	}
4036	}
4037	if (across)
4038	std::cout<<std::endl;
4039	}
4040	} else if (type==FULLGENERALQUERY) {
4041	std::cout<<"Full list of commands is:"<<std::endl;
4042	int maxAcross=5;
4043	int limits[]={1,51,101,151,201,251,301,351,401};
4044	std::vector<std::string> types;
4045	types.push_back("Double parameters:");
4046	types.push_back("Branch and Cut double parameters:");
4047	types.push_back("Integer parameters:");
4048	types.push_back("Branch and Cut integer parameters:");
4049	types.push_back("Keyword parameters:");
4050	types.push_back("Branch and Cut keyword parameters:");
4051	types.push_back("Actions or string parameters:");
4052	types.push_back("Branch and Cut actions:");
4053	int iType;
4054	for (iType=0;iType<8;iType++) {
4055	int across=0;
4056	std::cout<<types[iType]<<" ";
4057	for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4058	int type = parameters_[iParam].type();
4059	if (type>=limits[iType]
4060	&&type<limits[iType+1]) {
4061	if (!across)
4062	std::cout<<" ";
4063	std::cout<<parameters_[iParam].matchName()<<" ";
4064	across++;
4065	if (across==maxAcross) {
4066	std::cout<<std::endl;
4067	across=0;
4068	}
4069	}
4070	}
4071	if (across)
4072	std::cout<<std::endl;
4073	}
4074	} else if (type<101) {
4075	// get next field as double
4076	double value = CoinReadGetDoubleField(argc,argv,&valid);
4077	if (!valid) {
4078	if (type<51) {
4079	parameters_[iParam].setDoubleParameter(lpSolver,value);
4080	} else if (type<81) {
4081	parameters_[iParam].setDoubleParameter(model_,value);
4082	} else {
4083	parameters_[iParam].setDoubleParameter(lpSolver,value);
4084	switch(type) {
4085	case DJFIX:
4086	djFix=value;
4087	if (goodModel&&djFix<1.0e20) {
4088	// do some fixing
4089	clpSolver = dynamic_cast< OsiClpSolverInterface*> (model_.solver());
4090	clpSolver->initialSolve();
4091	lpSolver = clpSolver->getModelPtr();
4092	int numberColumns = lpSolver->numberColumns();
4093	int i;
4094	const char * type = lpSolver->integerInformation();
4095	double * lower = lpSolver->columnLower();
4096	double * upper = lpSolver->columnUpper();
4097	double * solution = lpSolver->primalColumnSolution();
4098	double * dj = lpSolver->dualColumnSolution();
4099	int numberFixed=0;
4100	double dextra4 = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
4101	if (dextra4)
4102	printf("Multiple for continuous dj fixing is %g\n",dextra4);
4103	for (i=0;i<numberColumns;i++) {
4104	double djValue = dj[i];
4105	if (!type[i])
4106	djValue *= dextra4;
4107	if (type[i]\|\|dextra4) {
4108	double value = solution[i];
4109	if (value<lower[i]+1.0e-5&&djValue>djFix) {
4110	solution[i]=lower[i];
4111	upper[i]=lower[i];
4112	numberFixed++;
4113	} else if (value>upper[i]-1.0e-5&&djValue<-djFix) {
4114	solution[i]=upper[i];
4115	lower[i]=upper[i];
4116	numberFixed++;
4117	}
4118	}
4119	}
4120	sprintf(generalPrint,"%d columns fixed\n",numberFixed);
4121	generalMessageHandler->message(CLP_GENERAL,generalMessages)
4122	<< generalPrint
4123	<<CoinMessageEol;
4124	}
4125	break;
4126	case GAPRATIO:
4127	gapRatio=value;
4128	break;
4129	case TIGHTENFACTOR:
4130	tightenFactor=value;
4131	if(!complicatedInteger)
4132	defaultSettings=false; // user knows what she is doing
4133	break;
4134	default:
4135	break;
4136	}
4137	}
4138	} else if (valid==1) {
4139	std::cout<<" is illegal for double parameter "<<parameters_[iParam].name()<<" value remains "<<
4140	parameters_[iParam].doubleValue()<<std::endl;
4141	} else {
4142	std::cout<<parameters_[iParam].name()<<" has value "<<
4143	parameters_[iParam].doubleValue()<<std::endl;
4144	}
4145	} else if (type<201) {
4146	// get next field as int
4147	int value = CoinReadGetIntField(argc,argv,&valid);
4148	if (!valid) {
4149	if (type<151) {
4150	if (parameters_[iParam].type()==PRESOLVEPASS)
4151	preSolve = value;
4152	else if (parameters_[iParam].type()==IDIOT)
4153	doIdiot = value;
4154	else if (parameters_[iParam].type()==SPRINT)
4155	doSprint = value;
4156	else if (parameters_[iParam].type()==OUTPUTFORMAT)
4157	outputFormat = value;
4158	else if (parameters_[iParam].type()==SLPVALUE)
4159	slpValue = value;
4160	else if (parameters_[iParam].type()==CPP)
4161	cppValue = value;
4162	else if (parameters_[iParam].type()==PRESOLVEOPTIONS)
4163	presolveOptions = value;
4164	else if (parameters_[iParam].type()==PRINTOPTIONS)
4165	printOptions = value;
4166	else if (parameters_[iParam].type()==SUBSTITUTION)
4167	substitution = value;
4168	else if (parameters_[iParam].type()==DUALIZE)
4169	dualize = value;
4170	else if (parameters_[iParam].type()==PROCESSTUNE)
4171	tunePreProcess = value;
4172	else if (parameters_[iParam].type()==VERBOSE)
4173	verbose = value;
4174	parameters_[iParam].setIntParameter(lpSolver,value);
4175	} else {
4176	if (parameters_[iParam].type()==CUTPASS)
4177	cutPass = value;
4178	else if (parameters_[iParam].type()==CUTPASSINTREE)
4179	cutPassInTree = value;
4180	else if (parameters_[iParam].type()==STRONGBRANCHING\|\|
4181	parameters_[iParam].type()==NUMBERBEFORE)
4182	strongChanged=true;
4183	parameters_[iParam].setIntParameter(model_,value);
4184	}
4185	} else if (valid==1) {
4186	std::cout<<" is illegal for integer parameter "<<parameters_[iParam].name()<<" value remains "<<
4187	parameters_[iParam].intValue()<<std::endl;
4188	} else {
4189	std::cout<<parameters_[iParam].name()<<" has value "<<
4190	parameters_[iParam].intValue()<<std::endl;
4191	}
4192	} else if (type<301) {
4193	// one of several strings
4194	std::string value = CoinReadGetString(argc,argv);
4195	int action = parameters_[iParam].parameterOption(value);
4196	if (action<0) {
4197	if (value!="EOL") {
4198	// no match
4199	parameters_[iParam].printOptions();
4200	} else {
4201	// print current value
4202	std::cout<<parameters_[iParam].name()<<" has value "<<
4203	parameters_[iParam].currentOption()<<std::endl;
4204	}
4205	} else {
4206	parameters_[iParam].setCurrentOption(action,!noPrinting_);
4207	// for now hard wired
4208	switch (type) {
4209	case DIRECTION:
4210	if (action==0)
4211	lpSolver->setOptimizationDirection(1);
4212	else if (action==1)
4213	lpSolver->setOptimizationDirection(-1);
4214	else
4215	lpSolver->setOptimizationDirection(0);
4216	break;
4217	case DUALPIVOT:
4218	if (action==0) {
4219	ClpDualRowSteepest steep(3);
4220	lpSolver->setDualRowPivotAlgorithm(steep);
4221	} else if (action==1) {
4222	//ClpDualRowDantzig dantzig;
4223	ClpDualRowSteepest dantzig(5);
4224	lpSolver->setDualRowPivotAlgorithm(dantzig);
4225	} else if (action==2) {
4226	// partial steep
4227	ClpDualRowSteepest steep(2);
4228	lpSolver->setDualRowPivotAlgorithm(steep);
4229	} else {
4230	ClpDualRowSteepest steep;
4231	lpSolver->setDualRowPivotAlgorithm(steep);
4232	}
4233	break;
4234	case PRIMALPIVOT:
4235	if (action==0) {
4236	ClpPrimalColumnSteepest steep(3);
4237	lpSolver->setPrimalColumnPivotAlgorithm(steep);
4238	} else if (action==1) {
4239	ClpPrimalColumnSteepest steep(0);
4240	lpSolver->setPrimalColumnPivotAlgorithm(steep);
4241	} else if (action==2) {
4242	ClpPrimalColumnDantzig dantzig;
4243	lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
4244	} else if (action==3) {
4245	ClpPrimalColumnSteepest steep(4);
4246	lpSolver->setPrimalColumnPivotAlgorithm(steep);
4247	} else if (action==4) {
4248	ClpPrimalColumnSteepest steep(1);
4249	lpSolver->setPrimalColumnPivotAlgorithm(steep);
4250	} else if (action==5) {
4251	ClpPrimalColumnSteepest steep(2);
4252	lpSolver->setPrimalColumnPivotAlgorithm(steep);
4253	} else if (action==6) {
4254	ClpPrimalColumnSteepest steep(10);
4255	lpSolver->setPrimalColumnPivotAlgorithm(steep);
4256	}
4257	break;
4258	case SCALING:
4259	lpSolver->scaling(action);
4260	solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
4261	doScaling = action;
4262	break;
4263	case AUTOSCALE:
4264	lpSolver->setAutomaticScaling(action!=0);
4265	break;
4266	case SPARSEFACTOR:
4267	lpSolver->setSparseFactorization((1-action)!=0);
4268	break;
4269	case BIASLU:
4270	lpSolver->factorization()->setBiasLU(action);
4271	break;
4272	case PERTURBATION:
4273	if (action==0)
4274	lpSolver->setPerturbation(50);
4275	else
4276	lpSolver->setPerturbation(100);
4277	break;
4278	case ERRORSALLOWED:
4279	allowImportErrors = action;
4280	break;
4281	case INTPRINT:
4282	printMode=action;
4283	break;
4284	//case ALGORITHM:
4285	//algorithm = action;
4286	//defaultSettings=false; // user knows what she is doing
4287	//abort();
4288	//break;
4289	case KEEPNAMES:
4290	keepImportNames = 1-action;
4291	break;
4292	case PRESOLVE:
4293	if (action==0)
4294	preSolve = 5;
4295	else if (action==1)
4296	preSolve=0;
4297	else if (action==2)
4298	preSolve=10;
4299	else
4300	preSolveFile=true;
4301	break;
4302	case PFI:
4303	lpSolver->factorization()->setForrestTomlin(action==0);
4304	break;
4305	case CRASH:
4306	doCrash=action;
4307	break;
4308	case VECTOR:
4309	doVector=action;
4310	break;
4311	case MESSAGES:
4312	lpSolver->messageHandler()->setPrefix(action!=0);
4313	break;
4314	case CHOLESKY:
4315	choleskyType = action;
4316	break;
4317	case GAMMA:
4318	gamma=action;
4319	break;
4320	case BARRIERSCALE:
4321	scaleBarrier=action;
4322	break;
4323	case KKT:
4324	doKKT=action;
4325	break;
4326	case CROSSOVER:
4327	crossover=action;
4328	break;
4329	case SOS:
4330	doSOS=action;
4331	break;
4332	case GOMORYCUTS:
4333	defaultSettings=false; // user knows what she is doing
4334	gomoryAction = action;
4335	break;
4336	case PROBINGCUTS:
4337	defaultSettings=false; // user knows what she is doing
4338	probingAction = action;
4339	break;
4340	case KNAPSACKCUTS:
4341	defaultSettings=false; // user knows what she is doing
4342	knapsackAction = action;
4343	break;
4344	case REDSPLITCUTS:
4345	defaultSettings=false; // user knows what she is doing
4346	redsplitAction = action;
4347	break;
4348	case CLIQUECUTS:
4349	defaultSettings=false; // user knows what she is doing
4350	cliqueAction = action;
4351	break;
4352	case FLOWCUTS:
4353	defaultSettings=false; // user knows what she is doing
4354	flowAction = action;
4355	break;
4356	case MIXEDCUTS:
4357	defaultSettings=false; // user knows what she is doing
4358	mixedAction = action;
4359	break;
4360	case TWOMIRCUTS:
4361	defaultSettings=false; // user knows what she is doing
4362	twomirAction = action;
4363	break;
4364	case LANDPCUTS:
4365	defaultSettings=false; // user knows what she is doing
4366	landpAction = action;
4367	break;
4368	case RESIDCUTS:
4369	defaultSettings=false; // user knows what she is doing
4370	residualCapacityAction = action;
4371	break;
4372	case ROUNDING:
4373	defaultSettings=false; // user knows what she is doing
4374	break;
4375	case FPUMP:
4376	defaultSettings=false; // user knows what she is doing
4377	break;
4378	case RINS:
4379	break;
4380	case CUTSSTRATEGY:
4381	gomoryAction = action;
4382	probingAction = action;
4383	knapsackAction = action;
4384	cliqueAction = action;
4385	flowAction = action;
4386	mixedAction = action;
4387	twomirAction = action;
4388	//landpAction = action;
4389	parameters_[whichParam(GOMORYCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4390	parameters_[whichParam(PROBINGCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4391	parameters_[whichParam(KNAPSACKCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4392	parameters_[whichParam(CLIQUECUTS,numberParameters_,parameters_)].setCurrentOption(action);
4393	parameters_[whichParam(FLOWCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4394	parameters_[whichParam(MIXEDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4395	parameters_[whichParam(TWOMIRCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4396	if (!action) {
4397	redsplitAction = action;
4398	parameters_[whichParam(REDSPLITCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4399	landpAction = action;
4400	parameters_[whichParam(LANDPCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4401	residualCapacityAction = action;
4402	parameters_[whichParam(RESIDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4403	}
4404	break;
4405	case HEURISTICSTRATEGY:
4406	parameters_[whichParam(ROUNDING,numberParameters_,parameters_)].setCurrentOption(action);
4407	parameters_[whichParam(GREEDY,numberParameters_,parameters_)].setCurrentOption(action);
4408	parameters_[whichParam(COMBINE,numberParameters_,parameters_)].setCurrentOption(action);
4409	//parameters_[whichParam(LOCALTREE,numberParameters_,parameters_)].setCurrentOption(action);
4410	parameters_[whichParam(FPUMP,numberParameters_,parameters_)].setCurrentOption(action);
4411	break;
4412	case GREEDY:
4413	defaultSettings=false; // user knows what she is doing
4414	break;
4415	case COMBINE:
4416	defaultSettings=false; // user knows what she is doing
4417	break;
4418	case LOCALTREE:
4419	defaultSettings=false; // user knows what she is doing
4420	break;
4421	case COSTSTRATEGY:
4422	useCosts=action;
4423	break;
4424	case NODESTRATEGY:
4425	nodeStrategy=action;
4426	break;
4427	case PREPROCESS:
4428	preProcess = action;
4429	break;
4430	case USESOLUTION:
4431	useSolution = action;
4432	break;
4433	default:
4434	abort();
4435	}
4436	}
4437	} else {
4438	// action
4439	if (type==EXIT) {
4440	#ifdef COIN_HAS_ASL
4441	if(statusUserFunction_[0]) {
4442	if (info.numberIntegers\|\|info.numberBinary) {
4443	// integer
4444	} else {
4445	// linear
4446	}
4447	writeAmpl(&info);
4448	freeArrays2(&info);
4449	freeArgs(&info);
4450	}
4451	#endif
4452	break; // stop all
4453	}
4454	switch (type) {
4455	case DUALSIMPLEX:
4456	case PRIMALSIMPLEX:
4457	case SOLVECONTINUOUS:
4458	case BARRIER:
4459	if (goodModel) {
4460	double objScale =
4461	parameters_[whichParam(OBJSCALE2,numberParameters_,parameters_)].doubleValue();
4462	if (objScale!=1.0) {
4463	int iColumn;
4464	int numberColumns=lpSolver->numberColumns();
4465	double * dualColumnSolution =
4466	lpSolver->dualColumnSolution();
4467	ClpObjective * obj = lpSolver->objectiveAsObject();
4468	assert(dynamic_cast<ClpLinearObjective *> (obj));
4469	double offset;
4470	double * objective = obj->gradient(NULL,NULL,offset,true);
4471	for (iColumn=0;iColumn<numberColumns;iColumn++) {
4472	dualColumnSolution[iColumn] *= objScale;
4473	objective[iColumn] *= objScale;;
4474	}
4475	int iRow;
4476	int numberRows=lpSolver->numberRows();
4477	double * dualRowSolution =
4478	lpSolver->dualRowSolution();
4479	for (iRow=0;iRow<numberRows;iRow++)
4480	dualRowSolution[iRow] *= objScale;
4481	lpSolver->setObjectiveOffset(objScale*lpSolver->objectiveOffset());
4482	}
4483	ClpSolve::SolveType method;
4484	ClpSolve::PresolveType presolveType;
4485	ClpSimplex * model2 = lpSolver;
4486	if (dualize) {
4487	bool tryIt=true;
4488	double fractionColumn=1.0;
4489	double fractionRow=1.0;
4490	if (dualize==3) {
4491	dualize=1;
4492	int numberColumns=lpSolver->numberColumns();
4493	int numberRows=lpSolver->numberRows();
4494	if (numberRows<50000\|\|5*numberColumns>numberRows) {
4495	tryIt=false;
4496	} else {
4497	fractionColumn=0.1;
4498	fractionRow=0.1;
4499	}
4500	}
4501	if (tryIt) {
4502	model2 = ((ClpSimplexOther *) model2)->dualOfModel(fractionRow,fractionColumn);
4503	if (model2) {
4504	sprintf(generalPrint,"Dual of model has %d rows and %d columns",
4505	model2->numberRows(),model2->numberColumns());
4506	generalMessageHandler->message(CLP_GENERAL,generalMessages)
4507	<< generalPrint
4508	<<CoinMessageEol;
4509	model2->setOptimizationDirection(1.0);
4510	} else {
4511	model2 = lpSolver;
4512	dualize=0;
4513	}
4514	} else {
4515	dualize=0;
4516	}
4517	}
4518	if (noPrinting_)
4519	lpSolver->setLogLevel(0);
4520	ClpSolve solveOptions;
4521	solveOptions.setPresolveActions(presolveOptions);
4522	solveOptions.setSubstitution(substitution);
4523	if (preSolve!=5&&preSolve) {
4524	presolveType=ClpSolve::presolveNumber;
4525	if (preSolve<0) {
4526	preSolve = - preSolve;
4527	if (preSolve<=100) {
4528	presolveType=ClpSolve::presolveNumber;
4529	sprintf(generalPrint,"Doing %d presolve passes - picking up non-costed slacks",
4530	preSolve);
4531	generalMessageHandler->message(CLP_GENERAL,generalMessages)
4532	<< generalPrint
4533	<<CoinMessageEol;
4534	solveOptions.setDoSingletonColumn(true);
4535	} else {
4536	preSolve -=100;
4537	presolveType=ClpSolve::presolveNumberCost;
4538	sprintf(generalPrint,"Doing %d presolve passes - picking up costed slacks",
4539	preSolve);
4540	generalMessageHandler->message(CLP_GENERAL,generalMessages)
4541	<< generalPrint
4542	<<CoinMessageEol;
4543	}
4544	}
4545	} else if (preSolve) {
4546	presolveType=ClpSolve::presolveOn;
4547	} else {
4548	presolveType=ClpSolve::presolveOff;
4549	}
4550	solveOptions.setPresolveType(presolveType,preSolve);
4551	if (type==DUALSIMPLEX\|\|type==SOLVECONTINUOUS) {
4552	method=ClpSolve::useDual;
4553	} else if (type==PRIMALSIMPLEX) {
4554	method=ClpSolve::usePrimalorSprint;
4555	} else {
4556	method = ClpSolve::useBarrier;
4557	if (crossover==1) {
4558	method=ClpSolve::useBarrierNoCross;
4559	} else if (crossover==2) {
4560	ClpObjective * obj = lpSolver->objectiveAsObject();
4561	if (obj->type()>1) {
4562	method=ClpSolve::useBarrierNoCross;
4563	presolveType=ClpSolve::presolveOff;
4564	solveOptions.setPresolveType(presolveType,preSolve);
4565	}
4566	}
4567	}
4568	solveOptions.setSolveType(method);
4569	if(preSolveFile)
4570	presolveOptions \|= 0x40000000;
4571	solveOptions.setSpecialOption(4,presolveOptions);
4572	solveOptions.setSpecialOption(5,printOptions);
4573	if (doVector) {
4574	ClpMatrixBase * matrix = lpSolver->clpMatrix();
4575	if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
4576	ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
4577	clpMatrix->makeSpecialColumnCopy();
4578	}
4579	}
4580	if (method==ClpSolve::useDual) {
4581	// dual
4582	if (doCrash)
4583	solveOptions.setSpecialOption(0,1,doCrash); // crash
4584	else if (doIdiot)
4585	solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
4586	} else if (method==ClpSolve::usePrimalorSprint) {
4587	// primal
4588	// if slp turn everything off
4589	if (slpValue>0) {
4590	doCrash=false;
4591	doSprint=0;
4592	doIdiot=-1;
4593	solveOptions.setSpecialOption(1,10,slpValue); // slp
4594	method=ClpSolve::usePrimal;
4595	}
4596	if (doCrash) {
4597	solveOptions.setSpecialOption(1,1,doCrash); // crash
4598	} else if (doSprint>0) {
4599	// sprint overrides idiot
4600	solveOptions.setSpecialOption(1,3,doSprint); // sprint
4601	} else if (doIdiot>0) {
4602	solveOptions.setSpecialOption(1,2,doIdiot); // idiot
4603	} else if (slpValue<=0) {
4604	if (doIdiot==0) {
4605	if (doSprint==0)
4606	solveOptions.setSpecialOption(1,4); // all slack
4607	else
4608	solveOptions.setSpecialOption(1,9); // all slack or sprint
4609	} else {
4610	if (doSprint==0)
4611	solveOptions.setSpecialOption(1,8); // all slack or idiot
4612	else
4613	solveOptions.setSpecialOption(1,7); // initiative
4614	}
4615	}
4616	if (basisHasValues==-1)
4617	solveOptions.setSpecialOption(1,11); // switch off values
4618	} else if (method==ClpSolve::useBarrier\|\|method==ClpSolve::useBarrierNoCross) {
4619	int barrierOptions = choleskyType;
4620	if (scaleBarrier)
4621	barrierOptions \|= 8;
4622	if (doKKT)
4623	barrierOptions \|= 16;
4624	if (gamma)
4625	barrierOptions \|= 32*gamma;
4626	if (crossover==3)
4627	barrierOptions \|= 256; // try presolve in crossover
4628	solveOptions.setSpecialOption(4,barrierOptions);
4629	}
4630	model2->setMaximumSeconds(model_.getMaximumSeconds());
4631	#ifdef COIN_HAS_LINK
4632	OsiSolverInterface * coinSolver = model_.solver();
4633	OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
4634	if (!linkSolver) {
4635	model2->initialSolve(solveOptions);
4636	} else {
4637	// special solver
4638	int testOsiOptions = parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].intValue();
4639	double * solution = NULL;
4640	if (testOsiOptions<10) {
4641	solution = linkSolver->nonlinearSLP(slpValue>0 ? slpValue : 20 ,1.0e-5);
4642	} else if (testOsiOptions>=10) {
4643	CoinModel coinModel = *linkSolver->coinModel();
4644	ClpSimplex * tempModel = approximateSolution(coinModel,slpValue>0 ? slpValue : 50 ,1.0e-5,0);
4645	assert (tempModel);
4646	solution = CoinCopyOfArray(tempModel->primalColumnSolution(),coinModel.numberColumns());
4647	model2->setObjectiveValue(tempModel->objectiveValue());
4648	model2->setProblemStatus(tempModel->problemStatus());
4649	model2->setSecondaryStatus(tempModel->secondaryStatus());
4650	delete tempModel;
4651	}
4652	if (solution) {
4653	memcpy(model2->primalColumnSolution(),solution,
4654	CoinMin(model2->numberColumns(),linkSolver->coinModel()->numberColumns())*sizeof(double));
4655	delete [] solution;
4656	} else {
4657	printf("No nonlinear solution\n");
4658	}
4659	}
4660	#else
4661	model2->initialSolve(solveOptions);
4662	#endif
4663	{
4664	// map states
4665	/* clp status
4666	-1 - unknown e.g. before solve or if postSolve says not optimal
4667	0 - optimal
4668	1 - primal infeasible
4669	2 - dual infeasible
4670	3 - stopped on iterations or time
4671	4 - stopped due to errors
4672	5 - stopped by event handler (virtual int ClpEventHandler::event()) */
4673	/* cbc status
4674	-1 before branchAndBound
4675	0 finished - check isProvenOptimal or isProvenInfeasible to see if solution found
4676	(or check value of best solution)
4677	1 stopped - on maxnodes, maxsols, maxtime
4678	2 difficulties so run was abandoned
4679	(5 event user programmed event occurred) */
4680	/* clp secondary status of problem - may get extended
4681	0 - none
4682	1 - primal infeasible because dual limit reached OR probably primal
4683	infeasible but can't prove it (main status 4)
4684	2 - scaled problem optimal - unscaled problem has primal infeasibilities
4685	3 - scaled problem optimal - unscaled problem has dual infeasibilities
4686	4 - scaled problem optimal - unscaled problem has primal and dual infeasibilities
4687	5 - giving up in primal with flagged variables
4688	6 - failed due to empty problem check
4689	7 - postSolve says not optimal
4690	8 - failed due to bad element check
4691	9 - status was 3 and stopped on time
4692	100 up - translation of enum from ClpEventHandler
4693	*/
4694	/* cbc secondary status of problem
4695	-1 unset (status_ will also be -1)
4696	0 search completed with solution
4697	1 linear relaxation not feasible (or worse than cutoff)
4698	2 stopped on gap
4699	3 stopped on nodes
4700	4 stopped on time
4701	5 stopped on user event
4702	6 stopped on solutions
4703	7 linear relaxation unbounded
4704	*/
4705	int iStatus = model2->status();
4706	int iStatus2 = model2->secondaryStatus();
4707	if (iStatus==0) {
4708	iStatus2=0;
4709	} else if (iStatus==1) {
4710	iStatus=0;
4711	iStatus2=1; // say infeasible
4712	} else if (iStatus==2) {
4713	iStatus=0;
4714	iStatus2=7; // say unbounded
4715	} else if (iStatus==3) {
4716	iStatus=1;
4717	if (iStatus2==9)
4718	iStatus2=4;
4719	else
4720	iStatus2=3; // Use nodes - as closer than solutions
4721	} else if (iStatus==4) {
4722	iStatus=2; // difficulties
4723	iStatus2=0;
4724	}
4725	model_.setProblemStatus(iStatus);
4726	model_.setSecondaryStatus(iStatus2);
4727	assert (lpSolver==clpSolver->getModelPtr());
4728	assert (clpSolver==model_.solver());
4729	clpSolver->setWarmStart(NULL);
4730	// and in babModel if exists
4731	if (babModel_) {
4732	babModel_->setProblemStatus(iStatus);
4733	babModel_->setSecondaryStatus(iStatus2);
4734	}
4735	#ifdef NEW_STYLE_SOLVER
4736	int returnCode = callBack_->callBack(&model_,1);
4737	#else
4738	int returnCode=callBack(&model,1);
4739	#endif
4740	if (returnCode) {
4741	// exit if user wants
4742	#ifdef NEW_STYLE_SOLVER
4743	updateModel(model2,returnMode);
4744	#else
4745	delete babModel_;
4746	babModel_ = NULL;
4747	#endif
4748	return returnCode;
4749	}
4750	}
4751	basisHasValues=1;
4752	if (dualize) {
4753	int returnCode=((ClpSimplexOther *) lpSolver)->restoreFromDual(model2);
4754	if (model2->status()==3)
4755	returnCode=0;
4756	delete model2;
4757	if (returnCode&&dualize!=2)
4758	lpSolver->primal(1);
4759	model2=lpSolver;
4760	}
4761	#ifdef NEW_STYLE_SOLVER
4762	updateModel(model2,returnMode);
4763	for (iUser=0;iUser<numberUserFunctions_;iUser++) {
4764	if (statusUserFunction_[iUser])
4765	userFunction_[iUser]->exportSolution(this,1);
4766	}
4767	#endif
4768	#ifdef COIN_HAS_ASL
4769	if (statusUserFunction_[0]) {
4770	double value = model2->getObjValue()*model2->getObjSense();
4771	char buf[300];
4772	int pos=0;
4773	int iStat = model2->status();
4774	if (iStat==0) {
4775	pos += sprintf(buf+pos,"optimal," );
4776	} else if (iStat==1) {
4777	// infeasible
4778	pos += sprintf(buf+pos,"infeasible,");
4779	} else if (iStat==2) {
4780	// unbounded
4781	pos += sprintf(buf+pos,"unbounded,");
4782	} else if (iStat==3) {
4783	pos += sprintf(buf+pos,"stopped on iterations or time,");
4784	} else if (iStat==4) {
4785	iStat = 7;
4786	pos += sprintf(buf+pos,"stopped on difficulties,");
4787	} else if (iStat==5) {
4788	iStat = 3;
4789	pos += sprintf(buf+pos,"stopped on ctrl-c,");
4790	} else {
4791	pos += sprintf(buf+pos,"status unknown,");
4792	iStat=6;
4793	}
4794	info.problemStatus=iStat;
4795	info.objValue = value;
4796	pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
4797	value);
4798	sprintf(buf+pos,"\n%d iterations",
4799	model2->getIterationCount());
4800	free(info.primalSolution);
4801	int numberColumns=model2->numberColumns();
4802	info.primalSolution = (double ) malloc(numberColumnssizeof(double));
4803	CoinCopyN(model2->primalColumnSolution(),numberColumns,info.primalSolution);
4804	int numberRows = model2->numberRows();
4805	free(info.dualSolution);
4806	info.dualSolution = (double ) malloc(numberRowssizeof(double));
4807	CoinCopyN(model2->dualRowSolution(),numberRows,info.dualSolution);
4808	CoinWarmStartBasis * basis = model2->getBasis();
4809	free(info.rowStatus);
4810	info.rowStatus = (int ) malloc(numberRowssizeof(int));
4811	free(info.columnStatus);
4812	info.columnStatus = (int ) malloc(numberColumnssizeof(int));
4813	// Put basis in
4814	int i;
4815	// free,basic,ub,lb are 0,1,2,3
4816	for (i=0;i<numberRows;i++) {
4817	CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
4818	info.rowStatus[i]=status;
4819	}
4820	for (i=0;i<numberColumns;i++) {
4821	CoinWarmStartBasis::Status status = basis->getStructStatus(i);
4822	info.columnStatus[i]=status;
4823	}
4824	// put buffer into info
4825	strcpy(info.buffer,buf);
4826	delete basis;
4827	}
4828	#endif
4829	} else {
4830	#ifndef DISALLOW_PRINTING
4831	std::cout<<"** Current model not valid"<<std::endl;
4832	#endif
4833	}
4834	break;
4835	case STATISTICS:
4836	if (goodModel) {
4837	// If presolve on look at presolved
4838	bool deleteModel2=false;
4839	ClpSimplex * model2 = lpSolver;
4840	if (preSolve) {
4841	ClpPresolve pinfo;
4842	int presolveOptions2 = presolveOptions&~0x40000000;
4843	if ((presolveOptions2&0xffff)!=0)
4844	pinfo.setPresolveActions(presolveOptions2);
4845	pinfo.setSubstitution(substitution);
4846	if ((printOptions&1)!=0)
4847	pinfo.statistics();
4848	double presolveTolerance =
4849	parameters_[whichParam(PRESOLVETOLERANCE,numberParameters_,parameters_)].doubleValue();
4850	model2 =
4851	pinfo.presolvedModel(*lpSolver,presolveTolerance,
4852	true,preSolve);
4853	if (model2) {
4854	printf("Statistics for presolved model\n");
4855	deleteModel2=true;
4856	} else {
4857	printf("Presolved model looks infeasible - will use unpresolved\n");
4858	model2 = lpSolver;
4859	}
4860	} else {
4861	printf("Statistics for unpresolved model\n");
4862	model2 = lpSolver;
4863	}
4864	statistics(lpSolver,model2);
4865	if (deleteModel2)
4866	delete model2;
4867	} else {
4868	#ifndef DISALLOW_PRINTING
4869	std::cout<<"** Current model not valid"<<std::endl;
4870	#endif
4871	}
4872	break;
4873	case TIGHTEN:
4874	if (goodModel) {
4875	int numberInfeasibilities = lpSolver->tightenPrimalBounds();
4876	if (numberInfeasibilities)
4877	std::cout<<"** Analysis indicates model infeasible"<<std::endl;
4878	} else {
4879	#ifndef DISALLOW_PRINTING
4880	std::cout<<"** Current model not valid"<<std::endl;
4881	#endif
4882	}
4883	break;
4884	case PLUSMINUS:
4885	if (goodModel) {
4886	ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
4887	ClpPackedMatrix* clpMatrix =
4888	dynamic_cast< ClpPackedMatrix*>(saveMatrix);
4889	if (clpMatrix) {
4890	ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
4891	if (newMatrix->getIndices()) {
4892	lpSolver->replaceMatrix(newMatrix);
4893	delete saveMatrix;
4894	std::cout<<"Matrix converted to +- one matrix"<<std::endl;
4895	} else {
4896	std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
4897	}
4898	} else {
4899	std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
4900	}
4901	} else {
4902	#ifndef DISALLOW_PRINTING
4903	std::cout<<"** Current model not valid"<<std::endl;
4904	#endif
4905	}
4906	break;
4907	case OUTDUPROWS:
4908	if (goodModel) {
4909	int numberRows = clpSolver->getNumRows();
4910	//int nOut = outDupRow(clpSolver);
4911	CglDuplicateRow dupcuts(clpSolver);
4912	storedCuts = dupcuts.outDuplicates(clpSolver)!=0;
4913	int nOut = numberRows-clpSolver->getNumRows();
4914	if (nOut&&!noPrinting_)
4915	sprintf(generalPrint,"%d rows eliminated",nOut);
4916	generalMessageHandler->message(CLP_GENERAL,generalMessages)
4917	<< generalPrint
4918	<<CoinMessageEol;
4919	} else {
4920	#ifndef DISALLOW_PRINTING
4921	std::cout<<"** Current model not valid"<<std::endl;
4922	#endif
4923	}
4924	break;
4925	case NETWORK:
4926	if (goodModel) {
4927	ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
4928	ClpPackedMatrix* clpMatrix =
4929	dynamic_cast< ClpPackedMatrix*>(saveMatrix);
4930	if (clpMatrix) {
4931	ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
4932	if (newMatrix->getIndices()) {
4933	lpSolver->replaceMatrix(newMatrix);
4934	delete saveMatrix;
4935	std::cout<<"Matrix converted to network matrix"<<std::endl;
4936	} else {
4937	std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
4938	}
4939	} else {
4940	std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
4941	}
4942	} else {
4943	#ifndef DISALLOW_PRINTING
4944	std::cout<<"** Current model not valid"<<std::endl;
4945	#endif
4946	}
4947	break;
4948	case DOHEURISTIC:
4949	if (goodModel) {
4950	int vubAction = parameters_[whichParam(VUBTRY,numberParameters_,parameters_)].intValue();
4951	if (vubAction!=-1) {
4952	// look at vubs
4953	// Just ones which affect >= extra3
4954	int extra3 = parameters_[whichParam(EXTRA3,numberParameters_,parameters_)].intValue();
4955	/* 3 is fraction of integer variables fixed (0.97)
4956	4 is fraction of all variables fixed (0.0)
4957	*/
4958	double dextra[5];
4959	dextra[1] = parameters_[whichParam(DEXTRA1,numberParameters_,parameters_)].doubleValue();
4960	dextra[2] = parameters_[whichParam(DEXTRA2,numberParameters_,parameters_)].doubleValue();
4961	dextra[3] = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
4962	dextra[4] = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
4963	if (!dextra[3])
4964	dextra[3] = 0.97;
4965	//OsiClpSolverInterface * newSolver =
4966	fixVubs(model_,extra3,vubAction,generalMessageHandler,
4967	debugValues,dextra);
4968	//assert (!newSolver);
4969	}
4970	// Actually do heuristics
4971	doHeuristics(&model_,2);
4972	if (model_.bestSolution()) {
4973	model_.setProblemStatus(1);
4974	model_.setSecondaryStatus(6);
4975	#ifdef COIN_HAS_ASL
4976	if (statusUserFunction_[0]) {
4977	double value = model_.getObjValue();
4978	char buf[300];
4979	int pos=0;
4980	pos += sprintf(buf+pos,"feasible,");
4981	info.problemStatus=0;
4982	info.objValue = value;
4983	pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
4984	value);
4985	sprintf(buf+pos,"\n0 iterations");
4986	free(info.primalSolution);
4987	int numberColumns=lpSolver->numberColumns();
4988	info.primalSolution = (double ) malloc(numberColumnssizeof(double));
4989	CoinCopyN(model_.bestSolution(),numberColumns,info.primalSolution);
4990	int numberRows = lpSolver->numberRows();
4991	free(info.dualSolution);
4992	info.dualSolution = (double ) malloc(numberRowssizeof(double));
4993	CoinZeroN(info.dualSolution,numberRows);
4994	CoinWarmStartBasis * basis = lpSolver->getBasis();
4995	free(info.rowStatus);
4996	info.rowStatus = (int ) malloc(numberRowssizeof(int));
4997	free(info.columnStatus);
4998	info.columnStatus = (int ) malloc(numberColumnssizeof(int));
4999	// Put basis in
5000	int i;
5001	// free,basic,ub,lb are 0,1,2,3
5002	for (i=0;i<numberRows;i++) {
5003	CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
5004	info.rowStatus[i]=status;
5005	}
5006	for (i=0;i<numberColumns;i++) {
5007	CoinWarmStartBasis::Status status = basis->getStructStatus(i);
5008	info.columnStatus[i]=status;
5009	}
5010	// put buffer into info
5011	strcpy(info.buffer,buf);
5012	delete basis;
5013	}
5014	#endif
5015	}
5016	#ifdef NEW_STYLE_SOLVER
5017	int returnCode = callBack_->callBack(&model_,6);
5018	#else
5019	int returnCode=callBack(&model,6);
5020	#endif
5021	if (returnCode) {
5022	// exit if user wants
5023	#ifdef NEW_STYLE_SOLVER
5024	updateModel(NULL,returnMode);
5025	#else
5026	delete babModel_;
5027	babModel_ = NULL;
5028	#endif
5029	return returnCode;
5030	}
5031	}
5032	break;
5033	case MIPLIB:
5034	// User can set options - main difference is lack of model and CglPreProcess
5035	goodModel=true;
5036	/*
5037	Run branch-and-cut. First set a few options -- node comparison, scaling.
5038	Print elapsed time at the end.
5039	*/
5040	case BAB: // branchAndBound
5041	case STRENGTHEN:
5042	if (goodModel) {
5043	bool miplib = type==MIPLIB;
5044	int logLevel = parameters_[slog].intValue();
5045	// Reduce printout
5046	if (logLevel<=1)
5047	model_.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
5048	{
5049	OsiSolverInterface * solver = model_.solver();
5050	OsiClpSolverInterface * si =
5051	dynamic_cast<OsiClpSolverInterface *>(solver) ;
5052	assert (si != NULL);
5053	si->getModelPtr()->scaling(doScaling);
5054	// See if quadratic
5055	#ifdef COIN_HAS_LINK
5056	if (!complicatedInteger) {
5057	ClpSimplex * lpSolver = si->getModelPtr();
5058	ClpQuadraticObjective * obj = (dynamic_cast< ClpQuadraticObjective*>(lpSolver->objectiveAsObject()));
5059	if (obj) {
5060	preProcess=0;
5061	int testOsiOptions = parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].intValue();
5062	parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].setIntValue(CoinMax(0,testOsiOptions));
5063	// create coin model
5064	coinModel = lpSolver->createCoinModel();
5065	assert (coinModel);
5066	// load from coin model
5067	OsiSolverLink solver1;
5068	OsiSolverInterface * solver2 = solver1.clone();
5069	model_.assignSolver(solver2,false);
5070	OsiSolverLink * si =
5071	dynamic_cast<OsiSolverLink *>(model_.solver()) ;
5072	assert (si != NULL);
5073	si->setDefaultMeshSize(0.001);
5074	// need some relative granularity
5075	si->setDefaultBound(100.0);
5076	double dextra3 = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
5077	if (dextra3)
5078	si->setDefaultMeshSize(dextra3);
5079	si->setDefaultBound(1000.0);
5080	si->setIntegerPriority(1000);
5081	si->setBiLinearPriority(10000);
5082	si->setSpecialOptions2(2+4+8);
5083	CoinModel * model2 = (CoinModel *) coinModel;
5084	si->load(*model2,true, parameters_[log].intValue());
5085	// redo
5086	solver = model_.solver();
5087	clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
5088	lpSolver = clpSolver->getModelPtr();
5089	clpSolver->messageHandler()->setLogLevel(0) ;
5090	testOsiParameters=0;
5091	complicatedInteger=2; // allow cuts
5092	OsiSolverInterface * coinSolver = model_.solver();
5093	OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
5094	if (linkSolver->quadraticModel()) {
5095	ClpSimplex * qp = linkSolver->quadraticModel();
5096	//linkSolver->nonlinearSLP(CoinMax(slpValue,10),1.0e-5);
5097	qp->nonlinearSLP(CoinMax(slpValue,40),1.0e-5);
5098	qp->primal(1);
5099	OsiSolverLinearizedQuadratic solver2(qp);
5100	const double * solution=NULL;
5101	// Reduce printout
5102	solver2.setHintParam(OsiDoReducePrint,true,OsiHintTry);
5103	CbcModel model2(solver2);
5104	// Now do requested saves and modifications
5105	CbcModel * cbcModel = & model2;
5106	OsiSolverInterface * osiModel = model2.solver();
5107	OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
5108	ClpSimplex * clpModel = osiclpModel->getModelPtr();
5109
5110	// Set changed values
5111
5112	CglProbing probing;
5113	probing.setMaxProbe(10);
5114	probing.setMaxLook(10);
5115	probing.setMaxElements(200);
5116	probing.setMaxProbeRoot(50);
5117	probing.setMaxLookRoot(10);
5118	probing.setRowCuts(3);
5119	probing.setUsingObjective(true);
5120	cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
5121	cbcModel->cutGenerator(0)->setTiming(true);
5122
5123	CglGomory gomory;
5124	gomory.setLimitAtRoot(512);
5125	cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
5126	cbcModel->cutGenerator(1)->setTiming(true);
5127
5128	CglKnapsackCover knapsackCover;
5129	cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
5130	cbcModel->cutGenerator(2)->setTiming(true);
5131
5132	CglRedSplit redSplit;
5133	cbcModel->addCutGenerator(&redSplit,-99,"RedSplit",true,false,false,-100,-1,-1);
5134	cbcModel->cutGenerator(3)->setTiming(true);
5135
5136	CglClique clique;
5137	clique.setStarCliqueReport(false);
5138	clique.setRowCliqueReport(false);
5139	clique.setMinViolation(0.1);
5140	cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
5141	cbcModel->cutGenerator(4)->setTiming(true);
5142
5143	CglMixedIntegerRounding2 mixedIntegerRounding2;
5144	cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
5145	cbcModel->cutGenerator(5)->setTiming(true);
5146
5147	CglFlowCover flowCover;
5148	cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
5149	cbcModel->cutGenerator(6)->setTiming(true);
5150
5151	CglTwomir twomir;
5152	twomir.setMaxElements(250);
5153	cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
5154	cbcModel->cutGenerator(7)->setTiming(true);
5155
5156	CbcHeuristicFPump heuristicFPump(*cbcModel);
5157	heuristicFPump.setWhen(13);
5158	heuristicFPump.setMaximumPasses(20);
5159	heuristicFPump.setMaximumRetries(7);
5160	heuristicFPump.setHeuristicName("feasibility pump");
5161	heuristicFPump.setInitialWeight(1);
5162	heuristicFPump.setFractionSmall(0.6);
5163	cbcModel->addHeuristic(&heuristicFPump);
5164
5165	CbcRounding rounding(*cbcModel);
5166	rounding.setHeuristicName("rounding");
5167	cbcModel->addHeuristic(&rounding);
5168
5169	CbcHeuristicLocal heuristicLocal(*cbcModel);
5170	heuristicLocal.setHeuristicName("combine solutions");
5171	heuristicLocal.setSearchType(1);
5172	heuristicLocal.setFractionSmall(0.6);
5173	cbcModel->addHeuristic(&heuristicLocal);
5174
5175	CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
5176	heuristicGreedyCover.setHeuristicName("greedy cover");
5177	cbcModel->addHeuristic(&heuristicGreedyCover);
5178
5179	CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
5180	heuristicGreedyEquality.setHeuristicName("greedy equality");
5181	cbcModel->addHeuristic(&heuristicGreedyEquality);
5182
5183	CbcCompareDefault compare;
5184	cbcModel->setNodeComparison(compare);
5185	cbcModel->setNumberBeforeTrust(5);
5186	cbcModel->setSpecialOptions(2);
5187	cbcModel->messageHandler()->setLogLevel(1);
5188	cbcModel->setMaximumCutPassesAtRoot(-100);
5189	cbcModel->setMaximumCutPasses(1);
5190	cbcModel->setMinimumDrop(0.05);
5191	// For branchAndBound this may help
5192	clpModel->defaultFactorizationFrequency();
5193	clpModel->setDualBound(1.0001e+08);
5194	clpModel->setPerturbation(50);
5195	osiclpModel->setSpecialOptions(193);
5196	osiclpModel->messageHandler()->setLogLevel(0);
5197	osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
5198	osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
5199	// You can save some time by switching off message building
5200	// clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
5201
5202	// Solve
5203
5204	cbcModel->initialSolve();
5205	if (clpModel->tightenPrimalBounds()!=0) {
5206	#ifndef DISALLOW_PRINTING
5207	std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
5208	#endif
5209	break;
5210	}
5211	clpModel->dual(); // clean up
5212	cbcModel->initialSolve();
5213	#ifdef CBC_THREAD
5214	int numberThreads =parameters_[whichParam(THREADS,numberParameters_,parameters_)].intValue();
5215	cbcModel->setNumberThreads(numberThreads%100);
5216	cbcModel->setThreadMode(numberThreads/100);
5217	#endif
5218	cbcModel->branchAndBound();
5219	OsiSolverLinearizedQuadratic * solver3 = dynamic_cast<OsiSolverLinearizedQuadratic *> (model2.solver());
5220	assert (solver3);
5221	solution = solver3->bestSolution();
5222	double bestObjectiveValue = solver3->bestObjectiveValue();
5223	linkSolver->setBestObjectiveValue(bestObjectiveValue);
5224	linkSolver->setBestSolution(solution,solver3->getNumCols());
5225	CbcHeuristicDynamic3 dynamic(model_);
5226	dynamic.setHeuristicName("dynamic pass thru");
5227	model_.addHeuristic(&dynamic);
5228	// if convex
5229	if ((linkSolver->specialOptions2()&4)!=0) {
5230	int numberColumns = coinModel->numberColumns();
5231	assert (linkSolver->objectiveVariable()==numberColumns);
5232	// add OA cut
5233	double offset;
5234	double * gradient = new double [numberColumns+1];
5235	memcpy(gradient,qp->objectiveAsObject()->gradient(qp,solution,offset,true,2),
5236	numberColumns*sizeof(double));
5237	double rhs = 0.0;
5238	int * column = new int[numberColumns+1];
5239	int n=0;
5240	for (int i=0;i<numberColumns;i++) {
5241	double value = gradient[i];
5242	if (fabs(value)>1.0e-12) {
5243	gradient[n]=value;
5244	rhs += value*solution[i];
5245	column[n++]=i;
5246	}
5247	}
5248	gradient[n]=-1.0;
5249	column[n++]=numberColumns;
5250	storedAmpl.addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
5251	delete [] gradient;
5252	delete [] column;
5253	}
5254	// could do three way branching round a) continuous b) best solution
5255	printf("obj %g\n",bestObjectiveValue);
5256	linkSolver->initialSolve();
5257	�