Main Lemma Repository
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

CubicSplineInterpolator.h 6.1KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 09/25/2013 08:20:14 AM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@lemmasoftware.org
  14. * @copyright Copyright (c) 2013,2018 Trevor Irons
  15. */
  16. #ifndef CUBICSPLINEINTERPOLATOR_INC
  17. #define CUBICSPLINEINTERPOLATOR_INC
  18. #include "LemmaObject.h"
  19. namespace Lemma {
  20. /**
  21. * \brief Simple struct to hold spline terms
  22. * \details This struct holds the actual knot and spline interpolation terms.
  23. */
  24. struct SplineSet{
  25. VectorXr a;
  26. VectorXr b;
  27. VectorXr c;
  28. VectorXr d;
  29. VectorXr x;
  30. SplineSet( ) {
  31. }
  32. SplineSet(const Index& n) {
  33. a = VectorXr::Zero(n+1);
  34. b = VectorXr::Zero(n);
  35. c = VectorXr::Zero(n+1);
  36. d = VectorXr::Zero(n);
  37. x = VectorXr::Zero(n+1);
  38. }
  39. };
  40. /**
  41. * \ingroup LemmaCore
  42. * \brief Real 1D Natural cubic spline interpolator.
  43. * \details Splines are fit between knots \f$j\f$ according to the forulae
  44. * \f[ S_j(x) = a_j + b_j(x - x_j) + c_j(x-x_j)^2 + d_j(x-x_y)^3 \f]
  45. * The spline must satisfy the following conditions
  46. * \f{eqnarray} {
  47. * S_i(x_i) & = & y_i = S_{i-1}(x_i), i = 1,..., n-1 \\
  48. * S'_i(x_i) & = & S'_{i-1}(x_i), i = 1,..., n-1 \\
  49. * S''_i(x_i) & = & S''_{i-1}(x_i), i = 1,..., n-1 \\
  50. * S''_0(x_0) & = & S''_{n-1}(x_n) = 0
  51. * \f}
  52. */
  53. class CubicSplineInterpolator : public LemmaObject {
  54. friend std::ostream &operator<<(std::ostream &stream,
  55. const CubicSplineInterpolator& ob);
  56. //struct ctor_key {};
  57. public:
  58. // ==================== LIFECYCLE =======================
  59. /** Default constructor */
  60. explicit CubicSplineInterpolator ( const ctor_key& );
  61. /** DeSerializing constructor, locked use factory DeSerialize method*/
  62. CubicSplineInterpolator ( const YAML::Node& node, const ctor_key& );
  63. /** Destructor use smart pointers to auto delete */
  64. virtual ~CubicSplineInterpolator ();
  65. /**
  66. * Factory method for generating concrete class.
  67. * @return a std::shared_ptr of type CubicSplineInterpolator
  68. */
  69. static std::shared_ptr<CubicSplineInterpolator> NewSP();
  70. /**
  71. * Uses YAML to serialize this object.
  72. * @return a YAML::Node
  73. */
  74. virtual YAML::Node Serialize() const;
  75. /**
  76. * Constructs an object from a YAML::Node.
  77. */
  78. static std::shared_ptr< CubicSplineInterpolator > DeSerialize(const YAML::Node& node);
  79. /**
  80. * Constructs an object from a string representation of a YAML::Node. This is primarily
  81. * used in Python wrapping
  82. */
  83. static std::shared_ptr<CubicSplineInterpolator> DeSerialize( const std::string& node ) {
  84. return CubicSplineInterpolator::DeSerialize(YAML::LoadFile(node));
  85. }
  86. // ==================== OPERATORS =======================
  87. // ==================== OPERATIONS =======================
  88. /** Sets the knots to use for interpolation.
  89. @param[in] x are the absissa values
  90. @param[in] y are the ordinate values
  91. */
  92. void SetKnots(const VectorXr& x, const VectorXr& y);
  93. /** Resets the knots to use for interpolation, when abscissa values haven't changed.
  94. @param[in] y are the ordinate values
  95. */
  96. void ResetKnotOrdinate( const VectorXr& y );
  97. /** Interpolate a monotonically increasing ordered set.
  98. @param[in] x are the interpolation abscissa points
  99. @return the ordinate values at x
  100. */
  101. VectorXr InterpolateOrderedSet(const VectorXr& x);
  102. /** integrates the spline from x0 to x1. Uses composite Simpson's rule and n is the number of segments
  103. * @param[in] x0 is left argument
  104. * @param[in] x1 is right argument
  105. * @param[in] n is the number of points, must be even
  106. */
  107. Real Integrate(const Real& x0, const Real& x1, const int& n);
  108. /** integrates using cubic spline values. Taken from AMRIRA P223F project code Leroi, which in turn was based on
  109. This is a modification of the FUNCTION PPVALU in the book
  110. "A PRACTICAL GUIDE TO SPLINES" by C. DE BOOR
  111. */
  112. Real Integrate(const Real& x0, const Real& x1);
  113. /** @returns the know abscissa values
  114. */
  115. VectorXr GetKnotAbscissa();
  116. /** @returns the know abscissa values
  117. */
  118. VectorXr GetKnotOrdinate();
  119. /** Interpolation at a single point.
  120. @param[in] x is the interpolation abscissa point
  121. @param[in] i is an optional index to start searching at. Defaults to zero
  122. @return the ordinate value at x
  123. */
  124. Real Interpolate(const Real& x, int& i);
  125. /** Interpolation at a single point.
  126. @param[in] x is the interpolation abscissa point
  127. @return the ordinate value at x
  128. */
  129. Real Interpolate(const Real& x);
  130. // ==================== ACCESS =======================
  131. // ==================== INQUIRY =======================
  132. /** Returns the name of the underlying class, similiar to Python's type */
  133. virtual inline std::string GetName() const {
  134. return CName;
  135. }
  136. protected:
  137. // ==================== OPERATIONS =======================
  138. /** Finds the interval of knots in spline to use for integration.
  139. */
  140. Index Interval(const Real& x);
  141. private:
  142. /** Copy */
  143. CubicSplineInterpolator( const CubicSplineInterpolator& ) = delete;
  144. /** ASCII string representation of the class name */
  145. static constexpr auto CName = "CubicSplineInterpolator";
  146. SplineSet Spline;
  147. Index ilo;
  148. int mflag;
  149. // ==================== DATA MEMBERS =========================
  150. }; // ----- end of class CubicSplineInterpolator -----
  151. } // ----- end of Lemma name -----
  152. #endif // ----- #ifndef CUBICSPLINEINTERPOLATOR_INC -----
  153. /* vim: set tabstop=4 expandtab: */
  154. /* vim: set filetype=cpp syntax=cpp.doxygen: */