source: trunk/ClpNonLinearCost.cpp @ 162

Last change on this file since 162 was 162, checked in by forrest, 17 years ago

Piecewise linear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.6 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
7#include "CoinIndexedVector.hpp"
8
9#include "ClpNonLinearCost.hpp"
10#include "ClpSimplex.hpp"
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpNonLinearCost::ClpNonLinearCost () :
20  changeCost_(0.0),
21  largestInfeasibility_(0.0),
22  sumInfeasibilities_(0.0),
23  numberRows_(0),
24  numberColumns_(0),
25  start_(NULL),
26  whichRange_(NULL),
27  offset_(NULL),
28  lower_(NULL),
29  cost_(NULL),
30  model_(NULL),
31  infeasible_(NULL),
32  numberInfeasibilities_(-1),
33  convex_(true),
34  bothWays_(false)
35{
36
37}
38/* Constructor from simplex.
39   This will just set up wasteful arrays for linear, but
40   later may do dual analysis and even finding duplicate columns
41*/
42ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model)
43{
44  model_ = model;
45  numberRows_ = model_->numberRows();
46  numberColumns_ = model_->numberColumns();
47  int numberTotal = numberRows_+numberColumns_;
48  convex_ = true;
49  bothWays_ = false;
50  start_ = new int [numberTotal+1];
51  whichRange_ = new int [numberTotal];
52  offset_ = new int [numberTotal];
53  memset(offset_,0,numberTotal*sizeof(int));
54
55  numberInfeasibilities_=0;
56  changeCost_=0.0;
57  double infeasibilityCost = model_->infeasibilityCost();
58  sumInfeasibilities_=0.0;
59  largestInfeasibility_=0.0;
60
61  // First see how much space we need
62  int put=0;
63
64  int iSequence;
65  double * upper = model_->upperRegion();
66  double * lower = model_->lowerRegion();
67  double * cost = model_->costRegion();
68
69  for (iSequence=0;iSequence<numberTotal;iSequence++) {
70    if (upper[iSequence]<1.0e20)
71      put++;
72    put += 3;
73  }
74
75  lower_ = new double [put];
76  cost_ = new double [put];
77  infeasible_ = new unsigned int[(put+31)>>5];
78  memset(infeasible_,0,((put+31)>>5)*sizeof(unsigned int));
79
80  put=0;
81
82  start_[0]=0;
83
84  for (iSequence=0;iSequence<numberTotal;iSequence++) {
85
86    lower_[put] = -COIN_DBL_MAX;
87    setInfeasible(put,true);
88
89    cost_[put++] = cost[iSequence]-infeasibilityCost;
90    whichRange_[iSequence]=put;
91    lower_[put] = lower[iSequence];
92    cost_[put++] = cost[iSequence];
93    lower_[put] = upper[iSequence];
94    cost_[put++] = cost[iSequence]+infeasibilityCost;
95    if (upper[iSequence]<1.0e20) {
96      lower_[put] = COIN_DBL_MAX;
97      setInfeasible(put-1,true);
98      cost_[put++] = 1.0e50;
99    }
100    start_[iSequence+1]=put;
101  }
102
103}
104ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model,const int * starts,
105                   const double * lowerNon, const double * costNon)
106{
107
108  // what about scaling? - only try without it initially
109  assert(!model->scalingFlag());
110  model_ = model;
111  numberRows_ = model_->numberRows();
112  numberColumns_ = model_->numberColumns();
113  int numberTotal = numberRows_+numberColumns_;
114  convex_ = true;
115  bothWays_ = true;
116  start_ = new int [numberTotal+1];
117  whichRange_ = new int [numberTotal];
118  offset_ = new int [numberTotal];
119  memset(offset_,0,numberTotal*sizeof(int));
120 
121  double whichWay = model_->optimizationDirection();
122  printf("Direction %g\n",whichWay);
123
124  numberInfeasibilities_=0;
125  changeCost_=0.0;
126  double infeasibilityCost = model_->infeasibilityCost();
127  largestInfeasibility_=0.0;
128  sumInfeasibilities_=0.0;
129
130  int iSequence;
131  assert (!model_->rowObjective());
132  double * cost = model_->objective();
133
134  // First see how much space we need
135  // Also set up feasible bounds
136  int put=starts[numberColumns_];
137
138  double * columnUpper = model_->columnUpper();
139  double * columnLower = model_->columnLower();
140  for (iSequence=0;iSequence<numberColumns_;iSequence++) {
141    if (columnLower[iSequence]>-1.0e20)
142      put++;
143    if (columnUpper[iSequence]<1.0e20)
144      put++;
145  }
146
147  double * rowUpper = model_->rowUpper();
148  double * rowLower = model_->rowLower();
149  for (iSequence=0;iSequence<numberRows_;iSequence++) {
150    if (rowLower[iSequence]>-1.0e20)
151      put++;
152    if (rowUpper[iSequence]<1.0e20)
153      put++;
154    put +=2;
155  }
156  lower_ = new double [put];
157  cost_ = new double [put];
158  infeasible_ = new unsigned int[(put+31)>>5];
159  memset(infeasible_,0,((put+31)>>5)*sizeof(unsigned int));
160
161  // now fill in
162  put=0;
163
164  start_[0]=0;
165  for (iSequence=0;iSequence<numberTotal;iSequence++) {
166    lower_[put] = -COIN_DBL_MAX;
167    whichRange_[iSequence]=put+1;
168    double thisCost;
169    double lowerValue;
170    double upperValue;
171    if (iSequence>=numberColumns_) {
172      // rows
173      lowerValue = rowLower[iSequence-numberColumns_];
174      upperValue = rowUpper[iSequence-numberColumns_];
175      if (lowerValue>-1.0e30) {
176        setInfeasible(put,true);
177        cost_[put++] = -infeasibilityCost;
178        lower_[put] = lowerValue;
179      }
180      cost_[put++] = 0.0;
181      thisCost = 0.0;
182    } else {
183      // columns - move costs and see if convex
184      lowerValue = columnLower[iSequence];
185      upperValue = columnUpper[iSequence];
186      if (lowerValue>-1.0e30) {
187        setInfeasible(put,true);
188        cost_[put++] = whichWay*cost[iSequence]-infeasibilityCost;
189        lower_[put] = lowerValue;
190      }
191      int iIndex = starts[iSequence];
192      int end = starts[iSequence+1];
193      assert (fabs(columnLower[iSequence]-lowerNon[iIndex])<1.0e-8);
194      thisCost = -COIN_DBL_MAX;
195      for (;iIndex<end;iIndex++) {
196        if (lowerNon[iIndex]<columnUpper[iSequence]-1.0e-8) {
197          lower_[put] = lowerNon[iIndex];
198          cost_[put++] = whichWay*costNon[iIndex];
199          // check convexity
200          if (whichWay*costNon[iIndex]<thisCost-1.0e-12)
201            convex_ = false;
202          thisCost = whichWay*costNon[iIndex];
203        } else {
204          break;
205        }
206      }
207    }
208    lower_[put] = upperValue;
209    setInfeasible(put,true);
210    cost_[put++] = thisCost+infeasibilityCost;
211    if (upperValue<1.0e20) {
212      lower_[put] = COIN_DBL_MAX;
213      cost_[put++] = 1.0e50;
214    }
215    int iFirst = start_[iSequence];
216    if (lower_[iFirst] != -COIN_DBL_MAX) {
217      setInfeasible(iFirst,true);
218      whichRange_[iSequence]=iFirst+1;
219    } else {
220      whichRange_[iSequence]=iFirst;
221    }
222    start_[iSequence+1]=put;
223  }
224  // can't handle non-convex at present
225  assert(convex_);
226}
227
228//-------------------------------------------------------------------
229// Copy constructor
230//-------------------------------------------------------------------
231ClpNonLinearCost::ClpNonLinearCost (const ClpNonLinearCost & rhs) :
232  changeCost_(0.0),
233  largestInfeasibility_(0.0),
234  sumInfeasibilities_(0.0),
235  numberRows_(rhs.numberRows_),
236  numberColumns_(rhs.numberColumns_),
237  start_(NULL),
238  whichRange_(NULL),
239  offset_(NULL),
240  lower_(NULL),
241  cost_(NULL),
242  model_(NULL),
243  infeasible_(NULL),
244  numberInfeasibilities_(-1),
245  convex_(true),
246  bothWays_(rhs.bothWays_)
247{ 
248  if (numberRows_) {
249    int numberTotal = numberRows_+numberColumns_;
250    start_ = new int [numberTotal+1];
251    memcpy(start_,rhs.start_,(numberTotal+1)*sizeof(int));
252    whichRange_ = new int [numberTotal];
253    memcpy(whichRange_,rhs.whichRange_,numberTotal*sizeof(int));
254    offset_ = new int [numberTotal];
255    memcpy(offset_,rhs.offset_,numberTotal*sizeof(int));
256    int numberEntries = start_[numberTotal];
257    lower_ = new double [numberEntries];
258    memcpy(lower_,rhs.lower_,numberEntries*sizeof(double));
259    cost_ = new double [numberEntries];
260    memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
261    model_ = rhs.model_;
262    numberInfeasibilities_=rhs.numberInfeasibilities_;
263    changeCost_ = rhs.changeCost_;
264    largestInfeasibility_ = rhs.largestInfeasibility_;
265    sumInfeasibilities_ = rhs.sumInfeasibilities_;
266    convex_ = rhs.convex_;
267    infeasible_ = new unsigned int[(numberEntries+31)>>5];
268    memcpy(infeasible_,rhs.infeasible_,
269           ((numberEntries+31)>>5)*sizeof(unsigned int));
270  }
271}
272
273//-------------------------------------------------------------------
274// Destructor
275//-------------------------------------------------------------------
276ClpNonLinearCost::~ClpNonLinearCost ()
277{
278  delete [] start_;
279  delete [] whichRange_;
280  delete [] offset_;
281  delete [] lower_;
282  delete [] cost_;
283  delete [] infeasible_;
284}
285
286//----------------------------------------------------------------
287// Assignment operator
288//-------------------------------------------------------------------
289ClpNonLinearCost &
290ClpNonLinearCost::operator=(const ClpNonLinearCost& rhs)
291{
292  if (this != &rhs) {
293    numberRows_ = rhs.numberRows_;
294    numberColumns_ = rhs.numberColumns_;
295    delete [] start_;
296    delete [] whichRange_;
297    delete [] offset_;
298    delete [] lower_;
299    delete []cost_;
300    delete [] infeasible_;
301    start_ = NULL;
302    whichRange_ = NULL;
303    lower_ = NULL;
304    cost_= NULL;
305    infeasible_=NULL;
306    if (numberRows_) {
307      int numberTotal = numberRows_+numberColumns_;
308      start_ = new int [numberTotal+1];
309      memcpy(start_,rhs.start_,(numberTotal+1)*sizeof(int));
310      whichRange_ = new int [numberTotal];
311      memcpy(whichRange_,rhs.whichRange_,numberTotal*sizeof(int));
312      offset_ = new int [numberTotal];
313      memcpy(offset_,rhs.offset_,numberTotal*sizeof(int));
314      int numberEntries = start_[numberTotal];
315      lower_ = new double [numberEntries];
316      memcpy(lower_,rhs.lower_,numberEntries*sizeof(double));
317      cost_ = new double [numberEntries];
318      memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
319      infeasible_ = new unsigned int[(numberEntries+31)>>5];
320      memcpy(infeasible_,rhs.infeasible_,
321             ((numberEntries+31)>>5)*sizeof(unsigned int));
322    }
323    model_ = rhs.model_;
324    numberInfeasibilities_=rhs.numberInfeasibilities_;
325    changeCost_ = rhs.changeCost_;
326    largestInfeasibility_ = rhs.largestInfeasibility_;
327    sumInfeasibilities_ = rhs.sumInfeasibilities_;
328    convex_ = rhs.convex_;
329    bothWays_ = rhs.bothWays_;
330  }
331  return *this;
332}
333
334// Changes infeasible costs and computes number and cost of infeas
335// We will need to re-think objective offsets later
336// We will also need a 2 bit per variable array for some
337// purpose which will come to me later
338void 
339ClpNonLinearCost::checkInfeasibilities(bool toNearest)
340{
341  numberInfeasibilities_=0;
342  double infeasibilityCost = model_->infeasibilityCost();
343  changeCost_=0.0;
344  largestInfeasibility_ = 0.0;
345  sumInfeasibilities_ = 0.0;
346  double primalTolerance = model_->currentPrimalTolerance();
347 
348  int iSequence;
349  double * solution = model_->solutionRegion();
350  double * upper = model_->upperRegion();
351  double * lower = model_->lowerRegion();
352  double * cost = model_->costRegion();
353   
354  // nonbasic should be at a valid bound
355  for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
356    double lowerValue;
357    double upperValue;
358    double value=solution[iSequence];
359    int iRange;
360    // get correct place
361    int start = start_[iSequence];
362    int end = start_[iSequence+1]-1;
363    // correct costs for this infeasibility weight
364    if (infeasible(start))
365      cost_[start] = cost_[start+1]-infeasibilityCost;
366    if (infeasible(end-1))
367      cost_[end-1] = cost_[end-2]+infeasibilityCost;
368    for (iRange=start; iRange<end;iRange++) {
369      if (value<lower_[iRange+1]+primalTolerance) {
370        // put in better range if infeasible
371        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
372          iRange++;
373        whichRange_[iSequence]=iRange;
374        break;
375      }
376    }
377    assert(iRange<end);
378    lowerValue = lower_[iRange];
379    upperValue = lower_[iRange+1];
380    ClpSimplex::Status status = model_->getStatus(iSequence);
381    if (upperValue==lowerValue) {
382      if (status != ClpSimplex::basic) 
383        model_->setStatus(iSequence,ClpSimplex::isFixed);
384    }
385    switch(status) {
386     
387    case ClpSimplex::basic:
388    case ClpSimplex::superBasic:
389      // iRange is in correct place
390      // slot in here
391      if (infeasible(iRange)) {
392        if (lower_[iRange]<-1.0e50) {
393          //cost_[iRange] = cost_[iRange+1]-infeasibilityCost;
394          // possibly below
395          lowerValue = lower_[iRange+1];
396          if (value<lowerValue-primalTolerance) {
397            value = lowerValue-value;
398            sumInfeasibilities_ += value;
399            largestInfeasibility_ = max(largestInfeasibility_,value);
400            changeCost_ -= lowerValue*
401              (cost_[iRange]-cost[iSequence]);
402            numberInfeasibilities_++;
403          }
404        } else {
405          //cost_[iRange] = cost_[iRange-1]+infeasibilityCost;
406          // possibly above
407          upperValue = lower_[iRange];
408          if (value>upperValue+primalTolerance) {
409            value = value-upperValue;
410            sumInfeasibilities_ += value;
411            largestInfeasibility_ = max(largestInfeasibility_,value);
412            changeCost_ -= upperValue*
413              (cost_[iRange]-cost[iSequence]);
414            numberInfeasibilities_++;
415          }
416        }
417      }
418      //lower[iSequence] = lower_[iRange];
419      //upper[iSequence] = lower_[iRange+1];
420      //cost[iSequence] = cost_[iRange];
421      break;
422    case ClpSimplex::isFree:
423      if (toNearest)
424        solution[iSequence] = 0.0;
425      break;
426    case ClpSimplex::atUpperBound:
427      if (!toNearest) {
428        // With increasing tolerances - we may be at wrong place
429        if (fabs(value-upperValue)>primalTolerance*1.0001) {
430          assert(fabs(value-lowerValue)<=primalTolerance*1.0001); 
431          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
432        }
433      } else {
434        if (fabs(value-upperValue)<=fabs(value-lowerValue)) {
435          solution[iSequence] = upperValue;
436        } else {
437          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
438          solution[iSequence] = lowerValue;
439        }
440      }
441      break;
442    case ClpSimplex::atLowerBound:
443      if (!toNearest) {
444        // With increasing tolerances - we may be at wrong place
445        if (fabs(value-lowerValue)>primalTolerance*1.0001) {
446          assert(fabs(value-upperValue)<=primalTolerance*1.0001); 
447          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
448        }
449      } else {
450        if (fabs(value-lowerValue)<=fabs(value-upperValue)) {
451          solution[iSequence] = lowerValue;
452        } else {
453          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
454          solution[iSequence] = upperValue;
455        }
456      }
457      break;
458    case ClpSimplex::isFixed:
459      break;
460    }
461    lower[iSequence] = lower_[iRange];
462    upper[iSequence] = lower_[iRange+1];
463    cost[iSequence] = cost_[iRange];
464  }
465}
466/* Goes through one bound for each variable.
467   If array[i]*multiplier>0 goes down, otherwise up.
468   The indices are row indices and need converting to sequences
469*/
470void 
471ClpNonLinearCost::goThru(int numberInArray, double multiplier,
472              const int * index, const double * array,
473                         double * rhs)
474{
475  assert (model_!=NULL);
476  const int * pivotVariable = model_->pivotVariable();
477  int i;
478  for (i=0;i<numberInArray;i++) {
479    int iRow = index[i];
480    int iPivot = pivotVariable[iRow];
481    double alpha = multiplier*array[iRow];
482    // get where in bound sequence
483    int iRange = whichRange_[iPivot];
484    iRange += offset_[iPivot]; //add temporary bias
485    double value = model_->solution(iPivot);
486    if (alpha>0.0) {
487      // down one
488      iRange--;
489      assert(iRange>=start_[iPivot]);
490      rhs[iRow] = value - lower_[iRange];
491    } else {
492      // up one
493      iRange++;
494      assert(iRange<start_[iPivot+1]-1);
495      rhs[iRow] = lower_[iRange+1] - value;
496    }
497    offset_[iPivot] = iRange - whichRange_[iPivot];
498  }
499}
500/* Takes off last iteration (i.e. offsets closer to 0)
501 */
502void 
503ClpNonLinearCost::goBack(int numberInArray, const int * index, 
504              double * rhs)
505{
506  assert (model_!=NULL);
507  const int * pivotVariable = model_->pivotVariable();
508  int i;
509  for (i=0;i<numberInArray;i++) {
510    int iRow = index[i];
511    int iPivot = pivotVariable[iRow];
512    // get where in bound sequence
513    int iRange = whichRange_[iPivot];
514    bool down;
515    // get closer to original
516    if (offset_[iPivot]>0) {
517      offset_[iPivot]--;
518      assert (offset_[iPivot]>=0);
519      down = false;
520    } else {
521      offset_[iPivot]++;
522      assert (offset_[iPivot]<=0);
523      down = true;
524    }
525    iRange += offset_[iPivot]; //add temporary bias
526    double value = model_->solution(iPivot);
527    if (down) {
528      // down one
529      assert(iRange>=start_[iPivot]);
530      rhs[iRow] = value - lower_[iRange];
531    } else {
532      // up one
533      assert(iRange<start_[iPivot+1]-1);
534      rhs[iRow] = lower_[iRange+1] - value;
535    }
536  }
537}
538void 
539ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
540{
541  assert (model_!=NULL);
542  const int * pivotVariable = model_->pivotVariable();
543  int i;
544  int number = update->getNumElements();
545  const int * index = update->getIndices();
546  for (i=0;i<number;i++) {
547    int iRow = index[i];
548    int iPivot = pivotVariable[iRow];
549    offset_[iPivot]=0;
550  }
551#ifdef CLP_DEBUG
552  for (i=0;i<numberRows_+numberColumns_;i++) 
553    assert(!offset_[i]);
554#endif
555}
556void 
557ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
558{
559  assert (model_!=NULL);
560  double primalTolerance = model_->currentPrimalTolerance();
561  const int * pivotVariable = model_->pivotVariable();
562  int i;
563  for (i=0;i<numberInArray;i++) {
564    int iRow = index[i];
565    int iPivot = pivotVariable[iRow];
566    // get where in bound sequence
567    int iRange;
568    double value = model_->solution(iPivot);
569    int start = start_[iPivot];
570    int end = start_[iPivot+1]-1;
571    for (iRange=start; iRange<end;iRange++) {
572      if (value<lower_[iRange+1]+primalTolerance) {
573        // put in better range
574        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
575          iRange++;
576        break;
577      }
578    }
579    assert(iRange<end);
580    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
581    double & lower = model_->lowerAddress(iPivot);
582    double & upper = model_->upperAddress(iPivot);
583    double & cost = model_->costAddress(iPivot);
584    whichRange_[iPivot]=iRange;
585    lower = lower_[iRange];
586    upper = lower_[iRange+1];
587    cost = cost_[iRange];
588  }
589}
590/* Puts back correct infeasible costs for each variable
591   The input indices are row indices and need converting to sequences
592   for costs.
593   On input array is empty (but indices exist).  On exit just
594   changed costs will be stored as normal CoinIndexedVector
595*/
596void 
597ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
598{
599  assert (model_!=NULL);
600  double primalTolerance = model_->currentPrimalTolerance();
601  const int * pivotVariable = model_->pivotVariable();
602  int number=0;
603  int * index = update->getIndices();
604  double * work = update->denseVector();
605  int i;
606  for (i=0;i<numberInArray;i++) {
607    int iRow = index[i];
608    int iPivot = pivotVariable[iRow];
609    // get where in bound sequence
610    int iRange;
611    double value = model_->solution(iPivot);
612    int start = start_[iPivot];
613    int end = start_[iPivot+1]-1;
614    for (iRange=start; iRange<end;iRange++) {
615      if (value<lower_[iRange+1]+primalTolerance) {
616        // put in better range
617        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
618          iRange++;
619        break;
620      }
621    }
622    assert(iRange<end);
623    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
624    int jRange = whichRange_[iPivot];
625    if (iRange!=jRange) {
626      // changed
627      work[iRow] = cost_[jRange]-cost_[iRange];
628      index[number++]=iRow;
629      double & lower = model_->lowerAddress(iPivot);
630      double & upper = model_->upperAddress(iPivot);
631      double & cost = model_->costAddress(iPivot);
632      whichRange_[iPivot]=iRange;
633      lower = lower_[iRange];
634      upper = lower_[iRange+1];
635      cost = cost_[iRange];
636    }
637  }
638  update->setNumElements(number);
639}
640/* Sets bounds and cost for one variable - returns change in cost*/
641double 
642ClpNonLinearCost::setOne(int iPivot, double value)
643{
644  assert (model_!=NULL);
645  double primalTolerance = model_->currentPrimalTolerance();
646  // get where in bound sequence
647  int iRange;
648  int start = start_[iPivot];
649  int end = start_[iPivot+1]-1;
650  if (!bothWays_) {
651    for (iRange=start; iRange<end;iRange++) {
652      if (value<lower_[iRange+1]+primalTolerance) {
653        // put in better range
654        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
655          iRange++;
656        break;
657      }
658    }
659  } else {
660    // leave in current if possible
661    iRange = whichRange_[iPivot];
662    if (value<lower_[iRange]-primalTolerance||value>lower_[iRange+1]+primalTolerance) {
663      for (iRange=start; iRange<end;iRange++) {
664        if (value<lower_[iRange+1]+primalTolerance) {
665          // put in better range
666          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
667            iRange++;
668          break;
669        }
670      }
671    }
672  }
673  assert(iRange<end);
674  whichRange_[iPivot]=iRange;
675  double & lower = model_->lowerAddress(iPivot);
676  double & upper = model_->upperAddress(iPivot);
677  double & cost = model_->costAddress(iPivot);
678  lower = lower_[iRange];
679  upper = lower_[iRange+1];
680  ClpSimplex::Status status = model_->getStatus(iPivot);
681  if (upper==lower) {
682    if (status != ClpSimplex::basic) {
683      model_->setStatus(iPivot,ClpSimplex::isFixed);
684      status = ClpSimplex::basic; // so will skip
685    }
686  }
687  switch(status) {
688     
689  case ClpSimplex::basic:
690  case ClpSimplex::superBasic:
691  case ClpSimplex::isFree:
692    break;
693  case ClpSimplex::atUpperBound:
694  case ClpSimplex::atLowerBound:
695  case ClpSimplex::isFixed:
696    // set correctly
697    if (fabs(value-lower)<=primalTolerance*1.001){
698      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
699    } else if (fabs(value-upper)<=primalTolerance*1.001){
700      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
701    } else {
702      // set superBasic
703      model_->setStatus(iPivot,ClpSimplex::superBasic);
704    }
705    break;
706  }
707  double difference = cost-cost_[iRange]; 
708  cost = cost_[iRange];
709  changeCost_ += value*difference;
710  return difference;
711}
712// Returns nearest bound
713double 
714ClpNonLinearCost::nearest(int sequence, double solutionValue)
715{
716  assert (model_!=NULL);
717  // get where in bound sequence
718  int iRange;
719  int start = start_[sequence];
720  int end = start_[sequence+1];
721  int jRange=-1;
722  double nearest=COIN_DBL_MAX;
723  for (iRange=start; iRange<end;iRange++) {
724    if (fabs(solutionValue-lower_[iRange])<nearest) {
725      jRange=iRange;
726      nearest=fabs(solutionValue-lower_[iRange]);
727    }
728  }
729  assert(jRange<end);
730  return lower_[jRange];
731}
732
Note: See TracBrowser for help on using the repository browser.