.. _program_listing_file_include_gwmodelpp_GTDR.h:

Program Listing for File GTDR.h
===============================

|exhale_lsh| :ref:`Return to documentation for file <file_include_gwmodelpp_GTDR.h>` (``include/gwmodelpp/GTDR.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #ifndef GTDR_H
   #define GTDR_H
   
   #include <vector>
   #include "armadillo_config.h"
   #include <gsl/gsl_vector.h>
   #include "SpatialAlgorithm.h"
   #include "spatialweight/SpatialWeight.h"
   #include "IRegressionAnalysis.h"
   #include "VariableForwardSelector.h"
   #include "IParallelizable.h"
   #include "IBandwidthSelectable.h"
   
   namespace gwm
   {
   
   class GTDR : public SpatialAlgorithm, public IRegressionAnalysis, public IVarialbeSelectable, public IParallelizable, public IParallelOpenmpEnabled
   {
   public:
       typedef arma::mat (GTDR::*PredictCalculator)(const arma::mat&, const arma::mat&, const arma::vec&); 
   
       typedef arma::mat (GTDR::*FitCalculator)(const arma::mat&, const arma::vec&, arma::mat&, arma::vec&, arma::vec&, arma::mat&);   
   
       enum BandwidthCriterionType
       {
           CV,     
           AIC     
       };
   
       typedef double (GTDR::*BandwidthCriterionCalculator)(const std::vector<BandwidthWeight*>&); 
   
       typedef double (GTDR::*IndepVarCriterionCalculator)(const std::vector<std::size_t>&); 
   
   public:
       static RegressionDiagnostic CalcDiagnostic(const arma::mat& x, const arma::vec& y, const arma::mat& betas, const arma::vec& shat);
   
       static arma::vec Fitted(const arma::mat& x, const arma::mat& betas)
       {
           return sum(betas % x, 1);
       }
   
       static double RSS(const arma::mat& x, const arma::mat& y, const arma::mat& betas)
       {
           arma::vec r = y - Fitted(x, betas);
           return sum(r % r);
       }
   
       static double AICc(const arma::mat& x, const arma::mat& y, const arma::mat& betas, const arma::vec& shat)
       {
           double ss = RSS(x, y, betas), n = (double)x.n_rows;
           return n * log(ss / n) + n * log(2 * arma::datum::pi) + n * ((n + shat(0)) / (n - 2 - shat(0)));
       }
   
   public:
   
       GTDR() {}
   
       GTDR(const arma::mat& x, const arma::vec& y, const arma::mat& coords, const std::vector<SpatialWeight>& spatialWeights, bool hasHatMatrix = true, bool hasIntercept = true) : SpatialAlgorithm(coords),
           mX(x),
           mY(y),
           mSpatialWeights(spatialWeights),
           mHasHatMatrix(hasHatMatrix),
           mHasIntercept(hasIntercept)
       {
       }
   
       virtual ~GTDR() {}
   
   public:
   
       const arma::mat& betas() const { return mBetas; }
   
       bool hasHatMatrix() const { return mHasHatMatrix; }
   
       void setHasHatMatrix(bool flag) { mHasHatMatrix = flag; }
   
       const std::vector<SpatialWeight>& spatialWeights() const { return mSpatialWeights; }
   
       void setSpatialWeights(const std::vector<SpatialWeight>& spatialWeights) { mSpatialWeights = spatialWeights; }
   
       bool enableBandwidthOptimize() { return mEnableBandwidthOptimize; }
   
       void setEnableBandwidthOptimize(bool flag) { mEnableBandwidthOptimize = flag; }
   
       double bandwidthOptimizeEps() const { return mBandwidthOptimizeEps; }
   
       void setBandwidthOptimizeEps(double value) { mBandwidthOptimizeEps = value; }
   
       std::size_t bandwidthOptimizeMaxIter() const { return mBandwidthOptimizeMaxIter; }
   
       void setBandwidthOptimizeMaxIter(std::size_t value) { mBandwidthOptimizeMaxIter = value; }
   
       double bandwidthOptimizeStep() const { return mBandwidthOptimizeStep; }
   
       void setBandwidthOptimizeStep(double value) { mBandwidthOptimizeStep = value; }
   
       BandwidthCriterionType bandwidthCriterionType() const { return mBandwidthCriterionType; }
   
       void setBandwidthCriterionType(const BandwidthCriterionType& type);
   
       bool enableIndpenVarSelect() const { return mEnableIndepVarSelect; }
   
       void setEnableIndepVarSelect(bool flag) { mEnableIndepVarSelect = flag; }
   
       double indepVarSelectThreshold() const { return mIndepVarSelectThreshold; }
   
       void setIndepVarSelectThreshold(double threshold) { mIndepVarSelectThreshold = threshold; }
   
       const VariablesCriterionList& indepVarCriterionList() const { return mIndepVarCriterionList; }
   
       const std::vector<std::size_t>& selectedIndepVars() const { return mSelectedIndepVars; }
   
       arma::mat betasSE() { return mBetasSE; }
   
       arma::vec sHat() { return mSHat; }
   
       arma::vec qDiag() { return mQDiag; }
   
       arma::mat s() { return mS; }
   
   public: // Algorithm
       bool isValid() override;
   
   public: // IRegressionAnalysis
       virtual const arma::vec& dependentVariable() const override { return mY; }
       virtual void setDependentVariable(const arma::vec& y) override { mY = y; }
   
       virtual const arma::mat& independentVariables() const override { return mX; }
       virtual void setIndependentVariables(const arma::mat& x) override { mX = x; }
   
       virtual bool hasIntercept() const override { return mHasIntercept; }
       virtual void setHasIntercept(const bool has) override { mHasIntercept = has; }
   
       virtual RegressionDiagnostic diagnostic() const override { return mDiagnostic; }
   
       virtual arma::mat predict(const arma::mat& locations) override { return arma::mat(locations.n_rows, mX.n_cols, arma::fill::zeros); }
   
       virtual arma::mat fit() override;
   
   public:  // IVariableSelectable
       Status getCriterion(const std::vector<std::size_t>& variables, double& criterion) override
       {
           criterion = (this->*mIndepVarCriterionFunction)(variables);
           return mStatus;
       }
   
       std::vector<std::size_t> selectedVariables() override
       {
           return mSelectedIndepVars;
       }
   
   public:  // IParallelOpenmpEnabled
       int parallelAbility() const override
       {
           return ParallelType::SerialOnly
           #ifdef ENABLE_OPENMP
           | ParallelType::OpenMP
           #endif        
           ;
       }
       
       ParallelType parallelType() const override
       {
           return mParallelType;
       }
       
       void setParallelType(const ParallelType& type) override;
   
       void setOmpThreadNum(const int threadNum) override
       {
           mOmpThreadNum = threadNum;
       }
   
   
   public:
   
       double bandwidthCriterion(const std::vector<BandwidthWeight*>& bandwidths)
       {
           return (this->*mBandwidthCriterionFunction)(bandwidths);
       }
   
   protected:
   
       arma::mat predictSerial(const arma::mat& locations, const arma::mat& x, const arma::vec& y);
   
       arma::mat fitSerial(const arma::mat& x, const arma::vec& y, arma::mat& betasSE, arma::vec& shat, arma::vec& qdiag, arma::mat& S);
   
       double bandwidthCriterionAICSerial(const std::vector<BandwidthWeight*>& bandwidths);
   
       double bandwidthCriterionCVSerial(const std::vector<BandwidthWeight*>& bandwidths);
   
       double indepVarCriterionSerial(const std::vector<std::size_t>& indepVars);
   
   #ifdef ENABLE_OPENMP
       arma::mat predictOmp(const arma::mat& locations, const arma::mat& x, const arma::vec& y);
   
       arma::mat fitOmp(const arma::mat& x, const arma::vec& y, arma::mat& betasSE, arma::vec& shat, arma::vec& qdiag, arma::mat& S);
   
       double bandwidthCriterionAICOmp(const std::vector<BandwidthWeight*>& bandwidths);
       
       double bandwidthCriterionCVOmp(const std::vector<BandwidthWeight*>& bandwidths);
       
       double indepVarCriterionOmp(const std::vector<std::size_t>& indepVars);
   #endif
   
   private:
   
       bool isStoreS()
       {
           return mHasHatMatrix && (mCoords.n_rows < 8192);
       }
   
   private:
   
       arma::mat mX;                       
       arma::vec mY;                       
       std::vector<SpatialWeight> mSpatialWeights; 
       arma::mat mBetas;                   
       bool mHasHatMatrix = true;          
       bool mHasIntercept = true;          
       RegressionDiagnostic mDiagnostic = RegressionDiagnostic();   
   
       PredictCalculator mPredictFunction = &GTDR::predictSerial;  
       FitCalculator mFitFunction = &GTDR::fitSerial;              
   
       bool mEnableBandwidthOptimize = false;  
       BandwidthCriterionType mBandwidthCriterionType = BandwidthCriterionType::CV;    
       BandwidthCriterionCalculator mBandwidthCriterionFunction = &GTDR::bandwidthCriterionCVSerial;   
       double mBandwidthOptimizeEps = 1e-6;    
       std::size_t mBandwidthOptimizeMaxIter = 100000; 
       double mBandwidthOptimizeStep = 0.01;   
       double mBandwidthLastCriterion = DBL_MAX;   
   
       bool mEnableIndepVarSelect = false;     
       double mIndepVarSelectThreshold = 3.0;  
       VariablesCriterionList mIndepVarCriterionList;  
       IndepVarCriterionCalculator mIndepVarCriterionFunction = &GTDR::indepVarCriterionSerial;    
       std::vector<std::size_t> mSelectedIndepVars;    
       std::size_t mIndepVarSelectionProgressTotal = 0; 
       std::size_t mIndepVarSelectionProgressCurrent = 0; 
   
       ParallelType mParallelType = ParallelType::SerialOnly;  
       int mOmpThreadNum = 8;  
   
       arma::mat mBetasSE; 
       arma::vec mSHat;    
       arma::vec mQDiag;   
       arma::mat mS;       
   };
   
   
   class GTDRBandwidthOptimizer
   {
   public:
   
       struct Parameter
       {
           GTDR* instance;     
           std::vector<BandwidthWeight*>* bandwidths;  
           arma::uword featureCount;   
       };
   
       static double criterion_function(const gsl_vector* bws, void* params);
   
       static std::string infoBandwidthCriterion(const std::vector<BandwidthWeight*>& weights)
       {
           std::size_t number = 1;
           std::vector<std::string> labels(weights.size());
           std::transform(weights.cbegin(), weights.cend(), labels.begin(), [&number](const BandwidthWeight* bw)
           {
               return std::to_string(number++) + ":" + (bw->adaptive() ? "adaptive" : "fixed");
           });
           return std::string(GWM_LOG_TAG_BANDWIDTH_CIRTERION) + strjoin(",", labels) + ",criterion";
       }
   
       static std::string infoBandwidthCriterion(const std::vector<BandwidthWeight*>& weights, const double criterion)
       {
           std::vector<std::string> labels(weights.size());
           std::transform(weights.cbegin(), weights.cend(), labels.begin(), [](const BandwidthWeight* bw)
           {
               return std::to_string(bw->bandwidth());
           });
           return std::string(GWM_LOG_TAG_BANDWIDTH_CIRTERION) + strjoin(",", labels) + "," + std::to_string(criterion);
       }
   
   public:
   
       explicit GTDRBandwidthOptimizer(const std::vector<BandwidthWeight*>& weights) : mBandwidths(weights) {}
   
       const int optimize(GTDR* instance, arma::uword featureCount, std::size_t maxIter, double eps, double step);
   
   private:
       std::vector<BandwidthWeight*> mBandwidths;  
   };
   
   }
   
   #endif  // GTDR_H