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

Last change on this file since 1769 was 1769, checked in by forrest, 10 years ago

changes for advanced use of Clp

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