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

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347
  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: GetName
  89. // Description: Class identifier
  90. //--------------------------------------------------------------------------------------
  91. inline std::string QWEKey::GetName ( ) const {
  92. return CName;
  93. } // ----- end of method QWEKey::GetName -----
  94. //--------------------------------------------------------------------------------------
  95. // Class: QWEKey
  96. // Method: Zgauss
  97. //--------------------------------------------------------------------------------------
  98. Complex QWEKey::Zgauss ( const int &ikk, const EMMODE &imode,
  99. const int &itype, const Real &rho,
  100. const Real &wavef, KernelEM1DBase *Kernel ) {
  101. return Textrap(Kernel->GetManagerIndex(), Tn(Kernel->GetManagerIndex())) ;
  102. } // ----- end of method QWEKey::Zgauss -----
  103. //--------------------------------------------------------------------------------------
  104. // Class: QWEKey
  105. // Method: ComputeRelated
  106. //--------------------------------------------------------------------------------------
  107. void QWEKey::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DBase> Kernel ) {
  108. return ;
  109. } // ----- end of method QWEKey::ComputeRelated -----
  110. //--------------------------------------------------------------------------------------
  111. // Class: QWEKey
  112. // Method: ComputeRelated
  113. //--------------------------------------------------------------------------------------
  114. void QWEKey::ComputeRelated ( const Real& rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec ) {
  115. return ;
  116. } // ----- end of method QWEKey::ComputeRelated -----
  117. //--------------------------------------------------------------------------------------
  118. // Class: QWEKey
  119. // Method: ComputeRelated
  120. //--------------------------------------------------------------------------------------
  121. void QWEKey::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManagerIn ) {
  122. KernelManager = KernelManagerIn; // OK becauase this is internal and we know what we are doing
  123. Lambda = Bx.array()/rho;
  124. Intervals = xInt.array()/rho;
  125. int nrel = (int)(KernelManager->GetSTLVector().size());
  126. Zans = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  127. QWE(rho);
  128. return ;
  129. } // ----- end of method QWEKey::ComputeRelated -----
  130. //--------------------------------------------------------------------------------------
  131. // Class: QWEKey
  132. // Method: GaussQuadWeights
  133. //--------------------------------------------------------------------------------------
  134. void QWEKey::GaussQuadWeights(const int& N) {
  135. VectorXr Nv = VectorXr::LinSpaced(N-1, 1, N-1);
  136. VectorXr beta = 0.5 / (1.-(2.*Nv.array()).pow(-2)).sqrt();
  137. MatrixXr T = MatrixXr::Zero(N,N);
  138. //std::cerr << "Eigen ERROR BELOW, QWEKey.cpp QWEKey::GaussQuadWeights, COMMENTED OUT ";
  139. T.bottomLeftCorner(N-1, N-1) = beta.asDiagonal();
  140. Eigen::SelfAdjointEigenSolver<MatrixXr> eig( T.selfadjointView< Eigen::Lower >() );
  141. GaussAbscissa = eig.eigenvalues();
  142. GaussWeights = 2.*eig.eigenvectors().row(0).array().pow(2);
  143. }
  144. //--------------------------------------------------------------------------------------
  145. // Class: QWEKey
  146. // Method: BesselWeights
  147. //--------------------------------------------------------------------------------------
  148. #ifdef HAVE_BOOST_SPECIAL_FUNCTIONS
  149. void QWEKey::BesselWeights ( const sZeroType& sType ) {
  150. GaussQuadWeights(nQuad); // TODO should this be moved out of initializer?
  151. std::vector<Real> bz;
  152. xInt = VectorXr(nIntervalsMax+1);
  153. xInt(0) = 1e-20;
  154. switch (sType) {
  155. case J0:
  156. boost::math::cyl_bessel_j_zero(0.0, 1, nIntervalsMax, std::back_inserter(bz));
  157. xInt.tail(nIntervalsMax) = VectorXr::Map(&bz[0], nIntervalsMax);
  158. break;
  159. case J1:
  160. boost::math::cyl_bessel_j_zero(1.0, 1, nIntervalsMax, std::back_inserter(bz));
  161. xInt.tail(nIntervalsMax) = VectorXr::Map(&bz[0], nIntervalsMax);
  162. break;
  163. case NPI:
  164. xInt << 1e-20, VectorXr::LinSpaced(nIntervalsMax, 1, nIntervalsMax).array() * PI;
  165. break;
  166. }
  167. VectorXr dx = ( xInt.tail(nIntervalsMax) - xInt.head(nIntervalsMax) ).array() / 2.;
  168. // x = GaussAbscissa
  169. // dx in every row GaussWeights+1 rows, cols = n
  170. // dx[0] dx[1] ... dx[N] Gw[0] Gw[0] ... ndX
  171. // dx[0] dx[1] ... dx[N] Gw[1]
  172. MatrixXr Bxm = (dx.transpose().replicate(GaussAbscissa.size(), 1)).eval().array() *
  173. ((GaussAbscissa.replicate(1, dx.size()).array() + 1.));
  174. Bxm.array() += xInt.head(Bxm.cols()).transpose().replicate( Bxm.rows(), 1 ).array();
  175. Bx = VectorXr::Map( &Bxm(0,0), Bxm.size() );
  176. BJ0 = VectorXr(Bx.size());
  177. BJ1 = VectorXr(Bx.size());
  178. int iw = 0;
  179. for (int ii=0; ii<Bx.size(); ++ii) {
  180. BJ0(ii) = boost::math::cyl_bessel_j(0, Bx(ii)) * GaussWeights(iw);
  181. BJ1(ii) = boost::math::cyl_bessel_j(1, Bx(ii)) * GaussWeights(iw);
  182. ++iw;
  183. if (iw == GaussWeights.size()) iw = 0;
  184. }
  185. return ;
  186. } // ----- end of method QWEKey::BesselWeights -----
  187. #else
  188. void QWEKey::BesselWeights ( const sZeroType& sType ) {
  189. std::cerr << "QWEKey requires boost special functions";
  190. }
  191. #endif
  192. //--------------------------------------------------------------------------------------
  193. // Class: QWEKey
  194. // Method: QWE
  195. //--------------------------------------------------------------------------------------
  196. void QWEKey::QWE ( const Real& rho ) {
  197. // TODO, is -1 needed here?
  198. int nTerms = nIntervalsMax - nDelay;// - 1;
  199. int nrel = (int)(KernelManager->GetSTLVector().size());
  200. // TODO GREMLINS LIVE IN HERE
  201. MatrixXcr prev = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  202. for (int i=0; i<nDelay; ++i) {
  203. getEyKernel(i, 0, rho);
  204. prev += Zans;
  205. }
  206. // Some of these are complex
  207. TS = MatrixXcr::Zero(nrel, nTerms);
  208. Tn = VectorXi::Zero(nrel);
  209. Textrap = MatrixXcr::Zero(nrel, nTerms);
  210. TrelErr = MatrixXr::Zero(nrel, nTerms);
  211. TabsErr = MatrixXr::Zero(nrel, nTerms);
  212. VectorXi Converged = VectorXi::Zero(nrel);
  213. // is nTerms right, 1 array shifting
  214. for (int i=nDelay; i<nTerms; ++i) {
  215. int n = i-nDelay;
  216. getEyKernel(i, 0, rho);
  217. for (int j=0; j<nrel; ++j) {
  218. if (!Converged(j)) { //continue;
  219. Tn(j) = n; // order of the expansion
  220. TS(j,n+1) = TS(j, n) + Zans(0, j); // working array for transformation
  221. /* Compute the Shanks transform using the Epsilon algorithm:
  222. Structured after Weniger (1989, p26) */
  223. /* TI - some kind ob bug here, shanks extrapolation doesn't buy much for TEM at least */
  224. Complex aux2(0);
  225. for (int k=n+1; k>0; --k) {
  226. Complex aux1 = aux2;
  227. aux2 = TS(j,k-1);
  228. Complex ddff = TS(j,k) - aux2;
  229. if (std::abs(ddff) < std::numeric_limits<Real>::min() ) {
  230. TS(j,k-1) = std::numeric_limits<Real>::max() ;
  231. } else {
  232. TS(j,k-1) = aux1 + 1./ddff;
  233. }
  234. }
  235. // The extrapolated result plus the prev integration term:
  236. Textrap(j,n) = TS(j, (n-1)%2)+prev(0, j);
  237. //Textrap(j,n) = TS(j, n%2 + 1)+prev(0, j);
  238. // Step 3: Analyze for convergence:
  239. if (n > 1) {
  240. TabsErr(j,n) = std::abs( Textrap(j, n) - Textrap(j, n-1));
  241. TrelErr(j,n) = TabsErr(j, n) / std::abs(Textrap(j, n)) ;
  242. Converged(j) = TrelErr(j,n) < RelTol + AbsTol/std::abs(Textrap(j,n));
  243. }
  244. }
  245. }
  246. if ( Converged.all() == 1 ) break;
  247. }
  248. // Trim up results
  249. // Clean up the T structure arrays? We can't really do this
  250. // because they are fixed size, maybe see how they are used and
  251. // init to zero. If they are only summed we are OK.
  252. /*
  253. for (int j = 0; j<nrel; ++j) {:nKernels
  254. n = Tn(j);
  255. T(j).extrap = T(j).extrap(1:n);
  256. T(j).relErr = T(j).relErr(1:n);
  257. T(j).absErr = T(j).absErr(1:n);
  258. }
  259. */
  260. return ;
  261. } // ----- end of method QWEKey::QWE -----
  262. //--------------------------------------------------------------------------------------
  263. // Class: QWEKey
  264. // Method: getEyKernel
  265. //--------------------------------------------------------------------------------------
  266. void QWEKey::getEyKernel ( const int& i, const int& idx, const Real& rho ) {
  267. int bidx = i*nQuad;
  268. int nrel = (int)(KernelManager->GetSTLVector().size());
  269. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork =
  270. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nQuad, nrel);
  271. for (int ik=0; ik<nQuad; ++ik) {
  272. KernelManager->ComputeReflectionCoeffs( Lambda(bidx+ik), idx, rho );
  273. for (int ir2=0; ir2<nrel; ++ir2) {
  274. // Zwork* needed due to sign convention (e^-jwt) of FT in filter weights
  275. Zwork(ik, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(Lambda(bidx+ik)));
  276. }
  277. }
  278. Real bma = (Intervals(i+1)-Intervals(i))/2;
  279. for (int ir2=0; ir2<nrel; ++ir2) {
  280. if (KernelManager->GetSTLVector()[ir2]->GetBesselOrder() == 0) {
  281. Zans(0, ir2) = bma * Zwork.col(ir2).dot( BJ0.segment(bidx, nQuad) ); // / rho;
  282. } else {
  283. Zans(0, ir2) = bma * Zwork.col(ir2).dot( BJ1.segment(bidx, nQuad) ); // / rho;
  284. }
  285. }
  286. // fcount += nQuad
  287. return ;
  288. } // ----- end of method QWEKey::getEyKernel -----
  289. void QWEKey::TestPrivate(const int& N) {
  290. //GaussQuadWeights(N);
  291. //std::cout << "abscissa\n" << GaussAbscissa << std::endl;
  292. //std::cout << "weights\n" << GaussWeights << std::endl;
  293. BesselWeights( J1 );
  294. //BesselZeros(0, N);
  295. std::cout << std::scientific;
  296. std::cout << "BJ0" << BJ0 << std::endl;
  297. std::cout << "BJ1" << BJ1 << std::endl;
  298. //std::cout << "Bess Zero\n" << xInt << std::endl;
  299. }
  300. } // ----- end of Lemma name -----