source: trunk/ClpNonLinearCost.cpp @ 193

Last change on this file since 193 was 193, checked in by forrest, 16 years ago

For Yiming

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.2 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(bool toNearest)
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 
375  int iSequence;
376  double * solution = model_->solutionRegion();
377  double * upper = model_->upperRegion();
378  double * lower = model_->lowerRegion();
379  double * cost = model_->costRegion();
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)>primalTolerance*1.0001) {
457          assert(fabs(value-lowerValue)<=primalTolerance*1.0001); 
458          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
459        }
460      } else {
461        if (fabs(value-upperValue)<=fabs(value-lowerValue)) {
462          solution[iSequence] = upperValue;
463        } else {
464          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
465          solution[iSequence] = lowerValue;
466        }
467      }
468      break;
469    case ClpSimplex::atLowerBound:
470      if (!toNearest) {
471        // With increasing tolerances - we may be at wrong place
472        if (fabs(value-lowerValue)>primalTolerance*1.0001) {
473          assert(fabs(value-upperValue)<=primalTolerance*1.0001); 
474          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
475        }
476      } else {
477        if (fabs(value-lowerValue)<=fabs(value-upperValue)) {
478          solution[iSequence] = lowerValue;
479        } else {
480          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
481          solution[iSequence] = upperValue;
482        }
483      }
484      break;
485    case ClpSimplex::isFixed:
486      break;
487    }
488    lower[iSequence] = lower_[iRange];
489    upper[iSequence] = lower_[iRange+1];
490    cost[iSequence] = cost_[iRange];
491  }
492}
493/* Goes through one bound for each variable.
494   If array[i]*multiplier>0 goes down, otherwise up.
495   The indices are row indices and need converting to sequences
496*/
497void 
498ClpNonLinearCost::goThru(int numberInArray, double multiplier,
499              const int * index, const double * array,
500                         double * rhs)
501{
502  assert (model_!=NULL);
503  const int * pivotVariable = model_->pivotVariable();
504  int i;
505  for (i=0;i<numberInArray;i++) {
506    int iRow = index[i];
507    int iPivot = pivotVariable[iRow];
508    double alpha = multiplier*array[iRow];
509    // get where in bound sequence
510    int iRange = whichRange_[iPivot];
511    iRange += offset_[iPivot]; //add temporary bias
512    double value = model_->solution(iPivot);
513    if (alpha>0.0) {
514      // down one
515      iRange--;
516      assert(iRange>=start_[iPivot]);
517      rhs[iRow] = value - lower_[iRange];
518    } else {
519      // up one
520      iRange++;
521      assert(iRange<start_[iPivot+1]-1);
522      rhs[iRow] = lower_[iRange+1] - value;
523    }
524    offset_[iPivot] = iRange - whichRange_[iPivot];
525  }
526}
527/* Takes off last iteration (i.e. offsets closer to 0)
528 */
529void 
530ClpNonLinearCost::goBack(int numberInArray, const int * index, 
531              double * rhs)
532{
533  assert (model_!=NULL);
534  const int * pivotVariable = model_->pivotVariable();
535  int i;
536  for (i=0;i<numberInArray;i++) {
537    int iRow = index[i];
538    int iPivot = pivotVariable[iRow];
539    // get where in bound sequence
540    int iRange = whichRange_[iPivot];
541    bool down;
542    // get closer to original
543    if (offset_[iPivot]>0) {
544      offset_[iPivot]--;
545      assert (offset_[iPivot]>=0);
546      down = false;
547    } else {
548      offset_[iPivot]++;
549      assert (offset_[iPivot]<=0);
550      down = true;
551    }
552    iRange += offset_[iPivot]; //add temporary bias
553    double value = model_->solution(iPivot);
554    if (down) {
555      // down one
556      assert(iRange>=start_[iPivot]);
557      rhs[iRow] = value - lower_[iRange];
558    } else {
559      // up one
560      assert(iRange<start_[iPivot+1]-1);
561      rhs[iRow] = lower_[iRange+1] - value;
562    }
563  }
564}
565void 
566ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
567{
568  assert (model_!=NULL);
569  const int * pivotVariable = model_->pivotVariable();
570  int i;
571  int number = update->getNumElements();
572  const int * index = update->getIndices();
573  for (i=0;i<number;i++) {
574    int iRow = index[i];
575    int iPivot = pivotVariable[iRow];
576    offset_[iPivot]=0;
577  }
578#ifdef CLP_DEBUG
579  for (i=0;i<numberRows_+numberColumns_;i++) 
580    assert(!offset_[i]);
581#endif
582}
583void 
584ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
585{
586  assert (model_!=NULL);
587  double primalTolerance = model_->currentPrimalTolerance();
588  const int * pivotVariable = model_->pivotVariable();
589  int i;
590  for (i=0;i<numberInArray;i++) {
591    int iRow = index[i];
592    int iPivot = pivotVariable[iRow];
593    // get where in bound sequence
594    int iRange;
595    int currentRange = whichRange_[iPivot];
596    double value = model_->solution(iPivot);
597    int start = start_[iPivot];
598    int end = start_[iPivot+1]-1;
599    for (iRange=start; iRange<end;iRange++) {
600      if (value<lower_[iRange+1]+primalTolerance) {
601        // put in better range
602        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
603          iRange++;
604        break;
605      }
606    }
607    assert(iRange<end);
608    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
609    double & lower = model_->lowerAddress(iPivot);
610    double & upper = model_->upperAddress(iPivot);
611    double & cost = model_->costAddress(iPivot);
612    whichRange_[iPivot]=iRange;
613    if (iRange!=currentRange) {
614      if (infeasible(iRange))
615        numberInfeasibilities_++;
616      if (infeasible(currentRange))
617        numberInfeasibilities_--;
618    }
619    lower = lower_[iRange];
620    upper = lower_[iRange+1];
621    cost = cost_[iRange];
622  }
623}
624/* Puts back correct infeasible costs for each variable
625   The input indices are row indices and need converting to sequences
626   for costs.
627   On input array is empty (but indices exist).  On exit just
628   changed costs will be stored as normal CoinIndexedVector
629*/
630void 
631ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
632{
633  assert (model_!=NULL);
634  double primalTolerance = model_->currentPrimalTolerance();
635  const int * pivotVariable = model_->pivotVariable();
636  int number=0;
637  int * index = update->getIndices();
638  double * work = update->denseVector();
639  int i;
640  for (i=0;i<numberInArray;i++) {
641    int iRow = index[i];
642    int iPivot = pivotVariable[iRow];
643    // get where in bound sequence
644    int iRange;
645    double value = model_->solution(iPivot);
646    int start = start_[iPivot];
647    int end = start_[iPivot+1]-1;
648    for (iRange=start; iRange<end;iRange++) {
649      if (value<lower_[iRange+1]+primalTolerance) {
650        // put in better range
651        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
652          iRange++;
653        break;
654      }
655    }
656    assert(iRange<end);
657    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
658    int jRange = whichRange_[iPivot];
659    if (iRange!=jRange) {
660      // changed
661      work[iRow] = cost_[jRange]-cost_[iRange];
662      index[number++]=iRow;
663      double & lower = model_->lowerAddress(iPivot);
664      double & upper = model_->upperAddress(iPivot);
665      double & cost = model_->costAddress(iPivot);
666      whichRange_[iPivot]=iRange;
667      if (infeasible(iRange))
668        numberInfeasibilities_++;
669      if (infeasible(jRange))
670        numberInfeasibilities_--;
671      lower = lower_[iRange];
672      upper = lower_[iRange+1];
673      cost = cost_[iRange];
674    }
675  }
676  update->setNumElements(number);
677}
678/* Sets bounds and cost for one variable - returns change in cost*/
679double 
680ClpNonLinearCost::setOne(int iPivot, double value)
681{
682  assert (model_!=NULL);
683  double primalTolerance = model_->currentPrimalTolerance();
684  // get where in bound sequence
685  int iRange;
686  int currentRange = whichRange_[iPivot];
687  int start = start_[iPivot];
688  int end = start_[iPivot+1]-1;
689  if (!bothWays_) {
690    for (iRange=start; iRange<end;iRange++) {
691      if (value<lower_[iRange+1]+primalTolerance) {
692        // put in better range
693        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
694          iRange++;
695        break;
696      }
697    }
698  } else {
699    // leave in current if possible
700    iRange = whichRange_[iPivot];
701    if (value<lower_[iRange]-primalTolerance||value>lower_[iRange+1]+primalTolerance) {
702      for (iRange=start; iRange<end;iRange++) {
703        if (value<lower_[iRange+1]+primalTolerance) {
704          // put in better range
705          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
706            iRange++;
707          break;
708        }
709      }
710    }
711  }
712  assert(iRange<end);
713  whichRange_[iPivot]=iRange;
714  if (iRange!=currentRange) {
715    if (infeasible(iRange))
716      numberInfeasibilities_++;
717    if (infeasible(currentRange))
718      numberInfeasibilities_--;
719  }
720  double & lower = model_->lowerAddress(iPivot);
721  double & upper = model_->upperAddress(iPivot);
722  double & cost = model_->costAddress(iPivot);
723  lower = lower_[iRange];
724  upper = lower_[iRange+1];
725  ClpSimplex::Status status = model_->getStatus(iPivot);
726  if (upper==lower) {
727    if (status != ClpSimplex::basic) {
728      model_->setStatus(iPivot,ClpSimplex::isFixed);
729      status = ClpSimplex::basic; // so will skip
730    }
731  }
732  switch(status) {
733     
734  case ClpSimplex::basic:
735  case ClpSimplex::superBasic:
736  case ClpSimplex::isFree:
737    break;
738  case ClpSimplex::atUpperBound:
739  case ClpSimplex::atLowerBound:
740  case ClpSimplex::isFixed:
741    // set correctly
742    if (fabs(value-lower)<=primalTolerance*1.001){
743      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
744    } else if (fabs(value-upper)<=primalTolerance*1.001){
745      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
746    } else {
747      // set superBasic
748      model_->setStatus(iPivot,ClpSimplex::superBasic);
749    }
750    break;
751  }
752  double difference = cost-cost_[iRange]; 
753  cost = cost_[iRange];
754  changeCost_ += value*difference;
755  return difference;
756}
757// Returns nearest bound
758double 
759ClpNonLinearCost::nearest(int sequence, double solutionValue)
760{
761  assert (model_!=NULL);
762  // get where in bound sequence
763  int iRange;
764  int start = start_[sequence];
765  int end = start_[sequence+1];
766  int jRange=-1;
767  double nearest=COIN_DBL_MAX;
768  for (iRange=start; iRange<end;iRange++) {
769    if (fabs(solutionValue-lower_[iRange])<nearest) {
770      jRange=iRange;
771      nearest=fabs(solutionValue-lower_[iRange]);
772    }
773  }
774  assert(jRange<end);
775  return lower_[jRange];
776}
777// Sets inside bounds (i.e. non infinite - used in QP
778void 
779ClpNonLinearCost::setBounds(int sequence, double lower, double upper)
780{
781  int start = start_[sequence];
782  assert(start_[sequence+1]==start+4);
783  lower_[start+1]=lower;
784  lower_[start+2]=upper;
785}
786
Note: See TracBrowser for help on using the repository browser.