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.

FHT.h 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324
  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. // Forward declarations
  135. /*
  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. */
  151. template < HANKELTRANSFORMTYPE Type >
  152. Complex FHT<Type>::Zgauss(const int& ii, const Lemma::EMMODE& mode, const int& jj, const Real& val,
  153. const Real& val2, Lemma::KernelEM1DBase* Kernel){
  154. // TODO, in 101 or 51 we never reach here!!
  155. //std::cout << "Zgauss " << std::endl;
  156. return this->Zans(0, Kernel->GetManagerIndex());
  157. }
  158. // Specialisations
  159. // Note that ANDERSON801, CHAVE, QWEKEY will throw errors as they are not consistent
  160. // part of this class
  161. template < HANKELTRANSFORMTYPE Type >
  162. Real FHT< Type >::GetABSER() {
  163. return WT(0,0)/WT(1,0);
  164. }
  165. /* specializations could provide slighly better performance by reducing divides */
  166. // template < >
  167. // Real FHT< FHTKEY201 >::GetABSER() {
  168. // return WT(0,0)/WT(1,0);
  169. // }
  170. //
  171. // template < >
  172. // Real FHT< FHTKEY101 >::GetABSER() {
  173. // return WT(0,0)/WT(1,0);
  174. // }
  175. //
  176. // template < >
  177. // Real FHT< FHTKEY51 >::GetABSER() {
  178. // return WT(0,0)/WT(1,0);
  179. // }
  180. //--------------------------------------------------------------------------------------
  181. // Class: FHT
  182. // Method: ComputeRelated
  183. //--------------------------------------------------------------------------------------
  184. template < HANKELTRANSFORMTYPE Type >
  185. void FHT<Type>::ComputeRelated ( const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager ) {
  186. int nrel = (int)(KernelManager->GetSTLVector().size());
  187. Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
  188. Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(1, nrel);
  189. Zwork.resize(WT.rows(), nrel);
  190. VectorXr lambda = WT.col(0).array()/rho;
  191. int NumFun = 0;
  192. int idx = 0;
  193. // Get Kernel values
  194. for (int ir=0; ir<lambda.size(); ++ir) {
  195. // irelated loop
  196. ++NumFun;
  197. KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
  198. for (int ir2=0; ir2<nrel; ++ir2) {
  199. // Zwork* needed due to sign convention of filter weights
  200. Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
  201. }
  202. }
  203. for (int ir2=0; ir2<nrel; ++ir2) {
  204. Zans(0, ir2) = Zwork.col(ir2).dot(WT.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
  205. }
  206. return ;
  207. } // ----- end of method FHT::ComputeRelated -----
  208. //--------------------------------------------------------------------------------------
  209. // Class: FHT
  210. // Method: ComputeLaggedRelated
  211. //--------------------------------------------------------------------------------------
  212. template < HANKELTRANSFORMTYPE Type >
  213. void FHT<Type>::ComputeLaggedRelated ( const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager ) {
  214. int nrel = (int)(KernelManager->GetSTLVector().size());
  215. Eigen::Matrix< Complex, Eigen::Dynamic, Eigen::Dynamic > Zwork;
  216. Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nlag, nrel);
  217. Zwork.resize(WT.rows()+nlag, nrel); // Zwork needs to be expanded to filter length + nlag
  218. // lambda needs to be expanded to include lagged results
  219. VectorXr lambda = (VectorXr(WT.rows()+nlag) << WT.col(0).array()/rho, VectorXr::Zero(nlag)).finished();
  220. for (int ilam =WT.rows(); ilam< nlag+WT.rows(); ++ilam) {
  221. lambda(ilam) = lambda(ilam-1)/GetABSER();
  222. }
  223. int NumFun = 0;
  224. int idx = 0;
  225. VectorXr Arg(nlag);
  226. Arg(nlag-1) = rho;
  227. for (int ilag=nlag-2; ilag>=0; --ilag) {
  228. Arg(ilag) = Arg(ilag+1) * GetABSER();
  229. }
  230. // Get Kernel values
  231. for (int ir=0; ir<lambda.size(); ++ir) {
  232. // irelated loop
  233. ++NumFun;
  234. KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
  235. for (int ir2=0; ir2<nrel; ++ir2) {
  236. Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
  237. }
  238. }
  239. // Inner product and scale
  240. int ilagr = nlag-1; // Zwork is in opposite order from Arg
  241. for (int ilag=0; ilag<nlag; ++ilag) {
  242. for (int ir2=0; ir2<nrel; ++ir2) {
  243. Zans(ilagr, ir2) = Zwork.col(ir2).segment(ilag,WT.rows()).dot( WT.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder()+1) ) / Arg(ilagr);
  244. }
  245. ilagr -= 1;
  246. }
  247. // make sure vectors are empty
  248. splineVecReal.clear();
  249. splineVecImag.clear();
  250. // Now do cubic spline
  251. for (int ii=0; ii<Zans.cols(); ++ii) {
  252. auto SplineR = CubicSplineInterpolator::NewSP();
  253. SplineR->SetKnots( Arg, Zans.col(ii).real() );
  254. splineVecReal.push_back(SplineR);
  255. auto SplineI = CubicSplineInterpolator::NewSP();
  256. SplineI->SetKnots( Arg, Zans.col(ii).imag() );
  257. splineVecImag.push_back(SplineI);
  258. }
  259. return ;
  260. } // ----- end of method FHT::ComputeLaggedRelated -----
  261. } // ----- end of namespace Lemma ----
  262. #endif // ----- #ifndef FHT_INC -----
  263. /* vim: set tabstop=4 expandtab: */
  264. /* vim: set filetype=cpp: */