source: stable/1.6/Clp/src/ClpNonLinearCost.cpp @ 1609

Last change on this file since 1609 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

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