source: trunk/Clp/src/AbcNonLinearCost.cpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 10 months ago

formatting

  • Property svn:keywords set to Id
File size: 32.8 KB
Line 
1/* $Id: AbcNonLinearCost.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include <iostream>
8#include <cassert>
9
10#include "CoinIndexedVector.hpp"
11
12#include "AbcSimplex.hpp"
13#include "CoinHelperFunctions.hpp"
14#include "AbcNonLinearCost.hpp"
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22AbcNonLinearCost::AbcNonLinearCost()
23  : changeCost_(0.0)
24  , feasibleCost_(0.0)
25  , infeasibilityWeight_(-1.0)
26  , largestInfeasibility_(0.0)
27  , sumInfeasibilities_(0.0)
28  , averageTheta_(0.0)
29  , numberRows_(0)
30  , numberColumns_(0)
31  , model_(NULL)
32  , numberInfeasibilities_(-1)
33  , status_(NULL)
34  , bound_(NULL)
35  , cost_(NULL)
36{
37}
38//#define VALIDATE
39#ifdef VALIDATE
40static double *saveLowerV = NULL;
41static double *saveUpperV = NULL;
42#ifdef NDEBUG
43Validate should not be set if no debug
44#endif
45#endif
46  /* Constructor from simplex.
47                                    This will just set up wasteful arrays for linear, but
48                                    later may do dual analysis and even finding duplicate columns
49                                 */
50  AbcNonLinearCost::AbcNonLinearCost(AbcSimplex * model)
51{
52  model_ = model;
53  numberRows_ = model_->numberRows();
54  numberColumns_ = model_->numberColumns();
55  // If gub then we need this extra
56  int numberTotal = numberRows_ + numberColumns_;
57  numberInfeasibilities_ = 0;
58  changeCost_ = 0.0;
59  feasibleCost_ = 0.0;
60  infeasibilityWeight_ = -1.0;
61  double *cost = model_->costRegion();
62  // check if all 0
63  int iSequence;
64  bool allZero = true;
65  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
66    if (cost[iSequence]) {
67      allZero = false;
68      break;
69    }
70  }
71  if (allZero)
72    model_->setInfeasibilityCost(1.0);
73  sumInfeasibilities_ = 0.0;
74  averageTheta_ = 0.0;
75  largestInfeasibility_ = 0.0;
76  bound_ = new double[numberTotal];
77  cost_ = new double[numberTotal];
78  status_ = new unsigned char[numberTotal];
79  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
80    bound_[iSequence] = 0.0;
81    cost_[iSequence] = cost[iSequence];
82    setInitialStatus(status_[iSequence]);
83  }
84}
85// Refresh - assuming regions OK
86void AbcNonLinearCost::refresh()
87{
88  int numberTotal = numberRows_ + numberColumns_;
89  numberInfeasibilities_ = 0;
90  sumInfeasibilities_ = 0.0;
91  largestInfeasibility_ = 0.0;
92  double infeasibilityCost = model_->infeasibilityCost();
93  double primalTolerance = model_->currentPrimalTolerance();
94  double *cost = model_->costRegion();
95  double *upper = model_->upperRegion();
96  double *lower = model_->lowerRegion();
97  double *solution = model_->solutionRegion();
98  for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
99    cost_[iSequence] = cost[iSequence];
100    double value = solution[iSequence];
101    double lowerValue = lower[iSequence];
102    double upperValue = upper[iSequence];
103    if (value - upperValue <= primalTolerance) {
104      if (value - lowerValue >= -primalTolerance) {
105        // feasible
106        status_[iSequence] = static_cast< unsigned char >(CLP_FEASIBLE | (CLP_SAME << 4));
107        bound_[iSequence] = 0.0;
108      } else {
109        // below
110        double infeasibility = lowerValue - value - primalTolerance;
111        sumInfeasibilities_ += infeasibility;
112        largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
113        cost[iSequence] -= infeasibilityCost;
114        numberInfeasibilities_++;
115        status_[iSequence] = static_cast< unsigned char >(CLP_BELOW_LOWER | (CLP_SAME << 4));
116        bound_[iSequence] = upperValue;
117        upper[iSequence] = lowerValue;
118        lower[iSequence] = -COIN_DBL_MAX;
119      }
120    } else {
121      // above
122      double infeasibility = value - upperValue - primalTolerance;
123      sumInfeasibilities_ += infeasibility;
124      largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
125      cost[iSequence] += infeasibilityCost;
126      numberInfeasibilities_++;
127      status_[iSequence] = static_cast< unsigned char >(CLP_ABOVE_UPPER | (CLP_SAME << 4));
128      bound_[iSequence] = lowerValue;
129      lower[iSequence] = upperValue;
130      upper[iSequence] = COIN_DBL_MAX;
131    }
132  }
133  //     checkInfeasibilities(model_->primalTolerance());
134}
135// Refresh - from original
136void AbcNonLinearCost::refreshFromPerturbed(double tolerance)
137{
138  // original costs and perturbed bounds
139  model_->copyFromSaved(32 + 2);
140  refresh();
141  //checkInfeasibilities(tolerance);
142}
143// Refreshes costs always makes row costs zero
144void AbcNonLinearCost::refreshCosts(const double *columnCosts)
145{
146  double *cost = model_->costRegion();
147  // zero row costs
148  memset(cost + numberColumns_, 0, numberRows_ * sizeof(double));
149  // copy column costs
150  CoinMemcpyN(columnCosts, numberColumns_, cost);
151  for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
152    cost_[iSequence] = cost[iSequence];
153  }
154}
155//-------------------------------------------------------------------
156// Copy constructor
157//-------------------------------------------------------------------
158AbcNonLinearCost::AbcNonLinearCost(const AbcNonLinearCost &rhs)
159  : changeCost_(0.0)
160  , feasibleCost_(0.0)
161  , infeasibilityWeight_(-1.0)
162  , largestInfeasibility_(0.0)
163  , sumInfeasibilities_(0.0)
164  , averageTheta_(0.0)
165  , numberRows_(rhs.numberRows_)
166  , numberColumns_(rhs.numberColumns_)
167  , model_(NULL)
168  , numberInfeasibilities_(-1)
169  , status_(NULL)
170  , bound_(NULL)
171  , cost_(NULL)
172{
173  if (numberRows_) {
174    int numberTotal = numberRows_ + numberColumns_;
175    model_ = rhs.model_;
176    numberInfeasibilities_ = rhs.numberInfeasibilities_;
177    changeCost_ = rhs.changeCost_;
178    feasibleCost_ = rhs.feasibleCost_;
179    infeasibilityWeight_ = rhs.infeasibilityWeight_;
180    largestInfeasibility_ = rhs.largestInfeasibility_;
181    sumInfeasibilities_ = rhs.sumInfeasibilities_;
182    averageTheta_ = rhs.averageTheta_;
183    bound_ = CoinCopyOfArray(rhs.bound_, numberTotal);
184    cost_ = CoinCopyOfArray(rhs.cost_, numberTotal);
185    status_ = CoinCopyOfArray(rhs.status_, numberTotal);
186  }
187}
188
189//-------------------------------------------------------------------
190// Destructor
191//-------------------------------------------------------------------
192AbcNonLinearCost::~AbcNonLinearCost()
193{
194  delete[] status_;
195  delete[] bound_;
196  delete[] cost_;
197}
198
199//----------------------------------------------------------------
200// Assignment operator
201//-------------------------------------------------------------------
202AbcNonLinearCost &
203AbcNonLinearCost::operator=(const AbcNonLinearCost &rhs)
204{
205  if (this != &rhs) {
206    numberRows_ = rhs.numberRows_;
207    numberColumns_ = rhs.numberColumns_;
208    delete[] status_;
209    delete[] bound_;
210    delete[] cost_;
211    status_ = NULL;
212    bound_ = NULL;
213    cost_ = NULL;
214    if (numberRows_) {
215      int numberTotal = numberRows_ + numberColumns_;
216      bound_ = CoinCopyOfArray(rhs.bound_, numberTotal);
217      cost_ = CoinCopyOfArray(rhs.cost_, numberTotal);
218      status_ = CoinCopyOfArray(rhs.status_, numberTotal);
219    }
220    model_ = rhs.model_;
221    numberInfeasibilities_ = rhs.numberInfeasibilities_;
222    changeCost_ = rhs.changeCost_;
223    feasibleCost_ = rhs.feasibleCost_;
224    infeasibilityWeight_ = rhs.infeasibilityWeight_;
225    largestInfeasibility_ = rhs.largestInfeasibility_;
226    sumInfeasibilities_ = rhs.sumInfeasibilities_;
227    averageTheta_ = rhs.averageTheta_;
228  }
229  return *this;
230}
231// Changes infeasible costs and computes number and cost of infeas
232// We will need to re-think objective offsets later
233// We will also need a 2 bit per variable array for some
234// purpose which will come to me later
235void AbcNonLinearCost::checkInfeasibilities(double oldTolerance)
236{
237  numberInfeasibilities_ = 0;
238  double infeasibilityCost = model_->infeasibilityCost();
239  changeCost_ = 0.0;
240  largestInfeasibility_ = 0.0;
241  sumInfeasibilities_ = 0.0;
242  double primalTolerance = model_->currentPrimalTolerance();
243  int iSequence;
244  double *solution = model_->solutionRegion();
245  double *upper = model_->upperRegion();
246  double *lower = model_->lowerRegion();
247  double *cost = model_->costRegion();
248  bool toNearest = oldTolerance <= 0.0;
249  feasibleCost_ = 0.0;
250  //bool checkCosts = (infeasibilityWeight_ != infeasibilityCost);
251  infeasibilityWeight_ = infeasibilityCost;
252  int numberTotal = numberColumns_ + numberRows_;
253  // nonbasic should be at a valid bound
254  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
255    double value = solution[iSequence];
256    unsigned char iStatus = status_[iSequence];
257    assert(currentStatus(iStatus) == CLP_SAME);
258    double lowerValue = lower[iSequence];
259    double upperValue = upper[iSequence];
260    double costValue = cost_[iSequence];
261    double trueCost = costValue;
262    int iWhere = originalStatus(iStatus);
263    if (iWhere == CLP_BELOW_LOWER) {
264      lowerValue = upperValue;
265      upperValue = bound_[iSequence];
266      costValue -= infeasibilityCost;
267    } else if (iWhere == CLP_ABOVE_UPPER) {
268      upperValue = lowerValue;
269      lowerValue = bound_[iSequence];
270      costValue += infeasibilityCost;
271    }
272    // get correct place
273    int newWhere = CLP_FEASIBLE;
274    AbcSimplex::Status status = model_->getInternalStatus(iSequence);
275    if (upperValue == lowerValue && status != AbcSimplex::isFixed) {
276      if (status != AbcSimplex::basic) {
277        model_->setInternalStatus(iSequence, AbcSimplex::isFixed);
278        status = AbcSimplex::isFixed;
279      }
280    }
281    switch (status) {
282
283    case AbcSimplex::basic:
284    case AbcSimplex::superBasic:
285      if (value - upperValue <= primalTolerance) {
286        if (value - lowerValue >= -primalTolerance) {
287          // feasible
288          //newWhere=CLP_FEASIBLE;
289        } else {
290          // below
291          newWhere = CLP_BELOW_LOWER;
292          assert(fabs(lowerValue) < 1.0e100);
293          double infeasibility = lowerValue - value - primalTolerance;
294          sumInfeasibilities_ += infeasibility;
295          largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
296          costValue = trueCost - infeasibilityCost;
297          changeCost_ -= lowerValue * (costValue - cost[iSequence]);
298          numberInfeasibilities_++;
299        }
300      } else {
301        // above
302        newWhere = CLP_ABOVE_UPPER;
303        double infeasibility = value - upperValue - primalTolerance;
304        sumInfeasibilities_ += infeasibility;
305        largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
306        costValue = trueCost + infeasibilityCost;
307        changeCost_ -= upperValue * (costValue - cost[iSequence]);
308        numberInfeasibilities_++;
309      }
310      break;
311    case AbcSimplex::isFree:
312      break;
313    case AbcSimplex::atUpperBound:
314      if (!toNearest) {
315        // With increasing tolerances - we may be at wrong place
316        if (fabs(value - upperValue) > oldTolerance * 1.0001) {
317          if (fabs(value - lowerValue) <= oldTolerance * 1.0001) {
318            if (fabs(value - lowerValue) > primalTolerance) {
319              solution[iSequence] = lowerValue;
320              value = lowerValue;
321            }
322            model_->setInternalStatus(iSequence, AbcSimplex::atLowerBound);
323          } else {
324            if (value < upperValue) {
325              if (value > lowerValue) {
326                model_->setInternalStatus(iSequence, AbcSimplex::superBasic);
327              } else {
328                // set to lower bound as infeasible
329                solution[iSequence] = lowerValue;
330                value = lowerValue;
331                model_->setInternalStatus(iSequence, AbcSimplex::atLowerBound);
332              }
333            } else {
334              // set to upper bound as infeasible
335              solution[iSequence] = upperValue;
336              value = upperValue;
337            }
338          }
339        } else if (fabs(value - upperValue) > primalTolerance) {
340          solution[iSequence] = upperValue;
341          value = upperValue;
342        }
343      } else {
344        // Set to nearest and make at bound
345        if (fabs(value - lowerValue) < fabs(value - upperValue)) {
346          solution[iSequence] = lowerValue;
347          value = lowerValue;
348          model_->setInternalStatus(iSequence, AbcSimplex::atLowerBound);
349        } else {
350          solution[iSequence] = upperValue;
351          value = upperValue;
352        }
353      }
354      break;
355    case AbcSimplex::atLowerBound:
356      if (!toNearest) {
357        // With increasing tolerances - we may be at wrong place
358        if (fabs(value - lowerValue) > oldTolerance * 1.0001) {
359          if (fabs(value - upperValue) <= oldTolerance * 1.0001) {
360            if (fabs(value - upperValue) > primalTolerance) {
361              solution[iSequence] = upperValue;
362              value = upperValue;
363            }
364            model_->setInternalStatus(iSequence, AbcSimplex::atUpperBound);
365          } else {
366            if (value < upperValue) {
367              if (value > lowerValue) {
368                model_->setInternalStatus(iSequence, AbcSimplex::superBasic);
369              } else {
370                // set to lower bound as infeasible
371                solution[iSequence] = lowerValue;
372                value = lowerValue;
373              }
374            } else {
375              // set to upper bound as infeasible
376              solution[iSequence] = upperValue;
377              value = upperValue;
378              model_->setInternalStatus(iSequence, AbcSimplex::atUpperBound);
379            }
380          }
381        } else if (fabs(value - lowerValue) > primalTolerance) {
382          solution[iSequence] = lowerValue;
383          value = lowerValue;
384        }
385      } else {
386        // Set to nearest and make at bound
387        if (fabs(value - lowerValue) < fabs(value - upperValue)) {
388          solution[iSequence] = lowerValue;
389          value = lowerValue;
390        } else {
391          solution[iSequence] = upperValue;
392          value = upperValue;
393          model_->setInternalStatus(iSequence, AbcSimplex::atUpperBound);
394        }
395      }
396      break;
397    case AbcSimplex::isFixed:
398      solution[iSequence] = lowerValue;
399      value = lowerValue;
400      break;
401    }
402    if (iWhere != newWhere) {
403      setOriginalStatus(status_[iSequence], newWhere);
404      if (newWhere == CLP_BELOW_LOWER) {
405        bound_[iSequence] = upperValue;
406        upperValue = lowerValue;
407        lowerValue = -COIN_DBL_MAX;
408        costValue = trueCost - infeasibilityCost;
409      } else if (newWhere == CLP_ABOVE_UPPER) {
410        bound_[iSequence] = lowerValue;
411        lowerValue = upperValue;
412        upperValue = COIN_DBL_MAX;
413        costValue = trueCost + infeasibilityCost;
414      } else {
415        costValue = trueCost;
416      }
417      lower[iSequence] = lowerValue;
418      upper[iSequence] = upperValue;
419    }
420    // always do as other things may change
421    cost[iSequence] = costValue;
422    feasibleCost_ += trueCost * value;
423  }
424  model_->moveToBasic(14); // all except solution
425}
426// Puts feasible bounds into lower and upper
427void AbcNonLinearCost::feasibleBounds()
428{
429  int iSequence;
430  double *upper = model_->upperRegion();
431  double *lower = model_->lowerRegion();
432  double *cost = model_->costRegion();
433  int numberTotal = numberColumns_ + numberRows_;
434  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
435    unsigned char iStatus = status_[iSequence];
436    assert(currentStatus(iStatus) == CLP_SAME);
437    double lowerValue = lower[iSequence];
438    double upperValue = upper[iSequence];
439    double costValue = cost_[iSequence];
440    int iWhere = originalStatus(iStatus);
441    if (iWhere == CLP_BELOW_LOWER) {
442      lowerValue = upperValue;
443      upperValue = bound_[iSequence];
444      assert(fabs(lowerValue) < 1.0e100);
445    } else if (iWhere == CLP_ABOVE_UPPER) {
446      upperValue = lowerValue;
447      lowerValue = bound_[iSequence];
448    }
449    setOriginalStatus(status_[iSequence], CLP_FEASIBLE);
450    lower[iSequence] = lowerValue;
451    upper[iSequence] = upperValue;
452    cost[iSequence] = costValue;
453  }
454}
455void AbcNonLinearCost::goBackAll(const CoinIndexedVector *update)
456{
457  assert(model_ != NULL);
458  const int *pivotVariable = model_->pivotVariable();
459  int number = update->getNumElements();
460  const int *index = update->getIndices();
461  for (int i = 0; i < number; i++) {
462    int iRow = index[i];
463    int iSequence = pivotVariable[iRow];
464    setSameStatus(status_[iSequence]);
465  }
466}
467void AbcNonLinearCost::checkInfeasibilities(int numberInArray, const int *index)
468{
469  assert(model_ != NULL);
470  double primalTolerance = model_->currentPrimalTolerance();
471  const int *pivotVariable = model_->pivotVariable();
472  double *upper = model_->upperRegion();
473  double *lower = model_->lowerRegion();
474  double *cost = model_->costRegion();
475  double *solutionBasic = model_->solutionBasic();
476  double *upperBasic = model_->upperBasic();
477  double *lowerBasic = model_->lowerBasic();
478  double *costBasic = model_->costBasic();
479  for (int i = 0; i < numberInArray; i++) {
480    int iRow = index[i];
481    int iSequence = pivotVariable[iRow];
482    double value = solutionBasic[iRow];
483    unsigned char iStatus = status_[iSequence];
484    assert(currentStatus(iStatus) == CLP_SAME);
485    double lowerValue = lowerBasic[iRow];
486    double upperValue = upperBasic[iRow];
487    double costValue = cost_[iSequence];
488    int iWhere = originalStatus(iStatus);
489    if (iWhere == CLP_BELOW_LOWER) {
490      lowerValue = upperValue;
491      upperValue = bound_[iSequence];
492      numberInfeasibilities_--;
493      assert(fabs(lowerValue) < 1.0e100);
494    } else if (iWhere == CLP_ABOVE_UPPER) {
495      upperValue = lowerValue;
496      lowerValue = bound_[iSequence];
497      numberInfeasibilities_--;
498    }
499    // get correct place
500    int newWhere = CLP_FEASIBLE;
501    if (value - upperValue <= primalTolerance) {
502      if (value - lowerValue >= -primalTolerance) {
503        // feasible
504        //newWhere=CLP_FEASIBLE;
505      } else {
506        // below
507        newWhere = CLP_BELOW_LOWER;
508        assert(fabs(lowerValue) < 1.0e100);
509        costValue -= infeasibilityWeight_;
510        numberInfeasibilities_++;
511      }
512    } else {
513      // above
514      newWhere = CLP_ABOVE_UPPER;
515      costValue += infeasibilityWeight_;
516      numberInfeasibilities_++;
517    }
518    if (iWhere != newWhere) {
519      setOriginalStatus(status_[iSequence], newWhere);
520      if (newWhere == CLP_BELOW_LOWER) {
521        bound_[iSequence] = upperValue;
522        upperValue = lowerValue;
523        lowerValue = -COIN_DBL_MAX;
524      } else if (newWhere == CLP_ABOVE_UPPER) {
525        bound_[iSequence] = lowerValue;
526        lowerValue = upperValue;
527        upperValue = COIN_DBL_MAX;
528      }
529      lower[iSequence] = lowerValue;
530      upper[iSequence] = upperValue;
531      cost[iSequence] = costValue;
532      lowerBasic[iRow] = lowerValue;
533      upperBasic[iRow] = upperValue;
534      costBasic[iRow] = costValue;
535    }
536  }
537}
538/* Puts back correct infeasible costs for each variable
539   The input indices are row indices and need converting to sequences
540   for costs.
541   On input array is empty (but indices exist).  On exit just
542   changed costs will be stored as normal CoinIndexedVector
543*/
544void AbcNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector *update)
545{
546  assert(model_ != NULL);
547  double primalTolerance = model_->currentPrimalTolerance();
548  const int *pivotVariable = model_->pivotVariable();
549  int number = 0;
550  int *index = update->getIndices();
551  double *work = update->denseVector();
552  double *upper = model_->upperRegion();
553  double *lower = model_->lowerRegion();
554  double *cost = model_->costRegion();
555  double *solutionBasic = model_->solutionBasic();
556  double *upperBasic = model_->upperBasic();
557  double *lowerBasic = model_->lowerBasic();
558  double *costBasic = model_->costBasic();
559  for (int i = 0; i < numberInArray; i++) {
560    int iRow = index[i];
561    int iSequence = pivotVariable[iRow];
562    double value = solutionBasic[iRow];
563    unsigned char iStatus = status_[iSequence];
564    assert(currentStatus(iStatus) == CLP_SAME);
565    double lowerValue = lowerBasic[iRow];
566    double upperValue = upperBasic[iRow];
567    double costValue = cost_[iSequence];
568    int iWhere = originalStatus(iStatus);
569    if (iWhere == CLP_BELOW_LOWER) {
570      lowerValue = upperValue;
571      upperValue = bound_[iSequence];
572      numberInfeasibilities_--;
573      assert(fabs(lowerValue) < 1.0e100);
574    } else if (iWhere == CLP_ABOVE_UPPER) {
575      upperValue = lowerValue;
576      lowerValue = bound_[iSequence];
577      numberInfeasibilities_--;
578    }
579    // get correct place
580    int newWhere = CLP_FEASIBLE;
581    if (value - upperValue <= primalTolerance) {
582      if (value - lowerValue >= -primalTolerance) {
583        // feasible
584        //newWhere=CLP_FEASIBLE;
585      } else {
586        // below
587        newWhere = CLP_BELOW_LOWER;
588        costValue -= infeasibilityWeight_;
589        numberInfeasibilities_++;
590        assert(fabs(lowerValue) < 1.0e100);
591      }
592    } else {
593      // above
594      newWhere = CLP_ABOVE_UPPER;
595      costValue += infeasibilityWeight_;
596      numberInfeasibilities_++;
597    }
598    if (iWhere != newWhere) {
599      work[iRow] = cost[iSequence] - costValue;
600      index[number++] = iRow;
601      setOriginalStatus(status_[iSequence], newWhere);
602      if (newWhere == CLP_BELOW_LOWER) {
603        bound_[iSequence] = upperValue;
604        upperValue = lowerValue;
605        lowerValue = -COIN_DBL_MAX;
606      } else if (newWhere == CLP_ABOVE_UPPER) {
607        bound_[iSequence] = lowerValue;
608        lowerValue = upperValue;
609        upperValue = COIN_DBL_MAX;
610      }
611      lower[iSequence] = lowerValue;
612      upper[iSequence] = upperValue;
613      cost[iSequence] = costValue;
614      lowerBasic[iRow] = lowerValue;
615      upperBasic[iRow] = upperValue;
616      costBasic[iRow] = costValue;
617    }
618  }
619  update->setNumElements(number);
620}
621/* Sets bounds and cost for one variable - returns change in cost*/
622double
623AbcNonLinearCost::setOne(int iSequence, double value)
624{
625  assert(model_ != NULL);
626  double primalTolerance = model_->currentPrimalTolerance();
627  // difference in cost
628  double difference = 0.0;
629  double *upper = model_->upperRegion();
630  double *lower = model_->lowerRegion();
631  double *cost = model_->costRegion();
632  unsigned char iStatus = status_[iSequence];
633  assert(currentStatus(iStatus) == CLP_SAME);
634  double lowerValue = lower[iSequence];
635  double upperValue = upper[iSequence];
636  double costValue = cost_[iSequence];
637  int iWhere = originalStatus(iStatus);
638  if (iWhere == CLP_BELOW_LOWER) {
639    lowerValue = upperValue;
640    upperValue = bound_[iSequence];
641    numberInfeasibilities_--;
642    assert(fabs(lowerValue) < 1.0e100);
643  } else if (iWhere == CLP_ABOVE_UPPER) {
644    upperValue = lowerValue;
645    lowerValue = bound_[iSequence];
646    numberInfeasibilities_--;
647  }
648  // get correct place
649  int newWhere = CLP_FEASIBLE;
650  if (value - upperValue <= primalTolerance) {
651    if (value - lowerValue >= -primalTolerance) {
652      // feasible
653      //newWhere=CLP_FEASIBLE;
654    } else {
655      // below
656      newWhere = CLP_BELOW_LOWER;
657      costValue -= infeasibilityWeight_;
658      numberInfeasibilities_++;
659      assert(fabs(lowerValue) < 1.0e100);
660    }
661  } else {
662    // above
663    newWhere = CLP_ABOVE_UPPER;
664    costValue += infeasibilityWeight_;
665    numberInfeasibilities_++;
666  }
667  if (iWhere != newWhere) {
668    difference = cost[iSequence] - costValue;
669    setOriginalStatus(status_[iSequence], newWhere);
670    if (newWhere == CLP_BELOW_LOWER) {
671      bound_[iSequence] = upperValue;
672      upperValue = lowerValue;
673      lowerValue = -COIN_DBL_MAX;
674    } else if (newWhere == CLP_ABOVE_UPPER) {
675      bound_[iSequence] = lowerValue;
676      lowerValue = upperValue;
677      upperValue = COIN_DBL_MAX;
678    }
679    lower[iSequence] = lowerValue;
680    upper[iSequence] = upperValue;
681    cost[iSequence] = costValue;
682  }
683  AbcSimplex::Status status = model_->getInternalStatus(iSequence);
684  if (upperValue == lowerValue) {
685    model_->setInternalStatus(iSequence, AbcSimplex::isFixed);
686  }
687  switch (status) {
688
689  case AbcSimplex::basic:
690  case AbcSimplex::superBasic:
691  case AbcSimplex::isFree:
692    break;
693  case AbcSimplex::atUpperBound:
694  case AbcSimplex::atLowerBound:
695  case AbcSimplex::isFixed:
696    // set correctly
697    if (fabs(value - lowerValue) <= primalTolerance * 1.001) {
698      model_->setInternalStatus(iSequence, AbcSimplex::atLowerBound);
699    } else if (fabs(value - upperValue) <= primalTolerance * 1.001) {
700      model_->setInternalStatus(iSequence, AbcSimplex::atUpperBound);
701    } else {
702      // set superBasic
703      model_->setInternalStatus(iSequence, AbcSimplex::superBasic);
704    }
705    break;
706  }
707  changeCost_ += value * difference;
708  return difference;
709}
710/* Sets bounds and cost for one variable - returns change in cost*/
711double
712AbcNonLinearCost::setOneBasic(int iRow, double value)
713{
714  assert(model_ != NULL);
715  int iSequence = model_->pivotVariable()[iRow];
716  double primalTolerance = model_->currentPrimalTolerance();
717  // difference in cost
718  double difference = 0.0;
719  double *upper = model_->upperRegion();
720  double *lower = model_->lowerRegion();
721  double *cost = model_->costRegion();
722  double *upperBasic = model_->upperBasic();
723  double *lowerBasic = model_->lowerBasic();
724  double *costBasic = model_->costBasic();
725  unsigned char iStatus = status_[iSequence];
726  assert(currentStatus(iStatus) == CLP_SAME);
727  double lowerValue = lowerBasic[iRow];
728  double upperValue = upperBasic[iRow];
729  double costValue = cost_[iSequence];
730  int iWhere = originalStatus(iStatus);
731  if (iWhere == CLP_BELOW_LOWER) {
732    lowerValue = upperValue;
733    upperValue = bound_[iSequence];
734    numberInfeasibilities_--;
735    assert(fabs(lowerValue) < 1.0e100);
736  } else if (iWhere == CLP_ABOVE_UPPER) {
737    upperValue = lowerValue;
738    lowerValue = bound_[iSequence];
739    numberInfeasibilities_--;
740  }
741  // get correct place
742  int newWhere = CLP_FEASIBLE;
743  if (value - upperValue <= primalTolerance) {
744    if (value - lowerValue >= -primalTolerance) {
745      // feasible
746      //newWhere=CLP_FEASIBLE;
747    } else {
748      // below
749      newWhere = CLP_BELOW_LOWER;
750      costValue -= infeasibilityWeight_;
751      numberInfeasibilities_++;
752      assert(fabs(lowerValue) < 1.0e100);
753    }
754  } else {
755    // above
756    newWhere = CLP_ABOVE_UPPER;
757    costValue += infeasibilityWeight_;
758    numberInfeasibilities_++;
759  }
760  if (iWhere != newWhere) {
761    difference = cost[iSequence] - costValue;
762    setOriginalStatus(status_[iSequence], newWhere);
763    if (newWhere == CLP_BELOW_LOWER) {
764      bound_[iSequence] = upperValue;
765      upperValue = lowerValue;
766      lowerValue = -COIN_DBL_MAX;
767    } else if (newWhere == CLP_ABOVE_UPPER) {
768      bound_[iSequence] = lowerValue;
769      lowerValue = upperValue;
770      upperValue = COIN_DBL_MAX;
771    }
772    lower[iSequence] = lowerValue;
773    upper[iSequence] = upperValue;
774    cost[iSequence] = costValue;
775    lowerBasic[iRow] = lowerValue;
776    upperBasic[iRow] = upperValue;
777    costBasic[iRow] = costValue;
778  }
779  changeCost_ += value * difference;
780  return difference;
781}
782/* Sets bounds and cost for outgoing variable
783   may change value
784   Returns direction */
785int AbcNonLinearCost::setOneOutgoing(int iRow, double &value)
786{
787  assert(model_ != NULL);
788  int iSequence = model_->pivotVariable()[iRow];
789  double primalTolerance = model_->currentPrimalTolerance();
790  // difference in cost
791  double difference = 0.0;
792  int direction = 0;
793  double *upper = model_->upperRegion();
794  double *lower = model_->lowerRegion();
795  double *cost = model_->costRegion();
796  double *upperBasic = model_->upperBasic();
797  double *lowerBasic = model_->lowerBasic();
798  unsigned char iStatus = status_[iSequence];
799  assert(currentStatus(iStatus) == CLP_SAME);
800  double lowerValue = lowerBasic[iRow];
801  double upperValue = upperBasic[iRow];
802  double costValue = cost_[iSequence];
803  // Set perceived direction out
804  if (value <= lowerValue + 1.001 * primalTolerance) {
805    direction = 1;
806  } else if (value >= upperValue - 1.001 * primalTolerance) {
807    direction = -1;
808  } else {
809    // odd
810    direction = 0;
811  }
812  int iWhere = originalStatus(iStatus);
813  if (iWhere == CLP_BELOW_LOWER) {
814    lowerValue = upperValue;
815    upperValue = bound_[iSequence];
816    numberInfeasibilities_--;
817    assert(fabs(lowerValue) < 1.0e100);
818  } else if (iWhere == CLP_ABOVE_UPPER) {
819    upperValue = lowerValue;
820    lowerValue = bound_[iSequence];
821    numberInfeasibilities_--;
822  }
823  // get correct place
824  // If fixed give benefit of doubt
825  if (lowerValue == upperValue)
826    value = lowerValue;
827  int newWhere = CLP_FEASIBLE;
828  if (value - upperValue <= primalTolerance) {
829    if (value - lowerValue >= -primalTolerance) {
830      // feasible
831      //newWhere=CLP_FEASIBLE;
832    } else {
833      // below
834      newWhere = CLP_BELOW_LOWER;
835      costValue -= infeasibilityWeight_;
836      numberInfeasibilities_++;
837      assert(fabs(lowerValue) < 1.0e100);
838    }
839  } else {
840    // above
841    newWhere = CLP_ABOVE_UPPER;
842    costValue += infeasibilityWeight_;
843    numberInfeasibilities_++;
844  }
845  if (iWhere != newWhere) {
846    difference = cost[iSequence] - costValue;
847    setOriginalStatus(status_[iSequence], newWhere);
848    if (newWhere == CLP_BELOW_LOWER) {
849      bound_[iSequence] = upperValue;
850      upper[iSequence] = lowerValue;
851      lower[iSequence] = -COIN_DBL_MAX;
852    } else if (newWhere == CLP_ABOVE_UPPER) {
853      bound_[iSequence] = lowerValue;
854      lower[iSequence] = upperValue;
855      upper[iSequence] = COIN_DBL_MAX;
856    } else {
857      lower[iSequence] = lowerValue;
858      upper[iSequence] = upperValue;
859    }
860    cost[iSequence] = costValue;
861  }
862  // set correctly
863  if (fabs(value - lowerValue) <= primalTolerance * 1.001) {
864    value = CoinMin(value, lowerValue + primalTolerance);
865  } else if (fabs(value - upperValue) <= primalTolerance * 1.001) {
866    value = CoinMax(value, upperValue - primalTolerance);
867  } else {
868    //printf("*** variable wandered off bound %g %g %g!\n",
869    //     lowerValue,value,upperValue);
870    if (value - lowerValue <= upperValue - value)
871      value = lowerValue + primalTolerance;
872    else
873      value = upperValue - primalTolerance;
874  }
875  changeCost_ += value * difference;
876  return direction;
877}
878// Returns nearest bound
879double
880AbcNonLinearCost::nearest(int iRow, double solutionValue)
881{
882  assert(model_ != NULL);
883  int iSequence = model_->pivotVariable()[iRow];
884  double nearest = 0.0;
885  const double *upperBasic = model_->upperBasic();
886  const double *lowerBasic = model_->lowerBasic();
887  double lowerValue = lowerBasic[iRow];
888  double upperValue = upperBasic[iRow];
889  int iWhere = originalStatus(status_[iSequence]);
890  if (iWhere == CLP_BELOW_LOWER) {
891    lowerValue = upperValue;
892    upperValue = bound_[iSequence];
893    assert(fabs(lowerValue) < 1.0e100);
894  } else if (iWhere == CLP_ABOVE_UPPER) {
895    upperValue = lowerValue;
896    lowerValue = bound_[iSequence];
897  }
898  if (fabs(solutionValue - lowerValue) < fabs(solutionValue - upperValue))
899    nearest = lowerValue;
900  else
901    nearest = upperValue;
902  return nearest;
903}
904/// Feasible cost with offset and direction (i.e. for reporting)
905double
906AbcNonLinearCost::feasibleReportCost() const
907{
908  double value;
909  model_->getDblParam(ClpObjOffset, value);
910  return (feasibleCost_ + model_->objectiveAsObject()->nonlinearOffset()) * model_->optimizationDirection() / (model_->objectiveScale() * model_->rhsScale()) - value;
911}
912// Get rid of real costs (just for moment)
913void AbcNonLinearCost::zapCosts()
914{
915}
916#ifdef VALIDATE
917// For debug
918void AbcNonLinearCost::validate()
919{
920  double primalTolerance = model_->currentPrimalTolerance();
921  int iSequence;
922  const double *solution = model_->solutionRegion();
923  const double *upper = model_->upperRegion();
924  const double *lower = model_->lowerRegion();
925  const double *cost = model_->costRegion();
926  double infeasibilityCost = model_->infeasibilityCost();
927  int numberTotal = numberRows_ + numberColumns_;
928  int numberInfeasibilities = 0;
929  double sumInfeasibilities = 0.0;
930
931  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
932    double value = solution[iSequence];
933    int iStatus = status_[iSequence];
934    assert(currentStatus(iStatus) == CLP_SAME);
935    double lowerValue = lower[iSequence];
936    double upperValue = upper[iSequence];
937    double costValue = cost_[iSequence];
938    int iWhere = originalStatus(iStatus);
939    if (iWhere == CLP_BELOW_LOWER) {
940      lowerValue = upperValue;
941      upperValue = bound_[iSequence];
942      assert(fabs(lowerValue) < 1.0e100);
943      costValue -= infeasibilityCost;
944      assert(value <= lowerValue - primalTolerance);
945      numberInfeasibilities++;
946      sumInfeasibilities += lowerValue - value - primalTolerance;
947      assert(model_->getInternalStatus(iSequence) == AbcSimplex::basic);
948    } else if (iWhere == CLP_ABOVE_UPPER) {
949      upperValue = lowerValue;
950      lowerValue = bound_[iSequence];
951      costValue += infeasibilityCost;
952      assert(value >= upperValue + primalTolerance);
953      numberInfeasibilities++;
954      sumInfeasibilities += value - upperValue - primalTolerance;
955      assert(model_->getInternalStatus(iSequence) == AbcSimplex::basic);
956    } else {
957      assert(value >= lowerValue - primalTolerance && value <= upperValue + primalTolerance);
958    }
959    assert(lowerValue == saveLowerV[iSequence]);
960    assert(upperValue == saveUpperV[iSequence]);
961    assert(costValue == cost[iSequence]);
962  }
963  if (numberInfeasibilities)
964    printf("JJ %d infeasibilities summing to %g\n",
965      numberInfeasibilities, sumInfeasibilities);
966}
967#endif
968
969/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
970*/
Note: See TracBrowser for help on using the repository browser.