From 8521405c611c353406ca02286d00d99fff44ca22 Mon Sep 17 00:00:00 2001 From: Egor Kotov Date: Mon, 20 Apr 2026 14:21:11 +0200 Subject: [PATCH] Implement Local Spatial Heteroscedasticity (LOSH) - Added native C++ implementation of the LOSH statistic matching spdep::LOSH. - Implemented Multi-threaded Empirical Permutations for LOSH inference. - Registered standard cluster IDs and placeholder labels. --- gda_sa.cpp | 26 ++++- gda_sa.h | 18 +++ sa/UniLOSH.cpp | 306 +++++++++++++++++++++++++++++++++++++++++++++++++ sa/UniLOSH.h | 61 ++++++++++ 4 files changed, 408 insertions(+), 3 deletions(-) create mode 100644 sa/UniLOSH.cpp create mode 100644 sa/UniLOSH.h diff --git a/gda_sa.cpp b/gda_sa.cpp index fa2c5d9..e58a2e9 100644 --- a/gda_sa.cpp +++ b/gda_sa.cpp @@ -10,6 +10,7 @@ #include "sa/UniGstar.h" #include "sa/UniJoinCount.h" #include "sa/UniLocalMoran.h" +#include "sa/UniLOSH.h" #include "sa/MultiGeary.h" #include "sa/MultiJoinCount.h" #include "sa/BatchLocalMoran.h" @@ -201,9 +202,9 @@ LISA *gda_localjoincount(GeoDaWeight *w, } LISA *gda_localmultijoincount(GeoDaWeight *w, - const std::vector > &data, - const std::vector > &undefs, - double significance_cutoff, int nCPUs, int perm, const std::string& perm_method, int last_seed) + const std::vector > &data, + const std::vector > &undefs, + double significance_cutoff, int nCPUs, int perm, const std::string& perm_method, int last_seed) { if (w == 0) return 0; @@ -214,6 +215,25 @@ LISA *gda_localmultijoincount(GeoDaWeight *w, return jc; } +LISA *gda_locallosh(GeoDaWeight *w, + const std::vector &data, + const std::vector &undefs, + double significance_cutoff, int nCPUs, int perm, const std::string& perm_method, int last_seed, + double a) +{ + if (w == 0) + return 0; + + int num_obs = w->num_obs; + + std::vector copy_undefs = undefs; + if (copy_undefs.empty()) + { + copy_undefs.resize(num_obs, false); + } + UniLOSH *losh = new UniLOSH(num_obs, w, data, copy_undefs, significance_cutoff, nCPUs, perm, perm_method, last_seed, a); + return losh; +} double gda_fdr(LISA *lisa, double current_p) { if (lisa == 0) diff --git a/gda_sa.h b/gda_sa.h index 4b5d58a..5866ae4 100644 --- a/gda_sa.h +++ b/gda_sa.h @@ -108,6 +108,24 @@ LISA *gda_localjoincount(GeoDaWeight *w, const std::vector &data, LISA *gda_localmultijoincount(GeoDaWeight *w, const std::vector > &data, const std::vector > &undefs, double significance_cutoff, int nCPUs, int permutations, const std::string& permutation_method, int last_seed_used); + +/** + * + * @param w + * @param data + * @param undefs + * @param significance_cutoff + * @param nCPUs + * @param permutations + * @param permutation_method + * @param last_seed_used + * @param a + * @return + */ +LISA *gda_locallosh(GeoDaWeight *w, const std::vector &data, + const std::vector &undefs, double significance_cutoff, + int nCPUs, int permutations, const std::string& permutation_method, int last_seed_used, + double a); /** * * @param w diff --git a/sa/UniLOSH.cpp b/sa/UniLOSH.cpp new file mode 100644 index 0000000..8b7c22b --- /dev/null +++ b/sa/UniLOSH.cpp @@ -0,0 +1,306 @@ +// +// Created by Agent on 2026-04-19. +// + +#include +#include "../GeoDaSet.h" +#include "../GenUtils.h" +#include "../weights/GeodaWeight.h" +#include "UniLOSH.h" + +UniLOSH::~UniLOSH() { + +} + +UniLOSH::UniLOSH(int num_obs, + GeoDaWeight *w, + const std::vector &_data, + const std::vector &_undefs, + double significance_cutoff, + int _nCPUs, int _perm, + const std::string& _permutation_method, uint64_t _last_seed, + double _a) +: LISA(num_obs, w, _undefs, significance_cutoff, _nCPUs, _perm, _permutation_method, _last_seed), + CLUSTER_NOT_SIG(0), + CLUSTER_HETEROGENEOUS(1), + CLUSTER_HOMOGENEOUS(2), + CLUSTER_UNDEFINED(5), + CLUSTER_NEIGHBORLESS(6), + data(_data), + a(_a), + mean_e(0) +{ + labels.push_back("Not significant"); // 0 + labels.push_back("Heterogeneous"); // 1 + labels.push_back("Homogeneous"); // 2 + labels.push_back(""); // 3 + labels.push_back(""); // 4 + labels.push_back("Undefined"); // 5 + labels.push_back("Isolated"); // 6 + + colors.push_back("#eeeeee"); // 0 + colors.push_back("#FF00FF"); // 1 + colors.push_back("#00FFFF"); // 2 + colors.push_back(""); // 3 + colors.push_back(""); // 4 + colors.push_back("#464646"); // 5 + colors.push_back("#999999"); // 6 + + local_mean.resize(num_obs, 0.0); + local_residuals.resize(num_obs, 0.0); + + if (_perm == 0) { + calc_significances = false; + } + + Run(); +} + +void UniLOSH::ComputeLoalSA() { + mean_e = 0.0; + int valid_obs = 0; + + // Pass 1: Compute spatially lagged means and local residuals + for (int i=0; iIsMasked(i) == false) { + lag_vec[i] = 0; + lisa_vec[i] = 0; + cluster_vec[i] = CLUSTER_UNDEFINED; + } else { + if (weights->GetNbrSize(i) == 0) { + cluster_vec[i] = CLUSTER_NEIGHBORLESS; + } else { + const std::vector& nbrs = weights->GetNeighbors(i); + const std::vector& nbr_w = weights->GetNeighborWeights(i); + int num_nbrs = nbrs.size(); + + double sp_lag = 0.0; + double w_sum = 0.0; + + for (int j=0; j 0 && row_standardize) { + sp_lag /= w_sum; + } + + local_mean[i] = sp_lag; + double diff = std::abs(data[i] - sp_lag); + if (a == 2.0) { + local_residuals[i] = diff * diff; + } else if (a == 1.0) { + local_residuals[i] = diff; + } else { + local_residuals[i] = std::pow(diff, a); + } + + mean_e += local_residuals[i]; + valid_obs++; + } + } + } + + if (valid_obs > 0) { + mean_e /= valid_obs; + } + if (mean_e == 0.0) { + mean_e = 1.0; + } + + // Pass 2: Compute spatially lagged residuals (Hi) + for (int i=0; iIsMasked(i) == false || weights->GetNbrSize(i) == 0) continue; + + const std::vector& nbrs = weights->GetNeighbors(i); + const std::vector& nbr_w = weights->GetNeighborWeights(i); + int num_nbrs = nbrs.size(); + + double lag_e = 0.0; + double w_sum = 0.0; + + for (int j=0; j 0 && row_standardize) { + lag_e /= w_sum; + } + + lag_vec[i] = lag_e; + lisa_vec[i] = lag_e / mean_e; + + // Initialize clusters (will be filtered by significance later in GetClusterIndicators) + if (lisa_vec[i] > 1.0) { + cluster_vec[i] = CLUSTER_HETEROGENEOUS; + } else { + cluster_vec[i] = CLUSTER_HOMOGENEOUS; + } + } +} + +void UniLOSH::CalcPseudoP_range(int obs_start, int obs_end, uint64_t seed_start) +{ + GeoDaSet workPermutation(num_obs); + int max_rand = num_obs - 1; + + for (int cnt = obs_start; cnt <= obs_end; cnt++) { + if (undefs[cnt] || weights->IsMasked(cnt) == false) { + sig_cat_vec[cnt] = 6; // undefined + continue; + } + + int numNeighbors = weights->GetNbrSize(cnt); + if (numNeighbors == 0) { + sig_cat_vec[cnt] = 5; // neighborless cat + continue; + } + + std::vector nbr_w = weights->GetNeighborWeights(cnt); + std::vector permutedSA(permutations, 0.0); + + for (int perm = 0; perm < permutations; perm++) { + int rand = 0, newRandom; + double rng_val; + while (rand < numNeighbors) { + rng_val = Gda::ThomasWangHashDouble(seed_start++) * max_rand; + newRandom = (int)(rng_val < 0.0 ? std::ceil(rng_val - 0.5) : std::floor(rng_val + 0.5)); + + if (newRandom != cnt && !workPermutation.Belongs(newRandom) && weights->GetNbrSize(newRandom) > 0) { + workPermutation.Push(newRandom); + rand++; + } + } + + double permutedLag = 0.0; + double w_sum = 0.0; + for (int cp = 0; cp < numNeighbors; cp++) { + int nb = workPermutation.Pop(); + if (!undefs[nb]) { + permutedLag += local_residuals[nb] * nbr_w[cp]; + w_sum += nbr_w[cp]; + } + } + + if (w_sum > 0 && row_standardize) { + permutedLag /= w_sum; + } + permutedSA[perm] = permutedLag / mean_e; + } + + uint64_t countLarger = CountLargerSA(cnt, permutedSA); + double _sigLocal = (countLarger + 1.0) / (permutations + 1); + + if (_sigLocal <= 0.0001) sig_cat_vec[cnt] = 4; + else if (_sigLocal <= 0.001) sig_cat_vec[cnt] = 3; + else if (_sigLocal <= 0.01) sig_cat_vec[cnt] = 2; + else if (_sigLocal <= 0.05) sig_cat_vec[cnt] = 1; + else sig_cat_vec[cnt] = 0; + + sig_local_vec[cnt] = _sigLocal; + } +} + +void UniLOSH::PermCalcPseudoP_range(int obs_start, int obs_end, uint64_t seed_start) +{ + for (int cnt = obs_start; cnt <= obs_end; cnt++) { + if (undefs[cnt] || weights->IsMasked(cnt) == false) { + sig_cat_vec[cnt] = 6; // undefined + continue; + } + + int numNeighbors = weights->GetNbrSize(cnt); + if (numNeighbors == 0) { + sig_cat_vec[cnt] = 5; // neighborless cat + continue; + } + + std::vector nbr_w = weights->GetNeighborWeights(cnt); + std::vector permutedSA(permutations, 0.0); + + for (size_t perm = 0; perm < (size_t)permutations; perm++) { + double permutedLag = 0.0; + double w_sum = 0.0; + int* permNeighbors = perm_table[perm]; + for (int cp = 0; cp < numNeighbors; cp++) { + int nb = permNeighbors[cp]; + if (nb >= cnt) nb = nb + 1; + if (!undefs[nb]) { + permutedLag += local_residuals[nb] * nbr_w[cp]; + w_sum += nbr_w[cp]; + } + } + + if (w_sum > 0 && row_standardize) { + permutedLag /= w_sum; + } + permutedSA[perm] = permutedLag / mean_e; + } + + uint64_t countLarger = CountLargerSA(cnt, permutedSA); + double _sigLocal = (countLarger + 1.0) / (permutations + 1); + + if (_sigLocal <= 0.0001) sig_cat_vec[cnt] = 4; + else if (_sigLocal <= 0.001) sig_cat_vec[cnt] = 3; + else if (_sigLocal <= 0.01) sig_cat_vec[cnt] = 2; + else if (_sigLocal <= 0.05) sig_cat_vec[cnt] = 1; + else sig_cat_vec[cnt] = 0; + + sig_local_vec[cnt] = _sigLocal; + } +} + +void UniLOSH::PermLocalSA(int cnt, int perm, const std::vector &permNeighbors, + std::vector& permutedSA) { + // Should not be called because we override CalcPseudoP_range +} + +void UniLOSH::PermLocalSA(int cnt, int perm, int numNeighbors, const int* permNeighbors, + std::vector& permutedSA) { + // Should not be called because we override PermCalcPseudoP_range +} + +uint64_t UniLOSH::CountLargerSA(int cnt, const std::vector& permutedSA) +{ + uint64_t countLarger = 0; + for (int i=0; i= lisa_vec[cnt]) { + countLarger += 1; + } + } + + if (permutations - countLarger <= countLarger) { + countLarger = permutations - countLarger; + } + return countLarger; +} + +std::vector UniLOSH::GetClusterIndicators() { + std::vector clusters(num_obs); + double cutoff = GetSignificanceCutoff(); + for (int i=0; i cutoff && + (const unsigned long)cluster_vec[i] != CLUSTER_UNDEFINED && + (const unsigned long)cluster_vec[i] != CLUSTER_NEIGHBORLESS) + { + clusters[i] = CLUSTER_NOT_SIG; + } else { + clusters[i] = cluster_vec[i]; + } + } + return clusters; +} + +std::vector UniLOSH::GetLocalMean() { + return local_mean; +} + +std::vector UniLOSH::GetLocalResiduals() { + return local_residuals; +} diff --git a/sa/UniLOSH.h b/sa/UniLOSH.h new file mode 100644 index 0000000..c528153 --- /dev/null +++ b/sa/UniLOSH.h @@ -0,0 +1,61 @@ +// +// Created by Agent on 2026-04-19. +// + +#ifndef GEODA_UNILOSH_H +#define GEODA_UNILOSH_H + +#include + +#include "LISA.h" + +class GeoDaWeight; + +class UniLOSH : public LISA { + + const unsigned long CLUSTER_NOT_SIG; + const unsigned long CLUSTER_HETEROGENEOUS; + const unsigned long CLUSTER_HOMOGENEOUS; + const unsigned long CLUSTER_UNDEFINED; + const unsigned long CLUSTER_NEIGHBORLESS; + +public: + UniLOSH(int num_obs, + GeoDaWeight* w, + const std::vector& data, + const std::vector& undefs, + double significance_cutoff, + int nCPUs, int permutations, + const std::string& _permutation_method, + uint64_t last_seed_used, + double a = 2.0); + + virtual ~UniLOSH(); + + virtual void ComputeLoalSA() ; + + virtual void CalcPseudoP_range(int obs_start, int obs_end, uint64_t seed_start); + + virtual void PermCalcPseudoP_range(int obs_start, int obs_end, uint64_t seed_start); + + virtual void PermLocalSA(int cnt, int perm, const std::vector &permNeighbors, std::vector& permutedSA); + + virtual void PermLocalSA(int cnt, int perm, int numNeighbors, const int* permNeighbors, std::vector& permutedSA); + + virtual uint64_t CountLargerSA(int cnt, const std::vector& permutedSA); + + virtual std::vector GetClusterIndicators(); + + std::vector GetLocalMean(); + std::vector GetLocalResiduals(); + +protected: + std::vector data; + double a; + double mean_e; + std::vector local_mean; + std::vector local_residuals; +}; + + +#endif //GEODA_UNILOSH_H