source: trunk/ClpNonLinearCost.cpp @ 448

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

compiler bug???? and assert in unitTest while Lou fixes presolve

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