implement gwm new feature

include/gwmodelpp/GWAverage.h

@@ -4,6 +4,7 @@
 #include "SpatialMonoscaleAlgorithm.h"
 #include "IMultivariableAnalysis.h"
 #include "IParallelizable.h"
+#include "IBandwidthSelectable.h"
 
 namespace gwm
 {
@@ -41,7 +42,7 @@ namespace gwm
  * - local interquartile ranges <- GWAverage::iqr()
  * - local quantile imbalances and coordinates <- GWAverage::qi()
  */
-class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis, public IParallelizable, public IParallelOpenmpEnabled
+class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysis, public IParallelizable, public IParallelOpenmpEnabled, public IBandwidthSelectable
 {
 public:
 
@@ -171,6 +172,11 @@ class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysi
     void setVariables(const arma::mat& x) override { mX = x; }
 
     void run() override;
+    // Implement IBandwidthSelectable
+    Status getCriterion(const std::unique_ptr<BandwidthWeight>& weight, double& criterion) override;
+    bool isAutoselectBandwidth() const { return mIsAutoselectBandwidth; }
+    void setIsAutoselectBandwidth(bool isAutoSelect) { mIsAutoselectBandwidth = isAutoSelect; }
+    const BandwidthCriterionList& bandwidthSelectionCriterionList() const { return mBandwidthSelectionCriterionList; }
 
     void calibration(const arma::mat& locations, const arma::mat& x);
 
@@ -239,6 +245,9 @@ class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysi
     arma::mat mIQR;           //!< \~english Local interquartile ranges \~chinese 局部分位距
     arma::mat mQI;            //!< \~english Local quantile imbalances and coordinates \~chinese 局部分位数不平衡度
 
+    bool mIsAutoselectBandwidth = false; //!< \~english Whether bandwidth autoselection is enabled \~chinese 是否启用带宽自动选择
+    BandwidthCriterionList mBandwidthSelectionCriterionList; //!< \~english List of bandwidth selection criterion values \~chinese 带宽优选指标值列表
+
     SummaryCalculator mSummaryFunction = &GWAverage::GWAverageSerial;  //!< \~english Calculator for summary statistics \~chinese 计算函数
     ParallelType mParallelType = ParallelType::SerialOnly;  //!< \~english Parallel type \~chinese 并行方法
     int mOmpThreadNum = 8;                                  //!< \~english Numbers of threads to be created while paralleling \~chinese 多线程所使用的线程数

src/gwmodelpp/GWAverage.cpp

@@ -1,4 +1,5 @@
 #include "GWAverage.h"
+#include "BandwidthSelector.h"
 #include <assert.h>
 #include <map>
 
@@ -67,6 +68,22 @@ void GWAverage::run()
     createDistanceParameter();
     GWM_LOG_STOP_RETURN(mStatus, void());
 
+    if (mIsAutoselectBandwidth)
+    {
+        GWM_LOG_STAGE("Bandwidth selection");
+        BandwidthWeight* bw0 = &mSpatialWeight.weight<BandwidthWeight>();
+        double lower = bw0->adaptive() ? 2.0 : mSpatialWeight.distance()->minDistance();
+        double upper = bw0->adaptive() ? (double)nRp : mSpatialWeight.distance()->maxDistance();
+
+        GWM_LOG_INFO(IBandwidthSelectable::infoBandwidthCriterion(*bw0));
+        BandwidthSelector selector(*bw0, lower, upper);
+        if (selector.optimize(this) == Status::Success)
+        {
+            mSpatialWeight.setWeight(selector.result());
+            mBandwidthSelectionCriterionList = selector.bandwidthCriterion();
+        }
+    }
+
     mLocalMean = mat(nRp, nVar, fill::zeros);
     mStandardDev = mat(nRp, nVar, fill::zeros);
     mLocalSkewness = mat(nRp, nVar, fill::zeros);
@@ -217,6 +234,77 @@ void GWAverage::setParallelType(const ParallelType &type)
     }
 }
 
+Status GWAverage::getCriterion(const std::unique_ptr<BandwidthWeight>& weight, double& criterion)
+{
+    if (!weight)
+    {
+        criterion = DBL_MAX;
+        return Status::Success;
+    }
+
+    uword n = mCoords.n_rows;
+    uword p = mX.n_cols;
+    if (n < 2 || p == 0)
+    {
+        criterion = DBL_MAX;
+        return Status::Success;
+    }
+
+    BandwidthWeight backup = mSpatialWeight.weight<BandwidthWeight>();
+    mSpatialWeight.setWeight(*weight);
+
+    double totalF = 0.0;
+    bool valid = false;
+
+    for (uword j = 0; j < p; j++)
+    {
+        mat weightedVals(n, n, fill::zeros);
+        for (uword i = 0; i < n; i++)
+        {
+            vec w = mSpatialWeight.weightVector(i);
+            if (w.n_rows != n)
+                continue;
+            weightedVals.col(i) = w % mX.col(j);
+        }
+
+        double yMean = accu(weightedVals) / double(n * n);
+        mat centered = weightedVals - yMean;
+        double SST = accu(centered % centered);
+        if (!isfinite(SST) || SST < 0.0)
+            continue;
+
+        rowvec groupMeans = mean(weightedVals, 0);
+        double SSE = 0.0;
+        for (uword i = 0; i < n; i++)
+        {
+            vec diff = weightedVals.col(i) - groupMeans(i);
+            SSE += accu(diff % diff);
+        }
+
+        double SSB = 0.0;
+        for (uword i = 0; i < n; i++)
+        {
+            double meanDiff = groupMeans(i) - yMean;
+            SSB += double(n) * meanDiff * meanDiff;
+        }
+        double MSB = SSB / double(n - 1);
+        double MSW = SSE / double(n * n - n);
+        if (!isfinite(MSB) || !isfinite(MSW) || MSW == 0.0)
+            continue;
+
+        double F = MSB / MSW;
+        if (!isfinite(F))
+            continue;
+
+        totalF += F;
+        valid = true;
+    }
+
+    criterion = valid ? -totalF : DBL_MAX;
+    mSpatialWeight.setWeight(backup);
+    return Status::Success;
+}
+
 void GWAverage::updateCalculator()
 {
     switch (mParallelType)