source: trunk/ClpNonLinearCost.cpp @ 2

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

Adding Clp to development branch

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