Lemma is an Electromagnetics API
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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327
  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 05/02/2018 09:46:38 PM
  11. * @author Trevor Irons (ti)
  12. * @email Trevor.Irons@utah.edu
  13. * @copyright Copyright (c) 2018, University of Utah
  14. * @copyright Copyright (c) 2018, Trevor Irons & Lemma Software, LLC
  15. */
  16. #ifndef FHT_INC
  17. #define FHT_INC
  18. #pragma once
  19. #include "HankelTransform.h"
  20. #include "CubicSplineInterpolator.h"
  21. namespace Lemma {
  22. /**
  23. \ingroup FDEM1D
  24. \brief Impliments lagged and related fast Hankel transform through
  25. digital filtering.
  26. \details A general Fast Hankel Transform routine which uses the digital
  27. filter apporach. Both lagged and related kernels are supported in
  28. order to minimize kernel function calls.
  29. This approach performs a complete sweep of the
  30. coefficients , for a variant that uses a longer filter which may
  31. be truncated, see FHTAnderson801.
  32. @see FHTAnderson801
  33. @see GQChave
  34. @see QWEKey
  35. */
  36. template < HANKELTRANSFORMTYPE Type >
  37. class FHT : public HankelTransform {
  38. friend std::ostream &operator<<(std::ostream &stream, const FHT<Type> &ob) {
  39. stream << ob.Serialize(); // << "\n";
  40. return stream;
  41. }
  42. public:
  43. // ==================== LIFECYCLE =======================
  44. /**
  45. * Default protected constructor, use NewSP methods to construct
  46. * @see FHT::NewSP
  47. */
  48. explicit FHT (const ctor_key& key ) : HankelTransform( key ) {
  49. }
  50. /**
  51. * Protected DeDerializing constructor, use factory DeSerialize method.
  52. * @see FHT::DeSerialize
  53. */
  54. FHT (const YAML::Node& node, const ctor_key& key) : HankelTransform(node, key) {
  55. }
  56. /** Default protected destructor, use smart pointers (std::shared_ptr) */
  57. ~FHT () {
  58. }
  59. /**
  60. * Factory method for generating concrete class.
  61. * @return a std::shared_ptr of type FHT
  62. */
  63. static std::shared_ptr< FHT > NewSP() {
  64. return std::make_shared< FHT >( ctor_key() );
  65. }
  66. /**
  67. * Uses YAML to serialize this object.
  68. * @return a YAML::Node
  69. * @see FHT::DeSerialize
  70. */
  71. YAML::Node Serialize() const {
  72. YAML::Node node = HankelTransform::Serialize();
  73. node.SetTag( this->GetName() ); // + enum2String(Type) );
  74. //node.SetTag( enum2String(Type) );
  75. //node["var"] = 0;
  76. return node;
  77. }
  78. /**
  79. * Constructs an FHT object from a YAML::Node.
  80. * @see FHT::Serialize
  81. */
  82. static std::shared_ptr<FHT> DeSerialize(const YAML::Node& node);
  83. // ==================== OPERATORS =======================
  84. // ==================== OPERATIONS =======================
  85. Complex Zgauss(const int&, const Lemma::EMMODE&, const int&, const Real&,
  86. const Real&, Lemma::KernelEM1DBase* Kernel);
  87. /// Computes related kernels, if applicable, otherwise this is
  88. /// just a dummy function.
  89. void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DBase> Kernel) {
  90. }
  91. void ComputeRelated(const Real& rho, std::vector< std::shared_ptr<KernelEM1DBase> > KernelVec) {
  92. }
  93. void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager);
  94. void ComputeLaggedRelated(const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager);
  95. // ==================== ACCESS =======================
  96. /**
  97. * @param[in] rho is the argument for lagged convolution evaluation from the
  98. * spline after calculation.
  99. */
  100. void SetLaggedArg(const Real& rho) {
  101. for (int i=0; i<Zans.cols(); ++ i) {
  102. Zans(0, i) = Complex( splineVecReal[i]->Interpolate(rho),
  103. splineVecImag[i]->Interpolate(rho) );
  104. }
  105. return ;
  106. }
  107. // ==================== INQUIRY =======================
  108. /**
  109. * @return filter asbscissa spacing
  110. */
  111. inline Real GetABSER();
  112. //{
  113. // return 0; //this->WT(0,0)/this->WT(1,0);
  114. //}
  115. /** Returns the name of the underlying class, similiar to Python's type */
  116. inline std::string GetName() const {
  117. return enum2String(Type); //this->CName;
  118. }
  119. protected:
  120. // ==================== LIFECYCLE =======================
  121. // ==================== DATA MEMBERS =========================
  122. private:
  123. // Filter Weights, these are specialized for each template type
  124. static const Eigen::Matrix<Real, Eigen::Dynamic, 3> WT;
  125. /// Spines for lagged convolutions (real part)
  126. std::vector <std::shared_ptr<CubicSplineInterpolator> > splineVecReal;
  127. /// Spines for lagged convolutions (imaginary part)
  128. std::vector < std::shared_ptr<CubicSplineInterpolator> > splineVecImag;
  129. /// Holds answer, dimensions are NumConv, and NumberRelated.
  130. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> Zans;
  131. /** ASCII string representation of the class name */
  132. //static constexpr auto CName = "FHT";
  133. }; // ----- end of class FHT ----
  134. // Clang wants forward declarations, MSVC doesn't
  135. #if defined( __clang__) || defined(__GNUC__) || defined(__GNUG__) || defined(__ICC) || defined(__INTEL_COMPILER)
  136. template<>
  137. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<FHTKEY201>::WT;
  138. template<>
  139. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<FHTKEY101>::WT;
  140. template<>
  141. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<FHTKEY51>::WT;
  142. template<>
  143. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<FHTKONG61>::WT;
  144. template<>
  145. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<FHTKONG121>::WT;
  146. template<>
  147. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<FHTKONG241>::WT;
  148. template<>
  149. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT<IRONS>::WT;
  150. // Clang wants generic declaration
  151. template < HANKELTRANSFORMTYPE Type >
  152. const Eigen::Matrix<Real, Eigen::Dynamic, 3> FHT< Type >::WT;
  153. #endif
  154. template < HANKELTRANSFORMTYPE Type >
  155. Complex FHT<Type>::Zgauss(const int& ii, const Lemma::EMMODE& mode, const int& jj, const Real& val,
  156. const Real& val2, Lemma::KernelEM1DBase* Kernel){
  157. // TODO, in 101 or 51 we never reach here!!
  158. //std::cout << "Zgauss " << std::endl;
  159. return this->Zans(0, Kernel->GetManagerIndex());
  160. }
  161. // Specialisations
  162. // Note that ANDERSON801, CHAVE, QWEKEY will throw errors as they are not consistent
  163. // part of this class
  164. template < HANKELTRANSFORMTYPE Type >
  165. Real FHT< Type >::GetABSER() {
  166. return WT(0,0)/WT(1,0);
  167. }
  168. /* specializations could provide slighly better performance by reducing divides */
  169. // template < >
  170. // Real FHT< FHTKEY201 >::GetABSER() {
  171. // return WT(0,0)/WT(1,0);
  172. // }
  173. //
  174. // template < >
  175. // Real FHT< FHTKEY101 >::GetABSER() {
  176. // return WT(0,0)/WT(1,0);
  177. // }
  178. //
  179. // template < >
  180. // Real FHT< FHTKEY51 >::GetABSER() {
  181. // return WT(0,0)/WT(1,0);
  182. // }
  183. //--------------------------------------------------------------------------------------
  184. // Class: FHT
  185. // Method: ComputeRelated
  186. //--------------------------------------------------------------------------------------
  187. template < HANKELTRANSFORMTYPE Type >
  188. void FHT<Type>::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager ) {
  189. int nrel = (int)(KernelManager->GetSTLVector().size());
  190. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
  191. Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  192. Zwork.resize(WT.rows(), nrel);
  193. VectorXr lambda = WT.col(0).array()/rho;
  194. int NumFun = 0;
  195. int idx = 0;
  196. // Get Kernel values
  197. for (int ir=0; ir<lambda.size(); ++ir) {
  198. // irelated loop
  199. ++NumFun;
  200. KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
  201. for (int ir2=0; ir2<nrel; ++ir2) {
  202. // Zwork* needed due to sign convention of filter weights
  203. Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
  204. }
  205. }
  206. for (int ir2=0; ir2<nrel; ++ir2) {
  207. Zans(0, ir2) = Zwork.col(ir2).dot(WT.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
  208. }
  209. return ;
  210. } // ----- end of method FHT::ComputeRelated -----
  211. //--------------------------------------------------------------------------------------
  212. // Class: FHT
  213. // Method: ComputeLaggedRelated
  214. //--------------------------------------------------------------------------------------
  215. template < HANKELTRANSFORMTYPE Type >
  216. void FHT<Type>::ComputeLaggedRelated ( const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager ) {
  217. int nrel = (int)(KernelManager->GetSTLVector().size());
  218. Eigen::Matrix< Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
  219. Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nlag, nrel);
  220. Zwork.resize(WT.rows()+nlag, nrel); // Zwork needs to be expanded to filter length + nlag
  221. // lambda needs to be expanded to include lagged results
  222. VectorXr lambda = (VectorXr(WT.rows()+nlag) << WT.col(0).array()/rho, VectorXr::Zero(nlag)).finished();
  223. for (Index ilam = WT.rows(); ilam< nlag+WT.rows(); ++ilam) {
  224. lambda(ilam) = lambda(ilam-1)/GetABSER();
  225. }
  226. int NumFun = 0;
  227. int idx = 0;
  228. VectorXr Arg(nlag);
  229. Arg(nlag-1) = rho;
  230. for (int ilag=nlag-2; ilag>=0; --ilag) {
  231. Arg(ilag) = Arg(ilag+1) * GetABSER();
  232. }
  233. // Get Kernel values
  234. for (int ir=0; ir<lambda.size(); ++ir) {
  235. // irelated loop
  236. ++NumFun;
  237. KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
  238. for (int ir2=0; ir2<nrel; ++ir2) {
  239. Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
  240. }
  241. }
  242. // Inner product and scale
  243. int ilagr = nlag-1; // Zwork is in opposite order from Arg
  244. for (int ilag=0; ilag<nlag; ++ilag) {
  245. for (int ir2=0; ir2<nrel; ++ir2) {
  246. Zans(ilagr, ir2) = Zwork.col(ir2).segment(ilag,WT.rows()).dot( WT.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder()+1) ) / Arg(ilagr);
  247. }
  248. ilagr -= 1;
  249. }
  250. // make sure vectors are empty
  251. splineVecReal.clear();
  252. splineVecImag.clear();
  253. // Now do cubic spline
  254. for (int ii=0; ii<Zans.cols(); ++ii) {
  255. auto SplineR = CubicSplineInterpolator::NewSP();
  256. SplineR->SetKnots( Arg, Zans.col(ii).real() );
  257. splineVecReal.push_back(SplineR);
  258. auto SplineI = CubicSplineInterpolator::NewSP();
  259. SplineI->SetKnots( Arg, Zans.col(ii).imag() );
  260. splineVecImag.push_back(SplineI);
  261. }
  262. return ;
  263. } // ----- end of method FHT::ComputeLaggedRelated -----
  264. } // ----- end of namespace Lemma ----
  265. #endif // ----- #ifndef FHT_INC -----
  266. /* vim: set tabstop=4 expandtab: */
  267. /* vim: set filetype=cpp: */