Program Listing for File GWRMultiscale.h

↰ Return to documentation for file (include/gwmodelpp/GWRMultiscale.h)

#ifndef GWRMULTISCALE_H
#define GWRMULTISCALE_H

#include <utility>
#include <string>
#include <initializer_list>
#include <optional>
#include "SpatialMultiscaleAlgorithm.h"
#include "spatialweight/SpatialWeight.h"
#include "IRegressionAnalysis.h"
#include "IBandwidthSelectable.h"
#include "BandwidthSelector.h"
#include "IParallelizable.h"


namespace gwm
{

#define GWM_LOG_TAG_MGWR_INITIAL_BW "#initial-bandwidth "
#define GWM_LOG_TAG_MGWR_BACKFITTING "#back-fitting "

#define GWM_LOG_MGWR_BACKFITTING(MESSAGE) { GWM_LOG_INFO((std::string(GWM_LOG_TAG_MGWR_BACKFITTING) + (MESSAGE))); }

class GWRMultiscale : public SpatialMultiscaleAlgorithm,
    public IBandwidthSelectable,
    public IRegressionAnalysis,
    public IParallelizable,
    public IParallelOpenmpEnabled,
    public IParallelCudaEnabled,
    public IParallelMpiEnabled
{
public:

    enum BandwidthInitilizeType
    {
        Null,
        Initial,
        Specified
    };

    static std::unordered_map<BandwidthInitilizeType,std::string> BandwidthInitilizeTypeNameMapper;

    enum BandwidthSelectionCriterionType
    {
        CV,
        AIC
    };

    static std::unordered_map<BandwidthSelectionCriterionType,std::string> BandwidthSelectionCriterionTypeNameMapper;

    enum BackFittingCriterionType
    {
        CVR,
        dCVR
    };

    static std::unordered_map<BackFittingCriterionType,std::string> BackFittingCriterionTypeNameMapper;

    typedef double (GWRMultiscale::*BandwidthSizeCriterionFunction)(BandwidthWeight*);

    typedef arma::vec (GWRMultiscale::*FitVarFunction)(const size_t);

    typedef arma::vec (GWRMultiscale::*FitVarCoreFunction)(const arma::vec&, const arma::vec&, const SpatialWeight&, arma::mat&);

    typedef arma::vec (GWRMultiscale::*FitVarCoreCVFunction)(const arma::vec&, const arma::vec&, const SpatialWeight&);

    typedef arma::vec (GWRMultiscale::*FitVarCoreSHatFunction)(const arma::vec&, const arma::vec&, const SpatialWeight&, arma::vec&);

private:

    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)));
    }

    bool isStoreS()
    {
        return mHasHatMatrix && (mCoords.n_rows < 8192);
    }

    static RegressionDiagnostic CalcDiagnostic(const arma::mat &x, const arma::vec &y, const arma::vec &shat, double RSS);

    void setInitSpatialWeight(const SpatialWeight &spatialWeight)
    {
        mInitSpatialWeight = spatialWeight;
    }

public:

    GWRMultiscale() {}

    GWRMultiscale(const arma::mat& x, const arma::vec& y, const arma::mat& coords, const std::vector<SpatialWeight>& spatialWeights)
        : SpatialMultiscaleAlgorithm(coords, spatialWeights)
    {
        mX = x;
        mY = y;
        if (spatialWeights.size() > 0)
            setInitSpatialWeight(spatialWeights[0]);
    }

    virtual ~GWRMultiscale() {}

public:

    void setGoldenUpperBounds(double value) { mGoldenUpperBounds = value; }

    void setGoldenLowerBounds(double value) { mGoldenLowerBounds = value; }

    const std::vector<BandwidthInitilizeType>& bandwidthInitilize() const { return GWRMultiscale::mBandwidthInitilize; }

    void setBandwidthInitilize(const std::vector<BandwidthInitilizeType> &bandwidthInitilize);

    const std::vector<BandwidthSelectionCriterionType>& bandwidthSelectionApproach() const { return GWRMultiscale::mBandwidthSelectionApproach; }

    void setBandwidthSelectionApproach(const std::vector<BandwidthSelectionCriterionType> &bandwidthSelectionApproach);

    const std::vector<bool>& preditorCentered() const { return mPreditorCentered; }

    void setPreditorCentered(const std::vector<bool> &preditorCentered) { mPreditorCentered = preditorCentered; }

    const std::vector<double>& bandwidthSelectThreshold() const { return mBandwidthSelectThreshold; }

    void setBandwidthSelectThreshold(const std::vector<double> &bandwidthSelectThreshold) { mBandwidthSelectThreshold = bandwidthSelectThreshold; }

    bool hasHatMatrix() const { return mHasHatMatrix; }

    void setHasHatMatrix(bool hasHatMatrix) { mHasHatMatrix = hasHatMatrix; }

    size_t bandwidthSelectRetryTimes() const { return (size_t)mBandwidthSelectRetryTimes; }

    void setBandwidthSelectRetryTimes(size_t bandwidthSelectRetryTimes) { mBandwidthSelectRetryTimes = (arma::uword)bandwidthSelectRetryTimes; }

    size_t maxIteration() const { return mMaxIteration; }

    void setMaxIteration(size_t maxIteration) { mMaxIteration = maxIteration; }

    BackFittingCriterionType criterionType() const { return mCriterionType; }

    void setCriterionType(const BackFittingCriterionType &criterionType) { mCriterionType = criterionType; }

    double criterionThreshold() const { return mCriterionThreshold; }

    void setCriterionThreshold(double criterionThreshold) { mCriterionThreshold = criterionThreshold; }

    int adaptiveLower() const { return mAdaptiveLower; }

    void setAdaptiveLower(int adaptiveLower) { mAdaptiveLower = adaptiveLower; }

    const arma::mat& betas() const { return mBetas; }

    BandwidthSizeCriterionFunction bandwidthSizeCriterionVar(BandwidthSelectionCriterionType type);

public:     // SpatialAlgorithm interface
    bool isValid() override;


public:     // SpatialMultiscaleAlgorithm interface
    virtual void setSpatialWeights(const std::vector<SpatialWeight> &spatialWeights) override;


public:     // IBandwidthSizeSelectable interface
    Status getCriterion(BandwidthWeight* weight, double& criterion) override
    {
        criterion = (this->*mBandwidthSizeCriterion)(weight);
        return mStatus;
    }


public:     // IRegressionAnalysis interface
    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; }

    arma::mat predict(const arma::mat& locations) override { return arma::mat(locations.n_rows, mX.n_cols, arma::fill::zeros); }

    arma::mat fit() override;

public:     // IParallelalbe interface
    int parallelAbility() const override
    {
        return ParallelType::SerialOnly
#ifdef ENABLE_OPENMP
            | ParallelType::OpenMP
#endif
#ifdef ENABLE_CUDA
            | ParallelType::CUDA
#endif
            ;
    }

    ParallelType parallelType() const override { return mParallelType; }

    void setParallelType(const ParallelType &type) override;


public:     // IOpenmpParallelable interface
    void setOmpThreadNum(const int threadNum) override { mOmpThreadNum = threadNum; }

public:     // IParallelCudaEnabled interface
    void setGPUId(const int gpuId) override { mGpuId = gpuId; }
    void setGroupSize(const std::size_t size) override { mGroupLength = size; }
    int workerId() override { return mWorkerId; }
    void setWorkerId(int id) override { mWorkerId = id; };
    void setWorkerNum(int size) override { mWorkerNum = size; };

protected:

    BandwidthWeight* bandwidth(size_t i)
    {
        return mSpatialWeights[i].weight<BandwidthWeight>();
    }

    arma::mat fitInitial();

    arma::mat backfitting(const arma::mat& betas0);

    arma::vec fitVarBase(const size_t var);

    double bandwidthSizeCriterionVarCVBase(BandwidthWeight* bandwidthWeight);

    double bandwidthSizeCriterionVarAICBase(BandwidthWeight* bandwidthWeight);

    arma::vec fitVarCoreSerial(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw, arma::mat& S);

    arma::vec fitVarCoreCVSerial(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw);

    arma::vec fitVarCoreSHatSerial(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw, arma::vec& shat);

#ifdef ENABLE_OPENMP
    arma::vec fitVarCoreOmp(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw, arma::mat& S);

    arma::vec fitVarCoreCVOmp(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw);

    arma::vec fitVarCoreSHatOmp(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw, arma::vec& shat);
#endif

#ifdef ENABLE_CUDA

    arma::vec fitVarCoreCuda(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw, arma::mat& S);

    arma::vec fitVarCoreCVCuda(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw);

    arma::vec fitVarCoreSHatCuda(const arma::vec& x, const arma::vec& y, const SpatialWeight& sw, arma::vec& shat);

#endif // ENABLE_CUDA

#ifdef ENABLE_MPI
    arma::vec fitVarMpi(const size_t var);
    double bandwidthSizeCriterionVarCVMpi(BandwidthWeight* bandwidthWeight);
    double bandwidthSizeCriterionVarAICMpi(BandwidthWeight* bandwidthWeight);
#endif // ENABLE_MPI

private:
    FitVarFunction mFitVar = &GWRMultiscale::fitVarBase;
    FitVarCoreFunction mFitVarCore = &GWRMultiscale::fitVarCoreSerial;
    FitVarCoreCVFunction mFitVarCoreCV = &GWRMultiscale::fitVarCoreCVSerial;
    FitVarCoreSHatFunction mFitVarCoreSHat = &GWRMultiscale::fitVarCoreSHatSerial;

    SpatialWeight mInitSpatialWeight;
    BandwidthSizeCriterionFunction mBandwidthSizeCriterion = &GWRMultiscale::bandwidthSizeCriterionVarCVBase;
    size_t mBandwidthSelectionCurrentIndex = 0;
    double mBandwidthLastCriterion = DBL_MAX;
    std::optional<double> mGoldenUpperBounds;
    std::optional<double> mGoldenLowerBounds;
    std::vector<double> mMaxDistances;
    std::vector<double> mMinDistances;

    std::vector<BandwidthInitilizeType> mBandwidthInitilize;
    std::vector<BandwidthSelectionCriterionType> mBandwidthSelectionApproach;
    std::vector<bool> mPreditorCentered;
    std::vector<double> mBandwidthSelectThreshold;
    arma::uword mBandwidthSelectRetryTimes = 5;
    size_t mMaxIteration = 500;
    BackFittingCriterionType mCriterionType = BackFittingCriterionType::dCVR;
    double mCriterionThreshold = 1e-6;
    int mAdaptiveLower = 10;

    bool mHasHatMatrix = true;

    arma::mat mX;
    arma::vec mY;
    arma::mat mBetas;
    arma::mat mBetasSE;
    arma::mat mBetasTV;
    bool mHasIntercept = true;

    arma::mat mS0;
    arma::cube mSArray;
    arma::cube mC;
    arma::mat mX0;
    arma::vec mY0;
    arma::vec mXi;
    arma::vec mYi;

    double mRSS0;

    RegressionDiagnostic mDiagnostic;

    ParallelType mParallelType = ParallelType::SerialOnly;
    int mOmpThreadNum = 8;
    size_t mGroupLength = 64;
    int mGpuId = 0;
    int mWorkerId = 0;
    int mWorkerNum = 1;
#ifdef ENABLE_MPI
    arma::uword mWorkRangeSize = 0;
#endif // ENABLE_MPI
    std::pair<arma::uword, arma::uword> mWorkRange = std::make_pair(arma::uword(0), arma::uword(0));

public:
    static int treeChildCount;
};

}

#endif // GWRMULTISCALE_H