From 4ff5fdd7d2a99c2c39a9ff4096641d48b0063dea Mon Sep 17 00:00:00 2001 From: Shaoqi Yang <1647246791@qq.com> Date: Sat, 17 Jan 2026 21:46:57 +0800 Subject: [PATCH 1/4] fix scores --- include/gwmodelpp/GWPCA.h | 21 ++++--- src/gwmodelpp/GWPCA.cpp | 14 +++-- test/CMakeFind/main.cpp | 10 ++-- test/testGWPCA.cpp | 117 +++++++++++++++++++++++++++++++++++++- 4 files changed, 140 insertions(+), 22 deletions(-) diff --git a/include/gwmodelpp/GWPCA.h b/include/gwmodelpp/GWPCA.h index 826d92ef..40d2d7f6 100644 --- a/include/gwmodelpp/GWPCA.h +++ b/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 @@ -86,7 +86,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,12 +105,13 @@ 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); } /** @@ -118,20 +119,22 @@ 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 solveSerial(const arma::mat& x, arma::cube& loadings, arma::mat& sdev); + arma::mat solveSerial(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& V, arma::vec & d); + void wpca(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 要保留的主成分数量 @@ -153,4 +156,4 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis } -#endif // GWPCA_H \ No newline at end of file +#endif // GWPCA_H diff --git a/src/gwmodelpp/GWPCA.cpp b/src/gwmodelpp/GWPCA.cpp index a2a954ac..5f61c9f1 100644 --- a/src/gwmodelpp/GWPCA.cpp +++ b/src/gwmodelpp/GWPCA.cpp @@ -10,30 +10,32 @@ 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, nDp, mK, arma::fill::zeros); for (uword i = 0; i < nDp; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = mSpatialWeight.weightVector(i); - mat V; + mat U, V; vec d; - wpca(x, w, V, d); + wpca(x, w, U, V, d); w0 = w; d_all.col(i) = d; for (int j = 0; j < mK; j++) { loadings.slice(j).row(i) = arma::trans(V.col(j)); + scores.slice(j).col(i) = U.col(j); } GWM_LOG_PROGRESS(i + 1, nDp); } @@ -44,9 +46,9 @@ mat GWPCA::solveSerial(const mat& x, cube& loadings, mat& sdev) 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); } diff --git a/test/CMakeFind/main.cpp b/test/CMakeFind/main.cpp index 1fb93361..e5939e88 100644 --- a/test/CMakeFind/main.cpp +++ b/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; } diff --git a/test/testGWPCA.cpp b/test/testGWPCA.cpp index 67d8bd8a..2c00b1a1 100644 --- a/test/testGWPCA.cpp +++ b/test/testGWPCA.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include "gwmodelpp/GWPCA.h" #include "gwmodelpp/spatialweight/CRSDistance.h" @@ -23,7 +24,6 @@ TEST_CASE("GWPCA: basic flow") { FAIL("Cannot load londonhp100 data."); } - CRSDistance distance(false); BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); SpatialWeight spatial(&bandwidth, &distance); @@ -78,6 +78,121 @@ TEST_CASE("GWPCA: basic flow") REQUIRE(approx_equal(loadings_pc2_q, loadings_pc2_q0, "absdiff", 1e-8)); } +TEST_CASE("GWPCA: summary statistics") +{ + mat londonhp100_coord, londonhp100_data; + vector londonhp100_fields; + if (!read_londonhp100(londonhp100_coord, londonhp100_data, londonhp100_fields)) + { + FAIL("Cannot load londonhp100 data."); + } + + CRSDistance distance(false); + BandwidthWeight bandwidth(15, true, BandwidthWeight::Gaussian); + SpatialWeight spatial(&bandwidth, &distance); + + mat x = londonhp100_data.cols(1, 3); + + GWPCA algorithm; + algorithm.setCoords(londonhp100_coord); + algorithm.setVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setKeepComponents(2); + REQUIRE_NOTHROW(algorithm.run()); + + vec p = {0.0, 0.25, 0.5, 0.75, 1.0}; + + mat localPV = algorithm.localPV(); + mat sdev = algorithm.sdev(); + mat localVariance = sdev % sdev; + + BandwidthWeight* bw_weight = algorithm.spatialWeight().weight(); + double bw_value = bw_weight->bandwidth(); + bool bw_adaptive = bw_weight->adaptive(); + BandwidthWeight::KernelFunctionType kernel_type = bw_weight->kernel(); + const char* bw_kernel_name = "Unknown"; + switch (kernel_type) + { + case BandwidthWeight::Gaussian: bw_kernel_name = "Gaussian"; break; + case BandwidthWeight::Exponential: bw_kernel_name = "Exponential"; break; + case BandwidthWeight::Bisquare: bw_kernel_name = "Bisquare"; break; + case BandwidthWeight::Tricube: bw_kernel_name = "Tricube"; break; + case BandwidthWeight::Boxcar: bw_kernel_name = "Boxcar"; break; + } + + printf("\n"); + printf("Summary of GWPCA information\n"); + printf("=====================================\n\n"); + printf("Bandwidth: %.0f (%s, %s kernel)\n", bw_value, bw_adaptive ? "adaptive" : "fixed", bw_kernel_name); + printf("\n"); + + printf("Local variance:\n"); + printf("%-15s %15s %15s\n", "", "Comp.1", "Comp.2"); + mat variance_q = quantile(localVariance, p, 0); + printf("%-15s %15.3f %15.3f\n", "Min", variance_q(0, 0), variance_q(0, 1)); + printf("%-15s %15.3f %15.3f\n", "1st Qu", variance_q(1, 0), variance_q(1, 1)); + printf("%-15s %15.3f %15.3f\n", "Median", variance_q(2, 0), variance_q(2, 1)); + printf("%-15s %15.3f %15.3f\n", "3rd Qu", variance_q(3, 0), variance_q(3, 1)); + printf("%-15s %15.3f %15.3f\n", "Max", variance_q(4, 0), variance_q(4, 1)); + printf("\n"); + + printf("Local Proportion of Variance:\n"); + printf("%-15s %15s %15s\n", "", "Comp.1", "Comp.2"); + mat pv_q = quantile(localPV, p, 0); + printf("%-15s %15.3f %15.3f\n", "Min", pv_q(0, 0), pv_q(0, 1)); + printf("%-15s %15.3f %15.3f\n", "1st Qu", pv_q(1, 0), pv_q(1, 1)); + printf("%-15s %15.3f %15.3f\n", "Median", pv_q(2, 0), pv_q(2, 1)); + printf("%-15s %15.3f %15.3f\n", "3rd Qu", pv_q(3, 0), pv_q(3, 1)); + printf("%-15s %15.3f %15.3f\n", "Max", pv_q(4, 0), pv_q(4, 1)); + + vec cumulative = sum(localPV, 1); + vec cum_q = quantile(cumulative, p); + printf("%-15s %15.3f %15s\n", "Cumulative", cum_q(0), ""); + printf("%-15s %15.3f %15s\n", "", cum_q(1), ""); + printf("%-15s %15.3f %15s\n", "", cum_q(2), ""); + printf("%-15s %15.3f %15s\n", "", cum_q(3), ""); + printf("%-15s %15.3f %15s\n", "", cum_q(4), ""); + printf("\n"); + + cube scores = algorithm.scores(); + if (scores.n_elem > 0) + { + printf("Local Scores:\n"); + printf("%-15s %15s %15s\n", "", "Comp.1", "Comp.2"); + + mat scores_pc1(scores.n_rows, scores.n_cols); + mat scores_pc2(scores.n_rows, scores.n_cols); + for (uword i = 0; i < scores.n_rows; i++) + { + for (uword j = 0; j < scores.n_cols; j++) + { + scores_pc1(i, j) = scores(i, j, 0); + if (scores.n_slices > 1) + { + scores_pc2(i, j) = scores(i, j, 1); + } + } + } + + vec scores_pc1_flat = vectorise(scores_pc1); + vec scores_pc2_flat = vectorise(scores_pc2); + + vec scores_pc1_q = quantile(scores_pc1_flat, p); + vec scores_pc2_q = quantile(scores_pc2_flat, p); + + printf("%-15s %15.3f %15.3f\n", "Min", scores_pc1_q(0), scores_pc2_q(0)); + printf("%-15s %15.3f %15.3f\n", "1st Qu", scores_pc1_q(1), scores_pc2_q(1)); + printf("%-15s %15.3f %15.3f\n", "Median", scores_pc1_q(2), scores_pc2_q(2)); + printf("%-15s %15.3f %15.3f\n", "3rd Qu", scores_pc1_q(3), scores_pc2_q(3)); + printf("%-15s %15.3f %15.3f\n", "Max", scores_pc1_q(4), scores_pc2_q(4)); + printf("\n"); + } + else + { + printf("Local Scores: (not yet implemented)\n\n"); + } +} + TEST_CASE("GWSS: cancel") { From 7423fad438cb93ee96166000f30bb1c238bf4d27 Mon Sep 17 00:00:00 2001 From: Shaoqi Yang <1647246791@qq.com> Date: Sat, 17 Jan 2026 22:55:46 +0800 Subject: [PATCH 2/4] fix robust --- include/gwmodelpp/GWPCA.h | 58 ++- src/gwmodelpp/GWPCA.cpp | 50 ++ test/gwmgwpcataskthread.cpp | 939 ++++++++++++++++++++++++++++++++++++ test/gwmgwpcataskthread.h | 223 +++++++++ 4 files changed, 1259 insertions(+), 11 deletions(-) create mode 100644 test/gwmgwpcataskthread.cpp create mode 100644 test/gwmgwpcataskthread.h diff --git a/include/gwmodelpp/GWPCA.h b/include/gwmodelpp/GWPCA.h index 40d2d7f6..ab2ce2ce 100644 --- a/include/gwmodelpp/GWPCA.h +++ b/include/gwmodelpp/GWPCA.h @@ -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 获取局部主成分值。 * @@ -115,7 +129,18 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis } /** - * @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 载荷矩阵 @@ -123,7 +148,7 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis * @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::cube& scores, 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的函数。 @@ -136,22 +161,33 @@ class GWPCA: public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis */ 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 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 }; } diff --git a/src/gwmodelpp/GWPCA.cpp b/src/gwmodelpp/GWPCA.cpp index 5f61c9f1..ce2845b8 100644 --- a/src/gwmodelpp/GWPCA.cpp +++ b/src/gwmodelpp/GWPCA.cpp @@ -53,6 +53,56 @@ void GWPCA::wpca(const mat& x, const vec& w, mat& U, mat& V, vec & d) 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()) diff --git a/test/gwmgwpcataskthread.cpp b/test/gwmgwpcataskthread.cpp new file mode 100644 index 00000000..fe8438ba --- /dev/null +++ b/test/gwmgwpcataskthread.cpp @@ -0,0 +1,939 @@ +#include "gwmgwpcataskthread.h" +#include +#include "TaskThread/gwmgeographicalweightedregressionalgorithm.h" +#include "gwmtaskthread.h" +#ifdef ENABLE_OpenMP +#include +#endif +#include +#include +#include +#include +#include +#include +#include + +int GwmGWPCATaskThread::treeChildCount = 0; + +GwmGWPCATaskThread::GwmGWPCATaskThread() : GwmSpatialMonoscaleAlgorithm(), + mAlgorithm(std::make_unique()) +{ + +} + +void GwmGWPCATaskThread::setCanceled(bool canceled) +{ + // mSelector.setCanceled(canceled); + // mSpatialWeight.distance()->setCanceled(canceled); + return GwmTaskThread::setCanceled(canceled); +} + +void GwmGWPCATaskThread::run() +{ + if(!checkCanceled()) + { + emit message(QString(tr("Setting data points ..."))); + initPoints(); + + emit message(QString(tr("Setting X and Y."))); + initXY(mX,mVariables); + + if(zscore()){ + emit message(QString(tr("Zscore normalizaiton..."))); + variableZscore(mX); + } + } + + if(mIsAutoselectBandwidth && !checkCanceled()) + { + emit message(QString(tr("Automatically selecting bandwidth ..."))); + emit tick(0, 0); + + if ((mSpatialWeight.distance()->type() == gwm::Distance::CRSDistance || + mSpatialWeight.distance()->type() == gwm::Distance::MinkwoskiDistance) && !checkCanceled()) + { + gwm::CRSDistance* d = static_cast(mSpatialWeight.distance()); + if(d) + { + d->makeParameter({ mDataPoints, mDataPoints }); + } + } + + gwm::BandwidthWeight* bandwidthWeight0 = mSpatialWeight.weight(); + if(!bandwidthWeight0) + { + qDebug() << "[GWPCA] ERROR: Cannot get bandwidth weight!"; + emit error(tr("Cannot get bandwidth weight for bandwidth selection.")); + return; + } + + double tmpMaxD = mSpatialWeight.distance()->maxDistance(); + double lower = bandwidthWeight0->adaptive() ? 2 : tmpMaxD / 5000; + double upper = bandwidthWeight0->adaptive() ? mDataPoints.n_rows : tmpMaxD; + + if(bandwidthWeight0->bandwidth() <= 0 || + (bandwidthWeight0->adaptive() && bandwidthWeight0->bandwidth() < lower) || + (!bandwidthWeight0->adaptive() && bandwidthWeight0->bandwidth() < lower)) + { + double initialBandwidth = bandwidthWeight0->adaptive() ? + std::max(2.0, std::min(20.0, (double)mDataPoints.n_rows * 0.1)) : + std::max(lower, tmpMaxD * 0.1); + bandwidthWeight0->setBandwidth(initialBandwidth); + } + + mSelector.setBandwidth(bandwidthWeight0); + mSelector.setLower(lower); + mSelector.setUpper(upper); + + try + { + gwm::BandwidthWeight* bandwidthWeight = mSelector.optimize(this); + if(bandwidthWeight && !checkCanceled()) + { + mSpatialWeight.setWeight(bandwidthWeight); + mSelector.setBandwidth(bandwidthWeight); + + gwm::BandwidthWeight* verifyBw = mSpatialWeight.weight(); + if(verifyBw && verifyBw->bandwidth() == 0) + { + qDebug() << "[GWPCA] ERROR: Bandwidth is 0 after setWeight! This may cause display issues."; + } + } + else if(!bandwidthWeight) + { + qDebug() << "[GWPCA] WARNING: Bandwidth optimization returned NULL"; + emit error(tr("Bandwidth optimization failed: no optimal bandwidth found.")); + return; + } + } + catch(const std::exception& e) + { + qDebug() << "[GWPCA] EXCEPTION during bandwidth optimization:" << e.what(); + emit error(QString(tr("Bandwidth optimization error: %1")).arg(e.what())); + return; + } + catch(...) + { + qDebug() << "[GWPCA] UNKNOWN EXCEPTION during bandwidth optimization"; + emit error(tr("Unknown error occurred during bandwidth optimization.")); + return; + } + } + + if(!checkCanceled()) + { + emit message(QString(tr("Principle components analyzing ..."))); + + try + { + if(Robust()) + { + emit message(QString(tr("Running Robust GWPCA ..."))); + mLocalPV = robustSolveSerial(mX, mLoadings, mSDev); + + if(checkCanceled()) + { + return; + } + + mVariance = mSDev % mSDev; + } + else + { + if(!mAlgorithm) + { + qDebug() << "[GWPCA] ERROR: mAlgorithm is NULL!"; + emit error(tr("GWPCA algorithm object is not initialized.")); + return; + } + + mAlgorithm->setCoords(mDataPoints); + mAlgorithm->setVariables(mX); + mAlgorithm->setSpatialWeight(mSpatialWeight); + mAlgorithm->setKeepComponents(mK); + + if(!mAlgorithm->isValid()) + { + qDebug() << "[GWPCA] ERROR: Algorithm configuration is invalid!"; + emit error(tr("GWPCA algorithm configuration is invalid.")); + return; + } + + mAlgorithm->setTelegram(std::make_unique(this)); + mAlgorithm->run(); + + auto status = mAlgorithm->status(); + if(status != gwm::Status::Success) + { + qDebug() << "[GWPCA] ERROR: Algorithm execution failed with status:" << static_cast(status); + emit error(tr("GWPCA algorithm execution failed.")); + return; + } + + mLocalPV = mAlgorithm->localPV(); + mLoadings = mAlgorithm->loadings(); + mSDev = mAlgorithm->sdev(); + mVariance = mSDev % mSDev; + + // Get scores from kernel library for comparison + if(scoresCal() && mDataPoints.n_rows <= 4096) + { + const cube& kernelScores = mAlgorithm->scores(); + mScoresFromKernel = kernelScores; + + qDebug() << "[GWPCA] Got scores from kernel library, dimensions:" << mScoresFromKernel.n_rows << "x" << mScoresFromKernel.n_cols << "x" << mScoresFromKernel.n_slices; + } + } + + if(scoresCal() && mDataPoints.n_rows <= 4096) + { + // Keep original calculateScores() result for comparison + calculateScores(); + + // Output comparison information + if(!Robust() && mScoresFromKernel.n_elem > 0) + { + qDebug() << "[GWPCA] ========== Scores Comparison Info =========="; + qDebug() << "[GWPCA] Original code scores dimensions:" << mScores.n_rows << "x" << mScores.n_cols << "x" << mScores.n_slices; + qDebug() << "[GWPCA] Kernel library scores dimensions:" << mScoresFromKernel.n_rows << "x" << mScoresFromKernel.n_cols << "x" << mScoresFromKernel.n_slices; + qDebug() << "[GWPCA] Note: Original code scores format is (nDp, mK, nDp), kernel library scores format is (nDp, nDp, mK)"; + qDebug() << "[GWPCA] Original code: mScores.slice(i) is the score matrix (nDp, mK) for point i"; + qDebug() << "[GWPCA] Kernel library: mScoresFromKernel.slice(j) is the score matrix (nDp, nDp) for component j"; + + int sampleSize = std::min(5, (int)mDataPoints.n_rows); + int sampleComps = std::min(mK, 3); + qDebug() << "[GWPCA] First" << sampleSize << "points score comparison:"; + for(int i = 0; i < sampleSize; i++) + { + qDebug() << "[GWPCA] Point" << i << "scores:"; + qDebug() << "[GWPCA] Original code (mScores.slice(" << i << ")):"; + for(int j = 0; j < sampleComps; j++) + { + vec scoreVec = mScores.slice(i).col(j); + QString scoreStr = "["; + for(uword k = 0; k < std::min((uword)5, scoreVec.n_elem); k++) + { + if(k > 0) scoreStr += ", "; + scoreStr += QString::number(scoreVec(k), 'f', 6); + } + if(scoreVec.n_elem > 5) scoreStr += ", ..."; + scoreStr += "]"; + qDebug() << "[GWPCA] Component" << (j+1) << ":" << scoreStr; + } + qDebug() << "[GWPCA] Kernel library (column" << i << "of each component):"; + for(int j = 0; j < sampleComps; j++) + { + vec scoreVec = mScoresFromKernel.slice(j).col(i); + QString scoreStr = "["; + for(uword k = 0; k < std::min((uword)5, scoreVec.n_elem); k++) + { + if(k > 0) scoreStr += ", "; + scoreStr += QString::number(scoreVec(k), 'f', 6); + } + if(scoreVec.n_elem > 5) scoreStr += ", ..."; + scoreStr += "]"; + qDebug() << "[GWPCA] Component" << (j+1) << ":" << scoreStr; + } + } + + qDebug() << "[GWPCA] Format converted comparison:"; + uword nDp = mDataPoints.n_rows; + cube kernelScoresConverted(nDp, mK, nDp, fill::zeros); + for(uword i = 0; i < nDp; i++) + { + for(int j = 0; j < mK; j++) + { + kernelScoresConverted.slice(i).col(j) = mScoresFromKernel.slice(j).col(i); + } + } + + for(int j = 0; j < sampleComps; j++) + { + double maxDiff = 0, meanDiff = 0; + uword diffCount = 0; + for(uword i = 0; i < nDp; i++) + { + vec diff = abs(mScores.slice(i).col(j) - kernelScoresConverted.slice(i).col(j)); + double localMax = diff.max(); + double localMean = mean(diff); + if(localMax > maxDiff) maxDiff = localMax; + meanDiff += localMean; + diffCount++; + } + meanDiff /= diffCount; + qDebug() << "[GWPCA] Component" << (j+1) << "score difference - Max diff:" << maxDiff << "Mean diff:" << meanDiff; + } + + qDebug() << "[GWPCA] ====================================="; + } + } + else + { + mScoresCal = false; + } + } + catch(const std::exception& e) + { + qDebug() << "[GWPCA] EXCEPTION caught:" << e.what(); + emit error(QString(tr("GWPCA execution error: %1")).arg(e.what())); + return; + } + catch(...) + { + qDebug() << "[GWPCA] UNKNOWN EXCEPTION caught"; + emit error(tr("Unknown error occurred during GWPCA execution.")); + return; + } + } + + QList win_var_PC1; + uvec iWinVar = index_max(mLoadings.slice(0), 1); + for(int i = 0; i < mDataPoints.n_rows && !checkCanceled(); i++) + { + win_var_PC1.append(mVariables.at(iWinVar(i)).name); + } + + if(!checkCanceled()) + { + CreateResultLayerData resultLayerData = { + qMakePair(QString("Comp.%1_PV"), mLocalPV), + qMakePair(QString("local_CP"), sum(mLocalPV, 1)) + }; + if(scoresCal()){ + for(int i = 0; i < mK; i++){ + resultLayerData += qMakePair(QString("Scores%1"), mScores.slice(i)); + } + }; + createResultLayer(resultLayerData,win_var_PC1); + + if(getPlot()){ + CreatePlotLayerData plotLayerData = {}; + vec vecLoadings(mDataPoints.n_rows * mVariables.size()); + QList var_pc; + int v = 0; + for(int i = 0; i < mDataPoints.n_rows; i++){ + for(int j=0; j < mVariables.size(); j++ ){ + vecLoadings(v) = mLoadings.slice(0)(i,j); + var_pc.append(mVariables.at(j).name); + v++; + } + } + plotLayerData += qMakePair(QStringLiteral("loadings"),vecLoadings); + createPlotLayer(plotLayerData, var_pc); + } + emit success(); + emit tick(100, 100); + } + if(checkCanceled()) + { + return; + } +} + +bool GwmGWPCATaskThread::isValid() +{ + gwm::BandwidthWeight* bandwidth = static_cast(mSpatialWeight.weight()); + if(bandwidth){ + if(!mIsAutoselectBandwidth) + { + if(bandwidth->adaptive()){ + if(bandwidth->bandwidth() <= mVariables.size()){ + return false; + } + } + } + if(mVariables.size() == 0){ + return false; + } + if(k()<=0 || k() > mVariables.size()){ + return false; + } + }else{ + return false; + } + return true; +} + +void GwmGWPCATaskThread::initPoints() +{ + if(!mDataLayer) + { + qDebug() << "[GWPCA::initPoints] ERROR: Data layer is NULL!"; + return; + } + + int nDp = mDataLayer->featureCount(); + mDataPoints = mat(nDp, 2, fill::zeros); + + QgsFeatureIterator iterator = mDataLayer->getFeatures(); + QgsFeature f; + for (int i = 0; iterator.nextFeature(f); i++) + { + QgsPointXY centroPoint = f.geometry().centroid().asPoint(); + mDataPoints(i, 0) = centroPoint.x(); + mDataPoints(i, 1) = centroPoint.y(); + } +} + +void GwmGWPCATaskThread::initXY(mat &x, const QList &indepVars) +{ + int nDp = mDataLayer->featureCount(), nVar = indepVars.size(); + + x = mat(nDp, nVar, fill::zeros); + + QgsFeatureIterator iterator = mDataLayer->getFeatures(); + QgsFeature f; + bool ok = false; + int errorCount = 0; + for (int i = 0; iterator.nextFeature(f); i++) + { + for (int k = 0; k < indepVars.size(); k++) + { + double vX = f.attribute(indepVars[k].name).toDouble(&ok); + if (ok) + { + x(i, k) = vX; + } + else + { + errorCount++; + if(errorCount <= 5) + { + qDebug() << "[GWPCA::initXY] WARNING: Cannot convert variable" << indepVars[k].name << "at row" << i << "to number"; + } + emit error(tr("Independent variable value cannot convert to a number. Set to 0.")); + } + } + } +} + +void GwmGWPCATaskThread::variableZscore(mat& x) +{ + mat xmean = mean(x); + mat xstd = stddev(x); + for (int k = 0; k < x.n_cols; k++) + { + x.col(k) = (x.col(k) - xmean(k))/xstd(k); + } +} + + +void GwmGWPCATaskThread::calculateScores() +{ + uword nDp = mDataPoints.n_rows, nVar = mVariables.size(); + mScores = cube(nDp, mK, nDp, fill::zeros); + + for(uword i = 0; i < nDp && !checkCanceled(); i++) + { + vec wt = mSpatialWeight.weightVector(i); + uvec positive = find(wt > 0); + vec newWt = wt.elem(positive); + mat newX = mX.rows(positive); + + if(newWt.n_rows <= 5) + { + break; + } + + mat V; + vec d; + if(!Robust()){ + wpca(newX, newWt, V, d); + }else{ + rwpca(newX, newWt, V, d); + } + + mat scorei(nDp, mK, fill::zeros); + for(int j = 0; j < mK && !checkCanceled(); j++) + { + mat score = newX.each_row() % trans(V.col(j)); + scorei.col(j) = sum(score, 1); + } + mScores.slice(i) = scorei; + emit tick(i + 1, nDp); + } +} + +void GwmGWPCATaskThread::wpca(const mat &x, const vec &wt, mat &V, vec &S) +{ + + mat xw = x.each_col() % wt; + mat centerized = (x.each_row() - sum(xw) / sum(wt)).each_col() % sqrt(wt); + //SVD + mat U; + svd(U,S,V,centerized); + //S即为R中的d + //V即为R中的v +} + +void GwmGWPCATaskThread::rwpca(const mat &x, const vec &wt, mat &V, vec &S) +{ + + mat mids = x; + mids = mids.each_row() - x.row((abs(wt - 0.5)).index_min()); + + mat score; + vec tsquared; + princomp(V, score, S, tsquared, mids.each_col() % wt); + +} + +mat GwmGWPCATaskThread::robustSolveSerial(const mat& x, cube& loadings, mat& sdev) +{ + int nDp = mDataPoints.n_rows, nVar = mX.n_cols; + + mat d_all(nVar, nDp, fill::zeros); + + loadings = cube(nDp, nVar, mK, fill::zeros); + + for(int i=0;idistance(i); + vec wt = mSpatialWeight.weightVector(i); + + uvec positive = find(wt > 0); + vec newWt = wt.elem(positive); + mat newX = x.rows(positive); + if(newWt.n_rows<=5) + { + break; + } + + mat V; + vec d; + rwpca(newX,newWt,V,d); + + mLatestWt = newWt; + d_all.col(i) = d; + + for(int j = 0; j < mK && !checkCanceled(); j++) + { + loadings.slice(j).row(i) = trans(V.col(j)); + } + emit tick(i, nDp); + } + + d_all = trans(d_all); + mat variance = (d_all / pow(sum(mLatestWt),0.5)) % (d_all / pow(sum(mLatestWt),0.5)); + + sdev = sqrt(variance); + + mat pv = variance.cols(0, mK - 1).each_col() % (1.0 / sum(variance,1)) * 100.0; + return pv; +} + +gwm::Status GwmGWPCATaskThread::getCriterion(gwm::BandwidthWeight* weight, double& criterion) +{ + if(checkCanceled()) + { + criterion = DBL_MAX; + return gwm::Status::Terminated; + } + + // 将内核库的BandwidthWeight转换为本地的GwmBandwidthWeight + GwmBandwidthWeight::KernelFunctionType kernelType = static_cast(weight->kernel()); + GwmBandwidthWeight localWeight(weight->bandwidth(), weight->adaptive(), kernelType); + + try + { + criterion = (this->*mBandwidthSelectCriterionFunction)(&localWeight); + return gwm::Status::Success; + } + catch(const std::exception& e) + { + qDebug() << "[GWPCA::getCriterion] EXCEPTION:" << e.what(); + criterion = DBL_MAX; + return gwm::Status::Success; + } + catch(...) + { + qDebug() << "[GWPCA::getCriterion] UNKNOWN EXCEPTION"; + criterion = DBL_MAX; + return gwm::Status::Success; + } +} + +void GwmGWPCATaskThread::setBandwidthSelectionCriterionType(const GwmGWPCATaskThread::BandwidthSelectionCriterionType &bandwidthSelectionCriterionType) +{ + mBandwidthSelectionCriterionType = bandwidthSelectionCriterionType; + QMap, BandwidthSelectCriterionFunction> mapper = { + std::make_pair(qMakePair(BandwidthSelectionCriterionType::CV, IParallelalbe::ParallelType::SerialOnly), &GwmGWPCATaskThread::bandwidthSizeCriterionCVSerial), + #ifdef ENABLE_OpenMP + std::make_pair(qMakePair(BandwidthSelectionCriterionType::CV, IParallelalbe::ParallelType::OpenMP), &GwmGWPCATaskThread::bandwidthSizeCriterionCVOmp), + #endif + //std::make_pair(qMakePair(BandwidthSelectionCriterionType::CV, IParallelalbe::ParallelType::CUDA), &GwmGWPCATaskThread::bandwidthSizeCriterionCVCuda), + //std::make_pair(qMakePair(BandwidthSelectionCriterionType::AIC, IParallelalbe::ParallelType::SerialOnly), &GwmGWPCATaskThread::bandwidthSizeCriterionAICSerial), + //std::make_pair(qMakePair(BandwidthSelectionCriterionType::AIC, IParallelalbe::ParallelType::OpenMP), &GwmGWPCATaskThread::bandwidthSizeCriterionAICOmp), + //std::make_pair(qMakePair(BandwidthSelectionCriterionType::AIC, IParallelalbe::ParallelType::CUDA), &GwmGWPCATaskThread::bandwidthSizeCriterionAICCuda) + }; + mBandwidthSelectCriterionFunction = mapper[qMakePair(bandwidthSelectionCriterionType, mParallelType)]; +} + +bool GwmGWPCATaskThread::isAutoselectBandwidth() const +{ + return mIsAutoselectBandwidth; +} + +void GwmGWPCATaskThread::setIsAutoselectBandwidth(bool isAutoselectBandwidth) +{ + mIsAutoselectBandwidth = isAutoselectBandwidth; +} + +void GwmGWPCATaskThread::setVariables(const QList &variables) +{ + mVariables = variables; +} + +void GwmGWPCATaskThread::setParallelType(const IParallelalbe::ParallelType &type) +{ + if(type & parallelAbility()) + { + mParallelType = type; + setBandwidthSelectionCriterionType(mBandwidthSelectionCriterionType); + + } +} + +void GwmGWPCATaskThread::createResultLayer(CreateResultLayerData data, QList winvar) +{ + QgsVectorLayer* srcLayer = mDataLayer; + QString layerFileName = QgsWkbTypes::displayString(srcLayer->wkbType()) + QStringLiteral("?"); + QString layerName = srcLayer->name(); + + + if(treeChildCount > 0) + { + if(!Robust()) layerName += QStringLiteral("_GWPCA") + "(" + QString::number(treeChildCount) + ")"; + else if(Robust()) layerName += QStringLiteral("_RGWPCA") + "(" + QString::number(treeChildCount) + ")"; + } else + { + if(!Robust()) layerName += QStringLiteral("_GWPCA"); + else if(Robust()) layerName += QStringLiteral("_RGWPCA"); + } + + treeChildCount++ ; + + +// if(!Robust()){ +// layerName += QStringLiteral("_GWPCA"); +// }else{ +// layerName += QStringLiteral("_RGWPCA"); +// } + + mResultLayer = new QgsVectorLayer(layerFileName, layerName, QStringLiteral("memory")); + mResultLayer->setCrs(srcLayer->crs()); + + + QgsFields fields; + for (QPair item : data) + { + QString title = item.first; + const mat& value = item.second; + if (value.n_cols > 1) + { + for (int k = 0; k < value.n_cols; k++) + { + QString fieldName = title.arg(k+1); + fields.append(QgsField(fieldName, QVariant::Double, QStringLiteral("double"))); + } + } + else + { + fields.append(QgsField(title, QVariant::Double, QStringLiteral("double"))); + } + } + fields.append(QgsField("win_var_PC1",QVariant::String,QStringLiteral("varchar"),255)); + mResultLayer->dataProvider()->addAttributes(fields.toList()); + mResultLayer->updateFields(); + + + mResultLayer->startEditing(); + QgsFeatureIterator iterator = srcLayer->getFeatures(); + QgsFeature f; + for (int i = 0; iterator.nextFeature(f); i++) + { + QgsFeature feature(fields); + feature.setGeometry(f.geometry()); + + int k = 0; + for (QPair item : data) + { + for (uword d = 0; d < item.second.n_cols; d++) + { + feature.setAttribute(k, item.second(i, d)); + k++; + } + } + feature.setAttribute("win_var_PC1",winvar[i]); + + mResultLayer->addFeature(feature); + } + mResultLayer->commitChanges(); +} + +void GwmGWPCATaskThread::createPlotLayer(CreatePlotLayerData data, QList varpc) +{ + QgsVectorLayer* srcLayer = mDataLayer; + QString layerFileName = QgsWkbTypes::displayString(srcLayer->wkbType()) + QStringLiteral("?"); + QString layerName = srcLayer->name(); + + + if(treeChildCount > 0) + { + layerName += QStringLiteral("_GlyphPlot") + "(" + QString::number(treeChildCount) + ")"; + } else + { + layerName += QStringLiteral("_GlyphPlot"); + } + + treeChildCount++ ; + + + if(!Robust()){ + layerName += QStringLiteral("_GWPCA"); + }else{ + layerName += QStringLiteral("_RGWPCA"); + } + + mPlotLayer = new QgsVectorLayer(QgsWkbTypes::displayString(QgsWkbTypes::LineString), layerName, QStringLiteral("memory")); + mPlotLayer->setCrs(srcLayer->crs()); + + QgsFields fields; +// for (QPair item : data) +// { +// QString title = item.first; +// const mat& value = item.second; +// if (value.n_cols > 1) +// { +// for (int k = 0; k < value.n_cols; k++) +// { +// QString fieldName = title.arg(k+1); +// fields.append(QgsField(fieldName, QVariant::Double, QStringLiteral("double"))); +// } +// } +// else +// { +// fields.append(QgsField(title, QVariant::Double, QStringLiteral("double"))); +// } +// } + fields.append(QgsField(QStringLiteral("loadings"), QVariant::Double, QStringLiteral("double"))); + fields.append(QgsField(QStringLiteral("var_name"),QVariant::String,QStringLiteral("varchar"),255)); + mPlotLayer->dataProvider()->addAttributes(fields.toList()); + mPlotLayer->updateFields(); + + mPlotLayer->startEditing(); + QgsFeatureIterator iterator = srcLayer->getFeatures(); + QgsFeature f; + double PI = 3.1415926535898; + int k = 0; + for (int i = 0; i item : data) +// { +// for (uword d = 0; d < item.second.n_cols; d++) +// { +// feature.setAttribute(k, item.second(i, d)); +// k++; +// } +// } + + for (QPair item : data) + { + for (uword d = 0; d < item.second.n_cols; d++) + { + + double x = sttPoint.x()+fabs(item.second(i, d)) * cos(k * angle)*10; + double y = sttPoint.y()+fabs(item.second(i, d)) * sin(k * angle)*10; +// std::cout< lineData ={}; + lineData+= sttPoint; + lineData+= QgsPointXY(x,y); +// std::cout<<"x: "<addFeature(feature); + mPlotLayer->addFeature(lineFeature); + } + mPlotLayer->commitChanges(); + qWarning("test %s", (QgsWkbTypes::displayString(mPlotLayer->wkbType())).toStdString().data()); +} + +double GwmGWPCATaskThread::bandwidthSizeCriterionCVSerial(GwmBandwidthWeight *weight) +{ + int mBandwidthCounter = 0; + int n = mX.n_rows; + int m = mX.n_cols; + double score = 0; + + if (mSpatialWeight.distance()->type() == gwm::Distance::CRSDistance || + mSpatialWeight.distance()->type() == gwm::Distance::MinkwoskiDistance) + { + gwm::CRSDistance* d = static_cast(mSpatialWeight.distance()); + if(d) + { + d->makeParameter({ mDataPoints, mDataPoints }); + } + } + + for (int i = 0; i < n && !checkCanceled(); i++) + { + vec distvi = mSpatialWeight.distance()->distance(i); + if(distvi.n_elem == 0) + { + qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] ERROR: Empty distance vector at point" << i; + score = DBL_MAX; + break; + } + + vec wt = weight->weight(distvi); + wt(i) = 0; + + uvec positive = find(wt > 0); + if(positive.n_elem == 0) + { + qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] WARNING: No positive weights at point" << i; + score = DBL_MAX; + break; + } + + vec newWt = wt.elem(positive); + mat newX = mX.rows(positive); + //判断length(newWt) + if(newWt.n_rows <= 1) + { + qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] WARNING: Insufficient points at" << i << ", n_rows:" << newWt.n_rows; + score = DBL_MAX; + break; + } + + mat V; + vec S; + try + { + wpca(newX, newWt, V, S); + if(V.n_cols < mK) + { + qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] WARNING: V.n_cols < mK at point" << i; + score = DBL_MAX; + break; + } + V = V.cols(0, mK - 1); + V = V * trans(V); + score = score + pow(sum(mX.row(i) - mX.row(i) * V),2); + } + catch(const std::exception& e) + { + qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] EXCEPTION at point" << i << ":" << e.what(); + score = DBL_MAX; + break; + } + + mBandwidthCounter++; + if (mBandwidthCounter < 10) + emit tick(mBandwidthCounter * 10 + i * 5 / n, 100); + } + return score; +} +#ifdef ENABLE_OpenMP +double GwmGWPCATaskThread::bandwidthSizeCriterionCVOmp(GwmBandwidthWeight *weight) +{ + int n = mX.n_rows; + int m = mX.n_cols; + double score = 0; + bool flag = true; + vec score_all(mOmpThreadNum, fill::zeros); + int current = 0; +#pragma omp parallel for num_threads(mOmpThreadNum) + for (int i = 0; i < n; i++) + { + if(flag && !checkCanceled()) + { + int thread = omp_get_thread_num(); + vec distvi = mSpatialWeight.distance()->distance(i); + vec wt = weight->weight(distvi); + wt(i) = 0; + uvec positive = find(wt > 0); + vec newWt = wt.elem(positive); + mat newX = mX.rows(positive); + //判断length(newWt) + if(newWt.n_rows <=1) + { + flag=false; + }else{ + mat V; + vec S; + // 带宽选择时总是使用普通GWPCA（wpca），无论是否Robust模式 + // 这样Robust GWPCA会使用与普通GWPCA相同的带宽 + wpca(newX, newWt, V, S); + V = V.cols(0, mK - 1); + V = V * trans(V); + score_all(thread) += pow(sum(mX.row(i) - mX.row(i) * V),2); + } + if(mSelector.counter<10) + emit tick(mSelector.counter * 10 + current * 10 / n, 100); + current++; + } + } + score = sum(score_all); + return score; +} +#endif +bool GwmGWPCATaskThread::zscore() const +{ + return mZscore; +} + +void GwmGWPCATaskThread::setZscore(bool zscore) +{ + mZscore = zscore; +} + +bool GwmGWPCATaskThread::scoresCal() const +{ + return mScoresCal; +} + +void GwmGWPCATaskThread::setScoresCal(bool scoresCal) +{ + mScoresCal = scoresCal; +} + +bool GwmGWPCATaskThread::Robust() const +{ + return mRobust; +} + +void GwmGWPCATaskThread::setRobust(bool robust) +{ + mRobust=robust; +} + +bool GwmGWPCATaskThread::getPlot() const +{ + return mPlot; +} + +void GwmGWPCATaskThread::setPlot(bool plot) +{ + mPlot=plot; +} + diff --git a/test/gwmgwpcataskthread.h b/test/gwmgwpcataskthread.h new file mode 100644 index 00000000..3cda47ac --- /dev/null +++ b/test/gwmgwpcataskthread.h @@ -0,0 +1,223 @@ +#ifndef GWMGWPCATASKTHREAD_H +#define GWMGWPCATASKTHREAD_H + +#include +#include "TaskThread/gwmspatialmonoscalealgorithm.h" +#include "TaskThread/imultivariableanalysis.h" +#include "TaskThread/iparallelable.h" + +#include "TaskThread/gwmbandwidthsizeselector.h" +#include + +class GwmGWPCATaskThread : public GwmSpatialMonoscaleAlgorithm, public IBandwidthSizeSelectable, public IGwmMultivariableAnalysis, public IOpenmpParallelable, public gwm::IBandwidthSelectable +{ + Q_OBJECT + enum BandwidthSelectionCriterionType + { + CV + }; + + typedef QList > CreateResultLayerData; + typedef QList > CreatePlotLayerData; + + typedef double (GwmGWPCATaskThread::*BandwidthSelectCriterionFunction)(GwmBandwidthWeight*); + +public: + GwmGWPCATaskThread(); + +public: + QList variables() const override + { + return QList(); + }; + void setVariables(const QList &variables) override; + void setVariables(const QList &&variables) + { + mVariables = variables; + }; + + bool isAutoselectBandwidth() const; + void setIsAutoselectBandwidth(bool isAutoselectBandwidth); + BandwidthSelectionCriterionType bandwidthSelectionCriterionType() const; + void setBandwidthSelectionCriterionType(const BandwidthSelectionCriterionType &bandwidthSelectionCriterionType); + + void setOmpThreadNum(const int threadNum) override + { + mOmpThreadNum = threadNum; + }; + + BandwidthCriterionList bandwidthSelectorCriterions() const + { + return mSelector.bandwidthCriterion(); + } + + double k() const; + void setK(double k); + + mat localPV() const; + + mat variance() const; + + mat sdev() const; + + cube loadings() const; + + cube scores() const; + +public: // QThread interface + void run() override; + + +public: // GwmTaskThread interface + QString name() const override { return tr("GWPCA"); } + + +public: // GwmSpatialMonoscaleAlgorithm interface + bool isValid() override; + + +public: // IParallelalbe interface + int parallelAbility() const override; + virtual ParallelType parallelType() const override; + virtual void setParallelType(const ParallelType& type) override; + +public: // IOpenmpParallelable interface + void setThreadNum(const int threadNum){}; + //void setOmpThreadNum(const int threadNum){}; + + bool zscore() const; + void setZscore(bool zscore); + + bool scoresCal() const; + void setScoresCal(bool scoresCal); + + bool Robust() const; + void setRobust(bool robust); + + bool getPlot() const; + void setPlot(bool plot); + + void setCanceled(bool canceled) override; + +private: + void initPoints(); + void initXY(mat& x, const QList& indepVars); + void variableZscore(mat& x); + void calculateScores(); + + void createResultLayer(CreateResultLayerData data,QList winvar); + void createPlotLayer(CreatePlotLayerData data, QList varpc); + std::unique_ptr mAlgorithm; + + // wpca函数保留，用于带宽选择（内核库的wpca是私有的，无法直接访问） + void wpca(const mat &x, const vec &wt, mat &V, vec &S); + + // Robust weighted PCA 函数 + void rwpca(const mat &x, const vec &wt, mat &V, vec &S); + + // Robust GWPCA 求解函数 + mat robustSolveSerial(const mat& x, cube& loadings, mat& sdev); + + double bandwidthSizeCriterionCVSerial(GwmBandwidthWeight* weight); +#ifdef ENABLE_OpenMP + double bandwidthSizeCriterionCVOmp(GwmBandwidthWeight* weight); +#endif + double criterion(GwmBandwidthWeight *weight) override + { + return (this->*mBandwidthSelectCriterionFunction)(weight); + } + +public: // gwm::IBandwidthSelectable interface + gwm::Status getCriterion(gwm::BandwidthWeight* weight, double& criterion) override; + +private: + QList mVariables; + mat mDataPoints; + + int mK = 2; + + mat mX; + vec mLatestWt; + + gwm::BandwidthSelector mSelector; + BandwidthSelectionCriterionType mBandwidthSelectionCriterionType = BandwidthSelectionCriterionType::CV; + BandwidthSelectCriterionFunction mBandwidthSelectCriterionFunction = &GwmGWPCATaskThread::bandwidthSizeCriterionCVSerial; + bool mIsAutoselectBandwidth = false; + + // IOpenmpParallelable interface + IParallelalbe::ParallelType mParallelType = IParallelalbe::ParallelType::SerialOnly; + int mOmpThreadNum = 8; + + mat mLocalPV; + mat mVariance; + mat mSDev; + cube mLoadings; + cube mScores; + cube mScoresFromKernel = cube(); // 从内核库获取的scores，用于对比 + + bool mZscore; + bool mScoresCal; + + +public: + static int treeChildCount; + +private: + //用户选择RobustGWPCA + bool mRobust; + + //用户选择GlyphPlot + bool mPlot; + +}; + +inline mat GwmGWPCATaskThread::variance() const +{ + return mVariance; +} + +inline mat GwmGWPCATaskThread::sdev() const +{ + return mSDev; +} + +inline int GwmGWPCATaskThread::parallelAbility() const +{ + return IParallelalbe::SerialOnly + #ifdef ENABLE_OpenMP + | IParallelalbe::OpenMP + #endif + ; +} + +inline IParallelalbe::ParallelType GwmGWPCATaskThread::parallelType() const +{ + return mParallelType; +} + +inline mat GwmGWPCATaskThread::localPV() const +{ + return mLocalPV; +} + +inline double GwmGWPCATaskThread::k() const +{ + return mK; +} + +inline cube GwmGWPCATaskThread::scores() const +{ + return mScores; +} + +inline cube GwmGWPCATaskThread::loadings() const +{ + return mLoadings; +} + +inline void GwmGWPCATaskThread::setK(double k) +{ + mK = k; +} + +#endif // GWMGWPCATASKTHREAD_H From 7111272a07051d188dda86e97faa7f9dfea71459 Mon Sep 17 00:00:00 2001 From: Shaoqi Yang <1647246791@qq.com> Date: Fri, 13 Mar 2026 22:43:43 +0800 Subject: [PATCH 3/4] fix calculation --- include/gwmodel.h | 3 ++- src/gwmodelpp/GWPCA.cpp | 37 +++++++++++++++++++++++++++++-------- 2 files changed, 31 insertions(+), 9 deletions(-) diff --git a/include/gwmodel.h b/include/gwmodel.h index 20966182..a1937616 100644 --- a/include/gwmodel.h +++ b/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 \ No newline at end of file +#endif // GWMODEL_H diff --git a/src/gwmodelpp/GWPCA.cpp b/src/gwmodelpp/GWPCA.cpp index ce2845b8..beae8598 100644 --- a/src/gwmodelpp/GWPCA.cpp +++ b/src/gwmodelpp/GWPCA.cpp @@ -12,7 +12,7 @@ void GWPCA::run() GWM_LOG_STAGE("Solving"); mLocalPV = pca(mX, mLoadings, mScores, mSDev); GWM_LOG_STOP_RETURN(mStatus, void()); - + mWinner = index_max(mLoadings.slice(0), 1); } @@ -20,28 +20,49 @@ 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, nDp, 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); + 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, U, 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)); - scores.slice(j).col(i) = U.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; } From 2c14e5b000a6068fa4aca24b2ccf93ee35fd97ce Mon Sep 17 00:00:00 2001 From: Shaoqi Yang <1647246791@qq.com> Date: Fri, 13 Mar 2026 23:07:06 +0800 Subject: [PATCH 4/4] remove useless test files --- test/gwmgwpcataskthread.cpp | 939 ------------------------------------ test/gwmgwpcataskthread.h | 223 --------- 2 files changed, 1162 deletions(-) delete mode 100644 test/gwmgwpcataskthread.cpp delete mode 100644 test/gwmgwpcataskthread.h diff --git a/test/gwmgwpcataskthread.cpp b/test/gwmgwpcataskthread.cpp deleted file mode 100644 index fe8438ba..00000000 --- a/test/gwmgwpcataskthread.cpp +++ /dev/null @@ -1,939 +0,0 @@ -#include "gwmgwpcataskthread.h" -#include -#include "TaskThread/gwmgeographicalweightedregressionalgorithm.h" -#include "gwmtaskthread.h" -#ifdef ENABLE_OpenMP -#include -#endif -#include -#include -#include -#include -#include -#include -#include - -int GwmGWPCATaskThread::treeChildCount = 0; - -GwmGWPCATaskThread::GwmGWPCATaskThread() : GwmSpatialMonoscaleAlgorithm(), - mAlgorithm(std::make_unique()) -{ - -} - -void GwmGWPCATaskThread::setCanceled(bool canceled) -{ - // mSelector.setCanceled(canceled); - // mSpatialWeight.distance()->setCanceled(canceled); - return GwmTaskThread::setCanceled(canceled); -} - -void GwmGWPCATaskThread::run() -{ - if(!checkCanceled()) - { - emit message(QString(tr("Setting data points ..."))); - initPoints(); - - emit message(QString(tr("Setting X and Y."))); - initXY(mX,mVariables); - - if(zscore()){ - emit message(QString(tr("Zscore normalizaiton..."))); - variableZscore(mX); - } - } - - if(mIsAutoselectBandwidth && !checkCanceled()) - { - emit message(QString(tr("Automatically selecting bandwidth ..."))); - emit tick(0, 0); - - if ((mSpatialWeight.distance()->type() == gwm::Distance::CRSDistance || - mSpatialWeight.distance()->type() == gwm::Distance::MinkwoskiDistance) && !checkCanceled()) - { - gwm::CRSDistance* d = static_cast(mSpatialWeight.distance()); - if(d) - { - d->makeParameter({ mDataPoints, mDataPoints }); - } - } - - gwm::BandwidthWeight* bandwidthWeight0 = mSpatialWeight.weight(); - if(!bandwidthWeight0) - { - qDebug() << "[GWPCA] ERROR: Cannot get bandwidth weight!"; - emit error(tr("Cannot get bandwidth weight for bandwidth selection.")); - return; - } - - double tmpMaxD = mSpatialWeight.distance()->maxDistance(); - double lower = bandwidthWeight0->adaptive() ? 2 : tmpMaxD / 5000; - double upper = bandwidthWeight0->adaptive() ? mDataPoints.n_rows : tmpMaxD; - - if(bandwidthWeight0->bandwidth() <= 0 || - (bandwidthWeight0->adaptive() && bandwidthWeight0->bandwidth() < lower) || - (!bandwidthWeight0->adaptive() && bandwidthWeight0->bandwidth() < lower)) - { - double initialBandwidth = bandwidthWeight0->adaptive() ? - std::max(2.0, std::min(20.0, (double)mDataPoints.n_rows * 0.1)) : - std::max(lower, tmpMaxD * 0.1); - bandwidthWeight0->setBandwidth(initialBandwidth); - } - - mSelector.setBandwidth(bandwidthWeight0); - mSelector.setLower(lower); - mSelector.setUpper(upper); - - try - { - gwm::BandwidthWeight* bandwidthWeight = mSelector.optimize(this); - if(bandwidthWeight && !checkCanceled()) - { - mSpatialWeight.setWeight(bandwidthWeight); - mSelector.setBandwidth(bandwidthWeight); - - gwm::BandwidthWeight* verifyBw = mSpatialWeight.weight(); - if(verifyBw && verifyBw->bandwidth() == 0) - { - qDebug() << "[GWPCA] ERROR: Bandwidth is 0 after setWeight! This may cause display issues."; - } - } - else if(!bandwidthWeight) - { - qDebug() << "[GWPCA] WARNING: Bandwidth optimization returned NULL"; - emit error(tr("Bandwidth optimization failed: no optimal bandwidth found.")); - return; - } - } - catch(const std::exception& e) - { - qDebug() << "[GWPCA] EXCEPTION during bandwidth optimization:" << e.what(); - emit error(QString(tr("Bandwidth optimization error: %1")).arg(e.what())); - return; - } - catch(...) - { - qDebug() << "[GWPCA] UNKNOWN EXCEPTION during bandwidth optimization"; - emit error(tr("Unknown error occurred during bandwidth optimization.")); - return; - } - } - - if(!checkCanceled()) - { - emit message(QString(tr("Principle components analyzing ..."))); - - try - { - if(Robust()) - { - emit message(QString(tr("Running Robust GWPCA ..."))); - mLocalPV = robustSolveSerial(mX, mLoadings, mSDev); - - if(checkCanceled()) - { - return; - } - - mVariance = mSDev % mSDev; - } - else - { - if(!mAlgorithm) - { - qDebug() << "[GWPCA] ERROR: mAlgorithm is NULL!"; - emit error(tr("GWPCA algorithm object is not initialized.")); - return; - } - - mAlgorithm->setCoords(mDataPoints); - mAlgorithm->setVariables(mX); - mAlgorithm->setSpatialWeight(mSpatialWeight); - mAlgorithm->setKeepComponents(mK); - - if(!mAlgorithm->isValid()) - { - qDebug() << "[GWPCA] ERROR: Algorithm configuration is invalid!"; - emit error(tr("GWPCA algorithm configuration is invalid.")); - return; - } - - mAlgorithm->setTelegram(std::make_unique(this)); - mAlgorithm->run(); - - auto status = mAlgorithm->status(); - if(status != gwm::Status::Success) - { - qDebug() << "[GWPCA] ERROR: Algorithm execution failed with status:" << static_cast(status); - emit error(tr("GWPCA algorithm execution failed.")); - return; - } - - mLocalPV = mAlgorithm->localPV(); - mLoadings = mAlgorithm->loadings(); - mSDev = mAlgorithm->sdev(); - mVariance = mSDev % mSDev; - - // Get scores from kernel library for comparison - if(scoresCal() && mDataPoints.n_rows <= 4096) - { - const cube& kernelScores = mAlgorithm->scores(); - mScoresFromKernel = kernelScores; - - qDebug() << "[GWPCA] Got scores from kernel library, dimensions:" << mScoresFromKernel.n_rows << "x" << mScoresFromKernel.n_cols << "x" << mScoresFromKernel.n_slices; - } - } - - if(scoresCal() && mDataPoints.n_rows <= 4096) - { - // Keep original calculateScores() result for comparison - calculateScores(); - - // Output comparison information - if(!Robust() && mScoresFromKernel.n_elem > 0) - { - qDebug() << "[GWPCA] ========== Scores Comparison Info =========="; - qDebug() << "[GWPCA] Original code scores dimensions:" << mScores.n_rows << "x" << mScores.n_cols << "x" << mScores.n_slices; - qDebug() << "[GWPCA] Kernel library scores dimensions:" << mScoresFromKernel.n_rows << "x" << mScoresFromKernel.n_cols << "x" << mScoresFromKernel.n_slices; - qDebug() << "[GWPCA] Note: Original code scores format is (nDp, mK, nDp), kernel library scores format is (nDp, nDp, mK)"; - qDebug() << "[GWPCA] Original code: mScores.slice(i) is the score matrix (nDp, mK) for point i"; - qDebug() << "[GWPCA] Kernel library: mScoresFromKernel.slice(j) is the score matrix (nDp, nDp) for component j"; - - int sampleSize = std::min(5, (int)mDataPoints.n_rows); - int sampleComps = std::min(mK, 3); - qDebug() << "[GWPCA] First" << sampleSize << "points score comparison:"; - for(int i = 0; i < sampleSize; i++) - { - qDebug() << "[GWPCA] Point" << i << "scores:"; - qDebug() << "[GWPCA] Original code (mScores.slice(" << i << ")):"; - for(int j = 0; j < sampleComps; j++) - { - vec scoreVec = mScores.slice(i).col(j); - QString scoreStr = "["; - for(uword k = 0; k < std::min((uword)5, scoreVec.n_elem); k++) - { - if(k > 0) scoreStr += ", "; - scoreStr += QString::number(scoreVec(k), 'f', 6); - } - if(scoreVec.n_elem > 5) scoreStr += ", ..."; - scoreStr += "]"; - qDebug() << "[GWPCA] Component" << (j+1) << ":" << scoreStr; - } - qDebug() << "[GWPCA] Kernel library (column" << i << "of each component):"; - for(int j = 0; j < sampleComps; j++) - { - vec scoreVec = mScoresFromKernel.slice(j).col(i); - QString scoreStr = "["; - for(uword k = 0; k < std::min((uword)5, scoreVec.n_elem); k++) - { - if(k > 0) scoreStr += ", "; - scoreStr += QString::number(scoreVec(k), 'f', 6); - } - if(scoreVec.n_elem > 5) scoreStr += ", ..."; - scoreStr += "]"; - qDebug() << "[GWPCA] Component" << (j+1) << ":" << scoreStr; - } - } - - qDebug() << "[GWPCA] Format converted comparison:"; - uword nDp = mDataPoints.n_rows; - cube kernelScoresConverted(nDp, mK, nDp, fill::zeros); - for(uword i = 0; i < nDp; i++) - { - for(int j = 0; j < mK; j++) - { - kernelScoresConverted.slice(i).col(j) = mScoresFromKernel.slice(j).col(i); - } - } - - for(int j = 0; j < sampleComps; j++) - { - double maxDiff = 0, meanDiff = 0; - uword diffCount = 0; - for(uword i = 0; i < nDp; i++) - { - vec diff = abs(mScores.slice(i).col(j) - kernelScoresConverted.slice(i).col(j)); - double localMax = diff.max(); - double localMean = mean(diff); - if(localMax > maxDiff) maxDiff = localMax; - meanDiff += localMean; - diffCount++; - } - meanDiff /= diffCount; - qDebug() << "[GWPCA] Component" << (j+1) << "score difference - Max diff:" << maxDiff << "Mean diff:" << meanDiff; - } - - qDebug() << "[GWPCA] ====================================="; - } - } - else - { - mScoresCal = false; - } - } - catch(const std::exception& e) - { - qDebug() << "[GWPCA] EXCEPTION caught:" << e.what(); - emit error(QString(tr("GWPCA execution error: %1")).arg(e.what())); - return; - } - catch(...) - { - qDebug() << "[GWPCA] UNKNOWN EXCEPTION caught"; - emit error(tr("Unknown error occurred during GWPCA execution.")); - return; - } - } - - QList win_var_PC1; - uvec iWinVar = index_max(mLoadings.slice(0), 1); - for(int i = 0; i < mDataPoints.n_rows && !checkCanceled(); i++) - { - win_var_PC1.append(mVariables.at(iWinVar(i)).name); - } - - if(!checkCanceled()) - { - CreateResultLayerData resultLayerData = { - qMakePair(QString("Comp.%1_PV"), mLocalPV), - qMakePair(QString("local_CP"), sum(mLocalPV, 1)) - }; - if(scoresCal()){ - for(int i = 0; i < mK; i++){ - resultLayerData += qMakePair(QString("Scores%1"), mScores.slice(i)); - } - }; - createResultLayer(resultLayerData,win_var_PC1); - - if(getPlot()){ - CreatePlotLayerData plotLayerData = {}; - vec vecLoadings(mDataPoints.n_rows * mVariables.size()); - QList var_pc; - int v = 0; - for(int i = 0; i < mDataPoints.n_rows; i++){ - for(int j=0; j < mVariables.size(); j++ ){ - vecLoadings(v) = mLoadings.slice(0)(i,j); - var_pc.append(mVariables.at(j).name); - v++; - } - } - plotLayerData += qMakePair(QStringLiteral("loadings"),vecLoadings); - createPlotLayer(plotLayerData, var_pc); - } - emit success(); - emit tick(100, 100); - } - if(checkCanceled()) - { - return; - } -} - -bool GwmGWPCATaskThread::isValid() -{ - gwm::BandwidthWeight* bandwidth = static_cast(mSpatialWeight.weight()); - if(bandwidth){ - if(!mIsAutoselectBandwidth) - { - if(bandwidth->adaptive()){ - if(bandwidth->bandwidth() <= mVariables.size()){ - return false; - } - } - } - if(mVariables.size() == 0){ - return false; - } - if(k()<=0 || k() > mVariables.size()){ - return false; - } - }else{ - return false; - } - return true; -} - -void GwmGWPCATaskThread::initPoints() -{ - if(!mDataLayer) - { - qDebug() << "[GWPCA::initPoints] ERROR: Data layer is NULL!"; - return; - } - - int nDp = mDataLayer->featureCount(); - mDataPoints = mat(nDp, 2, fill::zeros); - - QgsFeatureIterator iterator = mDataLayer->getFeatures(); - QgsFeature f; - for (int i = 0; iterator.nextFeature(f); i++) - { - QgsPointXY centroPoint = f.geometry().centroid().asPoint(); - mDataPoints(i, 0) = centroPoint.x(); - mDataPoints(i, 1) = centroPoint.y(); - } -} - -void GwmGWPCATaskThread::initXY(mat &x, const QList &indepVars) -{ - int nDp = mDataLayer->featureCount(), nVar = indepVars.size(); - - x = mat(nDp, nVar, fill::zeros); - - QgsFeatureIterator iterator = mDataLayer->getFeatures(); - QgsFeature f; - bool ok = false; - int errorCount = 0; - for (int i = 0; iterator.nextFeature(f); i++) - { - for (int k = 0; k < indepVars.size(); k++) - { - double vX = f.attribute(indepVars[k].name).toDouble(&ok); - if (ok) - { - x(i, k) = vX; - } - else - { - errorCount++; - if(errorCount <= 5) - { - qDebug() << "[GWPCA::initXY] WARNING: Cannot convert variable" << indepVars[k].name << "at row" << i << "to number"; - } - emit error(tr("Independent variable value cannot convert to a number. Set to 0.")); - } - } - } -} - -void GwmGWPCATaskThread::variableZscore(mat& x) -{ - mat xmean = mean(x); - mat xstd = stddev(x); - for (int k = 0; k < x.n_cols; k++) - { - x.col(k) = (x.col(k) - xmean(k))/xstd(k); - } -} - - -void GwmGWPCATaskThread::calculateScores() -{ - uword nDp = mDataPoints.n_rows, nVar = mVariables.size(); - mScores = cube(nDp, mK, nDp, fill::zeros); - - for(uword i = 0; i < nDp && !checkCanceled(); i++) - { - vec wt = mSpatialWeight.weightVector(i); - uvec positive = find(wt > 0); - vec newWt = wt.elem(positive); - mat newX = mX.rows(positive); - - if(newWt.n_rows <= 5) - { - break; - } - - mat V; - vec d; - if(!Robust()){ - wpca(newX, newWt, V, d); - }else{ - rwpca(newX, newWt, V, d); - } - - mat scorei(nDp, mK, fill::zeros); - for(int j = 0; j < mK && !checkCanceled(); j++) - { - mat score = newX.each_row() % trans(V.col(j)); - scorei.col(j) = sum(score, 1); - } - mScores.slice(i) = scorei; - emit tick(i + 1, nDp); - } -} - -void GwmGWPCATaskThread::wpca(const mat &x, const vec &wt, mat &V, vec &S) -{ - - mat xw = x.each_col() % wt; - mat centerized = (x.each_row() - sum(xw) / sum(wt)).each_col() % sqrt(wt); - //SVD - mat U; - svd(U,S,V,centerized); - //S即为R中的d - //V即为R中的v -} - -void GwmGWPCATaskThread::rwpca(const mat &x, const vec &wt, mat &V, vec &S) -{ - - mat mids = x; - mids = mids.each_row() - x.row((abs(wt - 0.5)).index_min()); - - mat score; - vec tsquared; - princomp(V, score, S, tsquared, mids.each_col() % wt); - -} - -mat GwmGWPCATaskThread::robustSolveSerial(const mat& x, cube& loadings, mat& sdev) -{ - int nDp = mDataPoints.n_rows, nVar = mX.n_cols; - - mat d_all(nVar, nDp, fill::zeros); - - loadings = cube(nDp, nVar, mK, fill::zeros); - - for(int i=0;idistance(i); - vec wt = mSpatialWeight.weightVector(i); - - uvec positive = find(wt > 0); - vec newWt = wt.elem(positive); - mat newX = x.rows(positive); - if(newWt.n_rows<=5) - { - break; - } - - mat V; - vec d; - rwpca(newX,newWt,V,d); - - mLatestWt = newWt; - d_all.col(i) = d; - - for(int j = 0; j < mK && !checkCanceled(); j++) - { - loadings.slice(j).row(i) = trans(V.col(j)); - } - emit tick(i, nDp); - } - - d_all = trans(d_all); - mat variance = (d_all / pow(sum(mLatestWt),0.5)) % (d_all / pow(sum(mLatestWt),0.5)); - - sdev = sqrt(variance); - - mat pv = variance.cols(0, mK - 1).each_col() % (1.0 / sum(variance,1)) * 100.0; - return pv; -} - -gwm::Status GwmGWPCATaskThread::getCriterion(gwm::BandwidthWeight* weight, double& criterion) -{ - if(checkCanceled()) - { - criterion = DBL_MAX; - return gwm::Status::Terminated; - } - - // 将内核库的BandwidthWeight转换为本地的GwmBandwidthWeight - GwmBandwidthWeight::KernelFunctionType kernelType = static_cast(weight->kernel()); - GwmBandwidthWeight localWeight(weight->bandwidth(), weight->adaptive(), kernelType); - - try - { - criterion = (this->*mBandwidthSelectCriterionFunction)(&localWeight); - return gwm::Status::Success; - } - catch(const std::exception& e) - { - qDebug() << "[GWPCA::getCriterion] EXCEPTION:" << e.what(); - criterion = DBL_MAX; - return gwm::Status::Success; - } - catch(...) - { - qDebug() << "[GWPCA::getCriterion] UNKNOWN EXCEPTION"; - criterion = DBL_MAX; - return gwm::Status::Success; - } -} - -void GwmGWPCATaskThread::setBandwidthSelectionCriterionType(const GwmGWPCATaskThread::BandwidthSelectionCriterionType &bandwidthSelectionCriterionType) -{ - mBandwidthSelectionCriterionType = bandwidthSelectionCriterionType; - QMap, BandwidthSelectCriterionFunction> mapper = { - std::make_pair(qMakePair(BandwidthSelectionCriterionType::CV, IParallelalbe::ParallelType::SerialOnly), &GwmGWPCATaskThread::bandwidthSizeCriterionCVSerial), - #ifdef ENABLE_OpenMP - std::make_pair(qMakePair(BandwidthSelectionCriterionType::CV, IParallelalbe::ParallelType::OpenMP), &GwmGWPCATaskThread::bandwidthSizeCriterionCVOmp), - #endif - //std::make_pair(qMakePair(BandwidthSelectionCriterionType::CV, IParallelalbe::ParallelType::CUDA), &GwmGWPCATaskThread::bandwidthSizeCriterionCVCuda), - //std::make_pair(qMakePair(BandwidthSelectionCriterionType::AIC, IParallelalbe::ParallelType::SerialOnly), &GwmGWPCATaskThread::bandwidthSizeCriterionAICSerial), - //std::make_pair(qMakePair(BandwidthSelectionCriterionType::AIC, IParallelalbe::ParallelType::OpenMP), &GwmGWPCATaskThread::bandwidthSizeCriterionAICOmp), - //std::make_pair(qMakePair(BandwidthSelectionCriterionType::AIC, IParallelalbe::ParallelType::CUDA), &GwmGWPCATaskThread::bandwidthSizeCriterionAICCuda) - }; - mBandwidthSelectCriterionFunction = mapper[qMakePair(bandwidthSelectionCriterionType, mParallelType)]; -} - -bool GwmGWPCATaskThread::isAutoselectBandwidth() const -{ - return mIsAutoselectBandwidth; -} - -void GwmGWPCATaskThread::setIsAutoselectBandwidth(bool isAutoselectBandwidth) -{ - mIsAutoselectBandwidth = isAutoselectBandwidth; -} - -void GwmGWPCATaskThread::setVariables(const QList &variables) -{ - mVariables = variables; -} - -void GwmGWPCATaskThread::setParallelType(const IParallelalbe::ParallelType &type) -{ - if(type & parallelAbility()) - { - mParallelType = type; - setBandwidthSelectionCriterionType(mBandwidthSelectionCriterionType); - - } -} - -void GwmGWPCATaskThread::createResultLayer(CreateResultLayerData data, QList winvar) -{ - QgsVectorLayer* srcLayer = mDataLayer; - QString layerFileName = QgsWkbTypes::displayString(srcLayer->wkbType()) + QStringLiteral("?"); - QString layerName = srcLayer->name(); - - - if(treeChildCount > 0) - { - if(!Robust()) layerName += QStringLiteral("_GWPCA") + "(" + QString::number(treeChildCount) + ")"; - else if(Robust()) layerName += QStringLiteral("_RGWPCA") + "(" + QString::number(treeChildCount) + ")"; - } else - { - if(!Robust()) layerName += QStringLiteral("_GWPCA"); - else if(Robust()) layerName += QStringLiteral("_RGWPCA"); - } - - treeChildCount++ ; - - -// if(!Robust()){ -// layerName += QStringLiteral("_GWPCA"); -// }else{ -// layerName += QStringLiteral("_RGWPCA"); -// } - - mResultLayer = new QgsVectorLayer(layerFileName, layerName, QStringLiteral("memory")); - mResultLayer->setCrs(srcLayer->crs()); - - - QgsFields fields; - for (QPair item : data) - { - QString title = item.first; - const mat& value = item.second; - if (value.n_cols > 1) - { - for (int k = 0; k < value.n_cols; k++) - { - QString fieldName = title.arg(k+1); - fields.append(QgsField(fieldName, QVariant::Double, QStringLiteral("double"))); - } - } - else - { - fields.append(QgsField(title, QVariant::Double, QStringLiteral("double"))); - } - } - fields.append(QgsField("win_var_PC1",QVariant::String,QStringLiteral("varchar"),255)); - mResultLayer->dataProvider()->addAttributes(fields.toList()); - mResultLayer->updateFields(); - - - mResultLayer->startEditing(); - QgsFeatureIterator iterator = srcLayer->getFeatures(); - QgsFeature f; - for (int i = 0; iterator.nextFeature(f); i++) - { - QgsFeature feature(fields); - feature.setGeometry(f.geometry()); - - int k = 0; - for (QPair item : data) - { - for (uword d = 0; d < item.second.n_cols; d++) - { - feature.setAttribute(k, item.second(i, d)); - k++; - } - } - feature.setAttribute("win_var_PC1",winvar[i]); - - mResultLayer->addFeature(feature); - } - mResultLayer->commitChanges(); -} - -void GwmGWPCATaskThread::createPlotLayer(CreatePlotLayerData data, QList varpc) -{ - QgsVectorLayer* srcLayer = mDataLayer; - QString layerFileName = QgsWkbTypes::displayString(srcLayer->wkbType()) + QStringLiteral("?"); - QString layerName = srcLayer->name(); - - - if(treeChildCount > 0) - { - layerName += QStringLiteral("_GlyphPlot") + "(" + QString::number(treeChildCount) + ")"; - } else - { - layerName += QStringLiteral("_GlyphPlot"); - } - - treeChildCount++ ; - - - if(!Robust()){ - layerName += QStringLiteral("_GWPCA"); - }else{ - layerName += QStringLiteral("_RGWPCA"); - } - - mPlotLayer = new QgsVectorLayer(QgsWkbTypes::displayString(QgsWkbTypes::LineString), layerName, QStringLiteral("memory")); - mPlotLayer->setCrs(srcLayer->crs()); - - QgsFields fields; -// for (QPair item : data) -// { -// QString title = item.first; -// const mat& value = item.second; -// if (value.n_cols > 1) -// { -// for (int k = 0; k < value.n_cols; k++) -// { -// QString fieldName = title.arg(k+1); -// fields.append(QgsField(fieldName, QVariant::Double, QStringLiteral("double"))); -// } -// } -// else -// { -// fields.append(QgsField(title, QVariant::Double, QStringLiteral("double"))); -// } -// } - fields.append(QgsField(QStringLiteral("loadings"), QVariant::Double, QStringLiteral("double"))); - fields.append(QgsField(QStringLiteral("var_name"),QVariant::String,QStringLiteral("varchar"),255)); - mPlotLayer->dataProvider()->addAttributes(fields.toList()); - mPlotLayer->updateFields(); - - mPlotLayer->startEditing(); - QgsFeatureIterator iterator = srcLayer->getFeatures(); - QgsFeature f; - double PI = 3.1415926535898; - int k = 0; - for (int i = 0; i item : data) -// { -// for (uword d = 0; d < item.second.n_cols; d++) -// { -// feature.setAttribute(k, item.second(i, d)); -// k++; -// } -// } - - for (QPair item : data) - { - for (uword d = 0; d < item.second.n_cols; d++) - { - - double x = sttPoint.x()+fabs(item.second(i, d)) * cos(k * angle)*10; - double y = sttPoint.y()+fabs(item.second(i, d)) * sin(k * angle)*10; -// std::cout< lineData ={}; - lineData+= sttPoint; - lineData+= QgsPointXY(x,y); -// std::cout<<"x: "<addFeature(feature); - mPlotLayer->addFeature(lineFeature); - } - mPlotLayer->commitChanges(); - qWarning("test %s", (QgsWkbTypes::displayString(mPlotLayer->wkbType())).toStdString().data()); -} - -double GwmGWPCATaskThread::bandwidthSizeCriterionCVSerial(GwmBandwidthWeight *weight) -{ - int mBandwidthCounter = 0; - int n = mX.n_rows; - int m = mX.n_cols; - double score = 0; - - if (mSpatialWeight.distance()->type() == gwm::Distance::CRSDistance || - mSpatialWeight.distance()->type() == gwm::Distance::MinkwoskiDistance) - { - gwm::CRSDistance* d = static_cast(mSpatialWeight.distance()); - if(d) - { - d->makeParameter({ mDataPoints, mDataPoints }); - } - } - - for (int i = 0; i < n && !checkCanceled(); i++) - { - vec distvi = mSpatialWeight.distance()->distance(i); - if(distvi.n_elem == 0) - { - qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] ERROR: Empty distance vector at point" << i; - score = DBL_MAX; - break; - } - - vec wt = weight->weight(distvi); - wt(i) = 0; - - uvec positive = find(wt > 0); - if(positive.n_elem == 0) - { - qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] WARNING: No positive weights at point" << i; - score = DBL_MAX; - break; - } - - vec newWt = wt.elem(positive); - mat newX = mX.rows(positive); - //判断length(newWt) - if(newWt.n_rows <= 1) - { - qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] WARNING: Insufficient points at" << i << ", n_rows:" << newWt.n_rows; - score = DBL_MAX; - break; - } - - mat V; - vec S; - try - { - wpca(newX, newWt, V, S); - if(V.n_cols < mK) - { - qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] WARNING: V.n_cols < mK at point" << i; - score = DBL_MAX; - break; - } - V = V.cols(0, mK - 1); - V = V * trans(V); - score = score + pow(sum(mX.row(i) - mX.row(i) * V),2); - } - catch(const std::exception& e) - { - qDebug() << "[GWPCA::bandwidthSizeCriterionCVSerial] EXCEPTION at point" << i << ":" << e.what(); - score = DBL_MAX; - break; - } - - mBandwidthCounter++; - if (mBandwidthCounter < 10) - emit tick(mBandwidthCounter * 10 + i * 5 / n, 100); - } - return score; -} -#ifdef ENABLE_OpenMP -double GwmGWPCATaskThread::bandwidthSizeCriterionCVOmp(GwmBandwidthWeight *weight) -{ - int n = mX.n_rows; - int m = mX.n_cols; - double score = 0; - bool flag = true; - vec score_all(mOmpThreadNum, fill::zeros); - int current = 0; -#pragma omp parallel for num_threads(mOmpThreadNum) - for (int i = 0; i < n; i++) - { - if(flag && !checkCanceled()) - { - int thread = omp_get_thread_num(); - vec distvi = mSpatialWeight.distance()->distance(i); - vec wt = weight->weight(distvi); - wt(i) = 0; - uvec positive = find(wt > 0); - vec newWt = wt.elem(positive); - mat newX = mX.rows(positive); - //判断length(newWt) - if(newWt.n_rows <=1) - { - flag=false; - }else{ - mat V; - vec S; - // 带宽选择时总是使用普通GWPCA（wpca），无论是否Robust模式 - // 这样Robust GWPCA会使用与普通GWPCA相同的带宽 - wpca(newX, newWt, V, S); - V = V.cols(0, mK - 1); - V = V * trans(V); - score_all(thread) += pow(sum(mX.row(i) - mX.row(i) * V),2); - } - if(mSelector.counter<10) - emit tick(mSelector.counter * 10 + current * 10 / n, 100); - current++; - } - } - score = sum(score_all); - return score; -} -#endif -bool GwmGWPCATaskThread::zscore() const -{ - return mZscore; -} - -void GwmGWPCATaskThread::setZscore(bool zscore) -{ - mZscore = zscore; -} - -bool GwmGWPCATaskThread::scoresCal() const -{ - return mScoresCal; -} - -void GwmGWPCATaskThread::setScoresCal(bool scoresCal) -{ - mScoresCal = scoresCal; -} - -bool GwmGWPCATaskThread::Robust() const -{ - return mRobust; -} - -void GwmGWPCATaskThread::setRobust(bool robust) -{ - mRobust=robust; -} - -bool GwmGWPCATaskThread::getPlot() const -{ - return mPlot; -} - -void GwmGWPCATaskThread::setPlot(bool plot) -{ - mPlot=plot; -} - diff --git a/test/gwmgwpcataskthread.h b/test/gwmgwpcataskthread.h deleted file mode 100644 index 3cda47ac..00000000 --- a/test/gwmgwpcataskthread.h +++ /dev/null @@ -1,223 +0,0 @@ -#ifndef GWMGWPCATASKTHREAD_H -#define GWMGWPCATASKTHREAD_H - -#include -#include "TaskThread/gwmspatialmonoscalealgorithm.h" -#include "TaskThread/imultivariableanalysis.h" -#include "TaskThread/iparallelable.h" - -#include "TaskThread/gwmbandwidthsizeselector.h" -#include - -class GwmGWPCATaskThread : public GwmSpatialMonoscaleAlgorithm, public IBandwidthSizeSelectable, public IGwmMultivariableAnalysis, public IOpenmpParallelable, public gwm::IBandwidthSelectable -{ - Q_OBJECT - enum BandwidthSelectionCriterionType - { - CV - }; - - typedef QList > CreateResultLayerData; - typedef QList > CreatePlotLayerData; - - typedef double (GwmGWPCATaskThread::*BandwidthSelectCriterionFunction)(GwmBandwidthWeight*); - -public: - GwmGWPCATaskThread(); - -public: - QList variables() const override - { - return QList(); - }; - void setVariables(const QList &variables) override; - void setVariables(const QList &&variables) - { - mVariables = variables; - }; - - bool isAutoselectBandwidth() const; - void setIsAutoselectBandwidth(bool isAutoselectBandwidth); - BandwidthSelectionCriterionType bandwidthSelectionCriterionType() const; - void setBandwidthSelectionCriterionType(const BandwidthSelectionCriterionType &bandwidthSelectionCriterionType); - - void setOmpThreadNum(const int threadNum) override - { - mOmpThreadNum = threadNum; - }; - - BandwidthCriterionList bandwidthSelectorCriterions() const - { - return mSelector.bandwidthCriterion(); - } - - double k() const; - void setK(double k); - - mat localPV() const; - - mat variance() const; - - mat sdev() const; - - cube loadings() const; - - cube scores() const; - -public: // QThread interface - void run() override; - - -public: // GwmTaskThread interface - QString name() const override { return tr("GWPCA"); } - - -public: // GwmSpatialMonoscaleAlgorithm interface - bool isValid() override; - - -public: // IParallelalbe interface - int parallelAbility() const override; - virtual ParallelType parallelType() const override; - virtual void setParallelType(const ParallelType& type) override; - -public: // IOpenmpParallelable interface - void setThreadNum(const int threadNum){}; - //void setOmpThreadNum(const int threadNum){}; - - bool zscore() const; - void setZscore(bool zscore); - - bool scoresCal() const; - void setScoresCal(bool scoresCal); - - bool Robust() const; - void setRobust(bool robust); - - bool getPlot() const; - void setPlot(bool plot); - - void setCanceled(bool canceled) override; - -private: - void initPoints(); - void initXY(mat& x, const QList& indepVars); - void variableZscore(mat& x); - void calculateScores(); - - void createResultLayer(CreateResultLayerData data,QList winvar); - void createPlotLayer(CreatePlotLayerData data, QList varpc); - std::unique_ptr mAlgorithm; - - // wpca函数保留，用于带宽选择（内核库的wpca是私有的，无法直接访问） - void wpca(const mat &x, const vec &wt, mat &V, vec &S); - - // Robust weighted PCA 函数 - void rwpca(const mat &x, const vec &wt, mat &V, vec &S); - - // Robust GWPCA 求解函数 - mat robustSolveSerial(const mat& x, cube& loadings, mat& sdev); - - double bandwidthSizeCriterionCVSerial(GwmBandwidthWeight* weight); -#ifdef ENABLE_OpenMP - double bandwidthSizeCriterionCVOmp(GwmBandwidthWeight* weight); -#endif - double criterion(GwmBandwidthWeight *weight) override - { - return (this->*mBandwidthSelectCriterionFunction)(weight); - } - -public: // gwm::IBandwidthSelectable interface - gwm::Status getCriterion(gwm::BandwidthWeight* weight, double& criterion) override; - -private: - QList mVariables; - mat mDataPoints; - - int mK = 2; - - mat mX; - vec mLatestWt; - - gwm::BandwidthSelector mSelector; - BandwidthSelectionCriterionType mBandwidthSelectionCriterionType = BandwidthSelectionCriterionType::CV; - BandwidthSelectCriterionFunction mBandwidthSelectCriterionFunction = &GwmGWPCATaskThread::bandwidthSizeCriterionCVSerial; - bool mIsAutoselectBandwidth = false; - - // IOpenmpParallelable interface - IParallelalbe::ParallelType mParallelType = IParallelalbe::ParallelType::SerialOnly; - int mOmpThreadNum = 8; - - mat mLocalPV; - mat mVariance; - mat mSDev; - cube mLoadings; - cube mScores; - cube mScoresFromKernel = cube(); // 从内核库获取的scores，用于对比 - - bool mZscore; - bool mScoresCal; - - -public: - static int treeChildCount; - -private: - //用户选择RobustGWPCA - bool mRobust; - - //用户选择GlyphPlot - bool mPlot; - -}; - -inline mat GwmGWPCATaskThread::variance() const -{ - return mVariance; -} - -inline mat GwmGWPCATaskThread::sdev() const -{ - return mSDev; -} - -inline int GwmGWPCATaskThread::parallelAbility() const -{ - return IParallelalbe::SerialOnly - #ifdef ENABLE_OpenMP - | IParallelalbe::OpenMP - #endif - ; -} - -inline IParallelalbe::ParallelType GwmGWPCATaskThread::parallelType() const -{ - return mParallelType; -} - -inline mat GwmGWPCATaskThread::localPV() const -{ - return mLocalPV; -} - -inline double GwmGWPCATaskThread::k() const -{ - return mK; -} - -inline cube GwmGWPCATaskThread::scores() const -{ - return mScores; -} - -inline cube GwmGWPCATaskThread::loadings() const -{ - return mLoadings; -} - -inline void GwmGWPCATaskThread::setK(double k) -{ - mK = k; -} - -#endif // GWMGWPCATASKTHREAD_H