source: branches/pre/ClpNonLinearCost.cpp @ 190

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

Stuff

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