Context Navigation

source: trunk/Cbc/src/CbcSOS.cpp @ 1433

Last change on this file since 1433 was 1432, checked in by bjarni, 10 years ago
Added extra return at end of each source file where needed, to remove possible linefeed conflicts (NightlyBuild? errors)
File size: 34.8 KB

Line
1	// Edwin 11/9/2009-- carved out of CbcBranchActual
2
3	#if defined(_MSC_VER)
4	// Turn off compiler warning about long names
5	# pragma warning(disable:4786)
6	#endif
7	#include <cassert>
8	#include <cstdlib>
9	#include <cmath>
10	#include <cfloat>
11	//#define CBC_DEBUG
12
13	#include "CoinTypes.hpp"
14	#include "OsiSolverInterface.hpp"
15	#include "OsiSolverBranch.hpp"
16	#include "CbcModel.hpp"
17	#include "CbcMessage.hpp"
18	#include "CbcSOS.hpp"
19	#include "CbcBranchActual.hpp"
20	#include "CoinSort.hpp"
21	#include "CoinError.hpp"
22
23	//##############################################################################
24
25	// Default Constructor
26	CbcSOS::CbcSOS ()
27	: CbcObject(),
28	members_(NULL),
29	weights_(NULL),
30	shadowEstimateDown_(1.0),
31	shadowEstimateUp_(1.0),
32	downDynamicPseudoRatio_(0.0),
33	upDynamicPseudoRatio_(0.0),
34	numberTimesDown_(0),
35	numberTimesUp_(0),
36	numberMembers_(0),
37	sosType_(-1),
38	integerValued_(false)
39	{
40	}
41
42	// Useful constructor (which are indices)
43	CbcSOS::CbcSOS (CbcModel * model, int numberMembers,
44	const int * which, const double * weights, int identifier, int type)
45	: CbcObject(model),
46	shadowEstimateDown_(1.0),
47	shadowEstimateUp_(1.0),
48	downDynamicPseudoRatio_(0.0),
49	upDynamicPseudoRatio_(0.0),
50	numberTimesDown_(0),
51	numberTimesUp_(0),
52	numberMembers_(numberMembers),
53	sosType_(type)
54	{
55	id_ = identifier;
56	integerValued_ = type == 1;
57	if (integerValued_) {
58	// check all members integer
59	OsiSolverInterface * solver = model->solver();
60	if (solver) {
61	for (int i = 0; i < numberMembers_; i++) {
62	if (!solver->isInteger(which[i]))
63	integerValued_ = false;
64	}
65	} else {
66	// can't tell
67	integerValued_ = false;
68	}
69	}
70	if (numberMembers_) {
71	members_ = new int[numberMembers_];
72	weights_ = new double[numberMembers_];
73	memcpy(members_, which, numberMembers_*sizeof(int));
74	if (weights) {
75	memcpy(weights_, weights, numberMembers_*sizeof(double));
76	} else {
77	for (int i = 0; i < numberMembers_; i++)
78	weights_[i] = i;
79	}
80	// sort so weights increasing
81	CoinSort_2(weights_, weights_ + numberMembers_, members_);
82	/*
83	Force all weights to be distinct; note that the separation enforced here
84	(1.0e-10) is not sufficien to pass the test in infeasibility().
85	*/
86
87	double last = -COIN_DBL_MAX;
88	int i;
89	for (i = 0; i < numberMembers_; i++) {
90	double possible = CoinMax(last + 1.0e-10, weights_[i]);
91	weights_[i] = possible;
92	last = possible;
93	}
94	} else {
95	members_ = NULL;
96	weights_ = NULL;
97	}
98	assert (sosType_ > 0 && sosType_ < 3);
99	}
100
101	// Copy constructor
102	CbcSOS::CbcSOS ( const CbcSOS & rhs)
103	: CbcObject(rhs)
104	{
105	shadowEstimateDown_ = rhs.shadowEstimateDown_;
106	shadowEstimateUp_ = rhs.shadowEstimateUp_;
107	downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
108	upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
109	numberTimesDown_ = rhs.numberTimesDown_;
110	numberTimesUp_ = rhs.numberTimesUp_;
111	numberMembers_ = rhs.numberMembers_;
112	sosType_ = rhs.sosType_;
113	integerValued_ = rhs.integerValued_;
114	if (numberMembers_) {
115	members_ = new int[numberMembers_];
116	weights_ = new double[numberMembers_];
117	memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
118	memcpy(weights_, rhs.weights_, numberMembers_*sizeof(double));
119	} else {
120	members_ = NULL;
121	weights_ = NULL;
122	}
123	}
124
125	// Clone
126	CbcObject *
127	CbcSOS::clone() const
128	{
129	return new CbcSOS(*this);
130	}
131
132	// Assignment operator
133	CbcSOS &
134	CbcSOS::operator=( const CbcSOS & rhs)
135	{
136	if (this != &rhs) {
137	CbcObject::operator=(rhs);
138	delete [] members_;
139	delete [] weights_;
140	shadowEstimateDown_ = rhs.shadowEstimateDown_;
141	shadowEstimateUp_ = rhs.shadowEstimateUp_;
142	downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
143	upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
144	numberTimesDown_ = rhs.numberTimesDown_;
145	numberTimesUp_ = rhs.numberTimesUp_;
146	numberMembers_ = rhs.numberMembers_;
147	sosType_ = rhs.sosType_;
148	integerValued_ = rhs.integerValued_;
149	if (numberMembers_) {
150	members_ = new int[numberMembers_];
151	weights_ = new double[numberMembers_];
152	memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
153	memcpy(weights_, rhs.weights_, numberMembers_*sizeof(double));
154	} else {
155	members_ = NULL;
156	weights_ = NULL;
157	}
158	}
159	return *this;
160	}
161
162	// Destructor
163	CbcSOS::~CbcSOS ()
164	{
165	delete [] members_;
166	delete [] weights_;
167	}
168	/*
169	Routine to calculate standard infeasibility of an SOS set and return a
170	preferred branching direction. This routine looks to have undergone
171	incomplete revision. There is vestigial code. preferredWay is unconditionally
172	set to 1. There used to be a comment `large is 0.5' but John removed it
173	at some point. Have to check to see if it no longer applies or if John
174	thought it provided too much information.
175	*/
176	double
177	CbcSOS::infeasibility(const OsiBranchingInformation * info,
178	int &preferredWay) const
179	{
180	int j;
181	int firstNonZero = -1;
182	int lastNonZero = -1;
183	OsiSolverInterface * solver = model_->solver();
184	const double * solution = model_->testSolution();
185	//const double * lower = solver->getColLower();
186	const double * upper = solver->getColUpper();
187	//double largestValue=0.0;
188	double integerTolerance =
189	model_->getDblParam(CbcModel::CbcIntegerTolerance);
190	double weight = 0.0;
191	double sum = 0.0;
192
193	// check bounds etc
194	double lastWeight = -1.0e100;
195	for (j = 0; j < numberMembers_; j++) {
196	int iColumn = members_[j];
197	/*
198	The value used here (1.0e-7) is larger than the value enforced in the
199	constructor.
200	*/
201
202	if (lastWeight >= weights_[j] - 1.0e-7)
203	throw CoinError("Weights too close together in SOS", "infeasibility", "CbcSOS");
204	double value = CoinMax(0.0, solution[iColumn]);
205	sum += value;
206	/*
207	If we're not making assumptions about integrality, why check integerTolerance
208	here? Convenient tolerance? Why not just check against the upper bound?
209
210	The calculation of weight looks to be a relic --- in the end, the value isn't
211	used to calculate either the return value or preferredWay.
212	*/
213	if (value > integerTolerance && upper[iColumn]) {
214	// Possibly due to scaling a fixed variable might slip through
215	if (value > upper[iColumn]) {
216	value = upper[iColumn];
217	// Could change to #ifdef CBC_DEBUG
218	#ifndef NDEBUG
219	if (model_->messageHandler()->logLevel() > 2)
220	printf("** Variable %d (%d) has value %g and upper bound of %g\n",
221	iColumn, j, value, upper[iColumn]);
222	#endif
223	}
224	weight += weights_[j] * value;
225	if (firstNonZero < 0)
226	firstNonZero = j;
227	lastNonZero = j;
228	}
229	}
230	/* ?? */
231	preferredWay = 1;
232	/*
233	SOS1 allows one nonzero; SOS2 allows two consecutive nonzeros. Infeasibility
234	is calculated as (.5)(range of nonzero values)/(number of members). So if
235	the first and last elements of the set are nonzero, we have maximum
236	infeasibility.
237	*/
238	if (lastNonZero - firstNonZero >= sosType_) {
239	// find where to branch
240	assert (sum > 0.0);
241	weight /= sum;
242	if (info->defaultDual_ >= 0.0 && info->usefulRegion_ && info->columnStart_) {
243	assert (sosType_ == 1);
244	int iWhere;
245	for (iWhere = firstNonZero; iWhere < lastNonZero - 1; iWhere++) {
246	if (weight < weights_[iWhere+1]) {
247	break;
248	}
249	}
250	int jColumnDown = members_[iWhere];
251	int jColumnUp = members_[iWhere+1];
252	int n = 0;
253	CoinBigIndex j;
254	double objMove = info->objective_[jColumnDown];
255	for (j = info->columnStart_[jColumnDown];
256	j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
257	double value = info->elementByColumn_[j];
258	int iRow = info->row_[j];
259	info->indexRegion_[n++] = iRow;
260	info->usefulRegion_[iRow] = value;
261	}
262	for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++) {
263	int jColumn = members_[iWhere];
264	double solValue = info->solution_[jColumn];
265	if (!solValue)
266	continue;
267	objMove -= info->objective_[jColumn] * solValue;
268	for (j = info->columnStart_[jColumn];
269	j < info->columnStart_[jColumn] + info->columnLength_[jColumn]; j++) {
270	double value = -info->elementByColumn_[j] * solValue;
271	int iRow = info->row_[j];
272	double oldValue = info->usefulRegion_[iRow];
273	if (!oldValue) {
274	info->indexRegion_[n++] = iRow;
275	} else {
276	value += oldValue;
277	if (!value)
278	value = 1.0e-100;
279	}
280	info->usefulRegion_[iRow] = value;
281	}
282	}
283	const double * pi = info->pi_;
284	const double * activity = info->rowActivity_;
285	const double * lower = info->rowLower_;
286	const double * upper = info->rowUpper_;
287	double tolerance = info->primalTolerance_;
288	double direction = info->direction_;
289	shadowEstimateDown_ = objMove * direction;
290	bool infeasible = false;
291	for (int k = 0; k < n; k++) {
292	int iRow = info->indexRegion_[k];
293	double movement = info->usefulRegion_[iRow];
294	// not this time info->usefulRegion_[iRow]=0.0;
295	if (lower[iRow] < -1.0e20)
296	assert (pi[iRow] <= 1.0e-3);
297	if (upper[iRow] > 1.0e20)
298	assert (pi[iRow] >= -1.0e-3);
299	double valueP = pi[iRow] * direction;
300	// if move makes infeasible then make at least default
301	double newValue = activity[iRow] + movement;
302	if (newValue > upper[iRow] + tolerance \|\| newValue < lower[iRow] - tolerance) {
303	shadowEstimateDown_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
304	infeasible = true;
305	}
306	}
307	if (shadowEstimateDown_ < info->integerTolerance_) {
308	if (!infeasible) {
309	shadowEstimateDown_ = 1.0e-10;
310	#ifdef COIN_DEVELOP
311	printf("zero pseudoShadowPrice\n");
312	#endif
313	} else
314	shadowEstimateDown_ = info->integerTolerance_;
315	}
316	// And other way
317	// take off
318	objMove -= info->objective_[jColumnDown];
319	for (j = info->columnStart_[jColumnDown];
320	j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
321	double value = -info->elementByColumn_[j];
322	int iRow = info->row_[j];
323	double oldValue = info->usefulRegion_[iRow];
324	if (!oldValue) {
325	info->indexRegion_[n++] = iRow;
326	} else {
327	value += oldValue;
328	if (!value)
329	value = 1.0e-100;
330	}
331	info->usefulRegion_[iRow] = value;
332	}
333	// add on
334	objMove += info->objective_[jColumnUp];
335	for (j = info->columnStart_[jColumnUp];
336	j < info->columnStart_[jColumnUp] + info->columnLength_[jColumnUp]; j++) {
337	double value = info->elementByColumn_[j];
338	int iRow = info->row_[j];
339	double oldValue = info->usefulRegion_[iRow];
340	if (!oldValue) {
341	info->indexRegion_[n++] = iRow;
342	} else {
343	value += oldValue;
344	if (!value)
345	value = 1.0e-100;
346	}
347	info->usefulRegion_[iRow] = value;
348	}
349	shadowEstimateUp_ = objMove * direction;
350	infeasible = false;
351	for (int k = 0; k < n; k++) {
352	int iRow = info->indexRegion_[k];
353	double movement = info->usefulRegion_[iRow];
354	info->usefulRegion_[iRow] = 0.0;
355	if (lower[iRow] < -1.0e20)
356	assert (pi[iRow] <= 1.0e-3);
357	if (upper[iRow] > 1.0e20)
358	assert (pi[iRow] >= -1.0e-3);
359	double valueP = pi[iRow] * direction;
360	// if move makes infeasible then make at least default
361	double newValue = activity[iRow] + movement;
362	if (newValue > upper[iRow] + tolerance \|\| newValue < lower[iRow] - tolerance) {
363	shadowEstimateUp_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
364	infeasible = true;
365	}
366	}
367	if (shadowEstimateUp_ < info->integerTolerance_) {
368	if (!infeasible) {
369	shadowEstimateUp_ = 1.0e-10;
370	#ifdef COIN_DEVELOP
371	printf("zero pseudoShadowPrice\n");
372	#endif
373	} else
374	shadowEstimateUp_ = info->integerTolerance_;
375	}
376	// adjust
377	double downCost = shadowEstimateDown_;
378	double upCost = shadowEstimateUp_;
379	if (numberTimesDown_)
380	downCost *= downDynamicPseudoRatio_ /
381	static_cast<double> (numberTimesDown_);
382	if (numberTimesUp_)
383	upCost *= upDynamicPseudoRatio_ /
384	static_cast<double> (numberTimesUp_);
385	#define WEIGHT_AFTER 0.7
386	#define WEIGHT_BEFORE 0.1
387	int stateOfSearch = model_->stateOfSearch() % 10;
388	double returnValue = 0.0;
389	double minValue = CoinMin(downCost, upCost);
390	double maxValue = CoinMax(downCost, upCost);
391	if (stateOfSearch <= 2) {
392	// no branching solution
393	returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
394	} else {
395	returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
396	}
397	#ifdef PRINT_SHADOW
398	printf("%d id - down %d %g up %d %g shadow %g, %g returned %g\n",
399	id_, numberTimesDown_, downDynamicPseudoRatio_,
400	numberTimesUp_, upDynamicPseudoRatio_, shadowEstimateDown_,
401	shadowEstimateUp_, returnValue);
402	#endif
403	return returnValue;
404	} else {
405	double value = lastNonZero - firstNonZero + 1;
406	value *= 0.5 / static_cast<double> (numberMembers_);
407	return value;
408	}
409	} else {
410	return 0.0; // satisfied
411	}
412	}
413
414	// This looks at solution and sets bounds to contain solution
415	void
416	CbcSOS::feasibleRegion()
417	{
418	int j;
419	int firstNonZero = -1;
420	int lastNonZero = -1;
421	OsiSolverInterface * solver = model_->solver();
422	const double * solution = model_->testSolution();
423	//const double * lower = solver->getColLower();
424	const double * upper = solver->getColUpper();
425	double integerTolerance =
426	model_->getDblParam(CbcModel::CbcIntegerTolerance);
427	double weight = 0.0;
428	double sum = 0.0;
429
430	for (j = 0; j < numberMembers_; j++) {
431	int iColumn = members_[j];
432	double value = CoinMax(0.0, solution[iColumn]);
433	sum += value;
434	if (value > integerTolerance && upper[iColumn]) {
435	weight += weights_[j] * value;
436	if (firstNonZero < 0)
437	firstNonZero = j;
438	lastNonZero = j;
439	}
440	}
441	assert (lastNonZero - firstNonZero < sosType_) ;
442	for (j = 0; j < firstNonZero; j++) {
443	int iColumn = members_[j];
444	solver->setColUpper(iColumn, 0.0);
445	}
446	for (j = lastNonZero + 1; j < numberMembers_; j++) {
447	int iColumn = members_[j];
448	solver->setColUpper(iColumn, 0.0);
449	}
450	}
451	// Redoes data when sequence numbers change
452	void
453	CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
454	{
455	model_ = model;
456	int n2 = 0;
457	for (int j = 0; j < numberMembers_; j++) {
458	int iColumn = members_[j];
459	int i;
460	for (i = 0; i < numberColumns; i++) {
461	if (originalColumns[i] == iColumn)
462	break;
463	}
464	if (i < numberColumns) {
465	members_[n2] = i;
466	weights_[n2++] = weights_[j];
467	}
468	}
469	if (n2 < numberMembers_) {
470	//printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
471	numberMembers_ = n2;
472	}
473	}
474	CbcBranchingObject *
475	CbcSOS::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /info/, int way)
476	{
477	int j;
478	const double * solution = model_->testSolution();
479	double integerTolerance =
480	model_->getDblParam(CbcModel::CbcIntegerTolerance);
481	//OsiSolverInterface * solver = model_->solver();
482	const double * upper = solver->getColUpper();
483	int firstNonFixed = -1;
484	int lastNonFixed = -1;
485	int firstNonZero = -1;
486	int lastNonZero = -1;
487	double weight = 0.0;
488	double sum = 0.0;
489	for (j = 0; j < numberMembers_; j++) {
490	int iColumn = members_[j];
491	if (upper[iColumn]) {
492	double value = CoinMax(0.0, solution[iColumn]);
493	sum += value;
494	if (firstNonFixed < 0)
495	firstNonFixed = j;
496	lastNonFixed = j;
497	if (value > integerTolerance) {
498	weight += weights_[j] * value;
499	if (firstNonZero < 0)
500	firstNonZero = j;
501	lastNonZero = j;
502	}
503	}
504	}
505	assert (lastNonZero - firstNonZero >= sosType_) ;
506	// find where to branch
507	assert (sum > 0.0);
508	weight /= sum;
509	int iWhere;
510	double separator = 0.0;
511	for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
512	if (weight < weights_[iWhere+1])
513	break;
514	if (sosType_ == 1) {
515	// SOS 1
516	separator = 0.5 * (weights_[iWhere] + weights_[iWhere+1]);
517	} else {
518	// SOS 2
519	if (iWhere == firstNonFixed)
520	iWhere++;;
521	if (iWhere == lastNonFixed - 1)
522	iWhere = lastNonFixed - 2;
523	separator = weights_[iWhere+1];
524	}
525	// create object
526	CbcBranchingObject * branch;
527	branch = new CbcSOSBranchingObject(model_, this, way, separator);
528	branch->setOriginalObject(this);
529	return branch;
530	}
531	/* Pass in information on branch just done and create CbcObjectUpdateData instance.
532	If object does not need data then backward pointer will be NULL.
533	Assumes can get information from solver */
534	CbcObjectUpdateData
535	CbcSOS::createUpdateInformation(const OsiSolverInterface * solver,
536	const CbcNode * node,
537	const CbcBranchingObject * branchingObject)
538	{
539	double originalValue = node->objectiveValue();
540	int originalUnsatisfied = node->numberUnsatisfied();
541	double objectiveValue = solver->getObjValue() * solver->getObjSense();
542	int unsatisfied = 0;
543	int i;
544	//might be base model - doesn't matter
545	int numberIntegers = model_->numberIntegers();;
546	const double * solution = solver->getColSolution();
547	//const double * lower = solver->getColLower();
548	//const double * upper = solver->getColUpper();
549	double change = CoinMax(0.0, objectiveValue - originalValue);
550	int iStatus;
551	if (solver->isProvenOptimal())
552	iStatus = 0; // optimal
553	else if (solver->isIterationLimitReached()
554	&& !solver->isDualObjectiveLimitReached())
555	iStatus = 2; // unknown
556	else
557	iStatus = 1; // infeasible
558
559	bool feasible = iStatus != 1;
560	if (feasible) {
561	double integerTolerance =
562	model_->getDblParam(CbcModel::CbcIntegerTolerance);
563	const int * integerVariable = model_->integerVariable();
564	for (i = 0; i < numberIntegers; i++) {
565	int j = integerVariable[i];
566	double value = solution[j];
567	double nearest = floor(value + 0.5);
568	if (fabs(value - nearest) > integerTolerance)
569	unsatisfied++;
570	}
571	}
572	int way = branchingObject->way();
573	way = - way; // because after branch so moved on
574	double value = branchingObject->value();
575	CbcObjectUpdateData newData (this, way,
576	change, iStatus,
577	originalUnsatisfied - unsatisfied, value);
578	newData.originalObjective_ = originalValue;
579	// Solvers know about direction
580	double direction = solver->getObjSense();
581	solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
582	newData.cutoff_ *= direction;
583	return newData;
584	}
585	// Update object by CbcObjectUpdateData
586	void
587	CbcSOS::updateInformation(const CbcObjectUpdateData & data)
588	{
589	bool feasible = data.status_ != 1;
590	int way = data.way_;
591	//double value = data.branchingValue_;
592	double originalValue = data.originalObjective_;
593	double change = data.change_;
594	if (way < 0) {
595	// down
596	if (!feasible) {
597	double distanceToCutoff = 0.0;
598	//double objectiveValue = model_->getCurrentMinimizationObjValue();
599	distanceToCutoff = model_->getCutoff() - originalValue;
600	if (distanceToCutoff < 1.0e20)
601	change = distanceToCutoff * 2.0;
602	else
603	change = (downDynamicPseudoRatio_ * shadowEstimateDown_ + 1.0e-3) * 10.0;
604	}
605	change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
606	#ifdef PRINT_SHADOW
607	if (numberTimesDown_)
608	printf("Updating id %d - down change %g (true %g) - ndown %d estimated change %g - raw shadow estimate %g\n",
609	id_, change, data.change_, numberTimesDown_, shadowEstimateDown_*
610	(downDynamicPseudoRatio_ / ((double) numberTimesDown_)),
611	shadowEstimateDown_);
612	else
613	printf("Updating id %d - down change %g (true %g) - shadow estimate %g\n",
614	id_, change, data.change_, shadowEstimateDown_);
615	#endif
616	numberTimesDown_++;
617	downDynamicPseudoRatio_ += change / shadowEstimateDown_;
618	} else {
619	// up
620	if (!feasible) {
621	double distanceToCutoff = 0.0;
622	//double objectiveValue = model_->getCurrentMinimizationObjValue();
623	distanceToCutoff = model_->getCutoff() - originalValue;
624	if (distanceToCutoff < 1.0e20)
625	change = distanceToCutoff * 2.0;
626	else
627	change = (upDynamicPseudoRatio_ * shadowEstimateUp_ + 1.0e-3) * 10.0;
628	}
629	change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
630	#ifdef PRINT_SHADOW
631	if (numberTimesUp_)
632	printf("Updating id %d - up change %g (true %g) - nup %d estimated change %g - raw shadow estimate %g\n",
633	id_, change, data.change_, numberTimesUp_, shadowEstimateUp_*
634	(upDynamicPseudoRatio_ / ((double) numberTimesUp_)),
635	shadowEstimateUp_);
636	else
637	printf("Updating id %d - up change %g (true %g) - shadow estimate %g\n",
638	id_, change, data.change_, shadowEstimateUp_);
639	#endif
640	numberTimesUp_++;
641	upDynamicPseudoRatio_ += change / shadowEstimateUp_;
642	}
643	}
644
645	/* Create an OsiSolverBranch object
646
647	This returns NULL if branch not represented by bound changes
648	*/
649	OsiSolverBranch *
650	CbcSOS::solverBranch() const
651	{
652	int j;
653	const double * solution = model_->testSolution();
654	double integerTolerance =
655	model_->getDblParam(CbcModel::CbcIntegerTolerance);
656	OsiSolverInterface * solver = model_->solver();
657	const double * upper = solver->getColUpper();
658	int firstNonFixed = -1;
659	int lastNonFixed = -1;
660	int firstNonZero = -1;
661	int lastNonZero = -1;
662	double weight = 0.0;
663	double sum = 0.0;
664	double * fix = new double[numberMembers_];
665	int * which = new int[numberMembers_];
666	for (j = 0; j < numberMembers_; j++) {
667	int iColumn = members_[j];
668	// fix all on one side or other (even if fixed)
669	fix[j] = 0.0;
670	which[j] = iColumn;
671	if (upper[iColumn]) {
672	double value = CoinMax(0.0, solution[iColumn]);
673	sum += value;
674	if (firstNonFixed < 0)
675	firstNonFixed = j;
676	lastNonFixed = j;
677	if (value > integerTolerance) {
678	weight += weights_[j] * value;
679	if (firstNonZero < 0)
680	firstNonZero = j;
681	lastNonZero = j;
682	}
683	}
684	}
685	assert (lastNonZero - firstNonZero >= sosType_) ;
686	// find where to branch
687	assert (sum > 0.0);
688	weight /= sum;
689	// down branch fixes ones above weight to 0
690	int iWhere;
691	int iDownStart = 0;
692	int iUpEnd = 0;
693	for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
694	if (weight < weights_[iWhere+1])
695	break;
696	if (sosType_ == 1) {
697	// SOS 1
698	iUpEnd = iWhere + 1;
699	iDownStart = iUpEnd;
700	} else {
701	// SOS 2
702	if (iWhere == firstNonFixed)
703	iWhere++;;
704	if (iWhere == lastNonFixed - 1)
705	iWhere = lastNonFixed - 2;
706	iUpEnd = iWhere + 1;
707	iDownStart = iUpEnd + 1;
708	}
709	//
710	OsiSolverBranch * branch = new OsiSolverBranch();
711	branch->addBranch(-1, 0, NULL, NULL, numberMembers_ - iDownStart, which + iDownStart, fix);
712	branch->addBranch(1, 0, NULL, NULL, iUpEnd, which, fix);
713	delete [] fix;
714	delete [] which;
715	return branch;
716	}
717	// Construct an OsiSOS object
718	OsiSOS *
719	CbcSOS::osiObject(const OsiSolverInterface * solver) const
720	{
721	OsiSOS * obj = new OsiSOS(solver, numberMembers_, members_, weights_, sosType_);
722	obj->setPriority(priority());
723	return obj;
724	}
725
726	// Default Constructor
727	CbcSOSBranchingObject::CbcSOSBranchingObject()
728	: CbcBranchingObject(),
729	firstNonzero_(-1),
730	lastNonzero_(-1)
731	{
732	set_ = NULL;
733	separator_ = 0.0;
734	}
735
736	// Useful constructor
737	CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
738	const CbcSOS * set,
739	int way ,
740	double separator)
741	: CbcBranchingObject(model, set->id(), way, 0.5)
742	{
743	set_ = set;
744	separator_ = separator;
745	computeNonzeroRange();
746	}
747
748	// Copy constructor
749	CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
750	: CbcBranchingObject(rhs),
751	firstNonzero_(rhs.firstNonzero_),
752	lastNonzero_(rhs.lastNonzero_)
753	{
754	set_ = rhs.set_;
755	separator_ = rhs.separator_;
756	}
757
758	// Assignment operator
759	CbcSOSBranchingObject &
760	CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject & rhs)
761	{
762	if (this != &rhs) {
763	CbcBranchingObject::operator=(rhs);
764	set_ = rhs.set_;
765	separator_ = rhs.separator_;
766	firstNonzero_ = rhs.firstNonzero_;
767	lastNonzero_ = rhs.lastNonzero_;
768	}
769	return *this;
770	}
771	CbcBranchingObject *
772	CbcSOSBranchingObject::clone() const
773	{
774	return (new CbcSOSBranchingObject(*this));
775	}
776
777
778	// Destructor
779	CbcSOSBranchingObject::~CbcSOSBranchingObject ()
780	{
781	}
782
783	void
784	CbcSOSBranchingObject::computeNonzeroRange()
785	{
786	const int numberMembers = set_->numberMembers();
787	const double * weights = set_->weights();
788	int i = 0;
789	if (way_ < 0) {
790	for ( i = 0; i < numberMembers; i++) {
791	if (weights[i] > separator_)
792	break;
793	}
794	assert (i < numberMembers);
795	firstNonzero_ = 0;
796	lastNonzero_ = i;
797	} else {
798	for ( i = 0; i < numberMembers; i++) {
799	if (weights[i] >= separator_)
800	break;
801	}
802	assert (i < numberMembers);
803	firstNonzero_ = i;
804	lastNonzero_ = numberMembers;
805	}
806	}
807
808	double
809	CbcSOSBranchingObject::branch()
810	{
811	decrementNumberBranchesLeft();
812	int numberMembers = set_->numberMembers();
813	const int * which = set_->members();
814	const double * weights = set_->weights();
815	OsiSolverInterface * solver = model_->solver();
816	//const double * lower = solver->getColLower();
817	//const double * upper = solver->getColUpper();
818	// *** for way - up means fix all those in down section
819	if (way_ < 0) {
820	int i;
821	for ( i = 0; i < numberMembers; i++) {
822	if (weights[i] > separator_)
823	break;
824	}
825	assert (i < numberMembers);
826	for (; i < numberMembers; i++)
827	solver->setColUpper(which[i], 0.0);
828	way_ = 1; // Swap direction
829	} else {
830	int i;
831	for ( i = 0; i < numberMembers; i++) {
832	if (weights[i] >= separator_)
833	break;
834	else
835	solver->setColUpper(which[i], 0.0);
836	}
837	assert (i < numberMembers);
838	way_ = -1; // Swap direction
839	}
840	computeNonzeroRange();
841	return 0.0;
842	}
843	/* Update bounds in solver as in 'branch' and update given bounds.
844	branchState is -1 for 'down' +1 for 'up' */
845	void
846	CbcSOSBranchingObject::fix(OsiSolverInterface * solver,
847	double * /lower/, double * upper,
848	int branchState) const
849	{
850	int numberMembers = set_->numberMembers();
851	const int * which = set_->members();
852	const double * weights = set_->weights();
853	//const double * lower = solver->getColLower();
854	//const double * upper = solver->getColUpper();
855	// *** for way - up means fix all those in down section
856	if (branchState < 0) {
857	int i;
858	for ( i = 0; i < numberMembers; i++) {
859	if (weights[i] > separator_)
860	break;
861	}
862	assert (i < numberMembers);
863	for (; i < numberMembers; i++) {
864	solver->setColUpper(which[i], 0.0);
865	upper[which[i]] = 0.0;
866	}
867	} else {
868	int i;
869	for ( i = 0; i < numberMembers; i++) {
870	if (weights[i] >= separator_) {
871	break;
872	} else {
873	solver->setColUpper(which[i], 0.0);
874	upper[which[i]] = 0.0;
875	}
876	}
877	assert (i < numberMembers);
878	}
879	}
880	// Print what would happen
881	void
882	CbcSOSBranchingObject::print()
883	{
884	int numberMembers = set_->numberMembers();
885	const int * which = set_->members();
886	const double * weights = set_->weights();
887	OsiSolverInterface * solver = model_->solver();
888	//const double * lower = solver->getColLower();
889	const double * upper = solver->getColUpper();
890	int first = numberMembers;
891	int last = -1;
892	int numberFixed = 0;
893	int numberOther = 0;
894	int i;
895	for ( i = 0; i < numberMembers; i++) {
896	double bound = upper[which[i]];
897	if (bound) {
898	first = CoinMin(first, i);
899	last = CoinMax(last, i);
900	}
901	}
902	// *** for way - up means fix all those in down section
903	if (way_ < 0) {
904	printf("SOS Down");
905	for ( i = 0; i < numberMembers; i++) {
906	double bound = upper[which[i]];
907	if (weights[i] > separator_)
908	break;
909	else if (bound)
910	numberOther++;
911	}
912	assert (i < numberMembers);
913	for (; i < numberMembers; i++) {
914	double bound = upper[which[i]];
915	if (bound)
916	numberFixed++;
917	}
918	} else {
919	printf("SOS Up");
920	for ( i = 0; i < numberMembers; i++) {
921	double bound = upper[which[i]];
922	if (weights[i] >= separator_)
923	break;
924	else if (bound)
925	numberFixed++;
926	}
927	assert (i < numberMembers);
928	for (; i < numberMembers; i++) {
929	double bound = upper[which[i]];
930	if (bound)
931	numberOther++;
932	}
933	}
934	printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
935	separator_, which[first], weights[first], which[last], weights[last], numberFixed, numberOther);
936	}
937
938	/** Compare the original object of \c this with the original object of \c
939	brObj. Assumes that there is an ordering of the original objects.
940	This method should be invoked only if \c this and brObj are of the same
941	type.
942	Return negative/0/positive depending on whether \c this is
943	smaller/same/larger than the argument.
944	*/
945	int
946	CbcSOSBranchingObject::compareOriginalObject
947	(const CbcBranchingObject* brObj) const
948	{
949	const CbcSOSBranchingObject* br =
950	dynamic_cast<const CbcSOSBranchingObject*>(brObj);
951	assert(br);
952	const CbcSOS* s0 = set_;
953	const CbcSOS* s1 = br->set_;
954	if (s0->sosType() != s1->sosType()) {
955	return s0->sosType() - s1->sosType();
956	}
957	if (s0->numberMembers() != s1->numberMembers()) {
958	return s0->numberMembers() - s1->numberMembers();
959	}
960	const int memberCmp = memcmp(s0->members(), s1->members(),
961	s0->numberMembers() * sizeof(int));
962	if (memberCmp != 0) {
963	return memberCmp;
964	}
965	return memcmp(s0->weights(), s1->weights(),
966	s0->numberMembers() * sizeof(double));
967	}
968
969	/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
970	same type and must have the same original object, but they may have
971	different feasible regions.
972	Return the appropriate CbcRangeCompare value (first argument being the
973	sub/superset if that's the case). In case of overlap (and if \c
974	replaceIfOverlap is true) replace the current branching object with one
975	whose feasible region is the overlap.
976	*/
977	CbcRangeCompare
978	CbcSOSBranchingObject::compareBranchingObject
979	(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
980	{
981	const CbcSOSBranchingObject* br =
982	dynamic_cast<const CbcSOSBranchingObject*>(brObj);
983	assert(br);
984	if (firstNonzero_ < br->firstNonzero_) {
985	if (lastNonzero_ >= br->lastNonzero_) {
986	return CbcRangeSuperset;
987	} else if (lastNonzero_ <= br->firstNonzero_) {
988	return CbcRangeDisjoint;
989	} else {
990	// overlap
991	if (replaceIfOverlap) {
992	firstNonzero_ = br->firstNonzero_;
993	}
994	return CbcRangeOverlap;
995	}
996	} else if (firstNonzero_ > br->firstNonzero_) {
997	if (lastNonzero_ <= br->lastNonzero_) {
998	return CbcRangeSubset;
999	} else if (firstNonzero_ >= br->lastNonzero_) {
1000	return CbcRangeDisjoint;
1001	} else {
1002	// overlap
1003	if (replaceIfOverlap) {
1004	lastNonzero_ = br->lastNonzero_;
1005	}
1006	return CbcRangeOverlap;
1007	}
1008	} else {
1009	if (lastNonzero_ == br->lastNonzero_) {
1010	return CbcRangeSame;
1011	}
1012	return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
1013	}
1014	return CbcRangeSame; // fake return
1015	}
1016

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

Download in other formats: