source: branches/pre/ClpNonLinearCost.cpp @ 180

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

new stuff

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