source: trunk/ClpNonLinearCost.cpp @ 284

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

Cosmetic printing when free variables

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