source: trunk/ClpNonLinearCost.cpp @ 114

Last change on this file since 114 was 114, checked in by forrest, 17 years ago

Re-order variables

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