From 7762b7851d6e94752756bc225db023b040751d1b Mon Sep 17 00:00:00 2001 From: OuGuangyu Date: Tue, 19 Mar 2024 17:31:58 +0800 Subject: [PATCH 1/4] fix r2_adj in basic flow --- include/gwmodelpp/GWRLocalCollinearity.h | 2 +- src/gwmodelpp/GWRLocalCollinearity.cpp | 47 ++++++++++++++++-------- test/testGWRLocalCollinearity.cpp | 8 ++-- 3 files changed, 36 insertions(+), 21 deletions(-) diff --git a/include/gwmodelpp/GWRLocalCollinearity.h b/include/gwmodelpp/GWRLocalCollinearity.h index ba0225c6..f4f32c1c 100644 --- a/include/gwmodelpp/GWRLocalCollinearity.h +++ b/include/gwmodelpp/GWRLocalCollinearity.h @@ -327,7 +327,7 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public bool mIsAutoselectBandwidth = false; double mTrS = 0; double mTrStS = 0; - arma::vec mSHat; + // arma::vec mSHat; FitCalculator mFitFunction = &GWRLocalCollinearity::fitSerial; PredictCalculator mPredictFunction = &GWRLocalCollinearity::predictSerial; diff --git a/src/gwmodelpp/GWRLocalCollinearity.cpp b/src/gwmodelpp/GWRLocalCollinearity.cpp index 785dc248..1f1e4b2a 100644 --- a/src/gwmodelpp/GWRLocalCollinearity.cpp +++ b/src/gwmodelpp/GWRLocalCollinearity.cpp @@ -16,14 +16,17 @@ RegressionDiagnostic GWRLocalCollinearity::CalcDiagnostic(const mat& x, const ve vec r = y - sum(betas % x, 1); double rss = sum(r % r); double n = (double)x.n_rows; - double AIC = n * log(rss / n) + n * log(2 * datum::pi) + n + shat(0); - double AICc = n * log(rss / n) + n * log(2 * datum::pi) + n * ((n + shat(0)) / (n - 2 - shat(0))); + // double AIC = n * log(rss / n) + n * log(2 * datum::pi) + n + shat(0);//这个计算结果不知道为什么不对。 + // double AICc = n * log(rss / n) + n * log(2 * datum::pi) + n * ((n + shat(0)) / (n - 2 - shat(0))); double edf = n - 2 * shat(0) + shat(1); double enp = 2 * shat(0) - shat(1); + double s2 = rss / (n - enp); + double AIC = n * (log(2 * datum::pi * s2) + 1) + 2 * (enp + 1); + double AICc = n * (log(2 * datum::pi * s2)) + n * ((1 + enp / n) / (1 - (enp + 2) / n)); double yss = sum((y - mean(y)) % (y - mean(y))); double r2 = 1 - rss / yss; double r2_adj = 1 - (1 - r2) * (n - 1) / (edf - 1); - return { rss, AIC, AICc, enp, edf, r2, r2_adj }; + return {rss, AIC, AICc, enp, edf, r2, r2_adj}; } GWRLocalCollinearity::GWRLocalCollinearity() @@ -76,14 +79,18 @@ mat GWRLocalCollinearity::fit() GWM_LOG_STAGE("Model Diagnostic"); vec mYHat = sum(mBetas % mX,1); vec mResidual = mY - mYHat; - mDiagnostic.RSS = sum(mResidual % mResidual); - mDiagnostic.ENP = 2*this->mTrS - this->mTrStS; - mDiagnostic.EDF = nDp - mDiagnostic.ENP; - double s2 = mDiagnostic.RSS / (nDp - mDiagnostic.ENP); - mDiagnostic.AIC = nDp * (log(2*M_PI*s2)+1) + 2*(mDiagnostic.ENP + 1); - mDiagnostic.AICc = nDp * (log(2*M_PI*s2)) + nDp*( (1+mDiagnostic.ENP/nDp) / (1-(mDiagnostic.ENP+2)/nDp) ); - mDiagnostic.RSquare = 1 - mDiagnostic.RSS/sum((mY - mean(mY)) % (mY - mean(mY))); - mDiagnostic.RSquareAdjust = 1 - (1 - mDiagnostic.RSquare)*(nDp - 1) / (mDiagnostic.EDF); + vec shat = vec(2, fill::zeros); + shat(0)=mTrS; + shat(1)=mTrStS; + mDiagnostic = CalcDiagnostic(mX, mY, mBetas, shat); + // mDiagnostic.RSS = sum(mResidual % mResidual); + // mDiagnostic.ENP = 2*this->mTrS - this->mTrStS; + // mDiagnostic.EDF = nDp - mDiagnostic.ENP; + // double s2 = mDiagnostic.RSS / (nDp - mDiagnostic.ENP); + // mDiagnostic.AIC = nDp * (log(2*M_PI*s2)+1) + 2*(mDiagnostic.ENP + 1); + // mDiagnostic.AICc = nDp * (log(2*M_PI*s2)) + nDp*( (1+mDiagnostic.ENP/nDp) / (1-(mDiagnostic.ENP+2)/nDp) ); + // mDiagnostic.RSquare = 1 - mDiagnostic.RSS/sum((mY - mean(mY)) % (mY - mean(mY))); + // mDiagnostic.RSquareAdjust = 1 - (1 - mDiagnostic.RSquare)*(nDp - 1) / (mDiagnostic.EDF);//这里少了一个-1 return mBetas; } @@ -334,6 +341,7 @@ mat GWRLocalCollinearity::fitSerial(const mat& x, const vec& y) vec localcn(nDp, fill::zeros); vec locallambda(nDp, fill::zeros); vec hatrow(nDp, fill::zeros); + vec shat = vec(2, fill::zeros); for(uword i=0;imTrS += hatrow(i); - this->mTrStS += sum(hatrow % hatrow); + try{ + mat xtwxinv = inv(xtwx); + rowvec hatrow = x1w.row(i) * xtwxinv * trans(x1w); + this->mTrS += hatrow(i); + this->mTrStS += sum(hatrow % hatrow); + } + catch (const exception& e) + { + GWM_LOG_ERROR(e.what()); + throw e; + } GWM_LOG_PROGRESS(i + 1, nDp); } return betas; diff --git a/test/testGWRLocalCollinearity.cpp b/test/testGWRLocalCollinearity.cpp index 5662cce7..cfa0b9d6 100644 --- a/test/testGWRLocalCollinearity.cpp +++ b/test/testGWRLocalCollinearity.cpp @@ -45,10 +45,10 @@ TEST_CASE("LocalCollinearityGWR: basic flow") REQUIRE_NOTHROW(algorithm.fit()); RegressionDiagnostic diagnostic = algorithm.diagnostic(); - /*REQUIRE_THAT(diagnostic.AIC, Catch::MatchersWithinAbs(2461.565456, 1e-6)); - REQUIRE_THAT(diagnostic.AICc, Catch::MatchersWithinAbs(2464.600255, 1e-6)); - REQUIRE_THAT(diagnostic.RSquare, Catch::MatchersWithinAbs(0.708010632044736, 1e-6)); - REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::MatchersWithinAbs(0.674975341723766, 1e-6));*/ + REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2461.5654565, 1e-8)); + REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2464.60025589, 1e-8)); + REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.708010632043, 1e-8)); + REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.674975341722, 1e-8)); REQUIRE(algorithm.hasIntercept() == true); } From 783384f070c3928a80261ba9cf760d4dd3b64c6c Mon Sep 17 00:00:00 2001 From: OuGuangyu Date: Tue, 19 Mar 2024 20:21:59 +0800 Subject: [PATCH 2/4] update test, all tests passed rebase master --- include/gwmodelpp/GWRLocalCollinearity.h | 25 ++-- src/gwmodelpp/GWRLocalCollinearity.cpp | 28 ++-- test/testGWRLocalCollinearity.cpp | 180 +++++++++++++---------- 3 files changed, 131 insertions(+), 102 deletions(-) diff --git a/include/gwmodelpp/GWRLocalCollinearity.h b/include/gwmodelpp/GWRLocalCollinearity.h index f4f32c1c..4621c53b 100644 --- a/include/gwmodelpp/GWRLocalCollinearity.h +++ b/include/gwmodelpp/GWRLocalCollinearity.h @@ -276,18 +276,19 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public double mBandwidthLastCriterion = DBL_MAX; //!< \~english Last criterion for bandwidth selection. \~chinese 上一次带宽优选的有效指标值。 - /** - * @brief \~english Get the CV. \~chinese 返回cv的函数 - * - * @param bw \~english \~chinese - * @param kernel \~english \~chinese - * @param adaptive \~english \~chinese - * @param lambda \~english \~chinese - * @param lambdaAdjust \~english \~chinese - * @param cnThresh \~english \~chinese - * @return double \~english \~chinese - */ - double LcrCV(double bw,arma::uword kernel, bool adaptive,double lambda,bool lambdaAdjust,double cnThresh); + // //这个函数并没有完成，计算local cv的话，直接放在了循环里面算，所以需求似乎也不大 + // /** + // * @brief \~english Get the CV. \~chinese 返回cv的函数 + // * + // * @param bw \~english \~chinese + // * @param kernel \~english \~chinese + // * @param adaptive \~english \~chinese + // * @param lambda \~english \~chinese + // * @param lambdaAdjust \~english \~chinese + // * @param cnThresh \~english \~chinese + // * @return double \~english \~chinese + // */ + // double LcrCV(double bw,arma::uword kernel, bool adaptive,double lambda,bool lambdaAdjust,double cnThresh); /** * @brief \~english Ridge linear regression. \~chinese 岭回归。 diff --git a/src/gwmodelpp/GWRLocalCollinearity.cpp b/src/gwmodelpp/GWRLocalCollinearity.cpp index 1f1e4b2a..c3a345f1 100644 --- a/src/gwmodelpp/GWRLocalCollinearity.cpp +++ b/src/gwmodelpp/GWRLocalCollinearity.cpp @@ -83,14 +83,6 @@ mat GWRLocalCollinearity::fit() shat(0)=mTrS; shat(1)=mTrStS; mDiagnostic = CalcDiagnostic(mX, mY, mBetas, shat); - // mDiagnostic.RSS = sum(mResidual % mResidual); - // mDiagnostic.ENP = 2*this->mTrS - this->mTrStS; - // mDiagnostic.EDF = nDp - mDiagnostic.ENP; - // double s2 = mDiagnostic.RSS / (nDp - mDiagnostic.ENP); - // mDiagnostic.AIC = nDp * (log(2*M_PI*s2)+1) + 2*(mDiagnostic.ENP + 1); - // mDiagnostic.AICc = nDp * (log(2*M_PI*s2)) + nDp*( (1+mDiagnostic.ENP/nDp) / (1-(mDiagnostic.ENP+2)/nDp) ); - // mDiagnostic.RSquare = 1 - mDiagnostic.RSS/sum((mY - mean(mY)) % (mY - mean(mY))); - // mDiagnostic.RSquareAdjust = 1 - (1 - mDiagnostic.RSquare)*(nDp - 1) / (mDiagnostic.EDF);//这里少了一个-1 return mBetas; } @@ -471,15 +463,21 @@ mat GWRLocalCollinearity::fitOmp(const mat& x, const vec& y) } betas.row(i) = trans( ridgelm(wi,locallambda(i)) ); //如果没有给regressionpoint - mat xm = x; mat xtw = trans(x % (wi * wispan1)); mat xtwx = xtw * x; - mat xtwxinv = inv(xtwx); - rowvec hatrow = x1w.row(i) * xtwxinv * trans(x1w); - shat_all(0, thread) += hatrow(i); - shat_all(1, thread) += sum(hatrow % hatrow); - this->mTrS += hatrow(i); - this->mTrStS += sum(hatrow % hatrow); + try{ + mat xtwxinv = inv(xtwx); + rowvec hatrow = x1w.row(i) * xtwxinv * trans(x1w); + shat_all(0, thread) += hatrow(i); + shat_all(1, thread) += sum(hatrow % hatrow); + this->mTrS += hatrow(i); + this->mTrStS += sum(hatrow % hatrow); + } + catch (const exception& e) + { + GWM_LOG_ERROR(e.what()); + throw e; + } GWM_LOG_PROGRESS(i + 1, nDp); } vec shat = sum(shat_all,1); diff --git a/test/testGWRLocalCollinearity.cpp b/test/testGWRLocalCollinearity.cpp index cfa0b9d6..658b9238 100644 --- a/test/testGWRLocalCollinearity.cpp +++ b/test/testGWRLocalCollinearity.cpp @@ -19,7 +19,7 @@ using namespace std; using namespace arma; using namespace gwm; -TEST_CASE("LocalCollinearityGWR: basic flow") +TEST_CASE("LocalCollinearityGWR") { mat londonhp100_coord, londonhp100_data; vector londonhp100_fields; @@ -28,96 +28,123 @@ TEST_CASE("LocalCollinearityGWR: basic flow") FAIL("Cannot load londonhp100 data."); } - CRSDistance distance(false); - BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); - SpatialWeight spatial(&bandwidth, &distance); + SECTION("adaptive bandwidth | no bandwidth optimization | no lambda adjust") + { - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); + CRSDistance distance(false); + BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); + SpatialWeight spatial(&bandwidth, &distance); - GWRLocalCollinearity algorithm; - algorithm.setCoords(londonhp100_coord); - algorithm.setDependentVariable(y); - algorithm.setIndependentVariables(x); - algorithm.setSpatialWeight(spatial); - algorithm.setHasHatMatrix(true); - //algorithm.setLambdaAdjust(true); - REQUIRE_NOTHROW(algorithm.fit()); + vec y = londonhp100_data.col(0); + mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); - RegressionDiagnostic diagnostic = algorithm.diagnostic(); - REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2461.5654565, 1e-8)); - REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2464.60025589, 1e-8)); - REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.708010632043, 1e-8)); - REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.674975341722, 1e-8)); + GWRLocalCollinearity algorithm; + algorithm.setCoords(londonhp100_coord); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setHasHatMatrix(true); - REQUIRE(algorithm.hasIntercept() == true); -} + REQUIRE_NOTHROW(algorithm.fit()); -TEST_CASE("LocalCollinearityGWR: adaptive bandwidth autoselection of with CV") -{ - mat londonhp100_coord, londonhp100_data; - vector londonhp100_fields; - if (!read_londonhp100(londonhp100_coord, londonhp100_data, londonhp100_fields)) - { - FAIL("Cannot load londonhp100 data."); + RegressionDiagnostic diagnostic = algorithm.diagnostic(); + REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2461.5654565, 1e-8)); + REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2464.60025589, 1e-8)); + REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.708010632043, 1e-8)); + REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.674975341722, 1e-8)); + + REQUIRE(algorithm.hasIntercept() == true); } - CRSDistance distance(false); - BandwidthWeight bandwidth(0, true, BandwidthWeight::Gaussian); - SpatialWeight spatial(&bandwidth, &distance); + SECTION("adaptive bandwidth | bandwidth optimization CV | no lambda adjust") + { + CRSDistance distance(false); + BandwidthWeight bandwidth(0, true, BandwidthWeight::Gaussian); + SpatialWeight spatial(&bandwidth, &distance); - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); + vec y = londonhp100_data.col(0); + mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); - GWRLocalCollinearity algorithm; - algorithm.setCoords(londonhp100_coord); - algorithm.setDependentVariable(y); - algorithm.setIndependentVariables(x); - algorithm.setSpatialWeight(spatial); - algorithm.setHasHatMatrix(true); + GWRLocalCollinearity algorithm; + algorithm.setCoords(londonhp100_coord); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setHasHatMatrix(true); - algorithm.setIsAutoselectBandwidth(true); - algorithm.setBandwidthSelectionCriterion(GWRLocalCollinearity::BandwidthSelectionCriterionType::CV); - - REQUIRE_NOTHROW(algorithm.fit()); + algorithm.setIsAutoselectBandwidth(true); + algorithm.setBandwidthSelectionCriterion(GWRLocalCollinearity::BandwidthSelectionCriterionType::CV); - size_t bw = (size_t)algorithm.spatialWeight().weight()->bandwidth(); - REQUIRE(bw == 67); -} + REQUIRE_NOTHROW(algorithm.fit()); -TEST_CASE("LocalCollinearityGWR: ") -{ - mat londonhp100_coord, londonhp100_data; - vector londonhp100_fields; - if (!read_londonhp100(londonhp100_coord, londonhp100_data, londonhp100_fields)) - { - FAIL("Cannot load londonhp100 data."); + size_t bw = (size_t)algorithm.spatialWeight().weight()->bandwidth(); + REQUIRE(bw == 67); } - CRSDistance distance(false); - BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); - SpatialWeight spatial(&bandwidth, &distance); + SECTION("adaptive bandwidth | no bandwidth optimization | lambda adjust | CnThresh=20 ") + { - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); + CRSDistance distance(false); + BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); + SpatialWeight spatial(&bandwidth, &distance); - GWRLocalCollinearity algorithm; - algorithm.setCoords(londonhp100_coord); - algorithm.setDependentVariable(y); - algorithm.setIndependentVariables(x); - algorithm.setSpatialWeight(spatial); - algorithm.setHasHatMatrix(true); - algorithm.setLambdaAdjust(true); - REQUIRE_NOTHROW(algorithm.fit()); + vec y = londonhp100_data.col(0); + mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); - RegressionDiagnostic diagnostic = algorithm.diagnostic(); - /*REQUIRE_THAT(diagnostic.AIC, Catch::MatchersWithinAbs(2461.565456, 1e-6)); - REQUIRE_THAT(diagnostic.AICc, Catch::MatchersWithinAbs(2464.600255, 1e-6)); - REQUIRE_THAT(diagnostic.RSquare, Catch::MatchersWithinAbs(0.708010632044736, 1e-6)); - REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::MatchersWithinAbs(0.674975341723766, 1e-6));*/ + GWRLocalCollinearity algorithm; + algorithm.setCoords(londonhp100_coord); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setHasHatMatrix(true); + algorithm.setLambdaAdjust(true); + // algorithm.setLambda(0); + algorithm.setCnThresh(20); + REQUIRE_NOTHROW(algorithm.fit()); + RegressionDiagnostic diagnostic = algorithm.diagnostic(); + REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2461.8623182524, 1e-8)); + REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2464.8971176381, 1e-8)); + REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.70714253941241, 1e-8)); + REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.67400903424089, 1e-8)); + } + #ifdef ENABLE_OPENMP + SECTION("adaptive bandwidth | no bandwidth optimization | lambda adjust | CnThresh=20 ") + { + + CRSDistance distance(false); + BandwidthWeight bandwidth(0, true, BandwidthWeight::Gaussian); + SpatialWeight spatial(&bandwidth, &distance); + + vec y = londonhp100_data.col(0); + mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); + + GWRLocalCollinearity algorithm; + algorithm.setCoords(londonhp100_coord); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setHasHatMatrix(true); + algorithm.setIsAutoselectBandwidth(true); + algorithm.setBandwidthSelectionCriterion(GWRLocalCollinearity::BandwidthSelectionCriterionType::CV); + algorithm.setParallelType(ParallelType::OpenMP); + algorithm.setOmpThreadNum(6); + REQUIRE_NOTHROW(algorithm.fit()); + + + double bw = algorithm.spatialWeight().weight()->bandwidth(); + REQUIRE(bw == 67.0); + + RegressionDiagnostic diagnostic = algorithm.diagnostic(); + REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2458.2472656218, 1e-8)); + REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2459.743757379, 1e-8)); + REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.6873384732363, 1e-8)); + REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.664362879709, 1e-8)); + } + #endif } +<<<<<<< HEAD /* #ifdef ENABLE_OPENMP TEST_CASE("LocalCollinearityGWR: multithread basic flow") @@ -159,6 +186,8 @@ TEST_CASE("LocalCollinearityGWR: multithread basic flow") REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::MatchersWithinAbs(0.66436287, 1e-6)); } #endif*/ +======= +>>>>>>> 59ec57b (update test, all tests passed) TEST_CASE("LcGWR: cancel") { @@ -179,8 +208,9 @@ TEST_CASE("LcGWR: cancel") const initializer_list parallel_list = { ParallelType::SerialOnly #ifdef ENABLE_OPENMP - , ParallelType::OpenMP -#endif // ENABLE_OPENMP + , + ParallelType::OpenMP +#endif // ENABLE_OPENMP }; auto parallel = GENERATE_REF(values(parallel_list)); @@ -188,7 +218,8 @@ TEST_CASE("LcGWR: cancel") { auto stage = GENERATE(as{}, "bandwidthSize", "fit"); auto progress = GENERATE(0, 10); - INFO("Settings: " << "Parallel:" << parallel << ", Stage:" << stage << ", " << progress); + INFO("Settings: " + << "Parallel:" << parallel << ", Stage:" << stage << ", " << progress); auto telegram = make_unique(stage, progress); GWRLocalCollinearity algorithm; @@ -207,5 +238,4 @@ TEST_CASE("LcGWR: cancel") REQUIRE_NOTHROW(algorithm.fit()); REQUIRE(algorithm.status() == Status::Terminated); } - } From b666436e238087aec4cb3f538f770ddf61e2c609 Mon Sep 17 00:00:00 2001 From: OuGuangyu Date: Tue, 3 Dec 2024 21:31:19 +0800 Subject: [PATCH 3/4] rebase master, add 1 test section --- include/gwmodelpp/GWRLocalCollinearity.h | 2 +- test/testGWRLocalCollinearity.cpp | 91 ++++++++---------------- 2 files changed, 31 insertions(+), 62 deletions(-) diff --git a/include/gwmodelpp/GWRLocalCollinearity.h b/include/gwmodelpp/GWRLocalCollinearity.h index 4621c53b..1dc79e1c 100644 --- a/include/gwmodelpp/GWRLocalCollinearity.h +++ b/include/gwmodelpp/GWRLocalCollinearity.h @@ -323,7 +323,7 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public double mLambda = 0; bool mLambdaAdjust = false; - double mCnThresh = 30; + double mCnThresh = 30; //maximum value for condition number, commonly set between 20 and 30 bool mHasHatMatrix = false; bool mIsAutoselectBandwidth = false; double mTrS = 0; diff --git a/test/testGWRLocalCollinearity.cpp b/test/testGWRLocalCollinearity.cpp index 658b9238..c7421b40 100644 --- a/test/testGWRLocalCollinearity.cpp +++ b/test/testGWRLocalCollinearity.cpp @@ -28,6 +28,9 @@ TEST_CASE("LocalCollinearityGWR") FAIL("Cannot load londonhp100 data."); } + vec y = londonhp100_data.col(0); + mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); + SECTION("adaptive bandwidth | no bandwidth optimization | no lambda adjust") { @@ -35,9 +38,6 @@ TEST_CASE("LocalCollinearityGWR") BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); SpatialWeight spatial(&bandwidth, &distance); - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); - GWRLocalCollinearity algorithm; algorithm.setCoords(londonhp100_coord); algorithm.setDependentVariable(y); @@ -62,9 +62,6 @@ TEST_CASE("LocalCollinearityGWR") BandwidthWeight bandwidth(0, true, BandwidthWeight::Gaussian); SpatialWeight spatial(&bandwidth, &distance); - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); - GWRLocalCollinearity algorithm; algorithm.setCoords(londonhp100_coord); algorithm.setDependentVariable(y); @@ -88,9 +85,6 @@ TEST_CASE("LocalCollinearityGWR") BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); SpatialWeight spatial(&bandwidth, &distance); - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); - GWRLocalCollinearity algorithm; algorithm.setCoords(londonhp100_coord); algorithm.setDependentVariable(y); @@ -108,17 +102,36 @@ TEST_CASE("LocalCollinearityGWR") REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.70714253941241, 1e-8)); REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.67400903424089, 1e-8)); } - - #ifdef ENABLE_OPENMP - SECTION("adaptive bandwidth | no bandwidth optimization | lambda adjust | CnThresh=20 ") + + SECTION("adaptive bandwidth | no bandwidth optimization | lambda 0.1 ") { CRSDistance distance(false); - BandwidthWeight bandwidth(0, true, BandwidthWeight::Gaussian); + BandwidthWeight bandwidth(36, true, BandwidthWeight::Gaussian); SpatialWeight spatial(&bandwidth, &distance); - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); + GWRLocalCollinearity algorithm; + algorithm.setCoords(londonhp100_coord); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setHasHatMatrix(true); + algorithm.setLambda(0.1); + REQUIRE_NOTHROW(algorithm.fit()); + + RegressionDiagnostic diagnostic = algorithm.diagnostic(); + REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2462.0038025123, 1e-8)); + REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2465.0386018980, 1e-8)); + REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.706727898945, 1e-8)); + REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.673547481901, 1e-8)); + } + +#ifdef ENABLE_OPENMP + SECTION("adaptive bandwidth | bandwidth optimization ") + { + CRSDistance distance(false); + BandwidthWeight bandwidth(0, true, BandwidthWeight::Gaussian); + SpatialWeight spatial(&bandwidth, &distance); GWRLocalCollinearity algorithm; algorithm.setCoords(londonhp100_coord); @@ -129,10 +142,9 @@ TEST_CASE("LocalCollinearityGWR") algorithm.setIsAutoselectBandwidth(true); algorithm.setBandwidthSelectionCriterion(GWRLocalCollinearity::BandwidthSelectionCriterionType::CV); algorithm.setParallelType(ParallelType::OpenMP); - algorithm.setOmpThreadNum(6); + algorithm.setOmpThreadNum(omp_get_num_threads()); REQUIRE_NOTHROW(algorithm.fit()); - double bw = algorithm.spatialWeight().weight()->bandwidth(); REQUIRE(bw == 67.0); @@ -142,52 +154,9 @@ TEST_CASE("LocalCollinearityGWR") REQUIRE_THAT(diagnostic.RSquare, Catch::Matchers::WithinAbs(0.6873384732363, 1e-8)); REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::Matchers::WithinAbs(0.664362879709, 1e-8)); } - #endif -} -<<<<<<< HEAD -/* -#ifdef ENABLE_OPENMP -TEST_CASE("LocalCollinearityGWR: multithread basic flow") -{ - 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(0, true, BandwidthWeight::Gaussian); - SpatialWeight spatial(&bandwidth, &distance); - - vec y = londonhp100_data.col(0); - mat x = join_rows(ones(londonhp100_coord.n_rows), londonhp100_data.cols(1, 3)); +#endif - GWRLocalCollinearity algorithm; - algorithm.setCoords(londonhp100_coord); - algorithm.setDependentVariable(y); - algorithm.setIndependentVariables(x); - algorithm.setSpatialWeight(spatial); - algorithm.setHasHatMatrix(true); - algorithm.setIsAutoselectBandwidth(true); - algorithm.setBandwidthSelectionCriterion(GWRLocalCollinearity::BandwidthSelectionCriterionType::CV); - algorithm.setParallelType(ParallelType::OpenMP); - algorithm.setOmpThreadNum(omp_get_num_threads()); - REQUIRE_NOTHROW(algorithm.fit()); - - - double bw = algorithm.spatialWeight().weight()->bandwidth(); - REQUIRE(bw == 67.0); - - RegressionDiagnostic diagnostic = algorithm.diagnostic(); - REQUIRE_THAT(diagnostic.AIC, Catch::MatchersWithinAbs(2458.2472656218, 1e-6)); - REQUIRE_THAT(diagnostic.AICc, Catch::MatchersWithinAbs(2459.743757379, 1e-6)); - REQUIRE_THAT(diagnostic.RSquare, Catch::MatchersWithinAbs(0.68733847, 1e-6)); - REQUIRE_THAT(diagnostic.RSquareAdjust, Catch::MatchersWithinAbs(0.66436287, 1e-6)); } -#endif*/ -======= ->>>>>>> 59ec57b (update test, all tests passed) TEST_CASE("LcGWR: cancel") { From 8f0962e3aca5c78b4504f2f6988a54cadaa614d9 Mon Sep 17 00:00:00 2001 From: OuGuangyu Date: Wed, 8 Jan 2025 23:19:07 +0800 Subject: [PATCH 4/4] add local_cn and local_lambda with test --- include/gwmodelpp/GWRLocalCollinearity.h | 18 +++++++++++++++--- src/gwmodelpp/GWRLocalCollinearity.cpp | 16 +++++++--------- test/testGWRLocalCollinearity.cpp | 15 +++++++++++++++ 3 files changed, 37 insertions(+), 12 deletions(-) diff --git a/include/gwmodelpp/GWRLocalCollinearity.h b/include/gwmodelpp/GWRLocalCollinearity.h index 1dc79e1c..92afc0ef 100644 --- a/include/gwmodelpp/GWRLocalCollinearity.h +++ b/include/gwmodelpp/GWRLocalCollinearity.h @@ -33,7 +33,7 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public }; typedef double (GWRLocalCollinearity::*BandwidthSelectionCriterionCalculator)(BandwidthWeight*); //!< \~english Calculator to get criterion for bandwidth optimization \~chinese 带宽优选指标值计算函数 - typedef arma::mat (GWRLocalCollinearity::*FitCalculator)(const arma::mat&, const arma::vec&); //!< \~english Calculator to predict \~chinese 用于预测的函数 + typedef arma::mat (GWRLocalCollinearity::*FitCalculator)(const arma::mat&, const arma::vec&, arma::vec&, arma::vec&); //!< \~english Calculator to predict \~chinese 用于预测的函数 typedef arma::mat (GWRLocalCollinearity::*PredictCalculator)(const arma::mat&, const arma::mat&, const arma::vec&); //!< \~english Calculator to predict \~chinese 用于预测的函数 /** @@ -150,6 +150,16 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public return mDiagnostic; } + arma::vec localCN() const + { + return mLocalCN; + } + + arma::vec localLambda() const + { + return mLocalLambda; + } + /** * @brief \~english Get whether bandwidth optimization is enabled. \~chinese 获取是否进行带宽优选。 * @@ -231,7 +241,7 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public * @param x \~english Independent variables \~chinese 自变量 * @param y \~english Dependent variables \~chinese 因变量 */ - arma::mat fitSerial(const arma::mat& x, const arma::vec& y); + arma::mat fitSerial(const arma::mat& x, const arma::vec& y, arma::vec& localcn, arma::vec& locallambda); /** * @brief \~english Multithreading implementation of fitting function. \~chinese 拟合函数的多线程实现。 @@ -239,7 +249,7 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public * @param x \~english Independent variables \~chinese 自变量 * @param y \~english Dependent variables \~chinese 因变量 */ - arma::mat fitOmp(const arma::mat& x, const arma::vec& y); + arma::mat fitOmp(const arma::mat& x, const arma::vec& y, arma::vec& localcn, arma::vec& locallambda); /** * @brief \~english Non-parallel implementation of prediction function. \~chinese 预测函数的非并行实现。 @@ -329,6 +339,8 @@ class GWRLocalCollinearity : public GWRBase, public IBandwidthSelectable, public double mTrS = 0; double mTrStS = 0; // arma::vec mSHat; + arma::vec mLocalCN; + arma::vec mLocalLambda; FitCalculator mFitFunction = &GWRLocalCollinearity::fitSerial; PredictCalculator mPredictFunction = &GWRLocalCollinearity::predictSerial; diff --git a/src/gwmodelpp/GWRLocalCollinearity.cpp b/src/gwmodelpp/GWRLocalCollinearity.cpp index c3a345f1..66fc8f99 100644 --- a/src/gwmodelpp/GWRLocalCollinearity.cpp +++ b/src/gwmodelpp/GWRLocalCollinearity.cpp @@ -16,8 +16,6 @@ RegressionDiagnostic GWRLocalCollinearity::CalcDiagnostic(const mat& x, const ve vec r = y - sum(betas % x, 1); double rss = sum(r % r); double n = (double)x.n_rows; - // double AIC = n * log(rss / n) + n * log(2 * datum::pi) + n + shat(0);//这个计算结果不知道为什么不对。 - // double AICc = n * log(rss / n) + n * log(2 * datum::pi) + n * ((n + shat(0)) / (n - 2 - shat(0))); double edf = n - 2 * shat(0) + shat(1); double enp = 2 * shat(0) - shat(1); double s2 = rss / (n - enp); @@ -73,7 +71,7 @@ mat GWRLocalCollinearity::fit() vec hatrow(nDp,fill::zeros); GWM_LOG_STAGE("Model fitting"); - mBetas = (this->*mFitFunction)(mX, mY); + mBetas = (this->*mFitFunction)(mX, mY, mLocalCN, mLocalLambda); GWM_LOG_STOP_RETURN(mStatus, mat(nDp, nVar, arma::fill::zeros)); GWM_LOG_STAGE("Model Diagnostic"); @@ -326,12 +324,12 @@ double GWRLocalCollinearity::bandwidthSizeCriterionCVOmp(BandwidthWeight* bandwi } #endif -mat GWRLocalCollinearity::fitSerial(const mat& x, const vec& y) +mat GWRLocalCollinearity::fitSerial(const mat& x, const vec& y, vec& localcn, vec& locallambda) { uword nDp = mCoords.n_rows, nVar = x.n_cols; mat betas(nDp, nVar, fill::zeros); - vec localcn(nDp, fill::zeros); - vec locallambda(nDp, fill::zeros); + localcn=vec(nDp, fill::zeros); + locallambda=vec(nDp, fill::zeros); vec hatrow(nDp, fill::zeros); vec shat = vec(2, fill::zeros); for(uword i=0;i()->bandwidth(); REQUIRE(bw == 67.0); + REQUIRE_THAT(algorithm.localCN().max(), Catch::Matchers::WithinAbs(50.634464389348, 1e-8)); + REQUIRE_THAT(algorithm.localCN().min(), Catch::Matchers::WithinAbs(45.520575900277, 1e-8)); + RegressionDiagnostic diagnostic = algorithm.diagnostic(); REQUIRE_THAT(diagnostic.AIC, Catch::Matchers::WithinAbs(2458.2472656218, 1e-8)); REQUIRE_THAT(diagnostic.AICc, Catch::Matchers::WithinAbs(2459.743757379, 1e-8));