source: trunk/Clp/src/ClpNonLinearCost.cpp

Last change on this file was 2385, checked in by unxusr, 9 months ago

formatting

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