work towards lagged Key201

Modules/FDEM1D/include/EMEarth1D.h

@@ -22,6 +22,9 @@
 //#include "KernelEM1DManager.h"
 
 #include "KernelEM1DSpec.h"
+
+#include "HankelTransformFactory.h"
+
 #include "GQChave.h"
 #include "FHTAnderson801.h"
 #include "FHTKey201.h"
@@ -174,7 +177,7 @@ namespace Lemma {
             /** Used internally, this is the innermost loop of the MakeCalc3,
              *  and CalculateWireAntennaField routines.
              */
-            void SolveLaggedTxRxPair(const int &irec, FHTAnderson801* Hankel,
+            void SolveLaggedTxRxPair(const int &irec, HankelTransform* Hankel,
                     const Real &wavef, const int &ifreq,
                     PolygonalWireAntenna* antenna);
 

Modules/FDEM1D/include/HankelTransform.h

@@ -69,6 +69,14 @@ namespace Lemma {
                 /** Returns the name of the underlying class, similiar to Python's type */
                 virtual inline std::string GetName() const = 0 ;
 
+
+                virtual Real GetABSER() { return 0; }
+
+                virtual void ComputeLaggedRelated( const Real& rhomax, const int& nlag, std::shared_ptr<KernelEM1DManager> Manager ) {
+                }
+
+                virtual void SetLaggedArg( const Real& rho ) {}
+
                 // ====================  DATA MEMBERS  =======================
 
             protected:

Modules/FDEM1D/include/HankelTransformFactory.h

@@ -0,0 +1,140 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      04/18/2018 11:59:00 AM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2018, University of Utah
+ * @copyright Copyright (c) 2018, Lemma Software, LLC
+ */
+
+
+#pragma once
+#include "LemmaObject.h"
+
+#include "FHTAnderson801.h"
+#include "FHTKey201.h"
+
+namespace Lemma {
+
+    /**
+     * \ingroup FDEM1D
+     * \brief  Factory generator of HankelTranform types
+     * \details
+     */
+    class HankelTransformFactory : public LemmaObject {
+
+        friend std::ostream &operator<<(std::ostream &stream, const HankelTransformFactory &ob);
+
+        protected:
+        /*
+         *  This key is used to lock the constructor. It is protected so that inhereted
+         *  classes also have the key to contruct their base class.
+         */
+        struct ctor_key {};
+
+        public:
+
+        // ====================  LIFECYCLE     =======================
+
+        /**
+         * Default constructor.
+         * @note This method is locked, and cannot be called directly.
+         *       The reason that the method is public is to enable the use
+         *       of make_shared whilst enforcing the use of shared_ptr,
+         *       in c++-17, this curiosity may be resolved.
+         * @see HankelTransformFactory::NewSP
+         */
+        explicit HankelTransformFactory ( const ctor_key& );
+
+        /**
+         * DeSerializing constructor.
+         * @note This method is locked, and cannot be called directly.
+         *       The reason that the method is public is to enable the use
+         *       of make_shared whilst enforcing the use of shared_ptr,
+         *       in c++-17, this curiosity may be resolved.
+         * @see HankelTransformFactory::DeSerialize
+         */
+        HankelTransformFactory ( const YAML::Node& node, const ctor_key& );
+
+        /**
+         * Default destructor.
+         * @note This method should never be called due to the mandated
+         *       use of smart pointers. It is necessary to keep the method
+         *       public in order to allow for the use of the more efficient
+         *       make_shared constructor.
+         */
+        virtual ~HankelTransformFactory ();
+
+        /**
+         *  Uses YAML to serialize this object.
+         *  @return a YAML::Node
+         *  @see HankelTransformFactory::DeSerialize
+         */
+        virtual YAML::Node Serialize() const;
+
+        /*
+         *  Factory method for generating concrete class.
+         *  @return a std::shared_ptr of type HankelTransformFactory
+         */
+        //static std::shared_ptr< HankelTransformFactory > NewSP();
+        static std::shared_ptr< HankelTransform > NewSP( const HANKELTRANSFORMTYPE Type) {
+            switch (Type) {
+                case ANDERSON801:
+                    break;
+            }
+            return FHTKey201::NewSP();
+            //return FHTAnderson801::NewSP();
+        }
+
+        /**
+         *   Constructs an HankelTransformFactory object from a YAML::Node.
+         *   @see HankelTransformFactory::Serialize
+         */
+        static std::shared_ptr<HankelTransformFactory> DeSerialize(const YAML::Node& node);
+
+        // ====================  OPERATORS     =======================
+
+        // ====================  OPERATIONS    =======================
+
+        // ====================  ACCESS        =======================
+
+        // ====================  INQUIRY       =======================
+        /**
+         *  Returns the name of the underlying class, similiar to Python's type
+         *  @return string of class name
+         */
+        virtual inline std::string GetName() const {
+            return CName;
+        }
+
+        protected:
+
+        // ====================  LIFECYCLE     =======================
+
+        /** Copy is disabled */
+        HankelTransformFactory( const HankelTransformFactory& ) = delete;
+
+        // ====================  DATA MEMBERS  =========================
+
+        private:
+
+        /** ASCII string representation of the class name */
+        static constexpr auto CName = "HankelTransformFactory";
+
+    }; // -----  end of class  HankelTransformFactory  -----
+}  // -----  end of namespace Lemma ----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */
+
+

Modules/FDEM1D/src/EMEarth1D.cpp

@@ -231,7 +231,7 @@ namespace Lemma {
         if (Antenna->GetName() == std::string("PolygonalWireAntenna") || Antenna->GetName() == std::string("TEMTransmitter") ) {
             icalc += 1;
             // Check to see if they are all on a plane? If so we can do this fast
-            if (Antenna->IsHorizontallyPlanar() && 1==2 && ( HankelType == ANDERSON801 || HankelType== FHTKEY201  )) {
+            if (Antenna->IsHorizontallyPlanar() && ( HankelType == ANDERSON801 || HankelType== FHTKEY201  )) {
                 #ifdef HAVE_BOOST_PROGRESS
                 if (progressbar) {
                     disp = new boost::progress_display( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
@@ -243,12 +243,7 @@ namespace Lemma {
                     #pragma omp parallel
                     {
                     #endif
-                    //if (HankelType == ANDERSON801) {
-                        auto Hankel = FHTAnderson801::NewSP();
-                    //}
-                    //else if(HankelType == FHTKEY201) {
-                    //    auto Hankel = FHTKey201::NewSP();
-                    //}
+                    auto Hankel = HankelTransformFactory::NewSP( HankelType );
                     #ifdef LEMMAUSEOMP
                     #pragma omp for schedule(static, 1)
                     #endif
@@ -767,7 +762,7 @@ namespace Lemma {
 //         //tDipole->UpdateFields( ifreq,  Hankel, wavef );
 //     }
 
-    void EMEarth1D::SolveLaggedTxRxPair(const int &irec, FHTAnderson801* Hankel,
+    void EMEarth1D::SolveLaggedTxRxPair(const int &irec, HankelTransform* Hankel,
                     const Real &wavef, const int &ifreq, PolygonalWireAntenna* antenna) {
 
         antenna->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
@@ -781,7 +776,6 @@ namespace Lemma {
             rhomin = std::min(rhomin, rho);
             rhomax = std::max(rhomax, rho);
         }
-        //std::cout << "rhomin\t" << rhomin << "\trhomax" << rhomax << std::endl;
 
         // Determine number of lagged convolutions to do
         // TODO, can Hankel2 adjust the lagg spacing safely?
@@ -811,7 +805,6 @@ namespace Lemma {
             Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
             //std::cout << " in Lagged " <<  rho << "\t" << rhomin << "\t" << rhomax << std::endl;
             Hankel->SetLaggedArg( rho );
-            //std::cout << "out Lagged" << std::endl;
             tDipole->UpdateFields( ifreq,  Hankel, wavef );
         }
         //std::cout << "Spline\n";

Modules/FDEM1D/src/FHTKey201.cpp

@@ -376,7 +376,6 @@ namespace Lemma {
     //      Method:  ComputeLaggedRelated
     //--------------------------------------------------------------------------------------
     void FHTKey201::ComputeLaggedRelated ( const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager ) {
-
         int nrel = (int)(KernelManager->GetSTLVector().size());
         Eigen::Matrix<Complex, 201, Eigen::Dynamic > Zwork;
         Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nlag, nrel);
@@ -385,6 +384,12 @@ namespace Lemma {
         int NumFun = 0;
         int idx = 0;
 
+        VectorXr Arg(nlag);
+        Arg(nlag-1) = rho;
+        for (int ilag=nlag-2; ilag>=0; --ilag) {
+            Arg(ilag) = Arg(ilag+1) * GetABSER();
+        }
+
         // Get Kernel values
         for (int ir=0; ir<lambda.size(); ++ir) {
             // irelated loop
@@ -402,8 +407,30 @@ namespace Lemma {
         // in the interests of making them as generic and reusable as possible. This approach requires slightly
         // more multiplies, but the same number of kernel evaluations, which is the expensive part.
         // Inner product and scale
-        for (int ir2=0; ir2<nrel; ++ir2) {
-            Zans(0, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
+
+        for (int ilag=0; ilag<nlag; ++ilag) {
+            for (int ir2=0; ir2<nrel; ++ir2) {
+                //Zans(ilag, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
+                Zans(ilag, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/Arg(ilag);
+            }
+        }
+
+        // make sure vectors are empty
+        splineVecReal.clear();
+        splineVecImag.clear();
+
+        // Now do cubic spline
+        // TODO Check that knots are set in right order, Eigen has reverse()
+        //std::cout << "Arg\n" << Arg << std::endl;
+        //std::cout << "Zans\n" << Zans.col(0) << std::endl;
+        for (int ii=0; ii<Zans.cols(); ++ii) {
+            auto Spline = CubicSplineInterpolator::NewSP();
+            Spline->SetKnots( Arg, Zans.col(ii).real() );
+            splineVecReal.push_back(Spline);
+
+            auto SplineI = CubicSplineInterpolator::NewSP();
+            SplineI->SetKnots( Arg, Zans.col(ii).imag() );
+            splineVecImag.push_back(SplineI);
         }
 
         return ;
@@ -415,7 +442,7 @@ namespace Lemma {
     //      Method:  GetABSER
     //--------------------------------------------------------------------------------------
     Real FHTKey201::GetABSER (  ) {
-        return WT201(1,0)/WT201(0,0);
+        return WT201(0,0)/WT201(1,0);
     }		// -----  end of method FHTKey201::GetABSER  -----
 
 
@@ -431,5 +458,4 @@ namespace Lemma {
         return ;
     }		// -----  end of method FHTKey201::SetLaggedArg  -----
 
-
 }		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/CubicSplineInterpolator.cpp

@@ -178,9 +178,10 @@ namespace Lemma {
         }
         --i;
 
-        //if ( x > Spline.x[i] ) {
-        //    throw std::runtime_error("CubicSplineInterpolator::Interpolate ATTEMPT TO INTERPOLATE PAST LAST KNOT");
-        //}
+//         if ( x > Spline.x[i] ) {
+//             std::cout << "DOOM\t" << x << "\t" << i << "\t" << Spline.x[i];
+//             throw std::runtime_error("CubicSplineInterpolator::Interpolate ATTEMPT TO INTERPOLATE PAST LAST KNOT");
+//         }
 
         return Spline.a[i] + Spline.b[i]*(x-Spline.x[i]) + Spline.c[i]*((x-Spline.x[i])*(x-Spline.x[i])) +
                Spline.d[i]*((x-Spline.x[i])*(x-Spline.x[i])*(x-Spline.x[i]) );