source: branches/devel-1/ClpNonLinearCost.cpp @ 10

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

Improving reliability using IBM Burlington problems (Primal)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.0 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include <iostream>
9
10#include "ClpNonLinearCost.hpp"
11#include "ClpSimplex.hpp"
12#include "OsiIndexedVector.hpp"
13
14//#############################################################################
15// Constructors / Destructor / Assignment
16//#############################################################################
17
18//-------------------------------------------------------------------
19// Default Constructor
20//-------------------------------------------------------------------
21ClpNonLinearCost::ClpNonLinearCost () :
22  numberRows_(0),
23  numberColumns_(0),
24  start_(NULL),
25  whichRange_(NULL),
26  offset_(NULL),
27  lower_(NULL),
28  cost_(NULL),
29  model_(NULL),
30  numberInfeasibilities_(-1),
31  changeCost_(0.0),
32  largestInfeasibility_(0.0),
33  sumInfeasibilities_(0.0),
34  convex_(true)
35{
36
37}
38/* Constructor from simplex.
39   This will just set up wasteful arrays for linear, but
40   later may do dual analysis and even finding duplicate columns
41*/
42ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model)
43{
44  model_ = model;
45  numberRows_ = model_->numberRows();
46  numberColumns_ = model_->numberColumns();
47  int numberTotal = numberRows_+numberColumns_;
48  convex_ = true;
49  start_ = new int [numberTotal+1];
50  whichRange_ = new int [numberTotal];
51  offset_ = new int [numberTotal];
52  memset(offset_,0,numberTotal*sizeof(int));
53
54  numberInfeasibilities_=0;
55  changeCost_=0.0;
56  double infeasibilityCost = model_->infeasibilityCost();
57  sumInfeasibilities_=0.0;
58  largestInfeasibility_=0.0;
59
60  // First see how much space we need
61  int put=0;
62
63  int iSequence;
64  double * upper = model_->upperRegion();
65  double * lower = model_->lowerRegion();
66  double * cost = model_->costRegion();
67
68  for (iSequence=0;iSequence<numberTotal;iSequence++) {
69    if (upper[iSequence]<1.0e20)
70      put++;
71    put += 3;
72  }
73
74  lower_ = new double [put];
75  cost_ = new double [put];
76
77  put=0;
78
79  start_[0]=0;
80
81  for (iSequence=0;iSequence<numberTotal;iSequence++) {
82
83    lower_[put] = -DBL_MAX;
84    cost_[put++] = cost[iSequence]-infeasibilityCost;
85    whichRange_[iSequence]=put;
86    lower_[put] = lower[iSequence];
87    cost_[put++] = cost[iSequence];
88    lower_[put] = upper[iSequence];
89    cost_[put++] = cost[iSequence]+infeasibilityCost;
90    if (upper[iSequence]<1.0e20) {
91      lower_[put] = DBL_MAX;
92      cost_[put++] = 1.0e50;
93    }
94    start_[iSequence+1]=put;
95  }
96
97}
98
99ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model,const int * starts,
100                   const double * lowerNon, const double * costNon)
101{
102
103  // what about scaling? - only try without it initially
104  assert(!model->scalingFlag());
105  model_ = model;
106  numberRows_ = model_->numberRows();
107  numberColumns_ = model_->numberColumns();
108  int numberTotal = numberRows_+numberColumns_;
109  convex_ = true;
110  start_ = new int [numberTotal+1];
111  whichRange_ = new int [numberTotal];
112  offset_ = new int [numberTotal];
113  memset(offset_,0,numberTotal*sizeof(int));
114
115  numberInfeasibilities_=0;
116  changeCost_=0.0;
117  double infeasibilityCost = model_->infeasibilityCost();
118  largestInfeasibility_=0.0;
119  sumInfeasibilities_=0.0;
120
121  int iSequence;
122  double * upper = model_->upperRegion();
123  double * lower = model_->lowerRegion();
124  double * cost = model_->costRegion();
125  // First see how much space we need - we know column part
126  // but not infeasibilities - see how much extra
127  // may be over-estimate
128  int put=starts[numberColumns_]-numberColumns_;
129
130  for (iSequence=0;iSequence<numberTotal;iSequence++) {
131    if (lower[iSequence]>-1.0e20)
132      if (upper[iSequence]<1.0e20)
133        put++;
134    put += 3;
135  }
136
137  lower_ = new double [put];
138  cost_ = new double [put];
139
140  // now fill in
141  put=0;
142
143  start_[0]=0;
144
145  for (iSequence=0;iSequence<numberTotal;iSequence++) {
146    lower_[put] = -DBL_MAX;
147    cost_[put++] = cost[iSequence]-infeasibilityCost;
148    whichRange_[iSequence]=put;
149    double thisCost;
150    if (iSequence>=numberColumns_) {
151      // rows
152      lower_[put] = lower[iSequence];
153      cost_[put++] = cost[iSequence];
154      thisCost = cost[iSequence];
155    } else {
156      // columns - move costs and see if convex
157      int iIndex = starts[iSequence];
158      int end = starts[iSequence+1];
159      assert (fabs(lower[iSequence]-lowerNon[iIndex])<1.0e-8);
160      thisCost = -DBL_MAX;
161      for (;iIndex<end;iIndex++) {
162        if (lowerNon[iIndex]<upper[iSequence]-1.0e-8) {
163          lower_[put] = lowerNon[iIndex];
164          cost_[put++] = costNon[iIndex];
165          // check convexity
166          if (costNon[iIndex]<thisCost-1.0e-12)
167            convex_ = false;
168          thisCost = costNon[iIndex];
169        } else {
170          break;
171        }
172      }
173    }
174    lower_[put] = upper[iSequence];
175    cost_[put++] = thisCost+infeasibilityCost;
176    if (upper[iSequence]<1.0e20) {
177      lower_[put] = DBL_MAX;
178      cost_[put++] = 1.0e50;
179    }
180    start_[iSequence+1]=put;
181  }
182  // can't handle non-convex at present
183  assert(convex_);
184}
185
186//-------------------------------------------------------------------
187// Copy constructor
188//-------------------------------------------------------------------
189ClpNonLinearCost::ClpNonLinearCost (const ClpNonLinearCost & rhs) :
190  numberRows_(rhs.numberRows_),
191  numberColumns_(rhs.numberColumns_),
192  start_(NULL),
193  whichRange_(NULL),
194  offset_(NULL),
195  lower_(NULL),
196  cost_(NULL),
197  model_(NULL),
198  numberInfeasibilities_(-1),
199  changeCost_(0.0),
200  largestInfeasibility_(0.0),
201  sumInfeasibilities_(0.0),
202  convex_(true)
203{ 
204  if (numberRows_) {
205    int numberTotal = numberRows_+numberColumns_;
206    start_ = new int [numberTotal+1];
207    memcpy(start_,rhs.start_,(numberTotal+1)*sizeof(int));
208    whichRange_ = new int [numberTotal];
209    memcpy(whichRange_,rhs.whichRange_,numberTotal*sizeof(int));
210    offset_ = new int [numberTotal];
211    memcpy(offset_,rhs.offset_,numberTotal*sizeof(int));
212    int numberEntries = start_[numberTotal];
213    lower_ = new double [numberEntries];
214    memcpy(lower_,rhs.lower_,numberEntries*sizeof(double));
215    cost_ = new double [numberEntries];
216    memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
217    model_ = rhs.model_;
218    numberInfeasibilities_=rhs.numberInfeasibilities_;
219    changeCost_ = rhs.changeCost_;
220    largestInfeasibility_ = rhs.largestInfeasibility_;
221    sumInfeasibilities_ = rhs.sumInfeasibilities_;
222    convex_ = rhs.convex_;
223  }
224}
225
226//-------------------------------------------------------------------
227// Destructor
228//-------------------------------------------------------------------
229ClpNonLinearCost::~ClpNonLinearCost ()
230{
231  delete [] start_;
232  delete [] whichRange_;
233  delete [] offset_;
234  delete [] lower_;
235  delete [] cost_;
236
237}
238
239//----------------------------------------------------------------
240// Assignment operator
241//-------------------------------------------------------------------
242ClpNonLinearCost &
243ClpNonLinearCost::operator=(const ClpNonLinearCost& rhs)
244{
245  if (this != &rhs) {
246    numberRows_ = rhs.numberRows_;
247    numberColumns_ = rhs.numberColumns_;
248    delete [] start_;
249    delete [] whichRange_;
250    delete [] offset_;
251    delete [] lower_;
252    delete []cost_;
253    start_ = NULL;
254    whichRange_ = NULL;
255    lower_ = NULL;
256    cost_= NULL;
257    if (numberRows_) {
258      int numberTotal = numberRows_+numberColumns_;
259      start_ = new int [numberTotal+1];
260      memcpy(start_,rhs.start_,(numberTotal+1)*sizeof(int));
261      whichRange_ = new int [numberTotal];
262      memcpy(whichRange_,rhs.whichRange_,numberTotal*sizeof(int));
263      offset_ = new int [numberTotal];
264      memcpy(offset_,rhs.offset_,numberTotal*sizeof(int));
265      int numberEntries = start_[numberTotal];
266      lower_ = new double [numberEntries];
267      memcpy(lower_,rhs.lower_,numberEntries*sizeof(double));
268      cost_ = new double [numberEntries];
269      memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
270    }
271    model_ = rhs.model_;
272    numberInfeasibilities_=rhs.numberInfeasibilities_;
273    changeCost_ = rhs.changeCost_;
274    largestInfeasibility_ = rhs.largestInfeasibility_;
275    sumInfeasibilities_ = rhs.sumInfeasibilities_;
276    convex_ = rhs.convex_;
277  }
278  return *this;
279}
280
281// Changes infeasible costs and computes number and cost of infeas
282// We will need to re-think objective offsets later
283// We will also need a 2 bit per variable array for some
284// purpose which will come to me later
285void 
286ClpNonLinearCost::checkInfeasibilities(bool toNearest)
287{
288  numberInfeasibilities_=0;
289  double infeasibilityCost = model_->infeasibilityCost();
290  changeCost_=0.0;
291  largestInfeasibility_ = 0.0;
292  sumInfeasibilities_ = 0.0;
293  double primalTolerance = model_->currentPrimalTolerance();
294 
295  int iSequence;
296  double * solution = model_->solutionRegion();
297  double * upper = model_->upperRegion();
298  double * lower = model_->lowerRegion();
299  double * cost = model_->costRegion();
300   
301  // nonbasic should be at a valid bound
302  for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
303    double lowerValue=lower[iSequence];
304    double upperValue=upper[iSequence];
305    double value=solution[iSequence];
306    int iRange;
307    // correct costs
308   
309    if (lowerValue>-1.0e20) {
310      iRange=start_[iSequence];
311      cost_[iRange] = cost[iSequence]-infeasibilityCost;
312    }
313    if (upperValue<1.0e20) {
314      iRange=start_[iSequence+1]-2;
315      cost_[iRange] = cost[iSequence]+infeasibilityCost;
316    }
317    // get correct place
318    int start = start_[iSequence];
319    int end = start_[iSequence+1]-1;
320    for (iRange=start; iRange<end;iRange++) {
321      if (value<lower_[iRange+1]+primalTolerance) {
322        // put in better range if infeasible
323        if (value>=lower_[iRange+1]-primalTolerance&&iRange==start) 
324          iRange++;
325        whichRange_[iSequence]=iRange;
326        break;
327      }
328    }
329    assert(iRange<end);
330    switch(model_->getStatus(iSequence)) {
331     
332    case ClpSimplex::basic:
333    case ClpSimplex::superBasic:
334      // iRange is in correct place
335      // slot in here
336      if (value<lower[iSequence]-primalTolerance) {
337        value = lower[iSequence]-value;
338        sumInfeasibilities_ += value;
339        largestInfeasibility_ = max(largestInfeasibility_,value);
340        changeCost_ -= lower[iSequence]*
341          (cost_[iRange]-cost[iSequence]);
342        numberInfeasibilities_++;
343      } else if (value>upper[iSequence]+primalTolerance) {
344        value = value-upper[iSequence];
345        sumInfeasibilities_ += value;
346        largestInfeasibility_ = max(largestInfeasibility_,value);
347        changeCost_ -= upper[iSequence]*
348          (cost_[iRange]-cost[iSequence]);
349        numberInfeasibilities_++;
350      }
351      lower[iSequence] = lower_[iRange];
352      upper[iSequence] = lower_[iRange+1];
353      cost[iSequence] = cost_[iRange];
354      break;
355    case ClpSimplex::isFree:
356      if (toNearest)
357        solution[iSequence] = 0.0;
358      break;
359    case ClpSimplex::atUpperBound:
360      if (!toNearest) {
361        assert(fabs(value-upperValue)<=primalTolerance*1.0001) ;
362      } else {
363        solution[iSequence] = upperValue;
364      }
365      break;
366    case ClpSimplex::atLowerBound:
367      if (!toNearest) {
368        assert(fabs(value-lowerValue)<=primalTolerance*1.0001); 
369      } else {
370        solution[iSequence] = lowerValue;
371      }
372      break;
373    }
374  }
375}
376/* Goes through one bound for each variable.
377   If array[i]*multiplier>0 goes down, otherwise up.
378   The indices are row indices and need converting to sequences
379*/
380void 
381ClpNonLinearCost::goThru(int numberInArray, double multiplier,
382              const int * index, const double * array,
383                         double * rhs)
384{
385  assert (model_!=NULL);
386  const int * pivotVariable = model_->pivotVariable();
387  int i;
388  for (i=0;i<numberInArray;i++) {
389    int iRow = index[i];
390    int iPivot = pivotVariable[iRow];
391    double alpha = multiplier*array[iRow];
392    // get where in bound sequence
393    int iRange = whichRange_[iPivot];
394    iRange += offset_[iPivot]; //add temporary bias
395    double value = model_->solution(iPivot);
396    if (alpha>0.0) {
397      // down one
398      iRange--;
399      assert(iRange>=start_[iPivot]);
400      rhs[iRow] = value - lower_[iRange];
401    } else {
402      // up one
403      iRange++;
404      assert(iRange<start_[iPivot+1]-1);
405      rhs[iRow] = lower_[iRange+1] - value;
406    }
407    offset_[iPivot] = iRange - whichRange_[iPivot];
408  }
409}
410/* Takes off last iteration (i.e. offsets closer to 0)
411 */
412void 
413ClpNonLinearCost::goBack(int numberInArray, const int * index, 
414              double * rhs)
415{
416  assert (model_!=NULL);
417  const int * pivotVariable = model_->pivotVariable();
418  int i;
419  for (i=0;i<numberInArray;i++) {
420    int iRow = index[i];
421    int iPivot = pivotVariable[iRow];
422    // get where in bound sequence
423    int iRange = whichRange_[iPivot];
424    bool down;
425    // get closer to original
426    if (offset_[iPivot]>0) {
427      offset_[iPivot]--;
428      assert (offset_[iPivot]>=0);
429      down = false;
430    } else {
431      offset_[iPivot]++;
432      assert (offset_[iPivot]<=0);
433      down = true;
434    }
435    iRange += offset_[iPivot]; //add temporary bias
436    double value = model_->solution(iPivot);
437    if (down) {
438      // down one
439      assert(iRange>=start_[iPivot]);
440      rhs[iRow] = value - lower_[iRange];
441    } else {
442      // up one
443      assert(iRange<start_[iPivot+1]-1);
444      rhs[iRow] = lower_[iRange+1] - value;
445    }
446  }
447}
448void 
449ClpNonLinearCost::goBackAll(const OsiIndexedVector * update)
450{
451  assert (model_!=NULL);
452  const int * pivotVariable = model_->pivotVariable();
453  int i;
454  int number = update->getNumElements();
455  const int * index = update->getIndices();
456  for (i=0;i<number;i++) {
457    int iRow = index[i];
458    int iPivot = pivotVariable[iRow];
459    offset_[iPivot]=0;
460  }
461#ifdef CLP_DEBUG
462  for (i=0;i<numberRows_+numberColumns_;i++) 
463    assert(!offset_[i]);
464#endif
465}
466void 
467ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
468{
469  assert (model_!=NULL);
470  double primalTolerance = model_->currentPrimalTolerance();
471  const int * pivotVariable = model_->pivotVariable();
472  int i;
473  for (i=0;i<numberInArray;i++) {
474    int iRow = index[i];
475    int iPivot = pivotVariable[iRow];
476    // get where in bound sequence
477    int iRange;
478    double value = model_->solution(iPivot);
479    int start = start_[iPivot];
480    int end = start_[iPivot+1]-1;
481    for (iRange=start; iRange<end;iRange++) {
482      if (value<lower_[iRange+1]+primalTolerance) {
483        // put in better range
484        if (value>=lower_[iRange+1]-primalTolerance&&iRange==start) 
485          iRange++;
486        break;
487      }
488    }
489    assert(iRange<end);
490    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
491    double & lower = model_->lowerAddress(iPivot);
492    double & upper = model_->upperAddress(iPivot);
493    double & cost = model_->costAddress(iPivot);
494    whichRange_[iPivot]=iRange;
495    lower = lower_[iRange];
496    upper = lower_[iRange+1];
497    cost = cost_[iRange];
498  }
499}
500/* Puts back correct infeasible costs for each variable
501   The input indices are row indices and need converting to sequences
502   for costs.
503   On input array is empty (but indices exist).  On exit just
504   changed costs will be stored as normal OsiIndexedVector
505*/
506void 
507ClpNonLinearCost::checkChanged(int numberInArray, OsiIndexedVector * update)
508{
509  assert (model_!=NULL);
510  double primalTolerance = model_->currentPrimalTolerance();
511  const int * pivotVariable = model_->pivotVariable();
512  int number=0;
513  int * index = update->getIndices();
514  double * work = update->denseVector();
515  int i;
516  for (i=0;i<numberInArray;i++) {
517    int iRow = index[i];
518    int iPivot = pivotVariable[iRow];
519    // get where in bound sequence
520    int iRange;
521    double value = model_->solution(iPivot);
522    int start = start_[iPivot];
523    int end = start_[iPivot+1]-1;
524    for (iRange=start; iRange<end;iRange++) {
525      if (value<lower_[iRange+1]+primalTolerance) {
526        // put in better range
527        if (value>=lower_[iRange+1]-primalTolerance&&iRange==start) 
528          iRange++;
529        break;
530      }
531    }
532    assert(iRange<end);
533    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
534    int jRange = whichRange_[iPivot];
535    if (iRange!=jRange) {
536      // changed
537      work[iRow] = cost_[jRange]-cost_[iRange];
538      index[number++]=iRow;
539      double & lower = model_->lowerAddress(iPivot);
540      double & upper = model_->upperAddress(iPivot);
541      double & cost = model_->costAddress(iPivot);
542      whichRange_[iPivot]=iRange;
543      lower = lower_[iRange];
544      upper = lower_[iRange+1];
545      cost = cost_[iRange];
546    }
547  }
548  update->setNumElements(number);
549}
550/* Sets bounds and cost for one variable - returns change in cost*/
551double 
552ClpNonLinearCost::setOne(int iPivot, double value)
553{
554  assert (model_!=NULL);
555  double primalTolerance = model_->currentPrimalTolerance();
556  // get where in bound sequence
557  int iRange;
558  int start = start_[iPivot];
559  int end = start_[iPivot+1]-1;
560  for (iRange=start; iRange<end;iRange++) {
561    if (value<lower_[iRange+1]+primalTolerance) {
562      // put in better range
563      if (value>=lower_[iRange+1]-primalTolerance&&iRange==start) 
564        iRange++;
565      break;
566    }
567  }
568  assert(iRange<end);
569  whichRange_[iPivot]=iRange;
570  double & lower = model_->lowerAddress(iPivot);
571  double & upper = model_->upperAddress(iPivot);
572  double & cost = model_->costAddress(iPivot);
573  lower = lower_[iRange];
574  upper = lower_[iRange+1];
575  switch(model_->getStatus(iPivot)) {
576     
577  case ClpSimplex::basic:
578  case ClpSimplex::superBasic:
579  case ClpSimplex::isFree:
580    break;
581  case ClpSimplex::atUpperBound:
582  case ClpSimplex::atLowerBound:
583    // set correctly
584    if (fabs(value-lower)<=primalTolerance*1.001) 
585      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
586    else if (fabs(value-upper)<=primalTolerance*1.001) 
587      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
588    else
589      assert(fabs(value-lower)<=primalTolerance*1.0001||
590             fabs(value-upper)<=primalTolerance*1.0001);
591    break;
592  }
593  if (upper-lower<1.0e-8)
594    model_->setFixed(iPivot);
595  else
596    model_->clearFixed(iPivot);
597  double difference = cost-cost_[iRange]; 
598  cost = cost_[iRange];
599  changeCost_ += value*difference;
600  return difference;
601}
602// Returns nearest bound
603double 
604ClpNonLinearCost::nearest(int sequence, double solutionValue)
605{
606  assert (model_!=NULL);
607  // get where in bound sequence
608  int iRange;
609  int start = start_[sequence];
610  int end = start_[sequence+1];
611  int jRange=-1;
612  double nearest=DBL_MAX;
613  for (iRange=start; iRange<end;iRange++) {
614    if (fabs(solutionValue-lower_[iRange])<nearest) {
615      jRange=iRange;
616      nearest=fabs(solutionValue-lower_[iRange]);
617    }
618  }
619  assert(jRange<end);
620  return lower_[jRange];
621}
622
Note: See TracBrowser for help on using the repository browser.