Gwaverage anova

include/gwmodelpp/GWAverage.h

@@ -3,6 +3,7 @@
 
 #include "SpatialMonoscaleAlgorithm.h"
 #include "IMultivariableAnalysis.h"
+#include "IBandwidthSelectable.h"
 #include "IParallelizable.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:
 
@@ -100,6 +101,21 @@ class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysi
      */
     void setQuantile(bool quantile) { mQuantile = quantile; }
 
+    /**
+     * @brief \~english Get whether bandwidth autoselect is enabled. \~chinese 获取是否启用带宽自动选择。
+     * 
+     * @return true \~english if autoselect bandwidth is enabled \~chinese 启用带宽自动选择
+     * @return false \~english if autoselect bandwidth is disabled \~chinese 不启用带宽自动选择
+     */
+    bool isAutoselectBandwidth() const { return mIsAutoselectBandwidth; }
+
+    /**
+     * @brief \~english Set whether bandwidth autoselect is enabled. \~chinese 设置是否启用带宽自动选择。
+     * 
+     * @param isAutoselect \~english whether autoselect bandwidth is enabled \~chinese 是否启用带宽自动选择
+     */
+    void setAutoselectBandwidth(bool isAutoselect) { mIsAutoselectBandwidth = isAutoselect; }
+
     /**
      * @brief \~english Get local mean on each sample. \~chinese 获取每个样本的局部均值。
      * 
@@ -171,6 +187,8 @@ class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysi
     void setVariables(const arma::mat& x) override { mX = x; }
 
     void run() override;
+    // Implement IBandwidthSelectable
+    Status getCriterion(BandwidthWeight* weight, double& criterion) override;
 
     void calibration(const arma::mat& locations, const arma::mat& x);
 
@@ -238,7 +256,8 @@ class GWAverage : public SpatialMonoscaleAlgorithm, public IMultivariableAnalysi
     arma::mat mLocalMedian;   //!< \~english Local medians \~chinese 局部中位数
     arma::mat mIQR;           //!< \~english Local interquartile ranges \~chinese 局部分位距
     arma::mat mQI;            //!< \~english Local quantile imbalances and coordinates \~chinese 局部分位数不平衡度
-
+    bool mIsAutoselectBandwidth = false; //!< \~english Whether bandwidth autoselect is enabled \~chinese 是否启用带宽自动选择
+    BandwidthCriterionList mBandwidthSelectionCriterionList; //!< \~english List of bandwidth selection criteria \~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/BandwidthSelector.cpp

@@ -1,12 +1,16 @@
 #include "BandwidthSelector.h"
+#include <iostream>
+#include <iomanip>
 
 using namespace std;
 using namespace gwm;
 
 Status BandwidthSelector::optimize(IBandwidthSelectable* instance)
 {
-    unique_ptr<BandwidthWeight> w1 { static_cast<BandwidthWeight*>(mBandwidth.clone().release()) };
-    unique_ptr<BandwidthWeight> w2 { static_cast<BandwidthWeight*>(mBandwidth.clone().release()) };
+    cerr << "[BandwidthSelector] optimize() starting: lower=" << mLower << ", upper=" << mUpper
+         << ", adaptive=" << std::boolalpha << mBandwidth->adaptive() << ", kernel=" << mBandwidth->kernel() << "\n";
+    BandwidthWeight* w1 = static_cast<BandwidthWeight*>(mBandwidth->clone());
+    BandwidthWeight* w2 = static_cast<BandwidthWeight*>(mBandwidth->clone());
     double xU = mUpper, xL = mLower;
     bool adaptBw = mBandwidth.adaptive();
     const double eps = 1e-4;
@@ -32,6 +36,8 @@ Status BandwidthSelector::optimize(IBandwidthSelectable* instance)
         mBandwidthCriterion[x1] = f1;
     if (f2 < DBL_MAX)
         mBandwidthCriterion[x2] = f2;
+    cerr << fixed << setprecision(6);
+    cerr << "[BandwidthSelector] initial x1=" << x1 << ", f1=" << f1 << "; x2=" << x2 << ", f2=" << f2 << "\n";
     double d1 = f2 - f1;
     double xopt = f1 < f2 ? x1 : x2;
     double ea = 100;
@@ -47,7 +53,8 @@ Status BandwidthSelector::optimize(IBandwidthSelectable* instance)
             w1->setBandwidth(x1);
             s1 = instance->getCriterion(w1, f1);
             if (f1 < DBL_MAX)
-                mBandwidthCriterion[x1] = f1;
+            mBandwidthCriterion[x1] = f1;
+            cerr << "[BandwidthSelector] iter=" << iter << " x1=" << x1 << " f1=" << f1 << " x2=" << x2 << " f2=" << f2 << " xopt=" << xopt << " d=" << d << "\n";
         }
         else
         {
@@ -63,7 +70,11 @@ Status BandwidthSelector::optimize(IBandwidthSelectable* instance)
         iter = iter + 1;
         xopt = (f1 < f2) ? x1 : x2;
         d1 = f2 - f1;
+        cerr << "[BandwidthSelector] iter=" << iter << " x1=" << x1 << " f1=" << f1 << " x2=" << x2 << " f2=" << f2 << " xopt=" << xopt << " d=" << d << "\n";
     }
+    delete w1;
+    delete w2;
+    cerr << "[BandwidthSelector] optimize completed: xopt=" << xopt << ", s1=" << static_cast<int>(s1) << ", s2=" << static_cast<int>(s2) << "\n";
     if (s1 == Status::Success && s2 == Status::Success)
     {
         mOptimisedBandwidth.setBandwidth(xopt);