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 
30 #include <opm/common/OpmLog/OpmLog.hpp>
33 
34 #include <algorithm>
35 #include <cassert>
36 #include <iostream>
37 #include <tuple>
38 #include <vector>
39 #include <sstream>
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 
134  resizeArrays_(x.size());
135  if (x.size() > 0) {
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  }
145 
149  template <class PointArray>
150  void setArrayOfPoints(size_t nSamples,
151  const PointArray& points,
152  bool sortInputs = true)
153  {
154  // a linear function with less than two sampling points? what an incredible
155  // bad idea!
156  assert(nSamples > 1);
157 
158  resizeArrays_(nSamples);
159  for (size_t i = 0; i < nSamples; ++i) {
160  xValues_[i] = points[i][0];
161  yValues_[i] = points[i][1];
162  }
163 
164  if (sortInputs)
165  sortInput_();
166  else if (xValues_[0] > xValues_[numSamples() - 1])
167  reverseSamplingPoints_();
168  }
169 
184  template <class XYContainer>
185  void setContainerOfTuples(const XYContainer& points,
186  bool sortInputs = true)
187  {
188  // a linear function with less than two sampling points? what an incredible
189  // bad idea!
190  assert(points.size() > 1);
191 
192  resizeArrays_(points.size());
193  typename XYContainer::const_iterator it = points.begin();
194  typename XYContainer::const_iterator endIt = points.end();
195  for (unsigned i = 0; it != endIt; ++i, ++it) {
196  xValues_[i] = std::get<0>(*it);
197  yValues_[i] = std::get<1>(*it);
198  }
199 
200  if (sortInputs)
201  sortInput_();
202  else if (xValues_[0] > xValues_[numSamples() - 1])
203  reverseSamplingPoints_();
204  }
205 
209  size_t numSamples() const
210  { return xValues_.size(); }
211 
215  Scalar xMin() const
216  { return xValues_[0]; }
217 
221  Scalar xMax() const
222  { return xValues_[numSamples() - 1]; }
223 
227  Scalar xAt(size_t i) const
228  { return xValues_[i]; }
229 
230  const std::vector<Scalar>& xValues() const
231  { return xValues_; }
232 
233  const std::vector<Scalar>& yValues() const
234  { return yValues_; }
235 
239  Scalar valueAt(size_t i) const
240  { return yValues_[i]; }
241 
245  template <class Evaluation>
246  bool applies(const Evaluation& x) const
247  { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
248 
258  template <class Evaluation>
259  Evaluation eval(const Evaluation& x, bool extrapolate = false) const
260  {
261  size_t segIdx = findSegmentIndex_(x, extrapolate);
262 
263  Scalar x0 = xValues_[segIdx];
264  Scalar x1 = xValues_[segIdx + 1];
265 
266  Scalar y0 = yValues_[segIdx];
267  Scalar y1 = yValues_[segIdx + 1];
268 
269  return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
270  }
271 
283  template <class Evaluation>
284  Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
285  {
286  unsigned segIdx = findSegmentIndex_(x, extrapolate);
287  return evalDerivative_(x, segIdx);
288  }
289 
304  template <class Evaluation>
305  Evaluation evalSecondDerivative(const Evaluation&, bool = false) const
306  { return 0.0; }
307 
322  template <class Evaluation>
323  Evaluation evalThirdDerivative(const Evaluation&, bool = false) const
324  { return 0.0; }
325 
334  int monotonic(Scalar x0, Scalar x1,
335  [[maybe_unused]] bool extrapolate = false) const
336  {
337  assert(x0 != x1);
338 
339  // make sure that x0 is smaller than x1
340  if (x0 > x1)
341  std::swap(x0, x1);
342 
343  assert(x0 < x1);
344 
345  int r = 3;
346  if (x0 < xMin()) {
347  assert(extrapolate);
348 
349  x0 = xMin();
350  };
351 
352  size_t i = findSegmentIndex_(x0, extrapolate);
353  if (xValues_[i + 1] >= x1) {
354  // interval is fully contained within a single function
355  // segment
356  updateMonotonicity_(i, r);
357  return r;
358  }
359 
360  // the first segment overlaps with the specified interval
361  // partially
362  updateMonotonicity_(i, r);
363  ++ i;
364 
365  // make sure that the segments which are completly in the
366  // interval [x0, x1] all exhibit the same monotonicity.
367  size_t iEnd = findSegmentIndex_(x1, extrapolate);
368  for (; i < iEnd - 1; ++i) {
369  updateMonotonicity_(i, r);
370  if (!r)
371  return 0;
372  }
373 
374  // if the user asked for a part of the function which is
375  // extrapolated, we need to check the slope at the function's
376  // endpoint
377  if (x1 > xMax()) {
378  assert(extrapolate);
379 
380  Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2);
381  if (m < 0)
382  return (r < 0 || r==3)?-1:0;
383  else if (m > 0)
384  return (r > 0 || r==3)?1:0;
385 
386  return r;
387  }
388 
389  // check for the last segment
390  updateMonotonicity_(iEnd, r);
391 
392  return r;
393  }
394 
399  int monotonic() const
400  { return monotonic(xMin(), xMax()); }
401 
418  void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os = std::cout) const
419  {
420  Scalar x0 = std::min(xi0, xi1);
421  Scalar x1 = std::max(xi0, xi1);
422  const int n = numSamples() - 1;
423  for (unsigned i = 0; i <= k; ++i) {
424  double x = i*(x1 - x0)/k + x0;
425  double y;
426  double dy_dx;
427  if (!applies(x)) {
428  if (x < xValues_[0]) {
429  dy_dx = evalDerivative(xValues_[0]);
430  y = (x - xValues_[0])*dy_dx + yValues_[0];
431  }
432  else if (x > xValues_[n]) {
433  dy_dx = evalDerivative(xValues_[n]);
434  y = (x - xValues_[n])*dy_dx + yValues_[n];
435  }
436  else {
437  throw std::runtime_error("The sampling points given to a function must be sorted by their x value!");
438  }
439  }
440  else {
441  y = eval(x);
442  dy_dx = evalDerivative(x);
443  }
444 
445  os << x << " " << y << " " << dy_dx << "\n";
446  }
447  }
448 
449  bool operator==(const Tabulated1DFunction<Scalar>& data) const {
450  return xValues_ == data.xValues_ &&
451  yValues_ == data.yValues_;
452  }
453 
454 private:
455  template <class Evaluation>
456  size_t findSegmentIndex_(const Evaluation& x, bool extrapolate = false) const
457  {
458  if (!isfinite(x)) {
459  std::ostringstream sstream;
460  sstream << "We can not search for extrapolation/interpolation segment in an 1D table for non-finite value " << getValue(x) << " .";
461  throw std::runtime_error(sstream.str());
462  }
463 
464  if (!extrapolate && !applies(x))
465  throw std::logic_error("Trying to evaluate a tabulated function outside of its range");
466 
467  // we need at least two sampling points!
468  if (numSamples() < 2) {
469  std::ostringstream sstream;
470  sstream << "We need at least two sampling points to do interpolation/extrapolation,"
471  "and the table only contains {} sampling points" << numSamples();
472  throw std::logic_error(sstream.str());
473  }
474 
475  if (x <= xValues_[1])
476  return 0;
477  else if (x >= xValues_[xValues_.size() - 2])
478  return xValues_.size() - 2;
479  else {
480  // bisection
481  size_t lowerIdx = 1;
482  size_t upperIdx = xValues_.size() - 2;
483  while (lowerIdx + 1 < upperIdx) {
484  size_t pivotIdx = (lowerIdx + upperIdx) / 2;
485  if (x < xValues_[pivotIdx])
486  upperIdx = pivotIdx;
487  else
488  lowerIdx = pivotIdx;
489  }
490 
491  if (xValues_[lowerIdx] > x || x > xValues_[lowerIdx + 1]) {
492  std::ostringstream sstream;
493  sstream << "Problematic interpolation/extrapolation segment is found for the input value " << Opm::getValue(x)
494  << "\nthe lower index of the found segment is " << lowerIdx << ", the size of the table is " << numSamples()
495  << ",\nand the end values of the found segment are " << xValues_[lowerIdx] << " and " << xValues_[lowerIdx + 1]
496  << ", respectively.";
497  std::ostringstream sstream2;
498  sstream2 << " Outputting the problematic table for more information(with *** marking the found segment):";
499  for (size_t i = 0; i < numSamples(); ++i) {
500  if (i % 10 == 0)
501  sstream2 << "\n";
502  if (i == lowerIdx)
503  sstream2 << " ***";
504  sstream2 << " " << xValues_[i];
505  if (i == lowerIdx + 1)
506  sstream2 << " ***";
507  }
508  sstream2<< "\n";
509  OpmLog::debug(sstream.str() + "\n" + sstream2.str());
510  throw std::runtime_error(sstream.str());
511  }
512  return lowerIdx;
513  }
514  }
515 
516  template <class Evaluation>
517  Evaluation evalDerivative_(const Evaluation& x, size_t segIdx) const
518  {
519  Scalar x0 = xValues_[segIdx];
520  Scalar x1 = xValues_[segIdx + 1];
521 
522  Scalar y0 = yValues_[segIdx];
523  Scalar y1 = yValues_[segIdx + 1];
524 
525  Evaluation ret = blank(x);
526  ret = (y1 - y0)/(x1 - x0);
527  return ret;
528  }
529 
530  // returns the monotonicity of a segment
531  //
532  // The return value have the following meaning:
533  //
534  // 3: function is constant within interval [x0, x1]
535  // 1: function is monotonously increasing in the specified interval
536  // 0: function is not monotonic in the specified interval
537  // -1: function is monotonously decreasing in the specified interval
538  int updateMonotonicity_(size_t i, int& r) const
539  {
540  if (yValues_[i] < yValues_[i + 1]) {
541  // monotonically increasing?
542  if (r == 3 || r == 1)
543  r = 1;
544  else
545  r = 0;
546  return 1;
547  }
548  else if (yValues_[i] > yValues_[i + 1]) {
549  // monotonically decreasing?
550  if (r == 3 || r == -1)
551  r = -1;
552  else
553  r = 0;
554  return -1;
555  }
556 
557  return 3;
558  }
559 
563  struct ComparatorX_
564  {
565  ComparatorX_(const std::vector<Scalar>& x)
566  : x_(x)
567  {}
568 
569  bool operator ()(size_t idxA, size_t idxB) const
570  { return x_.at(idxA) < x_.at(idxB); }
571 
572  const std::vector<Scalar>& x_;
573  };
574 
578  void sortInput_()
579  {
580  size_t n = numSamples();
581 
582  // create a vector containing 0...n-1
583  std::vector<unsigned> idxVector(n);
584  for (unsigned i = 0; i < n; ++i)
585  idxVector[i] = i;
586 
587  // sort the indices according to the x values of the sample
588  // points
589  ComparatorX_ cmp(xValues_);
590  std::sort(idxVector.begin(), idxVector.end(), cmp);
591 
592  // reorder the sample points
593  std::vector<Scalar> tmpX(n), tmpY(n);
594  for (size_t i = 0; i < idxVector.size(); ++ i) {
595  tmpX[i] = xValues_[idxVector[i]];
596  tmpY[i] = yValues_[idxVector[i]];
597  }
598  xValues_ = tmpX;
599  yValues_ = tmpY;
600  }
601 
606  void reverseSamplingPoints_()
607  {
608  // reverse the arrays
609  size_t n = numSamples();
610  for (size_t i = 0; i <= (n - 1)/2; ++i) {
611  std::swap(xValues_[i], xValues_[n - i - 1]);
612  std::swap(yValues_[i], yValues_[n - i - 1]);
613  }
614  }
615 
619  void resizeArrays_(size_t nSamples)
620  {
621  xValues_.resize(nSamples);
622  yValues_.resize(nSamples);
623  }
624 
625  std::vector<Scalar> xValues_;
626  std::vector<Scalar> yValues_;
627 };
628 } // namespace Opm
629 
630 #endif
Provides the opm-material specific exception classes.
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:259
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:209
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Tabulated1DFunction.hpp:284
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:150
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:221
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:81
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:215
Evaluation evalThirdDerivative(const Evaluation &, bool=false) const
Evaluate the function's third derivative at a given position.
Definition: Tabulated1DFunction.hpp:323
int monotonic(Scalar x0, Scalar x1, [[maybe_unused]] bool extrapolate=false) const
Returns 1 if the function is monotonically increasing, -1 if the function is mononously decreasing an...
Definition: Tabulated1DFunction.hpp:334
Evaluation evalSecondDerivative(const Evaluation &, bool=false) const
Evaluate the function's second derivative at a given position.
Definition: Tabulated1DFunction.hpp:305
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 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:418
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:185
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:227
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:93
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:55
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval.
Definition: Tabulated1DFunction.hpp:399
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
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:128
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:246
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:239