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