From 2ba66f8542f45d89b6f17ff015608fb1bfa0f223 Mon Sep 17 00:00:00 2001 From: HPDell Date: Sun, 23 Feb 2025 18:23:21 +0000 Subject: [PATCH 1/5] fix: min distance remove 0 --- src/gwmodelpp/spatialweight/CRSDistance.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gwmodelpp/spatialweight/CRSDistance.cpp b/src/gwmodelpp/spatialweight/CRSDistance.cpp index 0c00116f..102f5e4d 100644 --- a/src/gwmodelpp/spatialweight/CRSDistance.cpp +++ b/src/gwmodelpp/spatialweight/CRSDistance.cpp @@ -149,7 +149,8 @@ double CRSDistance::minDistance() double minD = DBL_MAX; for (uword i = 0; i < mParameter->total; i++) { - double d = min(mCalculator(mParameter->focusPoints.row(i), mParameter->dataPoints)); + vec ds = mCalculator(mParameter->focusPoints.row(i), mParameter->dataPoints); + double d = min(ds.elem(find(ds > 0.0))); minD = d < minD ? d : minD; } return minD; From ede8c62608d463b28bbe6c3b1dceb619eb85819f Mon Sep 17 00:00:00 2001 From: HPDell Date: Sun, 23 Feb 2025 21:18:36 +0000 Subject: [PATCH 2/5] edit: make mWorkRange as pair --- include/gwmodelpp/GWRBasic.h | 2 +- include/gwmodelpp/GWRMultiscale.h | 2 +- src/gwmodelpp/GWRBasic.cpp | 63 +++++++++++++------------------ src/gwmodelpp/GWRMultiscale.cpp | 59 +++++++++++++---------------- 4 files changed, 55 insertions(+), 71 deletions(-) diff --git a/include/gwmodelpp/GWRBasic.h b/include/gwmodelpp/GWRBasic.h index 30665d24..e333c491 100644 --- a/include/gwmodelpp/GWRBasic.h +++ b/include/gwmodelpp/GWRBasic.h @@ -763,7 +763,7 @@ class GWRBasic : public GWRBase, public IBandwidthSelectable, public IVarialbeSe int mWorkerId = 0; int mWorkerNum = 1; arma::uword mWorkRangeSize = 0; - std::optional> mWorkRange; + std::pair mWorkRange = std::make_pair(arma::uword(0), arma::uword(0)); arma::mat mBetasSE; //!< \~english Standard errors of coefficient estimates. \~chinese 回归系数估计值的标准差。 arma::vec mSHat; //!< \~english A vector of \f$tr(S)\f$ and \f$tr(SS^T)\f$. \~chinese 由 \f$tr(S)\f$ 和 \f$tr(SS^T)\f$ 组成的向量。 diff --git a/include/gwmodelpp/GWRMultiscale.h b/include/gwmodelpp/GWRMultiscale.h index af2e77f8..ec00750c 100644 --- a/include/gwmodelpp/GWRMultiscale.h +++ b/include/gwmodelpp/GWRMultiscale.h @@ -817,7 +817,7 @@ class GWRMultiscale : public SpatialMultiscaleAlgorithm, #ifdef ENABLE_MPI arma::uword mWorkRangeSize = 0; #endif // ENABLE_MPI - std::optional> mWorkRange; + std::pair mWorkRange = std::make_pair(arma::uword(0), arma::uword(0)); public: static int treeChildCount; //!< \~english \~chinese diff --git a/src/gwmodelpp/GWRBasic.cpp b/src/gwmodelpp/GWRBasic.cpp index 12720c74..8985f853 100644 --- a/src/gwmodelpp/GWRBasic.cpp +++ b/src/gwmodelpp/GWRBasic.cpp @@ -63,6 +63,8 @@ mat GWRBasic::fit() { GWM_LOG_STAGE("Initializing"); uword nDp = mCoords.n_rows, nVars = mX.n_cols; + mWorkRange = make_pair(uword(0), nDp); + #ifdef ENABLE_MPI // Sync x, y, coords with other processes if (mParallelType & ParallelType::MPI) @@ -250,12 +252,11 @@ mat GWRBasic::fitCoreSerial(const mat &x, const vec &y, const SpatialWeight &sw) mBetasSE = mat(nVar, nDp, fill::zeros); mSHat = vec(2, fill::zeros); mQDiag = vec(nDp, fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; mS = mat(isStoreS() ? rangeSize : 1, nDp, fill::zeros); mC = cube(nVar, nDp, isStoreC() ? rangeSize : 1, fill::zeros); - // cout << mWorkerId << " process work range: [" << workRange.first << "," << workRange.second << "]\n"; - for (uword i = workRange.first; i < workRange.second; i++) + // cout << mWorkerId << " process work range: [" << mWorkRange.first << "," << mWorkRange.second << "]\n"; + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = mSpatialWeight.weightVector(i); @@ -274,8 +275,8 @@ mat GWRBasic::fitCoreSerial(const mat &x, const vec &y, const SpatialWeight &sw) vec p = - si.t(); p(i) += 1.0; mQDiag += p % p; - mS.row(isStoreS() ? (i - workRange.first) : 0) = si; - mC.slice(isStoreC() ? (i - workRange.first) : 0) = ci; + mS.row(isStoreS() ? (i - mWorkRange.first) : 0) = si; + mC.slice(isStoreC() ? (i - mWorkRange.first) : 0) = ci; } catch (const exception& e) { @@ -292,8 +293,7 @@ mat GWRBasic::fitCoreCVSerial(const mat& x, const vec& y, const SpatialWeight& s { uword nDp = mCoords.n_rows, nVar = x.n_cols; mat betas(nVar, nDp, fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = sw.weightVector(i); @@ -320,8 +320,7 @@ mat GWRBasic::fitCoreSHatSerial(const mat& x, const vec& y, const SpatialWeight& uword nDp = mCoords.n_rows, nVar = x.n_cols; mat betas(nVar, nDp, fill::zeros); shat = vec(2, arma::fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = sw.weightVector(i); @@ -460,8 +459,7 @@ mat GWRBasic::fitCoreOmp(const mat& x, const vec& y, const SpatialWeight& sw) uword nDp = mCoords.n_rows, nVar = x.n_cols; mat betas(nVar, nDp, fill::zeros); mBetasSE = mat(nVar, nDp, fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; mS = mat(isStoreS() ? rangeSize : 1, nDp, fill::zeros); mC = cube(nVar, nDp, isStoreC() ? rangeSize : 1, fill::zeros); mat shat_all(2, mOmpThreadNum, fill::zeros); @@ -469,7 +467,7 @@ mat GWRBasic::fitCoreOmp(const mat& x, const vec& y, const SpatialWeight& sw) bool success = true; std::exception except; #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); if (success) @@ -491,8 +489,8 @@ mat GWRBasic::fitCoreOmp(const mat& x, const vec& y, const SpatialWeight& sw) vec p = - si.t(); p(i) += 1.0; qDiag_all.col(thread) += p % p; - mS.row(isStoreS() ? (i - workRange.first) : 0) = si; - mC.slice(isStoreC() ? (i - workRange.first) : 0) = ci; + mS.row(isStoreS() ? (i - mWorkRange.first) : 0) = si; + mC.slice(isStoreC() ? (i - mWorkRange.first) : 0) = ci; } catch (const exception& e) { @@ -518,9 +516,8 @@ arma::mat GWRBasic::fitCoreCVOmp(const arma::mat& x, const arma::vec& y, const S uword nDp = mCoords.n_rows, nVar = x.n_cols; mat betas(nVar, nDp, fill::zeros); bool flag = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); if (flag) @@ -551,9 +548,8 @@ arma::mat GWRBasic::fitCoreSHatOmp(const arma::mat& x, const arma::vec& y, const mat betas(nVar, nDp, fill::zeros); mat shat_all(2, mOmpThreadNum, fill::zeros); int flag = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); if (flag) @@ -597,8 +593,7 @@ mat GWRBasic::fitCoreCuda(const mat& x, const vec& y, const SpatialWeight& sw) mQDiag = vec(nDp, arma::fill::zeros); mS = mat(isStoreS() ? nDp : 1, nDp, arma::fill::zeros); mC = cube(nVar, nDp, isStoreC() ? nDp : 1, arma::fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); cumat u_xt(xt), u_y(y); cumat u_betas(nVar, nDp), u_betasSE(nVar, nDp); @@ -616,7 +611,7 @@ mat GWRBasic::fitCoreCuda(const mat& x, const vec& y, const SpatialWeight& sw) for (size_t i = 0; i < groups; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -676,8 +671,7 @@ arma::mat GWRBasic::predictCuda(const mat& locations, const mat& x, const vec& y { uword nRp = locations.n_rows, nDp = mCoords.n_rows, nVar = x.n_cols; mat betas(nVar, nRp, fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); cumat u_xt(x.t()), u_y(y); custride u_xtwx(nVar, nVar, mGroupLength), u_xtwy(nVar, 1, mGroupLength), u_xtwxI(nVar, nVar, mGroupLength); @@ -690,7 +684,7 @@ arma::mat GWRBasic::predictCuda(const mat& locations, const mat& x, const vec& y for (size_t i = 0; i < groups; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(mSpatialWeight.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -731,13 +725,12 @@ arma::mat GWRBasic::fitCoreCVCuda(const arma::mat& x, const arma::vec& y, const p_info = new int[mGroupLength]; checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength)); bool success = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); for (size_t i = 0; i < groups && success; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -782,13 +775,12 @@ arma::mat GWRBasic::fitCoreSHatCuda(const arma::mat& x, const arma::vec& y, cons p_info = new int[mGroupLength]; checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength)); bool success = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); for (size_t i = 0; i < groups && success; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -991,11 +983,10 @@ arma::mat gwm::GWRBasic::fitMpi() // } GWM_MPI_MASTER_END GWM_MPI_WORKER_BEGIN - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - // printf("%d process work range: [%lld, %lld]\n", mWorkerId, workRange.first, workRange.second); - // cout << mWorkerId << " process work range: [" << workRange.first << "," << workRange.second << "]\n"; - mat betas = mBetas.cols(workRange.first, workRange.second - 1); - mat betasSE = mBetasSE.cols(workRange.first, workRange.second - 1); + // printf("%d process work range: [%lld, %lld]\n", mWorkerId, mWorkRange.first, mWorkRange.second); + // cout << mWorkerId << " process work range: [" << mWorkRange.first << "," << mWorkRange.second << "]\n"; + mat betas = mBetas.cols(mWorkRange.first, mWorkRange.second - 1); + mat betasSE = mBetasSE.cols(mWorkRange.first, mWorkRange.second - 1); MPI_Send(betas.memptr(), betas.n_elem, MPI_DOUBLE, 0, int(GWRBasicFitMpiTags::Betas), MPI_COMM_WORLD); MPI_Send(betasSE.memptr(), betasSE.n_elem, MPI_DOUBLE, 0, int(GWRBasicFitMpiTags::BetasSE), MPI_COMM_WORLD); MPI_Send(mSHat.memptr(), mSHat.n_elem, MPI_DOUBLE, 0, int(GWRBasicFitMpiTags::SHat), MPI_COMM_WORLD); diff --git a/src/gwmodelpp/GWRMultiscale.cpp b/src/gwmodelpp/GWRMultiscale.cpp index 3e03f4a7..4ffa5b25 100644 --- a/src/gwmodelpp/GWRMultiscale.cpp +++ b/src/gwmodelpp/GWRMultiscale.cpp @@ -80,6 +80,8 @@ RegressionDiagnostic GWRMultiscale::CalcDiagnostic(const mat &x, const vec &y, c mat GWRMultiscale::fit() { uword nDp = mX.n_rows, nVar = mX.n_cols; + mWorkRange = make_pair(uword(0), nDp); + // ******************************** // Centering and scaling predictors // ******************************** @@ -126,6 +128,7 @@ mat GWRMultiscale::fit() case ParallelType::MPI_CUDA: gwr.setGroupSize(mGroupLength); gwr.setGPUId(mGpuId); + break; default: break; } @@ -146,6 +149,7 @@ mat GWRMultiscale::fit() createDistanceParameter(nVar); mMaxDistances.resize(nVar); mMinDistances.resize(nVar); + #ifdef ENABLE_MPI if (mParallelType & ParallelType::MPI) { @@ -226,8 +230,6 @@ mat GWRMultiscale::fit() GWM_LOG_STOP_RETURN(mStatus, mat(nDp, nVar, arma::fill::zeros)); } - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - GWM_LOG_STAGE("Initializing diagnostic matrices"); mat idm = eye(nVar, nVar); if (mHasHatMatrix) @@ -237,9 +239,9 @@ mat GWRMultiscale::fit() mC = gwr.c(); for (uword i = 0; i < nVar; ++i) { - for (uword j = workRange.first; j < workRange.second ; ++j) + for (uword j = mWorkRange.first; j < mWorkRange.second ; ++j) { - uword e = j - workRange.first; + uword e = j - mWorkRange.first; mSArray.slice(i).row(e) = mX(j, i) * (idm.row(i) * mC.slice(e)); } } @@ -258,7 +260,7 @@ mat GWRMultiscale::fit() if (mParallelType & ParallelType::MPI) { vec shati(2); - shati(0) = trace(mS0.cols(workRange.first, workRange.second - 1)); + shati(0) = trace(mS0.cols(mWorkRange.first, mWorkRange.second - 1)); shati(1) = trace(mS0 * mS0.t()); MPI_Reduce(shati.memptr(), shat.memptr(), 2, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Bcast(shat.memptr(), 2, MPI_DOUBLE, 0, MPI_COMM_WORLD); @@ -492,13 +494,12 @@ vec GWRMultiscale::fitVarCoreSerial(const vec &x, const vec &y, const SpatialWei mat betas(1, nDp, fill::zeros); bool success = true; std::exception except; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); if (mHasHatMatrix) { mat ci, si; - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; S = mat(rangeSize, nDp, fill::zeros); - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = sw.weightVector(i); @@ -511,7 +512,7 @@ vec GWRMultiscale::fitVarCoreSerial(const vec &x, const vec &y, const SpatialWei betas.col(i) = xtwx_inv * xtwy; ci = xtwx_inv * xtw; si = x(i) * ci; - S.row(i - workRange.first) = si; + S.row(i - mWorkRange.first) = si; } catch (const exception& e) { @@ -524,7 +525,7 @@ vec GWRMultiscale::fitVarCoreSerial(const vec &x, const vec &y, const SpatialWei } else { - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = sw.weightVector(i); @@ -556,8 +557,7 @@ vec GWRMultiscale::fitVarCoreCVSerial(const vec &x, const vec &y, const SpatialW { uword nDp = x.n_rows; vec beta(nDp, fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = sw.weightVector(i); @@ -584,8 +584,7 @@ vec GWRMultiscale::fitVarCoreSHatSerial(const vec &x, const vec &y, const Spatia uword nDp = x.n_rows; vec betas(nDp, fill::zeros); shat = vec(2, fill::zeros); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = sw.weightVector(i); @@ -618,13 +617,12 @@ vec GWRMultiscale::fitVarCoreOmp(const vec &x, const vec &y, const SpatialWeight mat betas(1, nDp, fill::zeros); bool success = true; std::exception except; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); if (mHasHatMatrix) { - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; S = mat(rangeSize, nDp, fill::zeros); #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); vec w = sw.weightVector(i); @@ -637,7 +635,7 @@ vec GWRMultiscale::fitVarCoreOmp(const vec &x, const vec &y, const SpatialWeight betas.col(i) = xtwx_inv * xtwy; mat ci = xtwx_inv * xtw; mat si = x(i) * ci; - S.row(i - workRange.first) = si; + S.row(i - mWorkRange.first) = si; } catch (const exception& e) { @@ -651,7 +649,7 @@ vec GWRMultiscale::fitVarCoreOmp(const vec &x, const vec &y, const SpatialWeight else { #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); vec w = sw.weightVector(i); @@ -684,9 +682,8 @@ vec GWRMultiscale::fitVarCoreCVOmp(const vec &x, const vec &y, const SpatialWeig uword nDp = mCoords.n_rows; vec beta(nDp, fill::zeros); bool flag = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); if (flag) @@ -717,9 +714,8 @@ vec GWRMultiscale::fitVarCoreSHatOmp(const vec &x, const vec &y, const SpatialWe vec betas(nDp, fill::zeros); mat shat_all(2, mOmpThreadNum, fill::zeros); bool flag = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); #pragma omp parallel for num_threads(mOmpThreadNum) - for (uword i = workRange.first; i < workRange.second; i++) + for (uword i = mWorkRange.first; i < mWorkRange.second; i++) { GWM_LOG_STOP_CONTINUE(mStatus); if (flag) @@ -764,8 +760,7 @@ vec GWRMultiscale::fitVarCoreCuda(const vec &x, const vec &y, const SpatialWeigh int *d_info, *p_info; p_info = new int[mGroupLength]; checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength)); - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); S = mat(mHasHatMatrix ? rangeSize : 1, nDp, fill::zeros); if (mHasHatMatrix) @@ -774,7 +769,7 @@ vec GWRMultiscale::fitVarCoreCuda(const vec &x, const vec &y, const SpatialWeigh for (size_t i = 0; i < groups; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -806,7 +801,7 @@ vec GWRMultiscale::fitVarCoreCuda(const vec &x, const vec &y, const SpatialWeigh for (size_t i = 0; i < groups; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -846,13 +841,12 @@ vec GWRMultiscale::fitVarCoreCVCuda(const vec &x, const vec &y, const SpatialWei p_info = new int[mGroupLength]; checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength)); bool success = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); for (size_t i = 0; i < groups && success; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); @@ -899,13 +893,12 @@ vec GWRMultiscale::fitVarCoreSHatCuda(const vec &x, const vec &y, const SpatialW p_info = new int[mGroupLength]; checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength)); bool success = true; - std::pair workRange = mWorkRange.value_or(make_pair(0, nDp)); - uword rangeSize = workRange.second - workRange.first; + uword rangeSize = mWorkRange.second - mWorkRange.first; size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1); for (size_t i = 0; i < groups && success; i++) { GWM_LOG_STOP_BREAK(mStatus); - size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength; + size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength; for (size_t j = 0, e = begin + j; j < length; j++, e++) { checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem())); From 68391622d5fe5ea63b2b5962c293b8a0142e602e Mon Sep 17 00:00:00 2001 From: HPDell Date: Sun, 23 Feb 2025 21:19:43 +0000 Subject: [PATCH 3/5] edit: move gwr fit process to fitInitial --- include/gwmodelpp/GWRMultiscale.h | 2 + src/gwmodelpp/GWRMultiscale.cpp | 94 +++++++++++++++++-------------- 2 files changed, 54 insertions(+), 42 deletions(-) diff --git a/include/gwmodelpp/GWRMultiscale.h b/include/gwmodelpp/GWRMultiscale.h index ec00750c..b3d0f5d1 100644 --- a/include/gwmodelpp/GWRMultiscale.h +++ b/include/gwmodelpp/GWRMultiscale.h @@ -645,6 +645,8 @@ class GWRMultiscale : public SpatialMultiscaleAlgorithm, return mSpatialWeights[i].weight(); } + arma::mat fitInitial(); + /** * \~english * @brief The backfitting algorithm. diff --git a/src/gwmodelpp/GWRMultiscale.cpp b/src/gwmodelpp/GWRMultiscale.cpp index 4ffa5b25..a281cdb4 100644 --- a/src/gwmodelpp/GWRMultiscale.cpp +++ b/src/gwmodelpp/GWRMultiscale.cpp @@ -101,48 +101,7 @@ mat GWRMultiscale::fit() // Calculate the initial beta0 from the above bandwidths // ***************************************************** GWM_LOG_STAGE("Calculating initial betas"); - GWRBasic gwr; - gwr.setCoords(mCoords); - gwr.setDependentVariable(mY); - gwr.setIndependentVariables(mX); - gwr.setSpatialWeight(mInitSpatialWeight); - gwr.setIsAutoselectBandwidth(true); - switch (mBandwidthSelectionApproach[0]) - { - case GWRMultiscale::BandwidthSelectionCriterionType::CV: - gwr.setBandwidthSelectionCriterion(GWRBasic::BandwidthSelectionCriterionType::CV); - break; - case GWRMultiscale::BandwidthSelectionCriterionType::AIC: - gwr.setBandwidthSelectionCriterion(GWRBasic::BandwidthSelectionCriterionType::AIC); - default: - break; - } - gwr.setParallelType(mParallelType); - switch (mParallelType) - { - case ParallelType::OpenMP: - case ParallelType::MPI_MP: - gwr.setOmpThreadNum(mOmpThreadNum); - break; - case ParallelType::CUDA: - case ParallelType::MPI_CUDA: - gwr.setGroupSize(mGroupLength); - gwr.setGPUId(mGpuId); - break; - default: - break; - } - if (mParallelType & ParallelType::MPI) - { - gwr.setWorkerId(mWorkerId); - gwr.setWorkerNum(mWorkerNum); - } - gwr.setStoreS(mHasHatMatrix); - gwr.setStoreC(mHasHatMatrix); - if (mGoldenLowerBounds.has_value()) gwr.setGoldenLowerBounds(mGoldenLowerBounds.value()); - if (mGoldenUpperBounds.has_value()) gwr.setGoldenUpperBounds(mGoldenUpperBounds.value()); - mat betas = gwr.fit(); - mBetasSE = gwr.betasSE(); + mat betas = fitInitial(); GWM_LOG_STOP_RETURN(mStatus, mat(nDp, nVar, arma::fill::zeros)); GWM_LOG_STAGE("Initializing"); @@ -292,6 +251,57 @@ mat GWRMultiscale::fit() return mBetas; } +arma::mat gwm::GWRMultiscale::fitInitial() +{ + GWRBasic gwr; + gwr.setCoords(mCoords); + gwr.setDependentVariable(mY); + gwr.setIndependentVariables(mX); + gwr.setSpatialWeight(mInitSpatialWeight); + gwr.setIsAutoselectBandwidth(true); + switch (mBandwidthSelectionApproach[0]) + { + case GWRMultiscale::BandwidthSelectionCriterionType::CV: + gwr.setBandwidthSelectionCriterion(GWRBasic::BandwidthSelectionCriterionType::CV); + break; + case GWRMultiscale::BandwidthSelectionCriterionType::AIC: + gwr.setBandwidthSelectionCriterion(GWRBasic::BandwidthSelectionCriterionType::AIC); + default: + break; + } + gwr.setParallelType(mParallelType); + switch (mParallelType) + { + case ParallelType::OpenMP: + case ParallelType::MPI_MP: + gwr.setOmpThreadNum(mOmpThreadNum); + break; + case ParallelType::CUDA: + case ParallelType::MPI_CUDA: + gwr.setGroupSize(mGroupLength); + gwr.setGPUId(mGpuId); + break; + default: + break; + } + if (mParallelType & ParallelType::MPI) + { + gwr.setWorkerId(mWorkerId); + gwr.setWorkerNum(mWorkerNum); + } + gwr.setStoreS(mHasHatMatrix); + gwr.setStoreC(mHasHatMatrix); + if (mGoldenLowerBounds.has_value()) gwr.setGoldenLowerBounds(mGoldenLowerBounds.value()); + if (mGoldenUpperBounds.has_value()) gwr.setGoldenUpperBounds(mGoldenUpperBounds.value()); + mat betas = gwr.fit(); + mBetasSE = gwr.betasSE(); + if (mHasHatMatrix) { + mS0 = gwr.s(); + mC = gwr.c(); + } + return betas; +} + mat GWRMultiscale::backfitting(const mat& betas0) { uword nDp = mCoords.n_rows, nVar = mX.n_cols; From a806a801a2de93e566e88daa2aa6e35f745ec026 Mon Sep 17 00:00:00 2001 From: HPDell Date: Sun, 23 Feb 2025 21:19:55 +0000 Subject: [PATCH 4/5] edit: handle min and max distance --- src/gwmodelpp/GWRMultiscale.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/gwmodelpp/GWRMultiscale.cpp b/src/gwmodelpp/GWRMultiscale.cpp index a281cdb4..116fc03d 100644 --- a/src/gwmodelpp/GWRMultiscale.cpp +++ b/src/gwmodelpp/GWRMultiscale.cpp @@ -150,6 +150,11 @@ mat GWRMultiscale::fit() } } #endif // ENABLE_CUDA + for (size_t i = 0; i < nVar; i++) + { + mMinDistances[i] = mSpatialWeights[i].distance()->minDistance(); + mMaxDistances[i] = mSpatialWeights[i].distance()->maxDistance(); + } GWM_LOG_STOP_RETURN(mStatus, mat(nDp, nVar, arma::fill::zeros)); // *********************** @@ -167,8 +172,6 @@ mat GWRMultiscale::fit() mXi = mX.col(i); BandwidthWeight* bw0 = bandwidth(i); bool adaptive = bw0->adaptive(); - mMaxDistances[i] = mSpatialWeights[i].distance()->maxDistance(); - mMinDistances[i] = mSpatialWeights[i].distance()->minDistance(); GWM_LOG_INFO(string(GWM_LOG_TAG_MGWR_INITIAL_BW) + to_string(i)); BandwidthSelector selector; @@ -193,9 +196,7 @@ mat GWRMultiscale::fit() mat idm = eye(nVar, nVar); if (mHasHatMatrix) { - mS0 = gwr.s(); mSArray = cube(mS0.n_rows, mS0.n_cols, nVar, fill::zeros); - mC = gwr.c(); for (uword i = 0; i < nVar; ++i) { for (uword j = mWorkRange.first; j < mWorkRange.second ; ++j) @@ -336,9 +337,8 @@ mat GWRMultiscale::backfitting(const mat& betas0) bool adaptive = bwi0->adaptive(); BandwidthSelector selector; selector.setBandwidth(bwi0); - double maxDist = mMaxDistances[i], minDist = mMinDistances[i]; - selector.setLower(mGoldenLowerBounds.value_or(adaptive ? mAdaptiveLower : minDist)); - selector.setUpper(mGoldenUpperBounds.value_or(adaptive ? mCoords.n_rows : maxDist)); + selector.setLower(mGoldenLowerBounds.value_or(adaptive ? mAdaptiveLower : mMinDistances[i])); + selector.setUpper(mGoldenUpperBounds.value_or(adaptive ? mCoords.n_rows : mMaxDistances[i])); BandwidthWeight* bwi = selector.optimize(this); double bwi0s = bwi0->bandwidth(), bwi1s = bwi->bandwidth(); vector vbs_args { From d74545308193dd1cefecc5094bfd1fa8210eaf0d Mon Sep 17 00:00:00 2001 From: HPDell Date: Mon, 24 Feb 2025 13:05:02 +0000 Subject: [PATCH 5/5] fix: init mWorkRange in GWRRobust --- src/gwmodelpp/GWRRobust.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/gwmodelpp/GWRRobust.cpp b/src/gwmodelpp/GWRRobust.cpp index 2dfb2926..dc8c8e1f 100644 --- a/src/gwmodelpp/GWRRobust.cpp +++ b/src/gwmodelpp/GWRRobust.cpp @@ -31,6 +31,7 @@ mat GWRRobust::fit() { GWM_LOG_STAGE("Initializing"); uword nDp = mCoords.n_rows, nVar = mX.n_cols; + mWorkRange = make_pair(uword(0), nDp); createDistanceParameter(); GWM_LOG_STOP_RETURN(mStatus, mat(nDp, nVar, arma::fill::zeros));