From 75465ca1648b6b81dae4a7af7ba743e0d6bec251 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Fri, 26 Jul 2024 10:23:47 +0200
Subject: [PATCH] moved small_kendall to statistics

---
 rogueviz/sag/functions.cpp |  2 +-
 rogueviz/som/measures.cpp  | 40 ++----------------------------------
 rogueviz/statistics.cpp    | 42 +++++++++++++++++++++++++++++++++++++-
 3 files changed, 44 insertions(+), 40 deletions(-)

diff --git a/rogueviz/sag/functions.cpp b/rogueviz/sag/functions.cpp
index e7420d2f..0cb1d2f1 100644
--- a/rogueviz/sag/functions.cpp
+++ b/rogueviz/sag/functions.cpp
@@ -182,7 +182,7 @@ void optimize_sag_loglik_logistic() {
 
 void optimize_sag_loglik_match() {
   if(state &~ SS_WEIGHTED) return;
-  lsq::leastsquare_solver<2> lsqs;
+  stats::leastsquare_solver<2> lsqs;
 
   for(auto& ei: sagedges) {
     ld y = sagdist[sagid[ei.i]][sagid[ei.j]];
diff --git a/rogueviz/som/measures.cpp b/rogueviz/som/measures.cpp
index 76c651f8..a403220e 100644
--- a/rogueviz/som/measures.cpp
+++ b/rogueviz/som/measures.cpp
@@ -2,47 +2,11 @@
 // Copyright (C) 2011-2022 Tehora and Zeno Rogue, see 'hyper.cpp' for details
 
 #include "kohonen.h"
+#include "../statistics.cpp"
 
 namespace rogueviz {
 namespace measures {
 
-double kendall(const vector<pair<int, int>>& allp) {
-  int maxo = 0, maxe = 0;
-  for(const auto& a: allp) maxo = max(maxo, a.first), maxe = max(maxe, a.second);
-  maxo++; maxe++;
-  
-  if(maxo >= 256 || maxe >= 256) {
-    println(hlog, "more than 256!\n");
-    exit(1);
-    }
-  int cnt[256][256];
-  for(int a=0; a<maxo; a++)
-  for(int b=0; b<maxe; b++)
-    cnt[a][b] = 0;
-  for(const auto& a: allp) cnt[a.first][a.second]++;
-
-  // int i1 = 0, i2 = 0;
-  int K = isize(allp);
-  // allp.emplace_back(maxo, maxe);
-  vector<int> counts(maxe, 0);
-  vector<int> totals(maxe);
-  double tau = 0;
-  for(int i=0; i<maxo; i++) {
-    totals[0] = 0;
-    for(int ii=1; ii<maxe; ii++)
-      totals[0] -= counts[ii];
-    for(int ii=1; ii<maxe; ii++)
-      totals[ii] = totals[ii-1] + counts[ii] + counts[ii-1];
-    
-    for(int b=0; b<maxe; b++) {
-      tau += totals[b] * 1. * cnt[i][b];
-      counts[b] += cnt[i][b];
-      }
-    }
-  ld par = (K * (K-1.) / 2);
-  return tau / par;
-  }
-
 vector<pair<int, int>> recreate_topology(const vector<int>& mapp, const vector<pair<int, int> >& edges) {
 
   auto cmapp = mapp;
@@ -117,7 +81,7 @@ ld evaluate_measure(manidata& emb, manidata& orig, vector<int>& mapp, vector<pai
     for(int i=0; i<No; i++) if(mapp[i] != -1)
     for(int j=0; j<i; j++) if(mapp[j] != -1)
       allp.emplace_back(diso[i][j], dise[mapp[i]][mapp[j]]);
-    energy = kendall(allp);
+    energy = stats::small_kendall<256, 256>(allp);
     }
   else if(k == 3) {
     vector<bool> empty(Ne, true);
diff --git a/rogueviz/statistics.cpp b/rogueviz/statistics.cpp
index 2d735b1c..17cd9cf7 100644
--- a/rogueviz/statistics.cpp
+++ b/rogueviz/statistics.cpp
@@ -1,9 +1,12 @@
+#ifndef _STATISTICS_CPP_
+#define _STATISTICS_CPP_
+
 #include <vector>
 #include <cstdio>
 #include <cmath>
 #include <array>
 
-namespace lsq {
+namespace stats {
 
 using namespace std;
 
@@ -98,4 +101,41 @@ template<size_t N> struct leastsquare_solver {
     }
   };
 
+template<size_t dim1, size_t dim2> double small_kendall(const vector<pair<int, int>>& allp) {
+  int maxo = 0, maxe = 0;
+  for(const auto& a: allp) maxo = max(maxo, a.first), maxe = max(maxe, a.second);
+  maxo++; maxe++;
+
+  if(maxo >= dim1 || maxe >= dim2)
+    throw hr::hr_exception("small_kendall limit exceeded");
+  int cnt[dim1][dim2];
+  for(int a=0; a<maxo; a++)
+  for(int b=0; b<maxe; b++)
+    cnt[a][b] = 0;
+  for(const auto& a: allp) cnt[a.first][a.second]++;
+
+  // int i1 = 0, i2 = 0;
+  int K = hr::isize(allp);
+  // allp.emplace_back(maxo, maxe);
+  vector<int> counts(maxe, 0);
+  vector<int> totals(maxe);
+  double tau = 0;
+  for(int i=0; i<maxo; i++) {
+    totals[0] = 0;
+    for(int ii=1; ii<maxe; ii++)
+      totals[0] -= counts[ii];
+    for(int ii=1; ii<maxe; ii++)
+      totals[ii] = totals[ii-1] + counts[ii] + counts[ii-1];
+
+    for(int b=0; b<maxe; b++) {
+      tau += totals[b] * 1. * cnt[i][b];
+      counts[b] += cnt[i][b];
+      }
+    }
+  double par = (K * (K-1.) / 2);
+  return tau / par;
+  }
+
+
 }
+#endif