Skip to content

Commit b15dcf8

Browse files
HPDellHPDell
authored
Fix: issues caused by workRange (#101)
* fix: min distance remove 0 * edit: make mWorkRange as pair * edit: move gwr fit process to fitInitial * edit: handle min and max distance * fix: init mWorkRange in GWRRobust
1 parent b698f2e commit b15dcf8

File tree

6 files changed

+118
-120
lines changed

6 files changed

+118
-120
lines changed

‎include/gwmodelpp/GWRBasic.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -763,7 +763,7 @@ class GWRBasic : public GWRBase, public IBandwidthSelectable, public IVarialbeSe
763763
int mWorkerId = 0;
764764
int mWorkerNum = 1;
765765
arma::uword mWorkRangeSize = 0;
766-
std::optional<std::pair<arma::uword, arma::uword>> mWorkRange;
766+
std::pair<arma::uword, arma::uword> mWorkRange = std::make_pair(arma::uword(0), arma::uword(0));
767767

768768
arma::mat mBetasSE; //!< \~english Standard errors of coefficient estimates. \~chinese 回归系数估计值的标准差。
769769
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$ 组成的向量。

‎include/gwmodelpp/GWRMultiscale.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -645,6 +645,8 @@ class GWRMultiscale : public SpatialMultiscaleAlgorithm,
645645
return mSpatialWeights[i].weight<BandwidthWeight>();
646646
}
647647

648+
arma::mat fitInitial();
649+
648650
/**
649651
* \~english
650652
* @brief The backfitting algorithm.
@@ -817,7 +819,7 @@ class GWRMultiscale : public SpatialMultiscaleAlgorithm,
817819
#ifdef ENABLE_MPI
818820
arma::uword mWorkRangeSize = 0;
819821
#endif // ENABLE_MPI
820-
std::optional<std::pair<arma::uword, arma::uword>> mWorkRange;
822+
std::pair<arma::uword, arma::uword> mWorkRange = std::make_pair(arma::uword(0), arma::uword(0));
821823

822824
public:
823825
static int treeChildCount; //!< \~english \~chinese

‎src/gwmodelpp/GWRBasic.cpp‎

Lines changed: 27 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,8 @@ mat GWRBasic::fit()
6363
{
6464
GWM_LOG_STAGE("Initializing");
6565
uword nDp = mCoords.n_rows, nVars = mX.n_cols;
66+
mWorkRange = make_pair(uword(0), nDp);
67+
6668
#ifdef ENABLE_MPI
6769
// Sync x, y, coords with other processes
6870
if (mParallelType & ParallelType::MPI)
@@ -250,12 +252,11 @@ mat GWRBasic::fitCoreSerial(const mat &x, const vec &y, const SpatialWeight &sw)
250252
mBetasSE = mat(nVar, nDp, fill::zeros);
251253
mSHat = vec(2, fill::zeros);
252254
mQDiag = vec(nDp, fill::zeros);
253-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
254-
uword rangeSize = workRange.second - workRange.first;
255+
uword rangeSize = mWorkRange.second - mWorkRange.first;
255256
mS = mat(isStoreS() ? rangeSize : 1, nDp, fill::zeros);
256257
mC = cube(nVar, nDp, isStoreC() ? rangeSize : 1, fill::zeros);
257-
// cout << mWorkerId << " process work range: [" << workRange.first << "," << workRange.second << "]\n";
258-
for (uword i = workRange.first; i < workRange.second; i++)
258+
// cout << mWorkerId << " process work range: [" << mWorkRange.first << "," << mWorkRange.second << "]\n";
259+
for (uword i = mWorkRange.first; i < mWorkRange.second; i++)
259260
{
260261
GWM_LOG_STOP_BREAK(mStatus);
261262
vec w = mSpatialWeight.weightVector(i);
@@ -274,8 +275,8 @@ mat GWRBasic::fitCoreSerial(const mat &x, const vec &y, const SpatialWeight &sw)
274275
vec p = - si.t();
275276
p(i) += 1.0;
276277
mQDiag += p % p;
277-
mS.row(isStoreS() ? (i - workRange.first) : 0) = si;
278-
mC.slice(isStoreC() ? (i - workRange.first) : 0) = ci;
278+
mS.row(isStoreS() ? (i - mWorkRange.first) : 0) = si;
279+
mC.slice(isStoreC() ? (i - mWorkRange.first) : 0) = ci;
279280
}
280281
catch (const exception& e)
281282
{
@@ -292,8 +293,7 @@ mat GWRBasic::fitCoreCVSerial(const mat& x, const vec& y, const SpatialWeight& s
292293
{
293294
uword nDp = mCoords.n_rows, nVar = x.n_cols;
294295
mat betas(nVar, nDp, fill::zeros);
295-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
296-
for (uword i = workRange.first; i < workRange.second; i++)
296+
for (uword i = mWorkRange.first; i < mWorkRange.second; i++)
297297
{
298298
GWM_LOG_STOP_BREAK(mStatus);
299299
vec w = sw.weightVector(i);
@@ -320,8 +320,7 @@ mat GWRBasic::fitCoreSHatSerial(const mat& x, const vec& y, const SpatialWeight&
320320
uword nDp = mCoords.n_rows, nVar = x.n_cols;
321321
mat betas(nVar, nDp, fill::zeros);
322322
shat = vec(2, arma::fill::zeros);
323-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
324-
for (uword i = workRange.first; i < workRange.second; i++)
323+
for (uword i = mWorkRange.first; i < mWorkRange.second; i++)
325324
{
326325
GWM_LOG_STOP_BREAK(mStatus);
327326
vec w = sw.weightVector(i);
@@ -460,16 +459,15 @@ mat GWRBasic::fitCoreOmp(const mat& x, const vec& y, const SpatialWeight& sw)
460459
uword nDp = mCoords.n_rows, nVar = x.n_cols;
461460
mat betas(nVar, nDp, fill::zeros);
462461
mBetasSE = mat(nVar, nDp, fill::zeros);
463-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
464-
uword rangeSize = workRange.second - workRange.first;
462+
uword rangeSize = mWorkRange.second - mWorkRange.first;
465463
mS = mat(isStoreS() ? rangeSize : 1, nDp, fill::zeros);
466464
mC = cube(nVar, nDp, isStoreC() ? rangeSize : 1, fill::zeros);
467465
mat shat_all(2, mOmpThreadNum, fill::zeros);
468466
mat qDiag_all(nDp, mOmpThreadNum, fill::zeros);
469467
bool success = true;
470468
std::exception except;
471469
#pragma omp parallel for num_threads(mOmpThreadNum)
472-
for (uword i = workRange.first; i < workRange.second; i++)
470+
for (uword i = mWorkRange.first; i < mWorkRange.second; i++)
473471
{
474472
GWM_LOG_STOP_CONTINUE(mStatus);
475473
if (success)
@@ -491,8 +489,8 @@ mat GWRBasic::fitCoreOmp(const mat& x, const vec& y, const SpatialWeight& sw)
491489
vec p = - si.t();
492490
p(i) += 1.0;
493491
qDiag_all.col(thread) += p % p;
494-
mS.row(isStoreS() ? (i - workRange.first) : 0) = si;
495-
mC.slice(isStoreC() ? (i - workRange.first) : 0) = ci;
492+
mS.row(isStoreS() ? (i - mWorkRange.first) : 0) = si;
493+
mC.slice(isStoreC() ? (i - mWorkRange.first) : 0) = ci;
496494
}
497495
catch (const exception& e)
498496
{
@@ -518,9 +516,8 @@ arma::mat GWRBasic::fitCoreCVOmp(const arma::mat& x, const arma::vec& y, const S
518516
uword nDp = mCoords.n_rows, nVar = x.n_cols;
519517
mat betas(nVar, nDp, fill::zeros);
520518
bool flag = true;
521-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
522519
#pragma omp parallel for num_threads(mOmpThreadNum)
523-
for (uword i = workRange.first; i < workRange.second; i++)
520+
for (uword i = mWorkRange.first; i < mWorkRange.second; i++)
524521
{
525522
GWM_LOG_STOP_CONTINUE(mStatus);
526523
if (flag)
@@ -551,9 +548,8 @@ arma::mat GWRBasic::fitCoreSHatOmp(const arma::mat& x, const arma::vec& y, const
551548
mat betas(nVar, nDp, fill::zeros);
552549
mat shat_all(2, mOmpThreadNum, fill::zeros);
553550
int flag = true;
554-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
555551
#pragma omp parallel for num_threads(mOmpThreadNum)
556-
for (uword i = workRange.first; i < workRange.second; i++)
552+
for (uword i = mWorkRange.first; i < mWorkRange.second; i++)
557553
{
558554
GWM_LOG_STOP_CONTINUE(mStatus);
559555
if (flag)
@@ -597,8 +593,7 @@ mat GWRBasic::fitCoreCuda(const mat& x, const vec& y, const SpatialWeight& sw)
597593
mQDiag = vec(nDp, arma::fill::zeros);
598594
mS = mat(isStoreS() ? nDp : 1, nDp, arma::fill::zeros);
599595
mC = cube(nVar, nDp, isStoreC() ? nDp : 1, arma::fill::zeros);
600-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
601-
uword rangeSize = workRange.second - workRange.first;
596+
uword rangeSize = mWorkRange.second - mWorkRange.first;
602597
size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1);
603598
cumat u_xt(xt), u_y(y);
604599
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)
616611
for (size_t i = 0; i < groups; i++)
617612
{
618613
GWM_LOG_STOP_BREAK(mStatus);
619-
size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength;
614+
size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength;
620615
for (size_t j = 0, e = begin + j; j < length; j++, e++)
621616
{
622617
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
676671
{
677672
uword nRp = locations.n_rows, nDp = mCoords.n_rows, nVar = x.n_cols;
678673
mat betas(nVar, nRp, fill::zeros);
679-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
680-
uword rangeSize = workRange.second - workRange.first;
674+
uword rangeSize = mWorkRange.second - mWorkRange.first;
681675
size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1);
682676
cumat u_xt(x.t()), u_y(y);
683677
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
690684
for (size_t i = 0; i < groups; i++)
691685
{
692686
GWM_LOG_STOP_BREAK(mStatus);
693-
size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength;
687+
size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength;
694688
for (size_t j = 0, e = begin + j; j < length; j++, e++)
695689
{
696690
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
731725
p_info = new int[mGroupLength];
732726
checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength));
733727
bool success = true;
734-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
735-
uword rangeSize = workRange.second - workRange.first;
728+
uword rangeSize = mWorkRange.second - mWorkRange.first;
736729
size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1);
737730
for (size_t i = 0; i < groups && success; i++)
738731
{
739732
GWM_LOG_STOP_BREAK(mStatus);
740-
size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength;
733+
size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength;
741734
for (size_t j = 0, e = begin + j; j < length; j++, e++)
742735
{
743736
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
782775
p_info = new int[mGroupLength];
783776
checkCudaErrors(cudaMalloc(&d_info, sizeof(int) * mGroupLength));
784777
bool success = true;
785-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
786-
uword rangeSize = workRange.second - workRange.first;
778+
uword rangeSize = mWorkRange.second - mWorkRange.first;
787779
size_t groups = rangeSize / mGroupLength + (rangeSize % mGroupLength == 0 ? 0 : 1);
788780
for (size_t i = 0; i < groups && success; i++)
789781
{
790782
GWM_LOG_STOP_BREAK(mStatus);
791-
size_t begin = workRange.first + i * mGroupLength, length = (begin + mGroupLength > workRange.second) ? (workRange.second - begin) : mGroupLength;
783+
size_t begin = mWorkRange.first + i * mGroupLength, length = (begin + mGroupLength > mWorkRange.second) ? (mWorkRange.second - begin) : mGroupLength;
792784
for (size_t j = 0, e = begin + j; j < length; j++, e++)
793785
{
794786
checkCudaErrors(sw.weightVector(e, u_dists.dmem(), u_weights.dmem()));
@@ -991,11 +983,10 @@ arma::mat gwm::GWRBasic::fitMpi()
991983
// }
992984
GWM_MPI_MASTER_END
993985
GWM_MPI_WORKER_BEGIN
994-
std::pair<uword, uword> workRange = mWorkRange.value_or(make_pair(0, nDp));
995-
// printf("%d process work range: [%lld, %lld]\n", mWorkerId, workRange.first, workRange.second);
996-
// cout << mWorkerId << " process work range: [" << workRange.first << "," << workRange.second << "]\n";
997-
mat betas = mBetas.cols(workRange.first, workRange.second - 1);
998-
mat betasSE = mBetasSE.cols(workRange.first, workRange.second - 1);
986+
// printf("%d process work range: [%lld, %lld]\n", mWorkerId, mWorkRange.first, mWorkRange.second);
987+
// cout << mWorkerId << " process work range: [" << mWorkRange.first << "," << mWorkRange.second << "]\n";
988+
mat betas = mBetas.cols(mWorkRange.first, mWorkRange.second - 1);
989+
mat betasSE = mBetasSE.cols(mWorkRange.first, mWorkRange.second - 1);
999990
MPI_Send(betas.memptr(), betas.n_elem, MPI_DOUBLE, 0, int(GWRBasicFitMpiTags::Betas), MPI_COMM_WORLD);
1000991
MPI_Send(betasSE.memptr(), betasSE.n_elem, MPI_DOUBLE, 0, int(GWRBasicFitMpiTags::BetasSE), MPI_COMM_WORLD);
1001992
MPI_Send(mSHat.memptr(), mSHat.n_elem, MPI_DOUBLE, 0, int(GWRBasicFitMpiTags::SHat), MPI_COMM_WORLD);

0 commit comments

Comments
 (0)