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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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