My Project
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 
33 
34 #include <ostream>
35 #include <vector>
36 #include <tuple>
37 
38 namespace Opm
39 {
89 template<class Scalar>
90 class Spline
91 {
93  typedef std::vector<Scalar> Vector;
94 
95 public:
101  enum SplineType {
102  Full,
103  Natural,
104  Periodic,
105  Monotonic
106  };
107 
114  { }
115 
126  Spline(Scalar x0, Scalar x1,
127  Scalar y0, Scalar y1,
128  Scalar m0, Scalar m1)
129  { set(x0, x1, y0, y1, m0, m1); }
130 
139  template <class ScalarArrayX, class ScalarArrayY>
140  Spline(size_t nSamples,
141  const ScalarArrayX& x,
142  const ScalarArrayY& y,
143  SplineType splineType = Natural,
144  bool sortInputs = true)
145  { this->setXYArrays(nSamples, x, y, splineType, sortInputs); }
146 
154  template <class PointArray>
155  Spline(size_t nSamples,
156  const PointArray& points,
157  SplineType splineType = Natural,
158  bool sortInputs = true)
159  { this->setArrayOfPoints(nSamples, points, splineType, sortInputs); }
160 
168  template <class ScalarContainer>
169  Spline(const ScalarContainer& x,
170  const ScalarContainer& y,
171  SplineType splineType = Natural,
172  bool sortInputs = true)
173  { this->setXYContainers(x, y, splineType, sortInputs); }
174 
181  template <class PointContainer>
182  Spline(const PointContainer& points,
183  SplineType splineType = Natural,
184  bool sortInputs = true)
185  { this->setContainerOfPoints(points, splineType, sortInputs); }
186 
197  template <class ScalarArray>
198  Spline(size_t nSamples,
199  const ScalarArray& x,
200  const ScalarArray& y,
201  Scalar m0,
202  Scalar m1,
203  bool sortInputs = true)
204  { this->setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
205 
215  template <class PointArray>
216  Spline(size_t nSamples,
217  const PointArray& points,
218  Scalar m0,
219  Scalar m1,
220  bool sortInputs = true)
221  { this->setArrayOfPoints(nSamples, points, m0, m1, sortInputs); }
222 
232  template <class ScalarContainerX, class ScalarContainerY>
233  Spline(const ScalarContainerX& x,
234  const ScalarContainerY& y,
235  Scalar m0,
236  Scalar m1,
237  bool sortInputs = true)
238  { this->setXYContainers(x, y, m0, m1, sortInputs); }
239 
248  template <class PointContainer>
249  Spline(const PointContainer& points,
250  Scalar m0,
251  Scalar m1,
252  bool sortInputs = true)
253  { this->setContainerOfPoints(points, m0, m1, sortInputs); }
254 
266  void set(Scalar x0, Scalar x1,
267  Scalar y0, Scalar y1,
268  Scalar m0, Scalar m1)
269  {
270  slopeVec_.resize(2);
271  xPos_.resize(2);
272  yPos_.resize(2);
273 
274  if (x0 > x1) {
275  xPos_[0] = x1;
276  xPos_[1] = x0;
277  yPos_[0] = y1;
278  yPos_[1] = y0;
279  }
280  else {
281  xPos_[0] = x0;
282  xPos_[1] = x1;
283  yPos_[0] = y0;
284  yPos_[1] = y1;
285  }
286 
287  slopeVec_[0] = m0;
288  slopeVec_[1] = m1;
289 
290  Matrix M(numSamples());
291  Vector d(numSamples());
292  Vector moments(numSamples());
293  this->makeFullSystem_(M, d, m0, m1);
294 
295  // solve for the moments
296  M.solve(moments, d);
297 
298  this->setSlopesFromMoments_(slopeVec_, moments);
299  }
300 
301 
305  // Full splines //
309 
321  template <class ScalarArrayX, class ScalarArrayY>
322  void setXYArrays(size_t nSamples,
323  const ScalarArrayX& x,
324  const ScalarArrayY& y,
325  Scalar m0, Scalar m1,
326  bool sortInputs = true)
327  {
328  assert(nSamples > 1);
329 
330  setNumSamples_(nSamples);
331  for (size_t i = 0; i < nSamples; ++i) {
332  xPos_[i] = x[i];
333  yPos_[i] = y[i];
334  }
335 
336  if (sortInputs)
337  sortInput_();
338  else if (xPos_[0] > xPos_[numSamples() - 1])
340 
341  makeFullSpline_(m0, m1);
342  }
343 
355  template <class ScalarContainerX, class ScalarContainerY>
356  void setXYContainers(const ScalarContainerX& x,
357  const ScalarContainerY& y,
358  Scalar m0, Scalar m1,
359  bool sortInputs = true)
360  {
361  assert(x.size() == y.size());
362  assert(x.size() > 1);
363 
364  setNumSamples_(x.size());
365 
366  std::copy(x.begin(), x.end(), xPos_.begin());
367  std::copy(y.begin(), y.end(), yPos_.begin());
368 
369  if (sortInputs)
370  sortInput_();
371  else if (xPos_[0] > xPos_[numSamples() - 1])
373 
374  makeFullSpline_(m0, m1);
375  }
376 
389  template <class PointArray>
390  void setArrayOfPoints(size_t nSamples,
391  const PointArray& points,
392  Scalar m0,
393  Scalar m1,
394  bool sortInputs = true)
395  {
396  // a spline with no or just one sampling points? what an
397  // incredible bad idea!
398  assert(nSamples > 1);
399 
400  setNumSamples_(nSamples);
401  for (size_t i = 0; i < nSamples; ++i) {
402  xPos_[i] = points[i][0];
403  yPos_[i] = points[i][1];
404  }
405 
406  if (sortInputs)
407  sortInput_();
408  else if (xPos_[0] > xPos_[numSamples() - 1])
410 
411  makeFullSpline_(m0, m1);
412  }
413 
426  template <class XYContainer>
427  void setContainerOfPoints(const XYContainer& points,
428  Scalar m0,
429  Scalar m1,
430  bool sortInputs = true)
431  {
432  // a spline with no or just one sampling points? what an
433  // incredible bad idea!
434  assert(points.size() > 1);
435 
436  setNumSamples_(points.size());
437  typename XYContainer::const_iterator it = points.begin();
438  typename XYContainer::const_iterator endIt = points.end();
439  for (size_t i = 0; it != endIt; ++i, ++it) {
440  xPos_[i] = (*it)[0];
441  yPos_[i] = (*it)[1];
442  }
443 
444  if (sortInputs)
445  sortInput_();
446  else if (xPos_[0] > xPos_[numSamples() - 1])
448 
449  // make a full spline
450  makeFullSpline_(m0, m1);
451  }
452 
468  template <class XYContainer>
469  void setContainerOfTuples(const XYContainer& points,
470  Scalar m0,
471  Scalar m1,
472  bool sortInputs = true)
473  {
474  // resize internal arrays
475  setNumSamples_(points.size());
476  typename XYContainer::const_iterator it = points.begin();
477  typename XYContainer::const_iterator endIt = points.end();
478  for (unsigned i = 0; it != endIt; ++i, ++it) {
479  xPos_[i] = std::get<0>(*it);
480  yPos_[i] = std::get<1>(*it);
481  }
482 
483  if (sortInputs)
484  sortInput_();
485  else if (xPos_[0] > xPos_[numSamples() - 1])
487 
488  // make a full spline
489  makeFullSpline_(m0, m1);
490  }
491 
495  // Natural/Periodic splines //
499 
509  template <class ScalarArrayX, class ScalarArrayY>
510  void setXYArrays(size_t nSamples,
511  const ScalarArrayX& x,
512  const ScalarArrayY& y,
513  SplineType splineType = Natural,
514  bool sortInputs = true)
515  {
516  assert(nSamples > 1);
517 
518  setNumSamples_(nSamples);
519  for (size_t i = 0; i < nSamples; ++i) {
520  xPos_[i] = x[i];
521  yPos_[i] = y[i];
522  }
523 
524  if (sortInputs)
525  sortInput_();
526  else if (xPos_[0] > xPos_[numSamples() - 1])
528 
529  if (splineType == Periodic)
531  else if (splineType == Natural)
533  else if (splineType == Monotonic)
534  this->makeMonotonicSpline_(slopeVec_);
535  else
536  throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
537  }
538 
550  template <class ScalarContainerX, class ScalarContainerY>
551  void setXYContainers(const ScalarContainerX& x,
552  const ScalarContainerY& y,
553  SplineType splineType = Natural,
554  bool sortInputs = true)
555  {
556  assert(x.size() == y.size());
557  assert(x.size() > 1);
558 
559  setNumSamples_(x.size());
560  std::copy(x.begin(), x.end(), xPos_.begin());
561  std::copy(y.begin(), y.end(), yPos_.begin());
562 
563  if (sortInputs)
564  sortInput_();
565  else if (xPos_[0] > xPos_[numSamples() - 1])
567 
568  if (splineType == Periodic)
570  else if (splineType == Natural)
572  else if (splineType == Monotonic)
573  this->makeMonotonicSpline_(slopeVec_);
574  else
575  throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
576  }
577 
590  template <class PointArray>
591  void setArrayOfPoints(size_t nSamples,
592  const PointArray& points,
593  SplineType splineType = Natural,
594  bool sortInputs = true)
595  {
596  // a spline with no or just one sampling points? what an
597  // incredible bad idea!
598  assert(nSamples > 1);
599 
600  setNumSamples_(nSamples);
601  for (size_t i = 0; i < nSamples; ++i) {
602  xPos_[i] = points[i][0];
603  yPos_[i] = points[i][1];
604  }
605 
606  if (sortInputs)
607  sortInput_();
608  else if (xPos_[0] > xPos_[numSamples() - 1])
610 
611  if (splineType == Periodic)
613  else if (splineType == Natural)
615  else if (splineType == Monotonic)
616  this->makeMonotonicSpline_(slopeVec_);
617  else
618  throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
619  }
620 
632  template <class XYContainer>
633  void setContainerOfPoints(const XYContainer& points,
634  SplineType splineType = Natural,
635  bool sortInputs = true)
636  {
637  // a spline with no or just one sampling points? what an
638  // incredible bad idea!
639  assert(points.size() > 1);
640 
641  setNumSamples_(points.size());
642  typename XYContainer::const_iterator it = points.begin();
643  typename XYContainer::const_iterator endIt = points.end();
644  for (size_t i = 0; it != endIt; ++ i, ++it) {
645  xPos_[i] = (*it)[0];
646  yPos_[i] = (*it)[1];
647  }
648 
649  if (sortInputs)
650  sortInput_();
651  else if (xPos_[0] > xPos_[numSamples() - 1])
653 
654  if (splineType == Periodic)
656  else if (splineType == Natural)
658  else if (splineType == Monotonic)
659  this->makeMonotonicSpline_(slopeVec_);
660  else
661  throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
662  }
663 
678  template <class XYContainer>
679  void setContainerOfTuples(const XYContainer& points,
680  SplineType splineType = Natural,
681  bool sortInputs = true)
682  {
683  // resize internal arrays
684  setNumSamples_(points.size());
685  typename XYContainer::const_iterator it = points.begin();
686  typename XYContainer::const_iterator endIt = points.end();
687  for (unsigned i = 0; it != endIt; ++i, ++it) {
688  xPos_[i] = std::get<0>(*it);
689  yPos_[i] = std::get<1>(*it);
690  }
691 
692  if (sortInputs)
693  sortInput_();
694  else if (xPos_[0] > xPos_[numSamples() - 1])
696 
697  if (splineType == Periodic)
699  else if (splineType == Natural)
701  else if (splineType == Monotonic)
702  this->makeMonotonicSpline_(slopeVec_);
703  else
704  throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
705  }
706 
710  template <class Evaluation>
711  bool applies(const Evaluation& x) const
712  { return x_(0) <= x && x <= x_(numSamples() - 1); }
713 
717  size_t numSamples() const
718  { return xPos_.size(); }
719 
723  Scalar xAt(size_t sampleIdx) const
724  { return x_(sampleIdx); }
725 
729  Scalar valueAt(size_t sampleIdx) const
730  { return y_(sampleIdx); }
731 
749  void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream& os = std::cout) const
750  {
751  Scalar x0 = std::min(xi0, xi1);
752  Scalar x1 = std::max(xi0, xi1);
753  const size_t n = numSamples() - 1;
754  for (size_t i = 0; i <= k; ++i) {
755  double x = i*(x1 - x0)/k + x0;
756  double x_p1 = x + (x1 - x0)/k;
757  double y;
758  double dy_dx;
759  double mono = 1;
760  if (!applies(x)) {
761  if (x < x_(0)) {
762  dy_dx = evalDerivative(x_(0));
763  y = (x - x_(0))*dy_dx + y_(0);
764  mono = (dy_dx>0)?1:-1;
765  }
766  else if (x > x_(n)) {
767  dy_dx = evalDerivative(x_(n));
768  y = (x - x_(n))*dy_dx + y_(n);
769  mono = (dy_dx>0)?1:-1;
770  }
771  else {
772  throw std::runtime_error("The sampling points given to a spline must be sorted by their x value!");
773  }
774  }
775  else {
776  y = eval(x);
777  dy_dx = evalDerivative(x);
778  mono = monotonic(x, x_p1, /*extrapolate=*/true);
779  }
780 
781  os << x << " " << y << " " << dy_dx << " " << mono << "\n";
782  }
783  }
784 
796  template <class Evaluation>
797  Evaluation eval(const Evaluation& x, bool extrapolate = false) const
798  {
799  if (!extrapolate && !applies(x))
800  throw NumericalIssue("Tried to evaluate a spline outside of its range");
801 
802  // handle extrapolation
803  if (extrapolate) {
804  if (x < xAt(0)) {
805  Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
806  Scalar y0 = y_(0);
807  return y0 + m*(x - xAt(0));
808  }
809  else if (x > xAt(static_cast<size_t>(static_cast<long int>(numSamples()) - 1))) {
810  Scalar m = evalDerivative_(xAt(static_cast<size_t>(numSamples() - 1)),
811  /*segmentIdx=*/static_cast<size_t>(numSamples()-2));
812  Scalar y0 = y_(static_cast<size_t>(numSamples() - 1));
813  return y0 + m*(x - xAt(static_cast<size_t>(numSamples() - 1)));
814  }
815  }
816 
817  return eval_(x, segmentIdx_(scalarValue(x)));
818  }
819 
832  template <class Evaluation>
833  Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
834  {
835  if (!extrapolate && !applies(x))
836  throw NumericalIssue("Tried to evaluate the derivative of a spline outside of its range");
837 
838  // handle extrapolation
839  if (extrapolate) {
840  if (x < xAt(0))
841  return evalDerivative_(xAt(0), /*segmentIdx=*/0);
842  else if (x > xAt(numSamples() - 1))
843  return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
844  }
845 
846  return evalDerivative_(x, segmentIdx_(scalarValue(x)));
847  }
848 
861  template <class Evaluation>
862  Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
863  {
864  if (!extrapolate && !applies(x))
865  throw NumericalIssue("Tried to evaluate the second derivative of a spline outside of its range");
866  else if (extrapolate)
867  return 0.0;
868 
869  return evalDerivative2_(x, segmentIdx_(scalarValue(x)));
870  }
871 
884  template <class Evaluation>
885  Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
886  {
887  if (!extrapolate && !applies(x))
888  throw NumericalIssue("Tried to evaluate the third derivative of a spline outside of its range");
889  else if (extrapolate)
890  return 0.0;
891 
892  return evalDerivative3_(x, segmentIdx_(scalarValue(x)));
893  }
894 
901  template <class Evaluation>
902  Evaluation intersect(const Evaluation& a,
903  const Evaluation& b,
904  const Evaluation& c,
905  const Evaluation& d) const
906  {
907  return intersectInterval(xAt(0), xAt(numSamples() - 1), a, b, c, d);
908  }
909 
916  template <class Evaluation>
917  Evaluation intersectInterval(Scalar x0, Scalar x1,
918  const Evaluation& a,
919  const Evaluation& b,
920  const Evaluation& c,
921  const Evaluation& d) const
922  {
923  assert(applies(x0) && applies(x1));
924 
925  Evaluation tmpSol[3], sol = 0;
926  size_t nSol = 0;
927  size_t iFirst = segmentIdx_(x0);
928  size_t iLast = segmentIdx_(x1);
929  for (size_t i = iFirst; i <= iLast; ++i)
930  {
931  size_t nCur = intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
932  if (nCur == 1)
933  sol = tmpSol[0];
934 
935  nSol += nCur;
936  if (nSol > 1) {
937  throw std::runtime_error("Spline has more than one intersection"); //<<a<<"x^3 + "<<b<"x^2 + "<<c<"x + "<<d);
938  }
939  }
940 
941  if (nSol != 1)
942  throw std::runtime_error("Spline has no intersection"); //<<a<"x^3 + " <<b<"x^2 + "<<c<"x + "<<d<<"!");
943 
944  return sol;
945  }
946 
955  int monotonic(Scalar x0, Scalar x1,
956  [[maybe_unused]] bool extrapolate = false) const
957  {
958  assert(std::abs(x0 - x1) > 1e-30);
959 
960  // make sure that x0 is smaller than x1
961  if (x0 > x1)
962  std::swap(x0, x1);
963 
964  assert(x0 < x1);
965 
966  int r = 3;
967  if (x0 < xAt(0)) {
968  assert(extrapolate);
969  Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
970  if (std::abs(m) < 1e-20)
971  r = (m < 0)?-1:1;
972  x0 = xAt(0);
973  };
974 
975  size_t i = segmentIdx_(x0);
976  if (x_(i + 1) >= x1) {
977  // interval is fully contained within a single spline
978  // segment
979  monotonic_(i, x0, x1, r);
980  return r;
981  }
982 
983  // the first segment overlaps with the specified interval
984  // partially
985  monotonic_(i, x0, x_(i+1), r);
986  ++ i;
987 
988  // make sure that the segments which are completly in the
989  // interval [x0, x1] all exhibit the same monotonicity.
990  size_t iEnd = segmentIdx_(x1);
991  for (; i < iEnd - 1; ++i) {
992  monotonic_(i, x_(i), x_(i + 1), r);
993  if (!r)
994  return 0;
995  }
996 
997  // if the user asked for a part of the spline which is
998  // extrapolated, we need to check the slope at the spline's
999  // endpoint
1000  if (x1 > xAt(numSamples() - 1)) {
1001  assert(extrapolate);
1002 
1003  Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
1004  if (m < 0)
1005  return (r < 0 || r==3)?-1:0;
1006  else if (m > 0)
1007  return (r > 0 || r==3)?1:0;
1008 
1009  return r;
1010  }
1011 
1012  // check for the last segment
1013  monotonic_(iEnd, x_(iEnd), x1, r);
1014 
1015  return r;
1016  }
1017 
1022  int monotonic() const
1023  { return monotonic(xAt(0), xAt(numSamples() - 1)); }
1024 
1025 protected:
1030  {
1031  ComparatorX_(const std::vector<Scalar>& x)
1032  : x_(x)
1033  {}
1034 
1035  bool operator ()(unsigned idxA, unsigned idxB) const
1036  { return x_.at(idxA) < x_.at(idxB); }
1037 
1038  const std::vector<Scalar>& x_;
1039  };
1040 
1044  void sortInput_()
1045  {
1046  size_t n = numSamples();
1047 
1048  // create a vector containing 0...n-1
1049  std::vector<unsigned> idxVector(n);
1050  for (unsigned i = 0; i < n; ++i)
1051  idxVector[i] = i;
1052 
1053  // sort the indices according to the x values of the sample
1054  // points
1055  ComparatorX_ cmp(xPos_);
1056  std::sort(idxVector.begin(), idxVector.end(), cmp);
1057 
1058  // reorder the sample points
1059  std::vector<Scalar> tmpX(n), tmpY(n);
1060  for (size_t i = 0; i < idxVector.size(); ++ i) {
1061  tmpX[i] = xPos_[idxVector[i]];
1062  tmpY[i] = yPos_[idxVector[i]];
1063  }
1064  xPos_ = tmpX;
1065  yPos_ = tmpY;
1066  }
1067 
1073  {
1074  // reverse the arrays
1075  size_t n = numSamples();
1076  for (unsigned i = 0; i <= (n - 1)/2; ++i) {
1077  std::swap(xPos_[i], xPos_[n - i - 1]);
1078  std::swap(yPos_[i], yPos_[n - i - 1]);
1079  }
1080  }
1081 
1085  void setNumSamples_(size_t nSamples)
1086  {
1087  xPos_.resize(nSamples);
1088  yPos_.resize(nSamples);
1089  slopeVec_.resize(nSamples);
1090  }
1091 
1097  void makeFullSpline_(Scalar m0, Scalar m1)
1098  {
1099  Matrix M(numSamples());
1100  std::vector<Scalar> d(numSamples());
1101  std::vector<Scalar> moments(numSamples());
1102 
1103  // create linear system of equations
1104  this->makeFullSystem_(M, d, m0, m1);
1105 
1106  // solve for the moments (-> second derivatives)
1107  M.solve(moments, d);
1108 
1109  // convert the moments to slopes at the sample points
1110  this->setSlopesFromMoments_(slopeVec_, moments);
1111  }
1112 
1119  {
1120  Matrix M(numSamples(), numSamples());
1121  Vector d(numSamples());
1122  Vector moments(numSamples());
1123 
1124  // create linear system of equations
1125  this->makeNaturalSystem_(M, d);
1126 
1127  // solve for the moments (-> second derivatives)
1128  M.solve(moments, d);
1129 
1130  // convert the moments to slopes at the sample points
1131  this->setSlopesFromMoments_(slopeVec_, moments);
1132  }
1133 
1140  {
1141  Matrix M(numSamples() - 1);
1142  Vector d(numSamples() - 1);
1143  Vector moments(numSamples() - 1);
1144 
1145  // create linear system of equations. This is a bit hacky,
1146  // because it assumes that std::vector internally stores its
1147  // data as a big C-style array, but it saves us from yet
1148  // another copy operation
1149  this->makePeriodicSystem_(M, d);
1150 
1151  // solve for the moments (-> second derivatives)
1152  M.solve(moments, d);
1153 
1154  moments.resize(numSamples());
1155  for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i) {
1156  unsigned ui = static_cast<unsigned>(i);
1157  moments[ui+1] = moments[ui];
1158  }
1159  moments[0] = moments[numSamples() - 1];
1160 
1161  // convert the moments to slopes at the sample points
1162  this->setSlopesFromMoments_(slopeVec_, moments);
1163  }
1164 
1171  template <class DestVector, class SourceVector>
1172  void assignSamplingPoints_(DestVector& destX,
1173  DestVector& destY,
1174  const SourceVector& srcX,
1175  const SourceVector& srcY,
1176  unsigned nSamples)
1177  {
1178  assert(nSamples >= 2);
1179 
1180  // copy sample points, make sure that the first x value is
1181  // smaller than the last one
1182  for (unsigned i = 0; i < nSamples; ++i) {
1183  unsigned idx = i;
1184  if (srcX[0] > srcX[nSamples - 1])
1185  idx = nSamples - i - 1;
1186  destX[i] = srcX[idx];
1187  destY[i] = srcY[idx];
1188  }
1189  }
1190 
1191  template <class DestVector, class ListIterator>
1192  void assignFromArrayList_(DestVector& destX,
1193  DestVector& destY,
1194  const ListIterator& srcBegin,
1195  const ListIterator& srcEnd,
1196  unsigned nSamples)
1197  {
1198  assert(nSamples >= 2);
1199 
1200  // find out wether the x values are in reverse order
1201  ListIterator it = srcBegin;
1202  ++it;
1203  bool reverse = false;
1204  if ((*srcBegin)[0] > (*it)[0])
1205  reverse = true;
1206  --it;
1207 
1208  // loop over all sampling points
1209  for (unsigned i = 0; it != srcEnd; ++i, ++it) {
1210  unsigned idx = i;
1211  if (reverse)
1212  idx = nSamples - i - 1;
1213  destX[idx] = (*it)[0];
1214  destY[idx] = (*it)[1];
1215  }
1216  }
1217 
1225  template <class DestVector, class ListIterator>
1226  void assignFromTupleList_(DestVector& destX,
1227  DestVector& destY,
1228  ListIterator srcBegin,
1229  ListIterator srcEnd,
1230  unsigned nSamples)
1231  {
1232  assert(nSamples >= 2);
1233 
1234  // copy sample points, make sure that the first x value is
1235  // smaller than the last one
1236 
1237  // find out wether the x values are in reverse order
1238  ListIterator it = srcBegin;
1239  ++it;
1240  bool reverse = false;
1241  if (std::get<0>(*srcBegin) > std::get<0>(*it))
1242  reverse = true;
1243  --it;
1244 
1245  // loop over all sampling points
1246  for (unsigned i = 0; it != srcEnd; ++i, ++it) {
1247  unsigned idx = i;
1248  if (reverse)
1249  idx = nSamples - i - 1;
1250  destX[idx] = std::get<0>(*it);
1251  destY[idx] = std::get<1>(*it);
1252  }
1253  }
1254 
1259  template <class Vector, class Matrix>
1260  void makeFullSystem_(Matrix& M, Vector& d, Scalar m0, Scalar m1)
1261  {
1262  makeNaturalSystem_(M, d);
1263 
1264  size_t n = numSamples() - 1;
1265  // first row
1266  M[0][1] = 1;
1267  d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0);
1268 
1269  // last row
1270  M[n][n - 1] = 1;
1271 
1272  // right hand side
1273  d[n] =
1274  6/h_(n)
1275  *
1276  (m1 - (y_(n) - y_(n - 1))/h_(n));
1277  }
1278 
1283  template <class Vector, class Matrix>
1284  void makeNaturalSystem_(Matrix& M, Vector& d)
1285  {
1286  M = 0.0;
1287 
1288  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1289  // Springer, 2005, p. 111
1290  size_t n = numSamples() - 1;
1291 
1292  // second to next to last rows
1293  for (size_t i = 1; i < n; ++i) {
1294  Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
1295  Scalar mu_i = 1 - lambda_i;
1296  Scalar d_i =
1297  6 / (h_(i) + h_(i + 1))
1298  *
1299  ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
1300 
1301  M[i][i-1] = mu_i;
1302  M[i][i] = 2;
1303  M[i][i + 1] = lambda_i;
1304  d[i] = d_i;
1305  };
1306 
1307  // See Stroer, equation (2.5.2.7)
1308  Scalar lambda_0 = 0;
1309  Scalar d_0 = 0;
1310 
1311  Scalar mu_n = 0;
1312  Scalar d_n = 0;
1313 
1314  // first row
1315  M[0][0] = 2;
1316  M[0][1] = lambda_0;
1317  d[0] = d_0;
1318 
1319  // last row
1320  M[n][n-1] = mu_n;
1321  M[n][n] = 2;
1322  d[n] = d_n;
1323  }
1324 
1329  template <class Matrix, class Vector>
1330  void makePeriodicSystem_(Matrix& M, Vector& d)
1331  {
1332  M = 0.0;
1333 
1334  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1335  // Springer, 2005, p. 111
1336  size_t n = numSamples() - 1;
1337 
1338  assert(M.rows() == n);
1339 
1340  // second to next to last rows
1341  for (size_t i = 2; i < n; ++i) {
1342  Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
1343  Scalar mu_i = 1 - lambda_i;
1344  Scalar d_i =
1345  6 / (h_(i) + h_(i + 1))
1346  *
1347  ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
1348 
1349  M[i-1][i-2] = mu_i;
1350  M[i-1][i-1] = 2;
1351  M[i-1][i] = lambda_i;
1352  d[i-1] = d_i;
1353  };
1354 
1355  Scalar lambda_n = h_(1) / (h_(n) + h_(1));
1356  Scalar lambda_1 = h_(2) / (h_(1) + h_(2));;
1357  Scalar mu_1 = 1 - lambda_1;
1358  Scalar mu_n = 1 - lambda_n;
1359 
1360  Scalar d_1 =
1361  6 / (h_(1) + h_(2))
1362  *
1363  ( (y_(2) - y_(1))/h_(2) - (y_(1) - y_(0))/h_(1));
1364  Scalar d_n =
1365  6 / (h_(n) + h_(1))
1366  *
1367  ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n));
1368 
1369 
1370  // first row
1371  M[0][0] = 2;
1372  M[0][1] = lambda_1;
1373  M[0][n-1] = mu_1;
1374  d[0] = d_1;
1375 
1376  // last row
1377  M[n-1][0] = lambda_n;
1378  M[n-1][n-2] = mu_n;
1379  M[n-1][n-1] = 2;
1380  d[n-1] = d_n;
1381  }
1382 
1391  template <class Vector>
1392  void makeMonotonicSpline_(Vector& slopes)
1393  {
1394  auto n = numSamples();
1395 
1396  // calculate the slopes of the secant lines
1397  std::vector<Scalar> delta(n);
1398  for (size_t k = 0; k < n - 1; ++k)
1399  delta[k] = (y_(k + 1) - y_(k))/(x_(k + 1) - x_(k));
1400 
1401  // calculate the "raw" slopes at the sample points
1402  for (size_t k = 1; k < n - 1; ++k)
1403  slopes[k] = (delta[k - 1] + delta[k])/2;
1404  slopes[0] = delta[0];
1405  slopes[n - 1] = delta[n - 2];
1406 
1407  // post-process the "raw" slopes at the sample points
1408  for (size_t k = 0; k < n - 1; ++k) {
1409  if (std::abs(delta[k]) < 1e-50) {
1410  // make the spline flat if the inputs are equal
1411  slopes[k] = 0;
1412  slopes[k + 1] = 0;
1413  ++ k;
1414  continue;
1415  }
1416  else {
1417  Scalar alpha = slopes[k] / delta[k];
1418  Scalar beta = slopes[k + 1] / delta[k];
1419 
1420  if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
1421  slopes[k] = 0;
1422  }
1423  // limit (alpha, beta) to a circle of radius 3
1424  else if (alpha*alpha + beta*beta > 3*3) {
1425  Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
1426  slopes[k] = tau*alpha*delta[k];
1427  slopes[k + 1] = tau*beta*delta[k];
1428  }
1429  }
1430  }
1431  }
1432 
1440  template <class MomentsVector, class SlopeVector>
1441  void setSlopesFromMoments_(SlopeVector& slopes, const MomentsVector& moments)
1442  {
1443  size_t n = numSamples();
1444 
1445  // evaluate slope at the rightmost point.
1446  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1447  // Springer, 2005, p. 109
1448  Scalar mRight;
1449 
1450  {
1451  Scalar h = this->h_(n - 1);
1452  Scalar x = h;
1453  //Scalar x_1 = 0;
1454 
1455  Scalar A =
1456  (y_(n - 1) - y_(n - 2))/h
1457  -
1458  h/6*(moments[n-1] - moments[n - 2]);
1459 
1460  mRight =
1461  //- moments[n - 2] * x_1*x_1 / (2 * h)
1462  //+
1463  moments[n - 1] * x*x / (2 * h)
1464  +
1465  A;
1466  }
1467 
1468  // evaluate the slope for the first n-1 sample points
1469  for (size_t i = 0; i < n - 1; ++ i) {
1470  // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1471  // Springer, 2005, p. 109
1472  Scalar h_i = this->h_(i + 1);
1473  //Scalar x_i = 0;
1474  Scalar x_i1 = h_i;
1475 
1476  Scalar A_i =
1477  (y_(i+1) - y_(i))/h_i
1478  -
1479  h_i/6*(moments[i+1] - moments[i]);
1480 
1481  slopes[i] =
1482  - moments[i] * x_i1*x_i1 / (2 * h_i)
1483  +
1484  //moments[i + 1] * x_i*x_i / (2 * h_i)
1485  //+
1486  A_i;
1487 
1488  }
1489  slopes[n - 1] = mRight;
1490  }
1491 
1492 
1493  // evaluate the spline at a given the position and given the
1494  // segment index
1495  template <class Evaluation>
1496  Evaluation eval_(const Evaluation& x, size_t i) const
1497  {
1498  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1499  Scalar delta = h_(i + 1);
1500  Evaluation t = (x - x_(i))/delta;
1501 
1502  return
1503  h00_(t) * y_(i)
1504  + h10_(t) * slope_(i)*delta
1505  + h01_(t) * y_(i + 1)
1506  + h11_(t) * slope_(i + 1)*delta;
1507  }
1508 
1509  // evaluate the derivative of a spline given the actual position
1510  // and the segment index
1511  template <class Evaluation>
1512  Evaluation evalDerivative_(const Evaluation& x, size_t i) const
1513  {
1514  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1515  Scalar delta = h_(i + 1);
1516  Evaluation t = (x - x_(i))/delta;
1517  Evaluation alpha = 1 / delta;
1518 
1519  return
1520  alpha *
1521  (h00_prime_(t) * y_(i)
1522  + h10_prime_(t) * slope_(i)*delta
1523  + h01_prime_(t) * y_(i + 1)
1524  + h11_prime_(t) * slope_(i + 1)*delta);
1525  }
1526 
1527  // evaluate the second derivative of a spline given the actual
1528  // position and the segment index
1529  template <class Evaluation>
1530  Evaluation evalDerivative2_(const Evaluation& x, size_t i) const
1531  {
1532  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1533  Scalar delta = h_(i + 1);
1534  Evaluation t = (x - x_(i))/delta;
1535  Evaluation alpha = 1 / delta;
1536 
1537  return
1538  alpha*alpha
1539  *(h00_prime2_(t) * y_(i)
1540  + h10_prime2_(t) * slope_(i)*delta
1541  + h01_prime2_(t) * y_(i + 1)
1542  + h11_prime2_(t) * slope_(i + 1)*delta);
1543  }
1544 
1545  // evaluate the third derivative of a spline given the actual
1546  // position and the segment index
1547  template <class Evaluation>
1548  Evaluation evalDerivative3_(const Evaluation& x, size_t i) const
1549  {
1550  // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1551  Scalar delta = h_(i + 1);
1552  Evaluation t = (x - x_(i))/delta;
1553  Evaluation alpha = 1 / delta;
1554 
1555  return
1556  alpha*alpha*alpha
1557  *(h00_prime3_(t)*y_(i)
1558  + h10_prime3_(t)*slope_(i)*delta
1559  + h01_prime3_(t)*y_(i + 1)
1560  + h11_prime3_(t)*slope_(i + 1)*delta);
1561  }
1562 
1563  // hermite basis functions
1564  template <class Evaluation>
1565  Evaluation h00_(const Evaluation& t) const
1566  { return (2*t - 3)*t*t + 1; }
1567 
1568  template <class Evaluation>
1569  Evaluation h10_(const Evaluation& t) const
1570  { return ((t - 2)*t + 1)*t; }
1571 
1572  template <class Evaluation>
1573  Evaluation h01_(const Evaluation& t) const
1574  { return (-2*t + 3)*t*t; }
1575 
1576  template <class Evaluation>
1577  Evaluation h11_(const Evaluation& t) const
1578  { return (t - 1)*t*t; }
1579 
1580  // first derivative of the hermite basis functions
1581  template <class Evaluation>
1582  Evaluation h00_prime_(const Evaluation& t) const
1583  { return (3*2*t - 2*3)*t; }
1584 
1585  template <class Evaluation>
1586  Evaluation h10_prime_(const Evaluation& t) const
1587  { return (3*t - 2*2)*t + 1; }
1588 
1589  template <class Evaluation>
1590  Evaluation h01_prime_(const Evaluation& t) const
1591  { return (-3*2*t + 2*3)*t; }
1592 
1593  template <class Evaluation>
1594  Evaluation h11_prime_(const Evaluation& t) const
1595  { return (3*t - 2)*t; }
1596 
1597  // second derivative of the hermite basis functions
1598  template <class Evaluation>
1599  Evaluation h00_prime2_(const Evaluation& t) const
1600  { return 2*3*2*t - 2*3; }
1601 
1602  template <class Evaluation>
1603  Evaluation h10_prime2_(const Evaluation& t) const
1604  { return 2*3*t - 2*2; }
1605 
1606  template <class Evaluation>
1607  Evaluation h01_prime2_(const Evaluation& t) const
1608  { return -2*3*2*t + 2*3; }
1609 
1610  template <class Evaluation>
1611  Evaluation h11_prime2_(const Evaluation& t) const
1612  { return 2*3*t - 2; }
1613 
1614  // third derivative of the hermite basis functions
1615  template <class Evaluation>
1616  Scalar h00_prime3_(const Evaluation&) const
1617  { return 2*3*2; }
1618 
1619  template <class Evaluation>
1620  Scalar h10_prime3_(const Evaluation&) const
1621  { return 2*3; }
1622 
1623  template <class Evaluation>
1624  Scalar h01_prime3_(const Evaluation&) const
1625  { return -2*3*2; }
1626 
1627  template <class Evaluation>
1628  Scalar h11_prime3_(const Evaluation&) const
1629  { return 2*3; }
1630 
1631  // returns the monotonicality of an interval of a spline segment
1632  //
1633  // The return value have the following meaning:
1634  //
1635  // 3: spline is constant within interval [x0, x1]
1636  // 1: spline is monotonously increasing in the specified interval
1637  // 0: spline is not monotonic (or constant) in the specified interval
1638  // -1: spline is monotonously decreasing in the specified interval
1639  int monotonic_(size_t i, Scalar x0, Scalar x1, int& r) const
1640  {
1641  // coefficients of derivative in monomial basis
1642  Scalar a = 3*a_(i);
1643  Scalar b = 2*b_(i);
1644  Scalar c = c_(i);
1645 
1646  if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
1647  return 3; // constant in interval, r stays unchanged!
1648 
1649  Scalar disc = b*b - 4*a*c;
1650  if (disc < 0) {
1651  // discriminant of derivative is smaller than 0, i.e. the
1652  // segment's derivative does not exhibit any extrema.
1653  if (x0*(x0*a + b) + c > 0) {
1654  r = (r==3 || r == 1)?1:0;
1655  return 1;
1656  }
1657  else {
1658  r = (r==3 || r == -1)?-1:0;
1659  return -1;
1660  }
1661  }
1662  disc = std::sqrt(disc);
1663  Scalar xE1 = (-b + disc)/(2*a);
1664  Scalar xE2 = (-b - disc)/(2*a);
1665 
1666  if (std::abs(disc) < 1e-30) {
1667  // saddle point -> no extrema
1668  if (std::abs(xE1 - x0) < 1e-30)
1669  // make sure that we're not picking the saddle point
1670  // to determine whether we're monotonically increasing
1671  // or decreasing
1672  x0 = x1;
1673  if (x0*(x0*a + b) + c > 0) {
1674  r = (r==3 || r == 1)?1:0;
1675  return 1;
1676  }
1677  else {
1678  r = (r==3 || r == -1)?-1:0;
1679  return -1;
1680  }
1681  };
1682  if ((x0 < xE1 && xE1 < x1) ||
1683  (x0 < xE2 && xE2 < x1))
1684  {
1685  // there is an extremum in the range (x0, x1)
1686  r = 0;
1687  return 0;
1688  }
1689  // no extremum in range (x0, x1)
1690  x0 = (x0 + x1)/2; // pick point in the middle of the interval
1691  // to avoid extrema on the boundaries
1692  if (x0*(x0*a + b) + c > 0) {
1693  r = (r==3 || r == 1)?1:0;
1694  return 1;
1695  }
1696  else {
1697  r = (r==3 || r == -1)?-1:0;
1698  return -1;
1699  }
1700  }
1701 
1706  template <class Evaluation>
1707  size_t intersectSegment_(Evaluation* sol,
1708  size_t segIdx,
1709  const Evaluation& a,
1710  const Evaluation& b,
1711  const Evaluation& c,
1712  const Evaluation& d,
1713  Scalar x0 = -1e30, Scalar x1 = 1e30) const
1714  {
1715  unsigned n =
1717  a_(segIdx) - a,
1718  b_(segIdx) - b,
1719  c_(segIdx) - c,
1720  d_(segIdx) - d);
1721  x0 = std::max(x_(segIdx), x0);
1722  x1 = std::min(x_(segIdx+1), x1);
1723 
1724  // filter the intersections outside of the specified interval
1725  size_t k = 0;
1726  for (unsigned j = 0; j < n; ++j) {
1727  if (x0 <= sol[j] && sol[j] <= x1) {
1728  sol[k] = sol[j];
1729  ++k;
1730  }
1731  }
1732  return k;
1733  }
1734 
1735  // find the segment index for a given x coordinate
1736  size_t segmentIdx_(Scalar x) const
1737  {
1738  // bisection
1739  size_t iLow = 0;
1740  size_t iHigh = numSamples() - 1;
1741 
1742  while (iLow + 1 < iHigh) {
1743  size_t i = (iLow + iHigh) / 2;
1744  if (x_(i) > x)
1745  iHigh = i;
1746  else
1747  iLow = i;
1748  };
1749  return iLow;
1750  }
1751 
1755  Scalar h_(size_t i) const
1756  {
1757  assert(x_(i) > x_(i-1)); // the sampling points must be given
1758  // in ascending order
1759  return x_(i) - x_(i - 1);
1760  }
1761 
1765  Scalar x_(size_t i) const
1766  { return xPos_[i]; }
1767 
1771  Scalar y_(size_t i) const
1772  { return yPos_[i]; }
1773 
1778  Scalar slope_(size_t i) const
1779  { return slopeVec_[i]; }
1780 
1781  // returns the coefficient in front of the x^3 term. In Stoer this
1782  // is delta.
1783  Scalar a_(size_t i) const
1784  { return evalDerivative3_(/*x=*/Scalar(0.0), i)/6.0; }
1785 
1786  // returns the coefficient in front of the x^2 term In Stoer this
1787  // is gamma.
1788  Scalar b_(size_t i) const
1789  { return evalDerivative2_(/*x=*/Scalar(0.0), i)/2.0; }
1790 
1791  // returns the coefficient in front of the x^1 term. In Stoer this
1792  // is beta.
1793  Scalar c_(size_t i) const
1794  { return evalDerivative_(/*x=*/Scalar(0.0), i); }
1795 
1796  // returns the coefficient in front of the x^0 term. In Stoer this
1797  // is alpha.
1798  Scalar d_(size_t i) const
1799  { return eval_(/*x=*/Scalar(0.0), i); }
1800 
1801  Vector xPos_;
1802  Vector yPos_;
1803  Vector slopeVec_;
1804 };
1805 }
1806 
1807 #endif
Provides the opm-material specific exception classes.
Provides free functions to invert polynomials of degree 1, 2 and 3.
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition: Exceptions.hpp:46
Class implementing cubic splines.
Definition: Spline.hpp:91
size_t numSamples() const
Return the number of (x, y) values.
Definition: Spline.hpp:717
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:1260
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:591
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:322
Scalar valueAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:729
void makeNaturalSpline_()
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1118
void sortInput_()
Sort the sample points in ascending order of their x value.
Definition: Spline.hpp:1044
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:140
Spline(const PointContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:249
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
Convert the moments at the sample points to slopes.
Definition: Spline.hpp:1441
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:469
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:1330
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Spline.hpp:711
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:169
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:126
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,...
Definition: Spline.hpp:902
Scalar y_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1771
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Spline.hpp:833
SplineType
The type of the spline to be created.
Definition: Spline.hpp:101
Spline(const PointContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:182
int monotonic(Scalar x0, Scalar x1, [[maybe_unused]] bool extrapolate=false) const
Returns 1 if the spline is monotonically increasing, -1 if the spline is mononously decreasing and 0 ...
Definition: Spline.hpp:955
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:1284
Scalar h_(size_t i) const
Returns x[i] - x[i - 1].
Definition: Spline.hpp:1755
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:266
Spline(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:233
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:633
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's third derivative at a given position.
Definition: Spline.hpp:885
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:749
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:198
void assignSamplingPoints_(DestVector &destX, DestVector &destY, const SourceVector &srcX, const SourceVector &srcY, unsigned nSamples)
Set the sampling point vectors.
Definition: Spline.hpp:1172
Scalar xAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:723
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the spline as interval.
Definition: Spline.hpp:1022
void assignFromTupleList_(DestVector &destX, DestVector &destY, ListIterator srcBegin, ListIterator srcEnd, unsigned nSamples)
Set the sampling points.
Definition: Spline.hpp:1226
Spline(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:216
Scalar slope_(size_t i) const
Returns the slope (i.e.
Definition: Spline.hpp:1778
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Spline.hpp:797
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:1707
Scalar x_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1765
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:390
void makeFullSpline_(Scalar m0, Scalar m1)
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1097
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's second derivative at a given position.
Definition: Spline.hpp:862
void makeMonotonicSpline_(Vector &slopes)
Create a monotonic spline from the already set sampling points.
Definition: Spline.hpp:1392
void setNumSamples_(size_t nSamples)
Resizes the internal vectors to store the sample points.
Definition: Spline.hpp:1085
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:510
Spline()
Default constructor for a spline.
Definition: Spline.hpp:113
void reverseSamplingPoints_()
Reverse order of the elements in the arrays which contain the sampling points.
Definition: Spline.hpp:1072
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:679
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:551
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:356
void makePeriodicSpline_()
Create a periodic spline from the already set sampling points.
Definition: Spline.hpp:1139
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:427
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:917
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:155
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition: TridiagonalMatrix.hpp:51
size_t rows() const
Return the number of rows of the matrix.
Definition: TridiagonalMatrix.hpp:140
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition: TridiagonalMatrix.hpp:714
Helper class needed to sort the input sampling points.
Definition: Spline.hpp:1030