Context Navigation

source: trunk/Clp/src/ClpNonLinearCost.cpp @ 1769

Last change on this file since 1769 was 1769, checked in by forrest, 10 years ago
changes for advanced use of Clp
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 82.2 KB

Line
1	/* $Id: ClpNonLinearCost.cpp 1769 2011-07-26 09:31:51Z forrest $ */
2	// Copyright (C) 2002, 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 "CoinPragma.hpp"
7	#include <iostream>
8	#include <cassert>
9
10	#include "CoinIndexedVector.hpp"
11
12	#include "ClpSimplex.hpp"
13	#include "CoinHelperFunctions.hpp"
14	#include "ClpNonLinearCost.hpp"
15	//#############################################################################
16	// Constructors / Destructor / Assignment
17	//#############################################################################
18
19	//-------------------------------------------------------------------
20	// Default Constructor
21	//-------------------------------------------------------------------
22	ClpNonLinearCost::ClpNonLinearCost () :
23	changeCost_(0.0),
24	feasibleCost_(0.0),
25	infeasibilityWeight_(-1.0),
26	largestInfeasibility_(0.0),
27	sumInfeasibilities_(0.0),
28	averageTheta_(0.0),
29	numberRows_(0),
30	numberColumns_(0),
31	start_(NULL),
32	whichRange_(NULL),
33	offset_(NULL),
34	lower_(NULL),
35	cost_(NULL),
36	model_(NULL),
37	infeasible_(NULL),
38	numberInfeasibilities_(-1),
39	status_(NULL),
40	bound_(NULL),
41	cost2_(NULL),
42	method_(1),
43	convex_(true),
44	bothWays_(false)
45	{
46
47	}
48	//#define VALIDATE
49	#ifdef VALIDATE
50	static double * saveLowerV = NULL;
51	static double * saveUpperV = NULL;
52	#ifdef NDEBUG
53	Validate sgould not be set if no debug
54	#endif
55	#endif
56	/* Constructor from simplex.
57	This will just set up wasteful arrays for linear, but
58	later may do dual analysis and even finding duplicate columns
59	*/
60	ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model, int method)
61	{
62	method = 2;
63	model_ = model;
64	numberRows_ = model_->numberRows();
65	//if (numberRows_==402) {
66	//model_->setLogLevel(63);
67	//model_->setMaximumIterations(30000);
68	//}
69	numberColumns_ = model_->numberColumns();
70	// If gub then we need this extra
71	int numberExtra = model_->numberExtraRows();
72	if (numberExtra)
73	method = 1;
74	int numberTotal1 = numberRows_ + numberColumns_;
75	int numberTotal = numberTotal1 + numberExtra;
76	convex_ = true;
77	bothWays_ = false;
78	method_ = method;
79	numberInfeasibilities_ = 0;
80	changeCost_ = 0.0;
81	feasibleCost_ = 0.0;
82	infeasibilityWeight_ = -1.0;
83	double * cost = model_->costRegion();
84	// check if all 0
85	int iSequence;
86	bool allZero = true;
87	for (iSequence = 0; iSequence < numberTotal1; iSequence++) {
88	if (cost[iSequence]) {
89	allZero = false;
90	break;
91	}
92	}
93	if (allZero&&model_->clpMatrix()->type()<15)
94	model_->setInfeasibilityCost(1.0);
95	double infeasibilityCost = model_->infeasibilityCost();
96	sumInfeasibilities_ = 0.0;
97	averageTheta_ = 0.0;
98	largestInfeasibility_ = 0.0;
99	// All arrays NULL to start
100	status_ = NULL;
101	bound_ = NULL;
102	cost2_ = NULL;
103	start_ = NULL;
104	whichRange_ = NULL;
105	offset_ = NULL;
106	lower_ = NULL;
107	cost_ = NULL;
108	infeasible_ = NULL;
109
110	double * upper = model_->upperRegion();
111	double * lower = model_->lowerRegion();
112
113	// See how we are storing things
114	bool always4 = (model_->clpMatrix()->
115	generalExpanded(model_, 10, iSequence) != 0);
116	if (always4)
117	method_ = 1;
118	if (CLP_METHOD1) {
119	start_ = new int [numberTotal+1];
120	whichRange_ = new int [numberTotal];
121	offset_ = new int [numberTotal];
122	memset(offset_, 0, numberTotal * sizeof(int));
123
124
125	// First see how much space we need
126	int put = 0;
127
128	// For quadratic we need -inf,0,0,+inf
129	for (iSequence = 0; iSequence < numberTotal1; iSequence++) {
130	if (!always4) {
131	if (lower[iSequence] > -COIN_DBL_MAX)
132	put++;
133	if (upper[iSequence] < COIN_DBL_MAX)
134	put++;
135	put += 2;
136	} else {
137	put += 4;
138	}
139	}
140
141	// and for extra
142	put += 4 * numberExtra;
143	#ifndef NDEBUG
144	int kPut = put;
145	#endif
146	lower_ = new double [put];
147	cost_ = new double [put];
148	infeasible_ = new unsigned int[(put+31)>>5];
149	memset(infeasible_, 0, ((put + 31) >> 5)*sizeof(unsigned int));
150
151	put = 0;
152
153	start_[0] = 0;
154
155	for (iSequence = 0; iSequence < numberTotal1; iSequence++) {
156	if (!always4) {
157	if (lower[iSequence] > -COIN_DBL_MAX) {
158	lower_[put] = -COIN_DBL_MAX;
159	setInfeasible(put, true);
160	cost_[put++] = cost[iSequence] - infeasibilityCost;
161	}
162	whichRange_[iSequence] = put;
163	lower_[put] = lower[iSequence];
164	cost_[put++] = cost[iSequence];
165	lower_[put] = upper[iSequence];
166	cost_[put++] = cost[iSequence] + infeasibilityCost;
167	if (upper[iSequence] < COIN_DBL_MAX) {
168	lower_[put] = COIN_DBL_MAX;
169	setInfeasible(put - 1, true);
170	cost_[put++] = 1.0e50;
171	}
172	} else {
173	lower_[put] = -COIN_DBL_MAX;
174	setInfeasible(put, true);
175	cost_[put++] = cost[iSequence] - infeasibilityCost;
176	whichRange_[iSequence] = put;
177	lower_[put] = lower[iSequence];
178	cost_[put++] = cost[iSequence];
179	lower_[put] = upper[iSequence];
180	cost_[put++] = cost[iSequence] + infeasibilityCost;
181	lower_[put] = COIN_DBL_MAX;
182	setInfeasible(put - 1, true);
183	cost_[put++] = 1.0e50;
184	}
185	start_[iSequence+1] = put;
186	}
187	for (; iSequence < numberTotal; iSequence++) {
188
189	lower_[put] = -COIN_DBL_MAX;
190	setInfeasible(put, true);
191	put++;
192	whichRange_[iSequence] = put;
193	lower_[put] = 0.0;
194	cost_[put++] = 0.0;
195	lower_[put] = 0.0;
196	cost_[put++] = 0.0;
197	lower_[put] = COIN_DBL_MAX;
198	setInfeasible(put - 1, true);
199	cost_[put++] = 1.0e50;
200	start_[iSequence+1] = put;
201	}
202	assert (put <= kPut);
203	}
204	#ifdef FAST_CLPNON
205	// See how we are storing things
206	CoinAssert (model_->clpMatrix()->
207	generalExpanded(model_, 10, iSequence) == 0);
208	#endif
209	if (CLP_METHOD2) {
210	assert (!numberExtra);
211	bound_ = new double[numberTotal];
212	cost2_ = new double[numberTotal];
213	status_ = new unsigned char[numberTotal];
214	#ifdef VALIDATE
215	delete [] saveLowerV;
216	saveLowerV = CoinCopyOfArray(model_->lowerRegion(), numberTotal);
217	delete [] saveUpperV;
218	saveUpperV = CoinCopyOfArray(model_->upperRegion(), numberTotal);
219	#endif
220	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
221	bound_[iSequence] = 0.0;
222	cost2_[iSequence] = cost[iSequence];
223	setInitialStatus(status_[iSequence]);
224	}
225	}
226	}
227	// Refresh - assuming regions OK
228	void
229	ClpNonLinearCost::refresh()
230	{
231	int numberTotal = numberRows_ + numberColumns_;
232	numberInfeasibilities_ = 0;
233	sumInfeasibilities_ = 0.0;
234	largestInfeasibility_ = 0.0;
235	double infeasibilityCost = model_->infeasibilityCost();
236	double primalTolerance = model_->currentPrimalTolerance();
237	double * cost = model_->costRegion();
238	double * upper = model_->upperRegion();
239	double * lower = model_->lowerRegion();
240	double * solution = model_->solutionRegion();
241	for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
242	cost2_[iSequence] = cost[iSequence];
243	double value = solution[iSequence];
244	double lowerValue = lower[iSequence];
245	double upperValue = upper[iSequence];
246	if (value - upperValue <= primalTolerance) {
247	if (value - lowerValue >= -primalTolerance) {
248	// feasible
249	status_[iSequence] = static_cast<unsigned char>(CLP_FEASIBLE \| (CLP_SAME << 4));
250	bound_[iSequence] = 0.0;
251	} else {
252	// below
253	double infeasibility = lowerValue - value - primalTolerance;
254	sumInfeasibilities_ += infeasibility;
255	largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
256	cost[iSequence] -= infeasibilityCost;
257	numberInfeasibilities_++;
258	status_[iSequence] = static_cast<unsigned char>(CLP_BELOW_LOWER \| (CLP_SAME << 4));
259	bound_[iSequence] = upperValue;
260	upper[iSequence] = lowerValue;
261	lower[iSequence] = -COIN_DBL_MAX;
262	}
263	} else {
264	// above
265	double infeasibility = value - upperValue - primalTolerance;
266	sumInfeasibilities_ += infeasibility;
267	largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
268	cost[iSequence] += infeasibilityCost;
269	numberInfeasibilities_++;
270	status_[iSequence] = static_cast<unsigned char>(CLP_ABOVE_UPPER \| (CLP_SAME << 4));
271	bound_[iSequence] = lowerValue;
272	lower[iSequence] = upperValue;
273	upper[iSequence] = COIN_DBL_MAX;
274	}
275	}
276	// checkInfeasibilities(model_->primalTolerance());
277
278	}
279	// Refreshes costs always makes row costs zero
280	void
281	ClpNonLinearCost::refreshCosts(const double * columnCosts)
282	{
283	double * cost = model_->costRegion();
284	// zero row costs
285	memset(cost + numberColumns_, 0, numberRows_ * sizeof(double));
286	// copy column costs
287	CoinMemcpyN(columnCosts, numberColumns_, cost);
288	if ((method_ & 1) != 0) {
289	for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
290	int start = start_[iSequence];
291	int end = start_[iSequence+1] - 1;
292	double thisFeasibleCost = cost[iSequence];
293	if (infeasible(start)) {
294	cost_[start] = thisFeasibleCost - infeasibilityWeight_;
295	cost_[start+1] = thisFeasibleCost;
296	} else {
297	cost_[start] = thisFeasibleCost;
298	}
299	if (infeasible(end - 1)) {
300	cost_[end-1] = thisFeasibleCost + infeasibilityWeight_;
301	}
302	}
303	}
304	if (CLP_METHOD2) {
305	for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
306	cost2_[iSequence] = cost[iSequence];
307	}
308	}
309	}
310	ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model, const int * starts,
311	const double * lowerNon, const double * costNon)
312	{
313	#ifndef FAST_CLPNON
314	// what about scaling? - only try without it initially
315	assert(!model->scalingFlag());
316	model_ = model;
317	numberRows_ = model_->numberRows();
318	numberColumns_ = model_->numberColumns();
319	int numberTotal = numberRows_ + numberColumns_;
320	convex_ = true;
321	bothWays_ = true;
322	start_ = new int [numberTotal+1];
323	whichRange_ = new int [numberTotal];
324	offset_ = new int [numberTotal];
325	memset(offset_, 0, numberTotal * sizeof(int));
326
327	double whichWay = model_->optimizationDirection();
328	COIN_DETAIL_PRINT(printf("Direction %g\n", whichWay));
329
330	numberInfeasibilities_ = 0;
331	changeCost_ = 0.0;
332	feasibleCost_ = 0.0;
333	double infeasibilityCost = model_->infeasibilityCost();
334	infeasibilityWeight_ = infeasibilityCost;;
335	largestInfeasibility_ = 0.0;
336	sumInfeasibilities_ = 0.0;
337
338	int iSequence;
339	assert (!model_->rowObjective());
340	double * cost = model_->objective();
341
342	// First see how much space we need
343	// Also set up feasible bounds
344	int put = starts[numberColumns_];
345
346	double * columnUpper = model_->columnUpper();
347	double * columnLower = model_->columnLower();
348	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
349	if (columnLower[iSequence] > -1.0e20)
350	put++;
351	if (columnUpper[iSequence] < 1.0e20)
352	put++;
353	}
354
355	double * rowUpper = model_->rowUpper();
356	double * rowLower = model_->rowLower();
357	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
358	if (rowLower[iSequence] > -1.0e20)
359	put++;
360	if (rowUpper[iSequence] < 1.0e20)
361	put++;
362	put += 2;
363	}
364	lower_ = new double [put];
365	cost_ = new double [put];
366	infeasible_ = new unsigned int[(put+31)>>5];
367	memset(infeasible_, 0, ((put + 31) >> 5)*sizeof(unsigned int));
368
369	// now fill in
370	put = 0;
371
372	start_[0] = 0;
373	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
374	lower_[put] = -COIN_DBL_MAX;
375	whichRange_[iSequence] = put + 1;
376	double thisCost;
377	double lowerValue;
378	double upperValue;
379	if (iSequence >= numberColumns_) {
380	// rows
381	lowerValue = rowLower[iSequence-numberColumns_];
382	upperValue = rowUpper[iSequence-numberColumns_];
383	if (lowerValue > -1.0e30) {
384	setInfeasible(put, true);
385	cost_[put++] = -infeasibilityCost;
386	lower_[put] = lowerValue;
387	}
388	cost_[put++] = 0.0;
389	thisCost = 0.0;
390	} else {
391	// columns - move costs and see if convex
392	lowerValue = columnLower[iSequence];
393	upperValue = columnUpper[iSequence];
394	if (lowerValue > -1.0e30) {
395	setInfeasible(put, true);
396	cost_[put++] = whichWay * cost[iSequence] - infeasibilityCost;
397	lower_[put] = lowerValue;
398	}
399	int iIndex = starts[iSequence];
400	int end = starts[iSequence+1];
401	assert (fabs(columnLower[iSequence] - lowerNon[iIndex]) < 1.0e-8);
402	thisCost = -COIN_DBL_MAX;
403	for (; iIndex < end; iIndex++) {
404	if (lowerNon[iIndex] < columnUpper[iSequence] - 1.0e-8) {
405	lower_[put] = lowerNon[iIndex];
406	cost_[put++] = whichWay * costNon[iIndex];
407	// check convexity
408	if (whichWay * costNon[iIndex] < thisCost - 1.0e-12)
409	convex_ = false;
410	thisCost = whichWay * costNon[iIndex];
411	} else {
412	break;
413	}
414	}
415	}
416	lower_[put] = upperValue;
417	setInfeasible(put, true);
418	cost_[put++] = thisCost + infeasibilityCost;
419	if (upperValue < 1.0e20) {
420	lower_[put] = COIN_DBL_MAX;
421	cost_[put++] = 1.0e50;
422	}
423	int iFirst = start_[iSequence];
424	if (lower_[iFirst] != -COIN_DBL_MAX) {
425	setInfeasible(iFirst, true);
426	whichRange_[iSequence] = iFirst + 1;
427	} else {
428	whichRange_[iSequence] = iFirst;
429	}
430	start_[iSequence+1] = put;
431	}
432	// can't handle non-convex at present
433	assert(convex_);
434	status_ = NULL;
435	bound_ = NULL;
436	cost2_ = NULL;
437	method_ = 1;
438	#else
439	printf("recompile ClpNonLinearCost without FAST_CLPNON\n");
440	abort();
441	#endif
442	}
443	//-------------------------------------------------------------------
444	// Copy constructor
445	//-------------------------------------------------------------------
446	ClpNonLinearCost::ClpNonLinearCost (const ClpNonLinearCost & rhs) :
447	changeCost_(0.0),
448	feasibleCost_(0.0),
449	infeasibilityWeight_(-1.0),
450	largestInfeasibility_(0.0),
451	sumInfeasibilities_(0.0),
452	averageTheta_(0.0),
453	numberRows_(rhs.numberRows_),
454	numberColumns_(rhs.numberColumns_),
455	start_(NULL),
456	whichRange_(NULL),
457	offset_(NULL),
458	lower_(NULL),
459	cost_(NULL),
460	model_(NULL),
461	infeasible_(NULL),
462	numberInfeasibilities_(-1),
463	status_(NULL),
464	bound_(NULL),
465	cost2_(NULL),
466	method_(rhs.method_),
467	convex_(true),
468	bothWays_(rhs.bothWays_)
469	{
470	if (numberRows_) {
471	int numberTotal = numberRows_ + numberColumns_;
472	model_ = rhs.model_;
473	numberInfeasibilities_ = rhs.numberInfeasibilities_;
474	changeCost_ = rhs.changeCost_;
475	feasibleCost_ = rhs.feasibleCost_;
476	infeasibilityWeight_ = rhs.infeasibilityWeight_;
477	largestInfeasibility_ = rhs.largestInfeasibility_;
478	sumInfeasibilities_ = rhs.sumInfeasibilities_;
479	averageTheta_ = rhs.averageTheta_;
480	convex_ = rhs.convex_;
481	if (CLP_METHOD1) {
482	start_ = new int [numberTotal+1];
483	CoinMemcpyN(rhs.start_, (numberTotal + 1), start_);
484	whichRange_ = new int [numberTotal];
485	CoinMemcpyN(rhs.whichRange_, numberTotal, whichRange_);
486	offset_ = new int [numberTotal];
487	CoinMemcpyN(rhs.offset_, numberTotal, offset_);
488	int numberEntries = start_[numberTotal];
489	lower_ = new double [numberEntries];
490	CoinMemcpyN(rhs.lower_, numberEntries, lower_);
491	cost_ = new double [numberEntries];
492	CoinMemcpyN(rhs.cost_, numberEntries, cost_);
493	infeasible_ = new unsigned int[(numberEntries+31)>>5];
494	CoinMemcpyN(rhs.infeasible_, ((numberEntries + 31) >> 5), infeasible_);
495	}
496	if (CLP_METHOD2) {
497	bound_ = CoinCopyOfArray(rhs.bound_, numberTotal);
498	cost2_ = CoinCopyOfArray(rhs.cost2_, numberTotal);
499	status_ = CoinCopyOfArray(rhs.status_, numberTotal);
500	}
501	}
502	}
503
504	//-------------------------------------------------------------------
505	// Destructor
506	//-------------------------------------------------------------------
507	ClpNonLinearCost::~ClpNonLinearCost ()
508	{
509	delete [] start_;
510	delete [] whichRange_;
511	delete [] offset_;
512	delete [] lower_;
513	delete [] cost_;
514	delete [] infeasible_;
515	delete [] status_;
516	delete [] bound_;
517	delete [] cost2_;
518	}
519
520	//----------------------------------------------------------------
521	// Assignment operator
522	//-------------------------------------------------------------------
523	ClpNonLinearCost &
524	ClpNonLinearCost::operator=(const ClpNonLinearCost& rhs)
525	{
526	if (this != &rhs) {
527	numberRows_ = rhs.numberRows_;
528	numberColumns_ = rhs.numberColumns_;
529	delete [] start_;
530	delete [] whichRange_;
531	delete [] offset_;
532	delete [] lower_;
533	delete [] cost_;
534	delete [] infeasible_;
535	delete [] status_;
536	delete [] bound_;
537	delete [] cost2_;
538	start_ = NULL;
539	whichRange_ = NULL;
540	lower_ = NULL;
541	cost_ = NULL;
542	infeasible_ = NULL;
543	status_ = NULL;
544	bound_ = NULL;
545	cost2_ = NULL;
546	method_ = rhs.method_;
547	if (numberRows_) {
548	int numberTotal = numberRows_ + numberColumns_;
549	if (CLP_METHOD1) {
550	start_ = new int [numberTotal+1];
551	CoinMemcpyN(rhs.start_, (numberTotal + 1), start_);
552	whichRange_ = new int [numberTotal];
553	CoinMemcpyN(rhs.whichRange_, numberTotal, whichRange_);
554	offset_ = new int [numberTotal];
555	CoinMemcpyN(rhs.offset_, numberTotal, offset_);
556	int numberEntries = start_[numberTotal];
557	lower_ = new double [numberEntries];
558	CoinMemcpyN(rhs.lower_, numberEntries, lower_);
559	cost_ = new double [numberEntries];
560	CoinMemcpyN(rhs.cost_, numberEntries, cost_);
561	infeasible_ = new unsigned int[(numberEntries+31)>>5];
562	CoinMemcpyN(rhs.infeasible_, ((numberEntries + 31) >> 5), infeasible_);
563	}
564	if (CLP_METHOD2) {
565	bound_ = CoinCopyOfArray(rhs.bound_, numberTotal);
566	cost2_ = CoinCopyOfArray(rhs.cost2_, numberTotal);
567	status_ = CoinCopyOfArray(rhs.status_, numberTotal);
568	}
569	}
570	model_ = rhs.model_;
571	numberInfeasibilities_ = rhs.numberInfeasibilities_;
572	changeCost_ = rhs.changeCost_;
573	feasibleCost_ = rhs.feasibleCost_;
574	infeasibilityWeight_ = rhs.infeasibilityWeight_;
575	largestInfeasibility_ = rhs.largestInfeasibility_;
576	sumInfeasibilities_ = rhs.sumInfeasibilities_;
577	averageTheta_ = rhs.averageTheta_;
578	convex_ = rhs.convex_;
579	bothWays_ = rhs.bothWays_;
580	}
581	return *this;
582	}
583	// Changes infeasible costs and computes number and cost of infeas
584	// We will need to re-think objective offsets later
585	// We will also need a 2 bit per variable array for some
586	// purpose which will come to me later
587	void
588	ClpNonLinearCost::checkInfeasibilities(double oldTolerance)
589	{
590	numberInfeasibilities_ = 0;
591	double infeasibilityCost = model_->infeasibilityCost();
592	changeCost_ = 0.0;
593	largestInfeasibility_ = 0.0;
594	sumInfeasibilities_ = 0.0;
595	double primalTolerance = model_->currentPrimalTolerance();
596	int iSequence;
597	double * solution = model_->solutionRegion();
598	double * upper = model_->upperRegion();
599	double * lower = model_->lowerRegion();
600	double * cost = model_->costRegion();
601	bool toNearest = oldTolerance <= 0.0;
602	feasibleCost_ = 0.0;
603	//bool checkCosts = (infeasibilityWeight_ != infeasibilityCost);
604	infeasibilityWeight_ = infeasibilityCost;
605	int numberTotal = numberColumns_ + numberRows_;
606	//#define NONLIN_DEBUG
607	#ifdef NONLIN_DEBUG
608	double * saveSolution = NULL;
609	double * saveLower = NULL;
610	double * saveUpper = NULL;
611	unsigned char * saveStatus = NULL;
612	if (method_ == 3) {
613	// Save solution as we will be checking
614	saveSolution = CoinCopyOfArray(solution, numberTotal);
615	saveLower = CoinCopyOfArray(lower, numberTotal);
616	saveUpper = CoinCopyOfArray(upper, numberTotal);
617	saveStatus = CoinCopyOfArray(model_->statusArray(), numberTotal);
618	}
619	#else
620	assert (method_ != 3);
621	#endif
622	if (CLP_METHOD1) {
623	// nonbasic should be at a valid bound
624	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
625	double lowerValue;
626	double upperValue;
627	double value = solution[iSequence];
628	int iRange;
629	// get correct place
630	int start = start_[iSequence];
631	int end = start_[iSequence+1] - 1;
632	// correct costs for this infeasibility weight
633	// If free then true cost will be first
634	double thisFeasibleCost = cost_[start];
635	if (infeasible(start)) {
636	thisFeasibleCost = cost_[start+1];
637	cost_[start] = thisFeasibleCost - infeasibilityCost;
638	}
639	if (infeasible(end - 1)) {
640	thisFeasibleCost = cost_[end-2];
641	cost_[end-1] = thisFeasibleCost + infeasibilityCost;
642	}
643	for (iRange = start; iRange < end; iRange++) {
644	if (value < lower_[iRange+1] + primalTolerance) {
645	// put in better range if infeasible
646	if (value >= lower_[iRange+1] - primalTolerance && infeasible(iRange) && iRange == start)
647	iRange++;
648	whichRange_[iSequence] = iRange;
649	break;
650	}
651	}
652	assert(iRange < end);
653	lowerValue = lower_[iRange];
654	upperValue = lower_[iRange+1];
655	ClpSimplex::Status status = model_->getStatus(iSequence);
656	if (upperValue == lowerValue && status != ClpSimplex::isFixed) {
657	if (status != ClpSimplex::basic) {
658	model_->setStatus(iSequence, ClpSimplex::isFixed);
659	status = ClpSimplex::isFixed;
660	}
661	}
662	//#define PRINT_DETAIL7 2
663	#if PRINT_DETAIL7>1
664	printf("NL %d sol %g bounds %g %g\n",
665	iSequence,solution[iSequence],
666	lowerValue,upperValue);
667	#endif
668	switch(status) {
669
670	case ClpSimplex::basic:
671	case ClpSimplex::superBasic:
672	// iRange is in correct place
673	// slot in here
674	if (infeasible(iRange)) {
675	if (lower_[iRange] < -1.0e50) {
676	//cost_[iRange] = cost_[iRange+1]-infeasibilityCost;
677	// possibly below
678	lowerValue = lower_[iRange+1];
679	if (value - lowerValue < -primalTolerance) {
680	value = lowerValue - value - primalTolerance;
681	#ifndef NDEBUG
682	if(value > 1.0e15)
683	printf("nonlincostb %d %g %g %g\n",
684	iSequence, lowerValue, solution[iSequence], lower_[iRange+2]);
685	#endif
686	#if PRINT_DETAIL7
687	printf("**NL %d sol %g below %g\n",
688	iSequence,solution[iSequence],
689	lowerValue);
690	#endif
691	sumInfeasibilities_ += value;
692	largestInfeasibility_ = CoinMax(largestInfeasibility_, value);
693	changeCost_ -= lowerValue *
694	(cost_[iRange] - cost[iSequence]);
695	numberInfeasibilities_++;
696	}
697	} else {
698	//cost_[iRange] = cost_[iRange-1]+infeasibilityCost;
699	// possibly above
700	upperValue = lower_[iRange];
701	if (value - upperValue > primalTolerance) {
702	value = value - upperValue - primalTolerance;
703	#ifndef NDEBUG
704	if(value > 1.0e15)
705	printf("nonlincostu %d %g %g %g\n",
706	iSequence, lower_[iRange-1], solution[iSequence], upperValue);
707	#endif
708	#if PRINT_DETAIL7
709	printf("**NL %d sol %g above %g\n",
710	iSequence,solution[iSequence],
711	upperValue);
712	#endif
713	sumInfeasibilities_ += value;
714	largestInfeasibility_ = CoinMax(largestInfeasibility_, value);
715	changeCost_ -= upperValue *
716	(cost_[iRange] - cost[iSequence]);
717	numberInfeasibilities_++;
718	}
719	}
720	}
721	//lower[iSequence] = lower_[iRange];
722	//upper[iSequence] = lower_[iRange+1];
723	//cost[iSequence] = cost_[iRange];
724	break;
725	case ClpSimplex::isFree:
726	//if (toNearest)
727	//solution[iSequence] = 0.0;
728	break;
729	case ClpSimplex::atUpperBound:
730	if (!toNearest) {
731	// With increasing tolerances - we may be at wrong place
732	if (fabs(value - upperValue) > oldTolerance * 1.0001) {
733	if (fabs(value - lowerValue) <= oldTolerance * 1.0001) {
734	if (fabs(value - lowerValue) > primalTolerance)
735	solution[iSequence] = lowerValue;
736	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
737	} else {
738	model_->setStatus(iSequence, ClpSimplex::superBasic);
739	}
740	} else if (fabs(value - upperValue) > primalTolerance) {
741	solution[iSequence] = upperValue;
742	}
743	} else {
744	// Set to nearest and make at upper bound
745	int kRange;
746	iRange = -1;
747	double nearest = COIN_DBL_MAX;
748	for (kRange = start; kRange < end; kRange++) {
749	if (fabs(lower_[kRange] - value) < nearest) {
750	nearest = fabs(lower_[kRange] - value);
751	iRange = kRange;
752	}
753	}
754	assert (iRange >= 0);
755	iRange--;
756	whichRange_[iSequence] = iRange;
757	solution[iSequence] = lower_[iRange+1];
758	}
759	break;
760	case ClpSimplex::atLowerBound:
761	if (!toNearest) {
762	// With increasing tolerances - we may be at wrong place
763	// below stops compiler error with gcc 3.2!!!
764	if (iSequence == -119)
765	printf("ZZ %g %g %g %g\n", lowerValue, value, upperValue, oldTolerance);
766	if (fabs(value - lowerValue) > oldTolerance * 1.0001) {
767	if (fabs(value - upperValue) <= oldTolerance * 1.0001) {
768	if (fabs(value - upperValue) > primalTolerance)
769	solution[iSequence] = upperValue;
770	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
771	} else {
772	model_->setStatus(iSequence, ClpSimplex::superBasic);
773	}
774	} else if (fabs(value - lowerValue) > primalTolerance) {
775	solution[iSequence] = lowerValue;
776	}
777	} else {
778	// Set to nearest and make at lower bound
779	int kRange;
780	iRange = -1;
781	double nearest = COIN_DBL_MAX;
782	for (kRange = start; kRange < end; kRange++) {
783	if (fabs(lower_[kRange] - value) < nearest) {
784	nearest = fabs(lower_[kRange] - value);
785	iRange = kRange;
786	}
787	}
788	assert (iRange >= 0);
789	whichRange_[iSequence] = iRange;
790	solution[iSequence] = lower_[iRange];
791	}
792	break;
793	case ClpSimplex::isFixed:
794	if (toNearest) {
795	// Set to true fixed
796	for (iRange = start; iRange < end; iRange++) {
797	if (lower_[iRange] == lower_[iRange+1])
798	break;
799	}
800	if (iRange == end) {
801	// Odd - but make sensible
802	// Set to nearest and make at lower bound
803	int kRange;
804	iRange = -1;
805	double nearest = COIN_DBL_MAX;
806	for (kRange = start; kRange < end; kRange++) {
807	if (fabs(lower_[kRange] - value) < nearest) {
808	nearest = fabs(lower_[kRange] - value);
809	iRange = kRange;
810	}
811	}
812	assert (iRange >= 0);
813	whichRange_[iSequence] = iRange;
814	if (lower_[iRange] != lower_[iRange+1])
815	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
816	else
817	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
818	}
819	solution[iSequence] = lower_[iRange];
820	}
821	break;
822	}
823	lower[iSequence] = lower_[iRange];
824	upper[iSequence] = lower_[iRange+1];
825	cost[iSequence] = cost_[iRange];
826	feasibleCost_ += thisFeasibleCost * solution[iSequence];
827	//assert (iRange==whichRange_[iSequence]);
828	}
829	}
830	#ifdef NONLIN_DEBUG
831	double saveCost = feasibleCost_;
832	if (method_ == 3) {
833	feasibleCost_ = 0.0;
834	// Put back solution as we will be checking
835	unsigned char * statusA = model_->statusArray();
836	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
837	double value = solution[iSequence];
838	solution[iSequence] = saveSolution[iSequence];
839	saveSolution[iSequence] = value;
840	value = lower[iSequence];
841	lower[iSequence] = saveLower[iSequence];
842	saveLower[iSequence] = value;
843	value = upper[iSequence];
844	upper[iSequence] = saveUpper[iSequence];
845	saveUpper[iSequence] = value;
846	unsigned char value2 = statusA[iSequence];
847	statusA[iSequence] = saveStatus[iSequence];
848	saveStatus[iSequence] = value2;
849	}
850	}
851	#endif
852	if (CLP_METHOD2) {
853	//#define CLP_NON_JUST_BASIC
854	#ifndef CLP_NON_JUST_BASIC
855	// nonbasic should be at a valid bound
856	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
857	#else
858	const int * pivotVariable = model_->pivotVariable();
859	for (int i=0;i<numberRows_;i++) {
860	int iSequence = pivotVariable[i];
861	#endif
862	double value = solution[iSequence];
863	unsigned char iStatus = status_[iSequence];
864	assert (currentStatus(iStatus) == CLP_SAME);
865	double lowerValue = lower[iSequence];
866	double upperValue = upper[iSequence];
867	double costValue = cost2_[iSequence];
868	double trueCost = costValue;
869	int iWhere = originalStatus(iStatus);
870	if (iWhere == CLP_BELOW_LOWER) {
871	lowerValue = upperValue;
872	upperValue = bound_[iSequence];
873	costValue -= infeasibilityCost;
874	} else if (iWhere == CLP_ABOVE_UPPER) {
875	upperValue = lowerValue;
876	lowerValue = bound_[iSequence];
877	costValue += infeasibilityCost;
878	}
879	// get correct place
880	int newWhere = CLP_FEASIBLE;
881	ClpSimplex::Status status = model_->getStatus(iSequence);
882	if (upperValue == lowerValue && status != ClpSimplex::isFixed) {
883	if (status != ClpSimplex::basic) {
884	model_->setStatus(iSequence, ClpSimplex::isFixed);
885	status = ClpSimplex::isFixed;
886	}
887	}
888	switch(status) {
889
890	case ClpSimplex::basic:
891	case ClpSimplex::superBasic:
892	if (value - upperValue <= primalTolerance) {
893	if (value - lowerValue >= -primalTolerance) {
894	// feasible
895	//newWhere=CLP_FEASIBLE;
896	} else {
897	// below
898	newWhere = CLP_BELOW_LOWER;
899	assert (fabs(lowerValue) < 1.0e100);
900	double infeasibility = lowerValue - value - primalTolerance;
901	sumInfeasibilities_ += infeasibility;
902	largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
903	costValue = trueCost - infeasibilityCost;
904	changeCost_ -= lowerValue * (costValue - cost[iSequence]);
905	numberInfeasibilities_++;
906	}
907	} else {
908	// above
909	newWhere = CLP_ABOVE_UPPER;
910	double infeasibility = value - upperValue - primalTolerance;
911	sumInfeasibilities_ += infeasibility;
912	largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
913	costValue = trueCost + infeasibilityCost;
914	changeCost_ -= upperValue * (costValue - cost[iSequence]);
915	numberInfeasibilities_++;
916	}
917	break;
918	case ClpSimplex::isFree:
919	break;
920	case ClpSimplex::atUpperBound:
921	if (!toNearest) {
922	// With increasing tolerances - we may be at wrong place
923	if (fabs(value - upperValue) > oldTolerance * 1.0001) {
924	if (fabs(value - lowerValue) <= oldTolerance * 1.0001) {
925	if (fabs(value - lowerValue) > primalTolerance) {
926	solution[iSequence] = lowerValue;
927	value = lowerValue;
928	}
929	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
930	} else {
931	if (value < upperValue) {
932	if (value > lowerValue) {
933	model_->setStatus(iSequence, ClpSimplex::superBasic);
934	} else {
935	// set to lower bound as infeasible
936	solution[iSequence] = lowerValue;
937	value = lowerValue;
938	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
939	}
940	} else {
941	// set to upper bound as infeasible
942	solution[iSequence] = upperValue;
943	value = upperValue;
944	}
945	}
946	} else if (fabs(value - upperValue) > primalTolerance) {
947	solution[iSequence] = upperValue;
948	value = upperValue;
949	}
950	} else {
951	// Set to nearest and make at bound
952	if (fabs(value - lowerValue) < fabs(value - upperValue)) {
953	solution[iSequence] = lowerValue;
954	value = lowerValue;
955	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
956	} else {
957	solution[iSequence] = upperValue;
958	value = upperValue;
959	}
960	}
961	break;
962	case ClpSimplex::atLowerBound:
963	if (!toNearest) {
964	// With increasing tolerances - we may be at wrong place
965	if (fabs(value - lowerValue) > oldTolerance * 1.0001) {
966	if (fabs(value - upperValue) <= oldTolerance * 1.0001) {
967	if (fabs(value - upperValue) > primalTolerance) {
968	solution[iSequence] = upperValue;
969	value = upperValue;
970	}
971	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
972	} else {
973	if (value < upperValue) {
974	if (value > lowerValue) {
975	model_->setStatus(iSequence, ClpSimplex::superBasic);
976	} else {
977	// set to lower bound as infeasible
978	solution[iSequence] = lowerValue;
979	value = lowerValue;
980	}
981	} else {
982	// set to upper bound as infeasible
983	solution[iSequence] = upperValue;
984	value = upperValue;
985	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
986	}
987	}
988	} else if (fabs(value - lowerValue) > primalTolerance) {
989	solution[iSequence] = lowerValue;
990	value = lowerValue;
991	}
992	} else {
993	// Set to nearest and make at bound
994	if (fabs(value - lowerValue) < fabs(value - upperValue)) {
995	solution[iSequence] = lowerValue;
996	value = lowerValue;
997	} else {
998	solution[iSequence] = upperValue;
999	value = upperValue;
1000	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
1001	}
1002	}
1003	break;
1004	case ClpSimplex::isFixed:
1005	solution[iSequence] = lowerValue;
1006	value = lowerValue;
1007	break;
1008	}
1009	#ifdef NONLIN_DEBUG
1010	double lo = saveLower[iSequence];
1011	double up = saveUpper[iSequence];
1012	double cc = cost[iSequence];
1013	unsigned char ss = saveStatus[iSequence];
1014	unsigned char snow = model_->statusArray()[iSequence];
1015	#endif
1016	if (iWhere != newWhere) {
1017	setOriginalStatus(status_[iSequence], newWhere);
1018	if (newWhere == CLP_BELOW_LOWER) {
1019	bound_[iSequence] = upperValue;
1020	upperValue = lowerValue;
1021	lowerValue = -COIN_DBL_MAX;
1022	costValue = trueCost - infeasibilityCost;
1023	} else if (newWhere == CLP_ABOVE_UPPER) {
1024	bound_[iSequence] = lowerValue;
1025	lowerValue = upperValue;
1026	upperValue = COIN_DBL_MAX;
1027	costValue = trueCost + infeasibilityCost;
1028	} else {
1029	costValue = trueCost;
1030	}
1031	lower[iSequence] = lowerValue;
1032	upper[iSequence] = upperValue;
1033	}
1034	// always do as other things may change
1035	cost[iSequence] = costValue;
1036	#ifdef NONLIN_DEBUG
1037	if (method_ == 3) {
1038	assert (ss == snow);
1039	assert (cc == cost[iSequence]);
1040	assert (lo == lower[iSequence]);
1041	assert (up == upper[iSequence]);
1042	assert (value == saveSolution[iSequence]);
1043	}
1044	#endif
1045	feasibleCost_ += trueCost * value;
1046	}
1047	}
1048	#ifdef NONLIN_DEBUG
1049	if (method_ == 3)
1050	assert (fabs(saveCost - feasibleCost_) < 0.001 * (1.0 + CoinMax(fabs(saveCost), fabs(feasibleCost_))));
1051	delete [] saveSolution;
1052	delete [] saveLower;
1053	delete [] saveUpper;
1054	delete [] saveStatus;
1055	#endif
1056	//feasibleCost_ /= (model_->rhsScale()*model_->objScale());
1057	}
1058	// Puts feasible bounds into lower and upper
1059	void
1060	ClpNonLinearCost::feasibleBounds()
1061	{
1062	if (CLP_METHOD2) {
1063	int iSequence;
1064	double * upper = model_->upperRegion();
1065	double * lower = model_->lowerRegion();
1066	double * cost = model_->costRegion();
1067	int numberTotal = numberColumns_ + numberRows_;
1068	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1069	unsigned char iStatus = status_[iSequence];
1070	assert (currentStatus(iStatus) == CLP_SAME);
1071	double lowerValue = lower[iSequence];
1072	double upperValue = upper[iSequence];
1073	double costValue = cost2_[iSequence];
1074	int iWhere = originalStatus(iStatus);
1075	if (iWhere == CLP_BELOW_LOWER) {
1076	lowerValue = upperValue;
1077	upperValue = bound_[iSequence];
1078	assert (fabs(lowerValue) < 1.0e100);
1079	} else if (iWhere == CLP_ABOVE_UPPER) {
1080	upperValue = lowerValue;
1081	lowerValue = bound_[iSequence];
1082	}
1083	setOriginalStatus(status_[iSequence], CLP_FEASIBLE);
1084	lower[iSequence] = lowerValue;
1085	upper[iSequence] = upperValue;
1086	cost[iSequence] = costValue;
1087	}
1088	}
1089	}
1090	/* Goes through one bound for each variable.
1091	If array[i]*multiplier>0 goes down, otherwise up.
1092	The indices are row indices and need converting to sequences
1093	*/
1094	void
1095	ClpNonLinearCost::goThru(int numberInArray, double multiplier,
1096	const int * index, const double * array,
1097	double * rhs)
1098	{
1099	assert (model_ != NULL);
1100	abort();
1101	const int * pivotVariable = model_->pivotVariable();
1102	if (CLP_METHOD1) {
1103	for (int i = 0; i < numberInArray; i++) {
1104	int iRow = index[i];
1105	int iSequence = pivotVariable[iRow];
1106	double alpha = multiplier * array[iRow];
1107	// get where in bound sequence
1108	int iRange = whichRange_[iSequence];
1109	iRange += offset_[iSequence]; //add temporary bias
1110	double value = model_->solution(iSequence);
1111	if (alpha > 0.0) {
1112	// down one
1113	iRange--;
1114	assert(iRange >= start_[iSequence]);
1115	rhs[iRow] = value - lower_[iRange];
1116	} else {
1117	// up one
1118	iRange++;
1119	assert(iRange < start_[iSequence+1] - 1);
1120	rhs[iRow] = lower_[iRange+1] - value;
1121	}
1122	offset_[iSequence] = iRange - whichRange_[iSequence];
1123	}
1124	}
1125	#ifdef NONLIN_DEBUG
1126	double * saveRhs = NULL;
1127	if (method_ == 3) {
1128	int numberRows = model_->numberRows();
1129	saveRhs = CoinCopyOfArray(rhs, numberRows);
1130	}
1131	#endif
1132	if (CLP_METHOD2) {
1133	const double * solution = model_->solutionRegion();
1134	const double * upper = model_->upperRegion();
1135	const double * lower = model_->lowerRegion();
1136	for (int i = 0; i < numberInArray; i++) {
1137	int iRow = index[i];
1138	int iSequence = pivotVariable[iRow];
1139	double alpha = multiplier * array[iRow];
1140	double value = solution[iSequence];
1141	unsigned char iStatus = status_[iSequence];
1142	int iWhere = currentStatus(iStatus);
1143	if (iWhere == CLP_SAME)
1144	iWhere = originalStatus(iStatus);
1145	if (iWhere == CLP_FEASIBLE) {
1146	if (alpha > 0.0) {
1147	// going below
1148	iWhere = CLP_BELOW_LOWER;
1149	rhs[iRow] = value - lower[iSequence];
1150	} else {
1151	// going above
1152	iWhere = CLP_ABOVE_UPPER;
1153	rhs[iRow] = upper[iSequence] - value;
1154	}
1155	} else if(iWhere == CLP_BELOW_LOWER) {
1156	assert (alpha < 0);
1157	// going feasible
1158	iWhere = CLP_FEASIBLE;
1159	rhs[iRow] = upper[iSequence] - value;
1160	} else {
1161	assert (iWhere == CLP_ABOVE_UPPER);
1162	// going feasible
1163	iWhere = CLP_FEASIBLE;
1164	rhs[iRow] = value - lower[iSequence];
1165	}
1166	#ifdef NONLIN_DEBUG
1167	if (method_ == 3)
1168	assert (rhs[iRow] == saveRhs[iRow]);
1169	#endif
1170	setCurrentStatus(status_[iSequence], iWhere);
1171	}
1172	}
1173	#ifdef NONLIN_DEBUG
1174	delete [] saveRhs;
1175	#endif
1176	}
1177	/* Takes off last iteration (i.e. offsets closer to 0)
1178	*/
1179	void
1180	ClpNonLinearCost::goBack(int numberInArray, const int * index,
1181	double * rhs)
1182	{
1183	assert (model_ != NULL);
1184	abort();
1185	const int * pivotVariable = model_->pivotVariable();
1186	if (CLP_METHOD1) {
1187	for (int i = 0; i < numberInArray; i++) {
1188	int iRow = index[i];
1189	int iSequence = pivotVariable[iRow];
1190	// get where in bound sequence
1191	int iRange = whichRange_[iSequence];
1192	// get closer to original
1193	if (offset_[iSequence] > 0) {
1194	offset_[iSequence]--;
1195	assert (offset_[iSequence] >= 0);
1196	iRange += offset_[iSequence]; //add temporary bias
1197	double value = model_->solution(iSequence);
1198	// up one
1199	assert(iRange < start_[iSequence+1] - 1);
1200	rhs[iRow] = lower_[iRange+1] - value; // was earlier lower_[iRange]
1201	} else {
1202	offset_[iSequence]++;
1203	assert (offset_[iSequence] <= 0);
1204	iRange += offset_[iSequence]; //add temporary bias
1205	double value = model_->solution(iSequence);
1206	// down one
1207	assert(iRange >= start_[iSequence]);
1208	rhs[iRow] = value - lower_[iRange]; // was earlier lower_[iRange+1]
1209	}
1210	}
1211	}
1212	#ifdef NONLIN_DEBUG
1213	double * saveRhs = NULL;
1214	if (method_ == 3) {
1215	int numberRows = model_->numberRows();
1216	saveRhs = CoinCopyOfArray(rhs, numberRows);
1217	}
1218	#endif
1219	if (CLP_METHOD2) {
1220	const double * solution = model_->solutionRegion();
1221	const double * upper = model_->upperRegion();
1222	const double * lower = model_->lowerRegion();
1223	for (int i = 0; i < numberInArray; i++) {
1224	int iRow = index[i];
1225	int iSequence = pivotVariable[iRow];
1226	double value = solution[iSequence];
1227	unsigned char iStatus = status_[iSequence];
1228	int iWhere = currentStatus(iStatus);
1229	int original = originalStatus(iStatus);
1230	assert (iWhere != CLP_SAME);
1231	if (iWhere == CLP_FEASIBLE) {
1232	if (original == CLP_BELOW_LOWER) {
1233	// going below
1234	iWhere = CLP_BELOW_LOWER;
1235	rhs[iRow] = lower[iSequence] - value;
1236	} else {
1237	// going above
1238	iWhere = CLP_ABOVE_UPPER;
1239	rhs[iRow] = value - upper[iSequence];
1240	}
1241	} else if(iWhere == CLP_BELOW_LOWER) {
1242	// going feasible
1243	iWhere = CLP_FEASIBLE;
1244	rhs[iRow] = value - upper[iSequence];
1245	} else {
1246	// going feasible
1247	iWhere = CLP_FEASIBLE;
1248	rhs[iRow] = lower[iSequence] - value;
1249	}
1250	#ifdef NONLIN_DEBUG
1251	if (method_ == 3)
1252	assert (rhs[iRow] == saveRhs[iRow]);
1253	#endif
1254	setCurrentStatus(status_[iSequence], iWhere);
1255	}
1256	}
1257	#ifdef NONLIN_DEBUG
1258	delete [] saveRhs;
1259	#endif
1260	}
1261	void
1262	ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
1263	{
1264	assert (model_ != NULL);
1265	const int * pivotVariable = model_->pivotVariable();
1266	int number = update->getNumElements();
1267	const int * index = update->getIndices();
1268	if (CLP_METHOD1) {
1269	for (int i = 0; i < number; i++) {
1270	int iRow = index[i];
1271	int iSequence = pivotVariable[iRow];
1272	offset_[iSequence] = 0;
1273	}
1274	#ifdef CLP_DEBUG
1275	for (i = 0; i < numberRows_ + numberColumns_; i++)
1276	assert(!offset_[i]);
1277	#endif
1278	}
1279	if (CLP_METHOD2) {
1280	for (int i = 0; i < number; i++) {
1281	int iRow = index[i];
1282	int iSequence = pivotVariable[iRow];
1283	setSameStatus(status_[iSequence]);
1284	}
1285	}
1286	}
1287	void
1288	ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
1289	{
1290	assert (model_ != NULL);
1291	double primalTolerance = model_->currentPrimalTolerance();
1292	const int * pivotVariable = model_->pivotVariable();
1293	if (CLP_METHOD1) {
1294	for (int i = 0; i < numberInArray; i++) {
1295	int iRow = index[i];
1296	int iSequence = pivotVariable[iRow];
1297	// get where in bound sequence
1298	int iRange;
1299	int currentRange = whichRange_[iSequence];
1300	double value = model_->solution(iSequence);
1301	int start = start_[iSequence];
1302	int end = start_[iSequence+1] - 1;
1303	for (iRange = start; iRange < end; iRange++) {
1304	if (value < lower_[iRange+1] + primalTolerance) {
1305	// put in better range
1306	if (value >= lower_[iRange+1] - primalTolerance && infeasible(iRange) && iRange == start)
1307	iRange++;
1308	break;
1309	}
1310	}
1311	assert(iRange < end);
1312	assert(model_->getStatus(iSequence) == ClpSimplex::basic);
1313	double & lower = model_->lowerAddress(iSequence);
1314	double & upper = model_->upperAddress(iSequence);
1315	double & cost = model_->costAddress(iSequence);
1316	whichRange_[iSequence] = iRange;
1317	if (iRange != currentRange) {
1318	if (infeasible(iRange))
1319	numberInfeasibilities_++;
1320	if (infeasible(currentRange))
1321	numberInfeasibilities_--;
1322	}
1323	lower = lower_[iRange];
1324	upper = lower_[iRange+1];
1325	cost = cost_[iRange];
1326	}
1327	}
1328	if (CLP_METHOD2) {
1329	double * solution = model_->solutionRegion();
1330	double * upper = model_->upperRegion();
1331	double * lower = model_->lowerRegion();
1332	double * cost = model_->costRegion();
1333	for (int i = 0; i < numberInArray; i++) {
1334	int iRow = index[i];
1335	int iSequence = pivotVariable[iRow];
1336	double value = solution[iSequence];
1337	unsigned char iStatus = status_[iSequence];
1338	assert (currentStatus(iStatus) == CLP_SAME);
1339	double lowerValue = lower[iSequence];
1340	double upperValue = upper[iSequence];
1341	double costValue = cost2_[iSequence];
1342	int iWhere = originalStatus(iStatus);
1343	if (iWhere == CLP_BELOW_LOWER) {
1344	lowerValue = upperValue;
1345	upperValue = bound_[iSequence];
1346	numberInfeasibilities_--;
1347	assert (fabs(lowerValue) < 1.0e100);
1348	} else if (iWhere == CLP_ABOVE_UPPER) {
1349	upperValue = lowerValue;
1350	lowerValue = bound_[iSequence];
1351	numberInfeasibilities_--;
1352	}
1353	// get correct place
1354	int newWhere = CLP_FEASIBLE;
1355	if (value - upperValue <= primalTolerance) {
1356	if (value - lowerValue >= -primalTolerance) {
1357	// feasible
1358	//newWhere=CLP_FEASIBLE;
1359	} else {
1360	// below
1361	newWhere = CLP_BELOW_LOWER;
1362	assert (fabs(lowerValue) < 1.0e100);
1363	costValue -= infeasibilityWeight_;
1364	numberInfeasibilities_++;
1365	}
1366	} else {
1367	// above
1368	newWhere = CLP_ABOVE_UPPER;
1369	costValue += infeasibilityWeight_;
1370	numberInfeasibilities_++;
1371	}
1372	if (iWhere != newWhere) {
1373	setOriginalStatus(status_[iSequence], newWhere);
1374	if (newWhere == CLP_BELOW_LOWER) {
1375	bound_[iSequence] = upperValue;
1376	upperValue = lowerValue;
1377	lowerValue = -COIN_DBL_MAX;
1378	} else if (newWhere == CLP_ABOVE_UPPER) {
1379	bound_[iSequence] = lowerValue;
1380	lowerValue = upperValue;
1381	upperValue = COIN_DBL_MAX;
1382	}
1383	lower[iSequence] = lowerValue;
1384	upper[iSequence] = upperValue;
1385	cost[iSequence] = costValue;
1386	}
1387	}
1388	}
1389	}
1390	/* Puts back correct infeasible costs for each variable
1391	The input indices are row indices and need converting to sequences
1392	for costs.
1393	On input array is empty (but indices exist). On exit just
1394	changed costs will be stored as normal CoinIndexedVector
1395	*/
1396	void
1397	ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
1398	{
1399	assert (model_ != NULL);
1400	double primalTolerance = model_->currentPrimalTolerance();
1401	const int * pivotVariable = model_->pivotVariable();
1402	int number = 0;
1403	int * index = update->getIndices();
1404	double * work = update->denseVector();
1405	if (CLP_METHOD1) {
1406	for (int i = 0; i < numberInArray; i++) {
1407	int iRow = index[i];
1408	int iSequence = pivotVariable[iRow];
1409	// get where in bound sequence
1410	int iRange;
1411	double value = model_->solution(iSequence);
1412	int start = start_[iSequence];
1413	int end = start_[iSequence+1] - 1;
1414	for (iRange = start; iRange < end; iRange++) {
1415	if (value < lower_[iRange+1] + primalTolerance) {
1416	// put in better range
1417	if (value >= lower_[iRange+1] - primalTolerance && infeasible(iRange) && iRange == start)
1418	iRange++;
1419	break;
1420	}
1421	}
1422	assert(iRange < end);
1423	assert(model_->getStatus(iSequence) == ClpSimplex::basic);
1424	int jRange = whichRange_[iSequence];
1425	if (iRange != jRange) {
1426	// changed
1427	work[iRow] = cost_[jRange] - cost_[iRange];
1428	index[number++] = iRow;
1429	double & lower = model_->lowerAddress(iSequence);
1430	double & upper = model_->upperAddress(iSequence);
1431	double & cost = model_->costAddress(iSequence);
1432	whichRange_[iSequence] = iRange;
1433	if (infeasible(iRange))
1434	numberInfeasibilities_++;
1435	if (infeasible(jRange))
1436	numberInfeasibilities_--;
1437	lower = lower_[iRange];
1438	upper = lower_[iRange+1];
1439	cost = cost_[iRange];
1440	}
1441	}
1442	}
1443	if (CLP_METHOD2) {
1444	double * solution = model_->solutionRegion();
1445	double * upper = model_->upperRegion();
1446	double * lower = model_->lowerRegion();
1447	double * cost = model_->costRegion();
1448	for (int i = 0; i < numberInArray; i++) {
1449	int iRow = index[i];
1450	int iSequence = pivotVariable[iRow];
1451	double value = solution[iSequence];
1452	unsigned char iStatus = status_[iSequence];
1453	assert (currentStatus(iStatus) == CLP_SAME);
1454	double lowerValue = lower[iSequence];
1455	double upperValue = upper[iSequence];
1456	double costValue = cost2_[iSequence];
1457	int iWhere = originalStatus(iStatus);
1458	if (iWhere == CLP_BELOW_LOWER) {
1459	lowerValue = upperValue;
1460	upperValue = bound_[iSequence];
1461	numberInfeasibilities_--;
1462	assert (fabs(lowerValue) < 1.0e100);
1463	} else if (iWhere == CLP_ABOVE_UPPER) {
1464	upperValue = lowerValue;
1465	lowerValue = bound_[iSequence];
1466	numberInfeasibilities_--;
1467	}
1468	// get correct place
1469	int newWhere = CLP_FEASIBLE;
1470	if (value - upperValue <= primalTolerance) {
1471	if (value - lowerValue >= -primalTolerance) {
1472	// feasible
1473	//newWhere=CLP_FEASIBLE;
1474	} else {
1475	// below
1476	newWhere = CLP_BELOW_LOWER;
1477	costValue -= infeasibilityWeight_;
1478	numberInfeasibilities_++;
1479	assert (fabs(lowerValue) < 1.0e100);
1480	}
1481	} else {
1482	// above
1483	newWhere = CLP_ABOVE_UPPER;
1484	costValue += infeasibilityWeight_;
1485	numberInfeasibilities_++;
1486	}
1487	if (iWhere != newWhere) {
1488	work[iRow] = cost[iSequence] - costValue;
1489	index[number++] = iRow;
1490	setOriginalStatus(status_[iSequence], newWhere);
1491	if (newWhere == CLP_BELOW_LOWER) {
1492	bound_[iSequence] = upperValue;
1493	upperValue = lowerValue;
1494	lowerValue = -COIN_DBL_MAX;
1495	} else if (newWhere == CLP_ABOVE_UPPER) {
1496	bound_[iSequence] = lowerValue;
1497	lowerValue = upperValue;
1498	upperValue = COIN_DBL_MAX;
1499	}
1500	lower[iSequence] = lowerValue;
1501	upper[iSequence] = upperValue;
1502	cost[iSequence] = costValue;
1503	}
1504	}
1505	}
1506	update->setNumElements(number);
1507	}
1508	/* Sets bounds and cost for one variable - returns change in cost*/
1509	double
1510	ClpNonLinearCost::setOne(int iSequence, double value)
1511	{
1512	assert (model_ != NULL);
1513	double primalTolerance = model_->currentPrimalTolerance();
1514	// difference in cost
1515	double difference = 0.0;
1516	if (CLP_METHOD1) {
1517	// get where in bound sequence
1518	int iRange;
1519	int currentRange = whichRange_[iSequence];
1520	int start = start_[iSequence];
1521	int end = start_[iSequence+1] - 1;
1522	if (!bothWays_) {
1523	// If fixed try and get feasible
1524	if (lower_[start+1] == lower_[start+2] && fabs(value - lower_[start+1]) < 1.001 * primalTolerance) {
1525	iRange = start + 1;
1526	} else {
1527	for (iRange = start; iRange < end; iRange++) {
1528	if (value <= lower_[iRange+1] + primalTolerance) {
1529	// put in better range
1530	if (value >= lower_[iRange+1] - primalTolerance && infeasible(iRange) && iRange == start)
1531	iRange++;
1532	break;
1533	}
1534	}
1535	}
1536	} else {
1537	// leave in current if possible
1538	iRange = whichRange_[iSequence];
1539	if (value < lower_[iRange] - primalTolerance \|\| value > lower_[iRange+1] + primalTolerance) {
1540	for (iRange = start; iRange < end; iRange++) {
1541	if (value < lower_[iRange+1] + primalTolerance) {
1542	// put in better range
1543	if (value >= lower_[iRange+1] - primalTolerance && infeasible(iRange) && iRange == start)
1544	iRange++;
1545	break;
1546	}
1547	}
1548	}
1549	}
1550	assert(iRange < end);
1551	whichRange_[iSequence] = iRange;
1552	if (iRange != currentRange) {
1553	if (infeasible(iRange))
1554	numberInfeasibilities_++;
1555	if (infeasible(currentRange))
1556	numberInfeasibilities_--;
1557	}
1558	double & lower = model_->lowerAddress(iSequence);
1559	double & upper = model_->upperAddress(iSequence);
1560	double & cost = model_->costAddress(iSequence);
1561	lower = lower_[iRange];
1562	upper = lower_[iRange+1];
1563	ClpSimplex::Status status = model_->getStatus(iSequence);
1564	if (upper == lower) {
1565	if (status != ClpSimplex::basic) {
1566	model_->setStatus(iSequence, ClpSimplex::isFixed);
1567	status = ClpSimplex::basic; // so will skip
1568	}
1569	}
1570	switch(status) {
1571
1572	case ClpSimplex::basic:
1573	case ClpSimplex::superBasic:
1574	case ClpSimplex::isFree:
1575	break;
1576	case ClpSimplex::atUpperBound:
1577	case ClpSimplex::atLowerBound:
1578	case ClpSimplex::isFixed:
1579	// set correctly
1580	if (fabs(value - lower) <= primalTolerance * 1.001) {
1581	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
1582	} else if (fabs(value - upper) <= primalTolerance * 1.001) {
1583	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
1584	} else {
1585	// set superBasic
1586	model_->setStatus(iSequence, ClpSimplex::superBasic);
1587	}
1588	break;
1589	}
1590	difference = cost - cost_[iRange];
1591	cost = cost_[iRange];
1592	}
1593	if (CLP_METHOD2) {
1594	double * upper = model_->upperRegion();
1595	double * lower = model_->lowerRegion();
1596	double * cost = model_->costRegion();
1597	unsigned char iStatus = status_[iSequence];
1598	assert (currentStatus(iStatus) == CLP_SAME);
1599	double lowerValue = lower[iSequence];
1600	double upperValue = upper[iSequence];
1601	double costValue = cost2_[iSequence];
1602	int iWhere = originalStatus(iStatus);
1603	if (iWhere == CLP_BELOW_LOWER) {
1604	lowerValue = upperValue;
1605	upperValue = bound_[iSequence];
1606	numberInfeasibilities_--;
1607	assert (fabs(lowerValue) < 1.0e100);
1608	} else if (iWhere == CLP_ABOVE_UPPER) {
1609	upperValue = lowerValue;
1610	lowerValue = bound_[iSequence];
1611	numberInfeasibilities_--;
1612	}
1613	// get correct place
1614	int newWhere = CLP_FEASIBLE;
1615	if (value - upperValue <= primalTolerance) {
1616	if (value - lowerValue >= -primalTolerance) {
1617	// feasible
1618	//newWhere=CLP_FEASIBLE;
1619	} else {
1620	// below
1621	newWhere = CLP_BELOW_LOWER;
1622	costValue -= infeasibilityWeight_;
1623	numberInfeasibilities_++;
1624	assert (fabs(lowerValue) < 1.0e100);
1625	}
1626	} else {
1627	// above
1628	newWhere = CLP_ABOVE_UPPER;
1629	costValue += infeasibilityWeight_;
1630	numberInfeasibilities_++;
1631	}
1632	if (iWhere != newWhere) {
1633	difference = cost[iSequence] - costValue;
1634	setOriginalStatus(status_[iSequence], newWhere);
1635	if (newWhere == CLP_BELOW_LOWER) {
1636	bound_[iSequence] = upperValue;
1637	upperValue = lowerValue;
1638	lowerValue = -COIN_DBL_MAX;
1639	} else if (newWhere == CLP_ABOVE_UPPER) {
1640	bound_[iSequence] = lowerValue;
1641	lowerValue = upperValue;
1642	upperValue = COIN_DBL_MAX;
1643	}
1644	lower[iSequence] = lowerValue;
1645	upper[iSequence] = upperValue;
1646	cost[iSequence] = costValue;
1647	}
1648	ClpSimplex::Status status = model_->getStatus(iSequence);
1649	if (upperValue == lowerValue) {
1650	if (status != ClpSimplex::basic) {
1651	model_->setStatus(iSequence, ClpSimplex::isFixed);
1652	status = ClpSimplex::basic; // so will skip
1653	}
1654	}
1655	switch(status) {
1656
1657	case ClpSimplex::basic:
1658	case ClpSimplex::superBasic:
1659	case ClpSimplex::isFree:
1660	break;
1661	case ClpSimplex::atUpperBound:
1662	case ClpSimplex::atLowerBound:
1663	case ClpSimplex::isFixed:
1664	// set correctly
1665	if (fabs(value - lowerValue) <= primalTolerance * 1.001) {
1666	model_->setStatus(iSequence, ClpSimplex::atLowerBound);
1667	} else if (fabs(value - upperValue) <= primalTolerance * 1.001) {
1668	model_->setStatus(iSequence, ClpSimplex::atUpperBound);
1669	} else {
1670	// set superBasic
1671	model_->setStatus(iSequence, ClpSimplex::superBasic);
1672	}
1673	break;
1674	}
1675	}
1676	changeCost_ += value * difference;
1677	return difference;
1678	}
1679	/* Sets bounds and infeasible cost and true cost for one variable
1680	This is for gub and column generation etc */
1681	void
1682	ClpNonLinearCost::setOne(int sequence, double solutionValue, double lowerValue, double upperValue,
1683	double costValue)
1684	{
1685	if (CLP_METHOD1) {
1686	int iRange = -1;
1687	int start = start_[sequence];
1688	double infeasibilityCost = model_->infeasibilityCost();
1689	cost_[start] = costValue - infeasibilityCost;
1690	lower_[start+1] = lowerValue;
1691	cost_[start+1] = costValue;
1692	lower_[start+2] = upperValue;
1693	cost_[start+2] = costValue + infeasibilityCost;
1694	double primalTolerance = model_->currentPrimalTolerance();
1695	if (solutionValue - lowerValue >= -primalTolerance) {
1696	if (solutionValue - upperValue <= primalTolerance) {
1697	iRange = start + 1;
1698	} else {
1699	iRange = start + 2;
1700	}
1701	} else {
1702	iRange = start;
1703	}
1704	model_->costRegion()[sequence] = cost_[iRange];
1705	whichRange_[sequence] = iRange;
1706	}
1707	if (CLP_METHOD2) {
1708	bound_[sequence]=0.0;
1709	cost2_[sequence]=costValue;
1710	setInitialStatus(status_[sequence]);
1711	}
1712
1713	}
1714	/* Sets bounds and cost for outgoing variable
1715	may change value
1716	Returns direction */
1717	int
1718	ClpNonLinearCost::setOneOutgoing(int iSequence, double & value)
1719	{
1720	assert (model_ != NULL);
1721	double primalTolerance = model_->currentPrimalTolerance();
1722	// difference in cost
1723	double difference = 0.0;
1724	int direction = 0;
1725	if (CLP_METHOD1) {
1726	// get where in bound sequence
1727	int iRange;
1728	int currentRange = whichRange_[iSequence];
1729	int start = start_[iSequence];
1730	int end = start_[iSequence+1] - 1;
1731	// Set perceived direction out
1732	if (value <= lower_[currentRange] + 1.001 * primalTolerance) {
1733	direction = 1;
1734	} else if (value >= lower_[currentRange+1] - 1.001 * primalTolerance) {
1735	direction = -1;
1736	} else {
1737	// odd
1738	direction = 0;
1739	}
1740	// If fixed try and get feasible
1741	if (lower_[start+1] == lower_[start+2] && fabs(value - lower_[start+1]) < 1.001 * primalTolerance) {
1742	iRange = start + 1;
1743	} else {
1744	// See if exact
1745	for (iRange = start; iRange < end; iRange++) {
1746	if (value == lower_[iRange+1]) {
1747	// put in better range
1748	if (infeasible(iRange) && iRange == start)
1749	iRange++;
1750	break;
1751	}
1752	}
1753	if (iRange == end) {
1754	// not exact
1755	for (iRange = start; iRange < end; iRange++) {
1756	if (value <= lower_[iRange+1] + primalTolerance) {
1757	// put in better range
1758	if (value >= lower_[iRange+1] - primalTolerance && infeasible(iRange) && iRange == start)
1759	iRange++;
1760	break;
1761	}
1762	}
1763	}
1764	}
1765	assert(iRange < end);
1766	whichRange_[iSequence] = iRange;
1767	if (iRange != currentRange) {
1768	if (infeasible(iRange))
1769	numberInfeasibilities_++;
1770	if (infeasible(currentRange))
1771	numberInfeasibilities_--;
1772	}
1773	double & lower = model_->lowerAddress(iSequence);
1774	double & upper = model_->upperAddress(iSequence);
1775	double & cost = model_->costAddress(iSequence);
1776	lower = lower_[iRange];
1777	upper = lower_[iRange+1];
1778	if (upper == lower) {
1779	value = upper;
1780	} else {
1781	// set correctly
1782	if (fabs(value - lower) <= primalTolerance * 1.001) {
1783	value = CoinMin(value, lower + primalTolerance);
1784	} else if (fabs(value - upper) <= primalTolerance * 1.001) {
1785	value = CoinMax(value, upper - primalTolerance);
1786	} else {
1787	//printf("*** variable wandered off bound %g %g %g!\n",
1788	// lower,value,upper);
1789	if (value - lower <= upper - value)
1790	value = lower + primalTolerance;
1791	else
1792	value = upper - primalTolerance;
1793	}
1794	}
1795	difference = cost - cost_[iRange];
1796	cost = cost_[iRange];
1797	}
1798	if (CLP_METHOD2) {
1799	double * upper = model_->upperRegion();
1800	double * lower = model_->lowerRegion();
1801	double * cost = model_->costRegion();
1802	unsigned char iStatus = status_[iSequence];
1803	assert (currentStatus(iStatus) == CLP_SAME);
1804	double lowerValue = lower[iSequence];
1805	double upperValue = upper[iSequence];
1806	double costValue = cost2_[iSequence];
1807	// Set perceived direction out
1808	if (value <= lowerValue + 1.001 * primalTolerance) {
1809	direction = 1;
1810	} else if (value >= upperValue - 1.001 * primalTolerance) {
1811	direction = -1;
1812	} else {
1813	// odd
1814	direction = 0;
1815	}
1816	int iWhere = originalStatus(iStatus);
1817	if (iWhere == CLP_BELOW_LOWER) {
1818	lowerValue = upperValue;
1819	upperValue = bound_[iSequence];
1820	numberInfeasibilities_--;
1821	assert (fabs(lowerValue) < 1.0e100);
1822	} else if (iWhere == CLP_ABOVE_UPPER) {
1823	upperValue = lowerValue;
1824	lowerValue = bound_[iSequence];
1825	numberInfeasibilities_--;
1826	}
1827	// get correct place
1828	// If fixed give benefit of doubt
1829	if (lowerValue == upperValue)
1830	value = lowerValue;
1831	int newWhere = CLP_FEASIBLE;
1832	if (value - upperValue <= primalTolerance) {
1833	if (value - lowerValue >= -primalTolerance) {
1834	// feasible
1835	//newWhere=CLP_FEASIBLE;
1836	} else {
1837	// below
1838	newWhere = CLP_BELOW_LOWER;
1839	costValue -= infeasibilityWeight_;
1840	numberInfeasibilities_++;
1841	assert (fabs(lowerValue) < 1.0e100);
1842	}
1843	} else {
1844	// above
1845	newWhere = CLP_ABOVE_UPPER;
1846	costValue += infeasibilityWeight_;
1847	numberInfeasibilities_++;
1848	}
1849	if (iWhere != newWhere) {
1850	difference = cost[iSequence] - costValue;
1851	setOriginalStatus(status_[iSequence], newWhere);
1852	if (newWhere == CLP_BELOW_LOWER) {
1853	bound_[iSequence] = upperValue;
1854	upper[iSequence] = lowerValue;
1855	lower[iSequence] = -COIN_DBL_MAX;
1856	} else if (newWhere == CLP_ABOVE_UPPER) {
1857	bound_[iSequence] = lowerValue;
1858	lower[iSequence] = upperValue;
1859	upper[iSequence] = COIN_DBL_MAX;
1860	} else {
1861	lower[iSequence] = lowerValue;
1862	upper[iSequence] = upperValue;
1863	}
1864	cost[iSequence] = costValue;
1865	}
1866	// set correctly
1867	if (fabs(value - lowerValue) <= primalTolerance * 1.001) {
1868	value = CoinMin(value, lowerValue + primalTolerance);
1869	} else if (fabs(value - upperValue) <= primalTolerance * 1.001) {
1870	value = CoinMax(value, upperValue - primalTolerance);
1871	} else {
1872	//printf("*** variable wandered off bound %g %g %g!\n",
1873	// lowerValue,value,upperValue);
1874	if (value - lowerValue <= upperValue - value)
1875	value = lowerValue + primalTolerance;
1876	else
1877	value = upperValue - primalTolerance;
1878	}
1879	}
1880	changeCost_ += value * difference;
1881	return direction;
1882	}
1883	// Returns nearest bound
1884	double
1885	ClpNonLinearCost::nearest(int iSequence, double solutionValue)
1886	{
1887	assert (model_ != NULL);
1888	double nearest = 0.0;
1889	if (CLP_METHOD1) {
1890	// get where in bound sequence
1891	int iRange;
1892	int start = start_[iSequence];
1893	int end = start_[iSequence+1];
1894	int jRange = -1;
1895	nearest = COIN_DBL_MAX;
1896	for (iRange = start; iRange < end; iRange++) {
1897	if (fabs(solutionValue - lower_[iRange]) < nearest) {
1898	jRange = iRange;
1899	nearest = fabs(solutionValue - lower_[iRange]);
1900	}
1901	}
1902	assert(jRange < end);
1903	nearest = lower_[jRange];
1904	}
1905	if (CLP_METHOD2) {
1906	const double * upper = model_->upperRegion();
1907	const double * lower = model_->lowerRegion();
1908	double lowerValue = lower[iSequence];
1909	double upperValue = upper[iSequence];
1910	int iWhere = originalStatus(status_[iSequence]);
1911	if (iWhere == CLP_BELOW_LOWER) {
1912	lowerValue = upperValue;
1913	upperValue = bound_[iSequence];
1914	assert (fabs(lowerValue) < 1.0e100);
1915	} else if (iWhere == CLP_ABOVE_UPPER) {
1916	upperValue = lowerValue;
1917	lowerValue = bound_[iSequence];
1918	}
1919	if (fabs(solutionValue - lowerValue) < fabs(solutionValue - upperValue))
1920	nearest = lowerValue;
1921	else
1922	nearest = upperValue;
1923	}
1924	return nearest;
1925	}
1926	/// Feasible cost with offset and direction (i.e. for reporting)
1927	double
1928	ClpNonLinearCost::feasibleReportCost() const
1929	{
1930	double value;
1931	model_->getDblParam(ClpObjOffset, value);
1932	return (feasibleCost_ + model_->objectiveAsObject()->nonlinearOffset()) * model_->optimizationDirection() /
1933	(model_->objectiveScale() * model_->rhsScale()) - value;
1934	}
1935	// Get rid of real costs (just for moment)
1936	void
1937	ClpNonLinearCost::zapCosts()
1938	{
1939	int iSequence;
1940	double infeasibilityCost = model_->infeasibilityCost();
1941	// zero out all costs
1942	int numberTotal = numberColumns_ + numberRows_;
1943	if (CLP_METHOD1) {
1944	int n = start_[numberTotal];
1945	memset(cost_, 0, n * sizeof(double));
1946	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1947	int start = start_[iSequence];
1948	int end = start_[iSequence+1] - 1;
1949	// correct costs for this infeasibility weight
1950	if (infeasible(start)) {
1951	cost_[start] = -infeasibilityCost;
1952	}
1953	if (infeasible(end - 1)) {
1954	cost_[end-1] = infeasibilityCost;
1955	}
1956	}
1957	}
1958	if (CLP_METHOD2) {
1959	}
1960	}
1961	#ifdef VALIDATE
1962	// For debug
1963	void
1964	ClpNonLinearCost::validate()
1965	{
1966	double primalTolerance = model_->currentPrimalTolerance();
1967	int iSequence;
1968	const double * solution = model_->solutionRegion();
1969	const double * upper = model_->upperRegion();
1970	const double * lower = model_->lowerRegion();
1971	const double * cost = model_->costRegion();
1972	double infeasibilityCost = model_->infeasibilityCost();
1973	int numberTotal = numberRows_ + numberColumns_;
1974	int numberInfeasibilities = 0;
1975	double sumInfeasibilities = 0.0;
1976
1977	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1978	double value = solution[iSequence];
1979	int iStatus = status_[iSequence];
1980	assert (currentStatus(iStatus) == CLP_SAME);
1981	double lowerValue = lower[iSequence];
1982	double upperValue = upper[iSequence];
1983	double costValue = cost2_[iSequence];
1984	int iWhere = originalStatus(iStatus);
1985	if (iWhere == CLP_BELOW_LOWER) {
1986	lowerValue = upperValue;
1987	upperValue = bound_[iSequence];
1988	assert (fabs(lowerValue) < 1.0e100);
1989	costValue -= infeasibilityCost;
1990	assert (value <= lowerValue - primalTolerance);
1991	numberInfeasibilities++;
1992	sumInfeasibilities += lowerValue - value - primalTolerance;
1993	assert (model_->getStatus(iSequence) == ClpSimplex::basic);
1994	} else if (iWhere == CLP_ABOVE_UPPER) {
1995	upperValue = lowerValue;
1996	lowerValue = bound_[iSequence];
1997	costValue += infeasibilityCost;
1998	assert (value >= upperValue + primalTolerance);
1999	numberInfeasibilities++;
2000	sumInfeasibilities += value - upperValue - primalTolerance;
2001	assert (model_->getStatus(iSequence) == ClpSimplex::basic);
2002	} else {
2003	assert (value >= lowerValue - primalTolerance && value <= upperValue + primalTolerance);
2004	}
2005	assert (lowerValue == saveLowerV[iSequence]);
2006	assert (upperValue == saveUpperV[iSequence]);
2007	assert (costValue == cost[iSequence]);
2008	}
2009	if (numberInfeasibilities)
2010	printf("JJ %d infeasibilities summing to %g\n",
2011	numberInfeasibilities, sumInfeasibilities);
2012	}
2013	#endif

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

Download in other formats: