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

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