linux build script changes. Workign on debugging trace with matrix fu…

…nctions

include/_lanczos/lanczos.h

@@ -216,13 +216,19 @@ void sl_trace(
   const auto trace_f = [lanczos_degree, &sf, &estimates](int i, [[maybe_unused]] F* q, [[maybe_unused]] F* Q, F* nodes, F* weights){
     Eigen::Map< VectorF > nodes_v(nodes, lanczos_degree, 1);     // no-op
     Eigen::Map< VectorF > weights_v(weights, lanczos_degree, 1); // no-op
-    nodes_v.unaryExpr(sf);
+    for (int c = 0; c < lanczos_degree; ++c){
+      std::cout << node_v[c] << " -> " << sf(node_v[c]) << std::endl; 
+      node_v[c] = sf(node_v[c]);
+    }
+    // nodes_v.unaryExpr(sf);
     estimates[i] = (nodes_v * weights_v).sum();
   };
 
   // Type-erased function since the call is cheap
   std::function< bool (int) > early_stop; 
-  if (use_CLT){
+  if (atol == 0.0 && rtol == 0.0){
+    early_stop = [](int i) -> bool { return false; };
+  } else if (use_CLT){
     // Parameterize when to stop using either the CLT over the given confidence level or 
     // This runs in critical section of the SLQ, so we can depend sequential execution (but i will vary!)
     // See: https://math.stackexchange.com/questions/102978/incremental-computation-of-standard-deviation
@@ -239,7 +245,7 @@ void sl_trace(
       vr_est = L * vr_pre + denom * std::pow(estimates[i] - mu_pre, 2); // update sample variance
       mu_pre = mu_est;
       vr_pre = vr_est;
-      if (n < 2){
+      if (n < 3){
         return false; 
       } else {
         const auto sd_est = std::sqrt(vr_est);

include/_random_generator/random_concepts.h

@@ -42,7 +42,7 @@ struct Random64Engine : public Random64EngineConcept {
   void seed(std::seed_seq& S) { rng.seed(S); };
   uint64_t operator()() { return rng(); }
   size_t state_size() const {
-    if constexpr (Stateful64Engine< T >){ return T::state_size; } 
+    if constexpr (Stateful64Engine< T >){ return std::max(T::state_size, size_t(1)); } 
     return 1; 
   }
 };

src/primate/stats.py

@@ -13,7 +13,6 @@ def suggest_nv_trace(p: float, eps: float, f: str = "identity", dist: str = "rad
   else: 
     raise NotImplementedError("TODO")
 
-
 # def suggest_nv_mf(p: float, eps: float, max_f: float, k: int, e_min: float, e_max: float, n: int, dist: str = "rademacher"):
 #   k = e_max / e_min 
 #   rho = (np.sqrt(k) + 1)/(np.sqrt(k) - 1)

src/primate/trace.py

@@ -21,7 +21,7 @@ def sl_trace (
   orth: int = 0,
   confidence: float = 0.95,
   pdf: str = "rademacher",
-  rng: str = "pcg",
+  rng: str = "lcg",
   seed: int = -1,
   num_threads: int = 0,
   verbose: bool = False,
@@ -298,10 +298,6 @@ def sl_gauss(
 
 
 
-
-
-
-
 # def trace_est() -> int:
 #   return _trace.apply_smoothstep(1, 2)
 

tests/test_slq.py

@@ -71,12 +71,23 @@ def test_slq_fixed():
   assert np.isclose(A.trace() - tr_est, 0.0, atol=threshold)
 
 def test_slq_trace():
-  from primate.trace import sl_trace, _lanczos
+  from primate.trace import sl_trace
+  np.random.seed(1234)
+  n = 25
+  A = csc_array(symmetric(n), dtype=np.float32)
+  tr_est = sl_trace(A, maxiter = 200, num_threads=1, seed=-1)
+  assert len(tr_est) == 200
+  assert np.all(~np.isclose(tr_est, 0.0))
+  assert np.isclose(np.mean(tr_est), A.trace(), atol=1.0)
+
+def test_slq_trace_multithread():
+  from primate.trace import sl_trace
   np.random.seed(1234)
   n = 25
   A = csc_array(symmetric(n), dtype=np.float32)
-  tr_est = sl_trace(A, maxiter = 200, num_threads=1)
+  tr_est = sl_trace(A, maxiter = 200, atol=0.0, num_threads=6)
   assert len(tr_est) == 200
+  assert np.all(~np.isclose(tr_est, 0.0))
   assert np.isclose(np.mean(tr_est), A.trace(), atol=1.0)
 
 def test_slq_trace_clt_atol():
@@ -101,4 +112,13 @@ def test_slq_trace_clt_atol():
   assert converged_online == converged_ind, "SLQ not converging at correct index!"
 
 
+def test_slq_trace_f():
+  from primate.trace import sl_trace, _lanczos
+  np.random.seed(1234)
+  n = 30
+  A = csc_array(symmetric(n), dtype=np.float32)
+
+
+  np.isclose(np.mean(sl_trace(A, fun="numrank")), np.linalg.matrix_rank(A.todense())
+
 

todo/examples/xtrace.py

@@ -84,7 +84,7 @@ def xtrace(A, error_atol: float = 0.1, error_rtol: float = 1e-6, nv: int = 1, co
 		ns = min(A.shape[1] - Y.shape[1], ns)
 		NewOm = sample_isotropic(n, ns, method=method)
 
-		## If the sampled vectors have a lot of linear dependence  won't be (numerically) span a large enough subspace
+		## If the sampled vectors have a lot of linear dependence won't be (numerically) span a large enough subspace
 		## to permit sufficient exploration of the eigen-space, so we optionally re-sample based on a loose upper bound
 		## Based on: https://math.stackexchange.com/questions/1191503/conditioning-of-triangular-matrices
 		cond_numb_bound = np.inf

tools/cibw_linux.sh

@@ -1,20 +1,24 @@
 #!/usr/bin/env bash
 
-if [ -n "$(command -v yum)" ]; then 
-  dnf check-update
-  dnf update
-  yum install llvm-toolset
-  yum install openblas
+if [ -n "$(command -v yum)" ]; then
+  cat /etc/*-release 
+  yum update -y 
+  yum install centos-release-scl
+  yum install llvm-toolset-7
+  scl enable llvm-toolset-7
+  yum install -y openblas
   yum install -y python3.9
-  yum install python39-devel
+  yum install -y python39-devel
   alias python=python3.9
 elif [ -n "$(command -v apt-get)" ]; then 
+  cat /etc/*-release
   sudo apt-get update -y
   sudo apt-get install -y clang libomp-dev
   sudo apt-get update -y
   sudo apt-get install -y libopenblas-dev
   sudo apt-get install -y python3-dev
 elif [ -n "$(command -v apk)" ]; then 
+  cat /etc/*-release
   apk update
   apk add clang # looks like Alpine only supports up to clang 10 ?
   # apk add openmp