source: branches/pre/ClpNonLinearCost.cpp @ 213

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

Trying

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