123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567 |
- /* 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 -----
|