Context Navigation

source: stable/2.9/Cbc/src/CbcFathomDynamicProgramming.cpp @ 2127

Last change on this file since 2127 was 1886, checked in by stefan, 7 years ago
fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 34.2 KB

Line
1	/*
2	$Id: CbcFathomDynamicProgramming.cpp 1886 2013-04-06 18:13:35Z forrest $
3	*/
4	// Copyright (C) 2004, International Business Machines
5	// Corporation and others. All Rights Reserved.
6	// This code is licensed under the terms of the Eclipse Public License (EPL).
7
8	#if defined(_MSC_VER)
9	// Turn off compiler warning about long names
10	# pragma warning(disable:4786)
11	#endif
12	#include <cassert>
13	#include <cstdlib>
14	#include <cmath>
15	#include <cfloat>
16
17	#include "OsiSolverInterface.hpp"
18	#include "CbcModel.hpp"
19	#include "CbcMessage.hpp"
20	#include "CbcFathomDynamicProgramming.hpp"
21	#include "CoinHelperFunctions.hpp"
22	#include "CoinPackedMatrix.hpp"
23	#include "CoinSort.hpp"
24	// Default Constructor
25	CbcFathomDynamicProgramming::CbcFathomDynamicProgramming()
26	: CbcFathom(),
27	size_(0),
28	type_(-1),
29	cost_(NULL),
30	back_(NULL),
31	lookup_(NULL),
32	indices_(NULL),
33	numberActive_(0),
34	maximumSizeAllowed_(1000000),
35	startBit_(NULL),
36	numberBits_(NULL),
37	rhs_(NULL),
38	coefficients_(NULL),
39	target_(0),
40	numberNonOne_(0),
41	bitPattern_(0),
42	algorithm_(-1)
43	{
44
45	}
46
47	// Constructor from model
48	CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(CbcModel & model)
49	: CbcFathom(model),
50	cost_(NULL),
51	back_(NULL),
52	lookup_(NULL),
53	indices_(NULL),
54	numberActive_(0),
55	maximumSizeAllowed_(1000000),
56	startBit_(NULL),
57	numberBits_(NULL),
58	rhs_(NULL),
59	coefficients_(NULL),
60	target_(0),
61	numberNonOne_(0),
62	bitPattern_(0),
63	algorithm_(-1)
64	{
65	type_ = checkPossible();
66	}
67
68	// Destructor
69	CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming ()
70	{
71	gutsOfDelete();
72	}
73	// Does deleteions
74	void
75	CbcFathomDynamicProgramming::gutsOfDelete()
76	{
77	delete [] cost_;
78	delete [] back_;
79	delete [] lookup_;
80	delete [] indices_;
81	delete [] startBit_;
82	delete [] numberBits_;
83	delete [] rhs_;
84	delete [] coefficients_;
85	cost_ = NULL;
86	back_ = NULL;
87	lookup_ = NULL;
88	indices_ = NULL;
89	startBit_ = NULL;
90	numberBits_ = NULL;
91	rhs_ = NULL;
92	coefficients_ = NULL;
93	}
94	// Clone
95	CbcFathom *
96	CbcFathomDynamicProgramming::clone() const
97	{
98	return new CbcFathomDynamicProgramming(*this);
99	}
100
101	// Copy constructor
102	CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(const CbcFathomDynamicProgramming & rhs)
103	:
104	CbcFathom(rhs),
105	size_(rhs.size_),
106	type_(rhs.type_),
107	cost_(NULL),
108	back_(NULL),
109	lookup_(NULL),
110	indices_(NULL),
111	numberActive_(rhs.numberActive_),
112	maximumSizeAllowed_(rhs.maximumSizeAllowed_),
113	startBit_(NULL),
114	numberBits_(NULL),
115	rhs_(NULL),
116	coefficients_(NULL),
117	target_(rhs.target_),
118	numberNonOne_(rhs.numberNonOne_),
119	bitPattern_(rhs.bitPattern_),
120	algorithm_(rhs.algorithm_)
121	{
122	if (size_) {
123	cost_ = CoinCopyOfArray(rhs.cost_, size_);
124	back_ = CoinCopyOfArray(rhs.back_, size_);
125	int numberRows = model_->getNumRows();
126	lookup_ = CoinCopyOfArray(rhs.lookup_, numberRows);
127	startBit_ = CoinCopyOfArray(rhs.startBit_, numberActive_);
128	indices_ = CoinCopyOfArray(rhs.indices_, numberActive_);
129	numberBits_ = CoinCopyOfArray(rhs.numberBits_, numberActive_);
130	rhs_ = CoinCopyOfArray(rhs.rhs_, numberActive_);
131	coefficients_ = CoinCopyOfArray(rhs.coefficients_, numberActive_);
132	}
133	}
134	// Returns type
135	int
136	CbcFathomDynamicProgramming::checkPossible(int allowableSize)
137	{
138	algorithm_ = -1;
139	assert(model_->solver());
140	OsiSolverInterface * solver = model_->solver();
141	const CoinPackedMatrix * matrix = solver->getMatrixByCol();
142
143	int numberIntegers = model_->numberIntegers();
144	int numberColumns = solver->getNumCols();
145	size_ = 0;
146	if (numberIntegers != numberColumns)
147	return -1; // can't do dynamic programming
148
149	const double * lower = solver->getColLower();
150	const double * upper = solver->getColUpper();
151	const double * rowUpper = solver->getRowUpper();
152
153	int numberRows = model_->getNumRows();
154	int i;
155
156	// First check columns to see if possible
157	double * rhs = new double [numberRows];
158	CoinCopyN(rowUpper, numberRows, rhs);
159
160	// Column copy
161	const double * element = matrix->getElements();
162	const int * row = matrix->getIndices();
163	const CoinBigIndex * columnStart = matrix->getVectorStarts();
164	const int * columnLength = matrix->getVectorLengths();
165	bool bad = false;
166	/* It is just possible that we could say okay as
167	variables may get fixed but seems unlikely */
168	for (i = 0; i < numberColumns; i++) {
169	int j;
170	double lowerValue = lower[i];
171	assert (lowerValue == floor(lowerValue));
172	for (j = columnStart[i];
173	j < columnStart[i] + columnLength[i]; j++) {
174	int iRow = row[j];
175	double value = element[j];
176	if (upper[i] > lowerValue && (value <= 0.0 \|\| value != floor(value)))
177	bad = true;
178	if (lowerValue)
179	rhs[iRow] -= lowerValue * value;
180	}
181	}
182	// check possible (at present do not allow covering)
183	int numberActive = 0;
184	bool infeasible = false;
185	bool saveBad = bad;
186	for (i = 0; i < numberRows; i++) {
187	if (rhs[i] < 0)
188	infeasible = true;
189	else if (rhs[i] > 1.0e5 \|\| fabs(rhs[i] - floor(rhs[i] + 0.5)) > 1.0e-7)
190	bad = true;
191	else if (rhs[i] > 0.0)
192	numberActive++;
193	}
194	if (bad \|\| infeasible) {
195	delete [] rhs;
196	if (!saveBad && infeasible)
197	return -2;
198	else
199	return -1;
200	}
201	// check size of array needed
202	double size = 1.0;
203	double check = COIN_INT_MAX;
204	for (i = 0; i < numberRows; i++) {
205	int n = static_cast<int> (floor(rhs[i] + 0.5));
206	if (n) {
207	n++; // allow for 0,1... n
208	if (numberActive != 1) {
209	// power of 2
210	int iBit = 0;
211	int k = n;
212	k &= ~1;
213	while (k) {
214	iBit++;
215	k &= ~(1 << iBit);
216	}
217	// See if exact power
218	if (n != (1 << iBit)) {
219	// round up to next power of 2
220	n = 1 << (iBit + 1);
221	}
222	size *= n;
223	if (size >= check)
224	break;
225	} else {
226	size = n; // just one constraint
227	}
228	}
229	}
230	// set size needed
231	if (size >= check)
232	size_ = COIN_INT_MAX;
233	else
234	size_ = static_cast<int> (size);
235
236	int n01 = 0;
237	int nbadcoeff = 0;
238	// See if we can tighten bounds
239	for (i = 0; i < numberColumns; i++) {
240	int j;
241	double lowerValue = lower[i];
242	double gap = upper[i] - lowerValue;
243	for (j = columnStart[i];
244	j < columnStart[i] + columnLength[i]; j++) {
245	int iRow = row[j];
246	double value = element[j];
247	if (value != 1.0)
248	nbadcoeff++;
249	if (gap*value > rhs[iRow] + 1.0e-8)
250	gap = rhs[iRow] / value;
251	}
252	gap = lowerValue + floor(gap + 1.0e-7);
253	if (gap < upper[i])
254	solver->setColUpper(i, gap);
255	if (gap <= 1.0)
256	n01++;
257	}
258	if (allowableSize && size_ <= allowableSize) {
259	if (n01 == numberColumns && !nbadcoeff)
260	algorithm_ = 0; // easiest
261	else
262	algorithm_ = 1;
263	}
264	if (allowableSize && size_ <= allowableSize) {
265	numberActive_ = numberActive;
266	indices_ = new int [numberActive_];
267	cost_ = new double [size_];
268	CoinFillN(cost_, size_, COIN_DBL_MAX);
269	// but do nothing is okay
270	cost_[0] = 0.0;
271	back_ = new int[size_];
272	CoinFillN(back_, size_, -1);
273	startBit_ = new int[numberActive_];
274	numberBits_ = new int[numberActive_];
275	lookup_ = new int [numberRows];
276	rhs_ = new int [numberActive_];
277	numberActive = 0;
278	int kBit = 0;
279	for (i = 0; i < numberRows; i++) {
280	int n = static_cast<int> (floor(rhs[i] + 0.5));
281	if (n) {
282	lookup_[i] = numberActive;
283	rhs_[numberActive] = n;
284	startBit_[numberActive] = kBit;
285	n++; // allow for 0,1... n
286	int iBit = 0;
287	// power of 2
288	int k = n;
289	k &= ~1;
290	while (k) {
291	iBit++;
292	k &= ~(1 << iBit);
293	}
294	// See if exact power
295	if (n != (1 << iBit)) {
296	// round up to next power of 2
297	iBit++;
298	}
299	if (numberActive != 1) {
300	n = 1 << iBit;
301	size *= n;
302	if (size >= check)
303	break;
304	} else {
305	size = n; // just one constraint
306	}
307	numberBits_[numberActive++] = iBit;
308	kBit += iBit;
309	} else {
310	lookup_[i] = -1;
311	}
312	}
313	const double * rowLower = solver->getRowLower();
314	if (algorithm_ == 0) {
315	// rhs 1 and coefficients 1
316	// Get first possible solution for printing
317	target_ = -1;
318	int needed = 0;
319	int numberActive = 0;
320	for (i = 0; i < numberRows; i++) {
321	int newRow = lookup_[i];
322	if (newRow >= 0) {
323	if (rowLower[i] == rowUpper[i]) {
324	needed += 1 << numberActive;
325	numberActive++;
326	}
327	}
328	}
329	for (i = 0; i < size_; i++) {
330	if ((i&needed) == needed) {
331	break;
332	}
333	}
334	target_ = i;
335	} else {
336	coefficients_ = new int[numberActive_];
337	// If not too many general rhs then we can be more efficient
338	numberNonOne_ = 0;
339	for (i = 0; i < numberActive_; i++) {
340	if (rhs_[i] != 1)
341	numberNonOne_++;
342	}
343	if (numberNonOne_*2 < numberActive_) {
344	// put rhs >1 every second
345	int * permute = new int[numberActive_];
346	int * temp = new int[numberActive_];
347	// try different ways
348	int k = 0;
349	for (i = 0; i < numberRows; i++) {
350	int newRow = lookup_[i];
351	if (newRow >= 0 && rhs_[newRow] > 1) {
352	permute[newRow] = k;
353	k += 2;
354	}
355	}
356	// adjust so k points to last
357	k -= 2;
358	// and now rest
359	int k1 = 1;
360	for (i = 0; i < numberRows; i++) {
361	int newRow = lookup_[i];
362	if (newRow >= 0 && rhs_[newRow] == 1) {
363	permute[newRow] = k1;
364	k1++;
365	if (k1 <= k)
366	k1++;
367	}
368	}
369	for (i = 0; i < numberActive_; i++) {
370	int put = permute[i];
371	temp[put] = rhs_[i];
372	}
373	memcpy(rhs_, temp, numberActive_*sizeof(int));
374	for (i = 0; i < numberActive_; i++) {
375	int put = permute[i];
376	temp[put] = numberBits_[i];
377	}
378	memcpy(numberBits_, temp, numberActive_*sizeof(int));
379	k = 0;
380	for (i = 0; i < numberActive_; i++) {
381	startBit_[i] = k;
382	k += numberBits_[i];
383	}
384	for (i = 0; i < numberRows; i++) {
385	int newRow = lookup_[i];
386	if (newRow >= 0)
387	lookup_[i] = permute[newRow];
388	}
389	delete [] permute;
390	delete [] temp;
391	// mark new method
392	algorithm_ = 2;
393	}
394	// Get first possible solution for printing
395	target_ = -1;
396	int needed = 0;
397	int * lower2 = new int[numberActive_];
398	for (i = 0; i < numberRows; i++) {
399	int newRow = lookup_[i];
400	if (newRow >= 0) {
401	int gap = static_cast<int> (rowUpper[i] - CoinMax(0.0, rowLower[i]));
402	lower2[newRow] = rhs_[newRow] - gap;
403	int numberBits = numberBits_[newRow];
404	int startBit = startBit_[newRow];
405	if (numberBits == 1 && !gap) {
406	needed \|= 1 << startBit;
407	}
408	}
409	}
410	for (i = 0; i < size_; i++) {
411	if ((i&needed) == needed) {
412	// this one may do
413	bool good = true;
414	for (int kk = 0; kk < numberActive_; kk++) {
415	int numberBits = numberBits_[kk];
416	int startBit = startBit_[kk];
417	int size = 1 << numberBits;
418	int start = 1 << startBit;
419	int mask = start * (size - 1);
420	int level = (i & mask) >> startBit;
421	if (level < lower2[kk]) {
422	good = false;
423	break;
424	}
425	}
426	if (good) {
427	break;
428	}
429	}
430	}
431	delete [] lower2;
432	target_ = i;
433	}
434	}
435	delete [] rhs;
436	if (allowableSize && size_ > allowableSize) {
437	COIN_DETAIL_PRINT(printf("Too large - need %d entries x 8 bytes\n", size_));
438	return -1; // too big
439	} else {
440	return algorithm_;
441	}
442	}
443
444	// Resets stuff if model changes
445	void
446	CbcFathomDynamicProgramming::resetModel(CbcModel * model)
447	{
448	model_ = model;
449	type_ = checkPossible();
450	}
451	int
452	CbcFathomDynamicProgramming::fathom(double * & betterSolution)
453	{
454	int returnCode = 0;
455	int type = checkPossible(maximumSizeAllowed_);
456	assert (type != -1);
457	if (type == -2) {
458	// infeasible (so complete search done)
459	return 1;
460	}
461	if (algorithm_ >= 0) {
462	OsiSolverInterface * solver = model_->solver();
463	const double * lower = solver->getColLower();
464	const double * upper = solver->getColUpper();
465	const double * objective = solver->getObjCoefficients();
466	double direction = solver->getObjSense();
467	const CoinPackedMatrix * matrix = solver->getMatrixByCol();
468	// Column copy
469	const double * element = matrix->getElements();
470	const int * row = matrix->getIndices();
471	const CoinBigIndex * columnStart = matrix->getVectorStarts();
472	const int * columnLength = matrix->getVectorLengths();
473	const double * rowLower = solver->getRowLower();
474	const double * rowUpper = solver->getRowUpper();
475	int numberRows = model_->getNumRows();
476
477	int numberColumns = solver->getNumCols();
478	double offset;
479	solver->getDblParam(OsiObjOffset, offset);
480	double fixedObj = -offset;
481	int i;
482	// may be possible
483	double bestAtTarget = COIN_DBL_MAX;
484	for (i = 0; i < numberColumns; i++) {
485	if (size_ > 10000000 && (i % 100) == 0)
486	COIN_DETAIL_PRINT(printf("column %d\n", i));
487	double lowerValue = lower[i];
488	assert (lowerValue == floor(lowerValue));
489	double cost = direction * objective[i];
490	fixedObj += lowerValue * cost;
491	int gap = static_cast<int> (upper[i] - lowerValue);
492	CoinBigIndex start = columnStart[i];
493	tryColumn(columnLength[i], row + start, element + start, cost, gap);
494	if (cost_[target_] < bestAtTarget) {
495	if (model_->messageHandler()->logLevel() > 1)
496	printf("At column %d new best objective of %g\n", i, cost_[target_]);
497	bestAtTarget = cost_[target_];
498	}
499	}
500	returnCode = 1;
501	int needed = 0;
502	double bestValue = COIN_DBL_MAX;
503	int iBest = -1;
504	if (algorithm_ == 0) {
505	int numberActive = 0;
506	for (i = 0; i < numberRows; i++) {
507	int newRow = lookup_[i];
508	if (newRow >= 0) {
509	if (rowLower[i] == rowUpper[i]) {
510	needed += 1 << numberActive;
511	numberActive++;
512	}
513	}
514	}
515	for (i = 0; i < size_; i++) {
516	if ((i&needed) == needed) {
517	// this one will do
518	if (cost_[i] < bestValue) {
519	bestValue = cost_[i];
520	iBest = i;
521	}
522	}
523	}
524	} else {
525	int * lower = new int[numberActive_];
526	for (i = 0; i < numberRows; i++) {
527	int newRow = lookup_[i];
528	if (newRow >= 0) {
529	int gap = static_cast<int> (rowUpper[i] - CoinMax(0.0, rowLower[i]));
530	lower[newRow] = rhs_[newRow] - gap;
531	int numberBits = numberBits_[newRow];
532	int startBit = startBit_[newRow];
533	if (numberBits == 1 && !gap) {
534	needed \|= 1 << startBit;
535	}
536	}
537	}
538	for (i = 0; i < size_; i++) {
539	if ((i&needed) == needed) {
540	// this one may do
541	bool good = true;
542	for (int kk = 0; kk < numberActive_; kk++) {
543	int numberBits = numberBits_[kk];
544	int startBit = startBit_[kk];
545	int size = 1 << numberBits;
546	int start = 1 << startBit;
547	int mask = start * (size - 1);
548	int level = (i & mask) >> startBit;
549	if (level < lower[kk]) {
550	good = false;
551	break;
552	}
553	}
554	if (good && cost_[i] < bestValue) {
555	bestValue = cost_[i];
556	iBest = i;
557	}
558	}
559	}
560	delete [] lower;
561	}
562	if (bestValue < COIN_DBL_MAX) {
563	bestValue += fixedObj;
564	if (model_->messageHandler()->logLevel() > 1)
565	printf("Can get solution of %g\n", bestValue);
566	if (bestValue < model_->getMinimizationObjValue()) {
567	// set up solution
568	betterSolution = new double[numberColumns];
569	memcpy(betterSolution, lower, numberColumns*sizeof(double));
570	while (iBest > 0) {
571	int n = decodeBitPattern(iBest - back_[iBest], indices_, numberRows);
572	// Search for cheapest
573	double bestCost = COIN_DBL_MAX;
574	int iColumn = -1;
575	for (i = 0; i < numberColumns; i++) {
576	if (n == columnLength[i]) {
577	bool good = true;
578	for (int j = columnStart[i];
579	j < columnStart[i] + columnLength[i]; j++) {
580	int iRow = row[j];
581	double value = element[j];
582	int iValue = static_cast<int> (value);
583	if (iValue != indices_[iRow]) {
584	good = false;
585	break;
586	}
587	}
588	if (good && objective[i] < bestCost && betterSolution[i] < upper[i]) {
589	bestCost = objective[i];
590	iColumn = i;
591	}
592	}
593	}
594	assert (iColumn >= 0);
595	betterSolution[iColumn]++;
596	assert (betterSolution[iColumn] <= upper[iColumn]);
597	iBest = back_[iBest];
598	}
599	}
600	// paranoid check
601	double * rowActivity = new double [numberRows];
602	memset(rowActivity, 0, numberRows*sizeof(double));
603	for (i = 0; i < numberColumns; i++) {
604	int j;
605	double value = betterSolution[i];
606	if (value) {
607	for (j = columnStart[i];
608	j < columnStart[i] + columnLength[i]; j++) {
609	int iRow = row[j];
610	rowActivity[iRow] += value * element[j];
611	}
612	}
613	}
614	// check was feasible
615	bool feasible = true;
616	for (i = 0; i < numberRows; i++) {
617	if (rowActivity[i] < rowLower[i]) {
618	if (rowActivity[i] < rowLower[i] - 1.0e-8)
619	feasible = false;
620	} else if (rowActivity[i] > rowUpper[i]) {
621	if (rowActivity[i] > rowUpper[i] + 1.0e-8)
622	feasible = false;
623	}
624	}
625	if (feasible) {
626	if (model_->messageHandler()->logLevel() > 0)
627	printf("** good solution of %g by dynamic programming\n", bestValue);
628	}
629	delete [] rowActivity;
630	}
631	gutsOfDelete();
632	}
633	return returnCode;
634	}
635	/* Tries a column
636	returns true if was used in making any changes.
637	*/
638	bool
639	CbcFathomDynamicProgramming::tryColumn(int numberElements, const int * rows,
640	const double * coefficients, double cost,
641	int upper)
642	{
643	bool touched = false;
644	int n = 0;
645	if (algorithm_ == 0) {
646	for (int j = 0; j < numberElements; j++) {
647	int iRow = rows[j];
648	double value = coefficients[j];
649	int newRow = lookup_[iRow];
650	if (newRow < 0 \|\| value > rhs_[newRow]) {
651	n = 0;
652	break; //can't use
653	} else {
654	indices_[n++] = newRow;
655	}
656	}
657	if (n && upper) {
658	touched = addOneColumn0(n, indices_, cost);
659	}
660	} else {
661	for (int j = 0; j < numberElements; j++) {
662	int iRow = rows[j];
663	double value = coefficients[j];
664	int iValue = static_cast<int> (value);
665	int newRow = lookup_[iRow];
666	if (newRow < 0 \|\| iValue > rhs_[newRow]) {
667	n = 0;
668	break; //can't use
669	} else {
670	coefficients_[n] = iValue;
671	indices_[n++] = newRow;
672	if (upper*iValue > rhs_[newRow]) {
673	upper = rhs_[newRow] / iValue;
674	}
675	}
676	}
677	if (n) {
678	if (algorithm_ == 1) {
679	for (int k = 1; k <= upper; k++) {
680	bool t = addOneColumn1(n, indices_, coefficients_, cost);
681	if (t)
682	touched = true;
683	}
684	} else {
685	CoinSort_2(indices_, indices_ + n, coefficients_);
686	for (int k = 1; k <= upper; k++) {
687	bool t = addOneColumn1A(n, indices_, coefficients_, cost);
688	if (t)
689	touched = true;
690	}
691	}
692	}
693	}
694	return touched;
695	}
696	/* Adds one column if type 0,
697	returns true if was used in making any changes
698	*/
699	bool
700	CbcFathomDynamicProgramming::addOneColumn0(int numberElements, const int * rows,
701	double cost)
702	{
703	// build up mask
704	int mask = 0;
705	int i;
706	for (i = 0; i < numberElements; i++) {
707	int iRow = rows[i];
708	mask \|= 1 << iRow;
709	}
710	bitPattern_ = mask;
711	i = size_ - 1 - mask;
712	bool touched = false;
713	while (i >= 0) {
714	int kMask = i & mask;
715	if (kMask == 0) {
716	double thisCost = cost_[i];
717	if (thisCost != COIN_DBL_MAX) {
718	// possible
719	double newCost = thisCost + cost;
720	int next = i + mask;
721	if (cost_[next] > newCost) {
722	cost_[next] = newCost;
723	back_[next] = i;
724	touched = true;
725	}
726	}
727	i--;
728	} else {
729	// we can skip some
730	int k = (i&~mask);
731	#ifdef CBC_DEBUG
732	for (int j = i - 1; j > k; j--) {
733	int jMask = j & mask;
734	assert (jMask != 0);
735	}
736	#endif
737	i = k;
738	}
739	}
740	return touched;
741	}
742	/* Adds one attempt of one column of type 1,
743	returns true if was used in making any changes.
744	At present the user has to call it once for each possible value
745	*/
746	bool
747	CbcFathomDynamicProgramming::addOneColumn1(int numberElements, const int * rows,
748	const int * coefficients, double cost)
749	{
750	/* build up masks.
751	a) mask for 1 rhs
752	b) mask for addition
753	c) mask so adding will overflow
754	d) individual masks
755	*/
756	int mask1 = 0;
757	int maskAdd = 0;
758	int mask2 = 0;
759	int i;
760	int n2 = 0;
761	int mask[40];
762	int adjust[40];
763	assert (numberElements <= 40);
764	for (i = 0; i < numberElements; i++) {
765	int iRow = rows[i];
766	int numberBits = numberBits_[iRow];
767	int startBit = startBit_[iRow];
768	if (numberBits == 1) {
769	mask1 \|= 1 << startBit;
770	maskAdd \|= 1 << startBit;
771	mask2 \|= 1 << startBit;
772	} else {
773	int value = coefficients[i];
774	int size = 1 << numberBits;
775	int start = 1 << startBit;
776	assert (value < size);
777	maskAdd \|= start * value;
778	int gap = size - rhs_[iRow] - 1;
779	assert (gap >= 0);
780	int hi2 = rhs_[iRow] - value;
781	if (hi2 < size - 1)
782	hi2++;
783	adjust[n2] = start * hi2;
784	mask2 += start * gap;
785	mask[n2++] = start * (size - 1);
786	}
787	}
788	bitPattern_ = maskAdd;
789	i = size_ - 1 - maskAdd;
790	bool touched = false;
791	while (i >= 0) {
792	int kMask = i & mask1;
793	if (kMask == 0) {
794	bool good = true;
795	for (int kk = n2 - 1; kk >= 0; kk--) {
796	int iMask = mask[kk];
797	int jMask = iMask & mask2;
798	int kkMask = iMask & i;
799	kkMask += jMask;
800	if (kkMask > iMask) {
801	// we can skip some
802	int k = (i&~iMask);
803	k \|= adjust[kk];
804	#ifdef CBC_DEBUG
805	for (int j = i - 1; j > k; j--) {
806	int jMask = j & mask1;
807	if (jMask == 0) {
808	bool good = true;
809	for (int kk = n2 - 1; kk >= 0; kk--) {
810	int iMask = mask[kk];
811	int jMask = iMask & mask2;
812	int kkMask = iMask & i;
813	kkMask += jMask;
814	if (kkMask > iMask) {
815	good = false;
816	break;
817	}
818	}
819	assert (!good);
820	}
821	}
822	#endif
823	i = k;
824	good = false;
825	break;
826	}
827	}
828	if (good) {
829	double thisCost = cost_[i];
830	if (thisCost != COIN_DBL_MAX) {
831	// possible
832	double newCost = thisCost + cost;
833	int next = i + maskAdd;
834	if (cost_[next] > newCost) {
835	cost_[next] = newCost;
836	back_[next] = i;
837	touched = true;
838	}
839	}
840	}
841	i--;
842	} else {
843	// we can skip some
844	// we can skip some
845	int k = (i&~mask1);
846	#ifdef CBC_DEBUG
847	for (int j = i - 1; j > k; j--) {
848	int jMask = j & mask1;
849	assert (jMask != 0);
850	}
851	#endif
852	i = k;
853	}
854	}
855	return touched;
856	}
857	/* Adds one attempt of one column of type 1,
858	returns true if was used in making any changes.
859	At present the user has to call it once for each possible value
860	This version is when there are enough 1 rhs to do faster
861	*/
862	bool
863	CbcFathomDynamicProgramming::addOneColumn1A(int numberElements, const int * rows,
864	const int * coefficients, double cost)
865	{
866	/* build up masks.
867	a) mask for 1 rhs
868	b) mask for addition
869	c) mask so adding will overflow
870	d) mask for non 1 rhs
871	*/
872	int maskA = 0;
873	int maskAdd = 0;
874	int maskC = 0;
875	int maskD = 0;
876	int i;
877	for (i = 0; i < numberElements; i++) {
878	int iRow = rows[i];
879	int numberBits = numberBits_[iRow];
880	int startBit = startBit_[iRow];
881	if (numberBits == 1) {
882	maskA \|= 1 << startBit;
883	maskAdd \|= 1 << startBit;
884	} else {
885	int value = coefficients[i];
886	int size = 1 << numberBits;
887	int start = 1 << startBit;
888	assert (value < size);
889	maskAdd \|= start * value;
890	int gap = size - rhs_[iRow] + value - 1;
891	assert (gap > 0 && gap <= size - 1);
892	maskC \|= start * gap;
893	maskD \|= start * (size - 1);
894	}
895	}
896	bitPattern_ = maskAdd;
897	int maskDiff = maskD - maskC;
898	i = size_ - 1 - maskAdd;
899	bool touched = false;
900	if (!maskD) {
901	// Just ones
902	while (i >= 0) {
903	int kMask = i & maskA;
904	if (kMask == 0) {
905	double thisCost = cost_[i];
906	if (thisCost != COIN_DBL_MAX) {
907	// possible
908	double newCost = thisCost + cost;
909	int next = i + maskAdd;
910	if (cost_[next] > newCost) {
911	cost_[next] = newCost;
912	back_[next] = i;
913	touched = true;
914	}
915	}
916	i--;
917	} else {
918	// we can skip some
919	int k = (i&~maskA);
920	i = k;
921	}
922	}
923	} else {
924	// More general
925	while (i >= 0) {
926	int kMask = i & maskA;
927	if (kMask == 0) {
928	int added = i & maskD; // just bits belonging to non 1 rhs
929	added += maskC; // will overflow mask if bad
930	added &= (~maskD);
931	if (added == 0) {
932	double thisCost = cost_[i];
933	if (thisCost != COIN_DBL_MAX) {
934	// possible
935	double newCost = thisCost + cost;
936	int next = i + maskAdd;
937	if (cost_[next] > newCost) {
938	cost_[next] = newCost;
939	back_[next] = i;
940	touched = true;
941	}
942	}
943	i--;
944	} else {
945	// we can skip some
946	int k = i & ~ maskD; // clear all
947	// Put back enough - but only below where we are
948	int kk = (numberNonOne_ << 1) - 2;
949	assert (rhs_[kk] > 1);
950	int iMask = 0;
951	for (; kk >= 0; kk -= 2) {
952	iMask = 1 << startBit_[kk+1];
953	if ((added&iMask) != 0) {
954	iMask--;
955	break;
956	}
957	}
958	assert (kk >= 0);
959	iMask &= maskDiff;
960	k \|= iMask;
961	assert (k < i);
962	i = k;
963	}
964	} else {
965	// we can skip some
966	int k = (i&~maskA);
967	i = k;
968	}
969	}
970	}
971	return touched;
972	}
973	// update model
974	void CbcFathomDynamicProgramming::setModel(CbcModel * model)
975	{
976	model_ = model;
977	type_ = checkPossible();
978	}
979	// Gets bit pattern from original column
980	int CbcFathomDynamicProgramming::bitPattern(int numberElements, const int * rows,
981	const int * coefficients)
982	{
983	int i;
984	int mask = 0;
985	switch (algorithm_) {
986	// just ones
987	case 0:
988	for (i = 0; i < numberElements; i++) {
989	int iRow = rows[i];
990	iRow = lookup_[iRow];
991	if (iRow >= 0)
992	mask \|= 1 << iRow;
993	}
994	break;
995	//
996	case 1:
997	case 2:
998	for (i = 0; i < numberElements; i++) {
999	int iRow = rows[i];
1000	iRow = lookup_[iRow];
1001	if (iRow >= 0) {
1002	int startBit = startBit_[iRow];
1003	int value = coefficients[i];
1004	int start = 1 << startBit;
1005	mask \|= start * value;
1006	}
1007	}
1008	break;
1009	}
1010	return mask;
1011	}
1012	// Fills in original column (dense) from bit pattern
1013	int CbcFathomDynamicProgramming::decodeBitPattern(int bitPattern,
1014	int * values,
1015	int numberRows)
1016	{
1017	int i;
1018	int n = 0;
1019	switch (algorithm_) {
1020	// just ones
1021	case 0:
1022	for (i = 0; i < numberRows; i++) {
1023	values[i] = 0;
1024	int iRow = lookup_[i];
1025	if (iRow >= 0) {
1026	if ((bitPattern&(1 << iRow)) != 0) {
1027	values[i] = 1;
1028	n++;
1029	}
1030	}
1031	}
1032	break;
1033	//
1034	case 1:
1035	case 2:
1036	for (i = 0; i < numberRows; i++) {
1037	values[i] = 0;
1038	int iRow = lookup_[i];
1039	if (iRow >= 0) {
1040	int startBit = startBit_[iRow];
1041	int numberBits = numberBits_[iRow];
1042	int iValue = bitPattern >> startBit;
1043	iValue &= ((1 << numberBits) - 1);
1044	if (iValue) {
1045	values[i] = iValue;
1046	n++;
1047	}
1048	}
1049	}
1050	break;
1051	}
1052	return n;
1053	}
1054
1055

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

Download in other formats: