Lemma/Modules/FDEM1D/include/DigitalFilterIntegratorAnde...

/* This file is part of Lemma, a geophysical modelling and inversion API */

/* 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
  @author   Trevor Irons
  @date     02/07/2011
  @version  $Id: digitalfilterintegrator.h 124 2014-02-07 04:26:26Z tirons $
 **/

#pragma once

#ifndef  DIGITALFILTERINTEGRATOR_INC
#define  DIGITALFILTERINTEGRATOR_INC

#include "IntegrationKernel.h"

namespace Lemma {

    // ===================================================================
    //  Class:  DigitalFilterIntegratorAnderson
    /**
      \brief   Reimplimentation of Walt Anderson's digital filtering
               algorithms which are public domain.
      \details Walt Anderson wrote several routines for digital filtering
               that reused a lot of boilerplate code. The original fortran
               codes are available at:
        ftp://ftpext.usgs.gov/pub/cr/co/denver/musette/pub/anderson/
               This is an abstract class. Concrete classes must define the
               a function to make arguments.
               Template specilizations are found in the implimentation file
               digitalfilterintegrator.cpp. This was necessary to support
               the design paradigm of Lemma.
     */
    // ===================================================================
    template <typename Scalar>
    class DigitalFilterIntegratorAnderson : public LemmaObject {

        /** Prints out basic info about the class, Complex implimentation */

        template <typename Scalar2>
        friend std::ostream &operator<<(std::ostream &stream,
                        const DigitalFilterIntegratorAnderson<Scalar2> &ob);

        public:

            // ====================  LIFECYCLE     =======================

            /** Default protected constructor. */
            explicit DigitalFilterIntegratorAnderson (const ctor_key& key);

            /** Default protected constructor. */
            ~DigitalFilterIntegratorAnderson ();

            /**
             *  @return a YAML::Node serial object
             */
            YAML::Node Serialize() const;

            /** Returns the name of the underlying class, similiar to Python's type */
            virtual std::string GetName() const {
                return this->CName;
            }

            // ====================  OPERATORS     =======================

            // ====================  OPERATIONS    =======================

            /** Sets the number of desired convolutions. From Anderson, where
                this value is called NB:
                NUMBER OF LAGGED CONVOLUTIONS DESIRED (NB.GE.1).  USE
                NB=1 IF B=BMIN=BMAX (I.E.,NO LAGGING DESIRED).  USE
                NB>1 IF B IS LAGGED IN (BMIN,BMAX), WHERE
                BMIN=BMAX*DEXP(-.1D0*(NB-1)) DOES NOT UNDERFLOW THE DEXP
                RANGE.  THE B-LAGGED SPACING IS .1D0 IN LOG-SPACE. FOR
                CONVENIENCE IN SPLINE INTERPOLATION LATER, EACH B IN
                (BMIN,BMAX) IS RETURNED IN ARRAY ARG(I),I=1,NB, WHERE
                ARG(I+1)/ARG(I)=DEXP(.1D0) FOR ALL I.  IF BMAX>BMIN>0 IS
                GIVEN, THEN AN EFFECTIVE VALUE OF NB IS DETERMINED AS
                NB=DINT(10.*DLOG(BMAX/BMIN))+I, WHERE I>1 IS RECOMMENDED,
                PARTICULARLY IF USING SUBSEQUENT SPLINE INTERPOLATION FOR
                A DIFFERENT B-SPACING THAN USED IN THE DIGITAL FILTER. IF
                SPLINE INTERPOLATION IS TO BE USED LATER, IT IS GENERALLY
                BEST TO USE DLOG(ARG(I)) INSTEAD OF ARG(I) -VS- DANS(I),
                FOR I=1,NB.
            */
            void SetNumConv(const int& i);

            /** Evaluates related integrals */
            void ComputeRelated(const Real& rho,
                IntegrationKernel<Scalar>* Kernel);

            /** Evaluates the integral
            @param[in]tol  REQUESTED TRUNCATION TOLERANCE AT BOTH FILTER TAILS
              FOR ADAPTIVE CONVOLUTION FOR ALL NB TRANSFORMS.  THE
              TRUNCATION CRITERION IS ESTABLISHED DURING CONVOLUTION IN
              A FIXED ABSCISSA RANGE (USING WEIGHTS 426-461) OF THE
              COSINE FILTER AS THE MAXIMUM ABSOLUTE CONVOLVED PRODUCT
              TIMES TOL.  THE CONVOLUTION SUMMATION IS TERMINATED
              ON EITHER SIDE OF THE FIXED RANGE WHENEVER THE ABSOLUTE
              PRODUCT .LE. THE TRUNCATION CRITERION.  IN GENERAL, A
              DECREASING TOLERANCE WILL PRODUCE HIGHER ACCURACY SINCE
              MORE FILTER WEIGHTS ARE USED (UNLESS EXPONENT UNDERFLOW
              OCCURS IN THE TRANSFORM INPUT FUNCTION EVALUATION).
              ONE MAY SET TOL=0.D0 TO OBTAIN MAXIMUM ACCURACY FOR ALL
              NB DOUBLE-PRECISION COSINE TRANSFORMS IN DANS(NB).
              HOWEVER, THE ACTUAL RELATIVE ERRORS CANNOT BE EXPECTED TO
              BE SMALLER THAN ABOUT .1D-12 REGARDLESS OF THE TOLERANCE
              VALUE USED, SINCE DOUBLE-PRECISION FILTER WEIGHTS AND A
              DOUBLE-PRECISION FUNCTION ARE USED. IN ANY EVENT,
              ONE SHOULD ALWAYS CHOOSE TOL<<DESIRED RELATIVE ERROR.
           ** ACCURACY WARNING **  SOME HIGHLY OSCILLATORY FUNCTIONS
              FUN(G) AND (OR) LIMITING CASES OF B NEAR MACHINE-ZERO
              (OR INFINITY) SHOULD BE AVOIDED, OTHERWISE UNSATISFACTORY
              RESULTS (E.G., RELATIVE & ABSOLUTE ERRORS>>TOL) MAY OCCUR.
              @param[in] ntol  NUMBER OF CONSECUTIVE TIMES THE TRUNCATION CRITERION (TOL)
              IS TO BE MET AT EITHER FILTER TAIL BEFORE FILTER
              TRUNCATION OCCURS.  NTOL=1 SHOULD BE USED FOR INPUT
              FUNCTIONS THAT DO NOT HAVE MANY ZEROS IN (0,INFINITY). FOR
              OSCILLATORY FUNCTIONS WITH MANY ZEROS, NTOL>1 MAY BE USED
              TO INSURE A PREMATURE CUTOFF DOES NOT OCCUR FOR TRUNCATION
              (SEE USE OF ITOL,NTOL,TOL IN THE CODE BELOW).
            @param[in] rho integration argument argument
            */
            void Compute(const Real& rho, const int& ntol, const Real& tol);

            /** Attaches the integration kernel */
            void AttachKernel(IntegrationKernel<Scalar> *Kernel);

            /** Detaches the integration kernel */
            void DetachKernel( );

            // ====================  ACCESS        =======================

            /** Returns the Vector of Abscissa arguments
             */
            VectorXr  GetAbscissaArguments();

            Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic > GetAnswer();

            /** Returns the number of evaluations of kernel */
            int GetNumFun();
            // ====================  INQUIRY       =======================

        protected:

            // ====================  LIFECYCLE     =======================

            // ====================  DATA MEMBERS  =========================

            /** A rewrite of Anderson's Pseudo-subroutine. */
            void StoreRetreive(const int& idx, const int& lag, const Real& y,
                        Scalar& Zsum, const int& jrel, Scalar& C);

            // ====================  OPERATIONS    ==============================
            /** Sets the filter weights as prescribed by Anderson. This is called
                automatically when the SetBesselOrder(int ) is called. This is
                a protected function and cannot be called directly.
             */
            virtual void SetFilterWeights()=0;

            /** Computes the absolute maximum of the scalar */
            Scalar AbsMax(const Scalar& C, const Scalar& Cmax);

            // ====================  DATA MEMBERS  ==============================

            /** The hankel transform wavenumber embedded in the integral */
            Real Lambda;

            /** Number of times kernel was evaluated */
            int NumFun;

            /** Number of lagged convolutions
                must be greater or equal to 1
                It is set automatically in the @see Compute function so
                that  \f$ \rho \exp\left( -.1*(\mathtt{NumConv} -1) \right) \f$
                does not underflow the exponent range
             */
            int NumConv;

            /** Number of related kernels */
            int NumRel;

            /** Used as base for filter abscissa generation */
            Real ABSCISSA;

            /** Also used in abscissa generation \f$ ABSE = \exp{.1} \f$ */
            const Real ABSE;

            /** Also used in abscissa generation \f$ ABSER = 1 / \exp{.1} \f$ */
            const Real ABSER;

            /** Counter for calculated */
            int icount;

            /** Index of special case low side convolution filter weights */
            int ilow;

            /** Index of special case right side convolution filter weights */
            int ihi;

            /** Kernel Calculator */
            IntegrationKernel<Scalar> *IntKernel;

            /** Holds answer, dimensions are NumConv, and NumberRelated. */
            Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Ans;

            /** Filter weight coefficients.
             */
            VectorXr  FilterWeights;

            /** Stored abscissa arguments
             */
            VectorXr  Arg;

            /** Work from Anderson
             */
            Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Work;

            /// Key used internally
            Eigen::Matrix<int, Eigen::Dynamic, 1> Key;

        private:
            /** ASCII string representation of the class name */
            static constexpr auto CName = "DigitalFilterIntegratorAnderson.h";

    }; // -----  end of class  DigitalFilterIntegratorAnderson  -----


    /////////////////////////////////////////
    /////////////////////////////////////////
    //      Templated implimentation       //
    /////////////////////////////////////////
    /////////////////////////////////////////


    template <typename Scalar>
    std::ostream &operator<<(std::ostream &stream, const DigitalFilterIntegratorAnderson<Scalar> &ob) {
        stream << ob.Serialize()  << "\n";
        return stream;
    }

    // ====================  LIFECYCLE     =======================

    template <typename Scalar>
    DigitalFilterIntegratorAnderson<Scalar>::
            DigitalFilterIntegratorAnderson(const ctor_key& key) :
        LemmaObject(key), Lambda(0), NumFun(0), NumConv(0), NumRel(0),
        ABSCISSA(0),
        ABSE(1.10517091807564762),   //   exp(.1)
        ABSER(0.904837418035959573), // 1/exp(.1)
        ilow(0), ihi(0), IntKernel(NULL) {
    }

    template <typename Scalar>
    YAML::Node DigitalFilterIntegratorAnderson<Scalar>::Serialize() const {
        YAML::Node node = LemmaObject::Serialize();
        node.SetTag( GetName() );
        //node["Type"] = enum2String(Type);
        //node["Location"] = Location;
        //node["Phat"] = Phat;
        //node["Freqs"] = Freqs;
        //node["Phase"] = Phase;
        //node["Moment"] = Moment;
        return node;
    }

    template <typename Scalar>
    DigitalFilterIntegratorAnderson<Scalar>::~DigitalFilterIntegratorAnderson( ) {
    }

    template <typename Scalar>
    Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic >
        DigitalFilterIntegratorAnderson<Scalar>::GetAnswer() {
        return Ans;
    }

    // ====================  OPERATIONS    =======================

    template <typename Scalar>
    void DigitalFilterIntegratorAnderson<Scalar>::SetNumConv(const int& i) {
        this->NumConv = i;
    }

    template <typename Scalar>
    void DigitalFilterIntegratorAnderson<Scalar>::
            AttachKernel(IntegrationKernel<Scalar> *ck) {
		if (this->IntKernel == ck) return;
        if (this->IntKernel != NULL) {
			this->IntKernel->DetachFrom(this);
		}
		ck->AttachTo(this);
		this->IntKernel = ck;
	}


    template <typename Scalar>
	void DigitalFilterIntegratorAnderson<Scalar>::DetachKernel( ) {
		if (this->IntKernel != NULL) {
			this->IntKernel->DetachFrom(this);
		}
		this->IntKernel = NULL;
	}

//     template < >
//     inline Complex DigitalFilterIntegratorAnderson<Complex>::AbsMax(const Complex& C,
//         const Complex& Cmax) {
// 		return Complex(std::max(std::abs(real(C)), std::real(Cmax)),
// 					   std::max(std::abs(imag(C)), std::imag(Cmax)) );
//     }
//     template < >
//     Real DigitalFilterIntegratorAnderson<Real>::AbsMax(const Real& C,
//         const Real& Cmax) {
//         return std::max(C, Cmax);
//     }

    template <typename Scalar>
    VectorXr  DigitalFilterIntegratorAnderson<Scalar>::GetAbscissaArguments() {
        return this->Arg;
    }

//     ///////////////////////////////////////////
// 	// Computes the transform
//
//     template < >
//     void DigitalFilterIntegratorAnderson<Real>::Compute(const Real& rho,
//         const int& ntol, const Real& tol) {
//
// 		Real y1 = this->ABSCISSA/rho;
// 		this->Key.setZero();
//
// 		// Check to make sure everything is set
// 		if (rho<=0) {
// 			throw std::runtime_error("In DigitalFilterIntegratorAnderson Argument rho < 0.");
// 		}
//
// 		if (this->NumConv<1) {
// 			throw std::runtime_error("In DigitalFilterIntegratorAnderson NumConv is less than 1.");
// 		}
//
// 		if (this->IntKernel == NULL) {
// 			throw std::runtime_error("In DigitalFilterIntegratorAnderson Unset Kernel Calculator");
// 		}
//
// 		Arg = VectorXr::Zero(this->NumConv);
// 		Real y = std::pow(rho*ABSER, this->NumConv-1);
//
// 		if (y <= 0) {
// 			std::cerr << "Exponent Underflow Error";
// 			throw std::underflow_error("Exponent underflow");
// 		}
//
// 		this->Work.resize(Eigen::NoChange, this->IntKernel->GetNumRel());
//
//         // Main Loop
// 		int itol = 0;
// 		int none = 0;
// 		this->NumFun = 0;
// 		int idx = 0;
// 		int istore = 0;
// 		Real Sum(0.);
// 		Real Cmax(0.);
// 		Real C(0.);
// 		Ans.resize(this->NumConv, this->IntKernel->GetNumRel());
// 		Ans.setZero();
//  		// 1010 Loop
//  		for (int ilag=0; ilag < this->NumConv; ++ilag) {
//   			istore = this->NumConv - 1 - ilag;
//   			if (ilag > 0) y1 *= ABSE;
//   			Arg(istore) = ABSCISSA/y1;
//
//  			// 1000 Loop
//  			for (int irel=0; irel < this->IntKernel->GetNumRel(); ++irel) {
// 				//this->SetBesselOrder(this->Ckernel->GetBesselOrder(irel));
//   				none = 0;
//   				itol = ntol;
//   				Sum = Real(0);
//   				Cmax = Real(0);
//   				y = y1;
// 				// Begin right side convolution at weight 298
// 				// counting from 0
//   				idx = ilow;
//   				y *= ABSE;
// 				// Code Block 20 in Anderson
// 				do {
//  					this->StoreRetreive(idx, ilag, y, Sum, irel, C);
//                     Cmax = AbsMax(C, Cmax);
// 					++idx;
// 					y *= ABSE;
// 				} while (idx <= ihi);
// 				//if (real(Cmax) == 0 && imag(Cmax) == 0) none = 1;
// 				if (Cmax == Real(0)) none = 1;
// 				Cmax *= tol;
// 				// Code Block 30 in Anderson
// 				do {
//  					this->StoreRetreive(idx, ilag, y, Sum, irel, C);
// 					if (std::abs(C) <= Cmax) {
// 						--itol;
// 						if (itol < 0 || idx > FilterWeights.size()-1) break;
// 					} else {
// 						itol = ntol;
// 					}
// 					++idx;
// 					y *= ABSE;
// 				} while (idx < FilterWeights.size());
// 				itol = ntol;
// 				y = y1;
// 				// Code Block 60 in Anderson
// 				idx = ilow-1;
// 				do {
//  					this->StoreRetreive(idx, ilag, y, Sum, irel, C);
// 					if (std::abs(C) <= Cmax &&  none == 0) {
// 						--itol;
// 						if (itol < 0 || idx < 0) break;
// 					} else {
// 						itol = ntol;
// 					}
// 					--idx;
// 					y *= ABSER;
// 				} while (idx>=0);
//  				Ans(istore, irel) = Sum/Arg(istore);
//   			}  // End of 1000 loop
//   		}      // End of 1010 loop
// 	}
//
//
//     template < >
// 	void DigitalFilterIntegratorAnderson<Complex>::Compute(const Real &rho,
//             const int& ntol, const Real &tol) {
//
// 		Real y1 = this->ABSCISSA/rho;
//
// 		this->Key.setZero();
//
// 		// Check to make sure everything is set
// 		if (rho<=0) {
// 			throw std::runtime_error("In DigitalFilterIntegratorAnderson Argument rho < 0.");
// 		}
//
// 		if (this->NumConv<1) {
// 			throw std::runtime_error("In DigitalFilterIntegratorAnderson NumConv is less than 1.");
// 		}
//
// 		if (this->IntKernel == NULL) {
// 			throw std::runtime_error("In DigitalFilterIntegratorAnderson Unset Kernel Calculator");
// 		}
//
// 		Arg = VectorXr::Zero(this->NumConv);
// 		Real y = std::pow(rho*ABSER, this->NumConv-1);
//
// 		if (y <= 0) {
// 			std::cerr << "Exponent Underflow Error";
// 			throw std::underflow_error("Exponent underflow");
// 		}
//
// 		this->Work.resize(Eigen::NoChange, this->IntKernel->GetNumRel());
//
//         // Main Loop
// 		int itol = 0;
// 		int none = 0;
// 		this->NumFun = 0;
// 		int idx = 0;
// 		int istore = 0;
// 		Complex Zsum(0.);
// 		Complex Cmax(0.);
// 		Complex C(0.);
// 		Ans.resize(this->NumConv, this->IntKernel->GetNumRel());
// 		Ans.setZero();
//  		// 1010 Loop
//  		for (int ilag=0; ilag < this->NumConv; ++ilag) {
//   			istore = this->NumConv - 1 - ilag;
//   			if (ilag > 0) y1 *= ABSE;
//   			Arg(istore) = ABSCISSA/y1;
//  			// 1000 Loop
//  			for (int irel=0; irel < this->IntKernel->GetNumRel(); ++irel) {
// 				//this->SetBesselOrder(this->Ckernel->GetBesselOrder(irel));
//   				none = 0;
//   				itol = ntol;
//   				Zsum = Complex(0);
//   				Cmax = Complex(0);
//   				y = y1;
// 				// Begin right side convolution at weight 298
// 				// counting from 0
//   				idx = ilow;
//   				y *= ABSE;
// 				// Code Block 20 in Anderson
// 				do {
//  					this->StoreRetreive(idx, ilag, y, Zsum, irel, C);
//                     Cmax = AbsMax(C, Cmax);
// 					++idx;
// 					y *= ABSE;
// 				} while (idx <= ihi);
// 				//if (real(Cmax) == 0 && imag(Cmax) == 0) none = 1;
// 				if (Cmax == Complex(0)) none = 1;
// 				Cmax *= tol;
// 				// Code Block 30 in Anderson
// 				do {
//  					this->StoreRetreive(idx, ilag, y, Zsum, irel, C);
// 					if ( std::abs(real(C)) <= real(Cmax) &&
// 						 std::abs(imag(C)) <= imag(Cmax) ) {
// 						--itol;
// 						if (itol < 0 || idx > FilterWeights.size()-1) break;
// 					} else {
// 						itol = ntol;
// 					}
// 					++idx;
// 					y *= ABSE;
// 				} while (idx < FilterWeights.size());
// 				itol = ntol;
// 				y = y1;
// 				// Code Block 60 in Anderson
// 				idx = ilow-1;
// 				do {
//  					this->StoreRetreive(idx, ilag, y, Zsum, irel, C);
// 					if ( std::abs(real(C)) <= real(Cmax) &&
// 						 std::abs(imag(C)) <= imag(Cmax) &&
// 						 none == 0 ) {
// 						--itol;
// 						if (itol < 0 || idx < 0) break;
// 					} else {
// 						itol = ntol;
// 					}
// 					--idx;
// 					y *= ABSER;
// 				} while (idx>=0);
//  				Ans(istore, irel) = Zsum/Arg(istore);
//   			}  // End of 1000 loop
//   		}      // End of 1010 loop
// 	}

    template <typename Scalar>
    int DigitalFilterIntegratorAnderson<Scalar>::GetNumFun() {
        return NumFun;
    }

    // generic rewrite of store-retreive 'pseudo-subroutine'
    template <typename Scalar>
    void DigitalFilterIntegratorAnderson<Scalar>::StoreRetreive(const int &idx,
            const int& lag, const Real& y, Scalar& Sum,
            const int& jrel, Scalar& C) {

		int look = idx+lag;
 		int	iq = look/FilterWeights.size();
 		int	ir = look%FilterWeights.size();
		int iroll = iq*(FilterWeights.size()-1);

		if(this->Key(ir) <= iroll) {
			this->Key(ir) = iroll + ir;
			this->Lambda = y;
 			++this->NumFun;
 			for (int ir2=0; ir2<IntKernel->GetNumRel(); ++ir2) {
 			    this->Work(ir, ir2) = this->IntKernel->Argument(this->Lambda, ir2);
 			}
		}

		C = this->Work(ir, jrel) * this->FilterWeights(idx);
		Sum += C;
		return;

	}

}		// -----  end of Lemma  name  -----

#endif   // ----- #ifndef DIGITALFILTERINTEGRATOR_INC  -----