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.

QWEKey.cpp 15KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332
  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. /* Original code is port of algorithm published by Key2011
  9. %------------------------------------------------------------------%
  10. % Copyright (c) 2012 by the Society of Exploration Geophysicists. %
  11. % For more information, go to http://software.seg.org/2012/0003 . %
  12. % You must read and accept usage terms at: %
  13. % http://software.seg.org/disclaimer.txt before use. %
  14. %------------------------------------------------------------------%
  15. */
  16. /**
  17. * @file
  18. * @date 02/12/2014 10:28:15 AM
  19. * @author Trevor Irons (ti)
  20. * @email Trevor.Irons@lemmasoftware.org
  21. * @copyright Copyright (c) 2014, Trevor Irons
  22. */
  23. #include "QWEKey.h"
  24. namespace Lemma {
  25. // ==================== FRIEND METHODS =====================
  26. std::ostream &operator<<(std::ostream &stream, const QWEKey &ob) {
  27. stream << ob.Serialize() << "\n---\n"; // End of doc ---
  28. return stream;
  29. }
  30. // ==================== LIFECYCLE =======================
  31. //--------------------------------------------------------------------------------------
  32. // Class: QWEKey
  33. // Method: QWEKey
  34. // Description: constructor (protected)
  35. //--------------------------------------------------------------------------------------
  36. //
  37. QWEKey::QWEKey (const ctor_key& ) : HankelTransform( ), RelTol(1e-12), AbsTol(1e-32), nQuad(61), nDelay(1),
  38. //QWEKey::QWEKey (const std::string& name) : HankelTransform(name), RelTol(1e-38), AbsTol(1e-48), nQuad(39), nDelay(5),
  39. nIntervalsMax(40) {
  40. BesselWeights( J0 ); // TODO experiment with zero weight (J0, J1) options, should be static one time method
  41. } // ----- end of method QWEKey::QWEKey (constructor) -----
  42. //--------------------------------------------------------------------------------------
  43. // Class: QWEKey
  44. // Method: QWEKey
  45. // Description: constructor (locked)
  46. //--------------------------------------------------------------------------------------
  47. QWEKey::QWEKey( const YAML::Node& node, const ctor_key& ) : HankelTransform(node) {
  48. }
  49. //--------------------------------------------------------------------------------------
  50. // Class: QWEKey
  51. // Method: New()
  52. // Description: public constructor
  53. //--------------------------------------------------------------------------------------
  54. std::shared_ptr<QWEKey> QWEKey::NewSP() {
  55. return std::make_shared<QWEKey>( ctor_key() );
  56. }
  57. //--------------------------------------------------------------------------------------
  58. // Class: QWEKey
  59. // Method: ~QWEKey
  60. // Description: destructor (protected)
  61. //--------------------------------------------------------------------------------------
  62. QWEKey::~QWEKey () {
  63. } // ----- end of method QWEKey::~QWEKey (destructor) -----
  64. //--------------------------------------------------------------------------------------
  65. // Class: QWEKey
  66. // Method: DeSerialize
  67. // Description: Factory method, converts YAML node into object
  68. //--------------------------------------------------------------------------------------
  69. std::shared_ptr<QWEKey> QWEKey::DeSerialize( const YAML::Node& node ) {
  70. if (node.Tag() != "QWEKey") {
  71. throw DeSerializeTypeMismatch( "QWEKey", node.Tag());
  72. }
  73. return std::make_shared<QWEKey> ( node, ctor_key() );
  74. }
  75. //--------------------------------------------------------------------------------------
  76. // Class: QWEKey
  77. // Method: Serialize
  78. // Description: Converts object into Serialized version
  79. //--------------------------------------------------------------------------------------
  80. YAML::Node QWEKey::Serialize() const {
  81. YAML::Node node = HankelTransform::Serialize();
  82. node.SetTag( GetName() );
  83. //node["LayerConductivity"] = LayerConductivity;
  84. return node;
  85. }
  86. //--------------------------------------------------------------------------------------
  87. // Class: QWEKey
  88. // Method: Zgauss
  89. //--------------------------------------------------------------------------------------
  90. Complex QWEKey::Zgauss ( const int &ikk, const EMMODE &imode,
  91. const int &itype, const Real &rho,
  92. const Real &wavef, KernelEM1DBase *Kernel ) {
  93. return Textrap(Kernel->GetManagerIndex(), Tn(Kernel->GetManagerIndex())) ;
  94. } // ----- end of method QWEKey::Zgauss -----
  95. //--------------------------------------------------------------------------------------
  96. // Class: QWEKey
  97. // Method: ComputeRelated
  98. //--------------------------------------------------------------------------------------
  99. void QWEKey::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DBase> Kernel ) {
  100. return ;
  101. } // ----- end of method QWEKey::ComputeRelated -----
  102. //--------------------------------------------------------------------------------------
  103. // Class: QWEKey
  104. // Method: ComputeRelated
  105. //--------------------------------------------------------------------------------------
  106. void QWEKey::ComputeRelated ( const Real& rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec ) {
  107. return ;
  108. } // ----- end of method QWEKey::ComputeRelated -----
  109. //--------------------------------------------------------------------------------------
  110. // Class: QWEKey
  111. // Method: ComputeRelated
  112. //--------------------------------------------------------------------------------------
  113. void QWEKey::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManagerIn ) {
  114. KernelManager = KernelManagerIn; // OK becauase this is internal and we know what we are doing
  115. Lambda = Bx.array()/rho;
  116. Intervals = xInt.array()/rho;
  117. int nrel = (int)(KernelManager->GetSTLVector().size());
  118. Zans = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  119. QWE(rho);
  120. return ;
  121. } // ----- end of method QWEKey::ComputeRelated -----
  122. //--------------------------------------------------------------------------------------
  123. // Class: QWEKey
  124. // Method: GaussQuadWeights
  125. //--------------------------------------------------------------------------------------
  126. void QWEKey::GaussQuadWeights(const int& N) {
  127. VectorXr Nv = VectorXr::LinSpaced(N-1, 1, N-1);
  128. VectorXr beta = 0.5 / (1.-(2.*Nv.array()).pow(-2)).sqrt();
  129. MatrixXr T = MatrixXr::Zero(N,N);
  130. //std::cerr << "Eigen ERROR BELOW, QWEKey.cpp QWEKey::GaussQuadWeights, COMMENTED OUT ";
  131. T.bottomLeftCorner(N-1, N-1) = beta.asDiagonal();
  132. Eigen::SelfAdjointEigenSolver<MatrixXr> eig( T.selfadjointView< Eigen::Lower >() );
  133. GaussAbscissa = eig.eigenvalues();
  134. GaussWeights = 2.*eig.eigenvectors().row(0).array().pow(2);
  135. }
  136. //--------------------------------------------------------------------------------------
  137. // Class: QWEKey
  138. // Method: BesselWeights
  139. //--------------------------------------------------------------------------------------
  140. void QWEKey::BesselWeights ( const sZeroType& sType ) {
  141. GaussQuadWeights(nQuad); // TODO should this be moved out of initializer?
  142. std::vector<Real> bz;
  143. xInt = VectorXr(nIntervalsMax+1);
  144. xInt(0) = 1e-20;
  145. switch (sType) {
  146. case J0:
  147. boost::math::cyl_bessel_j_zero(0.0, 1, nIntervalsMax, std::back_inserter(bz));
  148. xInt.tail(nIntervalsMax) = VectorXr::Map(&bz[0], nIntervalsMax);
  149. break;
  150. case J1:
  151. boost::math::cyl_bessel_j_zero(1.0, 1, nIntervalsMax, std::back_inserter(bz));
  152. xInt.tail(nIntervalsMax) = VectorXr::Map(&bz[0], nIntervalsMax);
  153. break;
  154. case NPI:
  155. xInt << 1e-20, VectorXr::LinSpaced(nIntervalsMax, 1, nIntervalsMax).array() * PI;
  156. break;
  157. }
  158. VectorXr dx = ( xInt.tail(nIntervalsMax) - xInt.head(nIntervalsMax) ).array() / 2.;
  159. // x = GaussAbscissa
  160. // dx in every row GaussWeights+1 rows, cols = n
  161. // dx[0] dx[1] ... dx[N] Gw[0] Gw[0] ... ndX
  162. // dx[0] dx[1] ... dx[N] Gw[1]
  163. MatrixXr Bxm = (dx.transpose().replicate(GaussAbscissa.size(), 1)).eval().array() *
  164. ((GaussAbscissa.replicate(1, dx.size()).array() + 1.));
  165. Bxm.array() += xInt.head(Bxm.cols()).transpose().replicate( Bxm.rows(), 1 ).array();
  166. Bx = VectorXr::Map( &Bxm(0,0), Bxm.size() );
  167. BJ0 = VectorXr(Bx.size());
  168. BJ1 = VectorXr(Bx.size());
  169. int iw = 0;
  170. for (int ii=0; ii<Bx.size(); ++ii) {
  171. BJ0(ii) = boost::math::cyl_bessel_j(0, Bx(ii)) * GaussWeights(iw);
  172. BJ1(ii) = boost::math::cyl_bessel_j(1, Bx(ii)) * GaussWeights(iw);
  173. ++iw;
  174. if (iw == GaussWeights.size()) iw = 0;
  175. }
  176. return ;
  177. } // ----- end of method QWEKey::BesselWeights -----
  178. //--------------------------------------------------------------------------------------
  179. // Class: QWEKey
  180. // Method: QWE
  181. //--------------------------------------------------------------------------------------
  182. void QWEKey::QWE ( const Real& rho ) {
  183. // TODO, is -1 needed here?
  184. int nTerms = nIntervalsMax - nDelay;// - 1;
  185. int nrel = (int)(KernelManager->GetSTLVector().size());
  186. // TODO GREMLINS LIVE IN HERE
  187. MatrixXcr prev = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  188. for (int i=0; i<nDelay; ++i) {
  189. getEyKernel(i, 0, rho);
  190. prev += Zans;
  191. }
  192. // Some of these are complex
  193. TS = MatrixXcr::Zero(nrel, nTerms);
  194. Tn = VectorXi::Zero(nrel);
  195. Textrap = MatrixXcr::Zero(nrel, nTerms);
  196. TrelErr = MatrixXr::Zero(nrel, nTerms);
  197. TabsErr = MatrixXr::Zero(nrel, nTerms);
  198. VectorXi Converged = VectorXi::Zero(nrel);
  199. // is nTerms right, 1 array shifting
  200. for (int i=nDelay; i<nTerms; ++i) {
  201. int n = i-nDelay;
  202. getEyKernel(i, 0, rho);
  203. for (int j=0; j<nrel; ++j) {
  204. if (!Converged(j)) { //continue;
  205. Tn(j) = n; // order of the expansion
  206. TS(j,n+1) = TS(j, n) + Zans(0, j); // working array for transformation
  207. /* Compute the Shanks transform using the Epsilon algorithm:
  208. Structured after Weniger (1989, p26) */
  209. /* TI - some kind ob bug here, shanks extrapolation doesn't buy much for TEM at least */
  210. Complex aux2(0);
  211. for (int k=n+1; k>0; --k) {
  212. Complex aux1 = aux2;
  213. aux2 = TS(j,k-1);
  214. Complex ddff = TS(j,k) - aux2;
  215. if (std::abs(ddff) < std::numeric_limits<Real>::min() ) {
  216. TS(j,k-1) = std::numeric_limits<Real>::max() ;
  217. } else {
  218. TS(j,k-1) = aux1 + 1./ddff;
  219. }
  220. }
  221. // The extrapolated result plus the prev integration term:
  222. Textrap(j,n) = TS(j, (n-1)%2)+prev(0, j);
  223. //Textrap(j,n) = TS(j, n%2 + 1)+prev(0, j);
  224. // Step 3: Analyze for convergence:
  225. if (n > 1) {
  226. TabsErr(j,n) = std::abs( Textrap(j, n) - Textrap(j, n-1));
  227. TrelErr(j,n) = TabsErr(j, n) / std::abs(Textrap(j, n)) ;
  228. Converged(j) = TrelErr(j,n) < RelTol + AbsTol/std::abs(Textrap(j,n));
  229. }
  230. }
  231. }
  232. if ( Converged.all() == 1 ) break;
  233. }
  234. // Trim up results
  235. // Clean up the T structure arrays? We can't really do this
  236. // because they are fixed size, maybe see how they are used and
  237. // init to zero. If they are only summed we are OK.
  238. /*
  239. for (int j = 0; j<nrel; ++j) {:nKernels
  240. n = Tn(j);
  241. T(j).extrap = T(j).extrap(1:n);
  242. T(j).relErr = T(j).relErr(1:n);
  243. T(j).absErr = T(j).absErr(1:n);
  244. }
  245. */
  246. return ;
  247. } // ----- end of method QWEKey::QWE -----
  248. //--------------------------------------------------------------------------------------
  249. // Class: QWEKey
  250. // Method: getEyKernel
  251. //--------------------------------------------------------------------------------------
  252. void QWEKey::getEyKernel ( const int& i, const int& idx, const Real& rho ) {
  253. int bidx = i*nQuad;
  254. int nrel = (int)(KernelManager->GetSTLVector().size());
  255. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork =
  256. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nQuad, nrel);
  257. for (int ik=0; ik<nQuad; ++ik) {
  258. KernelManager->ComputeReflectionCoeffs( Lambda(bidx+ik), idx, rho );
  259. for (int ir2=0; ir2<nrel; ++ir2) {
  260. // Zwork* needed due to sign convention (e^-jwt) of FT in filter weights
  261. Zwork(ik, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(Lambda(bidx+ik)));
  262. }
  263. }
  264. Real bma = (Intervals(i+1)-Intervals(i))/2;
  265. for (int ir2=0; ir2<nrel; ++ir2) {
  266. if (KernelManager->GetSTLVector()[ir2]->GetBesselOrder() == 0) {
  267. Zans(0, ir2) = bma * Zwork.col(ir2).dot( BJ0.segment(bidx, nQuad) ); // / rho;
  268. } else {
  269. Zans(0, ir2) = bma * Zwork.col(ir2).dot( BJ1.segment(bidx, nQuad) ); // / rho;
  270. }
  271. }
  272. // fcount += nQuad
  273. return ;
  274. } // ----- end of method QWEKey::getEyKernel -----
  275. void QWEKey::TestPrivate(const int& N) {
  276. //GaussQuadWeights(N);
  277. //std::cout << "abscissa\n" << GaussAbscissa << std::endl;
  278. //std::cout << "weights\n" << GaussWeights << std::endl;
  279. BesselWeights( J1 );
  280. //BesselZeros(0, N);
  281. std::cout << std::scientific;
  282. std::cout << "BJ0" << BJ0 << std::endl;
  283. std::cout << "BJ1" << BJ1 << std::endl;
  284. //std::cout << "Bess Zero\n" << xInt << std::endl;
  285. }
  286. } // ----- end of Lemma name -----