diff --git a/include/_lanczos/lanczos.h b/include/_lanczos/lanczos.h index 493baa3..4ce5063 100644 --- a/include/_lanczos/lanczos.h +++ b/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); diff --git a/include/_random_generator/random_concepts.h b/include/_random_generator/random_concepts.h index dc7a973..0a41644 100644 --- a/include/_random_generator/random_concepts.h +++ b/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; } }; diff --git a/src/primate/stats.py b/src/primate/stats.py index d8fdcd2..d63e277 100644 --- a/src/primate/stats.py +++ b/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) diff --git a/src/primate/trace.py b/src/primate/trace.py index 9a555a6..e5bddcd 100644 --- a/src/primate/trace.py +++ b/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) diff --git a/tests/test_slq.py b/tests/test_slq.py index 7cb543b..ba7fbcc 100644 --- a/tests/test_slq.py +++ b/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()) + diff --git a/todo/examples/xtrace.py b/todo/examples/xtrace.py index 8cf18d8..a582ef7 100644 --- a/todo/examples/xtrace.py +++ b/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 diff --git a/todo/examples/xtrace_pythran.py b/todo/examples/xtrace_pythran.py deleted file mode 100644 index e69de29..0000000 diff --git a/tools/cibw_linux.sh b/tools/cibw_linux.sh index 1a3d14f..94b93c4 100644 --- a/tools/cibw_linux.sh +++ b/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