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

Last change on this file since 2030 was 2030, checked in by forrest, 5 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

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