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

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

Still trying to make more rugged on IBM Burlington problems

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