source: trunk/ClpNonLinearCost.cpp @ 170

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

Forget to update include/ClpModel.hpp

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.7 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    double value = model_->solution(iPivot);
596    int start = start_[iPivot];
597    int end = start_[iPivot+1]-1;
598    for (iRange=start; iRange<end;iRange++) {
599      if (value<lower_[iRange+1]+primalTolerance) {
600        // put in better range
601        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
602          iRange++;
603        break;
604      }
605    }
606    assert(iRange<end);
607    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
608    double & lower = model_->lowerAddress(iPivot);
609    double & upper = model_->upperAddress(iPivot);
610    double & cost = model_->costAddress(iPivot);
611    whichRange_[iPivot]=iRange;
612    lower = lower_[iRange];
613    upper = lower_[iRange+1];
614    cost = cost_[iRange];
615  }
616}
617/* Puts back correct infeasible costs for each variable
618   The input indices are row indices and need converting to sequences
619   for costs.
620   On input array is empty (but indices exist).  On exit just
621   changed costs will be stored as normal CoinIndexedVector
622*/
623void 
624ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
625{
626  assert (model_!=NULL);
627  double primalTolerance = model_->currentPrimalTolerance();
628  const int * pivotVariable = model_->pivotVariable();
629  int number=0;
630  int * index = update->getIndices();
631  double * work = update->denseVector();
632  int i;
633  for (i=0;i<numberInArray;i++) {
634    int iRow = index[i];
635    int iPivot = pivotVariable[iRow];
636    // get where in bound sequence
637    int iRange;
638    double value = model_->solution(iPivot);
639    int start = start_[iPivot];
640    int end = start_[iPivot+1]-1;
641    for (iRange=start; iRange<end;iRange++) {
642      if (value<lower_[iRange+1]+primalTolerance) {
643        // put in better range
644        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
645          iRange++;
646        break;
647      }
648    }
649    assert(iRange<end);
650    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
651    int jRange = whichRange_[iPivot];
652    if (iRange!=jRange) {
653      // changed
654      work[iRow] = cost_[jRange]-cost_[iRange];
655      index[number++]=iRow;
656      double & lower = model_->lowerAddress(iPivot);
657      double & upper = model_->upperAddress(iPivot);
658      double & cost = model_->costAddress(iPivot);
659      whichRange_[iPivot]=iRange;
660      lower = lower_[iRange];
661      upper = lower_[iRange+1];
662      cost = cost_[iRange];
663    }
664  }
665  update->setNumElements(number);
666}
667/* Sets bounds and cost for one variable - returns change in cost*/
668double 
669ClpNonLinearCost::setOne(int iPivot, double value)
670{
671  assert (model_!=NULL);
672  double primalTolerance = model_->currentPrimalTolerance();
673  // get where in bound sequence
674  int iRange;
675  int start = start_[iPivot];
676  int end = start_[iPivot+1]-1;
677  if (!bothWays_) {
678    for (iRange=start; iRange<end;iRange++) {
679      if (value<lower_[iRange+1]+primalTolerance) {
680        // put in better range
681        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
682          iRange++;
683        break;
684      }
685    }
686  } else {
687    // leave in current if possible
688    iRange = whichRange_[iPivot];
689    if (value<lower_[iRange]-primalTolerance||value>lower_[iRange+1]+primalTolerance) {
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    }
699  }
700  assert(iRange<end);
701  whichRange_[iPivot]=iRange;
702  double & lower = model_->lowerAddress(iPivot);
703  double & upper = model_->upperAddress(iPivot);
704  double & cost = model_->costAddress(iPivot);
705  lower = lower_[iRange];
706  upper = lower_[iRange+1];
707  ClpSimplex::Status status = model_->getStatus(iPivot);
708  if (upper==lower) {
709    if (status != ClpSimplex::basic) {
710      model_->setStatus(iPivot,ClpSimplex::isFixed);
711      status = ClpSimplex::basic; // so will skip
712    }
713  }
714  switch(status) {
715     
716  case ClpSimplex::basic:
717  case ClpSimplex::superBasic:
718  case ClpSimplex::isFree:
719    break;
720  case ClpSimplex::atUpperBound:
721  case ClpSimplex::atLowerBound:
722  case ClpSimplex::isFixed:
723    // set correctly
724    if (fabs(value-lower)<=primalTolerance*1.001){
725      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
726    } else if (fabs(value-upper)<=primalTolerance*1.001){
727      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
728    } else {
729      // set superBasic
730      model_->setStatus(iPivot,ClpSimplex::superBasic);
731    }
732    break;
733  }
734  double difference = cost-cost_[iRange]; 
735  cost = cost_[iRange];
736  changeCost_ += value*difference;
737  return difference;
738}
739// Returns nearest bound
740double 
741ClpNonLinearCost::nearest(int sequence, double solutionValue)
742{
743  assert (model_!=NULL);
744  // get where in bound sequence
745  int iRange;
746  int start = start_[sequence];
747  int end = start_[sequence+1];
748  int jRange=-1;
749  double nearest=COIN_DBL_MAX;
750  for (iRange=start; iRange<end;iRange++) {
751    if (fabs(solutionValue-lower_[iRange])<nearest) {
752      jRange=iRange;
753      nearest=fabs(solutionValue-lower_[iRange]);
754    }
755  }
756  assert(jRange<end);
757  return lower_[jRange];
758}
759// Sets inside bounds (i.e. non infinite - used in QP
760void 
761ClpNonLinearCost::setBounds(int sequence, double lower, double upper)
762{
763  int start = start_[sequence];
764  assert(start_[sequence+1]==start+4);
765  lower_[start+1]=lower;
766  lower_[start+2]=upper;
767}
768
Note: See TracBrowser for help on using the repository browser.