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

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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