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

Last change on this file since 1676 was 1676, checked in by forrest, 8 years ago

For GUB

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 80.1 KB
Line 
1/* $Id: ClpNonLinearCost.cpp 1676 2011-01-20 17:00:28Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include <iostream>
8#include <cassert>
9
10#include "CoinIndexedVector.hpp"
11
12#include "ClpSimplex.hpp"
13#include "CoinHelperFunctions.hpp"
14#include "ClpNonLinearCost.hpp"
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
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)
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     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
Note: See TracBrowser for help on using the repository browser.