Context Navigation

source: stable/2.9/Cbc/src/CbcTreeLocal.cpp @ 2180

Last change on this file since 2180 was 1839, checked in by forrest, 7 years ago
multiple root solvers, stronger strong branching and cutoff as constraint
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 65.3 KB

Line
1	/* $Id: CbcTreeLocal.cpp 1839 2013-01-16 18:41:25Z forrest $ */
2	// Copyright (C) 2004, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6	#include "CbcModel.hpp"
7	#include "CbcNode.hpp"
8	#include "CbcTreeLocal.hpp"
9	#include "CoinPackedMatrix.hpp"
10	#include "CoinTime.hpp"
11	#include "OsiRowCutDebugger.hpp"
12	#include <cassert>
13	#ifdef JJF_ZERO
14	// gdb doesn't always put breakpoints in this virtual function
15	// just stick xxxxxx() where you want to start
16	static void xxxxxx()
17	{
18	printf("break\n");
19	}
20	#endif
21	CbcTreeLocal::CbcTreeLocal()
22	: localNode_(NULL),
23	bestSolution_(NULL),
24	savedSolution_(NULL),
25	saveNumberSolutions_(0),
26	model_(NULL),
27	originalLower_(NULL),
28	originalUpper_(NULL),
29	range_(0),
30	typeCuts_(-1),
31	maxDiversification_(0),
32	diversification_(0),
33	nextStrong_(false),
34	rhs_(0.0),
35	savedGap_(0.0),
36	bestCutoff_(0.0),
37	timeLimit_(0),
38	startTime_(0),
39	nodeLimit_(0),
40	startNode_(-1),
41	searchType_(-1),
42	refine_(false)
43	{
44
45	}
46	/* Constructor with solution.
47	range is upper bound on difference from given solution.
48	maxDiversification is maximum number of diversifications to try
49	timeLimit is seconds in subTree
50	nodeLimit is nodes in subTree
51	*/
52	CbcTreeLocal::CbcTreeLocal(CbcModel * model, const double * solution ,
53	int range, int typeCuts, int maxDiversification,
54	int timeLimit, int nodeLimit, bool refine)
55	: localNode_(NULL),
56	bestSolution_(NULL),
57	savedSolution_(NULL),
58	saveNumberSolutions_(0),
59	model_(model),
60	originalLower_(NULL),
61	originalUpper_(NULL),
62	range_(range),
63	typeCuts_(typeCuts),
64	maxDiversification_(maxDiversification),
65	diversification_(0),
66	nextStrong_(false),
67	rhs_(0.0),
68	savedGap_(0.0),
69	bestCutoff_(0.0),
70	timeLimit_(timeLimit),
71	startTime_(0),
72	nodeLimit_(nodeLimit),
73	startNode_(-1),
74	searchType_(-1),
75	refine_(refine)
76	{
77
78	OsiSolverInterface * solver = model_->solver();
79	const double * lower = solver->getColLower();
80	const double * upper = solver->getColUpper();
81	//const double * solution = solver->getColSolution();
82	//const double * objective = solver->getObjCoefficients();
83	double primalTolerance;
84	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
85
86	// Get increment
87	model_->analyzeObjective();
88
89	{
90	// needed to sync cutoffs
91	double value ;
92	solver->getDblParam(OsiDualObjectiveLimit, value) ;
93	model_->setCutoff(value * solver->getObjSense());
94	}
95	bestCutoff_ = model_->getCutoff();
96	// save current gap
97	savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
98
99	// make sure integers found
100	model_->findIntegers(false);
101	int numberIntegers = model_->numberIntegers();
102	const int * integerVariable = model_->integerVariable();
103	int i;
104	double direction = solver->getObjSense();
105	double newSolutionValue = 1.0e50;
106	if (solution) {
107	// copy solution
108	solver->setColSolution(solution);
109	newSolutionValue = direction * solver->getObjValue();
110	}
111	originalLower_ = new double [numberIntegers];
112	originalUpper_ = new double [numberIntegers];
113	bool all01 = true;
114	int number01 = 0;
115	for (i = 0; i < numberIntegers; i++) {
116	int iColumn = integerVariable[i];
117	originalLower_[i] = lower[iColumn];
118	originalUpper_[i] = upper[iColumn];
119	if (upper[iColumn] - lower[iColumn] > 1.5)
120	all01 = false;
121	else if (upper[iColumn] - lower[iColumn] == 1.0)
122	number01++;
123	}
124	if (all01 && !typeCuts_)
125	typeCuts_ = 1; // may as well so we don't have to deal with refine
126	if (!number01 && !typeCuts_) {
127	if (model_->messageHandler()->logLevel() > 1)
128	printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
129	typeCuts_ = -1;
130	} else {
131	if (model_->messageHandler()->logLevel() > 1) {
132	std::string type;
133	if (all01) {
134	printf("%d 0-1 variables normal local cuts\n",
135	number01);
136	} else if (typeCuts_) {
137	printf("%d 0-1 variables, %d other - general integer local cuts\n",
138	number01, numberIntegers - number01);
139	} else {
140	printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
141	number01, numberIntegers - number01);
142	}
143	printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
144	maxDiversification_, range_, timeLimit_, nodeLimit_);
145	}
146	}
147	int numberColumns = model_->getNumCols();
148	savedSolution_ = new double [numberColumns];
149	memset(savedSolution_, 0, numberColumns*sizeof(double));
150	if (solution) {
151	rhs_ = range_;
152	// Check feasible
153	int goodSolution = createCut(solution, cut_);
154	if (goodSolution >= 0) {
155	for (i = 0; i < numberIntegers; i++) {
156	int iColumn = integerVariable[i];
157	double value = floor(solution[iColumn] + 0.5);
158	// fix so setBestSolution will work
159	solver->setColLower(iColumn, value);
160	solver->setColUpper(iColumn, value);
161	}
162	model_->reserveCurrentSolution();
163	// Create cut and get total gap
164	if (newSolutionValue < bestCutoff_) {
165	model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
166	bestCutoff_ = model_->getCutoff();
167	// save as best solution
168	memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
169	}
170	for (i = 0; i < numberIntegers; i++) {
171	int iColumn = integerVariable[i];
172	// restore bounds
173	solver->setColLower(iColumn, originalLower_[i]);
174	solver->setColUpper(iColumn, originalUpper_[i]);
175	}
176	// make sure can't stop on gap
177	model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
178	} else {
179	model_ = NULL;
180	}
181	} else {
182	// no solution
183	rhs_ = 1.0e50;
184	// make sure can't stop on gap
185	model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
186	}
187	}
188	CbcTreeLocal::~CbcTreeLocal()
189	{
190	delete [] originalLower_;
191	delete [] originalUpper_;
192	delete [] bestSolution_;
193	delete [] savedSolution_;
194	delete localNode_;
195	}
196	// Copy constructor
197	CbcTreeLocal::CbcTreeLocal ( const CbcTreeLocal & rhs)
198	: CbcTree(rhs),
199	saveNumberSolutions_(rhs.saveNumberSolutions_),
200	model_(rhs.model_),
201	range_(rhs.range_),
202	typeCuts_(rhs.typeCuts_),
203	maxDiversification_(rhs.maxDiversification_),
204	diversification_(rhs.diversification_),
205	nextStrong_(rhs.nextStrong_),
206	rhs_(rhs.rhs_),
207	savedGap_(rhs.savedGap_),
208	bestCutoff_(rhs.bestCutoff_),
209	timeLimit_(rhs.timeLimit_),
210	startTime_(rhs.startTime_),
211	nodeLimit_(rhs.nodeLimit_),
212	startNode_(rhs.startNode_),
213	searchType_(rhs.searchType_),
214	refine_(rhs.refine_)
215	{
216	cut_ = rhs.cut_;
217	fixedCut_ = rhs.fixedCut_;
218	if (rhs.localNode_)
219	localNode_ = new CbcNode(*rhs.localNode_);
220	else
221	localNode_ = NULL;
222	if (rhs.originalLower_) {
223	int numberIntegers = model_->numberIntegers();
224	originalLower_ = new double [numberIntegers];
225	memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
226	originalUpper_ = new double [numberIntegers];
227	memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
228	} else {
229	originalLower_ = NULL;
230	originalUpper_ = NULL;
231	}
232	if (rhs.bestSolution_) {
233	int numberColumns = model_->getNumCols();
234	bestSolution_ = new double [numberColumns];
235	memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
236	} else {
237	bestSolution_ = NULL;
238	}
239	if (rhs.savedSolution_) {
240	int numberColumns = model_->getNumCols();
241	savedSolution_ = new double [numberColumns];
242	memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
243	} else {
244	savedSolution_ = NULL;
245	}
246	}
247	//----------------------------------------------------------------
248	// Assignment operator
249	//-------------------------------------------------------------------
250	CbcTreeLocal &
251	CbcTreeLocal::operator=(const CbcTreeLocal & rhs)
252	{
253	if (this != &rhs) {
254	CbcTree::operator=(rhs);
255	saveNumberSolutions_ = rhs.saveNumberSolutions_;
256	cut_ = rhs.cut_;
257	fixedCut_ = rhs.fixedCut_;
258	delete localNode_;
259	if (rhs.localNode_)
260	localNode_ = new CbcNode(*rhs.localNode_);
261	else
262	localNode_ = NULL;
263	model_ = rhs.model_;
264	range_ = rhs.range_;
265	typeCuts_ = rhs.typeCuts_;
266	maxDiversification_ = rhs.maxDiversification_;
267	diversification_ = rhs.diversification_;
268	nextStrong_ = rhs.nextStrong_;
269	rhs_ = rhs.rhs_;
270	savedGap_ = rhs.savedGap_;
271	bestCutoff_ = rhs.bestCutoff_;
272	timeLimit_ = rhs.timeLimit_;
273	startTime_ = rhs.startTime_;
274	nodeLimit_ = rhs.nodeLimit_;
275	startNode_ = rhs.startNode_;
276	searchType_ = rhs.searchType_;
277	refine_ = rhs.refine_;
278	delete [] originalLower_;
279	delete [] originalUpper_;
280	if (rhs.originalLower_) {
281	int numberIntegers = model_->numberIntegers();
282	originalLower_ = new double [numberIntegers];
283	memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
284	originalUpper_ = new double [numberIntegers];
285	memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
286	} else {
287	originalLower_ = NULL;
288	originalUpper_ = NULL;
289	}
290	delete [] bestSolution_;
291	if (rhs.bestSolution_) {
292	int numberColumns = model_->getNumCols();
293	bestSolution_ = new double [numberColumns];
294	memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
295	} else {
296	bestSolution_ = NULL;
297	}
298	delete [] savedSolution_;
299	if (rhs.savedSolution_) {
300	int numberColumns = model_->getNumCols();
301	savedSolution_ = new double [numberColumns];
302	memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
303	} else {
304	savedSolution_ = NULL;
305	}
306	}
307	return *this;
308	}
309	// Clone
310	CbcTree *
311	CbcTreeLocal::clone() const
312	{
313	return new CbcTreeLocal(*this);
314	}
315	// Pass in solution (so can be used after heuristic)
316	void
317	CbcTreeLocal::passInSolution(const double * solution, double solutionValue)
318	{
319	int numberColumns = model_->getNumCols();
320	delete [] savedSolution_;
321	savedSolution_ = new double [numberColumns];
322	memcpy(savedSolution_, solution, numberColumns*sizeof(double));
323	rhs_ = range_;
324	// Check feasible
325	int goodSolution = createCut(solution, cut_);
326	if (goodSolution >= 0) {
327	bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
328	} else {
329	model_ = NULL;
330	}
331	}
332	// Return the top node of the heap
333	CbcNode *
334	CbcTreeLocal::top() const
335	{
336	#ifdef CBC_DEBUG
337	int smallest = 9999999;
338	int largest = -1;
339	double smallestD = 1.0e30;
340	double largestD = -1.0e30;
341	int n = nodes_.size();
342	for (int i = 0; i < n; i++) {
343	int nn = nodes_[i]->nodeInfo()->nodeNumber();
344	double dd = nodes_[i]->objectiveValue();
345	largest = CoinMax(largest, nn);
346	smallest = CoinMin(smallest, nn);
347	largestD = CoinMax(largestD, dd);
348	smallestD = CoinMin(smallestD, dd);
349	}
350	if (model_->messageHandler()->logLevel() > 1) {
351	printf("smallest %d, largest %d, top %d\n", smallest, largest,
352	nodes_.front()->nodeInfo()->nodeNumber());
353	printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
354	}
355	#endif
356	return nodes_.front();
357	}
358
359	// Add a node to the heap
360	void
361	CbcTreeLocal::push(CbcNode * x)
362	{
363	if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
364	startNode_ = model_->getNodeCount();
365	// save copy of node
366	localNode_ = new CbcNode(*x);
367
368	if (cut_.row().getNumElements()) {
369	// Add to global cuts
370	// we came in with solution
371	model_->makeGlobalCut(cut_);
372	if (model_->messageHandler()->logLevel() > 1)
373	printf("initial cut - rhs %g %g\n",
374	cut_.lb(), cut_.ub());
375	searchType_ = 1;
376	} else {
377	// stop on first solution
378	searchType_ = 0;
379	}
380	startTime_ = static_cast<int> (CoinCpuTime());
381	saveNumberSolutions_ = model_->getSolutionCount();
382	}
383	nodes_.push_back(x);
384	#ifdef CBC_DEBUG
385	if (model_->messageHandler()->logLevel() > 0)
386	printf("pushing node onto heap %d %x %x\n",
387	x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
388	#endif
389	std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
390	}
391
392	// Remove the top node from the heap
393	void
394	CbcTreeLocal::pop()
395	{
396	std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
397	nodes_.pop_back();
398	}
399	// Test if empty - does work if so
400	bool
401	CbcTreeLocal::empty()
402	{
403	if (typeCuts_ < 0)
404	return !nodes_.size();
405	/* state -
406	0 iterating
407	1 subtree finished optimal solution for subtree found
408	2 subtree finished and no solution found
409	3 subtree exiting and solution found
410	4 subtree exiting and no solution found
411	*/
412	int state = 0;
413	assert (searchType_ != 2);
414	if (searchType_) {
415	if (CoinCpuTime() - startTime_ > timeLimit_ \|\| model_->getNodeCount() - startNode_ >= nodeLimit_) {
416	state = 4;
417	}
418	} else {
419	if (model_->getSolutionCount() > saveNumberSolutions_) {
420	state = 4;
421	}
422	}
423	if (!nodes_.size())
424	state = 2;
425	if (!state) {
426	return false;
427	}
428	// Finished this phase
429	int numberColumns = model_->getNumCols();
430	if (model_->getSolutionCount() > saveNumberSolutions_) {
431	if (model_->getCutoff() < bestCutoff_) {
432	// Save solution
433	if (!bestSolution_)
434	bestSolution_ = new double [numberColumns];
435	memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
436	bestCutoff_ = model_->getCutoff();
437	}
438	state--;
439	}
440	// get rid of all nodes (safe even if already done)
441	double bestPossibleObjective;
442	cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
443
444	double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
445	if (model_->messageHandler()->logLevel() > 1)
446	printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
447	state,
448	model_->getNodeCount() - startNode_,
449	static_cast<int> (CoinCpuTime()) - startTime_,
450	model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
451	saveNumberSolutions_ = model_->getSolutionCount();
452	bool finished = false;
453	bool lastTry = false;
454	switch (state) {
455	case 1:
456	// solution found and subtree exhausted
457	if (rhs_ > 1.0e30) {
458	finished = true;
459	} else {
460	// find global cut and reverse
461	reverseCut(1);
462	searchType_ = 1; // first false
463	rhs_ = range_; // reset range
464	nextStrong_ = false;
465
466	// save best solution in this subtree
467	memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
468	}
469	break;
470	case 2:
471	// solution not found and subtree exhausted
472	if (rhs_ > 1.0e30) {
473	finished = true;
474	} else {
475	// find global cut and reverse
476	reverseCut(2);
477	searchType_ = 1; // first false
478	if (diversification_ < maxDiversification_) {
479	if (nextStrong_) {
480	diversification_++;
481	// cut is valid so don't model_->setCutoff(1.0e50);
482	searchType_ = 0;
483	}
484	nextStrong_ = true;
485	rhs_ += range_ / 2;
486	} else {
487	// This will be last try (may hit max time)
488	lastTry = true;
489	if (!maxDiversification_)
490	typeCuts_ = -1; // make sure can't start again
491	model_->setCutoff(bestCutoff_);
492	if (model_->messageHandler()->logLevel() > 1)
493	printf("Exiting local search with current set of cuts\n");
494	rhs_ = 1.0e100;
495	// Can now stop on gap
496	model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
497	}
498	}
499	break;
500	case 3:
501	// solution found and subtree not exhausted
502	if (rhs_ < 1.0e30) {
503	if (searchType_) {
504	if (!typeCuts_ && refine_ && searchType_ == 1) {
505	// We need to check we have best solution given these 0-1 values
506	OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
507	CbcModel * subModel = model_->subTreeModel(subSolver);
508	CbcTree normalTree;
509	subModel->passInTreeHandler(normalTree);
510	int numberIntegers = model_->numberIntegers();
511	const int * integerVariable = model_->integerVariable();
512	const double * solution = model_->bestSolution();
513	int i;
514	int numberColumns = model_->getNumCols();
515	for (i = 0; i < numberIntegers; i++) {
516	int iColumn = integerVariable[i];
517	double value = floor(solution[iColumn] + 0.5);
518	if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
519	continue; // skip as not 0-1
520	if (originalLower_[i] == originalUpper_[i])
521	continue;
522	subSolver->setColLower(iColumn, value);
523	subSolver->setColUpper(iColumn, value);
524	}
525	subSolver->initialSolve();
526	// We can copy cutoff
527	// But adjust
528	subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
529	subModel->setSolutionCount(0);
530	assert (subModel->isProvenOptimal());
531	if (!subModel->typePresolve()) {
532	subModel->branchAndBound();
533	if (subModel->status()) {
534	model_->incrementSubTreeStopped();
535	}
536	//printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
537	// subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
538	double newCutoff = subModel->getMinimizationObjValue() -
539	subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
540	if (subModel->getSolutionCount()) {
541	if (!subModel->status())
542	assert (subModel->isProvenOptimal());
543	memcpy(model_->bestSolution(), subModel->bestSolution(),
544	numberColumns*sizeof(double));
545	model_->setCutoff(newCutoff);
546	}
547	} else if (subModel->typePresolve() == 1) {
548	CbcModel * model2 = subModel->integerPresolve(true);
549	if (model2) {
550	// Do complete search
551	model2->branchAndBound();
552	// get back solution
553	subModel->originalModel(model2, false);
554	if (model2->status()) {
555	model_->incrementSubTreeStopped();
556	}
557	double newCutoff = model2->getMinimizationObjValue() -
558	model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
559	if (model2->getSolutionCount()) {
560	if (!model2->status())
561	assert (model2->isProvenOptimal());
562	memcpy(model_->bestSolution(), subModel->bestSolution(),
563	numberColumns*sizeof(double));
564	model_->setCutoff(newCutoff);
565	}
566	delete model2;
567	} else {
568	// infeasible - could just be - due to cutoff
569	}
570	} else {
571	// too dangerous at present
572	assert (subModel->typePresolve() != 2);
573	}
574	if (model_->getCutoff() < bestCutoff_) {
575	// Save solution
576	if (!bestSolution_)
577	bestSolution_ = new double [numberColumns];
578	memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
579	bestCutoff_ = model_->getCutoff();
580	}
581	delete subModel;
582	}
583	// we have done search to make sure best general solution
584	searchType_ = 1;
585	// Reverse cut weakly
586	reverseCut(3, rhs_);
587	} else {
588	searchType_ = 1;
589	// delete last cut
590	deleteCut(cut_);
591	}
592	} else {
593	searchType_ = 1;
594	}
595	// save best solution in this subtree
596	memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
597	nextStrong_ = false;
598	rhs_ = range_;
599	break;
600	case 4:
601	// solution not found and subtree not exhausted
602	if (maxDiversification_) {
603	if (nextStrong_) {
604	// Reverse cut weakly
605	reverseCut(4, rhs_);
606	model_->setCutoff(1.0e50);
607	diversification_++;
608	searchType_ = 0;
609	} else {
610	// delete last cut
611	deleteCut(cut_);
612	searchType_ = 1;
613	}
614	nextStrong_ = true;
615	rhs_ += range_ / 2;
616	} else {
617	// special case when using as heuristic
618	// Reverse cut weakly if lb -infinity
619	reverseCut(4, rhs_);
620	// This will be last try (may hit max time0
621	lastTry = true;
622	model_->setCutoff(bestCutoff_);
623	if (model_->messageHandler()->logLevel() > 1)
624	printf("Exiting local search with current set of cuts\n");
625	rhs_ = 1.0e100;
626	// Can now stop on gap
627	model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
628	typeCuts_ = -1;
629	}
630	break;
631	}
632	if (rhs_ < 1.0e30 \|\| lastTry) {
633	int goodSolution = createCut(savedSolution_, cut_);
634	if (goodSolution >= 0) {
635	// Add to global cuts
636	model_->makeGlobalCut(cut_);
637	CbcRowCuts * global = model_->globalCuts();
638	int n = global->sizeRowCuts();
639	OsiRowCut * rowCut = global->rowCutPtr(n - 1);
640	if (model_->messageHandler()->logLevel() > 1)
641	printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
642	n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
643	const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
644	if (debugger) {
645	if (debugger->invalidCut(*rowCut))
646	printf("ZZZZTree Global cut - cuts off optimal solution!\n");
647	}
648	for (int i = 0; i < n; i++) {
649	rowCut = global->rowCutPtr(i);
650	if (model_->messageHandler()->logLevel() > 0)
651	printf("%d - rhs %g %g\n",
652	i, rowCut->lb(), rowCut->ub());
653	}
654	}
655	// put back node
656	startTime_ = static_cast<int> (CoinCpuTime());
657	startNode_ = model_->getNodeCount();
658	if (localNode_) {
659	// save copy of node
660	CbcNode * localNode2 = new CbcNode(*localNode_);
661	// But localNode2 now owns cuts so swap
662	//printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
663	// localNode_,localNode_->nodeInfo());
664	nodes_.push_back(localNode_);
665	localNode_ = localNode2;
666	std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
667	}
668	}
669	return finished;
670	}
671	// We may have got an intelligent tree so give it one more chance
672	void
673	CbcTreeLocal::endSearch()
674	{
675	if (typeCuts_ >= 0) {
676	// copy best solution to model
677	int numberColumns = model_->getNumCols();
678	if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
679	memcpy(model_->bestSolution(), bestSolution_, numberColumns*sizeof(double));
680	model_->setCutoff(bestCutoff_);
681	// recompute objective value
682	const double * objCoef = model_->getObjCoefficients();
683	double objOffset = 0.0;
684	model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
685
686	// Compute dot product of objCoef and colSol and then adjust by offset
687	double objValue = -objOffset;
688	for ( int i = 0 ; i < numberColumns ; i++ )
689	objValue += objCoef[i] * bestSolution_[i];
690	model_->setMinimizationObjValue(objValue);
691	}
692	// Can now stop on gap
693	model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
694	}
695	}
696	// Create cut
697	int
698	CbcTreeLocal::createCut(const double * solution, OsiRowCut & rowCut)
699	{
700	if (rhs_ > 1.0e20)
701	return -1;
702	OsiSolverInterface * solver = model_->solver();
703	const double * rowLower = solver->getRowLower();
704	const double * rowUpper = solver->getRowUpper();
705	//const double * solution = solver->getColSolution();
706	//const double * objective = solver->getObjCoefficients();
707	double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
708	double primalTolerance;
709	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
710	// relax
711	primalTolerance *= 1000.0;
712
713	int numberRows = model_->getNumRows();
714
715	int numberIntegers = model_->numberIntegers();
716	const int * integerVariable = model_->integerVariable();
717	int i;
718
719	// Check feasible
720
721	double * rowActivity = new double[numberRows];
722	memset(rowActivity, 0, numberRows*sizeof(double));
723	solver->getMatrixByCol()->times(solution, rowActivity) ;
724	int goodSolution = 0;
725	// check was feasible
726	for (i = 0; i < numberRows; i++) {
727	if (rowActivity[i] < rowLower[i] - primalTolerance) {
728	goodSolution = -1;
729	} else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
730	goodSolution = -1;
731	}
732	}
733	delete [] rowActivity;
734	for (i = 0; i < numberIntegers; i++) {
735	int iColumn = integerVariable[i];
736	double value = solution[iColumn];
737	if (fabs(floor(value + 0.5) - value) > integerTolerance) {
738	goodSolution = -1;
739	}
740	}
741	// zap cut
742	if (goodSolution == 0) {
743	// Create cut and get total gap
744	CoinPackedVector cut;
745	double rhs = rhs_;
746	double maxValue = 0.0;
747	for (i = 0; i < numberIntegers; i++) {
748	int iColumn = integerVariable[i];
749	double value = floor(solution[iColumn] + 0.5);
750	/*
751	typeCuts_ == 0 restricts to binary, 1 allows general integer. But we're
752	still restricted to being up against a bound. Consider: the notion is that
753	the cut restricts us to a k-neighbourhood. For binary variables, this
754	amounts to k variables which change value. For general integer, we could
755	end up with a single variable sucking up all of k (hence mu --- the
756	variable must swing to its other bound to look like a movement of 1). For
757	variables in the middle of a range, we're talking about fabs(sol<j> - x<j>).
758	*/
759	if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
760	continue; // skip as not 0-1
761	if (originalLower_[i] == originalUpper_[i])
762	continue;
763	double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
764	if (value == originalLower_[i]) {
765	rhs += mu * originalLower_[i];
766	cut.insert(iColumn, 1.0);
767	maxValue += originalUpper_[i];
768	} else if (value == originalUpper_[i]) {
769	rhs -= mu * originalUpper_[i];
770	cut.insert(iColumn, -1.0);
771	maxValue += originalLower_[i];
772	}
773	}
774	if (maxValue < rhs - primalTolerance) {
775	if (model_->messageHandler()->logLevel() > 1)
776	printf("slack cut\n");
777	goodSolution = 1;
778	}
779	rowCut.setRow(cut);
780	rowCut.setLb(-COIN_DBL_MAX);
781	rowCut.setUb(rhs);
782	rowCut.setGloballyValid();
783	if (model_->messageHandler()->logLevel() > 1)
784	printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
785	#ifdef CBC_DEBUG
786	if (model_->messageHandler()->logLevel() > 0) {
787	int k;
788	for (k = 0; k < cut.getNumElements(); k++) {
789	printf("%i %g ", cut.getIndices()[k], cut.getElements()[k]);
790	if ((k + 1) % 5 == 0)
791	printf("\n");
792	}
793	if (k % 5 != 0)
794	printf("\n");
795	}
796	#endif
797	return goodSolution;
798	} else {
799	if (model_->messageHandler()->logLevel() > 1)
800	printf("Not a good solution\n");
801	return -1;
802	}
803	}
804	// Other side of last cut branch
805	void
806	CbcTreeLocal::reverseCut(int state, double bias)
807	{
808	// find global cut
809	CbcRowCuts * global = model_->globalCuts();
810	int n = global->sizeRowCuts();
811	int i;
812	OsiRowCut * rowCut = NULL;
813	for ( i = 0; i < n; i++) {
814	rowCut = global->rowCutPtr(i);
815	if (cut_ == *rowCut) {
816	break;
817	}
818	}
819	if (!rowCut) {
820	// must have got here in odd way e.g. strong branching
821	return;
822	}
823	if (rowCut->lb() > -1.0e10)
824	return;
825	// get smallest element
826	double smallest = COIN_DBL_MAX;
827	CoinPackedVector row = cut_.row();
828	for (int k = 0; k < row.getNumElements(); k++)
829	smallest = CoinMin(smallest, fabs(row.getElements()[k]));
830	if (!typeCuts_ && !refine_) {
831	// Reverse cut very very weakly
832	if (state > 2)
833	smallest = 0.0;
834	}
835	// replace by other way
836	if (model_->messageHandler()->logLevel() > 1)
837	printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
838	i, n, rowCut->lb(), rowCut->ub());
839	rowCut->setLb(rowCut->ub() + smallest - bias);
840	rowCut->setUb(COIN_DBL_MAX);
841	if (model_->messageHandler()->logLevel() > 1)
842	printf("new rhs %g %g, bias %g smallest %g ",
843	rowCut->lb(), rowCut->ub(), bias, smallest);
844	const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
845	if (debugger) {
846	if (debugger->invalidCut(*rowCut))
847	printf("ZZZZTree Global cut - cuts off optimal solution!\n");
848	}
849	}
850	// Delete last cut branch
851	void
852	CbcTreeLocal::deleteCut(OsiRowCut & cut)
853	{
854	// find global cut
855	CbcRowCuts * global = model_->globalCuts();
856	int n = global->sizeRowCuts();
857	int i;
858	OsiRowCut * rowCut = NULL;
859	for ( i = 0; i < n; i++) {
860	rowCut = global->rowCutPtr(i);
861	if (cut == *rowCut) {
862	break;
863	}
864	}
865	assert (i < n);
866	// delete last cut
867	if (model_->messageHandler()->logLevel() > 1)
868	printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
869	i, n, rowCut->lb(), rowCut->ub());
870	global->eraseRowCut(i);
871	}
872	// Create C++ lines to get to current state
873	void
874	CbcTreeLocal::generateCpp( FILE * fp)
875	{
876	CbcTreeLocal other;
877	fprintf(fp, "0#include \"CbcTreeLocal.hpp\"\n");
878	fprintf(fp, "5 CbcTreeLocal localTree(cbcModel,NULL);\n");
879	if (range_ != other.range_)
880	fprintf(fp, "5 localTree.setRange(%d);\n", range_);
881	if (typeCuts_ != other.typeCuts_)
882	fprintf(fp, "5 localTree.setTypeCuts(%d);\n", typeCuts_);
883	if (maxDiversification_ != other.maxDiversification_)
884	fprintf(fp, "5 localTree.setMaxDiversification(%d);\n", maxDiversification_);
885	if (timeLimit_ != other.timeLimit_)
886	fprintf(fp, "5 localTree.setTimeLimit(%d);\n", timeLimit_);
887	if (nodeLimit_ != other.nodeLimit_)
888	fprintf(fp, "5 localTree.setNodeLimit(%d);\n", nodeLimit_);
889	if (refine_ != other.refine_)
890	fprintf(fp, "5 localTree.setRefine(%s);\n", refine_ ? "true" : "false");
891	fprintf(fp, "5 cbcModel->passInTreeHandler(localTree);\n");
892	}
893
894
895	CbcTreeVariable::CbcTreeVariable()
896	: localNode_(NULL),
897	bestSolution_(NULL),
898	savedSolution_(NULL),
899	saveNumberSolutions_(0),
900	model_(NULL),
901	originalLower_(NULL),
902	originalUpper_(NULL),
903	range_(0),
904	typeCuts_(-1),
905	maxDiversification_(0),
906	diversification_(0),
907	nextStrong_(false),
908	rhs_(0.0),
909	savedGap_(0.0),
910	bestCutoff_(0.0),
911	timeLimit_(0),
912	startTime_(0),
913	nodeLimit_(0),
914	startNode_(-1),
915	searchType_(-1),
916	refine_(false)
917	{
918
919	}
920	/* Constructor with solution.
921	range is upper bound on difference from given solution.
922	maxDiversification is maximum number of diversifications to try
923	timeLimit is seconds in subTree
924	nodeLimit is nodes in subTree
925	*/
926	CbcTreeVariable::CbcTreeVariable(CbcModel * model, const double * solution ,
927	int range, int typeCuts, int maxDiversification,
928	int timeLimit, int nodeLimit, bool refine)
929	: localNode_(NULL),
930	bestSolution_(NULL),
931	savedSolution_(NULL),
932	saveNumberSolutions_(0),
933	model_(model),
934	originalLower_(NULL),
935	originalUpper_(NULL),
936	range_(range),
937	typeCuts_(typeCuts),
938	maxDiversification_(maxDiversification),
939	diversification_(0),
940	nextStrong_(false),
941	rhs_(0.0),
942	savedGap_(0.0),
943	bestCutoff_(0.0),
944	timeLimit_(timeLimit),
945	startTime_(0),
946	nodeLimit_(nodeLimit),
947	startNode_(-1),
948	searchType_(-1),
949	refine_(refine)
950	{
951
952	OsiSolverInterface * solver = model_->solver();
953	const double * lower = solver->getColLower();
954	const double * upper = solver->getColUpper();
955	//const double * solution = solver->getColSolution();
956	//const double * objective = solver->getObjCoefficients();
957	double primalTolerance;
958	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
959
960	// Get increment
961	model_->analyzeObjective();
962
963	{
964	// needed to sync cutoffs
965	double value ;
966	solver->getDblParam(OsiDualObjectiveLimit, value) ;
967	model_->setCutoff(value * solver->getObjSense());
968	}
969	bestCutoff_ = model_->getCutoff();
970	// save current gap
971	savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
972
973	// make sure integers found
974	model_->findIntegers(false);
975	int numberIntegers = model_->numberIntegers();
976	const int * integerVariable = model_->integerVariable();
977	int i;
978	double direction = solver->getObjSense();
979	double newSolutionValue = 1.0e50;
980	if (solution) {
981	// copy solution
982	solver->setColSolution(solution);
983	newSolutionValue = direction * solver->getObjValue();
984	}
985	originalLower_ = new double [numberIntegers];
986	originalUpper_ = new double [numberIntegers];
987	bool all01 = true;
988	int number01 = 0;
989	for (i = 0; i < numberIntegers; i++) {
990	int iColumn = integerVariable[i];
991	originalLower_[i] = lower[iColumn];
992	originalUpper_[i] = upper[iColumn];
993	if (upper[iColumn] - lower[iColumn] > 1.5)
994	all01 = false;
995	else if (upper[iColumn] - lower[iColumn] == 1.0)
996	number01++;
997	}
998	if (all01 && !typeCuts_)
999	typeCuts_ = 1; // may as well so we don't have to deal with refine
1000	if (!number01 && !typeCuts_) {
1001	if (model_->messageHandler()->logLevel() > 1)
1002	printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
1003	typeCuts_ = -1;
1004	} else {
1005	if (model_->messageHandler()->logLevel() > 1) {
1006	std::string type;
1007	if (all01) {
1008	printf("%d 0-1 variables normal local cuts\n",
1009	number01);
1010	} else if (typeCuts_) {
1011	printf("%d 0-1 variables, %d other - general integer local cuts\n",
1012	number01, numberIntegers - number01);
1013	} else {
1014	printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
1015	number01, numberIntegers - number01);
1016	}
1017	printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
1018	maxDiversification_, range_, timeLimit_, nodeLimit_);
1019	}
1020	}
1021	int numberColumns = model_->getNumCols();
1022	savedSolution_ = new double [numberColumns];
1023	memset(savedSolution_, 0, numberColumns*sizeof(double));
1024	if (solution) {
1025	rhs_ = range_;
1026	// Check feasible
1027	int goodSolution = createCut(solution, cut_);
1028	if (goodSolution >= 0) {
1029	for (i = 0; i < numberIntegers; i++) {
1030	int iColumn = integerVariable[i];
1031	double value = floor(solution[iColumn] + 0.5);
1032	// fix so setBestSolution will work
1033	solver->setColLower(iColumn, value);
1034	solver->setColUpper(iColumn, value);
1035	}
1036	model_->reserveCurrentSolution();
1037	// Create cut and get total gap
1038	if (newSolutionValue < bestCutoff_) {
1039	model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
1040	bestCutoff_ = model_->getCutoff();
1041	// save as best solution
1042	memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1043	}
1044	for (i = 0; i < numberIntegers; i++) {
1045	int iColumn = integerVariable[i];
1046	// restore bounds
1047	solver->setColLower(iColumn, originalLower_[i]);
1048	solver->setColUpper(iColumn, originalUpper_[i]);
1049	}
1050	// make sure can't stop on gap
1051	model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1052	} else {
1053	model_ = NULL;
1054	}
1055	} else {
1056	// no solution
1057	rhs_ = 1.0e50;
1058	// make sure can't stop on gap
1059	model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1060	}
1061	}
1062	CbcTreeVariable::~CbcTreeVariable()
1063	{
1064	delete [] originalLower_;
1065	delete [] originalUpper_;
1066	delete [] bestSolution_;
1067	delete [] savedSolution_;
1068	delete localNode_;
1069	}
1070	// Copy constructor
1071	CbcTreeVariable::CbcTreeVariable ( const CbcTreeVariable & rhs)
1072	: CbcTree(rhs),
1073	saveNumberSolutions_(rhs.saveNumberSolutions_),
1074	model_(rhs.model_),
1075	range_(rhs.range_),
1076	typeCuts_(rhs.typeCuts_),
1077	maxDiversification_(rhs.maxDiversification_),
1078	diversification_(rhs.diversification_),
1079	nextStrong_(rhs.nextStrong_),
1080	rhs_(rhs.rhs_),
1081	savedGap_(rhs.savedGap_),
1082	bestCutoff_(rhs.bestCutoff_),
1083	timeLimit_(rhs.timeLimit_),
1084	startTime_(rhs.startTime_),
1085	nodeLimit_(rhs.nodeLimit_),
1086	startNode_(rhs.startNode_),
1087	searchType_(rhs.searchType_),
1088	refine_(rhs.refine_)
1089	{
1090	cut_ = rhs.cut_;
1091	fixedCut_ = rhs.fixedCut_;
1092	if (rhs.localNode_)
1093	localNode_ = new CbcNode(*rhs.localNode_);
1094	else
1095	localNode_ = NULL;
1096	if (rhs.originalLower_) {
1097	int numberIntegers = model_->numberIntegers();
1098	originalLower_ = new double [numberIntegers];
1099	memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
1100	originalUpper_ = new double [numberIntegers];
1101	memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
1102	} else {
1103	originalLower_ = NULL;
1104	originalUpper_ = NULL;
1105	}
1106	if (rhs.bestSolution_) {
1107	int numberColumns = model_->getNumCols();
1108	bestSolution_ = new double [numberColumns];
1109	memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
1110	} else {
1111	bestSolution_ = NULL;
1112	}
1113	if (rhs.savedSolution_) {
1114	int numberColumns = model_->getNumCols();
1115	savedSolution_ = new double [numberColumns];
1116	memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
1117	} else {
1118	savedSolution_ = NULL;
1119	}
1120	}
1121	//----------------------------------------------------------------
1122	// Assignment operator
1123	//-------------------------------------------------------------------
1124	CbcTreeVariable &
1125	CbcTreeVariable::operator=(const CbcTreeVariable & rhs)
1126	{
1127	if (this != &rhs) {
1128	CbcTree::operator=(rhs);
1129	saveNumberSolutions_ = rhs.saveNumberSolutions_;
1130	cut_ = rhs.cut_;
1131	fixedCut_ = rhs.fixedCut_;
1132	delete localNode_;
1133	if (rhs.localNode_)
1134	localNode_ = new CbcNode(*rhs.localNode_);
1135	else
1136	localNode_ = NULL;
1137	model_ = rhs.model_;
1138	range_ = rhs.range_;
1139	typeCuts_ = rhs.typeCuts_;
1140	maxDiversification_ = rhs.maxDiversification_;
1141	diversification_ = rhs.diversification_;
1142	nextStrong_ = rhs.nextStrong_;
1143	rhs_ = rhs.rhs_;
1144	savedGap_ = rhs.savedGap_;
1145	bestCutoff_ = rhs.bestCutoff_;
1146	timeLimit_ = rhs.timeLimit_;
1147	startTime_ = rhs.startTime_;
1148	nodeLimit_ = rhs.nodeLimit_;
1149	startNode_ = rhs.startNode_;
1150	searchType_ = rhs.searchType_;
1151	refine_ = rhs.refine_;
1152	delete [] originalLower_;
1153	delete [] originalUpper_;
1154	if (rhs.originalLower_) {
1155	int numberIntegers = model_->numberIntegers();
1156	originalLower_ = new double [numberIntegers];
1157	memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
1158	originalUpper_ = new double [numberIntegers];
1159	memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
1160	} else {
1161	originalLower_ = NULL;
1162	originalUpper_ = NULL;
1163	}
1164	delete [] bestSolution_;
1165	if (rhs.bestSolution_) {
1166	int numberColumns = model_->getNumCols();
1167	bestSolution_ = new double [numberColumns];
1168	memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
1169	} else {
1170	bestSolution_ = NULL;
1171	}
1172	delete [] savedSolution_;
1173	if (rhs.savedSolution_) {
1174	int numberColumns = model_->getNumCols();
1175	savedSolution_ = new double [numberColumns];
1176	memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
1177	} else {
1178	savedSolution_ = NULL;
1179	}
1180	}
1181	return *this;
1182	}
1183	// Clone
1184	CbcTree *
1185	CbcTreeVariable::clone() const
1186	{
1187	return new CbcTreeVariable(*this);
1188	}
1189	// Pass in solution (so can be used after heuristic)
1190	void
1191	CbcTreeVariable::passInSolution(const double * solution, double solutionValue)
1192	{
1193	int numberColumns = model_->getNumCols();
1194	delete [] savedSolution_;
1195	savedSolution_ = new double [numberColumns];
1196	memcpy(savedSolution_, solution, numberColumns*sizeof(double));
1197	rhs_ = range_;
1198	// Check feasible
1199	int goodSolution = createCut(solution, cut_);
1200	if (goodSolution >= 0) {
1201	bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
1202	} else {
1203	model_ = NULL;
1204	}
1205	}
1206	// Return the top node of the heap
1207	CbcNode *
1208	CbcTreeVariable::top() const
1209	{
1210	#ifdef CBC_DEBUG
1211	int smallest = 9999999;
1212	int largest = -1;
1213	double smallestD = 1.0e30;
1214	double largestD = -1.0e30;
1215	int n = nodes_.size();
1216	for (int i = 0; i < n; i++) {
1217	int nn = nodes_[i]->nodeInfo()->nodeNumber();
1218	double dd = nodes_[i]->objectiveValue();
1219	largest = CoinMax(largest, nn);
1220	smallest = CoinMin(smallest, nn);
1221	largestD = CoinMax(largestD, dd);
1222	smallestD = CoinMin(smallestD, dd);
1223	}
1224	if (model_->messageHandler()->logLevel() > 1) {
1225	printf("smallest %d, largest %d, top %d\n", smallest, largest,
1226	nodes_.front()->nodeInfo()->nodeNumber());
1227	printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
1228	}
1229	#endif
1230	return nodes_.front();
1231	}
1232
1233	// Add a node to the heap
1234	void
1235	CbcTreeVariable::push(CbcNode * x)
1236	{
1237	if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
1238	startNode_ = model_->getNodeCount();
1239	// save copy of node
1240	localNode_ = new CbcNode(*x);
1241
1242	if (cut_.row().getNumElements()) {
1243	// Add to global cuts
1244	// we came in with solution
1245	model_->makeGlobalCut(cut_);
1246	if (model_->messageHandler()->logLevel() > 1)
1247	printf("initial cut - rhs %g %g\n",
1248	cut_.lb(), cut_.ub());
1249	searchType_ = 1;
1250	} else {
1251	// stop on first solution
1252	searchType_ = 0;
1253	}
1254	startTime_ = static_cast<int> (CoinCpuTime());
1255	saveNumberSolutions_ = model_->getSolutionCount();
1256	}
1257	nodes_.push_back(x);
1258	#ifdef CBC_DEBUG
1259	if (model_->messageHandler()->logLevel() > 0)
1260	printf("pushing node onto heap %d %x %x\n",
1261	x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
1262	#endif
1263	std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
1264	}
1265
1266	// Remove the top node from the heap
1267	void
1268	CbcTreeVariable::pop()
1269	{
1270	std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
1271	nodes_.pop_back();
1272	}
1273	// Test if empty - does work if so
1274	bool
1275	CbcTreeVariable::empty()
1276	{
1277	if (typeCuts_ < 0)
1278	return !nodes_.size();
1279	/* state -
1280	0 iterating
1281	1 subtree finished optimal solution for subtree found
1282	2 subtree finished and no solution found
1283	3 subtree exiting and solution found
1284	4 subtree exiting and no solution found
1285	*/
1286	int state = 0;
1287	assert (searchType_ != 2);
1288	if (searchType_) {
1289	if (CoinCpuTime() - startTime_ > timeLimit_ \|\| model_->getNodeCount() - startNode_ >= nodeLimit_) {
1290	state = 4;
1291	}
1292	} else {
1293	if (model_->getSolutionCount() > saveNumberSolutions_) {
1294	state = 4;
1295	}
1296	}
1297	if (!nodes_.size())
1298	state = 2;
1299	if (!state) {
1300	return false;
1301	}
1302	// Finished this phase
1303	int numberColumns = model_->getNumCols();
1304	if (model_->getSolutionCount() > saveNumberSolutions_) {
1305	if (model_->getCutoff() < bestCutoff_) {
1306	// Save solution
1307	if (!bestSolution_)
1308	bestSolution_ = new double [numberColumns];
1309	memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1310	bestCutoff_ = model_->getCutoff();
1311	}
1312	state--;
1313	}
1314	// get rid of all nodes (safe even if already done)
1315	double bestPossibleObjective;
1316	cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
1317
1318	double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
1319	if (model_->messageHandler()->logLevel() > 1)
1320	printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
1321	state,
1322	model_->getNodeCount() - startNode_,
1323	static_cast<int> (CoinCpuTime()) - startTime_,
1324	model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
1325	saveNumberSolutions_ = model_->getSolutionCount();
1326	bool finished = false;
1327	bool lastTry = false;
1328	switch (state) {
1329	case 1:
1330	// solution found and subtree exhausted
1331	if (rhs_ > 1.0e30) {
1332	finished = true;
1333	} else {
1334	// find global cut and reverse
1335	reverseCut(1);
1336	searchType_ = 1; // first false
1337	rhs_ = range_; // reset range
1338	nextStrong_ = false;
1339
1340	// save best solution in this subtree
1341	memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1342	}
1343	break;
1344	case 2:
1345	// solution not found and subtree exhausted
1346	if (rhs_ > 1.0e30) {
1347	finished = true;
1348	} else {
1349	// find global cut and reverse
1350	reverseCut(2);
1351	searchType_ = 1; // first false
1352	if (diversification_ < maxDiversification_) {
1353	if (nextStrong_) {
1354	diversification_++;
1355	// cut is valid so don't model_->setCutoff(1.0e50);
1356	searchType_ = 0;
1357	}
1358	nextStrong_ = true;
1359	rhs_ += range_ / 2;
1360	} else {
1361	// This will be last try (may hit max time)
1362	lastTry = true;
1363	if (!maxDiversification_)
1364	typeCuts_ = -1; // make sure can't start again
1365	model_->setCutoff(bestCutoff_);
1366	if (model_->messageHandler()->logLevel() > 1)
1367	printf("Exiting local search with current set of cuts\n");
1368	rhs_ = 1.0e100;
1369	// Can now stop on gap
1370	model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1371	}
1372	}
1373	break;
1374	case 3:
1375	// solution found and subtree not exhausted
1376	if (rhs_ < 1.0e30) {
1377	if (searchType_) {
1378	if (!typeCuts_ && refine_ && searchType_ == 1) {
1379	// We need to check we have best solution given these 0-1 values
1380	OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
1381	CbcModel * subModel = model_->subTreeModel(subSolver);
1382	CbcTree normalTree;
1383	subModel->passInTreeHandler(normalTree);
1384	int numberIntegers = model_->numberIntegers();
1385	const int * integerVariable = model_->integerVariable();
1386	const double * solution = model_->bestSolution();
1387	int i;
1388	int numberColumns = model_->getNumCols();
1389	for (i = 0; i < numberIntegers; i++) {
1390	int iColumn = integerVariable[i];
1391	double value = floor(solution[iColumn] + 0.5);
1392	if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1393	continue; // skip as not 0-1
1394	if (originalLower_[i] == originalUpper_[i])
1395	continue;
1396	subSolver->setColLower(iColumn, value);
1397	subSolver->setColUpper(iColumn, value);
1398	}
1399	subSolver->initialSolve();
1400	// We can copy cutoff
1401	// But adjust
1402	subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
1403	subModel->setSolutionCount(0);
1404	assert (subModel->isProvenOptimal());
1405	if (!subModel->typePresolve()) {
1406	subModel->branchAndBound();
1407	if (subModel->status()) {
1408	model_->incrementSubTreeStopped();
1409	}
1410	//printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
1411	// subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
1412	double newCutoff = subModel->getMinimizationObjValue() -
1413	subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
1414	if (subModel->getSolutionCount()) {
1415	if (!subModel->status())
1416	assert (subModel->isProvenOptimal());
1417	memcpy(model_->bestSolution(), subModel->bestSolution(),
1418	numberColumns*sizeof(double));
1419	model_->setCutoff(newCutoff);
1420	}
1421	} else if (subModel->typePresolve() == 1) {
1422	CbcModel * model2 = subModel->integerPresolve(true);
1423	if (model2) {
1424	// Do complete search
1425	model2->branchAndBound();
1426	// get back solution
1427	subModel->originalModel(model2, false);
1428	if (model2->status()) {
1429	model_->incrementSubTreeStopped();
1430	}
1431	double newCutoff = model2->getMinimizationObjValue() -
1432	model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
1433	if (model2->getSolutionCount()) {
1434	if (!model2->status())
1435	assert (model2->isProvenOptimal());
1436	memcpy(model_->bestSolution(), subModel->bestSolution(),
1437	numberColumns*sizeof(double));
1438	model_->setCutoff(newCutoff);
1439	}
1440	delete model2;
1441	} else {
1442	// infeasible - could just be - due to cutoff
1443	}
1444	} else {
1445	// too dangerous at present
1446	assert (subModel->typePresolve() != 2);
1447	}
1448	if (model_->getCutoff() < bestCutoff_) {
1449	// Save solution
1450	if (!bestSolution_)
1451	bestSolution_ = new double [numberColumns];
1452	memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1453	bestCutoff_ = model_->getCutoff();
1454	}
1455	delete subModel;
1456	}
1457	// we have done search to make sure best general solution
1458	searchType_ = 1;
1459	// Reverse cut weakly
1460	reverseCut(3, rhs_);
1461	} else {
1462	searchType_ = 1;
1463	// delete last cut
1464	deleteCut(cut_);
1465	}
1466	} else {
1467	searchType_ = 1;
1468	}
1469	// save best solution in this subtree
1470	memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1471	nextStrong_ = false;
1472	rhs_ = range_;
1473	break;
1474	case 4:
1475	// solution not found and subtree not exhausted
1476	if (maxDiversification_) {
1477	if (nextStrong_) {
1478	// Reverse cut weakly
1479	reverseCut(4, rhs_);
1480	model_->setCutoff(1.0e50);
1481	diversification_++;
1482	searchType_ = 0;
1483	} else {
1484	// delete last cut
1485	deleteCut(cut_);
1486	searchType_ = 1;
1487	}
1488	nextStrong_ = true;
1489	rhs_ += range_ / 2;
1490	} else {
1491	// special case when using as heuristic
1492	// Reverse cut weakly if lb -infinity
1493	reverseCut(4, rhs_);
1494	// This will be last try (may hit max time0
1495	lastTry = true;
1496	model_->setCutoff(bestCutoff_);
1497	if (model_->messageHandler()->logLevel() > 1)
1498	printf("Exiting local search with current set of cuts\n");
1499	rhs_ = 1.0e100;
1500	// Can now stop on gap
1501	model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1502	typeCuts_ = -1;
1503	}
1504	break;
1505	}
1506	if (rhs_ < 1.0e30 \|\| lastTry) {
1507	int goodSolution = createCut(savedSolution_, cut_);
1508	if (goodSolution >= 0) {
1509	// Add to global cuts
1510	model_->makeGlobalCut(cut_);
1511	CbcRowCuts * global = model_->globalCuts();
1512	int n = global->sizeRowCuts();
1513	OsiRowCut * rowCut = global->rowCutPtr(n - 1);
1514	if (model_->messageHandler()->logLevel() > 1)
1515	printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
1516	n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
1517	const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1518	if (debugger) {
1519	if (debugger->invalidCut(*rowCut))
1520	printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1521	}
1522	for (int i = 0; i < n; i++) {
1523	rowCut = global->rowCutPtr(i);
1524	if (model_->messageHandler()->logLevel() > 1)
1525	printf("%d - rhs %g %g\n",
1526	i, rowCut->lb(), rowCut->ub());
1527	}
1528	}
1529	// put back node
1530	startTime_ = static_cast<int> (CoinCpuTime());
1531	startNode_ = model_->getNodeCount();
1532	if (localNode_) {
1533	// save copy of node
1534	CbcNode * localNode2 = new CbcNode(*localNode_);
1535	// But localNode2 now owns cuts so swap
1536	//printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
1537	// localNode_,localNode_->nodeInfo());
1538	nodes_.push_back(localNode_);
1539	localNode_ = localNode2;
1540	std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
1541	}
1542	}
1543	return finished;
1544	}
1545	// We may have got an intelligent tree so give it one more chance
1546	void
1547	CbcTreeVariable::endSearch()
1548	{
1549	if (typeCuts_ >= 0) {
1550	// copy best solution to model
1551	int numberColumns = model_->getNumCols();
1552	if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
1553	memcpy(model_->bestSolution(), bestSolution_, numberColumns*sizeof(double));
1554	model_->setCutoff(bestCutoff_);
1555	// recompute objective value
1556	const double * objCoef = model_->getObjCoefficients();
1557	double objOffset = 0.0;
1558	model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
1559
1560	// Compute dot product of objCoef and colSol and then adjust by offset
1561	double objValue = -objOffset;
1562	for ( int i = 0 ; i < numberColumns ; i++ )
1563	objValue += objCoef[i] * bestSolution_[i];
1564	model_->setMinimizationObjValue(objValue);
1565	}
1566	// Can now stop on gap
1567	model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1568	}
1569	}
1570	// Create cut
1571	int
1572	CbcTreeVariable::createCut(const double * solution, OsiRowCut & rowCut)
1573	{
1574	if (rhs_ > 1.0e20)
1575	return -1;
1576	OsiSolverInterface * solver = model_->solver();
1577	const double * rowLower = solver->getRowLower();
1578	const double * rowUpper = solver->getRowUpper();
1579	//const double * solution = solver->getColSolution();
1580	//const double * objective = solver->getObjCoefficients();
1581	double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1582	double primalTolerance;
1583	solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1584	// relax
1585	primalTolerance *= 1000.0;
1586
1587	int numberRows = model_->getNumRows();
1588
1589	int numberIntegers = model_->numberIntegers();
1590	const int * integerVariable = model_->integerVariable();
1591	int i;
1592
1593	// Check feasible
1594
1595	double * rowActivity = new double[numberRows];
1596	memset(rowActivity, 0, numberRows*sizeof(double));
1597	solver->getMatrixByCol()->times(solution, rowActivity) ;
1598	int goodSolution = 0;
1599	// check was feasible
1600	for (i = 0; i < numberRows; i++) {
1601	if (rowActivity[i] < rowLower[i] - primalTolerance) {
1602	goodSolution = -1;
1603	} else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
1604	goodSolution = -1;
1605	}
1606	}
1607	delete [] rowActivity;
1608	for (i = 0; i < numberIntegers; i++) {
1609	int iColumn = integerVariable[i];
1610	double value = solution[iColumn];
1611	if (fabs(floor(value + 0.5) - value) > integerTolerance) {
1612	goodSolution = -1;
1613	}
1614	}
1615	// zap cut
1616	if (goodSolution == 0) {
1617	// Create cut and get total gap
1618	CoinPackedVector cut;
1619	double rhs = rhs_;
1620	double maxValue = 0.0;
1621	for (i = 0; i < numberIntegers; i++) {
1622	int iColumn = integerVariable[i];
1623	double value = floor(solution[iColumn] + 0.5);
1624	if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1625	continue; // skip as not 0-1
1626	if (originalLower_[i] == originalUpper_[i])
1627	continue;
1628	double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
1629	if (value == originalLower_[i]) {
1630	rhs += mu * originalLower_[i];
1631	cut.insert(iColumn, 1.0);
1632	maxValue += originalUpper_[i];
1633	} else if (value == originalUpper_[i]) {
1634	rhs -= mu * originalUpper_[i];
1635	cut.insert(iColumn, -1.0);
1636	maxValue += originalLower_[i];
1637	}
1638	}
1639	if (maxValue < rhs - primalTolerance) {
1640	if (model_->messageHandler()->logLevel() > 1)
1641	printf("slack cut\n");
1642	goodSolution = 1;
1643	}
1644	rowCut.setRow(cut);
1645	rowCut.setLb(-COIN_DBL_MAX);
1646	rowCut.setUb(rhs);
1647	rowCut.setGloballyValid();
1648	if (model_->messageHandler()->logLevel() > 1)
1649	printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
1650	#ifdef CBC_DEBUG
1651	if (model_->messageHandler()->logLevel() > 0) {
1652	int k;
1653	for (k = 0; k < cut.getNumElements(); k++) {
1654	printf("%i %g ", cut.getIndices()[k], cut.getElements()[k]);
1655	if ((k + 1) % 5 == 0)
1656	printf("\n");
1657	}
1658	if (k % 5 != 0)
1659	printf("\n");
1660	}
1661	#endif
1662	return goodSolution;
1663	} else {
1664	if (model_->messageHandler()->logLevel() > 1)
1665	printf("Not a good solution\n");
1666	return -1;
1667	}
1668	}
1669	// Other side of last cut branch
1670	void
1671	CbcTreeVariable::reverseCut(int state, double bias)
1672	{
1673	// find global cut
1674	CbcRowCuts * global = model_->globalCuts();
1675	int n = global->sizeRowCuts();
1676	int i;
1677	OsiRowCut * rowCut = NULL;
1678	for ( i = 0; i < n; i++) {
1679	rowCut = global->rowCutPtr(i);
1680	if (cut_ == *rowCut) {
1681	break;
1682	}
1683	}
1684	if (!rowCut) {
1685	// must have got here in odd way e.g. strong branching
1686	return;
1687	}
1688	if (rowCut->lb() > -1.0e10)
1689	return;
1690	// get smallest element
1691	double smallest = COIN_DBL_MAX;
1692	CoinPackedVector row = cut_.row();
1693	for (int k = 0; k < row.getNumElements(); k++)
1694	smallest = CoinMin(smallest, fabs(row.getElements()[k]));
1695	if (!typeCuts_ && !refine_) {
1696	// Reverse cut very very weakly
1697	if (state > 2)
1698	smallest = 0.0;
1699	}
1700	// replace by other way
1701	if (model_->messageHandler()->logLevel() > 1)
1702	printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
1703	i, n, rowCut->lb(), rowCut->ub());
1704	rowCut->setLb(rowCut->ub() + smallest - bias);
1705	rowCut->setUb(COIN_DBL_MAX);
1706	if (model_->messageHandler()->logLevel() > 1)
1707	printf("new rhs %g %g, bias %g smallest %g ",
1708	rowCut->lb(), rowCut->ub(), bias, smallest);
1709	const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1710	if (debugger) {
1711	if (debugger->invalidCut(*rowCut))
1712	printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1713	}
1714	}
1715	// Delete last cut branch
1716	void
1717	CbcTreeVariable::deleteCut(OsiRowCut & cut)
1718	{
1719	// find global cut
1720	CbcRowCuts * global = model_->globalCuts();
1721	int n = global->sizeRowCuts();
1722	int i;
1723	OsiRowCut * rowCut = NULL;
1724	for ( i = 0; i < n; i++) {
1725	rowCut = global->rowCutPtr(i);
1726	if (cut == *rowCut) {
1727	break;
1728	}
1729	}
1730	assert (i < n);
1731	// delete last cut
1732	if (model_->messageHandler()->logLevel() > 1)
1733	printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
1734	i, n, rowCut->lb(), rowCut->ub());
1735	global->eraseRowCut(i);
1736	}
1737	// Create C++ lines to get to current state
1738	void
1739	CbcTreeVariable::generateCpp( FILE * fp)
1740	{
1741	CbcTreeVariable other;
1742	fprintf(fp, "0#include \"CbcTreeVariable.hpp\"\n");
1743	fprintf(fp, "5 CbcTreeVariable variableTree(cbcModel,NULL);\n");
1744	if (range_ != other.range_)
1745	fprintf(fp, "5 variableTree.setRange(%d);\n", range_);
1746	if (typeCuts_ != other.typeCuts_)
1747	fprintf(fp, "5 variableTree.setTypeCuts(%d);\n", typeCuts_);
1748	if (maxDiversification_ != other.maxDiversification_)
1749	fprintf(fp, "5 variableTree.setMaxDiversification(%d);\n", maxDiversification_);
1750	if (timeLimit_ != other.timeLimit_)
1751	fprintf(fp, "5 variableTree.setTimeLimit(%d);\n", timeLimit_);
1752	if (nodeLimit_ != other.nodeLimit_)
1753	fprintf(fp, "5 variableTree.setNodeLimit(%d);\n", nodeLimit_);
1754	if (refine_ != other.refine_)
1755	fprintf(fp, "5 variableTree.setRefine(%s);\n", refine_ ? "true" : "false");
1756	fprintf(fp, "5 cbcModel->passInTreeHandler(variableTree);\n");
1757	}
1758
1759
1760

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

Download in other formats: