Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix: Local-collinearity GWR #96

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Fix: Local-collinearity GWR #96

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7762b78
fix r2_adj in basic flow
RoryuOu Mar 19, 2024
783384f
update test, all tests passed
RoryuOu Mar 19, 2024
b666436
rebase master, add 1 test section
RoryuOu Dec 3, 2024
8f0962e
add local_cn and local_lambda with test
RoryuOu Jan 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 30 additions & 17 deletions include/gwmodelpp/GWRLocalCollinearity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 用于预测的函数

/**
Expand Down Expand Up @@ -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 获取是否进行带宽优选。
*
Expand Down Expand Up @@ -231,15 +241,15 @@ 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 拟合函数的多线程实现。
*
* @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 预测函数的非并行实现。
Expand Down Expand Up @@ -276,18 +286,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 岭回归。
Expand Down Expand Up @@ -322,12 +333,14 @@ 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;
double mTrStS = 0;
arma::vec mSHat;
// arma::vec mSHat;
arma::vec mLocalCN;
arma::vec mLocalLambda;

FitCalculator mFitFunction = &GWRLocalCollinearity::fitSerial;
PredictCalculator mPredictFunction = &GWRLocalCollinearity::predictSerial;
Expand Down
71 changes: 41 additions & 30 deletions src/gwmodelpp/GWRLocalCollinearity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,15 @@ 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);
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()
Expand Down Expand Up @@ -70,20 +71,16 @@ 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");
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);

return mBetas;
}
Expand Down Expand Up @@ -327,13 +324,14 @@ 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<nDp ;i++)
{
GWM_LOG_STOP_BREAK(mStatus);
Expand Down Expand Up @@ -363,13 +361,20 @@ mat GWRLocalCollinearity::fitSerial(const mat& x, const vec& y)
}
betas.row(i) = trans(ridgelm(wi,locallambda(i)) );
//如果没有给regressionpoint
mat xm = x;
// 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);
this->mTrS += 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;
Expand Down Expand Up @@ -416,12 +421,12 @@ mat GWRLocalCollinearity::predictSerial(const arma::mat &locations, const mat& x
}

#ifdef ENABLE_OPENMP
mat GWRLocalCollinearity::fitOmp(const mat& x, const vec& y)
mat GWRLocalCollinearity::fitOmp(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);
mat shat_all(2, mOmpThreadNum, fill::zeros);
//int current = 0;
Expand Down Expand Up @@ -456,15 +461,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);
Expand Down
Loading