source: branches/pre/ClpNonLinearCost.cpp @ 192

Last change on this file since 192 was 192, checked in by forrest, 18 years ago

For Yiming

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