My Project
Tabulated1DFunction.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_TABULATED_1D_FUNCTION_HPP
28 #define OPM_TABULATED_1D_FUNCTION_HPP
29 
33 
34 #include <algorithm>
35 #include <cassert>
36 #include <iostream>
37 #include <tuple>
38 #include <vector>
39 
40 namespace Opm {
45 template <class Scalar>
47 {
48 public:
55  {}
56 
64  template <class ScalarArrayX, class ScalarArrayY>
65  Tabulated1DFunction(size_t nSamples,
66  const ScalarArrayX& x,
67  const ScalarArrayY& y,
68  bool sortInputs = true)
69  { this->setXYArrays(nSamples, x, y, sortInputs); }
70 
79  template <class ScalarContainer>
80  Tabulated1DFunction(const ScalarContainer& x,
81  const ScalarContainer& y,
82  bool sortInputs = true)
83  { this->setXYContainers(x, y, sortInputs); }
84 
91  template <class PointContainer>
92  Tabulated1DFunction(const PointContainer& points,
93  bool sortInputs = true)
94  { this->setContainerOfTuples(points, sortInputs); }
95 
101  template <class ScalarArrayX, class ScalarArrayY>
102  void setXYArrays(size_t nSamples,
103  const ScalarArrayX& x,
104  const ScalarArrayY& y,
105  bool sortInputs = true)
106  {
107  assert(nSamples > 1);
108 
109  resizeArrays_(nSamples);
110  for (size_t i = 0; i < nSamples; ++i) {
111  xValues_[i] = x[i];
112  yValues_[i] = y[i];
113  }
114 
115  if (sortInputs)
116  sortInput_();
117  else if (xValues_[0] > xValues_[numSamples() - 1])
118  reverseSamplingPoints_();
119  }
120 
126  template <class ScalarContainerX, class ScalarContainerY>
127  void setXYContainers(const ScalarContainerX& x,
128  const ScalarContainerY& y,
129  bool sortInputs = true)
130  {
131  assert(x.size() == y.size());
132 
133  resizeArrays_(x.size());
134  if (x.size() > 0) {
135  std::copy(x.begin(), x.end(), xValues_.begin());
136  std::copy(y.begin(), y.end(), yValues_.begin());
137 
138  if (sortInputs)
139  sortInput_();
140  else if (xValues_[0] > xValues_[numSamples() - 1])
141  reverseSamplingPoints_();
142  }
143  }
144 
148  template <class PointArray>
149  void setArrayOfPoints(size_t nSamples,
150  const PointArray& points,
151  bool sortInputs = true)
152  {
153  // a linear function with less than two sampling points? what an incredible
154  // bad idea!
155  assert(nSamples > 1);
156 
157  resizeArrays_(nSamples);
158  for (size_t i = 0; i < nSamples; ++i) {
159  xValues_[i] = points[i][0];
160  yValues_[i] = points[i][1];
161  }
162 
163  if (sortInputs)
164  sortInput_();
165  else if (xValues_[0] > xValues_[numSamples() - 1])
166  reverseSamplingPoints_();
167  }
168 
183  template <class XYContainer>
184  void setContainerOfTuples(const XYContainer& points,
185  bool sortInputs = true)
186  {
187  // a linear function with less than two sampling points? what an incredible
188  // bad idea!
189  assert(points.size() > 1);
190 
191  resizeArrays_(points.size());
192  typename XYContainer::const_iterator it = points.begin();
193  typename XYContainer::const_iterator endIt = points.end();
194  for (unsigned i = 0; it != endIt; ++i, ++it) {
195  xValues_[i] = std::get<0>(*it);
196  yValues_[i] = std::get<1>(*it);
197  }
198 
199  if (sortInputs)
200  sortInput_();
201  else if (xValues_[0] > xValues_[numSamples() - 1])
202  reverseSamplingPoints_();
203  }
204 
208  size_t numSamples() const
209  { return xValues_.size(); }
210 
214  Scalar xMin() const
215  { return xValues_[0]; }
216 
220  Scalar xMax() const
221  { return xValues_[numSamples() - 1]; }
222 
226  Scalar xAt(size_t i) const
227  { return xValues_[i]; }
228 
229  const std::vector<Scalar>& xValues() const
230  { return xValues_; }
231 
232  const std::vector<Scalar>& yValues() const
233  { return yValues_; }
234 
238  Scalar valueAt(size_t i) const
239  { return yValues_[i]; }
240 
244  template <class Evaluation>
245  bool applies(const Evaluation& x) const
246  { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
247 
257  template <class Evaluation>
258  Evaluation eval(const Evaluation& x, bool extrapolate = false) const
259  {
260  size_t segIdx = findSegmentIndex_(x, extrapolate);
261 
262  Scalar x0 = xValues_[segIdx];
263  Scalar x1 = xValues_[segIdx + 1];
264 
265  Scalar y0 = yValues_[segIdx];
266  Scalar y1 = yValues_[segIdx + 1];
267 
268  return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
269  }
270 
282  template <class Evaluation>
283  Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
284  {
285  unsigned segIdx = findSegmentIndex_(x, extrapolate);
286  return evalDerivative_(x, segIdx);
287  }
288 
303  template <class Evaluation>
304  Evaluation evalSecondDerivative(const Evaluation&, bool = false) const
305  { return 0.0; }
306 
321  template <class Evaluation>
322  Evaluation evalThirdDerivative(const Evaluation&, bool = false) const
323  { return 0.0; }
324 
333  int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED = false) const
334  {
335  assert(x0 != x1);
336 
337  // make sure that x0 is smaller than x1
338  if (x0 > x1)
339  std::swap(x0, x1);
340 
341  assert(x0 < x1);
342 
343  int r = 3;
344  if (x0 < xMin()) {
345  assert(extrapolate);
346 
347  x0 = xMin();
348  };
349 
350  size_t i = findSegmentIndex_(x0, extrapolate);
351  if (xValues_[i + 1] >= x1) {
352  // interval is fully contained within a single function
353  // segment
354  updateMonotonicity_(i, r);
355  return r;
356  }
357 
358  // the first segment overlaps with the specified interval
359  // partially
360  updateMonotonicity_(i, r);
361  ++ i;
362 
363  // make sure that the segments which are completly in the
364  // interval [x0, x1] all exhibit the same monotonicity.
365  size_t iEnd = findSegmentIndex_(x1, extrapolate);
366  for (; i < iEnd - 1; ++i) {
367  updateMonotonicity_(i, r);
368  if (!r)
369  return 0;
370  }
371 
372  // if the user asked for a part of the function which is
373  // extrapolated, we need to check the slope at the function's
374  // endpoint
375  if (x1 > xMax()) {
376  assert(extrapolate);
377 
378  Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2);
379  if (m < 0)
380  return (r < 0 || r==3)?-1:0;
381  else if (m > 0)
382  return (r > 0 || r==3)?1:0;
383 
384  return r;
385  }
386 
387  // check for the last segment
388  updateMonotonicity_(iEnd, r);
389 
390  return r;
391  }
392 
397  int monotonic() const
398  { return monotonic(xMin(), xMax()); }
399 
416  void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os = std::cout) const
417  {
418  Scalar x0 = std::min(xi0, xi1);
419  Scalar x1 = std::max(xi0, xi1);
420  const int n = numSamples() - 1;
421  for (unsigned i = 0; i <= k; ++i) {
422  double x = i*(x1 - x0)/k + x0;
423  double y;
424  double dy_dx;
425  if (!applies(x)) {
426  if (x < xValues_[0]) {
427  dy_dx = evalDerivative(xValues_[0]);
428  y = (x - xValues_[0])*dy_dx + yValues_[0];
429  }
430  else if (x > xValues_[n]) {
431  dy_dx = evalDerivative(xValues_[n]);
432  y = (x - xValues_[n])*dy_dx + yValues_[n];
433  }
434  else {
435  throw std::runtime_error("The sampling points given to a function must be sorted by their x value!");
436  }
437  }
438  else {
439  y = eval(x);
440  dy_dx = evalDerivative(x);
441  }
442 
443  os << x << " " << y << " " << dy_dx << "\n";
444  }
445  }
446 
447  bool operator==(const Tabulated1DFunction<Scalar>& data) const {
448  return xValues_ == data.xValues_ &&
449  yValues_ == data.yValues_;
450  }
451 
452 private:
453  template <class Evaluation>
454  size_t findSegmentIndex_(const Evaluation& x, bool extrapolate = false) const
455  {
456  if (!extrapolate && !applies(x))
457  throw NumericalIssue("Tried to evaluate a tabulated function outside of its range");
458 
459  // we need at least two sampling points!
460  assert(xValues_.size() >= 2);
461 
462  if (x <= xValues_[1])
463  return 0;
464  else if (x >= xValues_[xValues_.size() - 2])
465  return xValues_.size() - 2;
466  else {
467  // bisection
468  size_t lowerIdx = 1;
469  size_t upperIdx = xValues_.size() - 2;
470  while (lowerIdx + 1 < upperIdx) {
471  size_t pivotIdx = (lowerIdx + upperIdx) / 2;
472  if (x < xValues_[pivotIdx])
473  upperIdx = pivotIdx;
474  else
475  lowerIdx = pivotIdx;
476  }
477 
478  assert(xValues_[lowerIdx] <= x);
479  assert(x <= xValues_[lowerIdx + 1]);
480  return lowerIdx;
481  }
482  }
483 
484  template <class Evaluation>
485  Evaluation evalDerivative_(const Evaluation& x, size_t segIdx) const
486  {
487  Scalar x0 = xValues_[segIdx];
488  Scalar x1 = xValues_[segIdx + 1];
489 
490  Scalar y0 = yValues_[segIdx];
491  Scalar y1 = yValues_[segIdx + 1];
492 
493  Evaluation ret = blank(x);
494  ret = (y1 - y0)/(x1 - x0);
495  return ret;
496  }
497 
498  // returns the monotonicity of a segment
499  //
500  // The return value have the following meaning:
501  //
502  // 3: function is constant within interval [x0, x1]
503  // 1: function is monotonously increasing in the specified interval
504  // 0: function is not monotonic in the specified interval
505  // -1: function is monotonously decreasing in the specified interval
506  int updateMonotonicity_(size_t i, int& r) const
507  {
508  if (yValues_[i] < yValues_[i + 1]) {
509  // monotonically increasing?
510  if (r == 3 || r == 1)
511  r = 1;
512  else
513  r = 0;
514  return 1;
515  }
516  else if (yValues_[i] > yValues_[i + 1]) {
517  // monotonically decreasing?
518  if (r == 3 || r == -1)
519  r = -1;
520  else
521  r = 0;
522  return -1;
523  }
524 
525  return 3;
526  }
527 
531  struct ComparatorX_
532  {
533  ComparatorX_(const std::vector<Scalar>& x)
534  : x_(x)
535  {}
536 
537  bool operator ()(size_t idxA, size_t idxB) const
538  { return x_.at(idxA) < x_.at(idxB); }
539 
540  const std::vector<Scalar>& x_;
541  };
542 
546  void sortInput_()
547  {
548  size_t n = numSamples();
549 
550  // create a vector containing 0...n-1
551  std::vector<unsigned> idxVector(n);
552  for (unsigned i = 0; i < n; ++i)
553  idxVector[i] = i;
554 
555  // sort the indices according to the x values of the sample
556  // points
557  ComparatorX_ cmp(xValues_);
558  std::sort(idxVector.begin(), idxVector.end(), cmp);
559 
560  // reorder the sample points
561  std::vector<Scalar> tmpX(n), tmpY(n);
562  for (size_t i = 0; i < idxVector.size(); ++ i) {
563  tmpX[i] = xValues_[idxVector[i]];
564  tmpY[i] = yValues_[idxVector[i]];
565  }
566  xValues_ = tmpX;
567  yValues_ = tmpY;
568  }
569 
574  void reverseSamplingPoints_()
575  {
576  // reverse the arrays
577  size_t n = numSamples();
578  for (size_t i = 0; i <= (n - 1)/2; ++i) {
579  std::swap(xValues_[i], xValues_[n - i - 1]);
580  std::swap(yValues_[i], yValues_[n - i - 1]);
581  }
582  }
583 
587  void resizeArrays_(size_t nSamples)
588  {
589  xValues_.resize(nSamples);
590  yValues_.resize(nSamples);
591  }
592 
593  std::vector<Scalar> xValues_;
594  std::vector<Scalar> yValues_;
595 };
596 } // namespace Opm
597 
598 #endif
Provides the opm-material specific exception classes.
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Provides the OPM_UNUSED macro.
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:258
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:208
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Tabulated1DFunction.hpp:283
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:149
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:220
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:80
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:214
Evaluation evalThirdDerivative(const Evaluation &, bool=false) const
Evaluate the function's third derivative at a given position.
Definition: Tabulated1DFunction.hpp:322
int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED=false) const
Returns 1 if the function is monotonically increasing, -1 if the function is mononously decreasing an...
Definition: Tabulated1DFunction.hpp:333
Evaluation evalSecondDerivative(const Evaluation &, bool=false) const
Evaluate the function's second derivative at a given position.
Definition: Tabulated1DFunction.hpp:304
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:65
void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream &os=std::cout) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition: Tabulated1DFunction.hpp:416
void setContainerOfTuples(const XYContainer &points, bool sortInputs=true)
Set the sampling points of the piecewise linear function using a STL-compatible container of tuple-li...
Definition: Tabulated1DFunction.hpp:184
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:226
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:92
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:54
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval.
Definition: Tabulated1DFunction.hpp:397
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:102
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:127
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:245
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:238