Fix/gwpca

include/gwmodel.h

@@ -40,5 +40,6 @@
 #include "gwmodelpp/GWAverage.h"
 #include "gwmodelpp/GWCorrelation.h"
 #include "gwmodelpp/GWPCA.h"
+#include "gwmodelpp/GTWR.h" 
 
-#endif  // GWMODEL_H
+#endif  // GWMODEL_H

include/gwmodelpp/GWPCA.h

@@ -19,7 +19,7 @@ namespace gwm
 class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis
 {
 private:
-    typedef arma::mat (GWPCA::*Solver)(const arma::mat&, arma::cube&, arma::mat&);  //!< \~english Calculator to solve \~chinese 模型求解函数
+    typedef arma::mat (GWPCA::*Solver)(const arma::mat&, arma::cube&, arma::cube&, arma::mat&);  //!< \~english Calculator to solve \~chinese 模型求解函数
 
 public: // Constructors and Deconstructors
 
@@ -62,6 +62,20 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis
      */
     void setKeepComponents(int k) { mK = k; }
 
+    /**
+     * @brief \~english Get the Robust flag. \~chinese 获取是否使用鲁棒模式。
+     * 
+     * @return bool \~english Robust flag \~chinese 是否使用鲁棒模式
+     */
+    bool robust() { return mRobust; }
+
+    /**
+     * @brief \~english Set the Robust flag. \~chinese 设置是否使用鲁棒模式。
+     * 
+     * @param robust \~english Robust flag \~chinese 是否使用鲁棒模式
+     */
+    void setRobust(bool robust) { mRobust = robust; mSolver = robust ? &GWPCA::solveRobustSerial : &GWPCA::solveSerial; }
+
     /**
      * @brief \~english Get the Local Principle Values matrix. \~chinese 获取局部主成分值。
      * 
@@ -86,7 +100,7 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis
     /**
      * @brief \~english Get the Scores matrix. \~chinese 获取得分矩阵。
      * 
-     * @return arma::mat \~english Scores matrix \~chinese 得分矩阵
+     * @return arma::mat \~english Scores matrix \~chinese 得分矩阵1
      */
     const arma::cube& scores() { return mScores; }
 
@@ -105,52 +119,77 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis
      * 
      * @param x \~english Symmetric data matrix \~chinese 对称数据矩阵
      * @param loadings [out] \~english Out reference to loadings matrix \~chinese 载荷矩阵
+     * @param scores [out] \~english Out reference to scores matrix \~chinese 得分矩阵
      * @param sdev [out] \~english Out reference to standard deviation matrix \~chinese 标准差
      * @return arma::mat \~english Principle values matrix \~chinese 主成分值矩阵
      */
-    arma::mat pca(const arma::mat& x, arma::cube& loadings, arma::mat& sdev)
+    arma::mat pca(const arma::mat& x, arma::cube& loadings, arma::cube& scores, arma::mat& sdev)
     {
-        return (this->*mSolver)(x, loadings, sdev);
+        return (this->*mSolver)(x, loadings, scores, sdev);
     }
 
     /**
-     * @brief \~english Serial version of PCA funtion. \~chinese 单线程 PCA 函数。
+     * @brief \~english Serial version of PCA funtion. \~chinese 单线程 PCA 函数。1
+     * 
+     * @param x \~english Symmetric data matrix \~chinese 对称数据矩阵1
+     * @param loadings [out] \~english Out reference to loadings matrix \~chinese 载荷矩阵1
+     * @param scores [out] \~english Out reference to scores matrix \~chinese 得分矩阵1
+     * @param sdev [out] \~english Out reference to standard deviation matrix \~chinese 标准差1
+     * @return arma::mat \~english Principle values matrix \~chinese 主成分值矩阵1
+     */
+    arma::mat solveSerial(const arma::mat& x, arma::cube& loadings, arma::cube& scores, arma::mat& sdev);
+
+    /**
+     * @brief \~english Robust serial version of PCA function. \~chinese 鲁棒单线程 PCA 函数。
      * 
      * @param x \~english Symmetric data matrix \~chinese 对称数据矩阵
      * @param loadings [out] \~english Out reference to loadings matrix \~chinese 载荷矩阵
+     * @param scores [out] \~english Out reference to scores matrix \~chinese 得分矩阵
      * @param sdev [out] \~english Out reference to standard deviation matrix \~chinese 标准差
      * @return arma::mat \~english Principle values matrix \~chinese 主成分值矩阵
      */
-    arma::mat solveSerial(const arma::mat& x, arma::cube& loadings, arma::mat& sdev);
+    arma::mat solveRobustSerial(const arma::mat& x, arma::cube& loadings, arma::cube& scores, arma::mat& sdev);
 
     /**
      * @brief \~english Function to carry out weighted PCA. \~chinese 执行加权PCA的函数。
      * 
      * @param x \~english Symmetric data matrix \~chinese 对称数据矩阵
      * @param w \~english Weight vector \~chinese 权重向量
-     * @param V [out] \~english Right orthogonal matrix \~chinese 右边的正交矩阵
-     * @param d [out] \~english Rectangular diagonal matri \~chinese 矩形对角阵
+     * @param U [out] \~english Left orthogonal matrix (scores) \~chinese 左正交矩阵（得分）
+     * @param V [out] \~english Right orthogonal matrix (loadings) \~chinese 右正交矩阵（载荷）
+     * @param d [out] \~english Rectangular diagonal matrix \~chinese 矩形对角阵
+     */
+    void wpca(const arma::mat& x, const arma::vec& w, arma::mat& U, arma::mat& V, arma::vec & d);
+
+    /**
+     * @brief \~english Function to carry out robust weighted PCA. \~chinese 执行鲁棒加权PCA的函数。
+     * 
+     * @param x \~english Symmetric data matrix \~chinese 对称数据矩阵
+     * @param w \~english Weight vector \~chinese 权重向量
+     * @param U [out] \~english Left orthogonal matrix (scores) \~chinese 左正交矩阵（得分）
+     * @param V [out] \~english Right orthogonal matrix (loadings) \~chinese 右正交矩阵（载荷）
+     * @param d [out] \~english Rectangular diagonal matrix \~chinese 矩形对角阵
      */
-    void wpca(const arma::mat& x, const arma::vec& w, arma::mat& V, arma::vec & d);
+    void rwpca(const arma::mat& x, const arma::vec& w, arma::mat& U, arma::mat& V, arma::vec & d);
 
 private:    // Algorithm Parameters
     int mK = 2;  //!< \~english Number of components to be kept \~chinese 要保留的主成分数量
-    // bool mRobust = false;
+    bool mRobust = false;  //!< \~english Robust mode flag \~chinese 鲁棒模式标志
 
 private:    // Algorithm Results
-    arma::mat mLocalPV;               //!< \~english Local principle component values \~chinese 局部主成分值
-    arma::cube mLoadings;             //!< \~english Loadings for each component \~chinese 局部载荷矩阵
-    arma::mat mSDev;                  //!< \~english Standard Deviation \~chinese 标准差矩阵
-    arma::cube mScores;               //!< \~english Scores for each variable \~chinese 得分矩阵
-    arma::uvec mWinner;               //!< \~english Winner variable at each sample \~chinese 优胜变量索引值
+    arma::mat mLocalPV;               //!< \~english Local principle component values \~chinese 局部主成分值1
+    arma::cube mLoadings;             //!< \~english Loadings for each component \~chinese 局部载荷矩阵1
+    arma::mat mSDev;                  //!< \~english Standard Deviation \~chinese 标准差矩阵1
+    arma::cube mScores;               //!< \~english Scores for each variable \~chinese 得分矩阵1
+    arma::uvec mWinner;               //!< \~english Winner variable at each sample \~chinese 优胜变量索引值1
 
 private:    // Algorithm Runtime Variables
-    arma::mat mX;           //!< \~english Variable matrix \~chinese 变量矩阵
-    arma::vec mLatestWt;    //!< \~english Latest weigths \~chinese 最新的权重
+    arma::mat mX;           //!< \~english Variable matrix \~chinese 变量矩阵1
+    arma::vec mLatestWt;    //!< \~english Latest weigths \~chinese 最新的权重1
 
-    Solver mSolver = &GWPCA::solveSerial;   //!< \~english Calculator to solve \~chinese 模型求解函数
+    Solver mSolver = &GWPCA::solveSerial;   //!< \~english Calculator to solve \~chinese 模型求解函数1
 };
 
 }
 
-#endif  // GWPCA_H
+#endif  // GWPCA_H

src/gwmodelpp/GWPCA.cpp

@@ -10,47 +10,120 @@ void GWPCA::run()
     GWM_LOG_STOP_RETURN(mStatus, void());
 
     GWM_LOG_STAGE("Solving");
-    mLocalPV = pca(mX, mLoadings, mSDev);
+    mLocalPV = pca(mX, mLoadings, mScores, mSDev);
     GWM_LOG_STOP_RETURN(mStatus, void());
-    
+
     mWinner = index_max(mLoadings.slice(0), 1);
 }
 
-mat GWPCA::solveSerial(const mat& x, cube& loadings, mat& sdev)
+mat GWPCA::solveSerial(const mat& x, cube& loadings, cube& scores, mat& sdev)
 {
     uword nDp = mCoords.n_rows, nVar = mX.n_cols;
     mat d_all(nVar, nDp, arma::fill::zeros);
-    vec w0;
+
     loadings = cube(nDp, nVar, mK, arma::fill::zeros);
+    scores  = cube(nDp, mK, nDp, arma::fill::zeros);
+
     for (uword i = 0; i < nDp; i++)
     {
         GWM_LOG_STOP_BREAK(mStatus);
+
         vec w = mSpatialWeight.weightVector(i);
-        mat V;
+        uvec positive = find(w > 0);
+        vec newWt = w.elem(positive);
+        mat newX = x.rows(positive);
+        if (newWt.n_rows <= 5)
+        {
+            break;
+        }
+
+        mat U, V;
         vec d;
-        wpca(x, w, V, d);
-        w0 = w;
+        wpca(newX, newWt, U, V, d);
+
+        mLatestWt = newWt;
         d_all.col(i) = d;
+
         for (int j = 0; j < mK; j++)
         {
             loadings.slice(j).row(i) = arma::trans(V.col(j));
         }
+
+        mat scorei(nDp, mK, arma::fill::zeros);
+        for (int j = 0; j < mK; j++)
+        {
+            mat score = newX.each_row() % arma::trans(V.col(j));
+            scorei.col(j) = sum(score, 1);
+        }
+        scores.slice(i) = scorei;
+
         GWM_LOG_PROGRESS(i + 1, nDp);
     }
+
     d_all = trans(d_all);
-    mat variance = (d_all / sqrt(sum(w0))) % (d_all / sqrt(sum(w0)));
-    sdev = sqrt(variance);
+    mat variance = (d_all / pow(sum(mLatestWt), 0.5)) % (d_all / pow(sum(mLatestWt), 0.5));
+    sdev = arma::sqrt(variance);
     mat pv = variance.cols(0, mK - 1).each_col() % (1.0 / sum(variance, 1)) * 100.0;
     return pv;
 }
 
-void GWPCA::wpca(const mat& x, const vec& w, mat& V, vec & d)
+void GWPCA::wpca(const mat& x, const vec& w, mat& U, mat& V, vec & d)
 {
-    mat xw = x.each_col() % w, U;
+    mat xw = x.each_col() % w;
     mat centerized = (x.each_row() - sum(xw) / sum(w)).each_col() % sqrt(w);
     svd(U, d, V, centerized);
 }
 
+void GWPCA::rwpca(const mat& x, const vec& w, mat& U, mat& V, vec & d)
+{
+    mat mids = x;
+    uword medianIdx = (abs(w - 0.5)).index_min();
+    mids = mids.each_row() - x.row(medianIdx);
+    mat weighted = mids.each_col() % w;
+    mat score;
+    vec tsquared;
+    princomp(V, score, d, tsquared, weighted);
+    U = score;
+}
+
+mat GWPCA::solveRobustSerial(const mat& x, cube& loadings, cube& scores, mat& sdev)
+{
+    uword nDp = mCoords.n_rows, nVar = mX.n_cols;
+    mat d_all(nVar, nDp, arma::fill::zeros);
+    vec w0;
+    loadings = cube(nDp, nVar, mK, arma::fill::zeros);
+    scores = cube(nDp, nDp, mK, arma::fill::zeros);
+    for (uword i = 0; i < nDp; i++)
+    {
+        GWM_LOG_STOP_BREAK(mStatus);
+        vec w = mSpatialWeight.weightVector(i);
+        uvec positive = find(w > 0);
+        vec newWt = w.elem(positive);
+        mat newX = x.rows(positive);
+        if (newWt.n_rows <= 5)
+        {
+            continue;
+        }
+        mat U, V;
+        vec d;
+        rwpca(newX, newWt, U, V, d);
+        w0 = newWt;
+        d_all.col(i) = d;
+        for (int j = 0; j < mK; j++)
+        {
+            loadings.slice(j).row(i) = arma::trans(V.col(j));
+            mat scoreAll = x.each_row() % arma::trans(V.col(j));
+            scores.slice(j).col(i) = sum(scoreAll, 1);
+        }
+        GWM_LOG_PROGRESS(i + 1, nDp);
+    }
+    d_all = trans(d_all);
+    mat variance = (d_all / pow(sum(w0), 0.5)) % (d_all / pow(sum(w0), 0.5));
+    sdev = sqrt(variance);
+    mat pv = variance.cols(0, mK - 1).each_col() % (1.0 / sum(variance, 1)) * 100.0;
+    return pv;
+}
+
 bool GWPCA::isValid()
 {
     if (SpatialAlgorithm::isValid())

test/CMakeFind/main.cpp

@@ -6,14 +6,12 @@ using namespace arma;
 int main()
 {
     mat coords(100, 2, fill::randu);
-    mat x = join_rows(vec(100, fill::ones), mat(100, 2, fill::randu));
-    mat betas = mat(100, 3, fill::randu);
-    vec eps(100, fill::randu);
-    vec y = sum(x % betas, 1) + eps;
+    mat x = mat(100, 3, fill::randu);
     BandwidthWeight bw(36.0, true, BandwidthWeight::Gaussian);
     CRSDistance dist(false);
     SpatialWeight sw(&bw, &dist);
-    GWRBasic algorithm(x, y, coords, sw);
-    algorithm.fit();
+    GWPCA algorithm(x, coords, sw);
+    algorithm.setKeepComponents(2);
+    algorithm.run();
     return 0;
 }