source: trunk/ClpNonLinearCost.cpp @ 461

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

Trying to make faster

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