source: trunk/ClpNonLinearCost.cpp @ 225

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

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.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{
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          if (fabs(value-lowerValue)<=oldTolerance*1.0001) {
435            if  (fabs(value-lowerValue)>primalTolerance)
436              solution[iSequence]=lowerValue;
437            model_->setStatus(iSequence,ClpSimplex::atLowerBound);
438          } else {
439            model_->setStatus(iSequence,ClpSimplex::superBasic);
440          }
441        } else if  (fabs(value-upperValue)>primalTolerance) {
442          solution[iSequence]=upperValue;
443        }
444      } else {
445        // Set to nearest and make at upper bound
446        int kRange;
447        iRange=-1;
448        double nearest = COIN_DBL_MAX;
449        for (kRange=start; kRange<end;kRange++) {
450          if (fabs(lower_[kRange]-value)<nearest) {
451            nearest = fabs(lower_[kRange]-value);
452            iRange=kRange;
453          }
454        }
455        assert (iRange>=0);
456        iRange--;
457        solution[iSequence]=lower_[iRange+1];
458      }
459      break;
460    case ClpSimplex::atLowerBound:
461      if (!toNearest) {
462        // With increasing tolerances - we may be at wrong place
463        if (fabs(value-lowerValue)>oldTolerance*1.0001) {
464          if (fabs(value-upperValue)<=oldTolerance*1.0001) {
465            if  (fabs(value-upperValue)>primalTolerance)
466              solution[iSequence]=upperValue;
467            model_->setStatus(iSequence,ClpSimplex::atLowerBound);
468          } else {
469            model_->setStatus(iSequence,ClpSimplex::superBasic);
470          }
471        } else if  (fabs(value-lowerValue)>primalTolerance) {
472          solution[iSequence]=lowerValue;
473        }
474      } else {
475        // Set to nearest and make at lower bound
476        int kRange;
477        iRange=-1;
478        double nearest = COIN_DBL_MAX;
479        for (kRange=start; kRange<end;kRange++) {
480          if (fabs(lower_[kRange]-value)<nearest) {
481            nearest = fabs(lower_[kRange]-value);
482            iRange=kRange;
483          }
484        }
485        assert (iRange>=0);
486        solution[iSequence]=lower_[iRange];
487      }
488      break;
489    case ClpSimplex::isFixed:
490      if (toNearest) {
491        // Set to true fixed
492        for (iRange=start; iRange<end;iRange++) {
493          if (lower_[iRange]==lower_[iRange+1])
494            break;
495        }
496        solution[iSequence]=lower_[iRange];
497      }
498      break;
499    }
500    lower[iSequence] = lower_[iRange];
501    upper[iSequence] = lower_[iRange+1];
502    cost[iSequence] = cost_[iRange];
503  }
504}
505/* Goes through one bound for each variable.
506   If array[i]*multiplier>0 goes down, otherwise up.
507   The indices are row indices and need converting to sequences
508*/
509void 
510ClpNonLinearCost::goThru(int numberInArray, double multiplier,
511              const int * index, const double * array,
512                         double * rhs)
513{
514  assert (model_!=NULL);
515  const int * pivotVariable = model_->pivotVariable();
516  int i;
517  for (i=0;i<numberInArray;i++) {
518    int iRow = index[i];
519    int iPivot = pivotVariable[iRow];
520    double alpha = multiplier*array[iRow];
521    // get where in bound sequence
522    int iRange = whichRange_[iPivot];
523    iRange += offset_[iPivot]; //add temporary bias
524    double value = model_->solution(iPivot);
525    if (alpha>0.0) {
526      // down one
527      iRange--;
528      assert(iRange>=start_[iPivot]);
529      rhs[iRow] = value - lower_[iRange];
530    } else {
531      // up one
532      iRange++;
533      assert(iRange<start_[iPivot+1]-1);
534      rhs[iRow] = lower_[iRange+1] - value;
535    }
536    offset_[iPivot] = iRange - whichRange_[iPivot];
537  }
538}
539/* Takes off last iteration (i.e. offsets closer to 0)
540 */
541void 
542ClpNonLinearCost::goBack(int numberInArray, const int * index, 
543              double * rhs)
544{
545  assert (model_!=NULL);
546  const int * pivotVariable = model_->pivotVariable();
547  int i;
548  for (i=0;i<numberInArray;i++) {
549    int iRow = index[i];
550    int iPivot = pivotVariable[iRow];
551    // get where in bound sequence
552    int iRange = whichRange_[iPivot];
553    bool down;
554    // get closer to original
555    if (offset_[iPivot]>0) {
556      offset_[iPivot]--;
557      assert (offset_[iPivot]>=0);
558      down = false;
559    } else {
560      offset_[iPivot]++;
561      assert (offset_[iPivot]<=0);
562      down = true;
563    }
564    iRange += offset_[iPivot]; //add temporary bias
565    double value = model_->solution(iPivot);
566    if (down) {
567      // down one
568      assert(iRange>=start_[iPivot]);
569      rhs[iRow] = value - lower_[iRange];
570    } else {
571      // up one
572      assert(iRange<start_[iPivot+1]-1);
573      rhs[iRow] = lower_[iRange+1] - value;
574    }
575  }
576}
577void 
578ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
579{
580  assert (model_!=NULL);
581  const int * pivotVariable = model_->pivotVariable();
582  int i;
583  int number = update->getNumElements();
584  const int * index = update->getIndices();
585  for (i=0;i<number;i++) {
586    int iRow = index[i];
587    int iPivot = pivotVariable[iRow];
588    offset_[iPivot]=0;
589  }
590#ifdef CLP_DEBUG
591  for (i=0;i<numberRows_+numberColumns_;i++) 
592    assert(!offset_[i]);
593#endif
594}
595void 
596ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
597{
598  assert (model_!=NULL);
599  double primalTolerance = model_->currentPrimalTolerance();
600  const int * pivotVariable = model_->pivotVariable();
601  int i;
602  for (i=0;i<numberInArray;i++) {
603    int iRow = index[i];
604    int iPivot = pivotVariable[iRow];
605    // get where in bound sequence
606    int iRange;
607    int currentRange = whichRange_[iPivot];
608    double value = model_->solution(iPivot);
609    int start = start_[iPivot];
610    int end = start_[iPivot+1]-1;
611    for (iRange=start; iRange<end;iRange++) {
612      if (value<lower_[iRange+1]+primalTolerance) {
613        // put in better range
614        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
615          iRange++;
616        break;
617      }
618    }
619    assert(iRange<end);
620    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
621    double & lower = model_->lowerAddress(iPivot);
622    double & upper = model_->upperAddress(iPivot);
623    double & cost = model_->costAddress(iPivot);
624    whichRange_[iPivot]=iRange;
625    if (iRange!=currentRange) {
626      if (infeasible(iRange))
627        numberInfeasibilities_++;
628      if (infeasible(currentRange))
629        numberInfeasibilities_--;
630    }
631    lower = lower_[iRange];
632    upper = lower_[iRange+1];
633    cost = cost_[iRange];
634  }
635}
636/* Puts back correct infeasible costs for each variable
637   The input indices are row indices and need converting to sequences
638   for costs.
639   On input array is empty (but indices exist).  On exit just
640   changed costs will be stored as normal CoinIndexedVector
641*/
642void 
643ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
644{
645  assert (model_!=NULL);
646  double primalTolerance = model_->currentPrimalTolerance();
647  const int * pivotVariable = model_->pivotVariable();
648  int number=0;
649  int * index = update->getIndices();
650  double * work = update->denseVector();
651  int i;
652  for (i=0;i<numberInArray;i++) {
653    int iRow = index[i];
654    int iPivot = pivotVariable[iRow];
655    // get where in bound sequence
656    int iRange;
657    double value = model_->solution(iPivot);
658    int start = start_[iPivot];
659    int end = start_[iPivot+1]-1;
660    for (iRange=start; iRange<end;iRange++) {
661      if (value<lower_[iRange+1]+primalTolerance) {
662        // put in better range
663        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
664          iRange++;
665        break;
666      }
667    }
668    assert(iRange<end);
669    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
670    int jRange = whichRange_[iPivot];
671    if (iRange!=jRange) {
672      // changed
673      work[iRow] = cost_[jRange]-cost_[iRange];
674      index[number++]=iRow;
675      double & lower = model_->lowerAddress(iPivot);
676      double & upper = model_->upperAddress(iPivot);
677      double & cost = model_->costAddress(iPivot);
678      whichRange_[iPivot]=iRange;
679      if (infeasible(iRange))
680        numberInfeasibilities_++;
681      if (infeasible(jRange))
682        numberInfeasibilities_--;
683      lower = lower_[iRange];
684      upper = lower_[iRange+1];
685      cost = cost_[iRange];
686    }
687  }
688  update->setNumElements(number);
689}
690/* Sets bounds and cost for one variable - returns change in cost*/
691double 
692ClpNonLinearCost::setOne(int iPivot, double value)
693{
694  assert (model_!=NULL);
695  double primalTolerance = model_->currentPrimalTolerance();
696  // get where in bound sequence
697  int iRange;
698  int currentRange = whichRange_[iPivot];
699  int start = start_[iPivot];
700  int end = start_[iPivot+1]-1;
701  if (!bothWays_) {
702    // If fixed try and get feasible
703    if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
704      iRange =start+1;
705    } else {
706      for (iRange=start; iRange<end;iRange++) {
707        if (value<=lower_[iRange+1]+primalTolerance) {
708          // put in better range
709          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
710            iRange++;
711          break;
712        }
713      }
714    }
715  } else {
716    // leave in current if possible
717    iRange = whichRange_[iPivot];
718    if (value<lower_[iRange]-primalTolerance||value>lower_[iRange+1]+primalTolerance) {
719      for (iRange=start; iRange<end;iRange++) {
720        if (value<lower_[iRange+1]+primalTolerance) {
721          // put in better range
722          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
723            iRange++;
724          break;
725        }
726      }
727    }
728  }
729  assert(iRange<end);
730  whichRange_[iPivot]=iRange;
731  if (iRange!=currentRange) {
732    if (infeasible(iRange))
733      numberInfeasibilities_++;
734    if (infeasible(currentRange))
735      numberInfeasibilities_--;
736  }
737  double & lower = model_->lowerAddress(iPivot);
738  double & upper = model_->upperAddress(iPivot);
739  double & cost = model_->costAddress(iPivot);
740  lower = lower_[iRange];
741  upper = lower_[iRange+1];
742  ClpSimplex::Status status = model_->getStatus(iPivot);
743  if (upper==lower) {
744    if (status != ClpSimplex::basic) {
745      model_->setStatus(iPivot,ClpSimplex::isFixed);
746      status = ClpSimplex::basic; // so will skip
747    }
748  }
749  switch(status) {
750     
751  case ClpSimplex::basic:
752  case ClpSimplex::superBasic:
753  case ClpSimplex::isFree:
754    break;
755  case ClpSimplex::atUpperBound:
756  case ClpSimplex::atLowerBound:
757  case ClpSimplex::isFixed:
758    // set correctly
759    if (fabs(value-lower)<=primalTolerance*1.001){
760      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
761    } else if (fabs(value-upper)<=primalTolerance*1.001){
762      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
763    } else {
764      // set superBasic
765      model_->setStatus(iPivot,ClpSimplex::superBasic);
766    }
767    break;
768  }
769  double difference = cost-cost_[iRange]; 
770  cost = cost_[iRange];
771  changeCost_ += value*difference;
772  return difference;
773}
774/* Sets bounds and cost for outgoing variable
775   may change value
776   Returns direction */
777int 
778ClpNonLinearCost::setOneOutgoing(int iPivot, double & value)
779{
780  assert (model_!=NULL);
781  double primalTolerance = model_->currentPrimalTolerance();
782  // get where in bound sequence
783  int iRange;
784  int currentRange = whichRange_[iPivot];
785  int start = start_[iPivot];
786  int end = start_[iPivot+1]-1;
787  // Set perceived direction out
788  int direction;
789  if (value<=lower_[currentRange]+1.001*primalTolerance) {
790    direction=1;
791  } else if (value>=lower_[currentRange+1]-1.001*primalTolerance) {
792    direction=-1;
793  } else {
794    // odd
795    direction=0;
796  }
797  // If fixed try and get feasible
798  if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
799    iRange =start+1;
800  } else {
801    // See if exact
802    for (iRange=start; iRange<end;iRange++) {
803      if (value==lower_[iRange+1]) {
804        // put in better range
805        if (infeasible(iRange)&&iRange==start) 
806          iRange++;
807        break;
808      }
809    }
810    if (iRange==end) {
811      // not exact
812      for (iRange=start; iRange<end;iRange++) {
813        if (value<=lower_[iRange+1]+primalTolerance) {
814          // put in better range
815          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
816            iRange++;
817          break;
818        }
819      }
820    }
821  }
822  assert(iRange<end);
823  whichRange_[iPivot]=iRange;
824  if (iRange!=currentRange) {
825    if (infeasible(iRange))
826      numberInfeasibilities_++;
827    if (infeasible(currentRange))
828      numberInfeasibilities_--;
829  }
830  double & lower = model_->lowerAddress(iPivot);
831  double & upper = model_->upperAddress(iPivot);
832  double & cost = model_->costAddress(iPivot);
833  lower = lower_[iRange];
834  upper = lower_[iRange+1];
835  if (upper==lower) {
836    value=upper;
837  } else {
838    // set correctly
839    if (fabs(value-lower)<=primalTolerance*1.001){
840      value = min(value,lower+primalTolerance);
841    } else if (fabs(value-upper)<=primalTolerance*1.001){
842      value = max(value,upper-primalTolerance);
843    } else {
844      printf("*** variable wandered off bound %g %g %g!\n",
845             lower,value,upper);
846      if (value-lower<=upper-value) 
847        value = lower+primalTolerance;
848      else 
849        value = upper-primalTolerance;
850    }
851  }
852  double difference = cost-cost_[iRange]; 
853  cost = cost_[iRange];
854  changeCost_ += value*difference;
855  return direction;
856}
857// Returns nearest bound
858double 
859ClpNonLinearCost::nearest(int sequence, double solutionValue)
860{
861  assert (model_!=NULL);
862  // get where in bound sequence
863  int iRange;
864  int start = start_[sequence];
865  int end = start_[sequence+1];
866  int jRange=-1;
867  double nearest=COIN_DBL_MAX;
868  for (iRange=start; iRange<end;iRange++) {
869    if (fabs(solutionValue-lower_[iRange])<nearest) {
870      jRange=iRange;
871      nearest=fabs(solutionValue-lower_[iRange]);
872    }
873  }
874  assert(jRange<end);
875  return lower_[jRange];
876}
877
Note: See TracBrowser for help on using the repository browser.