Main Lemma Repository
Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.

CubicSplineInterpolator.cpp 14KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346
  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 02/07/2014 12:50:52 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@lemmasoftware.org
  14. * @copyright Copyright (c) 2014,2018 Trevor Irons
  15. */
  16. #include "CubicSplineInterpolator.h"
  17. #include<iostream>
  18. #include<fstream>
  19. #include<vector>
  20. #include<algorithm>
  21. #include<cmath>
  22. namespace Lemma {
  23. // ==================== FRIEND METHODS =====================
  24. std::ostream &operator << (std::ostream &stream, const CubicSplineInterpolator &ob) {
  25. stream << ob.Serialize() << "\n---\n"; // End of doc ---
  26. return stream;
  27. }
  28. // ==================== LIFECYCLE =======================
  29. //--------------------------------------------------------------------------------------
  30. // Class: CubicSplineInterpolator
  31. // Method: CubicSplineInterpolator
  32. // Description: constructor (locked with ctor_key)
  33. //--------------------------------------------------------------------------------------
  34. CubicSplineInterpolator::CubicSplineInterpolator ( const ctor_key& key ) : LemmaObject( key ) {
  35. } // ----- end of method CubicSplineInterpolator::CubicSplineInterpolator (constructor) -----
  36. //--------------------------------------------------------------------------------------
  37. // Class: CubicSplineInterpolator
  38. // Method: CubicSplineInterpolator
  39. // Description: DeSerializing constructor (locked with ctor_key)
  40. //--------------------------------------------------------------------------------------
  41. CubicSplineInterpolator::CubicSplineInterpolator (const YAML::Node& node, const ctor_key& key ) : LemmaObject(node, key) {
  42. } // ----- end of method CubicSplineInterpolator::CubicSplineInterpolator (constructor) -----
  43. //--------------------------------------------------------------------------------------
  44. // Class: CubicSplineInterpolator
  45. // Method: NewSP()
  46. // Description: public constructor
  47. //--------------------------------------------------------------------------------------
  48. std::shared_ptr<CubicSplineInterpolator> CubicSplineInterpolator::NewSP() {
  49. return std::make_shared<CubicSplineInterpolator>( ctor_key() );
  50. }
  51. //--------------------------------------------------------------------------------------
  52. // Class: CubicSplineInterpolator
  53. // Method: Serialize
  54. //--------------------------------------------------------------------------------------
  55. YAML::Node CubicSplineInterpolator::Serialize ( ) const {
  56. YAML::Node node = LemmaObject::Serialize();;
  57. node.SetTag( GetName() );
  58. // FILL IN CLASS SPECIFICS HERE
  59. return node;
  60. } // ----- end of method CubicSplineInterpolator::Serialize -----
  61. //--------------------------------------------------------------------------------------
  62. // Class: CubicSplineInterpolator
  63. // Method: DeSerialize
  64. //--------------------------------------------------------------------------------------
  65. std::shared_ptr<CubicSplineInterpolator> CubicSplineInterpolator::DeSerialize ( const YAML::Node& node ) {
  66. if (node.Tag() != "CubicSplineInterpolator") {
  67. throw DeSerializeTypeMismatch( "CubicSplineInterpolator", node.Tag());
  68. }
  69. return std::make_shared<CubicSplineInterpolator>( node, ctor_key() );
  70. }
  71. //--------------------------------------------------------------------------------------
  72. // Class: CubicSplineInterpolator
  73. // Method: ~CubicSplineInterpolator
  74. // Description: destructor (protected)
  75. //--------------------------------------------------------------------------------------
  76. CubicSplineInterpolator::~CubicSplineInterpolator () {
  77. } // ----- end of method CubicSplineInterpolator::~CubicSplineInterpolator (destructor) -----
  78. // ==================== OPERATIONS =======================
  79. //--------------------------------------------------------------------------------------
  80. // Class: CubicSplineInterpolator
  81. // Method: SetKnots
  82. //--------------------------------------------------------------------------------------
  83. void CubicSplineInterpolator::SetKnots ( const VectorXr& x, const VectorXr& y ) {
  84. //int n = x.size()-1;
  85. Index n = x.size()-1;
  86. Spline = SplineSet(n);
  87. Spline.a = y;
  88. Spline.x = x;
  89. VectorXr h = VectorXr::Zero(n);
  90. for(Index i=0; i<n; ++i)
  91. h(i) = Spline.x(i+1)-Spline.x(i);
  92. VectorXr alpha(n-1);
  93. for(Index i=1; i<n-1; ++i)
  94. alpha(i) = 3.*(Spline.a[i+1]-Spline.a[i])/h[i] - 3.*(Spline.a[i]-Spline.a[i-1])/h[i-1] ;
  95. VectorXr l = VectorXr::Zero(n+1);
  96. VectorXr mu = VectorXr::Zero(n+1);
  97. VectorXr z = VectorXr::Zero(n+1);
  98. l[0] = 1;
  99. mu[0] = 0;
  100. z[0] = 0;
  101. for(Index i = 1; i < n-1; ++i) {
  102. l[i] = 2 *(Spline.x[i+1]-Spline.x[i-1])-h[i-1]*mu[i-1];
  103. mu[i] = h[i]/l[i];
  104. z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i];
  105. }
  106. l[n] = 1;
  107. z[n] = 0;
  108. for(Index j = n-1; j >= 0; --j) {
  109. Spline.c[j] = z[j] - mu[j] * Spline.c[j+1];
  110. Spline.b[j] = (Spline.a[j+1]-Spline.a[j])/h[j]-h[j]*(Spline.c[j+1]+2*Spline.c[j])/3;
  111. Spline.d[j] = (Spline.c[j+1]-Spline.c[j])/3/h[j];
  112. }
  113. // On OSX, this causes a strange bug 'sometimes', alignment?
  114. //Spline.c = Spline.c.head(n);
  115. return;
  116. } // ----- end of method CubicSplineInterpolator::SetKnots -----
  117. //--------------------------------------------------------------------------------------
  118. // Class: CubicSplineInterplator
  119. // Method: ResetKnotOrdinate
  120. //--------------------------------------------------------------------------------------
  121. void CubicSplineInterpolator::ResetKnotOrdinate ( const VectorXr& y ) {
  122. VectorXr x = Spline.x;
  123. SetKnots(x, y);
  124. return ;
  125. } // ----- end of method CubicSplineInterplator::ResetKnotOrdinate -----
  126. //--------------------------------------------------------------------------------------
  127. // Class: CubicSplineInterpolator
  128. // Method: InterpolateOrderedSet
  129. //--------------------------------------------------------------------------------------
  130. VectorXr CubicSplineInterpolator::InterpolateOrderedSet ( const VectorXr& x ) {
  131. VectorXr y = VectorXr::Zero(x.size());
  132. int ii = 0;
  133. for (int iy=0; iy<y.size(); ++iy) {
  134. y[iy] = Interpolate(x[iy], ii);
  135. }
  136. return y;
  137. } // ----- end of method CubicSplineInterpolator::InterpolateOrderedSet -----
  138. //--------------------------------------------------------------------------------------
  139. // Class: CubicSplineInterpolator
  140. // Method: Interpolate
  141. //--------------------------------------------------------------------------------------
  142. Real CubicSplineInterpolator::Interpolate ( const Real& x, int& i) {
  143. // O(n) search, could do bisection, but if these are sorted, then this is quick
  144. while(Spline.x[i] < x && i<Spline.x.size()) {
  145. ++i;
  146. }
  147. --i;
  148. // if ( x > Spline.x[i] ) {
  149. // std::cout << "DOOM\t" << x << "\t" << i << "\t" << Spline.x[i];
  150. // throw std::runtime_error("CubicSplineInterpolator::Interpolate ATTEMPT TO INTERPOLATE PAST LAST KNOT");
  151. // }
  152. return Spline.a[i] + Spline.b[i]*(x-Spline.x[i]) + Spline.c[i]*((x-Spline.x[i])*(x-Spline.x[i])) +
  153. Spline.d[i]*((x-Spline.x[i])*(x-Spline.x[i])*(x-Spline.x[i]) );
  154. } // ----- end of method CubicSplineInterpolator::Interpolate -----
  155. //--------------------------------------------------------------------------------------
  156. // Class: CubicSplineInterpolator
  157. // Method: Interpolate
  158. //--------------------------------------------------------------------------------------
  159. Real CubicSplineInterpolator::Interpolate ( const Real& x ) {
  160. int ii(0);
  161. return Interpolate(x, ii);
  162. }
  163. //--------------------------------------------------------------------------------------
  164. // Class: CubicSplineInterpolator
  165. // Method: Integrate
  166. //--------------------------------------------------------------------------------------
  167. Real CubicSplineInterpolator::Integrate ( const Real& x0, const Real& x1, const int& n ) {
  168. assert(n > 0);
  169. // force n to be even?
  170. //n += (n % 2);
  171. Real S = Interpolate(x0) + Interpolate(x1);
  172. Real h = (x1 - x0) / Real(n);
  173. int ik = 0;
  174. for (int i=1; i<n; i+=2) {
  175. S += 4. * Interpolate(x0 + (Real)(i)*h, ik);
  176. }
  177. ik = 0;
  178. for (int i=2; i<n-1; i+=2) {
  179. S += 2. * Interpolate(x0 + (Real)(i)*h, ik);
  180. }
  181. return S * h / 3.;
  182. }
  183. //--------------------------------------------------------------------------------------
  184. // Class: CubicSplineInterpolator
  185. // Method: Integrate
  186. //--------------------------------------------------------------------------------------
  187. Real CubicSplineInterpolator::Integrate ( const Real& x0, const Real& x1 ) {
  188. Index i0 = Interval(x0);
  189. Index i1 = Interval(x1);
  190. Real h0 = x0 - Spline.x(i0);
  191. if (mflag == -1) h0 = 0;
  192. Real h1 = x1 - Spline.x(i1);
  193. Real cubint = (((Spline.d(i1)*h1/4.0 + Spline.c(i1) )*h1/3.0 +
  194. Spline.b(i1) )*h1/2.0 + Spline.a(i1) )*h1
  195. - (((Spline.d(i0)*h0/4.0 + Spline.c(i0) )*h0/3.0 +
  196. Spline.b(i0) )*h0/2.0 + Spline.a(i0) )*h0;
  197. // Include integrals over intervening intervals.
  198. if (i1 > i0) {
  199. for (Index i=i0; i<i1-1; ++i) {
  200. Real h = Spline.x(i+1) - Spline.x(i);
  201. cubint += (((Spline.d(i)*h/4.0 + Spline.c(i) )*h/3.0 +
  202. Spline.b(i))*h/2.0 + Spline.a(i) )*h;
  203. }
  204. }
  205. return cubint;
  206. } // ----- end of method CubicSplineInterpolator::Integrate -----
  207. //--------------------------------------------------------------------------------------
  208. // Class: CubicSplineInterpolator
  209. // Method: Interval
  210. //--------------------------------------------------------------------------------------
  211. Index CubicSplineInterpolator::Interval ( const Real& x ) {
  212. std::cerr << "ENTERING CubicSplineInterpolator::Inverval. Known bugs here" << std::endl;
  213. Index nx = Spline.x.size() - 2; // TODO check if this is correct or just -1
  214. // when x not in range
  215. if (x <= Spline.x(0) || nx <= 1 ) {
  216. mflag = -1;
  217. return 1;
  218. }
  219. if (x >= Spline.x(nx)) {
  220. mflag = 1;
  221. return nx;
  222. }
  223. mflag = 0;
  224. if (ilo >= nx) ilo = nx-1;
  225. Index ihi = ilo+1;
  226. // if x is already in the interval
  227. if ( x<Spline.x(ihi) && x >= Spline.x(ilo) ) {
  228. //std::cout << "TRIVIAL INTERVAL " << Spline.x(ilo) << "\t" << x << "\t" << Spline.x(ihi) << std::endl;
  229. return ilo;
  230. }
  231. if (x <= Spline.x(ilo)) { // decrease ilo to capture
  232. Index istep = 1;
  233. for (Index ix=1; ix<nx; ++ix) {
  234. ihi = ilo;
  235. ilo = ihi - istep;
  236. ilo = std::max(Index(1), ilo);
  237. if (x >= Spline.x(ilo) || ilo == 1) break;
  238. istep *= 2;
  239. }
  240. } else if (x >= Spline.x(ihi)) { // increase ihi to capture
  241. Index istep = 1;
  242. for (Index ix=1; ix<nx; ++ix) {
  243. ilo = ihi;
  244. ihi = ilo + istep;
  245. ihi = std::min(ihi, nx);
  246. if (x <= Spline.x(ihi) || ihi == nx) break;
  247. istep *= 2;
  248. }
  249. }
  250. // Now Spline.x(ilo) <= x < Spline.x(ihi) --> Narrow the interval.
  251. //std::cout << "WIDE INTERVAL " << Spline.x(ilo) << "\t" << x << "\t" << Spline.x(ihi) << std::endl;
  252. for (Index ix=1; ix<nx; ++ix) {
  253. Index middle = (ilo+ihi) / 2;
  254. if (middle == ilo) break;
  255. if (x < Spline.x(middle)) {
  256. ihi = middle;
  257. } else {
  258. ilo = middle;
  259. }
  260. }
  261. assert ( Spline.x(ilo) < x && x < Spline.x(ihi) );
  262. return ilo;
  263. } // ----- end of method CubicSplineInterpolator::Inverval -----
  264. //--------------------------------------------------------------------------------------
  265. // Class: CubicSplineInterplator
  266. // Method: GetAbscissa
  267. //--------------------------------------------------------------------------------------
  268. VectorXr CubicSplineInterpolator::GetKnotAbscissa ( ) {
  269. return Spline.x;
  270. } // ----- end of method CubicSplineInterplator::get_GetAbscissa -----
  271. //--------------------------------------------------------------------------------------
  272. // Class: CubicSplineInterpolator
  273. // Method: GetKnotOrdinate
  274. //--------------------------------------------------------------------------------------
  275. VectorXr CubicSplineInterpolator::GetKnotOrdinate ( ) {
  276. return Spline.a;
  277. } // ----- end of method CubicSplineInterpolator::GetKnotOrdinate -----
  278. } // ----- end of Lemma name -----
  279. /* vim: set tabstop=4 expandtab: */
  280. /* vim: set filetype=cpp: */