Context Navigation

CbcLinked.cpp @ 817

Last change on this file since 817 was 817, checked in by forrest, 13 years ago
so CbcLinked? will compile
File size: 244.5 KB

Line
1	// Copyright (C) 2006, International Business Machines
2	// Corporation and others. All Rights Reserved.
3	#include "CbcConfig.h"
4
5	#include "CoinTime.hpp"
6
7	#include "CoinHelperFunctions.hpp"
8	#include "CoinModel.hpp"
9	#include "ClpSimplex.hpp"
10	// returns jColumn (-2 if linear term, -1 if unknown) and coefficient
11	static
12	int decodeBit(char * phrase, char * & nextPhrase, double & coefficient, bool ifFirst, const CoinModel & model)
13	{
14	char * pos = phrase;
15	// may be leading - (or +)
16	char * pos2 = pos;
17	double value=1.0;
18	if (pos2=='-'\|\|pos2=='+')
19	pos2++;
20	// next terminator * or + or -
21	while (*pos2) {
22	if (pos2=='') {
23	break;
24	} else if (pos2=='-'\|\|pos2=='+') {
25	if (pos2==pos\|\|*(pos2-1)!='e')
26	break;
27	}
28	pos2++;
29	}
30	// if * must be number otherwise must be name
31	if (pos2=='') {
32	char * pos3 = pos;
33	while (pos3!=pos2) {
34	char x = *pos3;
35	pos3++;
36	assert ((x>='0'&&x<='9')\|\|x=='.'\|\|x=='+'\|\|x=='-'\|\|x=='e');
37	}
38	char saved = *pos2;
39	*pos2='\0';
40	value = atof(pos);
41	*pos2=saved;
42	// and down to next
43	pos2++;
44	pos=pos2;
45	while (*pos2) {
46	if (pos2=='-'\|\|pos2=='+')
47	break;
48	pos2++;
49	}
50	}
51	char saved = *pos2;
52	*pos2='\0';
53	// now name
54	// might have + or -
55	if (*pos=='+') {
56	pos++;
57	} else if (*pos=='-') {
58	pos++;
59	assert (value==1.0);
60	value = - value;
61	}
62	int jColumn = model.column(pos);
63	// must be column unless first when may be linear term
64	if (jColumn<0) {
65	if (ifFirst) {
66	char * pos3 = pos;
67	while (pos3!=pos2) {
68	char x = *pos3;
69	pos3++;
70	assert ((x>='0'&&x<='9')\|\|x=='.'\|\|x=='+'\|\|x=='-'\|\|x=='e');
71	}
72	assert(*pos2=='\0');
73	// keep possible -
74	value = value * atof(pos);
75	jColumn=-2;
76	} else {
77	// bad
78	*pos2=saved;
79	printf("bad nonlinear term %s\n",phrase);
80	abort();
81	}
82	}
83	*pos2=saved;
84	pos=pos2;
85	coefficient=value;
86	nextPhrase = pos;
87	return jColumn;
88	}
89	#include "ClpQuadraticObjective.hpp"
90	#include <cassert>
91	#if defined(_MSC_VER)
92	// Turn off compiler warning about long names
93	# pragma warning(disable:4786)
94	#endif
95	#include "CbcLinked.hpp"
96	#include "CoinIndexedVector.hpp"
97	#include "CoinMpsIO.hpp"
98	//#include "OsiSolverLink.hpp"
99	//#include "OsiBranchLink.hpp"
100	#include "ClpPackedMatrix.hpp"
101	#include "CoinTime.hpp"
102	#include "CbcModel.hpp"
103	#include "CbcCutGenerator.hpp"
104	#include "CglStored.hpp"
105	#include "CglPreProcess.hpp"
106	#include "CglGomory.hpp"
107	#include "CglProbing.hpp"
108	#include "CglKnapsackCover.hpp"
109	#include "CglRedSplit.hpp"
110	#include "CglClique.hpp"
111	#include "CglFlowCover.hpp"
112	#include "CglMixedIntegerRounding2.hpp"
113	#include "CglTwomir.hpp"
114	#include "CglDuplicateRow.hpp"
115	#include "CbcHeuristicFPump.hpp"
116	#include "CbcHeuristic.hpp"
117	#include "CbcHeuristicLocal.hpp"
118	#include "CbcHeuristicGreedy.hpp"
119	#include "ClpLinearObjective.hpp"
120	#include "CbcBranchActual.hpp"
121	#include "CbcCompareActual.hpp"
122	//#############################################################################
123	// Solve methods
124	//#############################################################################
125	void OsiSolverLink::initialSolve()
126	{
127	//writeMps("yy");
128	//exit(77);
129	specialOptions_ =0;
130	modelPtr_->setWhatsChanged(0);
131	if (numberVariables_) {
132	CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
133	// update all bounds before coefficients
134	for (int i=0;i<numberVariables_;i++ ) {
135	info_[i].updateBounds(modelPtr_);
136	}
137	int updated = updateCoefficients(modelPtr_,temp);
138	if (updated\|\|1) {
139	temp->removeGaps(1.0e-14);
140	ClpMatrixBase * save = modelPtr_->clpMatrix();
141	ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
142	assert (clpMatrix);
143	if (save->getNumRows()>temp->getNumRows()) {
144	// add in cuts
145	int numberRows = temp->getNumRows();
146	int * which = new int[numberRows];
147	for (int i=0;i<numberRows;i++)
148	which[i]=i;
149	save->deleteRows(numberRows,which);
150	delete [] which;
151	temp->bottomAppendPackedMatrix(*clpMatrix->matrix());
152	}
153	modelPtr_->replaceMatrix(temp,true);
154	} else {
155	delete temp;
156	}
157	}
158	if (0) {
159	const double * lower = getColLower();
160	const double * upper = getColUpper();
161	int n=0;
162	for (int i=84;i<84+16;i++) {
163	if (lower[i]+0.01<upper[i]) {
164	n++;
165	}
166	}
167	if (!n)
168	writeMps("sol_query");
169
170	}
171	//static int iPass=0;
172	//char temp[50];
173	//iPass++;
174	//sprintf(temp,"cc%d",iPass);
175	//writeMps(temp);
176	//writeMps("tight");
177	//exit(33);
178	//printf("wrote cc%d\n",iPass);
179	OsiClpSolverInterface::initialSolve();
180	int secondaryStatus = modelPtr_->secondaryStatus();
181	if (modelPtr_->status()==0&&(secondaryStatus==2\|\|secondaryStatus==4))
182	modelPtr_->cleanup(1);
183	//if (!isProvenOptimal())
184	//writeMps("yy");
185	if (isProvenOptimal()&&quadraticModel_&&modelPtr_->numberColumns()==quadraticModel_->numberColumns()) {
186	// see if qp can get better solution
187	const double * solution = modelPtr_->primalColumnSolution();
188	int numberColumns = modelPtr_->numberColumns();
189	bool satisfied=true;
190	for (int i=0;i<numberColumns;i++) {
191	if (isInteger(i)) {
192	double value = solution[i];
193	if (fabs(value-floor(value+0.5))>1.0e-6) {
194	satisfied=false;
195	break;
196	}
197	}
198	}
199	if (satisfied) {
200	ClpSimplex qpTemp(*quadraticModel_);
201	double * lower = qpTemp.columnLower();
202	double * upper = qpTemp.columnUpper();
203	double * lower2 = modelPtr_->columnLower();
204	double * upper2 = modelPtr_->columnUpper();
205	for (int i=0;i<numberColumns;i++) {
206	if (isInteger(i)) {
207	double value = floor(solution[i]+0.5);
208	lower[i]=value;
209	upper[i]=value;
210	} else {
211	lower[i]=lower2[i];
212	upper[i]=upper2[i];
213	}
214	}
215	//qpTemp.writeMps("bad.mps");
216	//modelPtr_->writeMps("bad2.mps");
217	//qpTemp.objectiveAsObject()->setActivated(0);
218	//qpTemp.primal();
219	//qpTemp.objectiveAsObject()->setActivated(1);
220	qpTemp.primal();
221	//assert (!qpTemp.problemStatus());
222	if (qpTemp.objectiveValue()<bestObjectiveValue_-1.0e-3&&!qpTemp.problemStatus()) {
223	delete [] bestSolution_;
224	bestSolution_ = CoinCopyOfArray(qpTemp.primalColumnSolution(),numberColumns);
225	bestObjectiveValue_ = qpTemp.objectiveValue();
226	printf("better qp objective of %g\n",bestObjectiveValue_);
227	// If model has stored then add cut (if convex)
228	if (cbcModel_&&(specialOptions2_&4)!=0) {
229	int numberGenerators = cbcModel_->numberCutGenerators();
230	int iGenerator;
231	cbcModel_->lockThread();
232	for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
233	CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
234	CglCutGenerator * gen = generator->generator();
235	CglStored * gen2 = dynamic_cast<CglStored *> (gen);
236	if (gen2) {
237	// add OA cut
238	double offset;
239	double * gradient = new double [numberColumns+1];
240	memcpy(gradient,qpTemp.objectiveAsObject()->gradient(&qpTemp,bestSolution_,offset,true,2),
241	numberColumns*sizeof(double));
242	// assume convex
243	double rhs = 0.0;
244	int * column = new int[numberColumns+1];
245	int n=0;
246	for (int i=0;i<numberColumns;i++) {
247	double value = gradient[i];
248	if (fabs(value)>1.0e-12) {
249	gradient[n]=value;
250	rhs += value*solution[i];
251	column[n++]=i;
252	}
253	}
254	gradient[n]=-1.0;
255	column[n++]=objectiveVariable_;
256	gen2->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
257	delete [] gradient;
258	delete [] column;
259	break;
260	}
261	}
262	cbcModel_->unlockThread();
263	}
264	}
265	}
266	}
267	}
268	//#define WRITE_MATRIX
269	#ifdef WRITE_MATRIX
270	static int xxxxxx=0;
271	#endif
272	//-----------------------------------------------------------------------------
273	void OsiSolverLink::resolve()
274	{
275	if (false) {
276	bool takeHint;
277	OsiHintStrength strength;
278	// Switch off printing if asked to
279	bool gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength));
280	assert (gotHint);
281	if (strength!=OsiHintIgnore&&takeHint) {
282	printf("no printing\n");
283	} else {
284	printf("printing\n");
285	}
286	}
287	specialOptions_ =0;
288	modelPtr_->setWhatsChanged(0);
289	bool allFixed=numberFix_>0;
290	bool feasible=true;
291	if (numberVariables_) {
292	CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
293	//bool best=true;
294	const double * lower = modelPtr_->columnLower();
295	const double * upper = modelPtr_->columnUpper();
296	// update all bounds before coefficients
297	for (int i=0;i<numberVariables_;i++ ) {
298	info_[i].updateBounds(modelPtr_);
299	int iColumn = info_[i].variable();
300	double lo = lower[iColumn];
301	double up = upper[iColumn];
302	if (up<lo)
303	feasible=false;
304	}
305	int updated=updateCoefficients(modelPtr_,temp);
306	if (updated) {
307	temp->removeGaps(1.0e-14);
308	ClpMatrixBase * save = modelPtr_->clpMatrix();
309	ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
310	assert (clpMatrix);
311	if (save->getNumRows()>temp->getNumRows()) {
312	// add in cuts
313	int numberRows = temp->getNumRows();
314	int * which = new int[numberRows];
315	for (int i=0;i<numberRows;i++)
316	which[i]=i;
317	CoinPackedMatrix * mat = clpMatrix->matrix();
318	// for debug
319	//mat = new CoinPackedMatrix(*mat);
320	mat->deleteRows(numberRows,which);
321	delete [] which;
322	temp->bottomAppendPackedMatrix(*mat);
323	temp->removeGaps(1.0e-14);
324	}
325	if (0) {
326	const CoinPackedMatrix * matrix = modelPtr_->matrix();
327	int numberColumns = matrix->getNumCols();
328	assert (numberColumns==temp->getNumCols());
329	const double * element1 = temp->getMutableElements();
330	const int * row1 = temp->getIndices();
331	const CoinBigIndex * columnStart1 = temp->getVectorStarts();
332	const int * columnLength1 = temp->getVectorLengths();
333	const double * element2 = matrix->getMutableElements();
334	const int * row2 = matrix->getIndices();
335	const CoinBigIndex * columnStart2 = matrix->getVectorStarts();
336	const int * columnLength2 = matrix->getVectorLengths();
337	for (int i=0;i<numberColumns;i++) {
338	assert (columnLength2[i]==columnLength1[i]);
339	int offset = columnStart2[i]-columnStart1[i];
340	for (int j=columnStart1[i];j<columnStart1[i]+columnLength1[i];j++) {
341	assert (row1[j]==row2[j+offset]);
342	assert (element1[j]==element2[j+offset]);
343	}
344	}
345	}
346	modelPtr_->replaceMatrix(temp,true);
347	} else {
348	delete temp;
349	}
350	}
351	#ifdef WRITE_MATRIX
352	{
353	xxxxxx++;
354	char temp[50];
355	sprintf(temp,"bb%d",xxxxxx);
356	writeMps(temp);
357	printf("wrote bb%d\n",xxxxxx);
358	}
359	#endif
360	if (0) {
361	const double * lower = getColLower();
362	const double * upper = getColUpper();
363	int n=0;
364	for (int i=60;i<64;i++) {
365	if (lower[i]) {
366	printf("%d bounds %g %g\n",i,lower[i],upper[i]);
367	n++;
368	}
369	}
370	if (n==1)
371	printf("just one?\n");
372	}
373	// check feasible
374	{
375	const double * lower = getColLower();
376	const double * upper = getColUpper();
377	int numberColumns = getNumCols();
378	for (int i=0;i<numberColumns;i++) {
379	if (lower[i]>upper[i]+1.0e-12) {
380	feasible = false;
381	break;
382	}
383	}
384	}
385	if (!feasible)
386	allFixed=false;
387	if ((specialOptions2_&1)==0)
388	allFixed=false;
389	int returnCode=-1;
390	// See if in strong branching
391	int maxIts = modelPtr_->maximumIterations();
392	if (feasible) {
393	if(maxIts>10000) {
394	// may do lots of work
395	if ((specialOptions2_&1)!=0) {
396	// see if fixed
397	const double * lower = modelPtr_->columnLower();
398	const double * upper = modelPtr_->columnUpper();
399	for (int i=0;i<numberFix_;i++ ) {
400	int iColumn = fixVariables_[i];
401	double lo = lower[iColumn];
402	double up = upper[iColumn];
403	if (up>lo) {
404	allFixed=false;
405	break;
406	}
407	}
408	returnCode=allFixed ? fathom(allFixed) : 0;
409	} else {
410	returnCode=0;
411	}
412	} else {
413	returnCode=0;
414	}
415	}
416	if (returnCode>=0) {
417	if (returnCode==0)
418	OsiClpSolverInterface::resolve();
419	if (!allFixed&&(specialOptions2_&1)!=0) {
420	const double * solution = getColSolution();
421	bool satisfied=true;
422	for (int i=0;i<numberVariables_;i++) {
423	int iColumn = info_[i].variable();
424	double value = solution[iColumn];
425	if (fabs(value-floor(value+0.5))>0.0001)
426	satisfied=false;
427	}
428	//if (satisfied)
429	//printf("satisfied but not fixed\n");
430	}
431	int satisfied=2;
432	const double * solution = getColSolution();
433	const double * lower = getColLower();
434	const double * upper = getColUpper();
435	int numberColumns2 = coinModel_.numberColumns();
436	for (int i=0;i<numberColumns2;i++) {
437	if (isInteger(i)) {
438	double value = solution[i];
439	if (fabs(value-floor(value+0.5))>1.0e-6) {
440	satisfied=0;
441	break;
442	} else if (upper[i]>lower[i]) {
443	satisfied=1;
444	}
445	}
446	}
447	if (isProvenOptimal()) {
448	//if (satisfied==2)
449	//printf("satisfied %d\n",satisfied);
450	if (satisfied&&(specialOptions2_&2)!=0) {
451	assert (quadraticModel_);
452	// look at true objective
453	double direction = modelPtr_->optimizationDirection();
454	assert (direction==1.0);
455	double value = - quadraticModel_->objectiveOffset();
456	const double * objective = quadraticModel_->objective();
457	int i;
458	for ( i=0;i<numberColumns2;i++)
459	value += solution[i]*objective[i];
460	// and now rest
461	for ( i =0;i<numberObjects_;i++) {
462	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
463	if (obj) {
464	value += obj->xyCoefficient(solution);
465	}
466	}
467	if (value<bestObjectiveValue_-1.0e-3) {
468	printf("obj of %g\n",value);
469	//modelPtr_->setDualObjectiveLimit(value);
470	delete [] bestSolution_;
471	bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols());
472	bestObjectiveValue_ = value;
473	if (maxIts<=10000&&cbcModel_) {
474	OsiSolverLink * solver2 = dynamic_cast<OsiSolverLink *> (cbcModel_->solver());
475	assert (solver2);
476	if (solver2!=this) {
477	// in strong branching - need to store in original solver
478	if (value<solver2->bestObjectiveValue_-1.0e-3) {
479	delete [] solver2->bestSolution_;
480	solver2->bestSolution_ = CoinCopyOfArray(bestSolution_,modelPtr_->getNumCols());
481	solver2->bestObjectiveValue_ = value;
482	}
483	}
484	}
485	// If model has stored then add cut (if convex)
486	if (cbcModel_&&(specialOptions2_&4)!=0&&quadraticModel_) {
487	int numberGenerators = cbcModel_->numberCutGenerators();
488	int iGenerator;
489	for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
490	CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
491	CglCutGenerator * gen = generator->generator();
492	CglStored * gen2 = dynamic_cast<CglStored *> (gen);
493	if (gen2) {
494	cbcModel_->lockThread();
495	// add OA cut
496	double offset=0.0;
497	int numberColumns = quadraticModel_->numberColumns();
498	double * gradient = new double [numberColumns+1];
499	// gradient from bilinear
500	int i;
501	CoinZeroN(gradient,numberColumns+1);
502	//const double * objective = modelPtr_->objective();
503	assert (objectiveRow_>=0);
504	const double * element = originalRowCopy_->getElements();
505	const int * column2 = originalRowCopy_->getIndices();
506	const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
507	//const int * rowLength = originalRowCopy_->getVectorLengths();
508	//int numberColumns2 = coinModel_.numberColumns();
509	for ( i=rowStart[objectiveRow_];i<rowStart[objectiveRow_+1];i++)
510	gradient[column2[i]] = element[i];
511	for ( i =0;i<numberObjects_;i++) {
512	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
513	if (obj) {
514	int xColumn = obj->xColumn();
515	int yColumn = obj->yColumn();
516	if (xColumn!=yColumn) {
517	double coefficient = /* 2.0* */obj->coefficient();
518	gradient[xColumn] += coefficient*solution[yColumn];
519	gradient[yColumn] += coefficient*solution[xColumn];
520	offset += coefficientsolution[xColumn]solution[yColumn];
521	} else {
522	double coefficient = obj->coefficient();
523	gradient[xColumn] += 2.0coefficientsolution[yColumn];
524	offset += coefficientsolution[xColumn]solution[yColumn];
525	}
526	}
527	}
528	// assume convex
529	double rhs = 0.0;
530	int * column = new int[numberColumns+1];
531	int n=0;
532	for (int i=0;i<numberColumns;i++) {
533	double value = gradient[i];
534	if (fabs(value)>1.0e-12) {
535	gradient[n]=value;
536	rhs += value*solution[i];
537	column[n++]=i;
538	}
539	}
540	gradient[n]=-1.0;
541	column[n++]=objectiveVariable_;
542	gen2->addCut(-COIN_DBL_MAX,offset+1.0e-4,n,column,gradient);
543	delete [] gradient;
544	delete [] column;
545	cbcModel_->unlockThread();
546	break;
547	}
548	}
549	}
550	}
551	} else if (satisfied==2) {
552	// is there anything left to do?
553	int i;
554	int numberContinuous=0;
555	double gap=0.0;
556	for ( i =0;i<numberObjects_;i++) {
557	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
558	if (obj) {
559	if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
560	numberContinuous++;
561	int xColumn = obj->xColumn();
562	double gapX = upper[xColumn]-lower[xColumn];
563	int yColumn = obj->yColumn();
564	double gapY = upper[yColumn]-lower[yColumn];
565	gap = CoinMax(gap,CoinMax(gapX,gapY));
566	}
567	}
568	}
569	if (numberContinuous&&0) {
570	// iterate to get solution and fathom node
571	int numberColumns2 = coinModel_.numberColumns();
572	double * lower2 = CoinCopyOfArray(getColLower(),numberColumns2);
573	double * upper2 = CoinCopyOfArray(getColUpper(),numberColumns2);
574	while (gap>defaultMeshSize_) {
575	gap *= 0.9;
576	const double * solution = getColSolution();
577	double newGap=0.0;
578	for ( i =0;i<numberObjects_;i++) {
579	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
580	if (obj&&(obj->branchingStrategy()&8)==0) {
581	if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
582	numberContinuous++;
583	// need to make sure new xy value in range
584	double xB[3];
585	double yB[3];
586	double xybar[4];
587	obj->getCoefficients(this,xB,yB,xybar);
588	//double xyTrue = obj->xyCoefficient(solution);
589	double xyLambda=0.0;
590	int firstLambda = obj->firstLambda();
591	for (int j=0;j<4;j++) {
592	xyLambda += solution[firstLambda+j]*xybar[j];
593	}
594	//printf ("x %d, y %d - true %g lambda %g\n",obj->xColumn(),obj->yColumn(),
595	// xyTrue,xyLambda);
596	int xColumn = obj->xColumn();
597	double gapX = upper[xColumn]-lower[xColumn];
598	int yColumn = obj->yColumn();
599	if (gapX>gap) {
600	double value = solution[xColumn];
601	double newLower = CoinMax(lower2[xColumn],value-0.5*gap);
602	double newUpper = CoinMin(upper2[xColumn],value+0.5*gap);
603	if (newUpper-newLower<0.99*gap) {
604	if (newLower==lower2[xColumn])
605	newUpper = CoinMin(upper2[xColumn],newLower+gap);
606	else if (newUpper==upper2[xColumn])
607	newLower = CoinMax(lower2[xColumn],newUpper-gap);
608	}
609	// see if problem
610	double lambda[4];
611	xB[0]=newLower;
612	xB[1]=newUpper;
613	xB[2]=value;
614	yB[2]=solution[yColumn];
615	xybar[0]=xB[0]*yB[0];
616	xybar[1]=xB[0]*yB[1];
617	xybar[2]=xB[1]*yB[0];
618	xybar[3]=xB[1]*yB[1];
619	double infeasibility=obj->computeLambdas(xB,yB,xybar,lambda);
620	assert (infeasibility<1.0e-9);
621	setColLower(xColumn,newLower);
622	setColUpper(xColumn,newUpper);
623	}
624	double gapY = upper[yColumn]-lower[yColumn];
625	if (gapY>gap) {
626	double value = solution[yColumn];
627	double newLower = CoinMax(lower2[yColumn],value-0.5*gap);
628	double newUpper = CoinMin(upper2[yColumn],value+0.5*gap);
629	if (newUpper-newLower<0.99*gap) {
630	if (newLower==lower2[yColumn])
631	newUpper = CoinMin(upper2[yColumn],newLower+gap);
632	else if (newUpper==upper2[yColumn])
633	newLower = CoinMax(lower2[yColumn],newUpper-gap);
634	}
635	// see if problem
636	double lambda[4];
637	yB[0]=newLower;
638	yB[1]=newUpper;
639	xybar[0]=xB[0]*yB[0];
640	xybar[1]=xB[0]*yB[1];
641	xybar[2]=xB[1]*yB[0];
642	xybar[3]=xB[1]*yB[1];
643	double infeasibility=obj->computeLambdas(xB,yB,xybar,lambda);
644	assert (infeasibility<1.0e-9);
645	setColLower(yColumn,newLower);
646	setColUpper(yColumn,newUpper);
647	}
648	newGap = CoinMax(newGap,CoinMax(gapX,gapY));
649	}
650	}
651	}
652	printf("solving with gap of %g\n",gap);
653	//OsiClpSolverInterface::resolve();
654	initialSolve();
655	if (!isProvenOptimal())
656	break;
657	}
658	delete [] lower2;
659	delete [] upper2;
660	//if (isProvenOptimal())
661	//writeMps("zz");
662	}
663	}
664	// ??? - try
665	// But skip if strong branching
666	CbcModel * cbcModel = (modelPtr_->maximumIterations()>10000) ? cbcModel_ : NULL;
667	if ((specialOptions2_&2)!=0) {
668	// If model has stored then add cut (if convex)
669	// off until I work out problem with ibell3a
670	if (cbcModel&&(specialOptions2_&4)!=0&&quadraticModel_) {
671	int numberGenerators = cbcModel_->numberCutGenerators();
672	int iGenerator;
673	for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
674	CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
675	CglCutGenerator * gen = generator->generator();
676	CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
677	if (gen2) {
678	double * solution2 = NULL;
679	int numberColumns = quadraticModel_->numberColumns();
680	int depth = cbcModel_->currentNode() ? cbcModel_->currentNode()->depth() : 0;
681	if (depth<5) {
682	ClpSimplex qpTemp(*quadraticModel_);
683	double * lower = qpTemp.columnLower();
684	double * upper = qpTemp.columnUpper();
685	double * lower2 = modelPtr_->columnLower();
686	double * upper2 = modelPtr_->columnUpper();
687	for (int i=0;i<numberColumns;i++) {
688	lower[i]=lower2[i];
689	upper[i]=upper2[i];
690	}
691	qpTemp.setLogLevel(modelPtr_->logLevel());
692	qpTemp.primal();
693	assert (!qpTemp.problemStatus());
694	if (qpTemp.objectiveValue()<bestObjectiveValue_-1.0e-3&&!qpTemp.problemStatus()) {
695	solution2 = CoinCopyOfArray(qpTemp.primalColumnSolution(),numberColumns);
696	} else {
697	printf("QP says expensive - kill\n");
698	modelPtr_->setProblemStatus(1);
699	modelPtr_->setObjectiveValue(COIN_DBL_MAX);
700	break;
701	}
702	}
703	const double * solution = getColSolution();
704	// add OA cut
705	doAOCuts(gen2, solution, solution);
706	if (solution2) {
707	doAOCuts(gen2, solution, solution2);
708	delete [] solution2;
709	}
710	break;
711	}
712	}
713	}
714	} else if (cbcModel&&(specialOptions2_&8)==8) {
715	// convex and nonlinear in constraints
716	int numberGenerators = cbcModel_->numberCutGenerators();
717	int iGenerator;
718	for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
719	CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
720	CglCutGenerator * gen = generator->generator();
721	CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
722	if (gen2) {
723	const double * solution = getColSolution();
724	const double * rowUpper = getRowUpper();
725	const double * rowLower = getRowLower();
726	const double * element = originalRowCopy_->getElements();
727	const int * column2 = originalRowCopy_->getIndices();
728	const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
729	//const int * rowLength = originalRowCopy_->getVectorLengths();
730	int numberColumns2 = CoinMax(coinModel_.numberColumns(),objectiveVariable_+1);
731	double * gradient = new double [numberColumns2];
732	int * column = new int[numberColumns2];
733	//const double * columnLower = modelPtr_->columnLower();
734	//const double * columnUpper = modelPtr_->columnUpper();
735	cbcModel_->lockThread();
736	for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
737	int iRow = rowNonLinear_[iNon];
738	bool convex = convex_[iNon]>0;
739	if (!convex_[iNon])
740	continue; // can't use this row
741	// add OA cuts
742	double offset=0.0;
743	// gradient from bilinear
744	int i;
745	CoinZeroN(gradient,numberColumns2);
746	//const double * objective = modelPtr_->objective();
747	for ( i=rowStart[iRow];i<rowStart[iRow+1];i++)
748	gradient[column2[i]] = element[i];
749	for ( i =startNonLinear_[iNon];i<startNonLinear_[iNon+1];i++) {
750	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
751	assert (obj);
752	int xColumn = obj->xColumn();
753	int yColumn = obj->yColumn();
754	if (xColumn!=yColumn) {
755	double coefficient = /* 2.0* */obj->coefficient();
756	gradient[xColumn] += coefficient*solution[yColumn];
757	gradient[yColumn] += coefficient*solution[xColumn];
758	offset += coefficientsolution[xColumn]solution[yColumn];
759	} else {
760	double coefficient = obj->coefficient();
761	gradient[xColumn] += 2.0coefficientsolution[yColumn];
762	offset += coefficientsolution[xColumn]solution[yColumn];
763	}
764	}
765	// assume convex
766	double rhs = 0.0;
767	int n=0;
768	for (int i=0;i<numberColumns2;i++) {
769	double value = gradient[i];
770	if (fabs(value)>1.0e-12) {
771	gradient[n]=value;
772	rhs += value*solution[i];
773	column[n++]=i;
774	}
775	}
776	if (iRow==objectiveRow_) {
777	gradient[n]=-1.0;
778	assert (objectiveVariable_>=0);
779	rhs -= solution[objectiveVariable_];
780	column[n++]=objectiveVariable_;
781	assert (convex);
782	} else if (convex) {
783	offset += rowUpper[iRow];
784	} else if (!convex) {
785	offset += rowLower[iRow];
786	}
787	if (convex&&rhs>offset+1.0e-5)
788	gen2->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
789	else if (!convex&&rhs<offset-1.0e-5)
790	gen2->addCut(offset-1.0e-7,COIN_DBL_MAX,n,column,gradient);
791	}
792	cbcModel_->unlockThread();
793	delete [] gradient;
794	delete [] column;
795	break;
796	}
797	}
798	}
799	}
800	} else {
801	modelPtr_->setProblemStatus(1);
802	modelPtr_->setObjectiveValue(COIN_DBL_MAX);
803	}
804	}
805	// Do OA cuts
806	int
807	OsiSolverLink::doAOCuts(CglTemporary * cutGen, const double * solution, const double * solution2)
808	{
809	cbcModel_->lockThread();
810	// add OA cut
811	double offset=0.0;
812	int numberColumns = quadraticModel_->numberColumns();
813	double * gradient = new double [numberColumns+1];
814	// gradient from bilinear
815	int i;
816	CoinZeroN(gradient,numberColumns+1);
817	//const double * objective = modelPtr_->objective();
818	assert (objectiveRow_>=0);
819	const double * element = originalRowCopy_->getElements();
820	const int * column2 = originalRowCopy_->getIndices();
821	const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
822	//const int * rowLength = originalRowCopy_->getVectorLengths();
823	//int numberColumns2 = coinModel_.numberColumns();
824	for ( i=rowStart[objectiveRow_];i<rowStart[objectiveRow_+1];i++)
825	gradient[column2[i]] = element[i];
826	//const double * columnLower = modelPtr_->columnLower();
827	//const double * columnUpper = modelPtr_->columnUpper();
828	for ( i =0;i<numberObjects_;i++) {
829	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
830	if (obj) {
831	int xColumn = obj->xColumn();
832	int yColumn = obj->yColumn();
833	if (xColumn!=yColumn) {
834	double coefficient = /* 2.0* */obj->coefficient();
835	gradient[xColumn] += coefficient*solution2[yColumn];
836	gradient[yColumn] += coefficient*solution2[xColumn];
837	offset += coefficientsolution2[xColumn]solution2[yColumn];
838	} else {
839	double coefficient = obj->coefficient();
840	gradient[xColumn] += 2.0coefficientsolution2[yColumn];
841	offset += coefficientsolution2[xColumn]solution2[yColumn];
842	}
843	}
844	}
845	// assume convex
846	double rhs = 0.0;
847	int * column = new int[numberColumns+1];
848	int n=0;
849	for (int i=0;i<numberColumns;i++) {
850	double value = gradient[i];
851	if (fabs(value)>1.0e-12) {
852	gradient[n]=value;
853	rhs += value*solution[i];
854	column[n++]=i;
855	}
856	}
857	gradient[n]=-1.0;
858	assert (objectiveVariable_>=0);
859	rhs -= solution[objectiveVariable_];
860	column[n++]=objectiveVariable_;
861	int returnCode=0;
862	if (rhs>offset+1.0e-5) {
863	cutGen->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
864	//printf("added cut with %d elements\n",n);
865	returnCode=1;
866	}
867	delete [] gradient;
868	delete [] column;
869	cbcModel_->unlockThread();
870	return returnCode;
871	}
872
873	//#############################################################################
874	// Constructors, destructors clone and assignment
875	//#############################################################################
876
877	//-------------------------------------------------------------------
878	// Default Constructor
879	//-------------------------------------------------------------------
880	OsiSolverLink::OsiSolverLink ()
881	: CbcOsiSolver()
882	{
883	gutsOfDestructor(true);
884	}
885	#if 0
886	/* returns
887	sequence of nonlinear or
888	-1 numeric
889	-2 not found
890	-3 too many terms
891	*/
892	static int getVariable(const CoinModel & model, char * expression,
893	int & linear)
894	{
895	int non=-1;
896	linear=-1;
897	if (strcmp(expression,"Numeric")) {
898	// function
899	char * first = strchr(expression,'*');
900	int numberColumns = model.numberColumns();
901	int j;
902	if (first) {
903	*first='\0';
904	for (j=0;j<numberColumns;j++) {
905	if (!strcmp(expression,model.columnName(j))) {
906	linear=j;
907	memmove(expression,first+1,strlen(first+1)+1);
908	break;
909	}
910	}
911	}
912	// find nonlinear
913	for (j=0;j<numberColumns;j++) {
914	const char * name = model.columnName(j);
915	first = strstr(expression,name);
916	if (first) {
917	if (first!=expression&&isalnum(*(first-1)))
918	continue; // not real match
919	first += strlen(name);
920	if (!isalnum(*first)) {
921	// match
922	non=j;
923	// but check no others
924	j++;
925	for (;j<numberColumns;j++) {
926	const char * name = model.columnName(j);
927	first = strstr(expression,name);
928	if (first) {
929	if (isalnum(*(first-1)))
930	continue; // not real match
931	first += strlen(name);
932	if (!isalnum(*first)) {
933	// match - ouch
934	non=-3;
935	break;
936	}
937	}
938	}
939	break;
940	}
941	}
942	}
943	if (non==-1)
944	non=-2;
945	}
946	return non;
947	}
948	#endif
949	/* This creates from a coinModel object
950
951	if errors.then number of sets is -1
952
953	This creates linked ordered sets information. It assumes -
954
955	for product terms syntax is yy*f(zz)
956	also just f(zz) is allowed
957	and even a constant
958
959	modelObject not const as may be changed as part of process.
960	*/
961	OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
962	: CbcOsiSolver()
963	{
964	gutsOfDestructor(true);
965	load(coinModel);
966	}
967	// need bounds
968	static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue,
969	CoinModel * model1, CoinModel * model2)
970	{
971	double lo = solver->getColLower()[column];
972	if (lo<-maximumValue) {
973	solver->setColLower(column,-maximumValue);
974	if (model1)
975	model1->setColLower(column,-maximumValue);
976	if (model2)
977	model2->setColLower(column,-maximumValue);
978	}
979	double up = solver->getColUpper()[column];
980	if (up>maximumValue) {
981	solver->setColUpper(column,maximumValue);
982	if (model1)
983	model1->setColUpper(column,maximumValue);
984	if (model2)
985	model2->setColUpper(column,maximumValue);
986	}
987	}
988	void OsiSolverLink::load ( CoinModel & coinModelOriginal, bool tightenBounds,int logLevel)
989	{
990	// first check and set up arrays
991	int numberColumns = coinModelOriginal.numberColumns();
992	int numberRows = coinModelOriginal.numberRows();
993	// List of nonlinear entries
994	int * which = new int[numberColumns];
995	numberVariables_=0;
996	//specialOptions2_=0;
997	int iColumn;
998	int numberErrors=0;
999	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1000	CoinModelLink triple=coinModelOriginal.firstInColumn(iColumn);
1001	bool linear=true;
1002	int n=0;
1003	// See if quadratic objective
1004	const char * expr = coinModelOriginal.getColumnObjectiveAsString(iColumn);
1005	if (strcmp(expr,"Numeric")) {
1006	linear=false;
1007	}
1008	while (triple.row()>=0) {
1009	int iRow = triple.row();
1010	const char * expr = coinModelOriginal.getElementAsString(iRow,iColumn);
1011	if (strcmp(expr,"Numeric")) {
1012	linear=false;
1013	}
1014	triple=coinModelOriginal.next(triple);
1015	n++;
1016	}
1017	if (!linear) {
1018	which[numberVariables_++]=iColumn;
1019	}
1020	}
1021	// return if nothing
1022	if (!numberVariables_) {
1023	delete [] which;
1024	coinModel_ = coinModelOriginal;
1025	int nInt=0;
1026	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1027	if (coinModel_.isInteger(iColumn))
1028	nInt++;
1029	}
1030	printf("There are %d integers\n",nInt);
1031	loadFromCoinModel(coinModelOriginal,true);
1032	OsiObject ** objects = new OsiObject * [nInt];
1033	nInt=0;
1034	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1035	if (coinModel_.isInteger(iColumn)) {
1036	objects[nInt] = new OsiSimpleInteger(this,iColumn);
1037	objects[nInt]->setPriority(integerPriority_);
1038	nInt++;
1039	}
1040	}
1041	addObjects(nInt,objects);
1042	int i;
1043	for (i=0;i<nInt;i++)
1044	delete objects[i];
1045	delete [] objects;
1046	return;
1047	} else {
1048	coinModel_ = coinModelOriginal;
1049	// arrays for tightening bounds
1050	int * freeRow = new int [numberRows];
1051	CoinZeroN(freeRow,numberRows);
1052	int * tryColumn = new int [numberColumns];
1053	CoinZeroN(tryColumn,numberColumns);
1054	int nBi=0;
1055	int numberQuadratic=0;
1056	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1057	// See if quadratic objective
1058	const char * expr = coinModel_.getColumnObjectiveAsString(iColumn);
1059	if (strcmp(expr,"Numeric")) {
1060	// check if valuex+-valuey....
1061	assert (strlen(expr)<20000);
1062	tryColumn[iColumn]=1;
1063	char temp[20000];
1064	strcpy(temp,expr);
1065	char * pos = temp;
1066	bool ifFirst=true;
1067	double linearTerm=0.0;
1068	while (*pos) {
1069	double value;
1070	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1071	// must be column unless first when may be linear term
1072	if (jColumn>=0) {
1073	tryColumn[jColumn]=1;
1074	numberQuadratic++;
1075	nBi++;
1076	} else if (jColumn==-2) {
1077	linearTerm = value;
1078	} else {
1079	printf("bad nonlinear term %s\n",temp);
1080	abort();
1081	}
1082	ifFirst=false;
1083	}
1084	coinModelOriginal.setObjective(iColumn,linearTerm);
1085	}
1086	}
1087	int iRow;
1088	int saveNBi=nBi;
1089	for (iRow=0;iRow<numberRows;iRow++) {
1090	CoinModelLink triple=coinModel_.firstInRow(iRow);
1091	while (triple.column()>=0) {
1092	int iColumn = triple.column();
1093	const char * el = coinModel_.getElementAsString(iRow,iColumn);
1094	if (strcmp("Numeric",el)) {
1095	// check if valuex+-valuey....
1096	assert (strlen(el)<20000);
1097	char temp[20000];
1098	strcpy(temp,el);
1099	char * pos = temp;
1100	bool ifFirst=true;
1101	double linearTerm=0.0;
1102	tryColumn[iColumn]=1;
1103	freeRow[iRow]=1;
1104	while (*pos) {
1105	double value;
1106	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1107	// must be column unless first when may be linear term
1108	if (jColumn>=0) {
1109	tryColumn[jColumn]=1;
1110	nBi++;
1111	} else if (jColumn==-2) {
1112	linearTerm = value;
1113	} else {
1114	printf("bad nonlinear term %s\n",temp);
1115	abort();
1116	}
1117	ifFirst=false;
1118	}
1119	coinModelOriginal.setElement(iRow,iColumn,linearTerm);
1120	}
1121	triple=coinModel_.next(triple);
1122	}
1123	}
1124	if (!nBi)
1125	exit(1);
1126	bool quadraticObjective=false;
1127	int nInt=0;
1128	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1129	if (coinModel_.isInteger(iColumn))
1130	nInt++;
1131	}
1132	printf("There are %d bilinear and %d integers\n",nBi,nInt);
1133	loadFromCoinModel(coinModelOriginal,true);
1134	CoinModel coinModel = coinModelOriginal;
1135	if (tightenBounds&&numberColumns<100) {
1136	// first fake bounds
1137	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1138	if (tryColumn[iColumn]) {
1139	fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
1140	}
1141	}
1142	ClpSimplex tempModel(*modelPtr_);
1143	int nDelete = 0;
1144	for (iRow=0;iRow<numberRows;iRow++) {
1145	if (freeRow[iRow])
1146	freeRow[nDelete++]=iRow;
1147	}
1148	tempModel.deleteRows(nDelete,freeRow);
1149	tempModel.setOptimizationDirection(1.0);
1150	if (logLevel<3) {
1151	tempModel.setLogLevel(0);
1152	tempModel.setSpecialOptions(32768);
1153	}
1154	double * objective = tempModel.objective();
1155	CoinZeroN(objective,numberColumns);
1156	// now up and down
1157	double * columnLower = modelPtr_->columnLower();
1158	double * columnUpper = modelPtr_->columnUpper();
1159	const double * solution = tempModel.primalColumnSolution();
1160	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1161	if (tryColumn[iColumn]) {
1162	objective[iColumn]=1.0;
1163	tempModel.primal(1);
1164	if (solution[iColumn]>columnLower[iColumn]+1.0e-3) {
1165	double value = solution[iColumn];
1166	if (coinModel_.isInteger(iColumn))
1167	value = ceil(value-0.9e-3);
1168	if (logLevel>1)
1169	printf("lower bound on %d changed from %g to %g\n",iColumn,columnLower[iColumn],value);
1170	columnLower[iColumn]=value;
1171	coinModel_.setColumnLower(iColumn,value);
1172	coinModel.setColumnLower(iColumn,value);
1173	}
1174	objective[iColumn]=-1.0;
1175	tempModel.primal(1);
1176	if (solution[iColumn]<columnUpper[iColumn]-1.0e-3) {
1177	double value = solution[iColumn];
1178	if (coinModel_.isInteger(iColumn))
1179	value = floor(value+0.9e-3);
1180	if (logLevel>1)
1181	printf("upper bound on %d changed from %g to %g\n",iColumn,columnUpper[iColumn],value);
1182	columnUpper[iColumn]=value;
1183	coinModel_.setColumnUpper(iColumn,value);
1184	coinModel.setColumnUpper(iColumn,value);
1185	}
1186	objective[iColumn]=0.0;
1187	}
1188	}
1189	}
1190	delete [] freeRow;
1191	delete [] tryColumn;
1192	CoinBigIndex * startQuadratic = NULL;
1193	int * columnQuadratic = NULL;
1194	double * elementQuadratic = NULL;
1195	if( saveNBi==nBi) {
1196	printf("all bilinearity in objective\n");
1197	specialOptions2_ \|= 2;
1198	quadraticObjective = true;
1199	// save copy as quadratic model
1200	quadraticModel_ = new ClpSimplex(*modelPtr_);
1201	startQuadratic = new CoinBigIndex [numberColumns+1];
1202	columnQuadratic = new int [numberQuadratic];
1203	elementQuadratic = new double [numberQuadratic];
1204	numberQuadratic=0;
1205	}
1206	//if (quadraticObjective\|\|((specialOptions2_&8)!=0&&saveNBi)) {
1207	if (saveNBi) {
1208	// add in objective as constraint
1209	objectiveVariable_= numberColumns;
1210	objectiveRow_ = coinModel.numberRows();
1211	coinModel.addColumn(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);
1212	int * column = new int[numberColumns+1];
1213	double * element = new double[numberColumns+1];
1214	double * objective = coinModel.objectiveArray();
1215	int n=0;
1216	for (int i=0;i<numberColumns;i++) {
1217	if (objective[i]) {
1218	column[n]=i;
1219	element[n++]=objective[i];
1220	objective[i]=0.0;
1221	}
1222	}
1223	column[n]=objectiveVariable_;
1224	element[n++]=-1.0;
1225	double offset = - coinModel.objectiveOffset();
1226	//assert (!offset); // get sign right if happens
1227	printf("***** offset %g\n",offset);
1228	coinModel.setObjectiveOffset(0.0);
1229	double lowerBound = -COIN_DBL_MAX;
1230	coinModel.addRow(n,column,element,lowerBound,offset);
1231	delete [] column;
1232	delete [] element;
1233	}
1234	OsiObject ** objects = new OsiObject * [nBi+nInt];
1235	char * marked = new char [numberColumns];
1236	memset(marked,0,numberColumns);
1237	// statistics I-I I-x x-x
1238	int stats[3]={0,0,0};
1239	double * sort = new double [nBi];
1240	nBi=nInt;
1241	const OsiObject justBi = const_cast<const OsiObject > (objects+nInt);
1242	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1243	if (quadraticObjective)
1244	startQuadratic[iColumn] = numberQuadratic;
1245	// See if quadratic objective
1246	const char * expr = coinModel_.getColumnObjectiveAsString(iColumn);
1247	if (strcmp(expr,"Numeric")) {
1248	// need bounds
1249	fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
1250	// valuexy
1251	char temp[20000];
1252	strcpy(temp,expr);
1253	char * pos = temp;
1254	bool ifFirst=true;
1255	while (*pos) {
1256	double value;
1257	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1258	// must be column unless first when may be linear term
1259	if (jColumn>=0) {
1260	if (quadraticObjective) {
1261	columnQuadratic[numberQuadratic]=jColumn;
1262	if (jColumn==iColumn)
1263	elementQuadratic[numberQuadratic++]=2.0*value; // convention
1264	else
1265	elementQuadratic[numberQuadratic++]=1.0*value; // convention
1266	}
1267	// need bounds
1268	fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_);
1269	double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
1270	if (meshI)
1271	marked[iColumn]=1;
1272	double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
1273	if (meshJ)
1274	marked[jColumn]=1;
1275	// stats etc
1276	if (meshI) {
1277	if (meshJ)
1278	stats[0]++;
1279	else
1280	stats[1]++;
1281	} else {
1282	if (meshJ)
1283	stats[1]++;
1284	else
1285	stats[2]++;
1286	}
1287	if (iColumn<=jColumn)
1288	sort[nBi-nInt]=iColumn + numberColumns*jColumn;
1289	else
1290	sort[nBi-nInt]=jColumn + numberColumns*iColumn;
1291	if (!meshJ&&!meshI) {
1292	meshI=defaultMeshSize_;
1293	meshJ=0.0;
1294	}
1295	OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,
1296	nBi-nInt,justBi);
1297	newObj->setPriority(biLinearPriority_);
1298	objects[nBi++] = newObj;
1299	} else if (jColumn==-2) {
1300	} else {
1301	printf("bad nonlinear term %s\n",temp);
1302	abort();
1303	}
1304	ifFirst=false;
1305	}
1306	}
1307	}
1308	// stats
1309	printf("There were %d I-I, %d I-x and %d x-x bilinear in objective\n",
1310	stats[0],stats[1],stats[2]);
1311	if (quadraticObjective) {
1312	startQuadratic[numberColumns] = numberQuadratic;
1313	quadraticModel_->loadQuadraticObjective(numberColumns,startQuadratic,
1314	columnQuadratic,elementQuadratic);
1315	delete [] startQuadratic;
1316	delete [] columnQuadratic;
1317	delete [] elementQuadratic;
1318	}
1319	for (iRow=0;iRow<numberRows;iRow++) {
1320	CoinModelLink triple=coinModel_.firstInRow(iRow);
1321	while (triple.column()>=0) {
1322	int iColumn = triple.column();
1323	const char * el = coinModel_.getElementAsString(iRow,iColumn);
1324	if (strcmp("Numeric",el)) {
1325	// need bounds
1326	fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
1327	// valuexy
1328	char temp[20000];
1329	strcpy(temp,el);
1330	char * pos = temp;
1331	bool ifFirst=true;
1332	while (*pos) {
1333	double value;
1334	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1335	// must be column unless first when may be linear term
1336	if (jColumn>=0) {
1337	// need bounds
1338	fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_);
1339	double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
1340	if (meshI)
1341	marked[iColumn]=1;
1342	double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
1343	if (meshJ)
1344	marked[jColumn]=1;
1345	// stats etc
1346	if (meshI) {
1347	if (meshJ)
1348	stats[0]++;
1349	else
1350	stats[1]++;
1351	} else {
1352	if (meshJ)
1353	stats[1]++;
1354	else
1355	stats[2]++;
1356	}
1357	if (iColumn<=jColumn)
1358	sort[nBi-nInt]=iColumn + numberColumns*jColumn;
1359	else
1360	sort[nBi-nInt]=jColumn + numberColumns*iColumn;
1361	if (!meshJ&&!meshI) {
1362	meshI=defaultMeshSize_;
1363	meshJ=0.0;
1364	}
1365	OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,iRow,value,meshI,meshJ,
1366	nBi-nInt,justBi);
1367	newObj->setPriority(biLinearPriority_);
1368	objects[nBi++] = newObj;
1369	} else if (jColumn==-2) {
1370	} else {
1371	printf("bad nonlinear term %s\n",temp);
1372	abort();
1373	}
1374	ifFirst=false;
1375	}
1376	}
1377	triple=coinModel_.next(triple);
1378	}
1379	}
1380	{
1381	// stats
1382	std::sort(sort,sort+nBi-nInt);
1383	int nDiff=0;
1384	double last=-1.0;
1385	for (int i=0;i<nBi-nInt;i++) {
1386	if (sort[i]!=last)
1387	nDiff++;
1388	last=sort[i];
1389	}
1390	delete [] sort;
1391	printf("There were %d I-I, %d I-x and %d x-x bilinear in total of which %d were duplicates\n",
1392	stats[0],stats[1],stats[2],nBi-nInt-nDiff);
1393	}
1394	// reload with all bilinear stuff
1395	loadFromCoinModel(coinModel,true);
1396	//exit(77);
1397	nInt=0;
1398	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1399	if (coinModel_.isInteger(iColumn)) {
1400	objects[nInt] = new OsiSimpleInteger(this,iColumn);
1401	if (marked[iColumn])
1402	objects[nInt]->setPriority(integerPriority_);
1403	else
1404	objects[nInt]->setPriority(integerPriority_);
1405	nInt++;
1406	}
1407	}
1408	nInt=nBi;
1409	delete [] marked;
1410	if (numberErrors) {
1411	// errors
1412	gutsOfDestructor();
1413	numberVariables_=-1;
1414	} else {
1415	addObjects(nInt,objects);
1416	int i;
1417	for (i=0;i<nInt;i++)
1418	delete objects[i];
1419	delete [] objects;
1420	// Now do dummy bound stuff
1421	matrix_ = new CoinPackedMatrix(*getMatrixByCol());
1422	info_ = new OsiLinkedBound [numberVariables_];
1423	for ( i=0;i<numberVariables_;i++) {
1424	info_[i] = OsiLinkedBound(this,which[i],0,NULL,NULL,NULL);
1425	}
1426	// Do row copy but just part
1427	int numberRows2 = objectiveRow_>=0 ? numberRows+1 : numberRows;
1428	int * whichRows = new int [numberRows2];
1429	int * whichColumns = new int [numberColumns];
1430	CoinIotaN(whichRows,numberRows2,0);
1431	CoinIotaN(whichColumns,numberColumns,0);
1432	originalRowCopy_ = new CoinPackedMatrix(*getMatrixByRow(),
1433	numberRows2,whichRows,
1434	numberColumns,whichColumns);
1435	delete [] whichColumns;
1436	numberNonLinearRows_=0;
1437	CoinZeroN(whichRows,numberRows2);
1438	for ( i =0;i<numberObjects_;i++) {
1439	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1440	if (obj) {
1441	int xyRow = obj->xyRow();
1442	assert (xyRow>=0&&xyRow<numberRows2); // even if obj we should move
1443	whichRows[xyRow]++;
1444	}
1445	}
1446	int * pos = new int [numberRows2];
1447	int n=0;
1448	for (i=0;i<numberRows2;i++) {
1449	if (whichRows[i]) {
1450	pos[numberNonLinearRows_]=n;
1451	n+=whichRows[i];
1452	whichRows[i]=numberNonLinearRows_;
1453	numberNonLinearRows_++;
1454	} else {
1455	whichRows[i]=-1;
1456	}
1457	}
1458	startNonLinear_ = new int [numberNonLinearRows_+1];
1459	memcpy(startNonLinear_,pos,numberNonLinearRows_*sizeof(int));
1460	startNonLinear_[numberNonLinearRows_]=n;
1461	rowNonLinear_ = new int [numberNonLinearRows_];
1462	convex_ = new int [numberNonLinearRows_];
1463	// do row numbers now
1464	numberNonLinearRows_=0;
1465	for (i=0;i<numberRows2;i++) {
1466	if (whichRows[i]>=0) {
1467	rowNonLinear_[numberNonLinearRows_++]=i;
1468	}
1469	}
1470	whichNonLinear_ = new int [n];
1471	for ( i =0;i<numberObjects_;i++) {
1472	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1473	if (obj) {
1474	int xyRow = obj->xyRow();
1475	int k=whichRows[xyRow];
1476	int put = pos[k];
1477	pos[k]++;
1478	whichNonLinear_[put]=i;
1479	}
1480	}
1481	delete [] pos;
1482	delete [] whichRows;
1483	analyzeObjects();
1484	}
1485	}
1486	// See if there are any quadratic bounds
1487	int nQ=0;
1488	const CoinPackedMatrix * rowCopy = getMatrixByRow();
1489	//const double * element = rowCopy->getElements();
1490	//const int * column = rowCopy->getIndices();
1491	//const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1492	const int * rowLength = rowCopy->getVectorLengths();
1493	const double * rowLower = getRowLower();
1494	const double * rowUpper = getRowUpper();
1495	for (int iObject =0;iObject<numberObjects_;iObject++) {
1496	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
1497	if (obj) {
1498	int xyRow = obj->xyRow();
1499	if (rowLength[xyRow]==4&&false) {
1500	// we have simple bound
1501	nQ++;
1502	double coefficient = obj->coefficient();
1503	double lo = rowLower[xyRow];
1504	double up = rowUpper[xyRow];
1505	if (coefficient!=1.0) {
1506	printf("*** double check code here\n");
1507	if (coefficient<0.0) {
1508	double temp = lo;
1509	lo = - up;
1510	up = - temp;
1511	coefficient = - coefficient;
1512	}
1513	if (lo>-1.0e20)
1514	lo /= coefficient;
1515	if (up <1.0e20)
1516	up /= coefficient;
1517	setRowLower(xyRow,lo);
1518	setRowUpper(xyRow,up);
1519	// we also need to change elements in matrix_
1520	}
1521	int type=0;
1522	if (lo==up) {
1523	// good news
1524	type=3;
1525	coefficient = lo;
1526	} else if (lo<-1.0e20) {
1527	assert (up<1.0e20);
1528	coefficient = up;
1529	type = 1;
1530	// can we make equality?
1531	} else if (up>1.0e20) {
1532	coefficient = lo;
1533	type = 2;
1534	// can we make equality?
1535	} else {
1536	// we would need extra code
1537	abort();
1538	}
1539	obj->setBoundType(type);
1540	obj->setCoefficient(coefficient);
1541	// can do better if integer?
1542	assert (!isInteger(obj->xColumn()));
1543	assert (!isInteger(obj->yColumn()));
1544	}
1545	}
1546	}
1547	delete [] which;
1548	if ((specialOptions2_&16)!=0)
1549	addTighterConstraints();
1550	}
1551	// Add reformulated bilinear constraints
1552	void
1553	OsiSolverLink::addTighterConstraints()
1554	{
1555	// This is first attempt - for now get working on trimloss
1556	int numberW=0;
1557	int * xW = new int[numberObjects_];
1558	int * yW = new int[numberObjects_];
1559	// Points to firstlambda
1560	int * wW = new int[numberObjects_];
1561	// Coefficient
1562	double * alphaW = new double[numberObjects_];
1563	// Objects
1564	OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_];
1565	int numberColumns = getNumCols();
1566	int firstLambda=numberColumns;
1567	// set up list (better to rethink and do properly as column ordered)
1568	int * list = new int[numberColumns];
1569	memset(list,0,numberColumns*sizeof(int));
1570	int i;
1571	for ( i =0;i<numberObjects_;i++) {
1572	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1573	if (obj) {
1574	//obj->setBranchingStrategy(4); // ***** temp
1575	objW[numberW]=obj;
1576	xW[numberW]=obj->xColumn();
1577	yW[numberW]=obj->yColumn();
1578	list[xW[numberW]]=1;
1579	list[yW[numberW]]=1;
1580	wW[numberW]=obj->firstLambda();
1581	firstLambda = CoinMin(firstLambda,obj->firstLambda());
1582	alphaW[numberW]=obj->coefficient();
1583	//assert (alphaW[numberW]==1.0); // fix when occurs
1584	numberW++;
1585	}
1586	}
1587	int nList = 0;
1588	for (i=0;i<numberColumns;i++) {
1589	if (list[i])
1590	list[nList++]=i;
1591	}
1592	// set up mark array
1593	char * mark = new char [firstLambda*firstLambda];
1594	memset(mark,0,firstLambda*firstLambda);
1595	for (i=0;i<numberW;i++) {
1596	int x = xW[i];
1597	int y = yW[i];
1598	mark[x*firstLambda+y]=1;
1599	mark[y*firstLambda+x]=1;
1600	}
1601	int numberRows2 = originalRowCopy_->getNumRows();
1602	int * addColumn = new int [numberColumns];
1603	double * addElement = new double [numberColumns];
1604	int * addW = new int [numberColumns];
1605	assert (objectiveRow_<0); // fix when occurs
1606	for (int iRow=0;iRow<numberRows2;iRow++) {
1607	for (int iList=0;iList<nList;iList++) {
1608	int kColumn = list[iList];
1609	const double * columnLower = getColLower();
1610	//const double * columnUpper = getColUpper();
1611	const double * rowLower = getRowLower();
1612	const double * rowUpper = getRowUpper();
1613	const CoinPackedMatrix * rowCopy = getMatrixByRow();
1614	const double * element = rowCopy->getElements();
1615	const int * column = rowCopy->getIndices();
1616	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1617	const int * rowLength = rowCopy->getVectorLengths();
1618	CoinBigIndex j;
1619	int numberElements = rowLength[iRow];
1620	int n=0;
1621	for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
1622	int iColumn = column[j];
1623	if (iColumn>=firstLambda) {
1624	// no good
1625	n=-1;
1626	break;
1627	}
1628	if (mark[iColumn*firstLambda+kColumn])
1629	n++;
1630	}
1631	if (n==numberElements) {
1632	printf("can add row %d\n",iRow);
1633	assert (columnLower[kColumn]>=0); // might be able to fix
1634	n=0;
1635	for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
1636	int xColumn=kColumn;
1637	int yColumn = column[j];
1638	int k;
1639	for (k=0;k<numberW;k++) {
1640	if ((xW[k]==yColumn&&yW[k]==xColumn)\|\|
1641	(yW[k]==yColumn&&xW[k]==xColumn))
1642	break;
1643	}
1644	assert (k<numberW);
1645	if (xW[k]!=xColumn) {
1646	int temp=xColumn;
1647	xColumn=yColumn;
1648	yColumn=temp;
1649	}
1650	addW[n/4]=k;
1651	int start = wW[k];
1652	double value = element[j];
1653	for (int kk=0;kk<4;kk++) {
1654	// Dummy value
1655	addElement[n]= value;
1656	addColumn[n++]=start+kk;
1657	}
1658	}
1659	addColumn[n++] = kColumn;
1660	double lo = rowLower[iRow];
1661	double up = rowUpper[iRow];
1662	if (lo>-1.0e20) {
1663	// and tell object
1664	for (j=0;j<n-1;j+=4) {
1665	int iObject = addW[j/4];
1666	objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]);
1667	}
1668	addElement[n-1]=-lo;
1669	if (lo==up)
1670	addRow(n,addColumn,addElement,0.0,0.0);
1671	else
1672	addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
1673	matrix_->appendRow(n,addColumn,addElement);
1674	}
1675	if (up<1.0e20&&up>lo) {
1676	// and tell object
1677	for (j=0;j<n-1;j+=4) {
1678	int iObject = addW[j/4];
1679	objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]);
1680	}
1681	addElement[n-1]=-up;
1682	addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
1683	matrix_->appendRow(n,addColumn,addElement);
1684	}
1685	}
1686	}
1687	}
1688	#if 0
1689	// possibly do bounds
1690	for (int iColumn=0;iColumn<firstLambda;iColumn++) {
1691	for (int iList=0;iList<nList;iList++) {
1692	int kColumn = list[iList];
1693	const double * columnLower = getColLower();
1694	const double * columnUpper = getColUpper();
1695	if (mark[iColumn*firstLambda+kColumn]) {
1696	printf("can add column %d\n",iColumn);
1697	assert (columnLower[kColumn]>=0); // might be able to fix
1698	int xColumn=kColumn;
1699	int yColumn = iColumn;
1700	int k;
1701	for (k=0;k<numberW;k++) {
1702	if ((xW[k]==yColumn&&yW[k]==xColumn)\|\|
1703	(yW[k]==yColumn&&xW[k]==xColumn))
1704	break;
1705	}
1706	assert (k<numberW);
1707	if (xW[k]!=xColumn) {
1708	int temp=xColumn;
1709	xColumn=yColumn;
1710	yColumn=temp;
1711	}
1712	int start = wW[k];
1713	int n=0;
1714	for (int kk=0;kk<4;kk++) {
1715	// Dummy value
1716	addElement[n]= 1.0e-19;
1717	addColumn[n++]=start+kk;
1718	}
1719	// Tell object about this
1720	objW[k]->addExtraRow(matrix_->getNumRows(),1.0);
1721	addColumn[n++] = kColumn;
1722	double lo = columnLower[iColumn];
1723	double up = columnUpper[iColumn];
1724	if (lo>-1.0e20) {
1725	addElement[n-1]=-lo;
1726	if (lo==up)
1727	addRow(n,addColumn,addElement,0.0,0.0);
1728	else
1729	addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
1730	matrix_->appendRow(n,addColumn,addElement);
1731	}
1732	if (up<1.0e20&&up>lo) {
1733	addElement[n-1]=-up;
1734	addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
1735	matrix_->appendRow(n,addColumn,addElement);
1736	}
1737	}
1738	}
1739	}
1740	#endif
1741	delete [] xW;
1742	delete [] yW;
1743	delete [] wW;
1744	delete [] alphaW;
1745	delete [] addColumn;
1746	delete [] addElement;
1747	delete [] addW;
1748	delete [] mark;
1749	delete [] list;
1750	delete [] objW;
1751	}
1752	// Set all biLinear priorities on x-x variables
1753	void
1754	OsiSolverLink::setBiLinearPriorities(int value,double meshSize)
1755	{
1756	OsiObject ** newObject = new OsiObject * [numberObjects_];
1757	int numberOdd=0;
1758	int i;
1759	for ( i =0;i<numberObjects_;i++) {
1760	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1761	if (obj) {
1762	if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
1763	double oldSatisfied = CoinMax(obj->xSatisfied(),
1764	obj->ySatisfied());
1765	OsiBiLinear * objNew = new OsiBiLinear(*obj);
1766	newObject[numberOdd++]=objNew;
1767	objNew->setXSatisfied(0.5*meshSize);
1768	obj->setXOtherSatisfied(0.5*meshSize);
1769	objNew->setXOtherSatisfied(oldSatisfied);
1770	objNew->setXMeshSize(meshSize);
1771	objNew->setYSatisfied(0.5*meshSize);
1772	obj->setYOtherSatisfied(0.5*meshSize);
1773	objNew->setYOtherSatisfied(oldSatisfied);
1774	objNew->setYMeshSize(meshSize);
1775	objNew->setXYSatisfied(0.25*meshSize);
1776	objNew->setPriority(value);
1777	objNew->setBranchingStrategy(8);
1778	}
1779	}
1780	}
1781	addObjects(numberOdd,newObject);
1782	for (i=0;i<numberOdd;i++)
1783	delete newObject[i];
1784	delete [] newObject;
1785	}
1786	/* Set options and priority on all or some biLinear variables
1787	1 - on I-I
1788	2 - on I-x
1789	4 - on x-x
1790	or combinations.
1791	-1 means leave (for priority value and strategy value)
1792	*/
1793	void
1794	OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue,
1795	int mode)
1796	{
1797	int i;
1798	for ( i =0;i<numberObjects_;i++) {
1799	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1800	if (obj) {
1801	bool change=false;
1802	if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0&&(mode&4)!=0)
1803	change=true;
1804	else if (((obj->xMeshSize()==1.0&&obj->yMeshSize()<1.0)\|\|
1805	(obj->xMeshSize()<1.0&&obj->yMeshSize()==1.0))&&(mode&2)!=0)
1806	change=true;
1807	else if (obj->xMeshSize()==1.0&&obj->yMeshSize()==1.0&&(mode&1)!=0)
1808	change=true;
1809	else if (obj->xMeshSize()>1.0\|\|obj->yMeshSize()>1.0)
1810	abort();
1811	if (change) {
1812	if (strategyValue>=0)
1813	obj->setBranchingStrategy(strategyValue);
1814	if (priorityValue>=0)
1815	obj->setPriority(priorityValue);
1816	}
1817	}
1818	}
1819	}
1820
1821	// Say convex (should work it out)
1822	void
1823	OsiSolverLink::sayConvex(bool convex)
1824	{
1825	specialOptions2_ \|= 4;
1826	if (convex_) {
1827	for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
1828	convex_[iNon]=convex ? 1 : -1;
1829	}
1830	}
1831	}
1832	// Set all mesh sizes on x-x variables
1833	void
1834	OsiSolverLink::setMeshSizes(double value)
1835	{
1836	int i;
1837	for ( i =0;i<numberObjects_;i++) {
1838	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1839	if (obj) {
1840	if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
1841	#if 0
1842	numberContinuous++;
1843	int xColumn = obj->xColumn();
1844	double gapX = upper[xColumn]-lower[xColumn];
1845	int yColumn = obj->yColumn();
1846	double gapY = upper[yColumn]-lower[yColumn];
1847	gap = CoinMax(gap,CoinMax(gapX,gapY));
1848	#endif
1849	obj->setMeshSizes(this,value,value);
1850	}
1851	}
1852	}
1853	}
1854	/* Solves nonlinear problem from CoinModel using SLP - may be used as crash
1855	for other algorithms when number of iterations small.
1856	Also exits if all problematical variables are changing
1857	less than deltaTolerance
1858	Returns solution array
1859	*/
1860	double *
1861	OsiSolverLink::nonlinearSLP(int numberPasses,double deltaTolerance)
1862	{
1863	if (!coinModel_.numberRows()) {
1864	printf("Model not set up or nonlinear arrays not created!\n");
1865	return NULL;
1866	}
1867	// first check and set up arrays
1868	int numberColumns = coinModel_.numberColumns();
1869	int numberRows = coinModel_.numberRows();
1870	char * markNonlinear = new char [numberColumns+numberRows];
1871	CoinZeroN(markNonlinear,numberColumns+numberRows);
1872	// List of nonlinear entries
1873	int * listNonLinearColumn = new int[numberColumns];
1874	// List of nonlinear constraints
1875	int * whichRow = new int [numberRows];
1876	CoinZeroN(whichRow,numberRows);
1877	int numberNonLinearColumns=0;
1878	int iColumn;
1879	CoinModel coinModel = coinModel_;
1880	//const CoinModelHash * stringArray = coinModel.stringArray();
1881	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1882	CoinModelLink triple=coinModel.firstInColumn(iColumn);
1883	bool linear=true;
1884	int n=0;
1885	// See if nonlinear objective
1886	const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
1887	if (strcmp(expr,"Numeric")) {
1888	linear=false;
1889	// try and see which columns
1890	assert (strlen(expr)<20000);
1891	char temp[20000];
1892	strcpy(temp,expr);
1893	char * pos = temp;
1894	bool ifFirst=true;
1895	double linearTerm=0.0;
1896	while (*pos) {
1897	double value;
1898	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1899	// must be column unless first when may be linear term
1900	if (jColumn>=0) {
1901	markNonlinear[jColumn]=1;
1902	} else if (jColumn==-2) {
1903	linearTerm = value;
1904	} else {
1905	printf("bad nonlinear term %s\n",temp);
1906	abort();
1907	}
1908	ifFirst=false;
1909	}
1910	}
1911	while (triple.row()>=0) {
1912	int iRow = triple.row();
1913	const char * expr = coinModel.getElementAsString(iRow,iColumn);
1914	if (strcmp(expr,"Numeric")) {
1915	linear=false;
1916	whichRow[iRow]++;
1917	// try and see which columns
1918	assert (strlen(expr)<20000);
1919	char temp[20000];
1920	strcpy(temp,expr);
1921	char * pos = temp;
1922	bool ifFirst=true;
1923	double linearTerm=0.0;
1924	while (*pos) {
1925	double value;
1926	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1927	// must be column unless first when may be linear term
1928	if (jColumn>=0) {
1929	markNonlinear[jColumn]=1;
1930	} else if (jColumn==-2) {
1931	linearTerm = value;
1932	} else {
1933	printf("bad nonlinear term %s\n",temp);
1934	abort();
1935	}
1936	ifFirst=false;
1937	}
1938	}
1939	triple=coinModel.next(triple);
1940	n++;
1941	}
1942	if (!linear) {
1943	markNonlinear[iColumn]=1;
1944	}
1945	}
1946	//int xxxx[]={3,2,0,4,3,0};
1947	//double initialSolution[6];
1948	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1949	if (markNonlinear[iColumn]) {
1950	// put in something
1951	double lower = coinModel.columnLower(iColumn);
1952	double upper = CoinMin(coinModel.columnUpper(iColumn),lower+1000.0);
1953	coinModel.associateElement(coinModel.columnName(iColumn),0.5*(lower+upper));
1954	//coinModel.associateElement(coinModel.columnName(iColumn),xxxx[iColumn]);
1955	listNonLinearColumn[numberNonLinearColumns++]=iColumn;
1956	//initialSolution[iColumn]=xxxx[iColumn];
1957	}
1958	}
1959	// if nothing just solve
1960	if (!numberNonLinearColumns) {
1961	delete [] listNonLinearColumn;
1962	delete [] whichRow;
1963	delete [] markNonlinear;
1964	ClpSimplex tempModel;
1965	tempModel.loadProblem(coinModel,true);
1966	tempModel.initialSolve();
1967	double * solution = CoinCopyOfArray(tempModel.getColSolution(),numberColumns);
1968	return solution;
1969	}
1970	// Create artificials
1971	ClpSimplex tempModel;
1972	tempModel.loadProblem(coinModel,true);
1973	const double * rowLower = tempModel.rowLower();
1974	const double * rowUpper = tempModel.rowUpper();
1975	bool takeAll=false;
1976	int iRow;
1977	int numberArtificials=0;
1978	for (iRow=0;iRow<numberRows;iRow++) {
1979	if (whichRow[iRow]\|\|takeAll) {
1980	if (rowLower[iRow]>-1.0e30)
1981	numberArtificials++;
1982	if (rowUpper[iRow]<1.0e30)
1983	numberArtificials++;
1984	}
1985	}
1986	CoinBigIndex * startArtificial = new CoinBigIndex [numberArtificials+1];
1987	int * rowArtificial = new int [numberArtificials];
1988	double * elementArtificial = new double [numberArtificials];
1989	double * objectiveArtificial = new double [numberArtificials];
1990	numberArtificials=0;
1991	startArtificial[0]=0;
1992	double artificialCost =1.0e9;
1993	for (iRow=0;iRow<numberRows;iRow++) {
1994	if (whichRow[iRow]\|\|takeAll) {
1995	if (rowLower[iRow]>-1.0e30) {
1996	rowArtificial[numberArtificials]=iRow;
1997	elementArtificial[numberArtificials]=1.0;
1998	objectiveArtificial[numberArtificials]=artificialCost;
1999	numberArtificials++;
2000	startArtificial[numberArtificials]=numberArtificials;
2001	}
2002	if (rowUpper[iRow]<1.0e30) {
2003	rowArtificial[numberArtificials]=iRow;
2004	elementArtificial[numberArtificials]=-1.0;
2005	objectiveArtificial[numberArtificials]=artificialCost;
2006	numberArtificials++;
2007	startArtificial[numberArtificials]=numberArtificials;
2008	}
2009	}
2010	}
2011	// Get first solution
2012	int numberColumnsSmall=numberColumns;
2013	ClpSimplex model;
2014	model.loadProblem(coinModel,true);
2015	model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
2016	startArtificial,rowArtificial,elementArtificial);
2017	double * columnLower = model.columnLower();
2018	double * columnUpper = model.columnUpper();
2019	double * trueLower = new double[numberNonLinearColumns];
2020	double * trueUpper = new double[numberNonLinearColumns];
2021	int jNon;
2022	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2023	iColumn=listNonLinearColumn[jNon];
2024	trueLower[jNon]=columnLower[iColumn];
2025	trueUpper[jNon]=columnUpper[iColumn];
2026	//columnLower[iColumn]=initialSolution[iColumn];
2027	//columnUpper[iColumn]=initialSolution[iColumn];
2028	}
2029	model.initialSolve();
2030	//model.writeMps("bad.mps");
2031	// redo number of columns
2032	numberColumns = model.numberColumns();
2033	int * last[3];
2034	double * solution = model.primalColumnSolution();
2035
2036	double * trust = new double[numberNonLinearColumns];
2037	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2038	iColumn=listNonLinearColumn[jNon];
2039	trust[jNon]=0.5;
2040	if (solution[iColumn]<trueLower[jNon])
2041	solution[iColumn]=trueLower[jNon];
2042	else if (solution[iColumn]>trueUpper[jNon])
2043	solution[iColumn]=trueUpper[jNon];
2044	}
2045	int iPass;
2046	double lastObjective=1.0e31;
2047	double * saveSolution = new double [numberColumns];
2048	double * saveRowSolution = new double [numberRows];
2049	memset(saveRowSolution,0,numberRows*sizeof(double));
2050	double * savePi = new double [numberRows];
2051	double * safeSolution = new double [numberColumns];
2052	unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
2053	double targetDrop=1.0e31;
2054	//double objectiveOffset;
2055	//model.getDblParam(ClpObjOffset,objectiveOffset);
2056	// 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
2057	for (iPass=0;iPass<3;iPass++) {
2058	last[iPass]=new int[numberNonLinearColumns];
2059	for (jNon=0;jNon<numberNonLinearColumns;jNon++)
2060	last[iPass][jNon]=0;
2061	}
2062	// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2063	int goodMove=-2;
2064	char * statusCheck = new char[numberColumns];
2065	double * changeRegion = new double [numberColumns];
2066	int logLevel=63;
2067	double dualTolerance = model.dualTolerance();
2068	double primalTolerance = model.primalTolerance();
2069	int lastGoodMove=1;
2070	for (iPass=0;iPass<numberPasses;iPass++) {
2071	lastGoodMove=goodMove;
2072	columnLower = model.columnLower();
2073	columnUpper = model.columnUpper();
2074	solution = model.primalColumnSolution();
2075	double * rowActivity = model.primalRowSolution();
2076	// redo objective
2077	ClpSimplex tempModel;
2078	// load new values
2079	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2080	iColumn=listNonLinearColumn[jNon];
2081	coinModel.associateElement(coinModel.columnName(iColumn),solution[iColumn]);
2082	}
2083	tempModel.loadProblem(coinModel);
2084	double objectiveOffset;
2085	tempModel.getDblParam(ClpObjOffset,objectiveOffset);
2086	double objValue=-objectiveOffset;
2087	const double * objective = tempModel.objective();
2088	for (iColumn=0;iColumn<numberColumnsSmall;iColumn++)
2089	objValue += solution[iColumn]*objective[iColumn];
2090	double * rowActivity2 = tempModel.primalRowSolution();
2091	const double * rowLower2 = tempModel.rowLower();
2092	const double * rowUpper2 = tempModel.rowUpper();
2093	memset(rowActivity2,0,numberRows*sizeof(double));
2094	tempModel.times(1.0,solution,rowActivity2);
2095	for (iRow=0;iRow<numberRows;iRow++) {
2096	if (rowActivity2[iRow]<rowLower2[iRow]-primalTolerance)
2097	objValue += (rowLower2[iRow]-rowActivity2[iRow]-primalTolerance)*artificialCost;
2098	else if (rowActivity2[iRow]>rowUpper2[iRow]+primalTolerance)
2099	objValue -= (rowUpper2[iRow]-rowActivity2[iRow]+primalTolerance)*artificialCost;
2100	}
2101	double theta=-1.0;
2102	double maxTheta=COIN_DBL_MAX;
2103	if (objValue<=lastObjective+1.0e-15*fabs(lastObjective)\|\|!iPass)
2104	goodMove=1;
2105	else
2106	goodMove=-1;
2107	//maxTheta=1.0;
2108	if (iPass) {
2109	int jNon=0;
2110	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2111	changeRegion[iColumn]=solution[iColumn]-saveSolution[iColumn];
2112	double alpha = changeRegion[iColumn];
2113	double oldValue = saveSolution[iColumn];
2114	if (markNonlinear[iColumn]==0) {
2115	// linear
2116	if (alpha<-1.0e-15) {
2117	// variable going towards lower bound
2118	double bound = columnLower[iColumn];
2119	oldValue -= bound;
2120	if (oldValue+maxTheta*alpha<0.0) {
2121	maxTheta = CoinMax(0.0,oldValue/(-alpha));
2122	}
2123	} else if (alpha>1.0e-15) {
2124	// variable going towards upper bound
2125	double bound = columnUpper[iColumn];
2126	oldValue = bound-oldValue;
2127	if (oldValue-maxTheta*alpha<0.0) {
2128	maxTheta = CoinMax(0.0,oldValue/alpha);
2129	}
2130	}
2131	} else {
2132	// nonlinear
2133	if (alpha<-1.0e-15) {
2134	// variable going towards lower bound
2135	double bound = trueLower[jNon];
2136	oldValue -= bound;
2137	if (oldValue+maxTheta*alpha<0.0) {
2138	maxTheta = CoinMax(0.0,oldValue/(-alpha));
2139	}
2140	} else if (alpha>1.0e-15) {
2141	// variable going towards upper bound
2142	double bound = trueUpper[jNon];
2143	oldValue = bound-oldValue;
2144	if (oldValue-maxTheta*alpha<0.0) {
2145	maxTheta = CoinMax(0.0,oldValue/alpha);
2146	}
2147	}
2148	jNon++;
2149	}
2150	}
2151	// make sure both accurate
2152	memset(rowActivity,0,numberRows*sizeof(double));
2153	model.times(1.0,solution,rowActivity);
2154	memset(saveRowSolution,0,numberRows*sizeof(double));
2155	model.times(1.0,saveSolution,saveRowSolution);
2156	for (int iRow=0;iRow<numberRows;iRow++) {
2157	double alpha =rowActivity[iRow]-saveRowSolution[iRow];
2158	double oldValue = saveRowSolution[iRow];
2159	if (alpha<-1.0e-15) {
2160	// variable going towards lower bound
2161	double bound = rowLower[iRow];
2162	oldValue -= bound;
2163	if (oldValue+maxTheta*alpha<0.0) {
2164	maxTheta = CoinMax(0.0,oldValue/(-alpha));
2165	}
2166	} else if (alpha>1.0e-15) {
2167	// variable going towards upper bound
2168	double bound = rowUpper[iRow];
2169	oldValue = bound-oldValue;
2170	if (oldValue-maxTheta*alpha<0.0) {
2171	maxTheta = CoinMax(0.0,oldValue/alpha);
2172	}
2173	}
2174	}
2175	} else {
2176	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2177	changeRegion[iColumn]=0.0;
2178	saveSolution[iColumn]=solution[iColumn];
2179	}
2180	memcpy(saveRowSolution,rowActivity,numberRows*sizeof(double));
2181	}
2182	if (goodMove>=0) {
2183	//theta = CoinMin(theta2,maxTheta);
2184	theta = maxTheta;
2185	if (theta>0.0&&theta<=1.0) {
2186	// update solution
2187	double lambda = 1.0-theta;
2188	for (iColumn=0;iColumn<numberColumns;iColumn++)
2189	solution[iColumn] = lambda * saveSolution[iColumn]
2190	+ theta * solution[iColumn];
2191	memset(rowActivity,0,numberRows*sizeof(double));
2192	model.times(1.0,solution,rowActivity);
2193	if (lambda>0.999) {
2194	memcpy(model.dualRowSolution(),savePi,numberRows*sizeof(double));
2195	memcpy(model.statusArray(),saveStatus,numberRows+numberColumns);
2196	}
2197	// redo rowActivity
2198	memset(rowActivity,0,numberRows*sizeof(double));
2199	model.times(1.0,solution,rowActivity);
2200	}
2201	}
2202	// load new values
2203	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2204	iColumn=listNonLinearColumn[jNon];
2205	coinModel.associateElement(coinModel.columnName(iColumn),solution[iColumn]);
2206	}
2207	double * sol2 = CoinCopyOfArray(model.primalColumnSolution(),numberColumns);
2208	unsigned char * status2 = CoinCopyOfArray(model.statusArray(),numberColumns);
2209	model.loadProblem(coinModel);
2210	model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
2211	startArtificial,rowArtificial,elementArtificial);
2212	memcpy(model.primalColumnSolution(),sol2,numberColumns*sizeof(double));
2213	memcpy(model.statusArray(),status2,numberColumns);
2214	delete [] sol2;
2215	delete [] status2;
2216	columnLower = model.columnLower();
2217	columnUpper = model.columnUpper();
2218	solution = model.primalColumnSolution();
2219	rowActivity = model.primalRowSolution();
2220	int * temp=last[2];
2221	last[2]=last[1];
2222	last[1]=last[0];
2223	last[0]=temp;
2224	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2225	iColumn=listNonLinearColumn[jNon];
2226	double change = solution[iColumn]-saveSolution[iColumn];
2227	if (change<-1.0e-5) {
2228	if (fabs(change+trust[jNon])<1.0e-5)
2229	temp[jNon]=-1;
2230	else
2231	temp[jNon]=-2;
2232	} else if(change>1.0e-5) {
2233	if (fabs(change-trust[jNon])<1.0e-5)
2234	temp[jNon]=1;
2235	else
2236	temp[jNon]=2;
2237	} else {
2238	temp[jNon]=0;
2239	}
2240	}
2241	// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2242	double maxDelta=0.0;
2243	if (goodMove>=0) {
2244	if (objValue<=lastObjective+1.0e-15*fabs(lastObjective))
2245	goodMove=1;
2246	else
2247	goodMove=0;
2248	} else {
2249	maxDelta=1.0e10;
2250	}
2251	double maxGap=0.0;
2252	int numberSmaller=0;
2253	int numberSmaller2=0;
2254	int numberLarger=0;
2255	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2256	iColumn=listNonLinearColumn[jNon];
2257	maxDelta = CoinMax(maxDelta,
2258	fabs(solution[iColumn]-saveSolution[iColumn]));
2259	if (goodMove>0) {
2260	if (last[0][jNon]*last[1][jNon]<0) {
2261	// halve
2262	trust[jNon] *= 0.5;
2263	numberSmaller2++;
2264	} else {
2265	if (last[0][jNon]==last[1][jNon]&&
2266	last[0][jNon]==last[2][jNon])
2267	trust[jNon] = CoinMin(1.5*trust[jNon],1.0e6);
2268	numberLarger++;
2269	}
2270	} else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
2271	trust[jNon] *= 0.2;
2272	numberSmaller++;
2273	}
2274	maxGap = CoinMax(maxGap,trust[jNon]);
2275	}
2276	#ifdef CLP_DEBUG
2277	if (logLevel&32)
2278	std::cout<<"largest gap is "<<maxGap<<" "
2279	<<numberSmaller+numberSmaller2<<" reduced ("
2280	<<numberSmaller<<" badMove ), "
2281	<<numberLarger<<" increased"<<std::endl;
2282	#endif
2283	if (iPass>10000) {
2284	for (jNon=0;jNon<numberNonLinearColumns;jNon++)
2285	trust[jNon] *=0.0001;
2286	}
2287	printf("last good %d goodMove %d\n",lastGoodMove,goodMove);
2288	if (goodMove>0) {
2289	double drop = lastObjective-objValue;
2290	printf("Pass %d, objective %g - drop %g maxDelta %g\n",iPass,objValue,drop,maxDelta);
2291	if (iPass>20&&drop<1.0e-12*fabs(objValue)&&lastGoodMove>0)
2292	drop=0.999e-4; // so will exit
2293	if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&theta<0.99999&&lastGoodMove>0) {
2294	if (logLevel>1)
2295	std::cout<<"Exiting as maxDelta < tolerance and small drop"<<std::endl;
2296	break;
2297	}
2298	} else if (!numberSmaller&&iPass>1) {
2299	if (logLevel>1)
2300	std::cout<<"Exiting as all gaps small"<<std::endl;
2301	break;
2302	}
2303	if (!iPass)
2304	goodMove=1;
2305	targetDrop=0.0;
2306	double * r = model.dualColumnSolution();
2307	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2308	iColumn=listNonLinearColumn[jNon];
2309	columnLower[iColumn]=CoinMax(solution[iColumn]
2310	-trust[jNon],
2311	trueLower[jNon]);
2312	columnUpper[iColumn]=CoinMin(solution[iColumn]
2313	+trust[jNon],
2314	trueUpper[jNon]);
2315	}
2316	if (iPass) {
2317	// get reduced costs
2318	model.matrix()->transposeTimes(savePi,
2319	model.dualColumnSolution());
2320	const double * objective = model.objective();
2321	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2322	iColumn=listNonLinearColumn[jNon];
2323	double dj = objective[iColumn]-r[iColumn];
2324	r[iColumn]=dj;
2325	if (dj<-dualTolerance)
2326	targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
2327	else if (dj>dualTolerance)
2328	targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
2329	}
2330	} else {
2331	memset(r,0,numberColumns*sizeof(double));
2332	}
2333	#if 0
2334	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2335	iColumn=listNonLinearColumn[jNon];
2336	if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
2337	columnLower[iColumn]=CoinMax(solution[iColumn],
2338	trueLower[jNon]);
2339	columnUpper[iColumn]=CoinMin(solution[iColumn]
2340	+trust[jNon],
2341	trueUpper[jNon]);
2342	} else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
2343	columnLower[iColumn]=CoinMax(solution[iColumn]
2344	-trust[jNon],
2345	trueLower[jNon]);
2346	columnUpper[iColumn]=CoinMin(solution[iColumn],
2347	trueUpper[jNon]);
2348	} else {
2349	columnLower[iColumn]=CoinMax(solution[iColumn]
2350	-trust[jNon],
2351	trueLower[jNon]);
2352	columnUpper[iColumn]=CoinMin(solution[iColumn]
2353	+trust[jNon],
2354	trueUpper[jNon]);
2355	}
2356	}
2357	#endif
2358	if (goodMove>0) {
2359	memcpy(saveSolution,solution,numberColumns*sizeof(double));
2360	memcpy(saveRowSolution,rowActivity,numberRows*sizeof(double));
2361	memcpy(savePi,model.dualRowSolution(),numberRows*sizeof(double));
2362	memcpy(saveStatus,model.statusArray(),numberRows+numberColumns);
2363
2364	#ifdef CLP_DEBUG
2365	if (logLevel&32)
2366	std::cout<<"Pass - "<<iPass
2367	<<", target drop is "<<targetDrop
2368	<<std::endl;
2369	#endif
2370	lastObjective = objValue;
2371	if (targetDrop<CoinMax(1.0e-8,CoinMin(1.0e-6,1.0e-6*fabs(objValue)))&&lastGoodMove&&iPass>3) {
2372	if (logLevel>1)
2373	printf("Exiting on target drop %g\n",targetDrop);
2374	break;
2375	}
2376	#ifdef CLP_DEBUG
2377	{
2378	double * r = model.dualColumnSolution();
2379	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2380	iColumn=listNonLinearColumn[jNon];
2381	if (logLevel&32)
2382	printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
2383	jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
2384	r[iColumn],statusCheck[iColumn],columnLower[iColumn],
2385	columnUpper[iColumn]);
2386	}
2387	}
2388	#endif
2389	model.scaling(false);
2390	model.primal(1);
2391	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2392	iColumn=listNonLinearColumn[jNon];
2393	printf("%d bounds etc %g %g %g\n",iColumn, columnLower[iColumn],solution[iColumn],columnUpper[iColumn]);
2394	}
2395	char temp[20];
2396	sprintf(temp,"pass%d.mps",iPass);
2397	//model.writeMps(temp);
2398	#ifdef CLP_DEBUG
2399	if (model.status()) {
2400	model.writeMps("xx.mps");
2401	}
2402	#endif
2403	if (model.status()==1) {
2404	// not feasible ! - backtrack and exit
2405	// use safe solution
2406	memcpy(solution,safeSolution,numberColumns*sizeof(double));
2407	memcpy(saveSolution,solution,numberColumns*sizeof(double));
2408	memset(rowActivity,0,numberRows*sizeof(double));
2409	model.times(1.0,solution,rowActivity);
2410	memcpy(saveRowSolution,rowActivity,numberRows*sizeof(double));
2411	memcpy(model.dualRowSolution(),savePi,numberRows*sizeof(double));
2412	memcpy(model.statusArray(),saveStatus,numberRows+numberColumns);
2413	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2414	iColumn=listNonLinearColumn[jNon];
2415	columnLower[iColumn]=CoinMax(solution[iColumn]
2416	-trust[jNon],
2417	trueLower[jNon]);
2418	columnUpper[iColumn]=CoinMin(solution[iColumn]
2419	+trust[jNon],
2420	trueUpper[jNon]);
2421	}
2422	break;
2423	} else {
2424	// save in case problems
2425	memcpy(safeSolution,solution,numberColumns*sizeof(double));
2426	}
2427	goodMove=1;
2428	} else {
2429	// bad pass - restore solution
2430	#ifdef CLP_DEBUG
2431	if (logLevel&32)
2432	printf("Backtracking\n");
2433	#endif
2434	// load old values
2435	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2436	iColumn=listNonLinearColumn[jNon];
2437	coinModel.associateElement(coinModel.columnName(iColumn),saveSolution[iColumn]);
2438	}
2439	model.loadProblem(coinModel);
2440	model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
2441	startArtificial,rowArtificial,elementArtificial);
2442	solution = model.primalColumnSolution();
2443	rowActivity = model.primalRowSolution();
2444	memcpy(solution,saveSolution,numberColumns*sizeof(double));
2445	memcpy(rowActivity,saveRowSolution,numberRows*sizeof(double));
2446	memcpy(model.dualRowSolution(),savePi,numberRows*sizeof(double));
2447	memcpy(model.statusArray(),saveStatus,numberRows+numberColumns);
2448	columnLower = model.columnLower();
2449	columnUpper = model.columnUpper();
2450	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2451	iColumn=listNonLinearColumn[jNon];
2452	columnLower[iColumn]=solution[iColumn];
2453	columnUpper[iColumn]=solution[iColumn];
2454	}
2455	model.primal(1);
2456	//model.writeMps("xx.mps");
2457	iPass--;
2458	goodMove=-1;
2459	}
2460	}
2461	// restore solution
2462	memcpy(solution,saveSolution,numberColumns*sizeof(double));
2463	delete [] statusCheck;
2464	delete [] savePi;
2465	delete [] saveStatus;
2466	// load new values
2467	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2468	iColumn=listNonLinearColumn[jNon];
2469	coinModel.associateElement(coinModel.columnName(iColumn),solution[iColumn]);
2470	}
2471	double * sol2 = CoinCopyOfArray(model.primalColumnSolution(),numberColumns);
2472	unsigned char * status2 = CoinCopyOfArray(model.statusArray(),numberColumns);
2473	model.loadProblem(coinModel);
2474	model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
2475	startArtificial,rowArtificial,elementArtificial);
2476	memcpy(model.primalColumnSolution(),sol2,numberColumns*sizeof(double));
2477	memcpy(model.statusArray(),status2,numberColumns);
2478	delete [] sol2;
2479	delete [] status2;
2480	columnLower = model.columnLower();
2481	columnUpper = model.columnUpper();
2482	solution = model.primalColumnSolution();
2483	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2484	iColumn=listNonLinearColumn[jNon];
2485	columnLower[iColumn]=CoinMax(solution[iColumn],
2486	trueLower[jNon]);
2487	columnUpper[iColumn]=CoinMin(solution[iColumn],
2488	trueUpper[jNon]);
2489	}
2490	model.primal(1);
2491	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2492	iColumn=listNonLinearColumn[jNon];
2493	columnLower[iColumn]= trueLower[jNon];
2494	columnUpper[iColumn]= trueUpper[jNon];
2495	}
2496	delete [] saveSolution;
2497	delete [] safeSolution;
2498	delete [] saveRowSolution;
2499	for (iPass=0;iPass<3;iPass++)
2500	delete [] last[iPass];
2501	delete [] trust;
2502	delete [] trueUpper;
2503	delete [] trueLower;
2504	delete [] changeRegion;
2505	delete [] startArtificial;
2506	delete [] rowArtificial;
2507	delete [] elementArtificial;
2508	delete [] objectiveArtificial;
2509	delete [] listNonLinearColumn;
2510	delete [] whichRow;
2511	delete [] markNonlinear;
2512	return CoinCopyOfArray(solution,coinModel.numberColumns());
2513	}
2514	/* Solve linearized quadratic objective branch and bound.
2515	Return cutoff and OA cut
2516	*/
2517	double
2518	OsiSolverLink::linearizedBAB(CglStored * cut)
2519	{
2520	double bestObjectiveValue=COIN_DBL_MAX;
2521	if (quadraticModel_) {
2522	ClpSimplex * qp = new ClpSimplex(*quadraticModel_);
2523	// bounds
2524	int numberColumns = qp->numberColumns();
2525	double * lower = qp->columnLower();
2526	double * upper = qp->columnUpper();
2527	const double * lower2 = getColLower();
2528	const double * upper2 = getColUpper();
2529	for (int i=0;i<numberColumns;i++) {
2530	lower[i] = CoinMax(lower[i],lower2[i]);
2531	upper[i] = CoinMin(upper[i],upper2[i]);
2532	}
2533	qp->nonlinearSLP(20,1.0e-5);
2534	qp->primal();
2535	OsiSolverLinearizedQuadratic solver2(qp);
2536	const double * solution=NULL;
2537	// Reduce printout
2538	solver2.setHintParam(OsiDoReducePrint,true,OsiHintTry);
2539	CbcModel model2(solver2);
2540	// Now do requested saves and modifications
2541	CbcModel * cbcModel = & model2;
2542	OsiSolverInterface * osiModel = model2.solver();
2543	OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
2544	ClpSimplex * clpModel = osiclpModel->getModelPtr();
2545
2546	// Set changed values
2547
2548	CglProbing probing;
2549	probing.setMaxProbe(10);
2550	probing.setMaxLook(10);
2551	probing.setMaxElements(200);
2552	probing.setMaxProbeRoot(50);
2553	probing.setMaxLookRoot(10);
2554	probing.setRowCuts(3);
2555	probing.setUsingObjective(true);
2556	cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
2557	cbcModel->cutGenerator(0)->setTiming(true);
2558
2559	CglGomory gomory;
2560	gomory.setLimitAtRoot(512);
2561	cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
2562	cbcModel->cutGenerator(1)->setTiming(true);
2563
2564	CglKnapsackCover knapsackCover;
2565	cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
2566	cbcModel->cutGenerator(2)->setTiming(true);
2567
2568	CglClique clique;
2569	clique.setStarCliqueReport(false);
2570	clique.setRowCliqueReport(false);
2571	clique.setMinViolation(0.1);
2572	cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
2573	cbcModel->cutGenerator(3)->setTiming(true);
2574
2575	CglMixedIntegerRounding2 mixedIntegerRounding2;
2576	cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
2577	cbcModel->cutGenerator(4)->setTiming(true);
2578
2579	CglFlowCover flowCover;
2580	cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
2581	cbcModel->cutGenerator(5)->setTiming(true);
2582
2583	CglTwomir twomir;
2584	twomir.setMaxElements(250);
2585	cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
2586	cbcModel->cutGenerator(6)->setTiming(true);
2587	// For now - switch off most heuristics (because CglPreProcess is bad with QP)
2588	#if 1
2589	CbcHeuristicFPump heuristicFPump(*cbcModel);
2590	heuristicFPump.setWhen(13);
2591	heuristicFPump.setMaximumPasses(20);
2592	heuristicFPump.setMaximumRetries(7);
2593	heuristicFPump.setAbsoluteIncrement(4332.64);
2594	cbcModel->addHeuristic(&heuristicFPump);
2595	heuristicFPump.setInitialWeight(1);
2596
2597	CbcHeuristicLocal heuristicLocal(*cbcModel);
2598	heuristicLocal.setSearchType(1);
2599	cbcModel->addHeuristic(&heuristicLocal);
2600
2601	CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
2602	cbcModel->addHeuristic(&heuristicGreedyCover);
2603
2604	CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
2605	cbcModel->addHeuristic(&heuristicGreedyEquality);
2606	#endif
2607
2608	CbcRounding rounding(*cbcModel);
2609	rounding.setHeuristicName("rounding");
2610	cbcModel->addHeuristic(&rounding);
2611
2612	cbcModel->setNumberBeforeTrust(5);
2613	cbcModel->setSpecialOptions(2);
2614	cbcModel->messageHandler()->setLogLevel(1);
2615	cbcModel->setMaximumCutPassesAtRoot(-100);
2616	cbcModel->setMaximumCutPasses(1);
2617	cbcModel->setMinimumDrop(0.05);
2618	// For branchAndBound this may help
2619	clpModel->defaultFactorizationFrequency();
2620	clpModel->setDualBound(1.0001e+08);
2621	clpModel->setPerturbation(50);
2622	osiclpModel->setSpecialOptions(193);
2623	osiclpModel->messageHandler()->setLogLevel(0);
2624	osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
2625	osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
2626	// You can save some time by switching off message building
2627	// clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
2628
2629	// Solve
2630
2631	cbcModel->initialSolve();
2632	if (clpModel->tightenPrimalBounds()!=0) {
2633	std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
2634	delete qp;
2635	return COIN_DBL_MAX;
2636	}
2637	clpModel->dual(); // clean up
2638	cbcModel->initialSolve();
2639	cbcModel->branchAndBound();
2640	OsiSolverLinearizedQuadratic * solver3 = dynamic_cast<OsiSolverLinearizedQuadratic *> (model2.solver());
2641	assert (solver3);
2642	solution = solver3->bestSolution();
2643	bestObjectiveValue = solver3->bestObjectiveValue();
2644	setBestObjectiveValue(bestObjectiveValue);
2645	setBestSolution(solution,solver3->getNumCols());
2646	// if convex
2647	if ((specialOptions2()&4)!=0) {
2648	if (cbcModel_)
2649	cbcModel_->lockThread();
2650	// add OA cut
2651	double offset;
2652	double * gradient = new double [numberColumns+1];
2653	memcpy(gradient,qp->objectiveAsObject()->gradient(qp,solution,offset,true,2),
2654	numberColumns*sizeof(double));
2655	double rhs = 0.0;
2656	int * column = new int[numberColumns+1];
2657	int n=0;
2658	for (int i=0;i<numberColumns;i++) {
2659	double value = gradient[i];
2660	if (fabs(value)>1.0e-12) {
2661	gradient[n]=value;
2662	rhs += value*solution[i];
2663	column[n++]=i;
2664	}
2665	}
2666	gradient[n]=-1.0;
2667	column[n++]=numberColumns;
2668	cut->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
2669	delete [] gradient;
2670	delete [] column;
2671	if (cbcModel_)
2672	cbcModel_->unlockThread();
2673	}
2674	delete qp;
2675	printf("obj %g\n",bestObjectiveValue);
2676	}
2677	return bestObjectiveValue;
2678	}
2679	/* Solves nonlinear problem from CoinModel using SLP - and then tries to get
2680	heuristic solution
2681	Returns solution array
2682	*/
2683	double *
2684	OsiSolverLink::heuristicSolution(int numberPasses,double deltaTolerance,int mode)
2685	{
2686	// get a solution
2687	CoinModel tempModel = coinModel_;
2688	ClpSimplex * temp = approximateSolution(tempModel, numberPasses, deltaTolerance);
2689	int numberColumns = coinModel_.numberColumns();
2690	double * solution = CoinCopyOfArray(temp->primalColumnSolution(),numberColumns);
2691	delete temp;
2692	if (mode==0) {
2693	return solution;
2694	} else if (mode==2) {
2695	const double * lower = getColLower();
2696	const double * upper = getColUpper();
2697	for (int iObject =0;iObject<numberObjects_;iObject++) {
2698	OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
2699	if (obj&&(obj->priority()<biLinearPriority_\|\|biLinearPriority_<=0)) {
2700	int iColumn = obj->columnNumber();
2701	double value = solution[iColumn];
2702	value = floor(value+0.5);
2703	if (fabs(value-solution[iColumn])>0.01) {
2704	setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0)));
2705	setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,1.0)));
2706	} else {
2707	// could fix to integer
2708	setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0)));
2709	setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,0.0)));
2710	}
2711	}
2712	}
2713	return solution;
2714	}
2715	OsiClpSolverInterface newSolver;
2716	if (mode==1) {
2717	// round all with priority < biLinearPriority_
2718	setFixedPriority(biLinearPriority_);
2719	// ? should we save and restore coin model
2720	tempModel = coinModel_;
2721	// solve modified problem
2722	char * mark = new char[numberColumns];
2723	memset(mark,0,numberColumns);
2724	for (int iObject =0;iObject<numberObjects_;iObject++) {
2725	OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
2726	if (obj&&obj->priority()<biLinearPriority_) {
2727	int iColumn = obj->columnNumber();
2728	double value = solution[iColumn];
2729	value = ceil(value-1.0e-7);
2730	tempModel.associateElement(coinModel_.columnName(iColumn),value);
2731	mark[iColumn]=1;
2732	}
2733	OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
2734	if (objB) {
2735	// if one or both continuous then fix one
2736	if (objB->xMeshSize()<1.0) {
2737	int xColumn = objB->xColumn();
2738	double value = solution[xColumn];
2739	tempModel.associateElement(coinModel_.columnName(xColumn),value);
2740	mark[xColumn]=1;
2741	} else if (objB->yMeshSize()<1.0) {
2742	int yColumn = objB->yColumn();
2743	double value = solution[yColumn];
2744	tempModel.associateElement(coinModel_.columnName(yColumn),value);
2745	mark[yColumn]=1;
2746	}
2747	}
2748	}
2749	CoinModel * reOrdered = tempModel.reorder(mark);
2750	assert (reOrdered);
2751	tempModel=*reOrdered;
2752	delete reOrdered;
2753	delete [] mark;
2754	newSolver.loadFromCoinModel(tempModel,true);
2755	for (int iObject =0;iObject<numberObjects_;iObject++) {
2756	OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
2757	if (obj&&obj->priority()<biLinearPriority_) {
2758	int iColumn = obj->columnNumber();
2759	double value = solution[iColumn];
2760	value = ceil(value-1.0e-7);
2761	newSolver.setColLower(iColumn,value);
2762	newSolver.setColUpper(iColumn,value);
2763	}
2764	OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
2765	if (objB) {
2766	// if one or both continuous then fix one
2767	if (objB->xMeshSize()<1.0) {
2768	int xColumn = objB->xColumn();
2769	double value = solution[xColumn];
2770	newSolver.setColLower(xColumn,value);
2771	newSolver.setColUpper(xColumn,value);
2772	} else if (objB->yMeshSize()<1.0) {
2773	int yColumn = objB->yColumn();
2774	double value = solution[yColumn];
2775	newSolver.setColLower(yColumn,value);
2776	newSolver.setColUpper(yColumn,value);
2777	}
2778	}
2779	}
2780	}
2781	CbcModel model(newSolver);
2782	// Now do requested saves and modifications
2783	CbcModel * cbcModel = & model;
2784	OsiSolverInterface * osiModel = model.solver();
2785	OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
2786	ClpSimplex * clpModel = osiclpModel->getModelPtr();
2787	CglProbing probing;
2788	probing.setMaxProbe(10);
2789	probing.setMaxLook(10);
2790	probing.setMaxElements(200);
2791	probing.setMaxProbeRoot(50);
2792	probing.setMaxLookRoot(10);
2793	probing.setRowCuts(3);
2794	probing.setRowCuts(0);
2795	probing.setUsingObjective(true);
2796	cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
2797
2798	CglGomory gomory;
2799	gomory.setLimitAtRoot(512);
2800	cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
2801
2802	CglKnapsackCover knapsackCover;
2803	cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
2804
2805	CglClique clique;
2806	clique.setStarCliqueReport(false);
2807	clique.setRowCliqueReport(false);
2808	clique.setMinViolation(0.1);
2809	cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
2810	CglMixedIntegerRounding2 mixedIntegerRounding2;
2811	cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
2812
2813	CglFlowCover flowCover;
2814	cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
2815
2816	CglTwomir twomir;
2817	twomir.setMaxElements(250);
2818	cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
2819	cbcModel->cutGenerator(6)->setTiming(true);
2820
2821	CbcHeuristicFPump heuristicFPump(*cbcModel);
2822	heuristicFPump.setWhen(1);
2823	heuristicFPump.setMaximumPasses(20);
2824	heuristicFPump.setDefaultRounding(0.5);
2825	cbcModel->addHeuristic(&heuristicFPump);
2826
2827	CbcRounding rounding(*cbcModel);
2828	cbcModel->addHeuristic(&rounding);
2829
2830	CbcHeuristicLocal heuristicLocal(*cbcModel);
2831	heuristicLocal.setSearchType(1);
2832	cbcModel->addHeuristic(&heuristicLocal);
2833
2834	CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
2835	cbcModel->addHeuristic(&heuristicGreedyCover);
2836
2837	CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
2838	cbcModel->addHeuristic(&heuristicGreedyEquality);
2839
2840	CbcCompareDefault compare;
2841	cbcModel->setNodeComparison(compare);
2842	cbcModel->setNumberBeforeTrust(5);
2843	cbcModel->setSpecialOptions(2);
2844	cbcModel->messageHandler()->setLogLevel(1);
2845	cbcModel->setMaximumCutPassesAtRoot(-100);
2846	cbcModel->setMaximumCutPasses(1);
2847	cbcModel->setMinimumDrop(0.05);
2848	clpModel->setNumberIterations(1);
2849	// For branchAndBound this may help
2850	clpModel->defaultFactorizationFrequency();
2851	clpModel->setDualBound(6.71523e+07);
2852	clpModel->setPerturbation(50);
2853	osiclpModel->setSpecialOptions(193);
2854	osiclpModel->messageHandler()->setLogLevel(0);
2855	osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
2856	osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
2857	// You can save some time by switching off message building
2858	// clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
2859	// Solve
2860
2861	cbcModel->initialSolve();
2862	//double cutoff = model_->getCutoff();
2863	if (!cbcModel_)
2864	cbcModel->setCutoff(1.0e50);
2865	else
2866	cbcModel->setCutoff(cbcModel_->getCutoff());
2867	// to change exits
2868	bool isFeasible=false;
2869	int saveLogLevel=clpModel->logLevel();
2870	clpModel->setLogLevel(0);
2871	int returnCode=0;
2872	if (clpModel->tightenPrimalBounds()!=0) {
2873	clpModel->setLogLevel(saveLogLevel);
2874	returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
2875	//clpModel->writeMps("infeas2.mps");
2876	} else {
2877	clpModel->setLogLevel(saveLogLevel);
2878	clpModel->dual(); // clean up
2879	// compute some things using problem size
2880	cbcModel->setMinimumDrop(min(5.0e-2,
2881	fabs(cbcModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
2882	if (cbcModel->getNumCols()<500)
2883	cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
2884	else if (cbcModel->getNumCols()<5000)
2885	cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop
2886	else
2887	cbcModel->setMaximumCutPassesAtRoot(20);
2888	cbcModel->setMaximumCutPasses(1);
2889	// Hand coded preprocessing
2890	CglPreProcess process;
2891	OsiSolverInterface * saveSolver=cbcModel->solver()->clone();
2892	// Tell solver we are in Branch and Cut
2893	saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
2894	// Default set of cut generators
2895	CglProbing generator1;
2896	generator1.setUsingObjective(true);
2897	generator1.setMaxPass(3);
2898	generator1.setMaxProbeRoot(saveSolver->getNumCols());
2899	generator1.setMaxElements(100);
2900	generator1.setMaxLookRoot(50);
2901	generator1.setRowCuts(3);
2902	// Add in generators
2903	process.addCutGenerator(&generator1);
2904	process.messageHandler()->setLogLevel(cbcModel->logLevel());
2905	OsiSolverInterface * solver2 =
2906	process.preProcessNonDefault(*saveSolver,0,10);
2907	// Tell solver we are not in Branch and Cut
2908	saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
2909	if (solver2)
2910	solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
2911	if (!solver2) {
2912	std::cout<<"Pre-processing says infeasible!"<<std::endl;
2913	delete saveSolver;
2914	returnCode=-1;
2915	} else {
2916	std::cout<<"processed model has "<<solver2->getNumRows()
2917	<<" rows, "<<solver2->getNumCols()
2918	<<" and "<<solver2->getNumElements()<<std::endl;
2919	// we have to keep solver2 so pass clone
2920	solver2 = solver2->clone();
2921	//solver2->writeMps("intmodel");
2922	cbcModel->assignSolver(solver2);
2923	cbcModel->initialSolve();
2924	cbcModel->branchAndBound();
2925	// For best solution
2926	int numberColumns = newSolver.getNumCols();
2927	if (cbcModel->getMinimizationObjValue()<1.0e50) {
2928	// post process
2929	process.postProcess(*cbcModel->solver());
2930	// Solution now back in saveSolver
2931	cbcModel->assignSolver(saveSolver);
2932	memcpy(cbcModel->bestSolution(),cbcModel->solver()->getColSolution(),
2933	numberColumns*sizeof(double));
2934	// put back in original solver
2935	newSolver.setColSolution(cbcModel->bestSolution());
2936	isFeasible=true;
2937	} else {
2938	delete saveSolver;
2939	}
2940	}
2941	}
2942	assert (!returnCode);
2943	abort();
2944	return solution;
2945	}
2946	// Analyze constraints to see which are convex (quadratic)
2947	void
2948	OsiSolverLink::analyzeObjects()
2949	{
2950	// space for starts
2951	int numberColumns = coinModel_.numberColumns();
2952	int * start = new int [numberColumns+1];
2953	const double * rowLower = getRowLower();
2954	const double * rowUpper = getRowUpper();
2955	for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
2956	int iRow = rowNonLinear_[iNon];
2957	int numberElements = startNonLinear_[iNon+1]-startNonLinear_[iNon];
2958	// triplet arrays
2959	int * iColumn = new int [2*numberElements+1];
2960	int * jColumn = new int [2*numberElements];
2961	double * element = new double [2*numberElements];
2962	int i;
2963	int n=0;
2964	for ( i =startNonLinear_[iNon];i<startNonLinear_[iNon+1];i++) {
2965	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
2966	assert (obj);
2967	int xColumn = obj->xColumn();
2968	int yColumn = obj->yColumn();
2969	double coefficient = obj->coefficient();
2970	if (xColumn!=yColumn) {
2971	iColumn[n]=xColumn;
2972	jColumn[n]=yColumn;
2973	element[n++]=coefficient;
2974	iColumn[n]=yColumn;
2975	jColumn[n]=xColumn;
2976	element[n++]=coefficient;
2977	} else {
2978	iColumn[n]=xColumn;
2979	jColumn[n]=xColumn;
2980	element[n++]=coefficient;
2981	}
2982	}
2983	// First sort in column order
2984	CoinSort_3(iColumn,iColumn+n,jColumn,element);
2985	// marker at end
2986	iColumn[n]=numberColumns;
2987	int lastI=iColumn[0];
2988	// compute starts
2989	start[0]=0;
2990	for (i=1;i<n+1;i++) {
2991	if (iColumn[i]!=lastI) {
2992	while (lastI<iColumn[i]) {
2993	start[lastI+1]=i;
2994	lastI++;
2995	}
2996	lastI=iColumn[i];
2997	}
2998	}
2999	// -1 unknown, 0 convex, 1 nonconvex
3000	int status=-1;
3001	int statusNegative=-1;
3002	int numberLong=0; // number with >2 elements
3003	for (int k=0;k<numberColumns;k++) {
3004	int first = start[k];
3005	int last = start[k+1];
3006	if (last>first) {
3007	int j;
3008	double diagonal = 0.0;
3009	int whichK=-1;
3010	for (j=first;j<last;j++) {
3011	if (jColumn[j]==k) {
3012	diagonal = element[j];
3013	status = diagonal >0 ? 0 : 1;
3014	statusNegative = diagonal <0 ? 0 : 1;
3015	whichK = (j==first) ? j+1 :j-1;
3016	break;
3017	}
3018	}
3019	if (last==first+1) {
3020	// just one entry
3021	if (!diagonal) {
3022	// one off diagonal - not positive semi definite
3023	status=1;
3024	statusNegative=1;
3025	}
3026	} else if (diagonal) {
3027	if (last==first+2) {
3028	// other column and element
3029	double otherElement=element[whichK];;
3030	int otherColumn = jColumn[whichK];
3031	double otherDiagonal=0.0;
3032	// check 2x2 determinant - unless past and 2 long
3033	if (otherColumn>i\|\|start[otherColumn+1]>start[otherColumn]+2) {
3034	for (j=start[otherColumn];j<start[otherColumn+1];j++) {
3035	if (jColumn[j]==otherColumn) {
3036	otherDiagonal = element[j];
3037	break;
3038	}
3039	}
3040	// determinant
3041	double determinant = diagonalotherDiagonal-otherElementotherElement;
3042	if (determinant<-1.0e-12) {
3043	// not positive semi definite
3044	status=1;
3045	statusNegative=1;
3046	} else if (start[otherColumn+1]>start[otherColumn]+2&&determinant<1.0e-12) {
3047	// not positive semi definite
3048	status=1;
3049	statusNegative=1;
3050	}
3051	}
3052	} else {
3053	numberLong++;
3054	}
3055	}
3056	}
3057	}
3058	if ((status==0\|\|statusNegative==0)&&numberLong) {
3059	// need to do more work
3060	//printf("Needs more work\n");
3061	}
3062	assert (status>0\|\|statusNegative>0);
3063	if (!status) {
3064	convex_[iNon]=1;
3065	// equality may be ok
3066	if (rowUpper[iRow]<1.0e20)
3067	specialOptions2_ \|= 8;
3068	else
3069	convex_[iNon]=0;
3070	} else if (!statusNegative) {
3071	convex_[iNon]=-1;
3072	// equality may be ok
3073	if (rowLower[iRow]>-1.0e20)
3074	specialOptions2_ \|= 8;
3075	else
3076	convex_[iNon]=0;
3077	} else {
3078	convex_[iNon]=0;
3079	}
3080	//printf("Convexity of row %d is %d\n",iRow,convex_[iNon]);
3081	delete [] iColumn;
3082	delete [] jColumn;
3083	delete [] element;
3084	}
3085	delete [] start;
3086	}
3087	//-------------------------------------------------------------------
3088	// Clone
3089	//-------------------------------------------------------------------
3090	OsiSolverInterface *
3091	OsiSolverLink::clone(bool copyData) const
3092	{
3093	assert (copyData);
3094	OsiSolverLink * newModel = new OsiSolverLink(*this);
3095	return newModel;
3096	}
3097
3098
3099	//-------------------------------------------------------------------
3100	// Copy constructor
3101	//-------------------------------------------------------------------
3102	OsiSolverLink::OsiSolverLink (
3103	const OsiSolverLink & rhs)
3104	: CbcOsiSolver(rhs)
3105	{
3106	gutsOfDestructor(true);
3107	gutsOfCopy(rhs);
3108	// something odd happens - try this
3109	OsiSolverInterface::operator=(rhs);
3110	}
3111
3112	//-------------------------------------------------------------------
3113	// Destructor
3114	//-------------------------------------------------------------------
3115	OsiSolverLink::~OsiSolverLink ()
3116	{
3117	gutsOfDestructor();
3118	}
3119
3120	//-------------------------------------------------------------------
3121	// Assignment operator
3122	//-------------------------------------------------------------------
3123	OsiSolverLink &
3124	OsiSolverLink::operator=(const OsiSolverLink& rhs)
3125	{
3126	if (this != &rhs) {
3127	gutsOfDestructor();
3128	CbcOsiSolver::operator=(rhs);
3129	gutsOfCopy(rhs);
3130	}
3131	return *this;
3132	}
3133	void
3134	OsiSolverLink::gutsOfDestructor(bool justNullify)
3135	{
3136	if (!justNullify) {
3137	delete matrix_;
3138	delete originalRowCopy_;
3139	delete [] info_;
3140	delete [] bestSolution_;
3141	delete quadraticModel_;
3142	delete [] startNonLinear_;
3143	delete [] rowNonLinear_;
3144	delete [] convex_;
3145	delete [] whichNonLinear_;
3146	delete [] fixVariables_;
3147	}
3148	matrix_ = NULL;
3149	originalRowCopy_ = NULL;
3150	quadraticModel_ = NULL;
3151	numberNonLinearRows_=0;
3152	startNonLinear_ = NULL;
3153	rowNonLinear_ = NULL;
3154	convex_ = NULL;
3155	whichNonLinear_ = NULL;
3156	info_ = NULL;
3157	fixVariables_=NULL;
3158	numberVariables_ = 0;
3159	specialOptions2_ = 0;
3160	objectiveRow_=-1;
3161	objectiveVariable_=-1;
3162	bestSolution_ = NULL;
3163	bestObjectiveValue_ =1.0e100;
3164	defaultMeshSize_ = 1.0e-4;
3165	defaultBound_ = 1.0e5;
3166	integerPriority_ = 1000;
3167	biLinearPriority_ = 10000;
3168	numberFix_=0;
3169	}
3170	void
3171	OsiSolverLink::gutsOfCopy(const OsiSolverLink & rhs)
3172	{
3173	coinModel_ = rhs.coinModel_;
3174	numberVariables_ = rhs.numberVariables_;
3175	numberNonLinearRows_ = rhs.numberNonLinearRows_;
3176	specialOptions2_ = rhs.specialOptions2_;
3177	objectiveRow_=rhs.objectiveRow_;
3178	objectiveVariable_=rhs.objectiveVariable_;
3179	bestObjectiveValue_ = rhs.bestObjectiveValue_;
3180	defaultMeshSize_ = rhs.defaultMeshSize_;
3181	defaultBound_ = rhs.defaultBound_;
3182	integerPriority_ = rhs.integerPriority_;
3183	biLinearPriority_ = rhs.biLinearPriority_;
3184	numberFix_ = rhs.numberFix_;
3185	if (numberVariables_) {
3186	if (rhs.matrix_)
3187	matrix_ = new CoinPackedMatrix(*rhs.matrix_);
3188	else
3189	matrix_=NULL;
3190	if (rhs.originalRowCopy_)
3191	originalRowCopy_ = new CoinPackedMatrix(*rhs.originalRowCopy_);
3192	else
3193	originalRowCopy_=NULL;
3194	info_ = new OsiLinkedBound [numberVariables_];
3195	for (int i=0;i<numberVariables_;i++) {
3196	info_[i] = OsiLinkedBound(rhs.info_[i]);
3197	}
3198	if (rhs.bestSolution_) {
3199	bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->getNumCols());
3200	} else {
3201	bestSolution_=NULL;
3202	}
3203	}
3204	if (numberNonLinearRows_) {
3205	startNonLinear_ = CoinCopyOfArray(rhs.startNonLinear_,numberNonLinearRows_+1);
3206	rowNonLinear_ = CoinCopyOfArray(rhs.rowNonLinear_,numberNonLinearRows_);
3207	convex_ = CoinCopyOfArray(rhs.convex_,numberNonLinearRows_);
3208	int numberEntries = startNonLinear_[numberNonLinearRows_];
3209	whichNonLinear_ = CoinCopyOfArray(rhs.whichNonLinear_,numberEntries);
3210	}
3211	if (rhs.quadraticModel_) {
3212	quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
3213	} else {
3214	quadraticModel_ = NULL;
3215	}
3216	fixVariables_ = CoinCopyOfArray(rhs.fixVariables_,numberFix_);
3217	}
3218	// Add a bound modifier
3219	void
3220	OsiSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected,
3221	double multiplier)
3222	{
3223	bool found=false;
3224	int i;
3225	for ( i=0;i<numberVariables_;i++) {
3226	if (info_[i].variable()==whichVariable) {
3227	found=true;
3228	break;
3229	}
3230	}
3231	if (!found) {
3232	// add in
3233	OsiLinkedBound * temp = new OsiLinkedBound [numberVariables_+1];
3234	for (int i=0;i<numberVariables_;i++)
3235	temp[i]= info_[i];
3236	delete [] info_;
3237	info_=temp;
3238	info_[numberVariables_++] = OsiLinkedBound(this,whichVariable,0,NULL,NULL,NULL);
3239	}
3240	info_[i].addBoundModifier(upperBoundAffected, useUpperBound,whichVariableAffected,multiplier);
3241	}
3242	// Update coefficients
3243	int
3244	OsiSolverLink::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
3245	{
3246	double * lower = solver->columnLower();
3247	double * upper = solver->columnUpper();
3248	double * objective = solver->objective();
3249	int numberChanged=0;
3250	for (int iObject =0;iObject<numberObjects_;iObject++) {
3251	OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
3252	if (obj) {
3253	numberChanged +=obj->updateCoefficients(lower,upper,objective,matrix,&basis_);
3254	}
3255	}
3256	return numberChanged;
3257	}
3258	// Set best solution found internally
3259	void
3260	OsiSolverLink::setBestSolution(const double * solution, int numberColumns)
3261	{
3262	delete [] bestSolution_;
3263	int numberColumnsThis = modelPtr_->numberColumns();
3264	bestSolution_ = new double [numberColumnsThis];
3265	CoinZeroN(bestSolution_,numberColumnsThis);
3266	memcpy(bestSolution_,solution,CoinMin(numberColumns,numberColumnsThis)*sizeof(double));
3267	}
3268	/* Two tier integer problem where when set of variables with priority
3269	less than this are fixed the problem becomes an easier integer problem
3270	*/
3271	void
3272	OsiSolverLink::setFixedPriority(int priorityValue)
3273	{
3274	delete [] fixVariables_;
3275	fixVariables_=NULL;
3276	numberFix_=0;
3277	int i;
3278	for ( i =0;i<numberObjects_;i++) {
3279	OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[i]);
3280	if (obj) {
3281	int iColumn = obj->columnNumber();
3282	assert (iColumn>=0);
3283	if (obj->priority()<priorityValue)
3284	numberFix_++;
3285	}
3286	}
3287	if (numberFix_) {
3288	specialOptions2_ \|= 1;
3289	fixVariables_ = new int [numberFix_];
3290	numberFix_=0;
3291	// need to make sure coinModel_ is correct
3292	int numberColumns = coinModel_.numberColumns();
3293	char * highPriority = new char [numberColumns];
3294	CoinZeroN(highPriority,numberColumns);
3295	for ( i =0;i<numberObjects_;i++) {
3296	OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[i]);
3297	if (obj) {
3298	int iColumn = obj->columnNumber();
3299	assert (iColumn>=0);
3300	if (iColumn<numberColumns) {
3301	if (obj->priority()<priorityValue) {
3302	object_[i]=new OsiSimpleFixedInteger(*obj);
3303	delete obj;
3304	fixVariables_[numberFix_++]=iColumn;
3305	highPriority[iColumn]=1;
3306	}
3307	}
3308	}
3309	}
3310	CoinModel * newModel = coinModel_.reorder(highPriority);
3311	if (newModel) {
3312	coinModel_ = * newModel;
3313	} else {
3314	printf("Unable to use priorities\n");
3315	delete [] fixVariables_;
3316	fixVariables_ = NULL;
3317	numberFix_=0;
3318	}
3319	delete newModel;
3320	delete [] highPriority;
3321	}
3322	}
3323	// Gets correct form for a quadratic row - user to delete
3324	CoinPackedMatrix *
3325	OsiSolverLink::quadraticRow(int rowNumber,double * linearRow) const
3326	{
3327	int numberColumns = coinModel_.numberColumns();
3328	int numberRows = coinModel_.numberRows();
3329	CoinZeroN(linearRow,numberColumns);
3330	int numberElements=0;
3331	assert (rowNumber>=0&&rowNumber<numberRows);
3332	CoinModelLink triple=coinModel_.firstInRow(rowNumber);
3333	while (triple.column()>=0) {
3334	int iColumn = triple.column();
3335	const char * expr = coinModel_.getElementAsString(rowNumber,iColumn);
3336	if (strcmp(expr,"Numeric")) {
3337	// try and see which columns
3338	assert (strlen(expr)<20000);
3339	char temp[20000];
3340	strcpy(temp,expr);
3341	char * pos = temp;
3342	bool ifFirst=true;
3343	while (*pos) {
3344	double value;
3345	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
3346	// must be column unless first when may be linear term
3347	if (jColumn>=0) {
3348	numberElements++;
3349	} else if (jColumn==-2) {
3350	linearRow[iColumn]=value;
3351	} else {
3352	printf("bad nonlinear term %s\n",temp);
3353	abort();
3354	}
3355	ifFirst=false;
3356	}
3357	} else {
3358	linearRow[iColumn]=coinModel_.getElement(rowNumber,iColumn);
3359	}
3360	triple=coinModel_.next(triple);
3361	}
3362	if (!numberElements) {
3363	return NULL;
3364	} else {
3365	int * column = new int[numberElements];
3366	int * column2 = new int[numberElements];
3367	double * element = new double[numberElements];
3368	numberElements=0;
3369	CoinModelLink triple=coinModel_.firstInRow(rowNumber);
3370	while (triple.column()>=0) {
3371	int iColumn = triple.column();
3372	const char * expr = coinModel_.getElementAsString(rowNumber,iColumn);
3373	if (strcmp(expr,"Numeric")) {
3374	// try and see which columns
3375	assert (strlen(expr)<20000);
3376	char temp[20000];
3377	strcpy(temp,expr);
3378	char * pos = temp;
3379	bool ifFirst=true;
3380	while (*pos) {
3381	double value;
3382	int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
3383	// must be column unless first when may be linear term
3384	if (jColumn>=0) {
3385	column[numberElements]=iColumn;
3386	column2[numberElements]=jColumn;
3387	element[numberElements++]=value;
3388	} else if (jColumn!=-2) {
3389	printf("bad nonlinear term %s\n",temp);
3390	abort();
3391	}
3392	ifFirst=false;
3393	}
3394	}
3395	triple=coinModel_.next(triple);
3396	}
3397	return new CoinPackedMatrix(true,column2,column,element,numberElements);
3398	}
3399	}
3400	/*
3401	Problem specific
3402	Returns -1 if node fathomed and no solution
3403	0 if did nothing
3404	1 if node fathomed and solution
3405	allFixed is true if all LinkedBound variables are fixed
3406	*/
3407	int
3408	OsiSolverLink::fathom(bool allFixed)
3409	{
3410	int returnCode=0;
3411	if (allFixed) {
3412	// solve anyway
3413	OsiClpSolverInterface::resolve();
3414	if (!isProvenOptimal()) {
3415	printf("cutoff before fathoming\n");
3416	return -1;
3417	}
3418	// all fixed so we can reformulate
3419	OsiClpSolverInterface newSolver;
3420	// set values
3421	const double * lower = modelPtr_->columnLower();
3422	const double * upper = modelPtr_->columnUpper();
3423	int i;
3424	for (i=0;i<numberFix_;i++ ) {
3425	int iColumn = fixVariables_[i];
3426	double lo = lower[iColumn];
3427	double up = upper[iColumn];
3428	assert (lo==up);
3429	//printf("column %d fixed to %g\n",iColumn,lo);
3430	coinModel_.associateElement(coinModel_.columnName(iColumn),lo);
3431	}
3432	newSolver.loadFromCoinModel(coinModel_,true);
3433	for (i=0;i<numberFix_;i++ ) {
3434	int iColumn = fixVariables_[i];
3435	newSolver.setColLower(iColumn,lower[iColumn]);
3436	newSolver.setColUpper(iColumn,lower[iColumn]);
3437	}
3438	// see if everything with objective fixed
3439	const double * objective = modelPtr_->objective();
3440	int numberColumns = newSolver.getNumCols();
3441	bool zeroObjective=true;
3442	double sum=0.0;
3443	for (i=0;i<numberColumns;i++) {
3444	if (upper[i]>lower[i]&&objective[i]) {
3445	zeroObjective=false;
3446	break;
3447	} else {
3448	sum += lower[i]*objective[i];
3449	}
3450	}
3451	int fake[]={5,4,3,2,0,0,0};
3452	bool onOptimalPath=true;
3453	for (i=0;i<7;i++) {
3454	if ((int) upper[i]!=fake[i])
3455	onOptimalPath=false;
3456	}
3457	if (onOptimalPath)
3458	printf("possible\n");
3459	if (zeroObjective) {
3460	// randomize objective
3461	ClpSimplex * clpModel = newSolver.getModelPtr();
3462	const double * element = clpModel->matrix()->getMutableElements();
3463	//const int * row = clpModel->matrix()->getIndices();
3464	const CoinBigIndex * columnStart = clpModel->matrix()->getVectorStarts();
3465	const int * columnLength = clpModel->matrix()->getVectorLengths();
3466	double * objective = clpModel->objective();
3467	for (i=0;i<numberColumns;i++) {
3468	if (clpModel->isInteger(i)) {
3469	double value=0.0;
3470	for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
3471	value += fabs(element[j]);
3472	}
3473	objective[i]=value;
3474	}
3475	}
3476	}
3477	//newSolver.writeMps("xx");
3478	CbcModel model(newSolver);
3479	// Now do requested saves and modifications
3480	CbcModel * cbcModel = & model;
3481	OsiSolverInterface * osiModel = model.solver();
3482	OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
3483	ClpSimplex * clpModel = osiclpModel->getModelPtr();
3484	CglProbing probing;
3485	probing.setMaxProbe(10);
3486	probing.setMaxLook(10);
3487	probing.setMaxElements(200);
3488	probing.setMaxProbeRoot(50);
3489	probing.setMaxLookRoot(10);
3490	probing.setRowCuts(3);
3491	probing.setRowCuts(0);
3492	probing.setUsingObjective(true);
3493	cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
3494
3495	CglGomory gomory;
3496	gomory.setLimitAtRoot(512);
3497	cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
3498
3499	CglKnapsackCover knapsackCover;
3500	cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
3501
3502	CglClique clique;
3503	clique.setStarCliqueReport(false);
3504	clique.setRowCliqueReport(false);
3505	clique.setMinViolation(0.1);
3506	cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
3507	CglMixedIntegerRounding2 mixedIntegerRounding2;
3508	cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
3509
3510	CglFlowCover flowCover;
3511	cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
3512
3513	CglTwomir twomir;
3514	twomir.setMaxElements(250);
3515	cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
3516	cbcModel->cutGenerator(6)->setTiming(true);
3517
3518	CbcHeuristicFPump heuristicFPump(*cbcModel);
3519	heuristicFPump.setWhen(1);
3520	heuristicFPump.setMaximumPasses(20);
3521	heuristicFPump.setDefaultRounding(0.5);
3522	cbcModel->addHeuristic(&heuristicFPump);
3523
3524	CbcRounding rounding(*cbcModel);
3525	cbcModel->addHeuristic(&rounding);
3526
3527	CbcHeuristicLocal heuristicLocal(*cbcModel);
3528	heuristicLocal.setSearchType(1);
3529	cbcModel->addHeuristic(&heuristicLocal);
3530
3531	CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
3532	cbcModel->addHeuristic(&heuristicGreedyCover);
3533
3534	CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
3535	cbcModel->addHeuristic(&heuristicGreedyEquality);
3536
3537	CbcCompareDefault compare;
3538	cbcModel->setNodeComparison(compare);
3539	cbcModel->setNumberBeforeTrust(5);
3540	cbcModel->setSpecialOptions(2);
3541	cbcModel->messageHandler()->setLogLevel(1);
3542	cbcModel->setMaximumCutPassesAtRoot(-100);
3543	cbcModel->setMaximumCutPasses(1);
3544	cbcModel->setMinimumDrop(0.05);
3545	clpModel->setNumberIterations(1);
3546	// For branchAndBound this may help
3547	clpModel->defaultFactorizationFrequency();
3548	clpModel->setDualBound(6.71523e+07);
3549	clpModel->setPerturbation(50);
3550	osiclpModel->setSpecialOptions(193);
3551	osiclpModel->messageHandler()->setLogLevel(0);
3552	osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
3553	osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
3554	// You can save some time by switching off message building
3555	// clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
3556	// Solve
3557
3558	cbcModel->initialSolve();
3559	//double cutoff = model_->getCutoff();
3560	if (zeroObjective\|\|!cbcModel_)
3561	cbcModel->setCutoff(1.0e50);
3562	else
3563	cbcModel->setCutoff(cbcModel_->getCutoff());
3564	// to change exits
3565	bool isFeasible=false;
3566	int saveLogLevel=clpModel->logLevel();
3567	clpModel->setLogLevel(0);
3568	if (clpModel->tightenPrimalBounds()!=0) {
3569	clpModel->setLogLevel(saveLogLevel);
3570	returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
3571	} else {
3572	clpModel->setLogLevel(saveLogLevel);
3573	clpModel->dual(); // clean up
3574	// compute some things using problem size
3575	cbcModel->setMinimumDrop(min(5.0e-2,
3576	fabs(cbcModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
3577	if (cbcModel->getNumCols()<500)
3578	cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
3579	else if (cbcModel->getNumCols()<5000)
3580	cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop
3581	else
3582	cbcModel->setMaximumCutPassesAtRoot(20);
3583	cbcModel->setMaximumCutPasses(1);
3584	// Hand coded preprocessing
3585	CglPreProcess process;
3586	OsiSolverInterface * saveSolver=cbcModel->solver()->clone();
3587	// Tell solver we are in Branch and Cut
3588	saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
3589	// Default set of cut generators
3590	CglProbing generator1;
3591	generator1.setUsingObjective(true);
3592	generator1.setMaxPass(3);
3593	generator1.setMaxProbeRoot(saveSolver->getNumCols());
3594	generator1.setMaxElements(100);
3595	generator1.setMaxLookRoot(50);
3596	generator1.setRowCuts(3);
3597	// Add in generators
3598	process.addCutGenerator(&generator1);
3599	process.messageHandler()->setLogLevel(cbcModel->logLevel());
3600	OsiSolverInterface * solver2 =
3601	process.preProcessNonDefault(*saveSolver,0,10);
3602	// Tell solver we are not in Branch and Cut
3603	saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
3604	if (solver2)
3605	solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
3606	if (!solver2) {
3607	std::cout<<"Pre-processing says infeasible!"<<std::endl;
3608	delete saveSolver;
3609	returnCode=-1;
3610	} else {
3611	std::cout<<"processed model has "<<solver2->getNumRows()
3612	<<" rows, "<<solver2->getNumCols()
3613	<<" and "<<solver2->getNumElements()<<std::endl;
3614	// we have to keep solver2 so pass clone
3615	solver2 = solver2->clone();
3616	//solver2->writeMps("intmodel");
3617	cbcModel->assignSolver(solver2);
3618	cbcModel->initialSolve();
3619	if (zeroObjective) {
3620	cbcModel->setMaximumSolutions(1); // just getting a solution
3621	#if 0
3622	OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (cbcModel->solver());
3623	ClpSimplex * clpModel = osiclpModel->getModelPtr();
3624	const double * element = clpModel->matrix()->getMutableElements();
3625	//const int * row = clpModel->matrix()->getIndices();
3626	const CoinBigIndex * columnStart = clpModel->matrix()->getVectorStarts();
3627	const int * columnLength = clpModel->matrix()->getVectorLengths();
3628	int n=clpModel->numberColumns();
3629	int * sort2 = new int[n];
3630	int * pri = new int[n];
3631	double * sort = new double[n];
3632	int i;
3633	int nint=0;
3634	for (i=0;i<n;i++) {
3635	if (clpModel->isInteger(i)) {
3636	double largest=0.0;
3637	for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
3638	largest = CoinMax(largest,fabs(element[j]));
3639	}
3640	sort2[nint]=nint;
3641	sort[nint++]=-largest;
3642	}
3643	}
3644	CoinSort_2(sort,sort+nint,sort2);
3645	int kpri=1;
3646	double last = sort[0];
3647	for (i=0;i<nint;i++) {
3648	if (sort[i]!=last) {
3649	kpri++;
3650	last=sort[i];
3651	}
3652	pri[sort2[i]]=kpri;
3653	}
3654	cbcModel->passInPriorities(pri,false);
3655	delete [] sort;
3656	delete [] sort2;
3657	delete [] pri;
3658	#endif
3659	}
3660	cbcModel->branchAndBound();
3661	// For best solution
3662	int numberColumns = newSolver.getNumCols();
3663	if (cbcModel->getMinimizationObjValue()<1.0e50) {
3664	// post process
3665	process.postProcess(*cbcModel->solver());
3666	// Solution now back in saveSolver
3667	cbcModel->assignSolver(saveSolver);
3668	memcpy(cbcModel->bestSolution(),cbcModel->solver()->getColSolution(),
3669	numberColumns*sizeof(double));
3670	// put back in original solver
3671	newSolver.setColSolution(cbcModel->bestSolution());
3672	isFeasible=true;
3673	} else {
3674	delete saveSolver;
3675	}
3676	}
3677	//const double * solution = newSolver.getColSolution();
3678	if (isFeasible&&cbcModel->getMinimizationObjValue()<1.0e50) {
3679	int numberColumns = this->getNumCols();
3680	int i;
3681	const double * solution = cbcModel->bestSolution();
3682	int numberColumns2 = newSolver.getNumCols();
3683	for (i=0;i<numberColumns2;i++) {
3684	double value = solution[i];
3685	assert (fabs(value-floor(value+0.5))<0.0001);
3686	value = floor(value+0.5);
3687	this->setColLower(i,value);
3688	this->setColUpper(i,value);
3689	}
3690	for (;i<numberColumns;i++) {
3691	this->setColLower(i,0.0);
3692	this->setColUpper(i,1.1);
3693	}
3694	// but take off cuts
3695	int numberRows = getNumRows();
3696	int numberRows2 = cbcModel_->continuousSolver()->getNumRows();
3697
3698	for (i=numberRows2;i<numberRows;i++)
3699	setRowBounds(i,-COIN_DBL_MAX,COIN_DBL_MAX);
3700	initialSolve();
3701	//if (!isProvenOptimal())
3702	//getModelPtr()->writeMps("bad.mps");
3703	if (isProvenOptimal()) {
3704	delete [] bestSolution_;
3705	bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols());
3706	bestObjectiveValue_ = modelPtr_->objectiveValue();
3707	printf("BB best value %g\n",bestObjectiveValue_);
3708	returnCode = 1;
3709	} else {
3710	printf("*** WHY BAD SOL\n");
3711	returnCode=-1;
3712	}
3713	} else {
3714	modelPtr_->setProblemStatus(1);
3715	modelPtr_->setObjectiveValue(COIN_DBL_MAX);
3716	returnCode = -1;
3717	}
3718	}
3719	}
3720	return returnCode;
3721	}
3722	//#############################################################################
3723	// Constructors, destructors and assignment
3724	//#############################################################################
3725
3726	//-------------------------------------------------------------------
3727	// Default Constructor
3728	//-------------------------------------------------------------------
3729	OsiLinkedBound::OsiLinkedBound ()
3730	{
3731	model_ = NULL;
3732	variable_ = -1;
3733	numberAffected_ = 0;
3734	maximumAffected_ = numberAffected_;
3735	affected_ = NULL;
3736	}
3737	// Useful Constructor
3738	OsiLinkedBound::OsiLinkedBound(OsiSolverInterface * model, int variable,
3739	int numberAffected, const int * positionL,
3740	const int * positionU, const double * multiplier)
3741	{
3742	model_ = model;
3743	variable_ = variable;
3744	numberAffected_ = 2*numberAffected;
3745	maximumAffected_ = numberAffected_;
3746	if (numberAffected_) {
3747	affected_ = new boundElementAction[numberAffected_];
3748	int n=0;
3749	for (int i=0;i<numberAffected;i++) {
3750	// LB
3751	boundElementAction action;
3752	action.affect=2;
3753	action.ubUsed=0;
3754	action.type=0;
3755	action.affected=positionL[i];
3756	action.multiplier=multiplier[i];
3757	affected_[n++]=action;
3758	// UB
3759	action.affect=2;
3760	action.ubUsed=1;
3761	action.type=0;
3762	action.affected=positionU[i];
3763	action.multiplier=multiplier[i];
3764	affected_[n++]=action;
3765	}
3766	} else {
3767	affected_ = NULL;
3768	}
3769	}
3770
3771	//-------------------------------------------------------------------
3772	// Copy constructor
3773	//-------------------------------------------------------------------
3774	OsiLinkedBound::OsiLinkedBound (
3775	const OsiLinkedBound & rhs)
3776	{
3777	model_ = rhs.model_;
3778	variable_ = rhs.variable_;
3779	numberAffected_ = rhs.numberAffected_;
3780	maximumAffected_ = rhs.maximumAffected_;
3781	if (numberAffected_) {
3782	affected_ = new boundElementAction[maximumAffected_];
3783	memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
3784	} else {
3785	affected_ = NULL;
3786	}
3787	}
3788
3789	//-------------------------------------------------------------------
3790	// Destructor
3791	//-------------------------------------------------------------------
3792	OsiLinkedBound::~OsiLinkedBound ()
3793	{
3794	delete [] affected_;
3795	}
3796
3797	//-------------------------------------------------------------------
3798	// Assignment operator
3799	//-------------------------------------------------------------------
3800	OsiLinkedBound &
3801	OsiLinkedBound::operator=(const OsiLinkedBound& rhs)
3802	{
3803	if (this != &rhs) {
3804	delete [] affected_;
3805	model_ = rhs.model_;
3806	variable_ = rhs.variable_;
3807	numberAffected_ = rhs.numberAffected_;
3808	maximumAffected_ = rhs.maximumAffected_;
3809	if (numberAffected_) {
3810	affected_ = new boundElementAction[maximumAffected_];
3811	memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
3812	} else {
3813	affected_ = NULL;
3814	}
3815	}
3816	return *this;
3817	}
3818	// Add a bound modifier
3819	void
3820	OsiLinkedBound::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable,
3821	double multiplier)
3822	{
3823	if (numberAffected_==maximumAffected_) {
3824	maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
3825	boundElementAction * temp = new boundElementAction[maximumAffected_];
3826	memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
3827	delete [] affected_;
3828	affected_ = temp;
3829	}
3830	boundElementAction action;
3831	action.affect=upperBoundAffected ? 1 : 0;
3832	action.ubUsed=useUpperBound ? 1 : 0;
3833	action.type=2;
3834	action.affected=whichVariable;
3835	action.multiplier=multiplier;
3836	affected_[numberAffected_++]=action;
3837
3838	}
3839	// Update other bounds
3840	void
3841	OsiLinkedBound::updateBounds(ClpSimplex * solver)
3842	{
3843	double * lower = solver->columnLower();
3844	double * upper = solver->columnUpper();
3845	double lo = lower[variable_];
3846	double up = upper[variable_];
3847	// printf("bounds for %d are %g and %g\n",variable_,lo,up);
3848	for (int j=0;j<numberAffected_;j++) {
3849	if (affected_[j].affect<2) {
3850	double multiplier = affected_[j].multiplier;
3851	assert (affected_[j].type==2);
3852	int iColumn = affected_[j].affected;
3853	double useValue = (affected_[j].ubUsed) ? up : lo;
3854	if (affected_[j].affect==0)
3855	lower[iColumn] = CoinMin(upper[iColumn],CoinMax(lower[iColumn],multiplier*useValue));
3856	else
3857	upper[iColumn] = CoinMax(lower[iColumn],CoinMin(upper[iColumn],multiplier*useValue));
3858	}
3859	}
3860	}
3861	#if 0
3862	// Add an element modifier
3863	void
3864	OsiLinkedBound::addCoefficientModifier(bool useUpperBound, int position,
3865	double multiplier)
3866	{
3867	if (numberAffected_==maximumAffected_) {
3868	maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
3869	boundElementAction * temp = new boundElementAction[maximumAffected_];
3870	memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
3871	delete [] affected_;
3872	affected_ = temp;
3873	}
3874	boundElementAction action;
3875	action.affect=2;
3876	action.ubUsed=useUpperBound ? 1 : 0;
3877	action.type=0;
3878	action.affected=position;
3879	action.multiplier=multiplier;
3880	affected_[numberAffected_++]=action;
3881
3882	}
3883	// Update coefficients
3884	void
3885	OsiLinkedBound::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
3886	{
3887	double * lower = solver->columnLower();
3888	double * upper = solver->columnUpper();
3889	double * element = matrix->getMutableElements();
3890	double lo = lower[variable_];
3891	double up = upper[variable_];
3892	// printf("bounds for %d are %g and %g\n",variable_,lo,up);
3893	for (int j=0;j<numberAffected_;j++) {
3894	if (affected_[j].affect==2) {
3895	double multiplier = affected_[j].multiplier;
3896	assert (affected_[j].type==0);
3897	int position = affected_[j].affected;
3898	//double old = element[position];
3899	if (affected_[j].ubUsed)
3900	element[position] = multiplier*up;
3901	else
3902	element[position] = multiplier*lo;
3903	//if ( old != element[position])
3904	//printf("change at %d from %g to %g\n",position,old,element[position]);
3905	}
3906	}
3907	}
3908	#endif
3909	// Default Constructor
3910	CbcHeuristicDynamic3::CbcHeuristicDynamic3()
3911	:CbcHeuristic()
3912	{
3913	}
3914
3915	// Constructor from model
3916	CbcHeuristicDynamic3::CbcHeuristicDynamic3(CbcModel & model)
3917	:CbcHeuristic(model)
3918	{
3919	}
3920
3921	// Destructor
3922	CbcHeuristicDynamic3::~CbcHeuristicDynamic3 ()
3923	{
3924	}
3925
3926	// Clone
3927	CbcHeuristic *
3928	CbcHeuristicDynamic3::clone() const
3929	{
3930	return new CbcHeuristicDynamic3(*this);
3931	}
3932
3933	// Copy constructor
3934	CbcHeuristicDynamic3::CbcHeuristicDynamic3(const CbcHeuristicDynamic3 & rhs)
3935	:
3936	CbcHeuristic(rhs)
3937	{
3938	}
3939
3940	// Returns 1 if solution, 0 if not
3941	int
3942	CbcHeuristicDynamic3::solution(double & solutionValue,
3943	double * betterSolution)
3944	{
3945	if (!model_)
3946	return 0;
3947	OsiSolverLink * clpSolver
3948	= dynamic_cast<OsiSolverLink *> (model_->solver());
3949	assert (clpSolver);
3950	double newSolutionValue = clpSolver->bestObjectiveValue();
3951	const double * solution = clpSolver->bestSolution();
3952	if (newSolutionValue<solutionValue&&solution) {
3953	int numberColumns = clpSolver->getNumCols();
3954	// new solution
3955	memcpy(betterSolution,solution,numberColumns*sizeof(double));
3956	solutionValue = newSolutionValue;
3957	return 1;
3958	} else {
3959	return 0;
3960	}
3961	}
3962	// update model
3963	void CbcHeuristicDynamic3::setModel(CbcModel * model)
3964	{
3965	model_ = model;
3966	}
3967	// Resets stuff if model changes
3968	void
3969	CbcHeuristicDynamic3::resetModel(CbcModel * model)
3970	{
3971	model_ = model;
3972	}
3973	#include <cassert>
3974	#include <cmath>
3975	#include <cfloat>
3976	//#define CBC_DEBUG
3977
3978	#include "OsiSolverInterface.hpp"
3979	//#include "OsiBranchLink.hpp"
3980	#include "CoinError.hpp"
3981	#include "CoinHelperFunctions.hpp"
3982	#include "CoinPackedMatrix.hpp"
3983	#include "CoinWarmStartBasis.hpp"
3984
3985	// Default Constructor
3986	OsiOldLink::OsiOldLink ()
3987	: OsiSOS(),
3988	numberLinks_(0)
3989	{
3990	}
3991
3992	// Useful constructor (which are indices)
3993	OsiOldLink::OsiOldLink (const OsiSolverInterface * solver, int numberMembers,
3994	int numberLinks, int first , const double * weights, int identifier)
3995	: OsiSOS(),
3996	numberLinks_(numberLinks)
3997	{
3998	numberMembers_ = numberMembers;
3999	members_ = NULL;
4000	sosType_ = 1;
4001	if (numberMembers_) {
4002	weights_ = new double[numberMembers_];
4003	members_ = new int[numberMembers_*numberLinks_];
4004	if (weights) {
4005	memcpy(weights_,weights,numberMembers_*sizeof(double));
4006	} else {
4007	for (int i=0;i<numberMembers_;i++)
4008	weights_[i]=i;
4009	}
4010	// weights must be increasing
4011	int i;
4012	double last=-COIN_DBL_MAX;
4013	for (i=0;i<numberMembers_;i++) {
4014	assert (weights_[i]>last+1.0e-12);
4015	last=weights_[i];
4016	}
4017	for (i=0;i<numberMembers_*numberLinks_;i++) {
4018	members_[i]=first+i;
4019	}
4020	} else {
4021	weights_ = NULL;
4022	}
4023	}
4024
4025	// Useful constructor (which are indices)
4026	OsiOldLink::OsiOldLink (const OsiSolverInterface * solver, int numberMembers,
4027	int numberLinks, int sosType, const int * which , const double * weights, int identifier)
4028	: OsiSOS(),
4029	numberLinks_(numberLinks)
4030	{
4031	numberMembers_ = numberMembers;
4032	members_ = NULL;
4033	sosType_ = 1;
4034	if (numberMembers_) {
4035	weights_ = new double[numberMembers_];
4036	members_ = new int[numberMembers_*numberLinks_];
4037	if (weights) {
4038	memcpy(weights_,weights,numberMembers_*sizeof(double));
4039	} else {
4040	for (int i=0;i<numberMembers_;i++)
4041	weights_[i]=i;
4042	}
4043	// weights must be increasing
4044	int i;
4045	double last=-COIN_DBL_MAX;
4046	for (i=0;i<numberMembers_;i++) {
4047	assert (weights_[i]>last+1.0e-12);
4048	last=weights_[i];
4049	}
4050	for (i=0;i<numberMembers_*numberLinks_;i++) {
4051	members_[i]= which[i];
4052	}
4053	} else {
4054	weights_ = NULL;
4055	}
4056	}
4057
4058	// Copy constructor
4059	OsiOldLink::OsiOldLink ( const OsiOldLink & rhs)
4060	:OsiSOS(rhs)
4061	{
4062	numberLinks_ = rhs.numberLinks_;
4063	if (numberMembers_) {
4064	delete [] members_;
4065	members_ = CoinCopyOfArray(rhs.members_,numberMembers_*numberLinks_);
4066	}
4067	}
4068
4069	// Clone
4070	OsiObject *
4071	OsiOldLink::clone() const
4072	{
4073	return new OsiOldLink(*this);
4074	}
4075
4076	// Assignment operator
4077	OsiOldLink &
4078	OsiOldLink::operator=( const OsiOldLink& rhs)
4079	{
4080	if (this!=&rhs) {
4081	OsiSOS::operator=(rhs);
4082	delete [] members_;
4083	numberLinks_ = rhs.numberLinks_;
4084	if (numberMembers_) {
4085	members_ = CoinCopyOfArray(rhs.members_,numberMembers_*numberLinks_);
4086	} else {
4087	members_ = NULL;
4088	}
4089	}
4090	return *this;
4091	}
4092
4093	// Destructor
4094	OsiOldLink::~OsiOldLink ()
4095	{
4096	}
4097
4098	// Infeasibility - large is 0.5
4099	double
4100	OsiOldLink::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
4101	{
4102	int j;
4103	int firstNonZero=-1;
4104	int lastNonZero = -1;
4105	const double * solution = info->solution_;
4106	//const double * lower = info->lower_;
4107	const double * upper = info->upper_;
4108	double integerTolerance = info->integerTolerance_;
4109	double weight = 0.0;
4110	double sum =0.0;
4111
4112	// check bounds etc
4113	double lastWeight=-1.0e100;
4114	int base=0;
4115	for (j=0;j<numberMembers_;j++) {
4116	for (int k=0;k<numberLinks_;k++) {
4117	int iColumn = members_[base+k];
4118	if (lastWeight>=weights_[j]-1.0e-7)
4119	throw CoinError("Weights too close together in OsiLink","infeasibility","OsiLink");
4120	lastWeight = weights_[j];
4121	double value = CoinMax(0.0,solution[iColumn]);
4122	sum += value;
4123	if (value>integerTolerance&&upper[iColumn]) {
4124	// Possibly due to scaling a fixed variable might slip through
4125	if (value>upper[iColumn]+1.0e-8) {
4126	#ifdef OSI_DEBUG
4127	printf("** Variable %d (%d) has value %g and upper bound of %g\n",
4128	iColumn,j,value,upper[iColumn]);
4129	#endif
4130	}
4131	value = CoinMin(value,upper[iColumn]);
4132	weight += weights_[j]*value;
4133	if (firstNonZero<0)
4134	firstNonZero=j;
4135	lastNonZero=j;
4136	}
4137	}
4138	base += numberLinks_;
4139	}
4140	double valueInfeasibility;
4141	whichWay=1;
4142	whichWay_=1;
4143	if (lastNonZero-firstNonZero>=sosType_) {
4144	// find where to branch
4145	assert (sum>0.0);
4146	weight /= sum;
4147	valueInfeasibility = lastNonZero-firstNonZero+1;
4148	valueInfeasibility *= 0.5/((double) numberMembers_);
4149	//#define DISTANCE
4150	#ifdef DISTANCE
4151	assert (sosType_==1); // code up
4152	/* may still be satisfied.
4153	For LOS type 2 we might wish to move coding around
4154	and keep initial info in model_ for speed
4155	*/
4156	int iWhere;
4157	bool possible=false;
4158	for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
4159	if (fabs(weight-weights_[iWhere])<1.0e-8) {
4160	possible=true;
4161	break;
4162	}
4163	}
4164	if (possible) {
4165	// One could move some of this (+ arrays) into model_
4166	const CoinPackedMatrix * matrix = solver->getMatrixByCol();
4167	const double * element = matrix->getMutableElements();
4168	const int * row = matrix->getIndices();
4169	const CoinBigIndex * columnStart = matrix->getVectorStarts();
4170	const int * columnLength = matrix->getVectorLengths();
4171	const double * rowSolution = solver->getRowActivity();
4172	const double * rowLower = solver->getRowLower();
4173	const double * rowUpper = solver->getRowUpper();
4174	int numberRows = matrix->getNumRows();
4175	double * array = new double [numberRows];
4176	CoinZeroN(array,numberRows);
4177	int * which = new int [numberRows];
4178	int n=0;
4179	int base=numberLinks_*firstNonZero;
4180	for (j=firstNonZero;j<=lastNonZero;j++) {
4181	for (int k=0;k<numberLinks_;k++) {
4182	int iColumn = members_[base+k];
4183	double value = CoinMax(0.0,solution[iColumn]);
4184	if (value>integerTolerance&&upper[iColumn]) {
4185	value = CoinMin(value,upper[iColumn]);
4186	for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
4187	int iRow = row[j];
4188	double a = array[iRow];
4189	if (a) {
4190	a += value*element[j];
4191	if (!a)
4192	a = 1.0e-100;
4193	} else {
4194	which[n++]=iRow;
4195	a=value*element[j];
4196	assert (a);
4197	}
4198	array[iRow]=a;
4199	}
4200	}
4201	}
4202	base += numberLinks_;
4203	}
4204	base=numberLinks_*iWhere;
4205	for (int k=0;k<numberLinks_;k++) {
4206	int iColumn = members_[base+k];
4207	const double value = 1.0;
4208	for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
4209	int iRow = row[j];
4210	double a = array[iRow];
4211	if (a) {
4212	a -= value*element[j];
4213	if (!a)
4214	a = 1.0e-100;
4215	} else {
4216	which[n++]=iRow;
4217	a=-value*element[j];
4218	assert (a);
4219	}
4220	array[iRow]=a;
4221	}
4222	}
4223	for (j=0;j<n;j++) {
4224	int iRow = which[j];
4225	// moving to point will increase row solution by this
4226	double distance = array[iRow];
4227	if (distance>1.0e-8) {
4228	if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
4229	possible=false;
4230	break;
4231	}
4232	} else if (distance<-1.0e-8) {
4233	if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
4234	possible=false;
4235	break;
4236	}
4237	}
4238	}
4239	for (j=0;j<n;j++)
4240	array[which[j]]=0.0;
4241	delete [] array;
4242	delete [] which;
4243	if (possible) {
4244	valueInfeasibility=0.0;
4245	printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere);
4246	}
4247	}
4248	#endif
4249	} else {
4250	valueInfeasibility = 0.0; // satisfied
4251	}
4252	infeasibility_=valueInfeasibility;
4253	otherInfeasibility_=1.0-valueInfeasibility;
4254	return valueInfeasibility;
4255	}
4256
4257	// This looks at solution and sets bounds to contain solution
4258	double
4259	OsiOldLink::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
4260	{
4261	int j;
4262	int firstNonZero=-1;
4263	int lastNonZero = -1;
4264	const double * solution = info->solution_;
4265	const double * upper = info->upper_;
4266	double integerTolerance = info->integerTolerance_;
4267	double weight = 0.0;
4268	double sum =0.0;
4269
4270	int base=0;
4271	for (j=0;j<numberMembers_;j++) {
4272	for (int k=0;k<numberLinks_;k++) {
4273	int iColumn = members_[base+k];
4274	double value = CoinMax(0.0,solution[iColumn]);
4275	sum += value;
4276	if (value>integerTolerance&&upper[iColumn]) {
4277	weight += weights_[j]*value;
4278	if (firstNonZero<0)
4279	firstNonZero=j;
4280	lastNonZero=j;
4281	}
4282	}
4283	base += numberLinks_;
4284	}
4285	#ifdef DISTANCE
4286	if (lastNonZero-firstNonZero>sosType_-1) {
4287	/* may still be satisfied.
4288	For LOS type 2 we might wish to move coding around
4289	and keep initial info in model_ for speed
4290	*/
4291	int iWhere;
4292	bool possible=false;
4293	for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
4294	if (fabs(weight-weights_[iWhere])<1.0e-8) {
4295	possible=true;
4296	break;
4297	}
4298	}
4299	if (possible) {
4300	// One could move some of this (+ arrays) into model_
4301	const CoinPackedMatrix * matrix = solver->getMatrixByCol();
4302	const double * element = matrix->getMutableElements();
4303	const int * row = matrix->getIndices();
4304	const CoinBigIndex * columnStart = matrix->getVectorStarts();
4305	const int * columnLength = matrix->getVectorLengths();
4306	const double * rowSolution = solver->getRowActivity();
4307	const double * rowLower = solver->getRowLower();
4308	const double * rowUpper = solver->getRowUpper();
4309	int numberRows = matrix->getNumRows();
4310	double * array = new double [numberRows];
4311	CoinZeroN(array,numberRows);
4312	int * which = new int [numberRows];
4313	int n=0;
4314	int base=numberLinks_*firstNonZero;
4315	for (j=firstNonZero;j<=lastNonZero;j++) {
4316	for (int k=0;k<numberLinks_;k++) {
4317	int iColumn = members_[base+k];
4318	double value = CoinMax(0.0,solution[iColumn]);
4319	if (value>integerTolerance&&upper[iColumn]) {
4320	value = CoinMin(value,upper[iColumn]);
4321	for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
4322	int iRow = row[j];
4323	double a = array[iRow];
4324	if (a) {
4325	a += value*element[j];
4326	if (!a)
4327	a = 1.0e-100;
4328	} else {
4329	which[n++]=iRow;
4330	a=value*element[j];
4331	assert (a);
4332	}
4333	array[iRow]=a;
4334	}
4335	}
4336	}
4337	base += numberLinks_;
4338	}
4339	base=numberLinks_*iWhere;
4340	for (int k=0;k<numberLinks_;k++) {
4341	int iColumn = members_[base+k];
4342	const double value = 1.0;
4343	for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
4344	int iRow = row[j];
4345	double a = array[iRow];
4346	if (a) {
4347	a -= value*element[j];
4348	if (!a)
4349	a = 1.0e-100;
4350	} else {
4351	which[n++]=iRow;
4352	a=-value*element[j];
4353	assert (a);
4354	}
4355	array[iRow]=a;
4356	}
4357	}
4358	for (j=0;j<n;j++) {
4359	int iRow = which[j];
4360	// moving to point will increase row solution by this
4361	double distance = array[iRow];
4362	if (distance>1.0e-8) {
4363	if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
4364	possible=false;
4365	break;
4366	}
4367	} else if (distance<-1.0e-8) {
4368	if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
4369	possible=false;
4370	break;
4371	}
4372	}
4373	}
4374	for (j=0;j<n;j++)
4375	array[which[j]]=0.0;
4376	delete [] array;
4377	delete [] which;
4378	if (possible) {
4379	printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere);
4380	firstNonZero=iWhere;
4381	lastNonZero=iWhere;
4382	}
4383	}
4384	}
4385	#else
4386	assert (lastNonZero-firstNonZero<sosType_) ;
4387	#endif
4388	base=0;
4389	for (j=0;j<firstNonZero;j++) {
4390	for (int k=0;k<numberLinks_;k++) {
4391	int iColumn = members_[base+k];
4392	solver->setColUpper(iColumn,0.0);
4393	}
4394	base += numberLinks_;
4395	}
4396	// skip
4397	base += numberLinks_;
4398	for (j=lastNonZero+1;j<numberMembers_;j++) {
4399	for (int k=0;k<numberLinks_;k++) {
4400	int iColumn = members_[base+k];
4401	solver->setColUpper(iColumn,0.0);
4402	}
4403	base += numberLinks_;
4404	}
4405	// go to coding as in OsiSOS
4406	abort();
4407	return -1.0;
4408	}
4409
4410	// Redoes data when sequence numbers change
4411	void
4412	OsiOldLink::resetSequenceEtc(int numberColumns, const int * originalColumns)
4413	{
4414	int n2=0;
4415	for (int j=0;j<numberMembers_*numberLinks_;j++) {
4416	int iColumn = members_[j];
4417	int i;
4418	for (i=0;i<numberColumns;i++) {
4419	if (originalColumns[i]==iColumn)
4420	break;
4421	}
4422	if (i<numberColumns) {
4423	members_[n2]=i;
4424	weights_[n2++]=weights_[j];
4425	}
4426	}
4427	if (n2<numberMembers_) {
4428	printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2/numberLinks_);
4429	numberMembers_=n2/numberLinks_;
4430	}
4431	}
4432
4433	// Creates a branching object
4434	OsiBranchingObject *
4435	OsiOldLink::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
4436	{
4437	int j;
4438	const double * solution = info->solution_;
4439	double tolerance = info->primalTolerance_;
4440	const double * upper = info->upper_;
4441	int firstNonFixed=-1;
4442	int lastNonFixed=-1;
4443	int firstNonZero=-1;
4444	int lastNonZero = -1;
4445	double weight = 0.0;
4446	double sum =0.0;
4447	int base=0;
4448	for (j=0;j<numberMembers_;j++) {
4449	for (int k=0;k<numberLinks_;k++) {
4450	int iColumn = members_[base+k];
4451	if (upper[iColumn]) {
4452	double value = CoinMax(0.0,solution[iColumn]);
4453	sum += value;
4454	if (firstNonFixed<0)
4455	firstNonFixed=j;
4456	lastNonFixed=j;
4457	if (value>tolerance) {
4458	weight += weights_[j]*value;
4459	if (firstNonZero<0)
4460	firstNonZero=j;
4461	lastNonZero=j;
4462	}
4463	}
4464	}
4465	base += numberLinks_;
4466	}
4467	assert (lastNonZero-firstNonZero>=sosType_) ;
4468	// find where to branch
4469	assert (sum>0.0);
4470	weight /= sum;
4471	int iWhere;
4472	double separator=0.0;
4473	for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
4474	if (weight<weights_[iWhere+1])
4475	break;
4476	if (sosType_==1) {
4477	// SOS 1
4478	separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
4479	} else {
4480	// SOS 2
4481	if (iWhere==firstNonFixed)
4482	iWhere++;;
4483	if (iWhere==lastNonFixed-1)
4484	iWhere = lastNonFixed-2;
4485	separator = weights_[iWhere+1];
4486	}
4487	// create object
4488	OsiBranchingObject * branch;
4489	branch = new OsiOldLinkBranchingObject(solver,this,way,separator);
4490	return branch;
4491	}
4492	OsiOldLinkBranchingObject::OsiOldLinkBranchingObject()
4493	:OsiSOSBranchingObject()
4494	{
4495	}
4496
4497	// Useful constructor
4498	OsiOldLinkBranchingObject::OsiOldLinkBranchingObject (OsiSolverInterface * solver,
4499	const OsiOldLink * set,
4500	int way ,
4501	double separator)
4502	:OsiSOSBranchingObject(solver,set,way,separator)
4503	{
4504	}
4505
4506	// Copy constructor
4507	OsiOldLinkBranchingObject::OsiOldLinkBranchingObject ( const OsiOldLinkBranchingObject & rhs) :OsiSOSBranchingObject(rhs)
4508	{
4509	}
4510
4511	// Assignment operator
4512	OsiOldLinkBranchingObject &
4513	OsiOldLinkBranchingObject::operator=( const OsiOldLinkBranchingObject& rhs)
4514	{
4515	if (this != &rhs) {
4516	OsiSOSBranchingObject::operator=(rhs);
4517	}
4518	return *this;
4519	}
4520	OsiBranchingObject *
4521	OsiOldLinkBranchingObject::clone() const
4522	{
4523	return (new OsiOldLinkBranchingObject(*this));
4524	}
4525
4526
4527	// Destructor
4528	OsiOldLinkBranchingObject::~OsiOldLinkBranchingObject ()
4529	{
4530	}
4531	double
4532	OsiOldLinkBranchingObject::branch(OsiSolverInterface * solver)
4533	{
4534	const OsiOldLink * set =
4535	dynamic_cast <const OsiOldLink *>(originalObject_) ;
4536	assert (set);
4537	int way = (!branchIndex_) ? (2firstBranch_-1) : -(2firstBranch_-1);
4538	branchIndex_++;
4539	int numberMembers = set->numberMembers();
4540	const int * which = set->members();
4541	const double * weights = set->weights();
4542	int numberLinks = set->numberLinks();
4543	//const double * lower = info->lower_;
4544	//const double * upper = solver->getColUpper();
4545	// *** for way - up means fix all those in down section
4546	if (way<0) {
4547	int i;
4548	for ( i=0;i<numberMembers;i++) {
4549	if (weights[i] > value_)
4550	break;
4551	}
4552	assert (i<numberMembers);
4553	int base=i*numberLinks;;
4554	for (;i<numberMembers;i++) {
4555	for (int k=0;k<numberLinks;k++) {
4556	int iColumn = which[base+k];
4557	solver->setColUpper(iColumn,0.0);
4558	}
4559	base += numberLinks;
4560	}
4561	} else {
4562	int i;
4563	int base=0;
4564	for ( i=0;i<numberMembers;i++) {
4565	if (weights[i] >= value_) {
4566	break;
4567	} else {
4568	for (int k=0;k<numberLinks;k++) {
4569	int iColumn = which[base+k];
4570	solver->setColUpper(iColumn,0.0);
4571	}
4572	base += numberLinks;
4573	}
4574	}
4575	assert (i<numberMembers);
4576	}
4577	return 0.0;
4578	}
4579	// Print what would happen
4580	void
4581	OsiOldLinkBranchingObject::print(const OsiSolverInterface * solver)
4582	{
4583	const OsiOldLink * set =
4584	dynamic_cast <const OsiOldLink *>(originalObject_) ;
4585	assert (set);
4586	int way = (!branchIndex_) ? (2firstBranch_-1) : -(2firstBranch_-1);
4587	int numberMembers = set->numberMembers();
4588	int numberLinks = set->numberLinks();
4589	const double * weights = set->weights();
4590	const int * which = set->members();
4591	const double * upper = solver->getColUpper();
4592	int first=numberMembers;
4593	int last=-1;
4594	int numberFixed=0;
4595	int numberOther=0;
4596	int i;
4597	int base=0;
4598	for ( i=0;i<numberMembers;i++) {
4599	for (int k=0;k<numberLinks;k++) {
4600	int iColumn = which[base+k];
4601	double bound = upper[iColumn];
4602	if (bound) {
4603	first = CoinMin(first,i);
4604	last = CoinMax(last,i);
4605	}
4606	}
4607	base += numberLinks;
4608	}
4609	// *** for way - up means fix all those in down section
4610	base=0;
4611	if (way<0) {
4612	printf("SOS Down");
4613	for ( i=0;i<numberMembers;i++) {
4614	if (weights[i] > value_)
4615	break;
4616	for (int k=0;k<numberLinks;k++) {
4617	int iColumn = which[base+k];
4618	double bound = upper[iColumn];
4619	if (bound)
4620	numberOther++;
4621	}
4622	base += numberLinks;
4623	}
4624	assert (i<numberMembers);
4625	for (;i<numberMembers;i++) {
4626	for (int k=0;k<numberLinks;k++) {
4627	int iColumn = which[base+k];
4628	double bound = upper[iColumn];
4629	if (bound)
4630	numberFixed++;
4631	}
4632	base += numberLinks;
4633	}
4634	} else {
4635	printf("SOS Up");
4636	for ( i=0;i<numberMembers;i++) {
4637	if (weights[i] >= value_)
4638	break;
4639	for (int k=0;k<numberLinks;k++) {
4640	int iColumn = which[base+k];
4641	double bound = upper[iColumn];
4642	if (bound)
4643	numberFixed++;
4644	}
4645	base += numberLinks;
4646	}
4647	assert (i<numberMembers);
4648	for (;i<numberMembers;i++) {
4649	for (int k=0;k<numberLinks;k++) {
4650	int iColumn = which[base+k];
4651	double bound = upper[iColumn];
4652	if (bound)
4653	numberOther++;
4654	}
4655	base += numberLinks;
4656	}
4657	}
4658	assert ((numberFixed%numberLinks)==0);
4659	assert ((numberOther%numberLinks)==0);
4660	printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
4661	value_,first,weights[first],last,weights[last],numberFixed/numberLinks,
4662	numberOther/numberLinks);
4663	}
4664	// Default Constructor
4665	OsiBiLinear::OsiBiLinear ()
4666	: OsiObject2(),
4667	coefficient_(0.0),
4668	xMeshSize_(0.0),
4669	yMeshSize_(0.0),
4670	xSatisfied_(1.0e-6),
4671	ySatisfied_(1.0e-6),
4672	xOtherSatisfied_(0.0),
4673	yOtherSatisfied_(0.0),
4674	xySatisfied_(1.0e-6),
4675	xyBranchValue_(0.0),
4676	xColumn_(-1),
4677	yColumn_(-1),
4678	firstLambda_(-1),
4679	branchingStrategy_(0),
4680	boundType_(0),
4681	xRow_(-1),
4682	yRow_(-1),
4683	xyRow_(-1),
4684	convexity_(-1),
4685	numberExtraRows_(0),
4686	multiplier_(NULL),
4687	extraRow_(NULL),
4688	chosen_(-1)
4689	{
4690	}
4691
4692	// Useful constructor
4693	OsiBiLinear::OsiBiLinear (OsiSolverInterface * solver, int xColumn,
4694	int yColumn, int xyRow, double coefficient,
4695	double xMesh, double yMesh,
4696	int numberExistingObjects,const OsiObject ** objects )
4697	: OsiObject2(),
4698	coefficient_(coefficient),
4699	xMeshSize_(xMesh),
4700	yMeshSize_(yMesh),
4701	xSatisfied_(1.0e-6),
4702	ySatisfied_(1.0e-6),
4703	xOtherSatisfied_(0.0),
4704	yOtherSatisfied_(0.0),
4705	xySatisfied_(1.0e-6),
4706	xyBranchValue_(0.0),
4707	xColumn_(xColumn),
4708	yColumn_(yColumn),
4709	firstLambda_(-1),
4710	branchingStrategy_(0),
4711	boundType_(0),
4712	xRow_(-1),
4713	yRow_(-1),
4714	xyRow_(xyRow),
4715	convexity_(-1),
4716	numberExtraRows_(0),
4717	multiplier_(NULL),
4718	extraRow_(NULL),
4719	chosen_(-1)
4720	{
4721	double columnLower[4];
4722	double columnUpper[4];
4723	double objective[4];
4724	double rowLower[3];
4725	double rowUpper[3];
4726	CoinBigIndex starts[5];
4727	int index[16];
4728	double element[16];
4729	int i;
4730	starts[0]=0;
4731	// rows
4732	int numberRows = solver->getNumRows();
4733	// convexity
4734	rowLower[0]=1.0;
4735	rowUpper[0]=1.0;
4736	convexity_ = numberRows;
4737	starts[1]=0;
4738	// x
4739	rowLower[1]=0.0;
4740	rowUpper[1]=0.0;
4741	index[0]=xColumn_;
4742	element[0]=-1.0;
4743	xRow_ = numberRows+1;
4744	starts[2]=1;
4745	int nAdd=2;
4746	if (xColumn_!=yColumn_) {
4747	rowLower[2]=0.0;
4748	rowUpper[2]=0.0;
4749	index[1]=yColumn;
4750	element[1]=-1.0;
4751	nAdd=3;
4752	yRow_ = numberRows+2;
4753	starts[3]=2;
4754	} else {
4755	yRow_=-1;
4756	branchingStrategy_=1;
4757	}
4758	// may be objective
4759	assert (xyRow_>=-1);
4760	solver->addRows(nAdd,starts,index,element,rowLower,rowUpper);
4761	int n=0;
4762	// order is LxLy, LxUy, UxLy and UxUy
4763	firstLambda_ = solver->getNumCols();
4764	// bit sloppy as theoretically could be infeasible but otherwise need to do more work
4765	double xB[2];
4766	double yB[2];
4767	const double * lower = solver->getColLower();
4768	const double * upper = solver->getColUpper();
4769	xB[0]=lower[xColumn_];
4770	xB[1]=upper[xColumn_];
4771	yB[0]=lower[yColumn_];
4772	yB[1]=upper[yColumn_];
4773	if (xMeshSize_!=floor(xMeshSize_)) {
4774	// not integral
4775	xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_);
4776	if (!yMeshSize_) {
4777	xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1])));
4778	}
4779	}
4780	if (yMeshSize_!=floor(yMeshSize_)) {
4781	// not integral
4782	ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_);
4783	if (!xMeshSize_) {
4784	xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1])));
4785	}
4786	}
4787	// adjust
4788	double distance;
4789	double steps;
4790	if (xMeshSize_) {
4791	distance = xB[1]-xB[0];
4792	steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
4793	distance = xB[0]+xMeshSize_*steps;
4794	if (fabs(xB[1]-distance)>xSatisfied_) {
4795	printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance);
4796	//double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
4797	//printf("xSatisfied increased to %g\n",newValue);
4798	//xSatisfied_ = newValue;
4799	//xB[1]=distance;
4800	//solver->setColUpper(xColumn_,distance);
4801	}
4802	}
4803	if (yMeshSize_) {
4804	distance = yB[1]-yB[0];
4805	steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
4806	distance = yB[0]+yMeshSize_*steps;
4807	if (fabs(yB[1]-distance)>ySatisfied_) {
4808	printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance);
4809	//double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
4810	//printf("ySatisfied increased to %g\n",newValue);
4811	//ySatisfied_ = newValue;
4812	//yB[1]=distance;
4813	//solver->setColUpper(yColumn_,distance);
4814	}
4815	}
4816	for (i=0;i<4;i++) {
4817	double x = (i<2) ? xB[0] : xB[1];
4818	double y = ((i&1)==0) ? yB[0] : yB[1];
4819	columnLower[i]=0.0;
4820	columnUpper[i]=2.0;
4821	objective[i]=0.0;
4822	double value;
4823	// xy
4824	value=coefficient_xy;
4825	if (xyRow_>=0) {
4826	if (fabs(value)<1.0e-19)
4827	value = 1.0e-19;
4828	element[n]=value;
4829	index[n++]=xyRow_;
4830	} else {
4831	objective[i]=value;
4832	}
4833	// convexity
4834	value=1.0;
4835	element[n]=value;
4836	index[n++]=0+numberRows;
4837	// x
4838	value=x;
4839	if (fabs(value)<1.0e-19)
4840	value = 1.0e-19;
4841	element[n]=value;
4842	index[n++]=1+numberRows;
4843	if (xColumn_!=yColumn_) {
4844	// y
4845	value=y;
4846	if (fabs(value)<1.0e-19)
4847	value = 1.0e-19;
4848	element[n]=value;
4849	index[n++]=2+numberRows;
4850	}
4851	starts[i+1]=n;
4852	}
4853	solver->addCols(4,starts,index,element,columnLower,columnUpper,objective);
4854	// At least one has to have a mesh
4855	if (!xMeshSize_&&(!yMeshSize_\|\|yRow_<0)) {
4856	printf("one of x and y must have a mesh size\n");
4857	abort();
4858	} else if (yRow_>=0) {
4859	if (!xMeshSize_)
4860	branchingStrategy_ = 2;
4861	else if (!yMeshSize_)
4862	branchingStrategy_ = 1;
4863	}
4864	// Now add constraints to link in x and or y to existing ones.
4865	bool xDone=false;
4866	bool yDone=false;
4867	// order is LxLy, LxUy, UxLy and UxUy
4868	for (i=numberExistingObjects-1;i>=0;i--) {
4869	const OsiObject * obj = objects[i];
4870	const OsiBiLinear * obj2 =
4871	dynamic_cast <const OsiBiLinear *>(obj) ;
4872	if (obj2) {
4873	if (xColumn_==obj2->xColumn_&&!xDone) {
4874	// make sure y equal
4875	double rhs=0.0;
4876	CoinBigIndex starts[2];
4877	int index[4];
4878	double element[4]= {1.0,1.0,-1.0,-1.0};
4879	starts[0]=0;
4880	starts[1]=4;
4881	index[0]=firstLambda_+0;
4882	index[1]=firstLambda_+1;
4883	index[2]=obj2->firstLambda_+0;
4884	index[3]=obj2->firstLambda_+1;
4885	solver->addRows(1,starts,index,element,&rhs,&rhs);
4886	xDone=true;
4887	}
4888	if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) {
4889	// make sure x equal
4890	double rhs=0.0;
4891	CoinBigIndex starts[2];
4892	int index[4];
4893	double element[4]= {1.0,1.0,-1.0,-1.0};
4894	starts[0]=0;
4895	starts[1]=4;
4896	index[0]=firstLambda_+0;
4897	index[1]=firstLambda_+2;
4898	index[2]=obj2->firstLambda_+0;
4899	index[3]=obj2->firstLambda_+2;
4900	solver->addRows(1,starts,index,element,&rhs,&rhs);
4901	yDone=true;
4902	}
4903	}
4904	}
4905	}
4906	// Set sizes and other stuff
4907	void
4908	OsiBiLinear::setMeshSizes(const OsiSolverInterface * solver, double x, double y)
4909	{
4910	xMeshSize_ = x;
4911	yMeshSize_ = y;
4912	double xB[2];
4913	double yB[2];
4914	const double * lower = solver->getColLower();
4915	const double * upper = solver->getColUpper();
4916	xB[0]=lower[xColumn_];
4917	xB[1]=upper[xColumn_];
4918	yB[0]=lower[yColumn_];
4919	yB[1]=upper[yColumn_];
4920	if (xMeshSize_!=floor(xMeshSize_)) {
4921	// not integral
4922	xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_);
4923	if (!yMeshSize_) {
4924	xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1])));
4925	}
4926	}
4927	if (yMeshSize_!=floor(yMeshSize_)) {
4928	// not integral
4929	ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_);
4930	if (!xMeshSize_) {
4931	xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1])));
4932	}
4933	}
4934	}
4935	// Useful constructor
4936	OsiBiLinear::OsiBiLinear (CoinModel * coinModel, int xColumn,
4937	int yColumn, int xyRow, double coefficient,
4938	double xMesh, double yMesh,
4939	int numberExistingObjects,const OsiObject ** objects )
4940	: OsiObject2(),
4941	coefficient_(coefficient),
4942	xMeshSize_(xMesh),
4943	yMeshSize_(yMesh),
4944	xSatisfied_(1.0e-6),
4945	ySatisfied_(1.0e-6),
4946	xOtherSatisfied_(0.0),
4947	yOtherSatisfied_(0.0),
4948	xySatisfied_(1.0e-6),
4949	xyBranchValue_(0.0),
4950	xColumn_(xColumn),
4951	yColumn_(yColumn),
4952	firstLambda_(-1),
4953	branchingStrategy_(0),
4954	boundType_(0),
4955	xRow_(-1),
4956	yRow_(-1),
4957	xyRow_(xyRow),
4958	convexity_(-1),
4959	numberExtraRows_(0),
4960	multiplier_(NULL),
4961	extraRow_(NULL),
4962	chosen_(-1)
4963	{
4964	double columnLower[4];
4965	double columnUpper[4];
4966	double objective[4];
4967	double rowLower[3];
4968	double rowUpper[3];
4969	CoinBigIndex starts[5];
4970	int index[16];
4971	double element[16];
4972	int i;
4973	starts[0]=0;
4974	// rows
4975	int numberRows = coinModel->numberRows();
4976	// convexity
4977	rowLower[0]=1.0;
4978	rowUpper[0]=1.0;
4979	convexity_ = numberRows;
4980	starts[1]=0;
4981	// x
4982	rowLower[1]=0.0;
4983	rowUpper[1]=0.0;
4984	index[0]=xColumn_;
4985	element[0]=-1.0;
4986	xRow_ = numberRows+1;
4987	starts[2]=1;
4988	int nAdd=2;
4989	if (xColumn_!=yColumn_) {
4990	rowLower[2]=0.0;
4991	rowUpper[2]=0.0;
4992	index[1]=yColumn;
4993	element[1]=-1.0;
4994	nAdd=3;
4995	yRow_ = numberRows+2;
4996	starts[3]=2;
4997	} else {
4998	yRow_=-1;
4999	branchingStrategy_=1;
5000	}
5001	// may be objective
5002	assert (xyRow_>=-1);
5003	for (i=0;i<nAdd;i++) {
5004	CoinBigIndex iStart = starts[i];
5005	coinModel->addRow(starts[i+1]-iStart,index+iStart,element+iStart,rowLower[i],rowUpper[i]);
5006	}
5007	int n=0;
5008	// order is LxLy, LxUy, UxLy and UxUy
5009	firstLambda_ = coinModel->numberColumns();
5010	// bit sloppy as theoretically could be infeasible but otherwise need to do more work
5011	double xB[2];
5012	double yB[2];
5013	const double * lower = coinModel->columnLowerArray();
5014	const double * upper = coinModel->columnUpperArray();
5015	xB[0]=lower[xColumn_];
5016	xB[1]=upper[xColumn_];
5017	yB[0]=lower[yColumn_];
5018	yB[1]=upper[yColumn_];
5019	if (xMeshSize_!=floor(xMeshSize_)) {
5020	// not integral
5021	xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_);
5022	if (!yMeshSize_) {
5023	xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1])));
5024	}
5025	}
5026	if (yMeshSize_!=floor(yMeshSize_)) {
5027	// not integral
5028	ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_);
5029	if (!xMeshSize_) {
5030	xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1])));
5031	}
5032	}
5033	// adjust
5034	double distance;
5035	double steps;
5036	if (xMeshSize_) {
5037	distance = xB[1]-xB[0];
5038	steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
5039	distance = xB[0]+xMeshSize_*steps;
5040	if (fabs(xB[1]-distance)>xSatisfied_) {
5041	printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance);
5042	//double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
5043	//printf("xSatisfied increased to %g\n",newValue);
5044	//xSatisfied_ = newValue;
5045	//xB[1]=distance;
5046	//coinModel->setColUpper(xColumn_,distance);
5047	}
5048	}
5049	if (yMeshSize_) {
5050	distance = yB[1]-yB[0];
5051	steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
5052	distance = yB[0]+yMeshSize_*steps;
5053	if (fabs(yB[1]-distance)>ySatisfied_) {
5054	printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance);
5055	//double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
5056	//printf("ySatisfied increased to %g\n",newValue);
5057	//ySatisfied_ = newValue;
5058	//yB[1]=distance;
5059	//coinModel->setColUpper(yColumn_,distance);
5060	}
5061	}
5062	for (i=0;i<4;i++) {
5063	double x = (i<2) ? xB[0] : xB[1];
5064	double y = ((i&1)==0) ? yB[0] : yB[1];
5065	columnLower[i]=0.0;
5066	columnUpper[i]=2.0;
5067	objective[i]=0.0;
5068	double value;
5069	// xy
5070	value=coefficient_xy;
5071	if (xyRow_>=0) {
5072	if (fabs(value)<1.0e-19)
5073	value = 1.0e-19;
5074	element[n]=value;
5075	index[n++]=xyRow_;
5076	} else {
5077	objective[i]=value;
5078	}
5079	// convexity
5080	value=1.0;
5081	element[n]=value;
5082	index[n++]=0+numberRows;
5083	// x
5084	value=x;
5085	if (fabs(value)<1.0e-19)
5086	value = 1.0e-19;
5087	element[n]=value;
5088	index[n++]=1+numberRows;
5089	if (xColumn_!=yColumn_) {
5090	// y
5091	value=y;
5092	if (fabs(value)<1.0e-19)
5093	value = 1.0e-19;
5094	element[n]=value;
5095	index[n++]=2+numberRows;
5096	}
5097	starts[i+1]=n;
5098	}
5099	for (i=0;i<4;i++) {
5100	CoinBigIndex iStart = starts[i];
5101	coinModel->addColumn(starts[i+1]-iStart,index+iStart,element+iStart,columnLower[i],
5102	columnUpper[i],objective[i]);
5103	}
5104	// At least one has to have a mesh
5105	if (!xMeshSize_&&(!yMeshSize_\|\|yRow_<0)) {
5106	printf("one of x and y must have a mesh size\n");
5107	abort();
5108	} else if (yRow_>=0) {
5109	if (!xMeshSize_)
5110	branchingStrategy_ = 2;
5111	else if (!yMeshSize_)
5112	branchingStrategy_ = 1;
5113	}
5114	// Now add constraints to link in x and or y to existing ones.
5115	bool xDone=false;
5116	bool yDone=false;
5117	// order is LxLy, LxUy, UxLy and UxUy
5118	for (i=numberExistingObjects-1;i>=0;i--) {
5119	const OsiObject * obj = objects[i];
5120	const OsiBiLinear * obj2 =
5121	dynamic_cast <const OsiBiLinear *>(obj) ;
5122	if (obj2) {
5123	if (xColumn_==obj2->xColumn_&&!xDone) {
5124	// make sure y equal
5125	double rhs=0.0;
5126	CoinBigIndex starts[2];
5127	int index[4];
5128	double element[4]= {1.0,1.0,-1.0,-1.0};
5129	starts[0]=0;
5130	starts[1]=4;
5131	index[0]=firstLambda_+0;
5132	index[1]=firstLambda_+1;
5133	index[2]=obj2->firstLambda_+0;
5134	index[3]=obj2->firstLambda_+1;
5135	coinModel->addRow(4,index,element,rhs,rhs);
5136	xDone=true;
5137	}
5138	if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) {
5139	// make sure x equal
5140	double rhs=0.0;
5141	CoinBigIndex starts[2];
5142	int index[4];
5143	double element[4]= {1.0,1.0,-1.0,-1.0};
5144	starts[0]=0;
5145	starts[1]=4;
5146	index[0]=firstLambda_+0;
5147	index[1]=firstLambda_+2;
5148	index[2]=obj2->firstLambda_+0;
5149	index[3]=obj2->firstLambda_+2;
5150	coinModel->addRow(4,index,element,rhs,rhs);
5151	yDone=true;
5152	}
5153	}
5154	}
5155	}
5156
5157	// Copy constructor
5158	OsiBiLinear::OsiBiLinear ( const OsiBiLinear & rhs)
5159	:OsiObject2(rhs),
5160	coefficient_(rhs.coefficient_),
5161	xMeshSize_(rhs.xMeshSize_),
5162	yMeshSize_(rhs.yMeshSize_),
5163	xSatisfied_(rhs.xSatisfied_),
5164	ySatisfied_(rhs.ySatisfied_),
5165	xOtherSatisfied_(rhs.xOtherSatisfied_),
5166	yOtherSatisfied_(rhs.yOtherSatisfied_),
5167	xySatisfied_(rhs.xySatisfied_),
5168	xyBranchValue_(rhs.xyBranchValue_),
5169	xColumn_(rhs.xColumn_),
5170	yColumn_(rhs.yColumn_),
5171	firstLambda_(rhs.firstLambda_),
5172	branchingStrategy_(rhs.branchingStrategy_),
5173	boundType_(rhs.boundType_),
5174	xRow_(rhs.xRow_),
5175	yRow_(rhs.yRow_),
5176	xyRow_(rhs.xyRow_),
5177	convexity_(rhs.convexity_),
5178	numberExtraRows_(rhs.numberExtraRows_),
5179	multiplier_(NULL),
5180	extraRow_(NULL),
5181	chosen_(rhs.chosen_)
5182	{
5183	if (numberExtraRows_) {
5184	multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
5185	extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
5186	}
5187	}
5188
5189	// Clone
5190	OsiObject *
5191	OsiBiLinear::clone() const
5192	{
5193	return new OsiBiLinear(*this);
5194	}
5195
5196	// Assignment operator
5197	OsiBiLinear &
5198	OsiBiLinear::operator=( const OsiBiLinear& rhs)
5199	{
5200	if (this!=&rhs) {
5201	OsiObject2::operator=(rhs);
5202	coefficient_ = rhs.coefficient_;
5203	xMeshSize_ = rhs.xMeshSize_;
5204	yMeshSize_ = rhs.yMeshSize_;
5205	xSatisfied_ = rhs.xSatisfied_;
5206	ySatisfied_ = rhs.ySatisfied_;
5207	xOtherSatisfied_ = rhs.xOtherSatisfied_;
5208	yOtherSatisfied_ = rhs.yOtherSatisfied_;
5209	xySatisfied_ = rhs.xySatisfied_;
5210	xyBranchValue_ = rhs.xyBranchValue_;
5211	xColumn_ = rhs.xColumn_;
5212	yColumn_ = rhs.yColumn_;
5213	firstLambda_ = rhs.firstLambda_;
5214	branchingStrategy_ = rhs.branchingStrategy_;
5215	boundType_ = rhs.boundType_;
5216	xRow_ = rhs.xRow_;
5217	yRow_ = rhs.yRow_;
5218	xyRow_ = rhs.xyRow_;
5219	convexity_ = rhs.convexity_;
5220	numberExtraRows_ = rhs.numberExtraRows_;
5221	delete [] multiplier_;
5222	delete [] extraRow_;
5223	if (numberExtraRows_) {
5224	multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
5225	extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
5226	} else {
5227	multiplier_ = NULL;
5228	extraRow_ = NULL;
5229	}
5230	chosen_ = rhs.chosen_;
5231	}
5232	return *this;
5233	}
5234
5235	// Destructor
5236	OsiBiLinear::~OsiBiLinear ()
5237	{
5238	delete [] multiplier_;
5239	delete [] extraRow_;
5240	}
5241	// Adds in data for extra row with variable coefficients
5242	void
5243	OsiBiLinear::addExtraRow(int row, double multiplier)
5244	{
5245	int * tempI = new int [numberExtraRows_+1];
5246	double * tempD = new double [numberExtraRows_+1];
5247	memcpy(tempI,extraRow_,numberExtraRows_*sizeof(int));
5248	memcpy(tempD,multiplier_,numberExtraRows_*sizeof(double));
5249	tempI[numberExtraRows_]=row;
5250	tempD[numberExtraRows_]=multiplier;
5251	if (numberExtraRows_)
5252	assert (row>tempI[numberExtraRows_-1]);
5253	numberExtraRows_++;
5254	delete [] extraRow_;
5255	extraRow_ = tempI;
5256	delete [] multiplier_;
5257	multiplier_ = tempD;
5258	}
5259	static bool testCoarse=true;
5260	// Infeasibility - large is 0.5
5261	double
5262	OsiBiLinear::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
5263	{
5264	// order is LxLy, LxUy, UxLy and UxUy
5265	double xB[2];
5266	double yB[2];
5267	xB[0]=info->lower_[xColumn_];
5268	xB[1]=info->upper_[xColumn_];
5269	yB[0]=info->lower_[yColumn_];
5270	yB[1]=info->upper_[yColumn_];
5271	#if 0
5272	if (info->lower_[1]<=43.0&&info->upper_[1]>=43.0) {
5273	if (info->lower_[4]<=49.0&&info->upper_[4]>=49.0) {
5274	if (info->lower_[2]<=16.0&&info->upper_[2]>=16.0) {
5275	if (info->lower_[3]<=19.0&&info->upper_[3]>=19.0) {
5276	printf("feas %g %g %g %g p %g t %g\n",
5277	info->solution_[1],
5278	info->solution_[2],
5279	info->solution_[3],
5280	info->solution_[4],
5281	info->solution_[0],
5282	info->solution_[5]);
5283	}
5284	}
5285	}
5286	}
5287	#endif
5288	double x = info->solution_[xColumn_];
5289	x = CoinMax(x,xB[0]);
5290	x = CoinMin(x,xB[1]);
5291	double y = info->solution_[yColumn_];
5292	y = CoinMax(y,yB[0]);
5293	y = CoinMin(y,yB[1]);
5294	int j;
5295	#ifndef NDEBUG
5296	double xLambda = 0.0;
5297	double yLambda = 0.0;
5298	if ((branchingStrategy_&4)==0) {
5299	for (j=0;j<4;j++) {
5300	int iX = j>>1;
5301	int iY = j&1;
5302	xLambda += xB[iX]*info->solution_[firstLambda_+j];
5303	if (yRow_>=0)
5304	yLambda += yB[iY]*info->solution_[firstLambda_+j];
5305	}
5306	} else {
5307	const double * element = info->elementByColumn_;
5308	const int * row = info->row_;
5309	const CoinBigIndex * columnStart = info->columnStart_;
5310	const int * columnLength = info->columnLength_;
5311	for (j=0;j<4;j++) {
5312	int iColumn = firstLambda_+j;
5313	int iStart = columnStart[iColumn];
5314	int iEnd = iStart + columnLength[iColumn];
5315	int k=iStart;
5316	double sol = info->solution_[iColumn];
5317	for (;k<iEnd;k++) {
5318	if (xRow_==row[k])
5319	xLambda += element[k]*sol;
5320	if (yRow_==row[k])
5321	yLambda += element[k]*sol;
5322	}
5323	}
5324	}
5325	assert (fabs(x-xLambda)<1.0e-1);
5326	if (yRow_>=0)
5327	assert (fabs(y-yLambda)<1.0e-1);
5328	#endif
5329	// If x or y not satisfied then branch on that
5330	double distance;
5331	double steps;
5332	bool xSatisfied;
5333	double xNew=xB[0];
5334	if (xMeshSize_) {
5335	if (x<0.5*(xB[0]+xB[1])) {
5336	distance = x-xB[0];
5337	steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
5338	xNew = xB[0]+steps*xMeshSize_;
5339	assert (xNew<=xB[1]+xSatisfied_);
5340	xSatisfied = (fabs(xNew-x)<xSatisfied_);
5341	} else {
5342	distance = xB[1]-x;
5343	steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
5344	xNew = xB[1]-steps*xMeshSize_;
5345	assert (xNew>=xB[0]-xSatisfied_);
5346	xSatisfied = (fabs(xNew-x)<xSatisfied_);
5347	}
5348	// but if first coarse grid then only if gap small
5349	if (testCoarse&&(branchingStrategy_&8)!=0&&xSatisfied&&
5350	xB[1]-xB[0]>=xMeshSize_) {
5351	// but allow if fine grid would allow
5352	if (fabs(xNew-x)>=xOtherSatisfied_&&fabs(yB[0]-y)>yOtherSatisfied_
5353	&&fabs(yB[1]-y)>yOtherSatisfied_) {
5354	xNew = 0.5*(xB[0]+xB[1]);
5355	x = xNew;
5356	xSatisfied=false;
5357	}
5358	}
5359	} else {
5360	xSatisfied=true;
5361	}
5362	bool ySatisfied;
5363	double yNew=yB[0];
5364	if (yMeshSize_) {
5365	if (y<0.5*(yB[0]+yB[1])) {
5366	distance = y-yB[0];
5367	steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
5368	yNew = yB[0]+steps*yMeshSize_;
5369	assert (yNew<=yB[1]+ySatisfied_);
5370	ySatisfied = (fabs(yNew-y)<ySatisfied_);
5371	} else {
5372	distance = yB[1]-y;
5373	steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
5374	yNew = yB[1]-steps*yMeshSize_;
5375	assert (yNew>=yB[0]-ySatisfied_);
5376	ySatisfied = (fabs(yNew-y)<ySatisfied_);
5377	}
5378	// but if first coarse grid then only if gap small
5379	if (testCoarse&&(branchingStrategy_&8)!=0&&ySatisfied&&
5380	yB[1]-yB[0]>=yMeshSize_) {
5381	// but allow if fine grid would allow
5382	if (fabs(yNew-y)>=yOtherSatisfied_&&fabs(xB[0]-x)>xOtherSatisfied_
5383	&&fabs(xB[1]-x)>xOtherSatisfied_) {
5384	yNew = 0.5*(yB[0]+yB[1]);
5385	y = yNew;
5386	ySatisfied=false;
5387	}
5388	}
5389	} else {
5390	ySatisfied=true;
5391	}
5392	/* There are several possibilities
5393	1 - one or both are unsatisfied and branching strategy tells us what to do
5394	2 - both are unsatisfied and branching strategy is 0
5395	3 - both are satisfied but xy is not
5396	3a one has bounds within satisfied_ - other does not
5397	(or neither have but branching strategy tells us what to do)
5398	3b neither do - and branching strategy does not tell us
5399	3c both do - treat as feasible knowing another copy of object will fix
5400	4 - both are satisfied and xy is satisfied - as 3c
5401	*/
5402	chosen_=-1;
5403	xyBranchValue_=COIN_DBL_MAX;
5404	whichWay_=0;
5405	double xyTrue = x*y;
5406	double xyLambda = 0.0;
5407	if ((branchingStrategy_&4)==0) {
5408