All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
31 #include <opm/common/ErrorMacros.hpp>
32 #include <opm/common/Exceptions.hpp>
33 #include <opm/common/Unused.hpp>
34 
35 #include <algorithm>
36 #include <cassert>
37 #include <iostream>
38 #include <tuple>
39 #include <vector>
40 
41 namespace Opm {
46 template <class Scalar>
48 {
49 public:
56  {}
57 
65  template <class ScalarArrayX, class ScalarArrayY>
66  Tabulated1DFunction(size_t nSamples,
67  const ScalarArrayX& x,
68  const ScalarArrayY& y,
69  bool sortInputs = true)
70  { this->setXYArrays(nSamples, x, y, sortInputs); }
71 
80  template <class ScalarContainer>
81  Tabulated1DFunction(const ScalarContainer& x,
82  const ScalarContainer& y,
83  bool sortInputs = true)
84  { this->setXYContainers(x, y, sortInputs); }
85 
92  template <class PointContainer>
93  Tabulated1DFunction(const PointContainer& points,
94  bool sortInputs = true)
95  { this->setContainerOfTuples(points, sortInputs); }
96 
102  template <class ScalarArrayX, class ScalarArrayY>
103  void setXYArrays(size_t nSamples,
104  const ScalarArrayX& x,
105  const ScalarArrayY& y,
106  bool sortInputs = true)
107  {
108  assert(nSamples > 1);
109 
110  resizeArrays_(nSamples);
111  for (size_t i = 0; i < nSamples; ++i) {
112  xValues_[i] = x[i];
113  yValues_[i] = y[i];
114  }
115 
116  if (sortInputs)
117  sortInput_();
118  else if (xValues_[0] > xValues_[numSamples() - 1])
119  reverseSamplingPoints_();
120  }
121 
127  template <class ScalarContainerX, class ScalarContainerY>
128  void setXYContainers(const ScalarContainerX& x,
129  const ScalarContainerY& y,
130  bool sortInputs = true)
131  {
132  assert(x.size() == y.size());
133  assert(x.size() > 1);
134 
135  resizeArrays_(x.size());
136  std::copy(x.begin(), x.end(), xValues_.begin());
137  std::copy(y.begin(), y.end(), yValues_.begin());
138 
139  if (sortInputs)
140  sortInput_();
141  else if (xValues_[0] > xValues_[numSamples() - 1])
142  reverseSamplingPoints_();
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 
232  Scalar valueAt(size_t i) const
233  { return yValues_[i]; }
234 
238  template <class Evaluation>
239  bool applies(const Evaluation& x) const
240  { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
241 
251  template <class Evaluation>
252  Evaluation eval(const Evaluation& x, bool extrapolate = false) const
253  {
254  if (!extrapolate && !applies(x))
255  OPM_THROW(Opm::NumericalProblem,
256  "Tried to evaluate a tabulated function outside of its range");
257 
258  size_t segIdx;
259 
260  if (extrapolate && x < xValues_.front())
261  segIdx = 0;
262  else if (extrapolate && x > xValues_.back())
263  segIdx = numSamples() - 2;
264  else
265  segIdx = findSegmentIndex_(Opm::scalarValue(x));
266 
267  Scalar x0 = xValues_[segIdx];
268  Scalar x1 = xValues_[segIdx + 1];
269 
270  Scalar y0 = yValues_[segIdx];
271  Scalar y1 = yValues_[segIdx + 1];
272 
273  return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
274  }
275 
287  template <class Evaluation>
288  Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
289  {
290  if (!extrapolate && !applies(x)) {
291  OPM_THROW(Opm::NumericalProblem,
292  "Tried to evaluate a derivative of a tabulated"
293  " function outside of its range");
294  }
295 
296  unsigned segIdx = findSegmentIndex_(Opm::scalarValue(x));
297  return evalDerivative_(x, segIdx);
298  }
299 
314  template <class Evaluation>
315  Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
316  {
317  if (!extrapolate && !applies(x)) {
318  OPM_THROW(Opm::NumericalProblem,
319  "Tried to evaluate a second derivative of a tabulated "
320  " function outside of its range");
321  }
322 
323  return 0.0;
324  }
325 
340  template <class Evaluation>
341  Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
342  {
343  if (!extrapolate && !applies(x)) {
344  OPM_THROW(Opm::NumericalProblem,
345  "Tried to evaluate a third derivative of a tabulated "
346  " function outside of its range");
347  }
348 
349  return 0.0;
350  }
351 
360  int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED = false) const
361  {
362  assert(x0 != x1);
363 
364  // make sure that x0 is smaller than x1
365  if (x0 > x1)
366  std::swap(x0, x1);
367 
368  assert(x0 < x1);
369 
370  int r = 3;
371  if (x0 < xMin()) {
372  assert(extrapolate);
373 
374  x0 = xMin();
375  };
376 
377  size_t i = findSegmentIndex_(x0);
378  if (xValues_[i + 1] >= x1) {
379  // interval is fully contained within a single function
380  // segment
381  updateMonotonicity_(i, r);
382  return r;
383  }
384 
385  // the first segment overlaps with the specified interval
386  // partially
387  updateMonotonicity_(i, r);
388  ++ i;
389 
390  // make sure that the segments which are completly in the
391  // interval [x0, x1] all exhibit the same monotonicity.
392  size_t iEnd = findSegmentIndex_(x1);
393  for (; i < iEnd - 1; ++i) {
394  updateMonotonicity_(i, r);
395  if (!r)
396  return 0;
397  }
398 
399  // if the user asked for a part of the function which is
400  // extrapolated, we need to check the slope at the function's
401  // endpoint
402  if (x1 > xMax()) {
403  assert(extrapolate);
404 
405  Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2);
406  if (m < 0)
407  return (r < 0 || r==3)?-1:0;
408  else if (m > 0)
409  return (r > 0 || r==3)?1:0;
410 
411  return r;
412  }
413 
414  // check for the last segment
415  updateMonotonicity_(iEnd, r);
416 
417  return r;
418  }
419 
424  int monotonic() const
425  { return monotonic(xMin(), xMax()); }
426 
443  void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os = std::cout) const
444  {
445  Scalar x0 = std::min(xi0, xi1);
446  Scalar x1 = std::max(xi0, xi1);
447  const int n = numSamples() - 1;
448  for (int i = 0; i <= k; ++i) {
449  double x = i*(x1 - x0)/k + x0;
450  double y;
451  double dy_dx;
452  if (!applies(x)) {
453  if (x < xValues_[0]) {
454  dy_dx = evalDerivative(xValues_[0]);
455  y = (x - xValues_[0])*dy_dx + yValues_[0];
456  }
457  else if (x > xValues_[n]) {
458  dy_dx = evalDerivative(xValues_[n]);
459  y = (x - xValues_[n])*dy_dx + yValues_[n];
460  }
461  else {
462  OPM_THROW(std::runtime_error,
463  "The sampling points given to a function must be sorted by their x value!");
464  }
465  }
466  else {
467  y = eval(x);
468  dy_dx = evalDerivative(x);
469  }
470 
471  os << x << " " << y << " " << dy_dx << "\n";
472  }
473  }
474 
475 private:
476  size_t findSegmentIndex_(Scalar x) const
477  {
478  // we need at least two sampling points!
479  assert(xValues_.size() >= 2);
480 
481  if (x <= xValues_[1])
482  return 0;
483  else if (x >= xValues_[xValues_.size() - 2])
484  return xValues_.size() - 2;
485  else {
486  // bisection
487  size_t segmentIdx = 1;
488  size_t upperIdx = xValues_.size() - 2;
489  while (segmentIdx + 1 < upperIdx) {
490  size_t pivotIdx = (segmentIdx + upperIdx) / 2;
491  if (x < xValues_[pivotIdx])
492  upperIdx = pivotIdx;
493  else
494  segmentIdx = pivotIdx;
495  }
496 
497  assert(xValues_[segmentIdx] <= x);
498  assert(x <= xValues_[segmentIdx + 1]);
499  return segmentIdx;
500  }
501  }
502 
503  template <class Evaluation>
504  Evaluation evalDerivative_(const Evaluation& x OPM_UNUSED, size_t segIdx) const
505  {
506  Scalar x0 = xValues_[segIdx];
507  Scalar x1 = xValues_[segIdx + 1];
508 
509  Scalar y0 = yValues_[segIdx];
510  Scalar y1 = yValues_[segIdx + 1];
511 
512  return (y1 - y0)/(x1 - x0);
513  }
514 
515  // returns the monotonicity of a segment
516  //
517  // The return value have the following meaning:
518  //
519  // 3: function is constant within interval [x0, x1]
520  // 1: function is monotonously increasing in the specified interval
521  // 0: function is not monotonic in the specified interval
522  // -1: function is monotonously decreasing in the specified interval
523  int updateMonotonicity_(size_t i, int& r) const
524  {
525  if (yValues_[i] < yValues_[i + 1]) {
526  // monotonically increasing?
527  if (r == 3 || r == 1)
528  r = 1;
529  else
530  r = 0;
531  return 1;
532  }
533  else if (yValues_[i] > yValues_[i + 1]) {
534  // monotonically decreasing?
535  if (r == 3 || r == -1)
536  r = -1;
537  else
538  r = 0;
539  return -1;
540  }
541 
542  return 3;
543  }
544 
548  struct ComparatorX_
549  {
550  ComparatorX_(const std::vector<Scalar>& x)
551  : x_(x)
552  {}
553 
554  bool operator ()(size_t idxA, size_t idxB) const
555  { return x_.at(idxA) < x_.at(idxB); }
556 
557  const std::vector<Scalar>& x_;
558  };
559 
563  void sortInput_()
564  {
565  size_t n = numSamples();
566 
567  // create a vector containing 0...n-1
568  std::vector<unsigned> idxVector(n);
569  for (unsigned i = 0; i < n; ++i)
570  idxVector[i] = i;
571 
572  // sort the indices according to the x values of the sample
573  // points
574  ComparatorX_ cmp(xValues_);
575  std::sort(idxVector.begin(), idxVector.end(), cmp);
576 
577  // reorder the sample points
578  std::vector<Scalar> tmpX(n), tmpY(n);
579  for (size_t i = 0; i < idxVector.size(); ++ i) {
580  tmpX[i] = xValues_[idxVector[i]];
581  tmpY[i] = yValues_[idxVector[i]];
582  }
583  xValues_ = tmpX;
584  yValues_ = tmpY;
585  }
586 
591  void reverseSamplingPoints_()
592  {
593  // reverse the arrays
594  size_t n = numSamples();
595  for (size_t i = 0; i <= (n - 1)/2; ++i) {
596  std::swap(xValues_[i], xValues_[n - i - 1]);
597  std::swap(yValues_[i], yValues_[n - i - 1]);
598  }
599  }
600 
604  void resizeArrays_(size_t nSamples)
605  {
606  xValues_.resize(nSamples);
607  yValues_.resize(nSamples);
608  }
609 
610  std::vector<Scalar> xValues_;
611  std::vector<Scalar> yValues_;
612 };
613 } // namespace Opm
614 
615 #endif
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval. ...
Definition: Tabulated1DFunction.hpp:424
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline&#39;s derivative at a given position.
Definition: Tabulated1DFunction.hpp:288
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:360
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:128
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:226
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the function&#39;s third derivative at a given position.
Definition: Tabulated1DFunction.hpp:341
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:103
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:55
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:239
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:220
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:93
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:214
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the function&#39;s second derivative at a given position.
Definition: Tabulated1DFunction.hpp:315
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:252
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:66
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
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:443
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:208
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:232
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:149
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:81