All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Spline.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_SPLINE_HPP
28 #define OPM_SPLINE_HPP
29 
32 #include <opm/common/ErrorMacros.hpp>
33 #include <opm/common/Exceptions.hpp>
34 #include <opm/common/Unused.hpp>
35 
36 #include <ostream>
37 #include <vector>
38 #include <tuple>
39 
40 namespace Opm
41 {
91 template<class Scalar>
92 class Spline
93 {
95  typedef std::vector<Scalar> Vector;
96 
97 public:
103  enum SplineType {
104  Full,
105  Natural,
106  Periodic,
107  Monotonic
108  };
109 
116  { }
117 
128  Spline(Scalar x0, Scalar x1,
129  Scalar y0, Scalar y1,
130  Scalar m0, Scalar m1)
131  { set(x0, x1, y0, y1, m0, m1); }
132 
141  template <class ScalarArrayX, class ScalarArrayY>
142  Spline(size_t nSamples,
143  const ScalarArrayX& x,
144  const ScalarArrayY& y,
145  SplineType splineType = Natural,
146  bool sortInputs = true)
147  { this->setXYArrays(nSamples, x, y, splineType, sortInputs); }
148 
156  template <class PointArray>
157  Spline(size_t nSamples,
158  const PointArray& points,
159  SplineType splineType = Natural,
160  bool sortInputs = true)
161  { this->setArrayOfPoints(nSamples, points, splineType, sortInputs); }
162 
170  template <class ScalarContainer>
171  Spline(const ScalarContainer& x,
172  const ScalarContainer& y,
173  SplineType splineType = Natural,
174  bool sortInputs = true)
175  { this->setXYContainers(x, y, splineType, sortInputs); }
176 
183  template <class PointContainer>
184  Spline(const PointContainer& points,
185  SplineType splineType = Natural,
186  bool sortInputs = true)
187  { this->setContainerOfPoints(points, splineType, sortInputs); }
188 
199  template <class ScalarArray>
200  Spline(size_t nSamples,
201  const ScalarArray& x,
202  const ScalarArray& y,
203  Scalar m0,
204  Scalar m1,
205  bool sortInputs = true)
206  { this->setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
207 
217  template <class PointArray>
218  Spline(size_t nSamples,
219  const PointArray& points,
220  Scalar m0,
221  Scalar m1,
222  bool sortInputs = true)
223  { this->setArrayOfPoints(nSamples, points, m0, m1, sortInputs); }
224 
234  template <class ScalarContainerX, class ScalarContainerY>
235  Spline(const ScalarContainerX& x,
236  const ScalarContainerY& y,
237  Scalar m0,
238  Scalar m1,
239  bool sortInputs = true)
240  { this->setXYContainers(x, y, m0, m1, sortInputs); }
241 
250  template <class PointContainer>
251  Spline(const PointContainer& points,
252  Scalar m0,
253  Scalar m1,
254  bool sortInputs = true)
255  { this->setContainerOfPoints(points, m0, m1, sortInputs); }
256 
268  void set(Scalar x0, Scalar x1,
269  Scalar y0, Scalar y1,
270  Scalar m0, Scalar m1)
271  {
272  slopeVec_.resize(2);
273  xPos_.resize(2);
274  yPos_.resize(2);
275 
276  if (x0 > x1) {
277  xPos_[0] = x1;
278  xPos_[1] = x0;
279  yPos_[0] = y1;
280  yPos_[1] = y0;
281  }
282  else {
283  xPos_[0] = x0;
284  xPos_[1] = x1;
285  yPos_[0] = y0;
286  yPos_[1] = y1;
287  }
288 
289  slopeVec_[0] = m0;
290  slopeVec_[1] = m1;
291 
292  Matrix M(numSamples());
293  Vector d(numSamples());
294  Vector moments(numSamples());
295  this->makeFullSystem_(M, d, m0, m1);
296 
297  // solve for the moments
298  M.solve(moments, d);
299 
300  this->setSlopesFromMoments_(slopeVec_, moments);
301  }
302 
303 
307  // Full splines //
311 
323  template <class ScalarArrayX, class ScalarArrayY>
324  void setXYArrays(size_t nSamples,
325  const ScalarArrayX& x,
326  const ScalarArrayY& y,
327  Scalar m0, Scalar m1,
328  bool sortInputs = true)
329  {
330  assert(nSamples > 1);
331 
332  setNumSamples_(nSamples);
333  for (size_t i = 0; i < nSamples; ++i) {
334  xPos_[i] = x[i];
335  yPos_[i] = y[i];
336  }
337 
338  if (sortInputs)
339  sortInput_();
340  else if (xPos_[0] > xPos_[numSamples() - 1])
342 
343  makeFullSpline_(m0, m1);
344  }
345 
357  template <class ScalarContainerX, class ScalarContainerY>
358  void setXYContainers(const ScalarContainerX& x,
359  const ScalarContainerY& y,
360  Scalar m0, Scalar m1,
361  bool sortInputs = true)
362  {
363  assert(x.size() == y.size());
364  assert(x.size() > 1);
365 
366  setNumSamples_(x.size());
367 
368  std::copy(x.begin(), x.end(), xPos_.begin());
369  std::copy(y.begin(), y.end(), yPos_.begin());
370 
371  if (sortInputs)
372  sortInput_();
373  else if (xPos_[0] > xPos_[numSamples() - 1])
375 
376  makeFullSpline_(m0, m1);
377  }
378 
391  template <class PointArray>
392  void setArrayOfPoints(size_t nSamples,
393  const PointArray& points,
394  Scalar m0,
395  Scalar m1,
396  bool sortInputs = true)
397  {
398  // a spline with no or just one sampling points? what an
399  // incredible bad idea!
400  assert(nSamples > 1);
401 
402  setNumSamples_(nSamples);
403  for (size_t i = 0; i < nSamples; ++i) {
404  xPos_[i] = points[i][0];
405  yPos_[i] = points[i][1];
406  }
407 
408  if (sortInputs)
409  sortInput_();
410  else if (xPos_[0] > xPos_[numSamples() - 1])
412 
413  makeFullSpline_(m0, m1);
414  }
415 
428  template <class XYContainer>
429  void setContainerOfPoints(const XYContainer& points,
430  Scalar m0,
431  Scalar m1,
432  bool sortInputs = true)
433  {
434  // a spline with no or just one sampling points? what an
435  // incredible bad idea!
436  assert(points.size() > 1);
437 
438  setNumSamples_(points.size());
439  typename XYContainer::const_iterator it = points.begin();
440  typename XYContainer::const_iterator endIt = points.end();
441  for (size_t i = 0; it != endIt; ++i, ++it) {
442  xPos_[i] = (*it)[0];
443  yPos_[i] = (*it)[1];
444  }
445 
446  if (sortInputs)
447  sortInput_();
448  else if (xPos_[0] > xPos_[numSamples() - 1])
450 
451  // make a full spline
452  makeFullSpline_(m0, m1);
453  }
454 
470  template <class XYContainer>
471  void setContainerOfTuples(const XYContainer& points,
472  Scalar m0,
473  Scalar m1,
474  bool sortInputs = true)
475  {
476  // resize internal arrays
477  setNumSamples_(points.size());
478  typename XYContainer::const_iterator it = points.begin();
479  typename XYContainer::const_iterator endIt = points.end();
480  for (unsigned i = 0; it != endIt; ++i, ++it) {
481  xPos_[i] = std::get<0>(*it);
482  yPos_[i] = std::get<1>(*it);
483  }
484 
485  if (sortInputs)
486  sortInput_();
487  else if (xPos_[0] > xPos_[numSamples() - 1])
489 
490  // make a full spline
491  makeFullSpline_(m0, m1);
492  }
493 
497  // Natural/Periodic splines //
501 
511  template <class ScalarArrayX, class ScalarArrayY>
512  void setXYArrays(size_t nSamples,
513  const ScalarArrayX& x,
514  const ScalarArrayY& y,
515  SplineType splineType = Natural,
516  bool sortInputs = true)
517  {
518  assert(nSamples > 1);
519 
520  setNumSamples_(nSamples);
521  for (size_t i = 0; i < nSamples; ++i) {
522  xPos_[i] = x[i];
523  yPos_[i] = y[i];
524  }
525 
526  if (sortInputs)
527  sortInput_();
528  else if (xPos_[0] > xPos_[numSamples() - 1])
530 
531  if (splineType == Periodic)
533  else if (splineType == Natural)
535  else if (splineType == Monotonic)
536  this->makeMonotonicSpline_(slopeVec_);
537  else
538  OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
539  }
540 
552  template <class ScalarContainerX, class ScalarContainerY>
553  void setXYContainers(const ScalarContainerX& x,
554  const ScalarContainerY& y,
555  SplineType splineType = Natural,
556  bool sortInputs = true)
557  {
558  assert(x.size() == y.size());
559  assert(x.size() > 1);
560 
561  setNumSamples_(x.size());
562  std::copy(x.begin(), x.end(), xPos_.begin());
563  std::copy(y.begin(), y.end(), yPos_.begin());
564 
565  if (sortInputs)
566  sortInput_();
567  else if (xPos_[0] > xPos_[numSamples() - 1])
569 
570  if (splineType == Periodic)
572  else if (splineType == Natural)
574  else if (splineType == Monotonic)
575  this->makeMonotonicSpline_(slopeVec_);
576  else
577  OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
578  }
579 
592  template <class PointArray>
593  void setArrayOfPoints(size_t nSamples,
594  const PointArray& points,
595  SplineType splineType = Natural,
596  bool sortInputs = true)
597  {
598  // a spline with no or just one sampling points? what an
599  // incredible bad idea!
600  assert(nSamples > 1);
601 
602  setNumSamples_(nSamples);
603  for (size_t i = 0; i < nSamples; ++i) {
604  xPos_[i] = points[i][0];
605  yPos_[i] = points[i][1];
606  }
607 
608  if (sortInputs)
609  sortInput_();
610  else if (xPos_[0] > xPos_[numSamples() - 1])
612 
613  if (splineType == Periodic)
615  else if (splineType == Natural)
617  else if (splineType == Monotonic)
618  this->makeMonotonicSpline_(slopeVec_);
619  else
620  OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
621  }
622 
634  template <class XYContainer>
635  void setContainerOfPoints(const XYContainer& points,
636  SplineType splineType = Natural,
637  bool sortInputs = true)
638  {
639  // a spline with no or just one sampling points? what an
640  // incredible bad idea!
641  assert(points.size() > 1);
642 
643  setNumSamples_(points.size());
644  typename XYContainer::const_iterator it = points.begin();
645  typename XYContainer::const_iterator endIt = points.end();
646  for (size_t i = 0; it != endIt; ++ i, ++it) {
647  xPos_[i] = (*it)[0];
648  yPos_[i] = (*it)[1];
649  }
650 
651  if (sortInputs)
652  sortInput_();
653  else if (xPos_[0] > xPos_[numSamples() - 1])
655 
656  if (splineType == Periodic)
658  else if (splineType == Natural)
660  else if (splineType == Monotonic)
661  this->makeMonotonicSpline_(slopeVec_);
662  else
663  OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
664  }
665 
680  template <class XYContainer>
681  void setContainerOfTuples(const XYContainer& points,
682  SplineType splineType = Natural,
683  bool sortInputs = true)
684  {
685  // resize internal arrays
686  setNumSamples_(points.size());
687  typename XYContainer::const_iterator it = points.begin();
688  typename XYContainer::const_iterator endIt = points.end();
689  for (unsigned i = 0; it != endIt; ++i, ++it) {
690  xPos_[i] = std::get<0>(*it);
691  yPos_[i] = std::get<1>(*it);
692  }
693 
694  if (sortInputs)
695  sortInput_();
696  else if (xPos_[0] > xPos_[numSamples() - 1])
698 
699  if (splineType == Periodic)
701  else if (splineType == Natural)
703  else if (splineType == Monotonic)
704  this->makeMonotonicSpline_(slopeVec_);
705  else
706  OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
707  }
708 
712  template <class Evaluation>
713  bool applies(const Evaluation& x) const
714  { return x_(0) <= x && x <= x_(numSamples() - 1); }
715 
719  size_t numSamples() const
720  { return xPos_.size(); }
721 
725  Scalar xAt(size_t sampleIdx) const
726  { return x_(sampleIdx); }
727 
731  Scalar valueAt(size_t sampleIdx) const
732  { return y_(sampleIdx); }
733 
751  void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream& os = std::cout) const
752  {
753  Scalar x0 = std::min(xi0, xi1);
754  Scalar x1 = std::max(xi0, xi1);
755  const size_t n = numSamples() - 1;
756  for (size_t i = 0; i <= k; ++i) {
757  double x = i*(x1 - x0)/k + x0;
758  double x_p1 = x + (x1 - x0)/k;
759  double y;
760  double dy_dx;
761  double mono = 1;
762  if (!applies(x)) {
763  if (x < x_(0)) {
764  dy_dx = evalDerivative(x_(0));
765  y = (x - x_(0))*dy_dx + y_(0);
766  mono = (dy_dx>0)?1:-1;
767  }
768  else if (x > x_(n)) {
769  dy_dx = evalDerivative(x_(n));
770  y = (x - x_(n))*dy_dx + y_(n);
771  mono = (dy_dx>0)?1:-1;
772  }
773  else {
774  OPM_THROW(std::runtime_error,
775  "The sampling points given to a spline must be sorted by their x value!");
776  }
777  }
778  else {
779  y = eval(x);
780  dy_dx = evalDerivative(x);
781  mono = monotonic(x, x_p1, /*extrapolate=*/true);
782  }
783 
784  os << x << " " << y << " " << dy_dx << " " << mono << "\n";
785  }
786  }
787 
799  template <class Evaluation>
800  Evaluation eval(const Evaluation& x, bool extrapolate = false) const
801  {
802  if (!extrapolate && !applies(x))
803  OPM_THROW(Opm::NumericalProblem,
804  "Tried to evaluate a spline outside of its range");
805 
806  // handle extrapolation
807  if (extrapolate) {
808  if (x < xAt(0)) {
809  Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
810  Scalar y0 = y_(0);
811  return y0 + m*(x - xAt(0));
812  }
813  else if (x > xAt(static_cast<size_t>(static_cast<long int>(numSamples()) - 1))) {
814  Scalar m = evalDerivative_(xAt(static_cast<size_t>(numSamples() - 1)),
815  /*segmentIdx=*/static_cast<size_t>(numSamples()-2));
816  Scalar y0 = y_(static_cast<size_t>(numSamples() - 1));
817  return y0 + m*(x - xAt(static_cast<size_t>(numSamples() - 1)));
818  }
819  }
820 
821  return eval_(x, segmentIdx_(Opm::scalarValue(x)));
822  }
823 
836  template <class Evaluation>
837  Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
838  {
839  if (!extrapolate && !applies(x))
840  OPM_THROW(Opm::NumericalProblem,
841  "Tried to evaluate the derivative of a spline outside of its range");
842 
843  // handle extrapolation
844  if (extrapolate) {
845  if (x < xAt(0))
846  return evalDerivative_(xAt(0), /*segmentIdx=*/0);
847  else if (x > xAt(numSamples() - 1))
848  return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
849  }
850 
851  return evalDerivative_(x, segmentIdx_(Opm::scalarValue(x)));
852  }
853 
866  template <class Evaluation>
867  Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
868  {
869  if (!extrapolate && !applies(x))
870  OPM_THROW(Opm::NumericalProblem,
871  "Tried to evaluate the second derivative of a spline outside of its range");
872  else if (extrapolate)
873  return 0.0;
874 
875  return evalDerivative2_(x, segmentIdx_(Opm::scalarValue(x)));
876  }
877 
890  template <class Evaluation>
891  Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
892  {
893  if (!extrapolate && !applies(x))
894  OPM_THROW(Opm::NumericalProblem,
895  "Tried to evaluate the third derivative of a spline outside of its range");
896  else if (extrapolate)
897  return 0.0;
898 
899  return evalDerivative3_(x, segmentIdx_(Opm::scalarValue(x)));
900  }
901 
908  template <class Evaluation>
909  Evaluation intersect(const Evaluation& a,
910  const Evaluation& b,
911  const Evaluation& c,
912  const Evaluation& d) const
913  {
914  return intersectInterval(xAt(0), xAt(numSamples() - 1), a, b, c, d);
915  }
916 
923  template <class Evaluation>
924  Evaluation intersectInterval(Scalar x0, Scalar x1,
925  const Evaluation& a,
926  const Evaluation& b,
927  const Evaluation& c,
928  const Evaluation& d) const
929  {
930  assert(applies(x0) && applies(x1));
931 
932  Evaluation tmpSol[3], sol = 0;
933  size_t nSol = 0;
934  size_t iFirst = segmentIdx_(x0);
935  size_t iLast = segmentIdx_(x1);
936  for (size_t i = iFirst; i <= iLast; ++i)
937  {
938  size_t nCur = intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
939  if (nCur == 1)
940  sol = tmpSol[0];
941 
942  nSol += nCur;
943  if (nSol > 1) {
944  OPM_THROW(std::runtime_error,
945  "Spline has more than one intersection"); //<<a<<"x^3 + "<<b<"x^2 + "<<c<"x + "<<d);
946  }
947  }
948 
949  if (nSol != 1)
950  OPM_THROW(std::runtime_error,
951  "Spline has no intersection"); //<<a<"x^3 + " <<b<"x^2 + "<<c<"x + "<<d<<"!");
952 
953  return sol;
954  }
955 
964  int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED = false) const
965  {
966  assert(std::abs(x0 - x1) > 1e-30);
967 
968  // make sure that x0 is smaller than x1
969  if (x0 > x1)
970  std::swap(x0, x1);
971 
972  assert(x0 < x1);
973 
974  int r = 3;
975  if (x0 < xAt(0)) {
976  assert(extrapolate);
977  Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
978  if (std::abs(m) < 1e-20)
979  r = (m < 0)?-1:1;
980  x0 = xAt(0);
981  };
982 
983  size_t i = segmentIdx_(x0);
984  if (x_(i + 1) >= x1) {
985  // interval is fully contained within a single spline
986  // segment
987  monotonic_(i, x0, x1, r);
988  return r;
989  }
990 
991  // the first segment overlaps with the specified interval
992  // partially
993  monotonic_(i, x0, x_(i+1), r);
994  ++ i;
995 
996  // make sure that the segments which are completly in the
997  // interval [x0, x1] all exhibit the same monotonicity.
998  size_t iEnd = segmentIdx_(x1);
999  for (; i < iEnd - 1; ++i) {
1000  monotonic_(i, x_(i), x_(i + 1), r);
1001  if (!r)
1002  return 0;
1003  }
1004 
1005  // if the user asked for a part of the spline which is
1006  // extrapolated, we need to check the slope at the spline's
1007  // endpoint
1008  if (x1 > xAt(numSamples() - 1)) {
1009  assert(extrapolate);
1010 
1011  Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
1012  if (m < 0)
1013  return (r < 0 || r==3)?-1:0;
1014  else if (m > 0)
1015  return (r > 0 || r==3)?1:0;
1016 
1017  return r;
1018  }
1019 
1020  // check for the last segment
1021  monotonic_(iEnd, x_(iEnd), x1, r);
1022 
1023  return r;
1024  }
1025 
1030  int monotonic() const
1031  { return monotonic(xAt(0), xAt(numSamples() - 1)); }
1032 
1033 protected:
1038  {
1039  ComparatorX_(const std::vector<Scalar>& x)
1040  : x_(x)
1041  {}
1042 
1043  bool operator ()(unsigned idxA, unsigned idxB) const
1044  { return x_.at(idxA) < x_.at(idxB); }
1045 
1046  const std::vector<Scalar>& x_;
1047  };
1048 
1052  void sortInput_()
1053  {
1054  size_t n = numSamples();
1055 
1056  // create a vector containing 0...n-1
1057  std::vector<unsigned> idxVector(n);
1058  for (unsigned i = 0; i < n; ++i)
1059  idxVector[i] = i;
1060 
1061  // sort the indices according to the x values of the sample
1062  // points
1063  ComparatorX_ cmp(xPos_);
1064  std::sort(idxVector.begin(), idxVector.end(), cmp);
1065 
1066  // reorder the sample points
1067  std::vector<Scalar> tmpX(n), tmpY(n);
1068  for (size_t i = 0; i < idxVector.size(); ++ i) {
1069  tmpX[i] = xPos_[idxVector[i]];
1070  tmpY[i] = yPos_[idxVector[i]];
1071  }
1072  xPos_ = tmpX;
1073  yPos_ = tmpY;
1074  }
1075 
1081  {
1082  // reverse the arrays
1083  size_t n = numSamples();
1084  for (unsigned i = 0; i <= (n - 1)/2; ++i) {
1085  std::swap(xPos_[i], xPos_[n - i - 1]);
1086  std::swap(yPos_[i], yPos_[n - i - 1]);
1087  }
1088  }
1089 
1093  void setNumSamples_(size_t nSamples)
1094  {
1095  xPos_.resize(nSamples);
1096  yPos_.resize(nSamples);
1097  slopeVec_.resize(nSamples);
1098  }
1099 
1105  void makeFullSpline_(Scalar m0, Scalar m1)
1106  {
1107  Matrix M(numSamples());
1108  std::vector<Scalar> d(numSamples());
1109  std::vector<Scalar> moments(numSamples());
1110 
1111  // create linear system of equations
1112  this->makeFullSystem_(M, d, m0, m1);
1113 
1114  // solve for the moments (-> second derivatives)
1115  M.solve(moments, d);
1116 
1117  // convert the moments to slopes at the sample points
1118  this->setSlopesFromMoments_(slopeVec_, moments);
1119  }
1120 
1127  {
1128  Matrix M(numSamples(), numSamples());
1129  Vector d(numSamples());
1130  Vector moments(numSamples());
1131 
1132  // create linear system of equations
1133  this->makeNaturalSystem_(M, d);
1134 
1135  // solve for the moments (-> second derivatives)
1136  M.solve(moments, d);
1137 
1138  // convert the moments to slopes at the sample points
1139  this->setSlopesFromMoments_(slopeVec_, moments);
1140  }
1141 
1148  {
1149  Matrix M(numSamples() - 1);
1150  Vector d(numSamples() - 1);
1151  Vector moments(numSamples() - 1);
1152 
1153  // create linear system of equations. This is a bit hacky,
1154  // because it assumes that std::vector internally stores its
1155  // data as a big C-style array, but it saves us from yet
1156  // another copy operation
1157  this->makePeriodicSystem_(M, d);
1158 
1159  // solve for the moments (-> second derivatives)
1160  M.solve(moments, d);
1161 
1162  moments.resize(numSamples());
1163  for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i) {
1164  unsigned ui = static_cast<unsigned>(i);
1165  moments[ui+1] = moments[ui];
1166  }
1167  moments[0] = moments[numSamples() - 1];
1168 
1169  // convert the moments to slopes at the sample points
1170  this->setSlopesFromMoments_(slopeVec_, moments);
1171  }
1172 
1179  template <class DestVector, class SourceVector>
1180  void assignSamplingPoints_(DestVector& destX,
1181  DestVector& destY,
1182  const SourceVector& srcX,
1183  const SourceVector& srcY,
1184  unsigned nSamples)
1185  {
1186  assert(nSamples >= 2);
1187 
1188  // copy sample points, make sure that the first x value is
1189  // smaller than the last one
1190  for (unsigned i = 0; i < nSamples; ++i) {
1191  unsigned idx = i;
1192  if (srcX[0] > srcX[nSamples - 1])
1193  idx = nSamples - i - 1;
1194  destX[i] = srcX[idx];
1195  destY[i] = srcY[idx];
1196  }
1197  }
1198 
1199  template <class DestVector, class ListIterator>
1200  void assignFromArrayList_(DestVector& destX,
1201  DestVector& destY,
1202  const ListIterator& srcBegin,
1203  const ListIterator& srcEnd,
1204  unsigned nSamples)
1205  {
1206  assert(nSamples >= 2);
1207 
1208  // find out wether the x values are in reverse order
1209  ListIterator it = srcBegin;
1210  ++it;
1211  bool reverse = false;
1212  if ((*srcBegin)[0] > (*it)[0])
1213  reverse = true;
1214  --it;
1215 
1216  // loop over all sampling points
1217  for (unsigned i = 0; it != srcEnd; ++i, ++it) {
1218  unsigned idx = i;
1219  if (reverse)
1220  idx = nSamples - i - 1;
1221  destX[i] = (*it)[0];
1222  destY[i] = (*it)[1];
1223  }
1224  }
1225 
1233  template <class DestVector, class ListIterator>
1234  void assignFromTupleList_(DestVector& destX,
1235  DestVector& destY,
1236  ListIterator srcBegin,
1237  ListIterator srcEnd,
1238  unsigned nSamples)
1239  {
1240  assert(nSamples >= 2);
1241 
1242  // copy sample points, make sure that the first x value is
1243  // smaller than the last one
1244 
1245  // find out wether the x values are in reverse order
1246  ListIterator it = srcBegin;
1247  ++it;
1248  bool reverse = false;
1249  if (std::get<0>(*srcBegin) > std::get<0>(*it))
1250  reverse = true;
1251  --it;
1252 
1253  // loop over all sampling points
1254  for (unsigned i = 0; it != srcEnd; ++i, ++it) {
1255  unsigned idx = i;
1256  if (reverse)
1257  idx = nSamples - i - 1;
1258  destX[i] = std::get<0>(*it);
1259  destY[i] = std::get<1>(*it);
1260  }
1261  }
1262 
1267  template <class Vector, class Matrix>
1268  void makeFullSystem_(Matrix& M, Vector& d, Scalar m0, Scalar m1)
1269  {
1270  makeNaturalSystem_(M, d);
1271 
1272  size_t n = numSamples() - 1;
1273  // first row
1274  M[0][1] = 1;
1275  d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0);
1276 
1277  // last row
1278  M[n][n - 1] = 1;
1279 
1280  // right hand side
1281  d[n] =
1282  6/h_(n)
1283  *
1284  (m1 - (y_(n) - y_(n - 1))/h_(n));
1285  }
1286 
1291  template <class Vector, class Matrix>
1292  void makeNaturalSystem_(Matrix& M, Vector& d)
1293  {
1294  M = 0.0;
1295 
1296  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1297  // Springer, 2005, p. 111
1298  size_t n = numSamples() - 1;
1299 
1300  // second to next to last rows
1301  for (size_t i = 1; i < n; ++i) {
1302  Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
1303  Scalar mu_i = 1 - lambda_i;
1304  Scalar d_i =
1305  6 / (h_(i) + h_(i + 1))
1306  *
1307  ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
1308 
1309  M[i][i-1] = mu_i;
1310  M[i][i] = 2;
1311  M[i][i + 1] = lambda_i;
1312  d[i] = d_i;
1313  };
1314 
1315  // See Stroer, equation (2.5.2.7)
1316  Scalar lambda_0 = 0;
1317  Scalar d_0 = 0;
1318 
1319  Scalar mu_n = 0;
1320  Scalar d_n = 0;
1321 
1322  // first row
1323  M[0][0] = 2;
1324  M[0][1] = lambda_0;
1325  d[0] = d_0;
1326 
1327  // last row
1328  M[n][n-1] = mu_n;
1329  M[n][n] = 2;
1330  d[n] = d_n;
1331  }
1332 
1337  template <class Matrix, class Vector>
1338  void makePeriodicSystem_(Matrix& M, Vector& d)
1339  {
1340  M = 0.0;
1341 
1342  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1343  // Springer, 2005, p. 111
1344  size_t n = numSamples() - 1;
1345 
1346  assert(M.rows() == n);
1347 
1348  // second to next to last rows
1349  for (size_t i = 2; i < n; ++i) {
1350  Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
1351  Scalar mu_i = 1 - lambda_i;
1352  Scalar d_i =
1353  6 / (h_(i) + h_(i + 1))
1354  *
1355  ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
1356 
1357  M[i-1][i-2] = mu_i;
1358  M[i-1][i-1] = 2;
1359  M[i-1][i] = lambda_i;
1360  d[i-1] = d_i;
1361  };
1362 
1363  Scalar lambda_n = h_(1) / (h_(n) + h_(1));
1364  Scalar lambda_1 = h_(2) / (h_(1) + h_(2));;
1365  Scalar mu_1 = 1 - lambda_1;
1366  Scalar mu_n = 1 - lambda_n;
1367 
1368  Scalar d_1 =
1369  6 / (h_(1) + h_(2))
1370  *
1371  ( (y_(2) - y_(1))/h_(2) - (y_(1) - y_(0))/h_(1));
1372  Scalar d_n =
1373  6 / (h_(n) + h_(1))
1374  *
1375  ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n));
1376 
1377 
1378  // first row
1379  M[0][0] = 2;
1380  M[0][1] = lambda_1;
1381  M[0][n-1] = mu_1;
1382  d[0] = d_1;
1383 
1384  // last row
1385  M[n-1][0] = lambda_n;
1386  M[n-1][n-2] = mu_n;
1387  M[n-1][n-1] = 2;
1388  d[n-1] = d_n;
1389  }
1390 
1399  template <class Vector>
1400  void makeMonotonicSpline_(Vector& slopes)
1401  {
1402  auto n = numSamples();
1403 
1404  // calculate the slopes of the secant lines
1405  std::vector<Scalar> delta(n);
1406  for (size_t k = 0; k < n - 1; ++k)
1407  delta[k] = (y_(k + 1) - y_(k))/(x_(k + 1) - x_(k));
1408 
1409  // calculate the "raw" slopes at the sample points
1410  for (size_t k = 1; k < n - 1; ++k)
1411  slopes[k] = (delta[k - 1] + delta[k])/2;
1412  slopes[0] = delta[0];
1413  slopes[n - 1] = delta[n - 2];
1414 
1415  // post-process the "raw" slopes at the sample points
1416  for (size_t k = 0; k < n - 1; ++k) {
1417  if (std::abs(delta[k]) < 1e-50) {
1418  // make the spline flat if the inputs are equal
1419  slopes[k] = 0;
1420  slopes[k + 1] = 0;
1421  ++ k;
1422  continue;
1423  }
1424  else {
1425  Scalar alpha = slopes[k] / delta[k];
1426  Scalar beta = slopes[k + 1] / delta[k];
1427 
1428  if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
1429  slopes[k] = 0;
1430  }
1431  // limit (alpha, beta) to a circle of radius 3
1432  else if (alpha*alpha + beta*beta > 3*3) {
1433  Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
1434  slopes[k] = tau*alpha*delta[k];
1435  slopes[k + 1] = tau*beta*delta[k];
1436  }
1437  }
1438  }
1439  }
1440 
1448  template <class MomentsVector, class SlopeVector>
1449  void setSlopesFromMoments_(SlopeVector& slopes, const MomentsVector& moments)
1450  {
1451  size_t n = numSamples();
1452 
1453  // evaluate slope at the rightmost point.
1454  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1455  // Springer, 2005, p. 109
1456  Scalar mRight;
1457 
1458  {
1459  Scalar h = this->h_(n - 1);
1460  Scalar x = h;
1461  //Scalar x_1 = 0;
1462 
1463  Scalar A =
1464  (y_(n - 1) - y_(n - 2))/h
1465  -
1466  h/6*(moments[n-1] - moments[n - 2]);
1467 
1468  mRight =
1469  //- moments[n - 2] * x_1*x_1 / (2 * h)
1470  //+
1471  moments[n - 1] * x*x / (2 * h)
1472  +
1473  A;
1474  }
1475 
1476  // evaluate the slope for the first n-1 sample points
1477  for (size_t i = 0; i < n - 1; ++ i) {
1478  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1479  // Springer, 2005, p. 109
1480  Scalar h_i = this->h_(i + 1);
1481  //Scalar x_i = 0;
1482  Scalar x_i1 = h_i;
1483 
1484  Scalar A_i =
1485  (y_(i+1) - y_(i))/h_i
1486  -
1487  h_i/6*(moments[i+1] - moments[i]);
1488 
1489  slopes[i] =
1490  - moments[i] * x_i1*x_i1 / (2 * h_i)
1491  +
1492  //moments[i + 1] * x_i*x_i / (2 * h_i)
1493  //+
1494  A_i;
1495 
1496  }
1497  slopes[n - 1] = mRight;
1498  }
1499 
1500 
1501  // evaluate the spline at a given the position and given the
1502  // segment index
1503  template <class Evaluation>
1504  Evaluation eval_(const Evaluation& x, size_t i) const
1505  {
1506  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1507  Scalar delta = h_(i + 1);
1508  Evaluation t = (x - x_(i))/delta;
1509 
1510  return
1511  h00_(t) * y_(i)
1512  + h10_(t) * slope_(i)*delta
1513  + h01_(t) * y_(i + 1)
1514  + h11_(t) * slope_(i + 1)*delta;
1515  }
1516 
1517  // evaluate the derivative of a spline given the actual position
1518  // and the segment index
1519  template <class Evaluation>
1520  Evaluation evalDerivative_(const Evaluation& x, size_t i) const
1521  {
1522  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1523  Scalar delta = h_(i + 1);
1524  Evaluation t = (x - x_(i))/delta;
1525  Evaluation alpha = 1 / delta;
1526 
1527  return
1528  alpha *
1529  (h00_prime_(t) * y_(i)
1530  + h10_prime_(t) * slope_(i)*delta
1531  + h01_prime_(t) * y_(i + 1)
1532  + h11_prime_(t) * slope_(i + 1)*delta);
1533  }
1534 
1535  // evaluate the second derivative of a spline given the actual
1536  // position and the segment index
1537  template <class Evaluation>
1538  Evaluation evalDerivative2_(const Evaluation& x, size_t i) const
1539  {
1540  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1541  Scalar delta = h_(i + 1);
1542  Evaluation t = (x - x_(i))/delta;
1543  Evaluation alpha = 1 / delta;
1544 
1545  return
1546  alpha*alpha
1547  *(h00_prime2_(t) * y_(i)
1548  + h10_prime2_(t) * slope_(i)*delta
1549  + h01_prime2_(t) * y_(i + 1)
1550  + h11_prime2_(t) * slope_(i + 1)*delta);
1551  }
1552 
1553  // evaluate the third derivative of a spline given the actual
1554  // position and the segment index
1555  template <class Evaluation>
1556  Evaluation evalDerivative3_(const Evaluation& x, size_t i) const
1557  {
1558  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1559  Scalar delta = h_(i + 1);
1560  Evaluation t = (x - x_(i))/delta;
1561  Evaluation alpha = 1 / delta;
1562 
1563  return
1564  alpha*alpha*alpha
1565  *(h00_prime3_(t)*y_(i)
1566  + h10_prime3_(t)*slope_(i)*delta
1567  + h01_prime3_(t)*y_(i + 1)
1568  + h11_prime3_(t)*slope_(i + 1)*delta);
1569  }
1570 
1571  // hermite basis functions
1572  template <class Evaluation>
1573  Evaluation h00_(const Evaluation& t) const
1574  { return (2*t - 3)*t*t + 1; }
1575 
1576  template <class Evaluation>
1577  Evaluation h10_(const Evaluation& t) const
1578  { return ((t - 2)*t + 1)*t; }
1579 
1580  template <class Evaluation>
1581  Evaluation h01_(const Evaluation& t) const
1582  { return (-2*t + 3)*t*t; }
1583 
1584  template <class Evaluation>
1585  Evaluation h11_(const Evaluation& t) const
1586  { return (t - 1)*t*t; }
1587 
1588  // first derivative of the hermite basis functions
1589  template <class Evaluation>
1590  Evaluation h00_prime_(const Evaluation& t) const
1591  { return (3*2*t - 2*3)*t; }
1592 
1593  template <class Evaluation>
1594  Evaluation h10_prime_(const Evaluation& t) const
1595  { return (3*t - 2*2)*t + 1; }
1596 
1597  template <class Evaluation>
1598  Evaluation h01_prime_(const Evaluation& t) const
1599  { return (-3*2*t + 2*3)*t; }
1600 
1601  template <class Evaluation>
1602  Evaluation h11_prime_(const Evaluation& t) const
1603  { return (3*t - 2)*t; }
1604 
1605  // second derivative of the hermite basis functions
1606  template <class Evaluation>
1607  Evaluation h00_prime2_(const Evaluation& t) const
1608  { return 2*3*2*t - 2*3; }
1609 
1610  template <class Evaluation>
1611  Evaluation h10_prime2_(const Evaluation& t) const
1612  { return 2*3*t - 2*2; }
1613 
1614  template <class Evaluation>
1615  Evaluation h01_prime2_(const Evaluation& t) const
1616  { return -2*3*2*t + 2*3; }
1617 
1618  template <class Evaluation>
1619  Evaluation h11_prime2_(const Evaluation& t) const
1620  { return 2*3*t - 2; }
1621 
1622  // third derivative of the hermite basis functions
1623  template <class Evaluation>
1624  Scalar h00_prime3_(const Evaluation& t OPM_UNUSED) const
1625  { return 2*3*2; }
1626 
1627  template <class Evaluation>
1628  Scalar h10_prime3_(const Evaluation& t OPM_UNUSED) const
1629  { return 2*3; }
1630 
1631  template <class Evaluation>
1632  Scalar h01_prime3_(const Evaluation& t OPM_UNUSED) const
1633  { return -2*3*2; }
1634 
1635  template <class Evaluation>
1636  Scalar h11_prime3_(const Evaluation& t OPM_UNUSED) const
1637  { return 2*3; }
1638 
1639  // returns the monotonicality of an interval of a spline segment
1640  //
1641  // The return value have the following meaning:
1642  //
1643  // 3: spline is constant within interval [x0, x1]
1644  // 1: spline is monotonously increasing in the specified interval
1645  // 0: spline is not monotonic (or constant) in the specified interval
1646  // -1: spline is monotonously decreasing in the specified interval
1647  int monotonic_(size_t i, Scalar x0, Scalar x1, int& r) const
1648  {
1649  // coefficients of derivative in monomial basis
1650  Scalar a = 3*a_(i);
1651  Scalar b = 2*b_(i);
1652  Scalar c = c_(i);
1653 
1654  if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
1655  return 3; // constant in interval, r stays unchanged!
1656 
1657  Scalar disc = b*b - 4*a*c;
1658  if (disc < 0) {
1659  // discriminant of derivative is smaller than 0, i.e. the
1660  // segment's derivative does not exhibit any extrema.
1661  if (x0*(x0*a + b) + c > 0) {
1662  r = (r==3 || r == 1)?1:0;
1663  return 1;
1664  }
1665  else {
1666  r = (r==3 || r == -1)?-1:0;
1667  return -1;
1668  }
1669  }
1670  disc = std::sqrt(disc);
1671  Scalar xE1 = (-b + disc)/(2*a);
1672  Scalar xE2 = (-b - disc)/(2*a);
1673 
1674  if (std::abs(disc) < 1e-30) {
1675  // saddle point -> no extrema
1676  if (std::abs(xE1 - x0) < 1e-30)
1677  // make sure that we're not picking the saddle point
1678  // to determine whether we're monotonically increasing
1679  // or decreasing
1680  x0 = x1;
1681  if (x0*(x0*a + b) + c > 0) {
1682  r = (r==3 || r == 1)?1:0;
1683  return 1;
1684  }
1685  else {
1686  r = (r==3 || r == -1)?-1:0;
1687  return -1;
1688  }
1689  };
1690  if ((x0 < xE1 && xE1 < x1) ||
1691  (x0 < xE2 && xE2 < x1))
1692  {
1693  // there is an extremum in the range (x0, x1)
1694  r = 0;
1695  return 0;
1696  }
1697  // no extremum in range (x0, x1)
1698  x0 = (x0 + x1)/2; // pick point in the middle of the interval
1699  // to avoid extrema on the boundaries
1700  if (x0*(x0*a + b) + c > 0) {
1701  r = (r==3 || r == 1)?1:0;
1702  return 1;
1703  }
1704  else {
1705  r = (r==3 || r == -1)?-1:0;
1706  return -1;
1707  }
1708  }
1709 
1714  template <class Evaluation>
1715  size_t intersectSegment_(Evaluation* sol,
1716  size_t segIdx,
1717  const Evaluation& a,
1718  const Evaluation& b,
1719  const Evaluation& c,
1720  const Evaluation& d,
1721  Scalar x0 = -1e30, Scalar x1 = 1e30) const
1722  {
1723  unsigned n =
1724  Opm::invertCubicPolynomial(sol,
1725  a_(segIdx) - a,
1726  b_(segIdx) - b,
1727  c_(segIdx) - c,
1728  d_(segIdx) - d);
1729  x0 = std::max(x_(segIdx), x0);
1730  x1 = std::min(x_(segIdx+1), x1);
1731 
1732  // filter the intersections outside of the specified interval
1733  size_t k = 0;
1734  for (unsigned j = 0; j < n; ++j) {
1735  if (x0 <= sol[j] && sol[j] <= x1) {
1736  sol[k] = sol[j];
1737  ++k;
1738  }
1739  }
1740  return k;
1741  }
1742 
1743  // find the segment index for a given x coordinate
1744  size_t segmentIdx_(Scalar x) const
1745  {
1746  // bisection
1747  size_t iLow = 0;
1748  size_t iHigh = numSamples() - 1;
1749 
1750  while (iLow + 1 < iHigh) {
1751  size_t i = (iLow + iHigh) / 2;
1752  if (x_(i) > x)
1753  iHigh = i;
1754  else
1755  iLow = i;
1756  };
1757  return iLow;
1758  }
1759 
1763  Scalar h_(size_t i) const
1764  {
1765  assert(x_(i) > x_(i-1)); // the sampling points must be given
1766  // in ascending order
1767  return x_(i) - x_(i - 1);
1768  }
1769 
1773  Scalar x_(size_t i) const
1774  { return xPos_[i]; }
1775 
1779  Scalar y_(size_t i) const
1780  { return yPos_[i]; }
1781 
1786  Scalar slope_(size_t i) const
1787  { return slopeVec_[i]; }
1788 
1789  // returns the coefficient in front of the x^3 term. In Stoer this
1790  // is delta.
1791  Scalar a_(size_t i) const
1792  { return evalDerivative3_(/*x=*/Scalar(0.0), i)/6.0; }
1793 
1794  // returns the coefficient in front of the x^2 term In Stoer this
1795  // is gamma.
1796  Scalar b_(size_t i) const
1797  { return evalDerivative2_(/*x=*/Scalar(0.0), i)/2.0; }
1798 
1799  // returns the coefficient in front of the x^1 term. In Stoer this
1800  // is beta.
1801  Scalar c_(size_t i) const
1802  { return evalDerivative_(/*x=*/Scalar(0.0), i); }
1803 
1804  // returns the coefficient in front of the x^0 term. In Stoer this
1805  // is alpha.
1806  Scalar d_(size_t i) const
1807  { return eval_(/*x=*/Scalar(0.0), i); }
1808 
1809  Vector xPos_;
1810  Vector yPos_;
1811  Vector slopeVec_;
1812 };
1813 }
1814 
1815 #endif
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using C-style arrays.
Definition: Spline.hpp:324
Spline(size_t nSamples, const ScalarArray &x, const ScalarArray &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:200
void assignSamplingPoints_(DestVector &destX, DestVector &destY, const SourceVector &srcX, const SourceVector &srcY, unsigned nSamples)
Set the sampling point vectors.
Definition: Spline.hpp:1180
void setNumSamples_(size_t nSamples)
Resizes the internal vectors to store the sample points.
Definition: Spline.hpp:1093
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline&#39;s derivative at a given position.
Definition: Spline.hpp:837
void makeMonotonicSpline_(Vector &slopes)
Create a monotonic spline from the already set sampling points.
Definition: Spline.hpp:1400
void sortInput_()
Sort the sample points in ascending order of their x value.
Definition: Spline.hpp:1052
SplineType
The type of the spline to be created.
Definition: Spline.hpp:103
Helper class needed to sort the input sampling points.
Definition: Spline.hpp:1037
void setArrayOfPoints(size_t nSamples, const PointArray &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a C-style array.
Definition: Spline.hpp:593
int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED=false) const
Returns 1 if the spline is monotonically increasing, -1 if the spline is mononously decreasing and 0 ...
Definition: Spline.hpp:964
void makeFullSpline_(Scalar m0, Scalar m1)
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1105
Spline()
Default constructor for a spline.
Definition: Spline.hpp:115
Provides free functions to invert polynomials of degree 1, 2 and 3.
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using STL-compatible containers...
Definition: Spline.hpp:358
Evaluation intersect(const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d) const
Find the intersections of the spline with a cubic polynomial in the whole interval, throws Opm::MathError exception if there is more or less than one solution.
Definition: Spline.hpp:909
Scalar slope_(size_t i) const
Returns the slope (i.e.
Definition: Spline.hpp:1786
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition: TridiagonalMatrix.hpp:714
size_t rows() const
Return the number of rows of the matrix.
Definition: TridiagonalMatrix.hpp:140
Spline(const PointContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:184
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points natural spline using C-style arrays.
Definition: Spline.hpp:512
void reverseSamplingPoints_()
Reverse order of the elements in the arrays which contain the sampling points.
Definition: Spline.hpp:1080
Spline(const ScalarContainer &x, const ScalarContainer &y, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:171
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using STL-compatible containers.
Definition: Spline.hpp:553
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Spline.hpp:800
Spline(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:142
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left...
size_t intersectSegment_(Evaluation *sol, size_t segIdx, const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d, Scalar x0=-1e30, Scalar x1=1e30) const
Find all the intersections of a segment of the spline with a cubic polynomial within a specified inte...
Definition: Spline.hpp:1715
Scalar xAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:725
void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream &os=std::cout) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition: Spline.hpp:751
void setContainerOfPoints(const XYContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a STL-compatible container of array-like objects...
Definition: Spline.hpp:635
Scalar valueAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:731
void makePeriodicSpline_()
Create a periodic spline from the already set sampling points.
Definition: Spline.hpp:1147
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the spline as interval.
Definition: Spline.hpp:1030
Evaluation intersectInterval(Scalar x0, Scalar x1, const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d) const
Find the intersections of the spline with a cubic polynomial in a sub-interval of the spline...
Definition: Spline.hpp:924
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left...
Definition: TridiagonalMatrix.hpp:50
size_t numSamples() const
Return the number of (x, y) values.
Definition: Spline.hpp:719
Spline(const PointContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:251
Scalar h_(size_t i) const
Returns x[i] - x[i - 1].
Definition: Spline.hpp:1763
Spline(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:235
Spline(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)
Convenience constructor for a full spline with just two sampling points.
Definition: Spline.hpp:128
void makeFullSystem_(Matrix &M, Vector &d, Scalar m0, Scalar m1)
Make the linear system of equations Mx = d which results in the moments of the full spline...
Definition: Spline.hpp:1268
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline&#39;s third derivative at a given position.
Definition: Spline.hpp:891
void setContainerOfTuples(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:471
void setContainerOfTuples(const XYContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a STL-compatible container of tuple-like objects...
Definition: Spline.hpp:681
void assignFromTupleList_(DestVector &destX, DestVector &destY, ListIterator srcBegin, ListIterator srcEnd, unsigned nSamples)
Set the sampling points.
Definition: Spline.hpp:1234
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Spline.hpp:713
Spline(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:218
void setArrayOfPoints(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a C-style array.
Definition: Spline.hpp:392
Spline(size_t nSamples, const PointArray &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:157
void set(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)
Set the sampling points and the boundary slopes of the spline with two sampling points.
Definition: Spline.hpp:268
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline&#39;s second derivative at a given position.
Definition: Spline.hpp:867
void makePeriodicSystem_(Matrix &M, Vector &d)
Make the linear system of equations Mx = d which results in the moments of the periodic spline...
Definition: Spline.hpp:1338
void makeNaturalSpline_()
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1126
Class implementing cubic splines.
Definition: Spline.hpp:92
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
Convert the moments at the sample points to slopes.
Definition: Spline.hpp:1449
Scalar y_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1779
void setContainerOfPoints(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:429
void makeNaturalSystem_(Matrix &M, Vector &d)
Make the linear system of equations Mx = d which results in the moments of the natural spline...
Definition: Spline.hpp:1292
Scalar x_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1773