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

Last change on this file since 27 was 19, checked in by ladanyi, 17 years ago

reordering Clp

  • 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 "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        assert(fabs(value-upperValue)<=primalTolerance*1.0001) ;
363      } else {
364        solution[iSequence] = upperValue;
365      }
366      break;
367    case ClpSimplex::atLowerBound:
368      if (!toNearest) {
369        assert(fabs(value-lowerValue)<=primalTolerance*1.0001); 
370      } else {
371        solution[iSequence] = lowerValue;
372      }
373      break;
374    }
375  }
376}
377/* Goes through one bound for each variable.
378   If array[i]*multiplier>0 goes down, otherwise up.
379   The indices are row indices and need converting to sequences
380*/
381void 
382ClpNonLinearCost::goThru(int numberInArray, double multiplier,
383              const int * index, const double * array,
384                         double * rhs)
385{
386  assert (model_!=NULL);
387  const int * pivotVariable = model_->pivotVariable();
388  int i;
389  for (i=0;i<numberInArray;i++) {
390    int iRow = index[i];
391    int iPivot = pivotVariable[iRow];
392    double alpha = multiplier*array[iRow];
393    // get where in bound sequence
394    int iRange = whichRange_[iPivot];
395    iRange += offset_[iPivot]; //add temporary bias
396    double value = model_->solution(iPivot);
397    if (alpha>0.0) {
398      // down one
399      iRange--;
400      assert(iRange>=start_[iPivot]);
401      rhs[iRow] = value - lower_[iRange];
402    } else {
403      // up one
404      iRange++;
405      assert(iRange<start_[iPivot+1]-1);
406      rhs[iRow] = lower_[iRange+1] - value;
407    }
408    offset_[iPivot] = iRange - whichRange_[iPivot];
409  }
410}
411/* Takes off last iteration (i.e. offsets closer to 0)
412 */
413void 
414ClpNonLinearCost::goBack(int numberInArray, const int * index, 
415              double * rhs)
416{
417  assert (model_!=NULL);
418  const int * pivotVariable = model_->pivotVariable();
419  int i;
420  for (i=0;i<numberInArray;i++) {
421    int iRow = index[i];
422    int iPivot = pivotVariable[iRow];
423    // get where in bound sequence
424    int iRange = whichRange_[iPivot];
425    bool down;
426    // get closer to original
427    if (offset_[iPivot]>0) {
428      offset_[iPivot]--;
429      assert (offset_[iPivot]>=0);
430      down = false;
431    } else {
432      offset_[iPivot]++;
433      assert (offset_[iPivot]<=0);
434      down = true;
435    }
436    iRange += offset_[iPivot]; //add temporary bias
437    double value = model_->solution(iPivot);
438    if (down) {
439      // down one
440      assert(iRange>=start_[iPivot]);
441      rhs[iRow] = value - lower_[iRange];
442    } else {
443      // up one
444      assert(iRange<start_[iPivot+1]-1);
445      rhs[iRow] = lower_[iRange+1] - value;
446    }
447  }
448}
449void 
450ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
451{
452  assert (model_!=NULL);
453  const int * pivotVariable = model_->pivotVariable();
454  int i;
455  int number = update->getNumElements();
456  const int * index = update->getIndices();
457  for (i=0;i<number;i++) {
458    int iRow = index[i];
459    int iPivot = pivotVariable[iRow];
460    offset_[iPivot]=0;
461  }
462#ifdef CLP_DEBUG
463  for (i=0;i<numberRows_+numberColumns_;i++) 
464    assert(!offset_[i]);
465#endif
466}
467void 
468ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
469{
470  assert (model_!=NULL);
471  double primalTolerance = model_->currentPrimalTolerance();
472  const int * pivotVariable = model_->pivotVariable();
473  int i;
474  for (i=0;i<numberInArray;i++) {
475    int iRow = index[i];
476    int iPivot = pivotVariable[iRow];
477    // get where in bound sequence
478    int iRange;
479    double value = model_->solution(iPivot);
480    int start = start_[iPivot];
481    int end = start_[iPivot+1]-1;
482    for (iRange=start; iRange<end;iRange++) {
483      if (value<lower_[iRange+1]+primalTolerance) {
484        // put in better range
485        if (value>=lower_[iRange+1]-primalTolerance&&iRange==start) 
486          iRange++;
487        break;
488      }
489    }
490    assert(iRange<end);
491    assert(model_->getStatus(iPivot)==ClpSimplex::basic);
492    double & lower = model_->lowerAddress(iPivot);
493    double & upper = model_->upperAddress(iPivot);
494    double & cost = model_->costAddress(iPivot);
495    whichRange_[iPivot]=iRange;
496    lower = lower_[iRange];
497    upper = lower_[iRange+1];
498    cost = cost_[iRange];
499  }
500}
501/* Puts back correct infeasible costs for each variable
502   The input indices are row indices and need converting to sequences
503   for costs.
504   On input array is empty (but indices exist).  On exit just
505   changed costs will be stored as normal CoinIndexedVector
506*/
507void 
508ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
509{
510  assert (model_!=NULL);
511  double primalTolerance = model_->currentPrimalTolerance();
512  const int * pivotVariable = model_->pivotVariable();
513  int number=0;
514  int * index = update->getIndices();
515  double * work = update->denseVector();
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    int jRange = whichRange_[iPivot];
536    if (iRange!=jRange) {
537      // changed
538      work[iRow] = cost_[jRange]-cost_[iRange];
539      index[number++]=iRow;
540      double & lower = model_->lowerAddress(iPivot);
541      double & upper = model_->upperAddress(iPivot);
542      double & cost = model_->costAddress(iPivot);
543      whichRange_[iPivot]=iRange;
544      lower = lower_[iRange];
545      upper = lower_[iRange+1];
546      cost = cost_[iRange];
547    }
548  }
549  update->setNumElements(number);
550}
551/* Sets bounds and cost for one variable - returns change in cost*/
552double 
553ClpNonLinearCost::setOne(int iPivot, double value)
554{
555  assert (model_!=NULL);
556  double primalTolerance = model_->currentPrimalTolerance();
557  // get where in bound sequence
558  int iRange;
559  int start = start_[iPivot];
560  int end = start_[iPivot+1]-1;
561  for (iRange=start; iRange<end;iRange++) {
562    if (value<lower_[iRange+1]+primalTolerance) {
563      // put in better range
564      if (value>=lower_[iRange+1]-primalTolerance&&iRange==start) 
565        iRange++;
566      break;
567    }
568  }
569  assert(iRange<end);
570  whichRange_[iPivot]=iRange;
571  double & lower = model_->lowerAddress(iPivot);
572  double & upper = model_->upperAddress(iPivot);
573  double & cost = model_->costAddress(iPivot);
574  lower = lower_[iRange];
575  upper = lower_[iRange+1];
576  switch(model_->getStatus(iPivot)) {
577     
578  case ClpSimplex::basic:
579  case ClpSimplex::superBasic:
580  case ClpSimplex::isFree:
581    break;
582  case ClpSimplex::atUpperBound:
583  case ClpSimplex::atLowerBound:
584    // set correctly
585    if (fabs(value-lower)<=primalTolerance*1.001) 
586      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
587    else if (fabs(value-upper)<=primalTolerance*1.001) 
588      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
589    else
590      assert(fabs(value-lower)<=primalTolerance*1.0001||
591             fabs(value-upper)<=primalTolerance*1.0001);
592    break;
593  }
594  if (upper-lower<1.0e-8)
595    model_->setFixed(iPivot);
596  else
597    model_->clearFixed(iPivot);
598  double difference = cost-cost_[iRange]; 
599  cost = cost_[iRange];
600  changeCost_ += value*difference;
601  return difference;
602}
603// Returns nearest bound
604double 
605ClpNonLinearCost::nearest(int sequence, double solutionValue)
606{
607  assert (model_!=NULL);
608  // get where in bound sequence
609  int iRange;
610  int start = start_[sequence];
611  int end = start_[sequence+1];
612  int jRange=-1;
613  double nearest=DBL_MAX;
614  for (iRange=start; iRange<end;iRange++) {
615    if (fabs(solutionValue-lower_[iRange])<nearest) {
616      jRange=iRange;
617      nearest=fabs(solutionValue-lower_[iRange]);
618    }
619  }
620  assert(jRange<end);
621  return lower_[jRange];
622}
623
Note: See TracBrowser for help on using the repository browser.