From 660893ee5dc81f0e293fd50ba1d5586c1d181de5 Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 08:28:05 -0700 Subject: [PATCH 1/9] Update .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index d011e93..208450c 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ .DS_Store __pycache__/ *.pyc +*.png +DiffusionGeometry/ \ No newline at end of file From 388e076083c6fdab846e68a8ef9ba5954b521c4d Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 09:20:33 -0700 Subject: [PATCH 2/9] Add Hodge parity harness and align circular mode defaults --- .gitignore | 5 +- benches/bench_dod.cpp | 4 +- benches/bench_pipelines.cpp | 5 +- include/igneous/data/topology.hpp | 420 ++++++++++++++++++----- include/igneous/ops/geometry.hpp | 108 ++++-- include/igneous/ops/hodge.hpp | 132 +++++-- include/igneous/ops/spectral.hpp | 154 ++++++--- notes/hodge/journal.md | 100 ++++++ notes/hodge/results/.gitkeep | 0 scripts/hodge/compare_hodge_outputs.py | 456 +++++++++++++++++++++++++ scripts/hodge/generate_input_torus.py | 36 ++ scripts/hodge/run_cpp_hodge.sh | 40 +++ scripts/hodge/run_parity_round.sh | 115 +++++++ scripts/hodge/run_reference_hodge.py | 174 ++++++++++ scripts/hodge/setup_reference_env.sh | 30 ++ src/main_diffusion.cpp | 2 +- src/main_hodge.cpp | 425 +++++++++++++++++++---- src/main_spectral.cpp | 2 +- tests/test_ops_hodge.cpp | 8 +- tests/test_ops_spectral_geometry.cpp | 3 +- tests/test_topology_diffusion.cpp | 6 +- 21 files changed, 1962 insertions(+), 263 deletions(-) create mode 100644 notes/hodge/journal.md create mode 100644 notes/hodge/results/.gitkeep create mode 100755 scripts/hodge/compare_hodge_outputs.py create mode 100755 scripts/hodge/generate_input_torus.py create mode 100755 scripts/hodge/run_cpp_hodge.sh create mode 100755 scripts/hodge/run_parity_round.sh create mode 100755 scripts/hodge/run_reference_hodge.py create mode 100755 scripts/hodge/setup_reference_env.sh diff --git a/.gitignore b/.gitignore index 208450c..e79e49a 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,7 @@ __pycache__/ *.pyc *.png -DiffusionGeometry/ \ No newline at end of file +DiffusionGeometry/ +notes/hodge/.venv_ref/ +notes/hodge/results/* +!notes/hodge/results/.gitkeep diff --git a/benches/bench_dod.cpp b/benches/bench_dod.cpp index 985e187..99b5ef2 100644 --- a/benches/bench_dod.cpp +++ b/benches/bench_dod.cpp @@ -83,7 +83,7 @@ static DiffusionMesh make_diffusion_cloud(size_t n_points) { } mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), 0.05f, 32}); + mesh.geometry.z_span(), 32}); return mesh; } @@ -127,7 +127,7 @@ static void bench_diffusion_build(benchmark::State &state) { for (auto _ : state) { mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), 0.05f, 32}); + mesh.geometry.z_span(), 32}); benchmark::DoNotOptimize(mesh.topology.markov_values.size()); } } diff --git a/benches/bench_pipelines.cpp b/benches/bench_pipelines.cpp index 507e660..f0bd410 100644 --- a/benches/bench_pipelines.cpp +++ b/benches/bench_pipelines.cpp @@ -83,9 +83,10 @@ int compute_max_y_vertex(const DiffusionMesh &mesh) { return max_y_idx; } -void build_diffusion_topology(DiffusionMesh &mesh, float bandwidth, int k_neighbors) { +void build_diffusion_topology(DiffusionMesh &mesh, float /*bandwidth*/, + int k_neighbors) { mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), bandwidth, k_neighbors}); + mesh.geometry.z_span(), k_neighbors}); } void bench_pipeline_diffusion_main(benchmark::State &state) { diff --git a/include/igneous/data/topology.hpp b/include/igneous/data/topology.hpp index dee967b..f4c3f5e 100644 --- a/include/igneous/data/topology.hpp +++ b/include/igneous/data/topology.hpp @@ -1,12 +1,14 @@ #pragma once #include +#include #include #include #include #include #include #include +#include #include #include #include @@ -318,18 +320,36 @@ struct DiffusionTopology { std::span x; std::span y; std::span z; - float bandwidth = 0.01f; int k_neighbors = 32; + int knn_bandwidth = 8; + float bandwidth_variability = -0.5f; + float c = 0.0f; + bool use_mean_centres = true; }; Eigen::VectorXf mu; Eigen::MatrixXf eigen_basis; + Eigen::VectorXf local_bandwidths; + Eigen::VectorXf symmetric_row_sums; + Eigen::MatrixXf immersion_coords; + bool use_mean_centres = true; + int knn_k = 0; + + // Dense kNN row storage used by the reference-style diffusion build. + std::vector knn_indices; + std::vector knn_distances; + std::vector knn_kernel; // Direct CSR view for hot diffusion operators. std::vector markov_row_offsets; std::vector markov_col_indices; std::vector markov_values; + // Symmetric kernel used to reproduce reference spectral basis construction. + std::vector symmetric_row_offsets; + std::vector symmetric_col_indices; + std::vector symmetric_values; + [[nodiscard]] size_t num_primitives() const { if (markov_row_offsets.empty()) { return 0; @@ -344,9 +364,20 @@ struct DiffusionTopology { core::gpu::invalidate_markov_cache(this); mu.resize(0); eigen_basis.resize(0, 0); + local_bandwidths.resize(0); + symmetric_row_sums.resize(0); + immersion_coords.resize(0, 0); + use_mean_centres = true; + knn_k = 0; + knn_indices.clear(); + knn_distances.clear(); + knn_kernel.clear(); markov_row_offsets.clear(); markov_col_indices.clear(); markov_values.clear(); + symmetric_row_offsets.clear(); + symmetric_col_indices.clear(); + symmetric_values.clear(); } struct PointCloudAdaptor { @@ -373,6 +404,88 @@ struct DiffusionTopology { nanoflann::L2_Simple_Adaptor, PointCloudAdaptor, 3>; + static std::vector build_epsilons() { + std::vector epsilons; + epsilons.reserve(80); + for (float exponent = -10.0f; exponent < 10.0f; exponent += 0.25f) { + epsilons.push_back(std::pow(2.0f, exponent)); + } + return epsilons; + } + + static std::pair tune_kernel(const std::vector &entries, + int n, int k, + const std::vector &epsilons) { + const int e_count = static_cast(epsilons.size()); + std::vector averages(static_cast(e_count), 0.0f); + const float inv_nk = 1.0f / static_cast(std::max(1, n * k)); + + for (int e = 0; e < e_count; ++e) { + const float epsilon = std::max(epsilons[static_cast(e)], 1e-12f); + float sum = 0.0f; + for (float value : entries) { + sum += std::exp(-value / epsilon); + } + averages[static_cast(e)] = sum * inv_nk; + } + + int best_idx = 0; + float best_criterion = -std::numeric_limits::infinity(); + for (int i = 0; i + 1 < e_count; ++i) { + const float avg0 = std::max(averages[static_cast(i)], 1e-30f); + const float avg1 = std::max(averages[static_cast(i + 1)], 1e-30f); + const float eps0 = std::max(epsilons[static_cast(i)], 1e-12f); + const float eps1 = std::max(epsilons[static_cast(i + 1)], 1e-12f); + const float denom = std::log(eps1) - std::log(eps0); + if (std::abs(denom) < 1e-12f) { + continue; + } + const float criterion = (std::log(avg1) - std::log(avg0)) / denom; + if (criterion > best_criterion) { + best_criterion = criterion; + best_idx = i; + } + } + + const float epsilon = + std::max(epsilons[static_cast(best_idx)], 1e-12f); + const float dim = 2.0f * best_criterion; + return {epsilon, dim}; + } + + static void compute_local_bandwidths(const std::vector &nbr_distances, + int n, int k, int knn_bandwidth, + Eigen::VectorXf &bandwidths_out) { + bandwidths_out.resize(n); + bandwidths_out.setOnes(); + if (k <= 1) { + return; + } + + const int end = std::max(2, std::min(knn_bandwidth, k)); + const int count = std::max(1, end - 1); + for (int i = 0; i < n; ++i) { + const size_t base = static_cast(i) * static_cast(k); + float accum = 0.0f; + for (int idx = 1; idx < end; ++idx) { + const float d = nbr_distances[base + static_cast(idx)]; + accum += d * d; + } + const float rms = std::sqrt(accum / static_cast(count)); + bandwidths_out[i] = std::max(rms, 1e-8f); + } + } + + static float vector_median(std::vector values) { + if (values.empty()) { + return 1.0f; + } + const size_t mid = values.size() / 2; + std::nth_element(values.begin(), values.begin() + static_cast(mid), + values.end()); + return values[mid]; + } + void build(Input input) { core::gpu::invalidate_markov_cache(this); const size_t n_verts = input.x.size(); @@ -381,127 +494,262 @@ struct DiffusionTopology { return; } + const int n = static_cast(n_verts); const int k_neighbors = std::max(1, std::min(input.k_neighbors, static_cast(n_verts))); - const float t_sq = std::max(input.bandwidth, 1e-8f); + const int knn_bandwidth = std::max(2, std::min(input.knn_bandwidth, k_neighbors)); + use_mean_centres = input.use_mean_centres; + knn_k = k_neighbors; PointCloudAdaptor adaptor{input.x, input.y, input.z}; KDTree tree(3, adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(32)); tree.buildIndex(); - markov_row_offsets.assign(n_verts + 1, 0); - mu.resize(static_cast(n_verts)); - mu.setZero(); const size_t row_stride = static_cast(k_neighbors); - const size_t row_capacity = n_verts * row_stride; + const size_t row_capacity = static_cast(n) * row_stride; - std::vector row_counts(n_verts, 0); - std::vector row_cols(row_capacity, 0); - std::vector row_values(row_capacity, 0.0f); + knn_indices.assign(row_capacity, 0); + knn_distances.assign(row_capacity, 0.0f); + knn_kernel.assign(row_capacity, 0.0f); core::parallel_for_index( - 0, static_cast(n_verts), + 0, n, [&](int row_idx) { const size_t i = static_cast(row_idx); const size_t base = i * row_stride; thread_local std::vector ret_index; thread_local std::vector out_dist_sqr; - if (ret_index.size() < static_cast(k_neighbors)) { - ret_index.resize(static_cast(k_neighbors)); + ret_index.assign(static_cast(k_neighbors), 0u); + out_dist_sqr.assign(static_cast(k_neighbors), 0.0f); + + const float query_pt[3] = {input.x[i], input.y[i], input.z[i]}; + const size_t found = tree.knnSearch(query_pt, static_cast(k_neighbors), + ret_index.data(), + out_dist_sqr.data()); + const size_t usable = std::min(found, static_cast(k_neighbors)); + for (size_t k = 0; k < usable; ++k) { + knn_indices[base + k] = static_cast(ret_index[k]); + knn_distances[base + k] = std::sqrt(std::max(0.0f, out_dist_sqr[k])); } - if (out_dist_sqr.size() < static_cast(k_neighbors)) { - out_dist_sqr.resize(static_cast(k_neighbors)); + for (size_t k = usable; k < static_cast(k_neighbors); ++k) { + knn_indices[base + k] = row_idx; + knn_distances[base + k] = 0.0f; } + }, + 64); - const float query_pt[3] = {input.x[i], input.y[i], input.z[i]}; - const size_t num_results = tree.knnSearch( - query_pt, static_cast(k_neighbors), ret_index.data(), - out_dist_sqr.data()); - - float row_mass = 0.0f; - for (size_t k = 0; k < num_results; ++k) { - const float k_val = std::exp(-out_dist_sqr[k] / t_sq); - if (k_val < 1e-8f) { - out_dist_sqr[k] = 0.0f; - continue; - } + Eigen::VectorXf bandwidths_A; + compute_local_bandwidths(knn_distances, n, k_neighbors, knn_bandwidth, + bandwidths_A); + + std::vector kernel_entries_A(row_capacity, 0.0f); + for (int i = 0; i < n; ++i) { + const size_t base = static_cast(i) * row_stride; + const float bw_i = std::max(bandwidths_A[i], 1e-8f); + for (int k = 0; k < k_neighbors; ++k) { + const int j = knn_indices[base + static_cast(k)]; + const float bw_j = std::max(bandwidths_A[j], 1e-8f); + const float dist = knn_distances[base + static_cast(k)]; + kernel_entries_A[base + static_cast(k)] = + (dist * dist) / (bw_j * bw_i); + } + } - out_dist_sqr[k] = k_val; - row_mass += k_val; - } + const std::vector epsilons = build_epsilons(); + const auto [epsilon_A, dim_A] = + tune_kernel(kernel_entries_A, n, k_neighbors, epsilons); + + constexpr float kPi = 3.14159265358979323846f; + const float kernel_A_norm = + std::pow(kPi * std::max(epsilon_A, 1e-12f), dim_A * 0.5f); + Eigen::VectorXf density_estimate_A(n); + density_estimate_A.setZero(); + for (int i = 0; i < n; ++i) { + const size_t base = static_cast(i) * row_stride; + float row_sum = 0.0f; + for (int k = 0; k < k_neighbors; ++k) { + row_sum += std::exp(-kernel_entries_A[base + static_cast(k)] / + std::max(epsilon_A, 1e-12f)); + } + const float bw = std::max(bandwidths_A[i], 1e-8f); + density_estimate_A[i] = + (row_sum / std::max(kernel_A_norm, 1e-12f)) / + (static_cast(n) * std::pow(bw, dim_A)); + density_estimate_A[i] = std::max(density_estimate_A[i], 1e-12f); + } - if (row_mass <= 1e-12f) { - row_cols[base] = static_cast(i); - row_values[base] = 1.0f; - row_counts[i] = 1; - mu[row_idx] = 1.0f; - return; - } + Eigen::VectorXf bandwidths_B(n); + for (int i = 0; i < n; ++i) { + bandwidths_B[i] = + std::pow(density_estimate_A[i], input.bandwidth_variability); + } + std::vector bw_values(static_cast(n)); + for (int i = 0; i < n; ++i) { + bw_values[static_cast(i)] = bandwidths_B[i]; + } + const float bw_median = std::max(vector_median(std::move(bw_values)), 1e-8f); + bandwidths_B /= bw_median; + + std::vector kernel_entries_B(row_capacity, 0.0f); + for (int i = 0; i < n; ++i) { + const size_t base = static_cast(i) * row_stride; + const float bw_i = std::max(bandwidths_B[i], 1e-8f); + for (int k = 0; k < k_neighbors; ++k) { + const int j = knn_indices[base + static_cast(k)]; + const float bw_j = std::max(bandwidths_B[j], 1e-8f); + const float dist = knn_distances[base + static_cast(k)]; + kernel_entries_B[base + static_cast(k)] = + (dist * dist) / (bw_j * bw_i); + } + } - const float inv_row_mass = 1.0f / row_mass; - int count = 0; - for (size_t k = 0; k < num_results; ++k) { - const float k_val = out_dist_sqr[k]; - if (k_val <= 0.0f) { - continue; - } + const auto [epsilon_B, dim_B] = + tune_kernel(kernel_entries_B, n, k_neighbors, epsilons); + + Eigen::VectorXf density_estimate_B(n); + density_estimate_B.setZero(); + for (int i = 0; i < n; ++i) { + const size_t base = static_cast(i) * row_stride; + float row_sum = 0.0f; + for (int k = 0; k < k_neighbors; ++k) { + row_sum += std::exp(-kernel_entries_B[base + static_cast(k)] / + std::max(epsilon_B, 1e-12f)); + } + const float bw = std::max(bandwidths_B[i], 1e-8f); + density_estimate_B[i] = row_sum / std::pow(bw, dim_B); + density_estimate_B[i] = std::max(density_estimate_B[i], 1e-12f); + } - row_cols[base + static_cast(count)] = - static_cast(ret_index[k]); - row_values[base + static_cast(count)] = k_val * inv_row_mass; - ++count; - } + const float alpha = 1.0f - (input.c * 0.5f) + + input.bandwidth_variability * ((dim_B * 0.5f) + 1.0f); + Eigen::VectorXf density_alpha(n); + for (int i = 0; i < n; ++i) { + density_alpha[i] = std::pow(density_estimate_B[i], -alpha); + } - if (count == 0) { - row_cols[base] = static_cast(i); - row_values[base] = 1.0f; - row_counts[i] = 1; - mu[row_idx] = 1.0f; - return; + markov_row_offsets.resize(static_cast(n) + 1); + for (int i = 0; i <= n; ++i) { + markov_row_offsets[static_cast(i)] = i * k_neighbors; + } + markov_col_indices = knn_indices; + markov_values.assign(row_capacity, 0.0f); + + for (int i = 0; i < n; ++i) { + const size_t base = static_cast(i) * row_stride; + float row_sum = 0.0f; + for (int k = 0; k < k_neighbors; ++k) { + const int j = knn_indices[base + static_cast(k)]; + const float kernel_B = + std::exp(-kernel_entries_B[base + static_cast(k)] / + std::max(epsilon_B, 1e-12f)); + const float val = kernel_B * density_alpha[i] * density_alpha[j]; + markov_values[base + static_cast(k)] = val; + row_sum += val; + } + + if (row_sum <= 1e-30f) { + int self_idx = -1; + for (int k = 0; k < k_neighbors; ++k) { + if (knn_indices[base + static_cast(k)] == i) { + self_idx = k; + break; } + } + if (self_idx < 0) { + self_idx = 0; + } + for (int k = 0; k < k_neighbors; ++k) { + markov_values[base + static_cast(k)] = + (k == self_idx) ? 1.0f : 0.0f; + } + knn_kernel[base + static_cast(self_idx)] = 1.0f; + for (int k = 0; k < k_neighbors; ++k) { + if (k != self_idx) { + knn_kernel[base + static_cast(k)] = 0.0f; + } + } + continue; + } - row_counts[i] = count; - mu[row_idx] = row_mass; - }, - 64); + const float inv_row_sum = 1.0f / row_sum; + for (int k = 0; k < k_neighbors; ++k) { + const size_t idx = base + static_cast(k); + markov_values[idx] *= inv_row_sum; + knn_kernel[idx] = markov_values[idx]; + } + } - int nnz = 0; - for (size_t i = 0; i < n_verts; ++i) { - markov_row_offsets[i] = nnz; - nnz += row_counts[i]; + std::vector> markov_triplets; + markov_triplets.reserve(row_capacity); + for (int i = 0; i < n; ++i) { + const int begin = markov_row_offsets[static_cast(i)]; + const int end = markov_row_offsets[static_cast(i) + 1]; + for (int idx = begin; idx < end; ++idx) { + markov_triplets.emplace_back(i, markov_col_indices[static_cast(idx)], + markov_values[static_cast(idx)]); + } } - markov_row_offsets[n_verts] = nnz; - markov_col_indices.resize(static_cast(nnz)); - markov_values.resize(static_cast(nnz)); + Eigen::SparseMatrix P_row(n, n); + P_row.setFromTriplets(markov_triplets.begin(), markov_triplets.end()); + P_row.makeCompressed(); - core::parallel_for_index( - 0, static_cast(n_verts), - [&](int row_idx) { - const size_t i = static_cast(row_idx); - const int begin = markov_row_offsets[i]; - const int count = row_counts[i]; - const size_t src_base = i * row_stride; - for (int k = 0; k < count; ++k) { - markov_col_indices[static_cast(begin + k)] = - row_cols[src_base + static_cast(k)]; - markov_values[static_cast(begin + k)] = - row_values[src_base + static_cast(k)]; - } - }, - 64); + Eigen::SparseMatrix Pt_row = P_row.transpose(); + Pt_row.makeCompressed(); + + P_row += Pt_row; + P_row *= 0.5f; + P_row.makeCompressed(); - const float mu_sum = mu.sum(); + const Eigen::VectorXf ones = Eigen::VectorXf::Ones(n); + symmetric_row_sums = P_row * ones; + symmetric_row_sums = + symmetric_row_sums.array().max(1e-12f).matrix(); + + const float mu_sum = symmetric_row_sums.sum(); if (mu_sum > 1e-12f) { - mu /= mu_sum; + mu = symmetric_row_sums / mu_sum; } else { - mu = Eigen::VectorXf::Constant(static_cast(n_verts), - 1.0f / static_cast(n_verts)); + mu = Eigen::VectorXf::Constant(n, 1.0f / static_cast(n)); + } + + const int *outer = P_row.outerIndexPtr(); + const int *inner = P_row.innerIndexPtr(); + const float *values = P_row.valuePtr(); + const int nnz_sym = P_row.nonZeros(); + symmetric_row_offsets.assign(outer, outer + (n + 1)); + symmetric_col_indices.assign(inner, inner + nnz_sym); + symmetric_values.assign(values, values + nnz_sym); + + local_bandwidths.resize(n); + for (int i = 0; i < n; ++i) { + const float bw = std::max(bandwidths_B[i], 1e-8f); + local_bandwidths[i] = (std::max(epsilon_B, 1e-12f) * bw * bw) / 4.0f; + } + + immersion_coords.resize(n, 3); + for (int i = 0; i < n; ++i) { + float rx = 0.0f; + float ry = 0.0f; + float rz = 0.0f; + const int begin = markov_row_offsets[static_cast(i)]; + const int end = markov_row_offsets[static_cast(i) + 1]; + for (int idx = begin; idx < end; ++idx) { + const int j = markov_col_indices[static_cast(idx)]; + const float w = markov_values[static_cast(idx)]; + rx += w * input.x[static_cast(j)]; + ry += w * input.y[static_cast(j)]; + rz += w * input.z[static_cast(j)]; + } + immersion_coords(i, 0) = rx; + immersion_coords(i, 1) = ry; + immersion_coords(i, 2) = rz; } if (std::getenv("IGNEOUS_BENCH_MODE") == nullptr) { - std::cout << "[Diffusion] Built connectivity for " << n_verts - << " points.\n"; + std::cout << "[Diffusion] Built reference-style kernel for " << n_verts + << " points (k=" << k_neighbors << ", knn_bw=" << knn_bandwidth + << ").\n"; } } }; diff --git a/include/igneous/ops/geometry.hpp b/include/igneous/ops/geometry.hpp index 226256f..5d280e0 100644 --- a/include/igneous/ops/geometry.hpp +++ b/include/igneous/ops/geometry.hpp @@ -34,25 +34,54 @@ void fill_coordinate_vectors(const MeshT &mesh, coords[d].resize(static_cast(n_verts)); } + if constexpr (requires { mesh.topology.immersion_coords; }) { + const auto &immersion = mesh.topology.immersion_coords; + if (immersion.rows() == static_cast(n_verts) && immersion.cols() >= 3) { + for (size_t i = 0; i < n_verts; ++i) { + const int idx = static_cast(i); + coords[0][idx] = immersion(idx, 0); + coords[1][idx] = immersion(idx, 1); + coords[2][idx] = immersion(idx, 2); + } + return; + } + } + + for (size_t i = 0; i < n_verts; ++i) { + const auto p = mesh.geometry.get_vec3(i); + const int idx = static_cast(i); + coords[0][idx] = p.x; + coords[1][idx] = p.y; + coords[2][idx] = p.z; + } +} + +template +void fill_data_coordinate_vectors(const MeshT &mesh, + std::array &coords) { + const size_t n_verts = mesh.geometry.num_points(); + for (int d = 0; d < 3; ++d) { + coords[d].resize(static_cast(n_verts)); + } + for (size_t i = 0; i < n_verts; ++i) { const auto p = mesh.geometry.get_vec3(i); - coords[0][static_cast(i)] = p.x; - coords[1][static_cast(i)] = p.y; - coords[2][static_cast(i)] = p.z; + const int idx = static_cast(i); + coords[0][idx] = p.x; + coords[1][idx] = p.y; + coords[2][idx] = p.z; } } template void carre_du_champ(const MeshT &mesh, Eigen::Ref f, - Eigen::Ref h, float bandwidth, + Eigen::Ref h, [[maybe_unused]] float bandwidth, Eigen::Ref gamma_out) { [[maybe_unused]] const int expected_size = static_cast(mesh.geometry.num_points()); assert(gamma_out.size() == expected_size); gamma_out.setZero(); - const float inv_2t = 1.0f / std::max(2.0f * bandwidth, 1e-8f); - if constexpr (requires { mesh.topology.markov_row_offsets; mesh.topology.markov_col_indices; @@ -64,25 +93,33 @@ void carre_du_champ(const MeshT &mesh, Eigen::Ref f, assert(row_offsets.size() == static_cast(expected_size) + 1); - const bool use_gpu = - core::compute_backend_from_env() == core::ComputeBackend::Gpu && - (core::gpu::gpu_force_enabled() || expected_size >= core::gpu::gpu_min_rows()); - if (use_gpu) { - if (core::gpu::carre_du_champ( - static_cast(&mesh.topology), - std::span(row_offsets.data(), row_offsets.size()), - std::span(col_indices.data(), col_indices.size()), - std::span(weights.data(), weights.size()), - std::span(f.data(), static_cast(expected_size)), - std::span(h.data(), static_cast(expected_size)), - inv_2t, - std::span(gamma_out.data(), static_cast(expected_size)))) { - return; - } - } - const float *f_data = f.data(); const float *h_data = h.data(); + Eigen::VectorXf means_f = Eigen::VectorXf::Zero(expected_size); + Eigen::VectorXf means_h = Eigen::VectorXf::Zero(expected_size); + const bool use_mean_centres = [&]() { + if constexpr (requires { mesh.topology.use_mean_centres; }) { + return mesh.topology.use_mean_centres; + } + return false; + }(); + + if (use_mean_centres) { + for (int i = 0; i < expected_size; ++i) { + const int begin = row_offsets[static_cast(i)]; + const int end = row_offsets[static_cast(i) + 1]; + const int *cols = col_indices.data() + begin; + const float *w = weights.data() + begin; + float mean_f = 0.0f; + float mean_h = 0.0f; + for (int idx = 0; idx < (end - begin); ++idx) { + mean_f += w[idx] * f_data[cols[idx]]; + mean_h += w[idx] * h_data[cols[idx]]; + } + means_f[i] = mean_f; + means_h[i] = mean_h; + } + } for (int i = 0; i < expected_size; ++i) { const int begin = row_offsets[static_cast(i)]; @@ -90,8 +127,8 @@ void carre_du_champ(const MeshT &mesh, Eigen::Ref f, const int count = end - begin; const int *cols = col_indices.data() + begin; const float *w = weights.data() + begin; - const float fi = f_data[i]; - const float hi = h_data[i]; + const float center_f = use_mean_centres ? means_f[i] : f_data[i]; + const float center_h = use_mean_centres ? means_h[i] : h_data[i]; float acc = 0.0f; int idx = 0; @@ -100,17 +137,23 @@ void carre_du_champ(const MeshT &mesh, Eigen::Ref f, const int j1 = cols[idx + 1]; const int j2 = cols[idx + 2]; const int j3 = cols[idx + 3]; - acc += w[idx + 0] * (f_data[j0] - fi) * (h_data[j0] - hi); - acc += w[idx + 1] * (f_data[j1] - fi) * (h_data[j1] - hi); - acc += w[idx + 2] * (f_data[j2] - fi) * (h_data[j2] - hi); - acc += w[idx + 3] * (f_data[j3] - fi) * (h_data[j3] - hi); + acc += w[idx + 0] * (f_data[j0] - center_f) * (h_data[j0] - center_h); + acc += w[idx + 1] * (f_data[j1] - center_f) * (h_data[j1] - center_h); + acc += w[idx + 2] * (f_data[j2] - center_f) * (h_data[j2] - center_h); + acc += w[idx + 3] * (f_data[j3] - center_f) * (h_data[j3] - center_h); } for (; idx < count; ++idx) { const int j = cols[idx]; - acc += w[idx] * (f_data[j] - fi) * (h_data[j] - hi); + acc += w[idx] * (f_data[j] - center_f) * (h_data[j] - center_h); } - gamma_out[i] = acc; + float denom = 2.0f; + if constexpr (requires { mesh.topology.local_bandwidths; }) { + if (mesh.topology.local_bandwidths.size() == expected_size) { + denom = 2.0f * std::max(mesh.topology.local_bandwidths[i], 1e-8f); + } + } + gamma_out[i] = acc / denom; } } else { const auto &P = mesh.topology.P; @@ -123,9 +166,8 @@ void carre_du_champ(const MeshT &mesh, Eigen::Ref f, gamma_out[i] += w * (f[j] - f[i]) * (h[j] - h[i]); } } + gamma_out *= 0.5f; } - - gamma_out *= inv_2t; } template diff --git a/include/igneous/ops/hodge.hpp b/include/igneous/ops/hodge.hpp index daf4a8f..b125b2f 100644 --- a/include/igneous/ops/hodge.hpp +++ b/include/igneous/ops/hodge.hpp @@ -2,10 +2,12 @@ #include #include +#include #include #include #include #include +#include #include #include @@ -157,18 +159,55 @@ inline Eigen::MatrixXf compute_hodge_laplacian_matrix( return L_down + E_up; } -inline auto compute_hodge_spectrum(const Eigen::MatrixXf &laplacian, - const Eigen::MatrixXf &mass_matrix) { - Eigen::GeneralizedSelfAdjointEigenSolver solver(laplacian, - mass_matrix); - return std::make_pair(solver.eigenvalues(), solver.eigenvectors()); +inline std::pair +compute_hodge_spectrum(const Eigen::MatrixXf &laplacian, + const Eigen::MatrixXf &mass_matrix, + float rcond = 1e-5f) { + Eigen::SelfAdjointEigenSolver mass_solver(mass_matrix); + if (mass_solver.info() != Eigen::Success) { + return std::make_pair(Eigen::VectorXf(), Eigen::MatrixXf()); + } + + const auto &mass_evals = mass_solver.eigenvalues(); + const auto &mass_evecs = mass_solver.eigenvectors(); + std::vector keep_indices; + keep_indices.reserve(static_cast(mass_evals.size())); + for (int i = 0; i < mass_evals.size(); ++i) { + if (mass_evals[i] > rcond) { + keep_indices.push_back(i); + } + } + + if (keep_indices.empty()) { + return std::make_pair(Eigen::VectorXf(), Eigen::MatrixXf::Zero(laplacian.rows(), 0)); + } + + const int k = static_cast(keep_indices.size()); + Eigen::MatrixXf phi(laplacian.rows(), k); + for (int col = 0; col < k; ++col) { + const int idx = keep_indices[static_cast(col)]; + const float inv_sqrt = 1.0f / std::sqrt(std::max(mass_evals[idx], 1e-12f)); + phi.col(col) = mass_evecs.col(idx) * inv_sqrt; + } + + const Eigen::MatrixXf restricted = phi.transpose() * laplacian * phi; + Eigen::SelfAdjointEigenSolver solver(restricted); + if (solver.info() != Eigen::Success) { + return std::make_pair(Eigen::VectorXf(), Eigen::MatrixXf::Zero(laplacian.rows(), 0)); + } + + const Eigen::VectorXf evals = solver.eigenvalues(); + const Eigen::MatrixXf evecs = phi * solver.eigenvectors(); + return std::make_pair(evals, evecs); } template Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, const Eigen::VectorXf &alpha_coeffs, float bandwidth, - float epsilon = 1e-3f) { + float lambda = 1.0f, + int positive_imag_mode = 0, + std::complex *selected_eval = nullptr) { const auto &U = mesh.topology.eigen_basis; const auto &mu = mesh.topology.mu; @@ -214,32 +253,87 @@ Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, }, 8); - Eigen::MatrixXf L_eps = X_op - epsilon * Eigen::MatrixXf::Identity(n0, n0); + Eigen::MatrixXf laplacian0_weak = Eigen::MatrixXf::Zero(n0, n0); + Eigen::VectorXf gamma_local(n_verts_i); + for (int i = 0; i < n0; ++i) { + for (int j = i; j < n0; ++j) { + carre_du_champ(mesh, U.col(i), U.col(j), bandwidth, gamma_local); + const float val = (gamma_local.array() * mu.array()).sum(); + laplacian0_weak(i, j) = val; + laplacian0_weak(j, i) = val; + } + } + + Eigen::MatrixXf function_gram = U.transpose() * (U.array().colwise() * mu.array()).matrix(); + Eigen::MatrixXf operator_weak = X_op - (lambda * laplacian0_weak); - Eigen::EigenSolver solver(L_eps); - const auto evals = solver.eigenvalues(); - const auto evecs = solver.eigenvectors(); + Eigen::GeneralizedEigenSolver solver(operator_weak, function_gram); + if (solver.info() != Eigen::Success) { + return Eigen::VectorXf::Zero(static_cast(n_verts)); + } + + Eigen::VectorXcf evals = solver.eigenvalues().cast>(); + Eigen::MatrixXcf evecs = solver.eigenvectors().cast>(); - int best_idx = 0; - float min_mag = std::numeric_limits::max(); + std::vector order(static_cast(n0)); for (int i = 0; i < n0; ++i) { - const float mag = std::abs(evals(i)); - if (mag > 1e-6f && mag < min_mag) { - min_mag = mag; - best_idx = i; + order[static_cast(i)] = i; + } + std::sort(order.begin(), order.end(), [&](int lhs, int rhs) { + const std::complex a = evals[lhs]; + const std::complex b = evals[rhs]; + const float abs_a = std::abs(a); + const float abs_b = std::abs(b); + if (std::abs(abs_a - abs_b) > 1e-7f) { + return abs_a < abs_b; + } + if (std::abs(a.real() - b.real()) > 1e-7f) { + return a.real() < b.real(); } + return a.imag() < b.imag(); + }); + + Eigen::VectorXcf sorted_evals(n0); + Eigen::MatrixXcf sorted_evecs(n0, n0); + for (int i = 0; i < n0; ++i) { + const int src = order[static_cast(i)]; + sorted_evals[i] = evals[src]; + sorted_evecs.col(i) = evecs.col(src); } - const Eigen::VectorXcf z_coeffs = evecs.col(best_idx); + std::vector positive_imag_indices; + for (int i = 0; i < n0; ++i) { + if (sorted_evals[i].imag() > 1e-6f) { + positive_imag_indices.push_back(i); + } + } + + if (positive_imag_indices.empty()) { + return Eigen::VectorXf::Zero(static_cast(n_verts)); + } + + const int mode = + std::clamp(positive_imag_mode, 0, + static_cast(positive_imag_indices.size()) - 1); + const int best_idx = positive_imag_indices[static_cast(mode)]; + if (selected_eval != nullptr) { + *selected_eval = sorted_evals[best_idx]; + } + + const Eigen::VectorXcf z_coeffs = sorted_evecs.col(best_idx); Eigen::VectorXf theta(static_cast(n_verts)); - constexpr float kPi = 3.14159265358979323846f; + constexpr float kTwoPi = 6.28318530717958647692f; for (size_t i = 0; i < n_verts; ++i) { std::complex zi(0.0f, 0.0f); for (int k = 0; k < n0; ++k) { zi += U(static_cast(i), k) * z_coeffs(k); } - theta[static_cast(i)] = (std::arg(zi) + kPi) / (2.0f * kPi); + float angle = std::atan2(zi.imag(), zi.real()); + if (angle < 0.0f) { + angle += static_cast(kTwoPi); + } + theta[static_cast(i)] = angle; } return theta; diff --git a/include/igneous/ops/spectral.hpp b/include/igneous/ops/spectral.hpp index 7e84243..0fd87d7 100644 --- a/include/igneous/ops/spectral.hpp +++ b/include/igneous/ops/spectral.hpp @@ -1,8 +1,8 @@ #pragma once #include #include -#include #include +#include #include #include #include @@ -116,6 +116,52 @@ class MarkovSymmetricCsrMatProd { int n_ = 0; }; +class SymmetricKernelCsrMatProd { +public: + using Scalar = float; + + SymmetricKernelCsrMatProd(const std::vector &row_offsets, + const std::vector &col_indices, + const std::vector &values, int n) + : row_offsets_(row_offsets), col_indices_(col_indices), values_(values), + n_(n) {} + + [[nodiscard]] int rows() const { return n_; } + [[nodiscard]] int cols() const { return n_; } + + void perform_op(const Scalar *x_in, Scalar *y_out) const { + core::parallel_for_index( + 0, n_, + [&](int i) { + const int begin = row_offsets_[static_cast(i)]; + const int end = row_offsets_[static_cast(i) + 1]; + const int count = end - begin; + const int *cols = col_indices_.data() + begin; + const float *vals = values_.data() + begin; + + float acc = 0.0f; + int k = 0; + for (; k + 3 < count; k += 4) { + acc += vals[k + 0] * x_in[cols[k + 0]]; + acc += vals[k + 1] * x_in[cols[k + 1]]; + acc += vals[k + 2] * x_in[cols[k + 2]]; + acc += vals[k + 3] * x_in[cols[k + 3]]; + } + for (; k < count; ++k) { + acc += vals[k] * x_in[cols[k]]; + } + y_out[i] = acc; + }, + 8192); + } + +private: + const std::vector &row_offsets_; + const std::vector &col_indices_; + const std::vector &values_; + int n_ = 0; +}; + // Computes the first k eigenvectors of the Markov Chain P. template void compute_eigenbasis(MeshT &mesh, int n_eigenvectors) { @@ -187,62 +233,76 @@ void compute_eigenbasis(MeshT &mesh, int n_eigenvectors) { } }; + if (n <= 1) { + mesh.topology.eigen_basis = + Eigen::MatrixXf::Ones(std::max(1, n), std::max(1, n)); + return; + } + if constexpr (requires { mesh.topology.markov_row_offsets; mesh.topology.markov_col_indices; mesh.topology.markov_values; }) { - const bool use_symmetric_solver = n >= 2200; - if (use_symmetric_solver) { - Eigen::VectorXf sqrt_mu(n); - Eigen::VectorXf inv_sqrt_mu(n); - for (int i = 0; i < n; ++i) { - const float mu_i = std::max(mesh.topology.mu[i], 1e-12f); - const float sqrt_mu_i = std::sqrt(mu_i); - sqrt_mu[i] = sqrt_mu_i; - inv_sqrt_mu[i] = 1.0f / sqrt_mu_i; - } + if constexpr (requires { + mesh.topology.symmetric_row_offsets; + mesh.topology.symmetric_col_indices; + mesh.topology.symmetric_values; + mesh.topology.symmetric_row_sums; + }) { + const auto &sym_row_offsets = mesh.topology.symmetric_row_offsets; + const auto &sym_col_indices = mesh.topology.symmetric_col_indices; + const auto &sym_values = mesh.topology.symmetric_values; + const auto &row_sums = mesh.topology.symmetric_row_sums; + if (!sym_row_offsets.empty() && + static_cast(sym_row_offsets.size()) == n + 1 && + row_sums.size() == n) { + const int k_eval = std::max(1, std::min(n_eigenvectors, std::max(1, n - 1))); + const int full_ncv = std::min(n, std::max(2 * k_eval + 1, 20)); + const bool try_compact = k_eval >= 32; + const int compact_ncv = + try_compact ? std::min(n, std::max(k_eval + 16, 20)) : full_ncv; + + SymmetricKernelCsrMatProd op(sym_row_offsets, sym_col_indices, sym_values, + n); + const auto solve_symmetric = [&](int ncv) -> bool { + Spectra::SymEigsSolver eigs(op, k_eval, ncv); + eigs.init(); + const int nconv = eigs.compute(Spectra::SortRule::LargestMagn); + if (eigs.info() != Spectra::CompInfo::Successful || nconv <= 0) { + return false; + } - MarkovSymmetricCsrMatProd sym_op(mesh.topology.markov_row_offsets, - mesh.topology.markov_col_indices, - mesh.topology.markov_values, sqrt_mu, - inv_sqrt_mu, n); - - const int full_ncv = std::min(n, std::max(2 * n_eigenvectors + 1, 20)); - const bool try_compact = n_eigenvectors >= 32; - const int compact_ncv = - try_compact ? std::min(n, std::max(n_eigenvectors + 16, 20)) : full_ncv; - - const auto solve_symmetric = [&](int ncv) -> bool { - Spectra::SymEigsSolver eigs( - sym_op, n_eigenvectors, ncv); - eigs.init(); - const int nconv = eigs.compute(Spectra::SortRule::LargestAlge); - if (eigs.info() != Spectra::CompInfo::Successful || nconv <= 0) { - return false; - } + Eigen::MatrixXf basis = eigs.eigenvectors(nconv); + Eigen::VectorXf inv_sqrt_rows(n); + for (int i = 0; i < n; ++i) { + inv_sqrt_rows[i] = 1.0f / std::sqrt(std::max(row_sums[i], 1e-12f)); + } + basis = basis.array().colwise() * inv_sqrt_rows.array(); + basis = basis.rowwise().reverse().eval(); + if (std::abs(basis(0, 0)) > 1e-12f) { + basis /= basis(0, 0); + } + mesh.topology.eigen_basis = basis; - mesh.topology.eigen_basis = eigs.eigenvectors(nconv); - mesh.topology.eigen_basis = - mesh.topology.eigen_basis.array().colwise() * inv_sqrt_mu.array(); + if (verbose) { + std::cout << "[Spectral] Converged! Found " << nconv + << " eigenvectors. Basis shape: " + << mesh.topology.eigen_basis.rows() << "x" + << mesh.topology.eigen_basis.cols() << "\n"; + } + return true; + }; - if (verbose) { - std::cout << "[Spectral] Converged! Found " << nconv - << " eigenvectors. Basis shape: " - << mesh.topology.eigen_basis.rows() << "x" - << mesh.topology.eigen_basis.cols() << "\n"; + if ((try_compact && solve_symmetric(compact_ncv)) || + solve_symmetric(full_ncv)) { + return; } - return true; - }; - - if ((try_compact && solve_symmetric(compact_ncv)) || - solve_symmetric(full_ncv)) { - return; - } - if (verbose) { - std::cerr - << "[Spectral] Symmetric solve failed, falling back to generic solver.\n"; + if (verbose) { + std::cerr << "[Spectral] Symmetric-kernel solve failed, falling back to " + "generic solver.\n"; + } } } diff --git a/notes/hodge/journal.md b/notes/hodge/journal.md new file mode 100644 index 0000000..1802af7 --- /dev/null +++ b/notes/hodge/journal.md @@ -0,0 +1,100 @@ +# Hodge Parity Journal + +Use one entry per parity hypothesis. + +## Entry Template +- Timestamp: +- Commit: +- Hypothesis: +- Files touched: +- Commands: +- Numeric parity deltas: +- Decision: `kept` | `rejected` | `deferred` +- Notes: + +## 2026-02-14 Baseline Harness + Gauge-Invariant Metrics +- Timestamp: 2026-02-14T16:09:21Z +- Commit: working tree (uncommitted) +- Hypothesis: Establish deterministic cross-language parity artifacts first, then compare harmonic forms in ambient space (not raw coefficient coordinates) to avoid basis-order artifacts. +- Files touched: + - `include/igneous/data/topology.hpp` + - `include/igneous/ops/geometry.hpp` + - `include/igneous/ops/spectral.hpp` + - `include/igneous/ops/hodge.hpp` + - `src/main_hodge.cpp` + - `src/main_diffusion.cpp` + - `src/main_spectral.cpp` + - `tests/test_topology_diffusion.cpp` + - `tests/test_ops_spectral_geometry.cpp` + - `tests/test_ops_hodge.cpp` + - `benches/bench_dod.cpp` + - `benches/bench_pipelines.cpp` + - `scripts/hodge/setup_reference_env.sh` + - `scripts/hodge/generate_input_torus.py` + - `scripts/hodge/run_reference_hodge.py` + - `scripts/hodge/run_cpp_hodge.sh` + - `scripts/hodge/compare_hodge_outputs.py` + - `scripts/hodge/run_parity_round.sh` + - `notes/hodge/journal.md` +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `scripts/hodge/run_parity_round.sh` + - `python3 scripts/hodge/compare_hodge_outputs.py --reference-dir notes/hodge/results/round_20260214-160921/reference --cpp-dir notes/hodge/results/round_20260214-160921/cpp --output-markdown notes/hodge/results/round_20260214-160921/report/parity_report.md --output-json notes/hodge/results/round_20260214-160921/report/parity_report.json --label \"Hodge parity round 20260214-160921\"` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-160921/report/parity_report.json` + - Composite score: `0.310740` + - Harmonic subspace max principal angle: `9.229825 deg` + - Harmonic Procrustes error: `0.129113` + - Circular correlation min: `0.098627` + - Circular wrapped P95 max: `2.936349 rad` +- Decision: `kept` +- Notes: + - Concrete structural mismatch identified: coefficient-space harmonic comparison is not gauge-invariant due basis ordering differences. + - Primary harmonic parity metrics now use `harmonic_ambient.csv` / `reference_harmonic_ambient.csv`. + +## 2026-02-14 Hypothesis A1: Normalized Symmetric-Kernel Eigensolve +- Timestamp: 2026-02-14T16:11:23Z +- Commit: working tree (uncommitted) +- Hypothesis: Solving eigenvectors on `D^{-1/2} K D^{-1/2}` instead of `K` will improve harmonic-form parity. +- Files touched: + - `include/igneous/ops/spectral.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-160921/report/parity_report.json scripts/hodge/run_parity_round.sh` + - `python3 scripts/hodge/compare_hodge_outputs.py --reference-dir notes/hodge/results/round_20260214-161123/reference --cpp-dir notes/hodge/results/round_20260214-161123/cpp --output-markdown notes/hodge/results/round_20260214-161123/report/parity_report.md --output-json notes/hodge/results/round_20260214-161123/report/parity_report.json --label \"Hodge parity round 20260214-161123\" --previous-json notes/hodge/results/round_20260214-160921/report/parity_report.json` +- Numeric parity deltas: + - Report: `notes/hodge/results/round_20260214-161123/report/parity_report.json` + - Composite: `0.249454` (`-19.72%` vs baseline) + - Harmonic subspace max principal angle: `4.066515 deg` (improved) + - Harmonic Procrustes error: `0.059640` (improved) + - Circular correlation min: `0.050558` (regressed) + - Circular wrapped P95 max: `2.964019 rad` (regressed) +- Decision: `rejected` +- Notes: + - Improvement was below the required `20%` keep threshold. + - Change discarded per commit gate policy. + +## 2026-02-14 Hypothesis D1: Circular Mode Selection Parity +- Timestamp: 2026-02-14T16:18:50Z +- Commit: working tree (uncommitted) +- Hypothesis: The reference workflow uses the first positive-imaginary circular mode for each harmonic form; setting second mode default to `0` (not `1`) should restore circular parity. +- Files touched: + - `src/main_hodge.cpp` + - `scripts/hodge/run_reference_hodge.py` + - `scripts/hodge/run_cpp_hodge.sh` + - `scripts/hodge/run_parity_round.sh` +- Commands: + - `cmake --build build -j8` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-160921/report/parity_report.json scripts/hodge/run_parity_round.sh` +- Numeric parity deltas: + - Report: `notes/hodge/results/round_20260214-161850/report/parity_report.json` + - Composite: `0.045839` (`-85.25%` vs baseline) + - Harmonic subspace max principal angle: `4.155913 deg` + - Harmonic Procrustes error: `0.060369` + - Circular correlation min: `0.995021` (major improvement) + - Circular wrapped P95 max: `0.170368 rad` (major improvement) +- Decision: `kept` +- Notes: + - This is a structural parity fix in eigenmode selection policy, not a numerical micro-optimization. + - Remaining gaps are now concentrated in harmonic-form thresholds and circular P95 tail. diff --git a/notes/hodge/results/.gitkeep b/notes/hodge/results/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/scripts/hodge/compare_hodge_outputs.py b/scripts/hodge/compare_hodge_outputs.py new file mode 100755 index 0000000..d2ac133 --- /dev/null +++ b/scripts/hodge/compare_hodge_outputs.py @@ -0,0 +1,456 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import csv +import json +import math +from pathlib import Path +from typing import Any + +import numpy as np + + +def load_csv_columns(path: Path) -> dict[str, np.ndarray]: + with path.open(newline="") as handle: + reader = csv.DictReader(handle) + if reader.fieldnames is None: + raise ValueError(f"CSV has no header: {path}") + + rows: dict[str, list[float]] = {field: [] for field in reader.fieldnames} + for row in reader: + for field in reader.fieldnames: + value = row.get(field, "") + rows[field].append(float(value) if value not in (None, "") else math.nan) + + return {k: np.asarray(v, dtype=np.float64) for k, v in rows.items()} + + +def sorted_prefixed_columns(columns: dict[str, np.ndarray], prefix: str) -> list[str]: + names = [name for name in columns.keys() if name.startswith(prefix)] + + def sort_key(name: str) -> tuple[int, str]: + suffix = name[len(prefix) :] + suffix = suffix.lstrip("_") + try: + return (int(suffix), name) + except ValueError: + return (10**9, name) + + return sorted(names, key=sort_key) + + +def to_xyz(columns: dict[str, np.ndarray]) -> np.ndarray: + required = ("x", "y", "z") + for field in required: + if field not in columns: + raise ValueError(f"Missing column '{field}'") + return np.column_stack([columns["x"], columns["y"], columns["z"]]) + + +def sorted_ambient_form_indices(columns: dict[str, np.ndarray]) -> list[int]: + indices: list[int] = [] + for name in columns.keys(): + if not name.startswith("form") or not name.endswith("_x"): + continue + middle = name[len("form") : -len("_x")] + try: + idx = int(middle) + except ValueError: + continue + if f"form{idx}_y" in columns and f"form{idx}_z" in columns: + indices.append(idx) + return sorted(set(indices)) + + +def harmonic_ambient_matrix( + columns: dict[str, np.ndarray], form_indices: list[int], n_points: int +) -> np.ndarray: + vectors: list[np.ndarray] = [] + for form_idx in form_indices: + vec = np.column_stack( + [ + columns[f"form{form_idx}_x"][:n_points], + columns[f"form{form_idx}_y"][:n_points], + columns[f"form{form_idx}_z"][:n_points], + ] + ).reshape(-1) + vectors.append(vec) + return np.column_stack(vectors) if vectors else np.zeros((n_points * 3, 0)) + + +def orthonormal_basis(matrix: np.ndarray, rcond: float) -> tuple[np.ndarray, int]: + if matrix.size == 0: + return np.zeros((matrix.shape[0], 0), dtype=np.float64), 0 + + u, s, _ = np.linalg.svd(matrix, full_matrices=False) + if s.size == 0: + return np.zeros((matrix.shape[0], 0), dtype=np.float64), 0 + + tol = max(rcond * float(s[0]), rcond) + rank = int(np.sum(s > tol)) + return u[:, :rank], rank + + +def principal_angles_deg(a: np.ndarray, b: np.ndarray, rcond: float) -> np.ndarray: + qa, rank_a = orthonormal_basis(a, rcond) + qb, rank_b = orthonormal_basis(b, rcond) + rank = min(rank_a, rank_b) + if rank == 0: + return np.asarray([90.0], dtype=np.float64) + + cross = qa[:, :rank].T @ qb[:, :rank] + singular_vals = np.linalg.svd(cross, compute_uv=False) + singular_vals = np.clip(singular_vals, -1.0, 1.0) + return np.degrees(np.arccos(singular_vals)) + + +def procrustes_relative_error(reference: np.ndarray, candidate: np.ndarray) -> float: + if reference.size == 0 or candidate.size == 0: + return float("nan") + + c = candidate.T @ reference + u, _, vt = np.linalg.svd(c, full_matrices=False) + r = u @ vt + aligned = candidate @ r + denom = max(float(np.linalg.norm(reference)), 1e-12) + return float(np.linalg.norm(aligned - reference) / denom) + + +def circular_mode_metrics(theta_ref: np.ndarray, theta_cpp: np.ndarray) -> dict[str, float | str]: + z_ref = np.exp(1.0j * theta_ref) + + candidates: list[tuple[str, np.ndarray]] = [ + ("direct", np.exp(1.0j * theta_cpp)), + ("conjugate", np.exp(-1.0j * theta_cpp)), + ] + + best: dict[str, float | str] | None = None + for orientation, z_raw in candidates: + corr_raw = np.vdot(z_raw, z_ref) + phase = float(np.angle(corr_raw)) if np.abs(corr_raw) > 0 else 0.0 + z_aligned = z_raw * np.exp(1.0j * phase) + + denom = max(float(np.linalg.norm(z_ref) * np.linalg.norm(z_aligned)), 1e-12) + corr_mag = float(np.abs(np.vdot(z_ref, z_aligned)) / denom) + + angle_diff = np.angle(z_ref * np.conj(z_aligned)) + abs_diff = np.abs(angle_diff) + + entry: dict[str, float | str] = { + "orientation": orientation, + "phase_shift_rad": phase, + "complex_correlation": corr_mag, + "wrapped_error_mean_rad": float(np.mean(abs_diff)), + "wrapped_error_p95_rad": float(np.quantile(abs_diff, 0.95)), + "wrapped_error_max_rad": float(np.max(abs_diff)), + } + + if best is None: + best = entry + continue + + if entry["complex_correlation"] > best["complex_correlation"]: + best = entry + elif entry["complex_correlation"] == best["complex_correlation"] and entry[ + "wrapped_error_p95_rad" + ] < best["wrapped_error_p95_rad"]: + best = entry + + assert best is not None + return best + + +def finite_or_default(value: float, default: float) -> float: + return value if np.isfinite(value) else default + + +def composite_score( + harmonic_subspace_max_angle_deg: float, + harmonic_procrustes_rel_error: float, + circular_complex_correlation_mean: float, + circular_wrapped_p95_rad_max: float, +) -> float: + # Lower is better. + angle_component = finite_or_default(harmonic_subspace_max_angle_deg / 90.0, 1.0) + proc_component = finite_or_default(harmonic_procrustes_rel_error, 1.0) + corr_component = finite_or_default(1.0 - circular_complex_correlation_mean, 1.0) + circ_component = finite_or_default(circular_wrapped_p95_rad_max / math.pi, 1.0) + + return float( + 0.35 * angle_component + + 0.35 * proc_component + + 0.15 * corr_component + + 0.15 * circ_component + ) + + +def safe_mean(values: list[float]) -> float: + if not values: + return float("nan") + return float(np.mean(np.asarray(values, dtype=np.float64))) + + +def render_markdown(payload: dict[str, Any]) -> str: + summary = payload["summary"] + circular = payload["circular_modes"] + gates = payload["gates"] + + lines: list[str] = [] + lines.append(f"# {payload['label']}") + lines.append("") + lines.append("## Inputs") + lines.append("") + lines.append(f"- Reference dir: `{payload['reference_dir']}`") + lines.append(f"- C++ dir: `{payload['cpp_dir']}`") + lines.append("") + lines.append("## Primary Metrics") + lines.append("") + lines.append(f"- Harmonic subspace max principal angle (deg): `{summary['harmonic_subspace_max_angle_deg']:.6f}`") + lines.append(f"- Harmonic Procrustes relative error: `{summary['harmonic_procrustes_rel_error']:.6f}`") + lines.append(f"- Circular complex correlation (min): `{summary['circular_complex_correlation_min']:.6f}`") + lines.append(f"- Circular complex correlation (mean): `{summary['circular_complex_correlation_mean']:.6f}`") + lines.append(f"- Circular wrapped angular P95 error max (rad): `{summary['circular_wrapped_p95_rad_max']:.6f}`") + lines.append(f"- Composite parity score (lower is better): `{summary['composite_score']:.6f}`") + lines.append("") + lines.append("## Circular Modes") + lines.append("") + lines.append("| Mode | Orientation | Correlation | P95 (rad) | Mean (rad) | Max (rad) |") + lines.append("| --- | --- | ---: | ---: | ---: | ---: |") + for row in circular: + lines.append( + f"| `{row['name']}` | `{row['orientation']}` | {row['complex_correlation']:.6f} | {row['wrapped_error_p95_rad']:.6f} | {row['wrapped_error_mean_rad']:.6f} | {row['wrapped_error_max_rad']:.6f} |" + ) + + lines.append("") + lines.append("## Gates") + lines.append("") + lines.append(f"- Final pass gate: `{'pass' if gates['final_pass'] else 'fail'}`") + lines.append(f"- Commit improvement gate: `{'pass' if gates['commit_gate_pass'] else 'fail'}`") + lines.append(f"- Commit gate details: `{gates['commit_gate_details']}`") + if gates.get("improvement_vs_previous") is not None: + lines.append( + f"- Composite improvement vs previous: `{100.0 * gates['improvement_vs_previous']:.2f}%`" + ) + + if payload["warnings"]: + lines.append("") + lines.append("## Warnings") + lines.append("") + for warning in payload["warnings"]: + lines.append(f"- {warning}") + + lines.append("") + return "\n".join(lines) + + +def main() -> int: + parser = argparse.ArgumentParser(description="Compare reference and C++ Hodge outputs") + parser.add_argument("--reference-dir", required=True) + parser.add_argument("--cpp-dir", required=True) + parser.add_argument("--output-markdown", required=True) + parser.add_argument("--output-json", required=True) + parser.add_argument("--label", default="Hodge parity comparison") + parser.add_argument("--previous-json", default="") + parser.add_argument("--gate-improvement", type=float, default=0.20) + parser.add_argument("--rcond", type=float, default=1e-10) + parser.add_argument("--final-max-angle-deg", type=float, default=3.0) + parser.add_argument("--final-max-procrustes", type=float, default=0.05) + parser.add_argument("--final-min-correlation", type=float, default=0.99) + parser.add_argument("--final-max-p95-rad", type=float, default=0.10) + args = parser.parse_args() + + reference_dir = Path(args.reference_dir) + cpp_dir = Path(args.cpp_dir) + output_markdown = Path(args.output_markdown) + output_json = Path(args.output_json) + + warnings: list[str] = [] + + ref_points_cols = load_csv_columns(reference_dir / "reference_points.csv") + cpp_points_cols = load_csv_columns(cpp_dir / "points.csv") + ref_points = to_xyz(ref_points_cols) + cpp_points = to_xyz(cpp_points_cols) + + n_points = min(ref_points.shape[0], cpp_points.shape[0]) + if ref_points.shape[0] != cpp_points.shape[0]: + warnings.append( + f"Point count mismatch: reference={ref_points.shape[0]}, cpp={cpp_points.shape[0]}; cropped to {n_points}." + ) + ref_points = ref_points[:n_points] + cpp_points = cpp_points[:n_points] + + point_diff = ref_points - cpp_points + point_rmse = float(np.sqrt(np.mean(point_diff**2))) + point_max_abs = float(np.max(np.abs(point_diff))) + if point_rmse > 1e-6: + warnings.append(f"Input point clouds differ (rmse={point_rmse:.3e}).") + + ref_ambient_cols = load_csv_columns(reference_dir / "reference_harmonic_ambient.csv") + cpp_ambient_cols = load_csv_columns(cpp_dir / "harmonic_ambient.csv") + ref_form_indices = sorted_ambient_form_indices(ref_ambient_cols) + cpp_form_indices = sorted_ambient_form_indices(cpp_ambient_cols) + common_form_indices = sorted(set(ref_form_indices).intersection(cpp_form_indices)) + if len(ref_form_indices) != len(cpp_form_indices): + warnings.append( + f"Harmonic ambient form count mismatch: reference={len(ref_form_indices)}, cpp={len(cpp_form_indices)}." + ) + if not common_form_indices: + raise ValueError("No common harmonic ambient form columns found.") + + ref_harmonic_ambient = harmonic_ambient_matrix( + ref_ambient_cols, common_form_indices, n_points + ) + cpp_harmonic_ambient = harmonic_ambient_matrix( + cpp_ambient_cols, common_form_indices, n_points + ) + + angles_deg = principal_angles_deg(ref_harmonic_ambient, cpp_harmonic_ambient, args.rcond) + max_angle_deg = float(np.max(angles_deg)) + procrustes_err = procrustes_relative_error(ref_harmonic_ambient, cpp_harmonic_ambient) + + # Secondary diagnostics in coefficient space (basis-order dependent). + ref_coeff_cols = load_csv_columns(reference_dir / "reference_harmonic_coeffs.csv") + cpp_coeff_cols = load_csv_columns(cpp_dir / "harmonic_coeffs.csv") + ref_form_cols = sorted_prefixed_columns(ref_coeff_cols, "form") + cpp_form_cols = sorted_prefixed_columns(cpp_coeff_cols, "form") + n_coeff_forms = min(len(ref_form_cols), len(cpp_form_cols)) + if n_coeff_forms == 0: + raise ValueError("No common harmonic coefficient form columns found.") + + ref_coeff = np.column_stack([ref_coeff_cols[name] for name in ref_form_cols[:n_coeff_forms]]) + cpp_coeff = np.column_stack([cpp_coeff_cols[name] for name in cpp_form_cols[:n_coeff_forms]]) + n_coeff_rows = min(ref_coeff.shape[0], cpp_coeff.shape[0]) + if ref_coeff.shape[0] != cpp_coeff.shape[0]: + warnings.append( + f"Harmonic coefficient row mismatch: reference={ref_coeff.shape[0]}, cpp={cpp_coeff.shape[0]}; cropped to {n_coeff_rows}." + ) + ref_coeff = ref_coeff[:n_coeff_rows] + cpp_coeff = cpp_coeff[:n_coeff_rows] + coeff_angles_deg = principal_angles_deg(ref_coeff, cpp_coeff, args.rcond) + coeff_procrustes_err = procrustes_relative_error(ref_coeff, cpp_coeff) + + ref_circ_cols = load_csv_columns(reference_dir / "reference_circular_coordinates.csv") + cpp_circ_cols = load_csv_columns(cpp_dir / "circular_coordinates.csv") + + ref_theta_cols = sorted_prefixed_columns(ref_circ_cols, "theta") + cpp_theta_cols = sorted_prefixed_columns(cpp_circ_cols, "theta") + n_theta = min(len(ref_theta_cols), len(cpp_theta_cols)) + if n_theta == 0: + raise ValueError("No common circular coordinate columns found.") + + circular_mode_rows: list[dict[str, Any]] = [] + circular_corr_values: list[float] = [] + circular_p95_values: list[float] = [] + + for idx in range(n_theta): + ref_col = ref_theta_cols[idx] + cpp_col = cpp_theta_cols[idx] + + theta_ref = ref_circ_cols[ref_col][:n_points] + theta_cpp = cpp_circ_cols[cpp_col][:n_points] + metrics = circular_mode_metrics(theta_ref, theta_cpp) + + row = { + "name": f"theta_{idx}", + "reference_column": ref_col, + "cpp_column": cpp_col, + **metrics, + } + circular_mode_rows.append(row) + circular_corr_values.append(float(metrics["complex_correlation"])) + circular_p95_values.append(float(metrics["wrapped_error_p95_rad"])) + + circular_corr_min = float(np.min(circular_corr_values)) + circular_corr_mean = safe_mean(circular_corr_values) + circular_p95_max = float(np.max(circular_p95_values)) + + score = composite_score( + harmonic_subspace_max_angle_deg=max_angle_deg, + harmonic_procrustes_rel_error=procrustes_err, + circular_complex_correlation_mean=circular_corr_mean, + circular_wrapped_p95_rad_max=circular_p95_max, + ) + + final_pass = ( + max_angle_deg <= args.final_max_angle_deg + and procrustes_err <= args.final_max_procrustes + and circular_corr_min >= args.final_min_correlation + and circular_p95_max <= args.final_max_p95_rad + ) + + improvement_vs_previous: float | None = None + commit_gate_pass = True + commit_gate_details = "baseline run (no previous composite score provided)" + + if args.previous_json: + previous_path = Path(args.previous_json) + if previous_path.exists(): + previous_payload = json.loads(previous_path.read_text()) + previous_score = float(previous_payload.get("summary", {}).get("composite_score", "nan")) + if np.isfinite(previous_score) and previous_score > 0.0: + improvement_vs_previous = (previous_score - score) / previous_score + commit_gate_pass = improvement_vs_previous >= args.gate_improvement + commit_gate_details = ( + f"improvement={100.0 * improvement_vs_previous:.2f}% (threshold={100.0 * args.gate_improvement:.2f}%)" + ) + else: + commit_gate_pass = False + commit_gate_details = "previous composite score missing or non-positive" + else: + commit_gate_pass = False + commit_gate_details = f"previous report not found: {previous_path}" + + payload: dict[str, Any] = { + "label": args.label, + "reference_dir": str(reference_dir), + "cpp_dir": str(cpp_dir), + "summary": { + "harmonic_subspace_max_angle_deg": max_angle_deg, + "harmonic_procrustes_rel_error": procrustes_err, + "circular_complex_correlation_min": circular_corr_min, + "circular_complex_correlation_mean": circular_corr_mean, + "circular_wrapped_p95_rad_max": circular_p95_max, + "composite_score": score, + }, + "harmonic_details": { + "n_forms_compared": len(common_form_indices), + "form_indices_compared": common_form_indices, + "ambient_principal_angles_deg": [float(v) for v in angles_deg.tolist()], + "ambient_procrustes_rel_error": procrustes_err, + "coeff_principal_angles_deg": [float(v) for v in coeff_angles_deg.tolist()], + "coeff_procrustes_rel_error": coeff_procrustes_err, + "n_coeff_rows_compared": n_coeff_rows, + }, + "circular_modes": circular_mode_rows, + "point_cloud": { + "n_points_compared": n_points, + "point_rmse": point_rmse, + "point_max_abs": point_max_abs, + }, + "gates": { + "final_pass": final_pass, + "commit_gate_pass": commit_gate_pass, + "commit_gate_details": commit_gate_details, + "improvement_vs_previous": improvement_vs_previous, + "thresholds": { + "final_max_angle_deg": args.final_max_angle_deg, + "final_max_procrustes": args.final_max_procrustes, + "final_min_correlation": args.final_min_correlation, + "final_max_p95_rad": args.final_max_p95_rad, + "commit_improvement": args.gate_improvement, + }, + }, + "warnings": warnings, + } + + output_json.parent.mkdir(parents=True, exist_ok=True) + output_markdown.parent.mkdir(parents=True, exist_ok=True) + output_json.write_text(json.dumps(payload, indent=2)) + output_markdown.write_text(render_markdown(payload)) + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/hodge/generate_input_torus.py b/scripts/hodge/generate_input_torus.py new file mode 100755 index 0000000..1f2c0fe --- /dev/null +++ b/scripts/hodge/generate_input_torus.py @@ -0,0 +1,36 @@ +#!/usr/bin/env python3 +import argparse +import csv +import pathlib +import numpy as np + + +def main() -> None: + parser = argparse.ArgumentParser(description="Generate deterministic torus point cloud CSV") + parser.add_argument("--output", required=True, help="Output CSV path") + parser.add_argument("--n", type=int, default=1000) + parser.add_argument("--major-radius", type=float, default=2.0) + parser.add_argument("--minor-radius", type=float, default=1.0) + parser.add_argument("--seed", type=int, default=0) + args = parser.parse_args() + + rng = np.random.default_rng(args.seed) + u = rng.uniform(0.0, 2.0 * np.pi, size=args.n) + v = rng.uniform(0.0, 2.0 * np.pi, size=args.n) + + x = (args.major_radius + args.minor_radius * np.cos(v)) * np.cos(u) + y = (args.major_radius + args.minor_radius * np.cos(v)) * np.sin(u) + z = args.minor_radius * np.sin(v) + + out_path = pathlib.Path(args.output) + out_path.parent.mkdir(parents=True, exist_ok=True) + + with out_path.open("w", newline="") as f: + writer = csv.writer(f) + writer.writerow(["x", "y", "z"]) + for row in np.stack([x, y, z], axis=1): + writer.writerow([float(row[0]), float(row[1]), float(row[2])]) + + +if __name__ == "__main__": + main() diff --git a/scripts/hodge/run_cpp_hodge.sh b/scripts/hodge/run_cpp_hodge.sh new file mode 100755 index 0000000..7d30a14 --- /dev/null +++ b/scripts/hodge/run_cpp_hodge.sh @@ -0,0 +1,40 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +ROOT_DIR="$(cd "${SCRIPT_DIR}/../.." && pwd)" + +INPUT_CSV="${1:?first arg must be input csv path}" +OUTPUT_DIR="${2:?second arg must be output dir path}" + +BUILD_DIR="${BUILD_DIR:-${ROOT_DIR}/build}" +EXE="${BUILD_DIR}/igneous-hodge" + +N_BASIS="${N_BASIS:-50}" +K_NEIGHBORS="${K_NEIGHBORS:-32}" +KNN_BANDWIDTH="${KNN_BANDWIDTH:-8}" +BANDWIDTH_VARIABILITY="${BANDWIDTH_VARIABILITY:--0.5}" +C_PARAM="${C_PARAM:-0.0}" +CIRCULAR_LAMBDA="${CIRCULAR_LAMBDA:-1.0}" +CIRCULAR_MODE_0="${CIRCULAR_MODE_0:-0}" +CIRCULAR_MODE_1="${CIRCULAR_MODE_1:-0}" + +if [[ ! -x "${EXE}" ]]; then + echo "Executable not found: ${EXE}" >&2 + exit 1 +fi + +mkdir -p "${OUTPUT_DIR}" + +"${EXE}" \ + --input-csv "${INPUT_CSV}" \ + --output-dir "${OUTPUT_DIR}" \ + --n-basis "${N_BASIS}" \ + --k-neighbors "${K_NEIGHBORS}" \ + --knn-bandwidth "${KNN_BANDWIDTH}" \ + --bandwidth-variability "${BANDWIDTH_VARIABILITY}" \ + --c "${C_PARAM}" \ + --circular-lambda "${CIRCULAR_LAMBDA}" \ + --circular-mode-0 "${CIRCULAR_MODE_0}" \ + --circular-mode-1 "${CIRCULAR_MODE_1}" \ + --no-ply diff --git a/scripts/hodge/run_parity_round.sh b/scripts/hodge/run_parity_round.sh new file mode 100755 index 0000000..247b7c9 --- /dev/null +++ b/scripts/hodge/run_parity_round.sh @@ -0,0 +1,115 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +ROOT_DIR="$(cd "${SCRIPT_DIR}/../.." && pwd)" + +RESULTS_ROOT="${RESULTS_ROOT:-${ROOT_DIR}/notes/hodge/results}" +ROUND_TAG="${ROUND_TAG:-$(date -u +%Y%m%d-%H%M%S)}" +ROUND_DIR="${RESULTS_ROOT}/round_${ROUND_TAG}" + +INPUT_DIR="${ROUND_DIR}/input" +REFERENCE_DIR="${ROUND_DIR}/reference" +CPP_DIR="${ROUND_DIR}/cpp" +REPORT_DIR="${ROUND_DIR}/report" +mkdir -p "${INPUT_DIR}" "${REFERENCE_DIR}" "${CPP_DIR}" "${REPORT_DIR}" + +VENV_PATH="${VENV_PATH:-${ROOT_DIR}/notes/hodge/.venv_ref}" +if [[ ! -x "${VENV_PATH}/bin/python" ]]; then + "${SCRIPT_DIR}/setup_reference_env.sh" "${VENV_PATH}" >/dev/null +fi +# shellcheck source=/dev/null +source "${VENV_PATH}/bin/activate" +PYTHON_BIN="${VENV_PATH}/bin/python" + +BUILD_DIR="${BUILD_DIR:-${ROOT_DIR}/build}" +if [[ ! -x "${BUILD_DIR}/igneous-hodge" ]]; then + cmake -S "${ROOT_DIR}" -B "${BUILD_DIR}" + cmake --build "${BUILD_DIR}" -j8 +fi + +N_POINTS="${N_POINTS:-1000}" +MAJOR_RADIUS="${MAJOR_RADIUS:-2.0}" +MINOR_RADIUS="${MINOR_RADIUS:-1.0}" +SEED="${SEED:-0}" + +N_FUNCTION_BASIS="${N_FUNCTION_BASIS:-50}" +N_COEFFICIENTS="${N_COEFFICIENTS:-50}" +K_NEIGHBORS="${K_NEIGHBORS:-32}" +KNN_BANDWIDTH="${KNN_BANDWIDTH:-8}" +BANDWIDTH_VARIABILITY="${BANDWIDTH_VARIABILITY:--0.5}" +C_PARAM="${C_PARAM:-0.0}" +CIRCULAR_LAMBDA="${CIRCULAR_LAMBDA:-1.0}" + +FORM_INDICES="${FORM_INDICES:-0,1}" +MODE_INDICES="${MODE_INDICES:-0,0}" + +INPUT_CSV="${INPUT_DIR}/torus.csv" +"${PYTHON_BIN}" "${SCRIPT_DIR}/generate_input_torus.py" \ + --output "${INPUT_CSV}" \ + --n "${N_POINTS}" \ + --major-radius "${MAJOR_RADIUS}" \ + --minor-radius "${MINOR_RADIUS}" \ + --seed "${SEED}" + +"${PYTHON_BIN}" "${SCRIPT_DIR}/run_reference_hodge.py" \ + --input-csv "${INPUT_CSV}" \ + --output-dir "${REFERENCE_DIR}" \ + --n-function-basis "${N_FUNCTION_BASIS}" \ + --n-coefficients "${N_COEFFICIENTS}" \ + --k-neighbors "${K_NEIGHBORS}" \ + --knn-bandwidth "${KNN_BANDWIDTH}" \ + --bandwidth-variability "${BANDWIDTH_VARIABILITY}" \ + --c "${C_PARAM}" \ + --circular-lambda "${CIRCULAR_LAMBDA}" \ + --form-indices "${FORM_INDICES}" \ + --mode-indices "${MODE_INDICES}" + +CIRCULAR_MODE_0="$(echo "${MODE_INDICES}" | cut -d',' -f1)" +CIRCULAR_MODE_1="$(echo "${MODE_INDICES}" | cut -d',' -f2)" + +N_BASIS="${N_FUNCTION_BASIS}" \ +K_NEIGHBORS="${K_NEIGHBORS}" \ +KNN_BANDWIDTH="${KNN_BANDWIDTH}" \ +BANDWIDTH_VARIABILITY="${BANDWIDTH_VARIABILITY}" \ +C_PARAM="${C_PARAM}" \ +CIRCULAR_LAMBDA="${CIRCULAR_LAMBDA}" \ +CIRCULAR_MODE_0="${CIRCULAR_MODE_0}" \ +CIRCULAR_MODE_1="${CIRCULAR_MODE_1}" \ +BUILD_DIR="${BUILD_DIR}" \ +"${SCRIPT_DIR}/run_cpp_hodge.sh" "${INPUT_CSV}" "${CPP_DIR}" + +COMPARE_ARGS=( + --reference-dir "${REFERENCE_DIR}" + --cpp-dir "${CPP_DIR}" + --output-markdown "${REPORT_DIR}/parity_report.md" + --output-json "${REPORT_DIR}/parity_report.json" + --label "Hodge parity round ${ROUND_TAG}" +) + +if [[ -n "${PREVIOUS_REPORT_JSON:-}" ]]; then + COMPARE_ARGS+=(--previous-json "${PREVIOUS_REPORT_JSON}") +fi + +"${PYTHON_BIN}" "${SCRIPT_DIR}/compare_hodge_outputs.py" "${COMPARE_ARGS[@]}" + +"${PYTHON_BIN}" - < np.ndarray: + data = np.genfromtxt(path, delimiter=",", names=True) + return np.column_stack([data["x"], data["y"], data["z"]]) + + +def write_table(path: pathlib.Path, header: List[str], rows: np.ndarray) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with path.open("w", newline="") as f: + writer = csv.writer(f) + writer.writerow(header) + for row in rows: + writer.writerow([float(v) for v in row]) + + +def circular_coordinates(form, lam: float, mode: int): + dg = form.dg + operator = form.sharp().operator - lam * dg.laplacian(0) + evals, efunctions = operator.spectrum() + + circular_indices = np.where(evals.imag > 0)[0] + if circular_indices.size == 0: + zeros = np.zeros(dg.n, dtype=float) + return zeros, complex(0.0, 0.0), -1 + + mode_idx = min(max(mode, 0), circular_indices.size - 1) + chosen_global = int(circular_indices[mode_idx]) + chosen_eval = evals[chosen_global] + + circular_funcs = efunctions[circular_indices].to_ambient() + chosen_func = circular_funcs[mode_idx] + angles = np.mod(np.arctan2(chosen_func.imag, chosen_func.real), 2.0 * np.pi) + return angles, chosen_eval, chosen_global + + +def main() -> None: + parser = argparse.ArgumentParser(description="Run reference DiffusionGeometry Hodge pipeline") + parser.add_argument("--input-csv", required=True) + parser.add_argument("--output-dir", required=True) + parser.add_argument("--n-function-basis", type=int, default=50) + parser.add_argument("--n-coefficients", type=int, default=50) + parser.add_argument("--k-neighbors", type=int, default=32) + parser.add_argument("--knn-bandwidth", type=int, default=8) + parser.add_argument("--bandwidth-variability", type=float, default=-0.5) + parser.add_argument("--c", type=float, default=0.0) + parser.add_argument("--circular-lambda", type=float, default=1.0) + parser.add_argument("--form-indices", default="0,1") + parser.add_argument("--mode-indices", default="0,0") + args = parser.parse_args() + + out_dir = pathlib.Path(args.output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + points = load_points(pathlib.Path(args.input_csv)) + + dg = DiffusionGeometry.from_point_cloud( + points, + n_function_basis=args.n_function_basis, + n_coefficients=args.n_coefficients, + knn_kernel=args.k_neighbors, + knn_bandwidth=args.knn_bandwidth, + bandwidth_variability=args.bandwidth_variability, + c=args.c, + use_mean_centres=True, + ) + + vals_1, forms_1 = dg.laplacian(1).spectrum() + + form_indices = [int(x.strip()) for x in args.form_indices.split(",") if x.strip()] + mode_indices = [int(x.strip()) for x in args.mode_indices.split(",") if x.strip()] + if len(mode_indices) < len(form_indices): + mode_indices.extend([0] * (len(form_indices) - len(mode_indices))) + + harmonic_coeffs = [] + harmonic_ambient = [] + circular_thetas = [] + circular_meta = [] + + for idx, form_idx in enumerate(form_indices): + form = forms_1[form_idx] + harmonic_coeffs.append(np.asarray(form.coeffs, dtype=float)) + ambient = np.asarray(form.to_ambient(), dtype=float) + harmonic_ambient.append(ambient) + + theta, eigval, eig_idx = circular_coordinates( + form, args.circular_lambda, mode_indices[idx] + ) + circular_thetas.append(theta) + circular_meta.append( + { + "name": f"theta_{idx}", + "form_index": int(form_idx), + "mode": int(mode_indices[idx]), + "selected_global_index": int(eig_idx), + "lambda": float(args.circular_lambda), + "eigenvalue_real": float(np.real(eigval)), + "eigenvalue_imag": float(np.imag(eigval)), + } + ) + + write_table( + out_dir / "reference_points.csv", + ["x", "y", "z"], + points, + ) + + spectrum_rows = np.column_stack([np.arange(vals_1.shape[0]), np.asarray(vals_1, dtype=float)]) + write_table( + out_dir / "reference_hodge_spectrum.csv", + ["mode", "lambda"], + spectrum_rows, + ) + + coeff_mat = np.column_stack(harmonic_coeffs) + coeff_rows = np.column_stack([np.arange(coeff_mat.shape[0]), coeff_mat]) + coeff_header = ["coeff_index"] + [f"form{i}" for i in range(coeff_mat.shape[1])] + write_table(out_dir / "reference_harmonic_coeffs.csv", coeff_header, coeff_rows) + + ambient_cols = [points] + ambient_header = ["x", "y", "z"] + for i, ambient in enumerate(harmonic_ambient): + ambient_cols.append(ambient) + ambient_header.extend([f"form{i}_x", f"form{i}_y", f"form{i}_z"]) + ambient_rows = np.column_stack(ambient_cols) + write_table(out_dir / "reference_harmonic_ambient.csv", ambient_header, ambient_rows) + + circular_rows = np.column_stack([points] + circular_thetas) + circular_header = ["x", "y", "z"] + [f"theta_{i}" for i in range(len(circular_thetas))] + write_table(out_dir / "reference_circular_coordinates.csv", circular_header, circular_rows) + + with (out_dir / "reference_circular_modes.csv").open("w", newline="") as f: + writer = csv.writer(f) + writer.writerow(["name", "form_index", "mode", "selected_global_index", "lambda", "eigenvalue_real", "eigenvalue_imag"]) + for row in circular_meta: + writer.writerow( + [ + row["name"], + row["form_index"], + row["mode"], + row["selected_global_index"], + row["lambda"], + row["eigenvalue_real"], + row["eigenvalue_imag"], + ] + ) + + metadata = { + "n_points": int(points.shape[0]), + "n_function_basis": int(args.n_function_basis), + "n_coefficients": int(args.n_coefficients), + "k_neighbors": int(args.k_neighbors), + "knn_bandwidth": int(args.knn_bandwidth), + "bandwidth_variability": float(args.bandwidth_variability), + "c": float(args.c), + "circular_lambda": float(args.circular_lambda), + "form_indices": form_indices, + "mode_indices": mode_indices, + } + (out_dir / "reference_metadata.json").write_text(json.dumps(metadata, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/scripts/hodge/setup_reference_env.sh b/scripts/hodge/setup_reference_env.sh new file mode 100755 index 0000000..541a451 --- /dev/null +++ b/scripts/hodge/setup_reference_env.sh @@ -0,0 +1,30 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +ROOT_DIR="$(cd "${SCRIPT_DIR}/../.." && pwd)" + +VENV_PATH="${1:-${ROOT_DIR}/notes/hodge/.venv_ref}" +PYTHON_BIN="${PYTHON_BIN:-/opt/homebrew/bin/python3.13}" + +if [[ ! -x "${PYTHON_BIN}" ]]; then + if command -v python3.13 >/dev/null 2>&1; then + PYTHON_BIN="$(command -v python3.13)" + else + echo "python3.13 not found. Install Python 3.13 or set PYTHON_BIN." >&2 + exit 1 + fi +fi + +if [[ ! -d "${VENV_PATH}" ]]; then + "${PYTHON_BIN}" -m venv "${VENV_PATH}" +fi + +# shellcheck source=/dev/null +source "${VENV_PATH}/bin/activate" + +python -m pip -q install --upgrade pip +python -m pip -q install numpy scipy scikit-learn opt_einsum +python -m pip -q install -e "${ROOT_DIR}/DiffusionGeometry" + +echo "${VENV_PATH}" diff --git a/src/main_diffusion.cpp b/src/main_diffusion.cpp index e19330a..793db18 100644 --- a/src/main_diffusion.cpp +++ b/src/main_diffusion.cpp @@ -42,7 +42,7 @@ int main(int argc, char **argv) { ops::normalize(mesh); mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), 0.005f, 32}); + mesh.geometry.z_span(), 32}); const int n = static_cast(mesh.geometry.num_points()); std::cout << "Markov Chain P: " << n << "x" << n << " (" diff --git a/src/main_hodge.cpp b/src/main_hodge.cpp index 762bf21..9509f4d 100644 --- a/src/main_hodge.cpp +++ b/src/main_hodge.cpp @@ -1,10 +1,13 @@ #include +#include #include #include #include #include #include #include +#include +#include #include #include @@ -17,65 +20,262 @@ using namespace igneous; using DiffusionMesh = data::Mesh; -static void generate_torus(DiffusionMesh &mesh, size_t n_points, float R, - float r) { +struct Config { + std::string input_csv; + std::string output_dir = "output_hodge"; + size_t n_points = 1000; + float major_radius = 2.0f; + float minor_radius = 1.0f; + uint32_t seed = 0; + + int n_basis = 50; + int k_neighbors = 32; + int knn_bandwidth = 8; + float bandwidth_variability = -0.5f; + float c = 0.0f; + + float circular_lambda = 1.0f; + int circular_mode_0 = 0; + int circular_mode_1 = 0; + + bool export_ply = true; +}; + +static void print_usage() { + std::cout << "Usage: ./build/igneous-hodge [options]\n" + << " --input-csv \n" + << " --output-dir \n" + << " --n-points \n" + << " --major-radius \n" + << " --minor-radius \n" + << " --seed \n" + << " --n-basis \n" + << " --k-neighbors \n" + << " --knn-bandwidth \n" + << " --bandwidth-variability \n" + << " --c \n" + << " --circular-lambda \n" + << " --circular-mode-0 \n" + << " --circular-mode-1 \n" + << " --no-ply\n"; +} + +static bool parse_args(int argc, char **argv, Config &cfg) { + for (int i = 1; i < argc; ++i) { + const std::string arg = argv[i]; + auto require_value = [&](const char *name) -> const char * { + if (i + 1 >= argc) { + std::cerr << "Missing value for " << name << "\n"; + std::exit(1); + } + ++i; + return argv[i]; + }; + + if (arg == "--help" || arg == "-h") { + print_usage(); + return false; + } + if (arg == "--input-csv") { + cfg.input_csv = require_value("--input-csv"); + continue; + } + if (arg == "--output-dir") { + cfg.output_dir = require_value("--output-dir"); + continue; + } + if (arg == "--n-points") { + cfg.n_points = static_cast(std::stoul(require_value("--n-points"))); + continue; + } + if (arg == "--major-radius") { + cfg.major_radius = std::stof(require_value("--major-radius")); + continue; + } + if (arg == "--minor-radius") { + cfg.minor_radius = std::stof(require_value("--minor-radius")); + continue; + } + if (arg == "--seed") { + cfg.seed = static_cast(std::stoul(require_value("--seed"))); + continue; + } + if (arg == "--n-basis") { + cfg.n_basis = std::stoi(require_value("--n-basis")); + continue; + } + if (arg == "--k-neighbors") { + cfg.k_neighbors = std::stoi(require_value("--k-neighbors")); + continue; + } + if (arg == "--knn-bandwidth") { + cfg.knn_bandwidth = std::stoi(require_value("--knn-bandwidth")); + continue; + } + if (arg == "--bandwidth-variability") { + cfg.bandwidth_variability = std::stof(require_value("--bandwidth-variability")); + continue; + } + if (arg == "--c") { + cfg.c = std::stof(require_value("--c")); + continue; + } + if (arg == "--circular-lambda") { + cfg.circular_lambda = std::stof(require_value("--circular-lambda")); + continue; + } + if (arg == "--circular-mode-0") { + cfg.circular_mode_0 = std::stoi(require_value("--circular-mode-0")); + continue; + } + if (arg == "--circular-mode-1") { + cfg.circular_mode_1 = std::stoi(require_value("--circular-mode-1")); + continue; + } + if (arg == "--no-ply") { + cfg.export_ply = false; + continue; + } + + std::cerr << "Unknown argument: " << arg << "\n"; + print_usage(); + return false; + } + + return true; +} + +static void generate_torus(DiffusionMesh &mesh, size_t n_points, float major_radius, + float minor_radius, uint32_t seed) { mesh.clear(); mesh.geometry.reserve(n_points); - std::mt19937 gen(42); - std::uniform_real_distribution dist(0.0f, 6.283185f); + std::mt19937 gen(seed); + std::uniform_real_distribution dist(0.0f, 6.28318530718f); for (size_t i = 0; i < n_points; ++i) { const float u = dist(gen); const float v = dist(gen); + const float x = (major_radius + minor_radius * std::cos(v)) * std::cos(u); + const float y = (major_radius + minor_radius * std::cos(v)) * std::sin(u); + const float z = minor_radius * std::sin(v); + mesh.geometry.push_point({x, y, z}); + } +} + +static bool load_point_cloud_csv(const std::string &filename, DiffusionMesh &mesh) { + std::ifstream file(filename); + if (!file.is_open()) { + std::cerr << "Failed to open input CSV: " << filename << "\n"; + return false; + } + + mesh.clear(); + std::string line; + while (std::getline(file, line)) { + if (line.empty()) { + continue; + } - const float x = (R + r * std::cos(v)) * std::cos(u); - const float y = (R + r * std::cos(v)) * std::sin(u); - const float z = r * std::sin(v); + for (char &c : line) { + if (c == ',') { + c = ' '; + } + } + std::istringstream iss(line); + float x = 0.0f; + float y = 0.0f; + float z = 0.0f; + if (!(iss >> x >> y >> z)) { + continue; + } mesh.geometry.push_point({x, y, z}); } + + return mesh.geometry.num_points() > 0; +} + +static void write_points_csv(const std::string &filename, const DiffusionMesh &mesh) { + std::ofstream file(filename); + file << "x,y,z\n"; + const size_t n = mesh.geometry.num_points(); + for (size_t i = 0; i < n; ++i) { + const auto p = mesh.geometry.get_vec3(i); + file << p.x << "," << p.y << "," << p.z << "\n"; + } +} + +static void write_spectrum_csv(const std::string &filename, + const Eigen::VectorXf &evals) { + std::ofstream file(filename); + file << "mode,lambda\n"; + for (int i = 0; i < evals.size(); ++i) { + file << i << "," << evals[i] << "\n"; + } +} + +static void write_harmonic_coeffs_csv(const std::string &filename, + const Eigen::MatrixXf &evecs, + int n_forms) { + std::ofstream file(filename); + file << "coeff_index"; + for (int form = 0; form < n_forms; ++form) { + file << ",form" << form; + } + file << "\n"; + + for (int i = 0; i < evecs.rows(); ++i) { + file << i; + for (int form = 0; form < n_forms; ++form) { + file << "," << evecs(i, form); + } + file << "\n"; + } } static std::vector -reconstruct_vector_field(const DiffusionMesh &mesh, - const Eigen::VectorXf &coeffs, float bandwidth) { +reconstruct_harmonic_ambient(const DiffusionMesh &mesh, + const Eigen::VectorXf &coeffs) { const size_t n_verts = mesh.geometry.num_points(); const int n_basis = mesh.topology.eigen_basis.cols(); + const auto &U = mesh.topology.eigen_basis; - std::array coords; - ops::fill_coordinate_vectors(mesh, coords); + std::array data_coords; + std::array immersion_coords; + ops::fill_data_coordinate_vectors(mesh, data_coords); + ops::fill_coordinate_vectors(mesh, immersion_coords); - std::array, 3> gamma{}; + std::array, 3> gamma_data_imm{}; for (int a = 0; a < 3; ++a) { for (int b = 0; b < 3; ++b) { - gamma[a][b].resize(static_cast(n_verts)); - ops::carre_du_champ(mesh, coords[a], coords[b], bandwidth, gamma[a][b]); + gamma_data_imm[a][b].resize(static_cast(n_verts)); + ops::carre_du_champ(mesh, data_coords[a], immersion_coords[b], 0.0f, + gamma_data_imm[a][b]); } } - std::vector field(n_verts, {0.0f, 0.0f, 0.0f}); - const auto &U = mesh.topology.eigen_basis; - + Eigen::MatrixXf alpha_mat(n_basis, 3); for (int k = 0; k < n_basis; ++k) { - for (int a = 0; a < 3; ++a) { - const float c = coeffs(k * 3 + a); - if (std::abs(c) < 1e-7f) { - continue; - } + alpha_mat(k, 0) = coeffs(k * 3 + 0); + alpha_mat(k, 1) = coeffs(k * 3 + 1); + alpha_mat(k, 2) = coeffs(k * 3 + 2); + } + const Eigen::MatrixXf q = U * alpha_mat; // [n x 3] - for (int b = 0; b < 3; ++b) { - for (size_t i = 0; i < n_verts; ++i) { - const float val = c * U(static_cast(i), k) * gamma[a][b][static_cast(i)]; - if (b == 0) - field[i].x += val; - else if (b == 1) - field[i].y += val; - else - field[i].z += val; - } - } - } + std::vector field(n_verts, {0.0f, 0.0f, 0.0f}); + for (size_t i = 0; i < n_verts; ++i) { + const int idx = static_cast(i); + const float vx = gamma_data_imm[0][0][idx] * q(idx, 0) + + gamma_data_imm[0][1][idx] * q(idx, 1) + + gamma_data_imm[0][2][idx] * q(idx, 2); + const float vy = gamma_data_imm[1][0][idx] * q(idx, 0) + + gamma_data_imm[1][1][idx] * q(idx, 1) + + gamma_data_imm[1][2][idx] * q(idx, 2); + const float vz = gamma_data_imm[2][0][idx] * q(idx, 0) + + gamma_data_imm[2][1][idx] * q(idx, 1) + + gamma_data_imm[2][2][idx] * q(idx, 2); + field[i] = {vx, vy, vz}; } return field; @@ -101,36 +301,139 @@ static void export_vector_field(const std::string &filename, } } -int main() { - const bool bench_mode = std::getenv("IGNEOUS_BENCH_MODE") != nullptr; - if (!bench_mode) { - std::filesystem::create_directory("output_hodge"); +static void write_harmonic_ambient_csv(const std::string &filename, + const DiffusionMesh &mesh, + const std::vector> &fields) { + std::ofstream file(filename); + file << "x,y,z"; + for (size_t form = 0; form < fields.size(); ++form) { + file << ",form" << form << "_x" + << ",form" << form << "_y" + << ",form" << form << "_z"; } + file << "\n"; + + const size_t n = mesh.geometry.num_points(); + for (size_t i = 0; i < n; ++i) { + const auto p = mesh.geometry.get_vec3(i); + file << p.x << "," << p.y << "," << p.z; + for (size_t form = 0; form < fields.size(); ++form) { + const auto v = fields[form][i]; + file << "," << v.x << "," << v.y << "," << v.z; + } + file << "\n"; + } +} + +static void write_circular_csv(const std::string &filename, + const DiffusionMesh &mesh, + const Eigen::VectorXf &theta_0, + const Eigen::VectorXf &theta_1) { + std::ofstream file(filename); + file << "x,y,z,theta_0,theta_1\n"; + const size_t n = mesh.geometry.num_points(); + for (size_t i = 0; i < n; ++i) { + const auto p = mesh.geometry.get_vec3(i); + const int idx = static_cast(i); + file << p.x << "," << p.y << "," << p.z << "," << theta_0[idx] << "," + << theta_1[idx] << "\n"; + } +} + +static void write_circular_modes_csv(const std::string &filename, + std::complex eval_0, + std::complex eval_1, + int mode_0, int mode_1, + float circular_lambda) { + std::ofstream file(filename); + file << "name,mode,lambda,eigenvalue_real,eigenvalue_imag\n"; + file << "theta_0," << mode_0 << "," << circular_lambda << "," << eval_0.real() + << "," << eval_0.imag() << "\n"; + file << "theta_1," << mode_1 << "," << circular_lambda << "," << eval_1.real() + << "," << eval_1.imag() << "\n"; +} + +int main(int argc, char **argv) { + Config cfg; + if (!parse_args(argc, argv, cfg)) { + return 0; + } + + const bool bench_mode = std::getenv("IGNEOUS_BENCH_MODE") != nullptr; + const bool export_ply = cfg.export_ply && !bench_mode; + + std::filesystem::create_directories(cfg.output_dir); DiffusionMesh mesh; - generate_torus(mesh, 4000, 2.0f, 0.8f); + if (!cfg.input_csv.empty()) { + if (!load_point_cloud_csv(cfg.input_csv, mesh)) { + return 1; + } + } else { + generate_torus(mesh, cfg.n_points, cfg.major_radius, cfg.minor_radius, cfg.seed); + } - const float bandwidth = 0.05f; - mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), bandwidth, 32}); + mesh.topology.build({mesh.geometry.x_span(), + mesh.geometry.y_span(), + mesh.geometry.z_span(), + cfg.k_neighbors, + cfg.knn_bandwidth, + cfg.bandwidth_variability, + cfg.c, + true}); - const int n_basis = 64; - ops::compute_eigenbasis(mesh, n_basis); + ops::compute_eigenbasis(mesh, cfg.n_basis); ops::GeometryWorkspace geom_ws; - const auto G = ops::compute_1form_gram_matrix(mesh, bandwidth, geom_ws); + const auto G = ops::compute_1form_gram_matrix(mesh, 0.0f, geom_ws); ops::HodgeWorkspace hodge_ws; - const auto D_weak = ops::compute_weak_exterior_derivative(mesh, bandwidth, hodge_ws); - const auto E_up = ops::compute_curl_energy_matrix(mesh, bandwidth, hodge_ws); + const auto D_weak = ops::compute_weak_exterior_derivative(mesh, 0.0f, hodge_ws); + const auto E_up = ops::compute_curl_energy_matrix(mesh, 0.0f, hodge_ws); const auto laplacian = ops::compute_hodge_laplacian_matrix(D_weak, E_up); auto [evals, evecs] = ops::compute_hodge_spectrum(laplacian, G); + if (evals.size() == 0 || evecs.cols() < 2) { + std::cerr << "Hodge spectrum solve failed.\n"; + return 1; + } - const auto theta_0 = ops::compute_circular_coordinates(mesh, evecs.col(0), bandwidth); - const auto theta_1 = ops::compute_circular_coordinates(mesh, evecs.col(1), bandwidth); + std::complex selected_eval_0(0.0f, 0.0f); + std::complex selected_eval_1(0.0f, 0.0f); + const auto theta_0 = ops::compute_circular_coordinates( + mesh, evecs.col(0), 0.0f, cfg.circular_lambda, cfg.circular_mode_0, + &selected_eval_0); + const auto theta_1 = ops::compute_circular_coordinates( + mesh, evecs.col(1), 0.0f, cfg.circular_lambda, cfg.circular_mode_1, + &selected_eval_1); - if (!bench_mode) { + std::cout << "HODGE SPECTRUM (first 12 modes)\n"; + for (int i = 0; i < 12 && i < evals.size(); ++i) { + std::cout << "Mode " << i << ": lambda = " << evals[i] << "\n"; + } + + write_points_csv(cfg.output_dir + "/points.csv", mesh); + write_spectrum_csv(cfg.output_dir + "/hodge_spectrum.csv", evals); + + const int export_forms = std::min(3, static_cast(evecs.cols())); + write_harmonic_coeffs_csv(cfg.output_dir + "/harmonic_coeffs.csv", evecs, + export_forms); + + std::vector> harmonic_fields; + harmonic_fields.reserve(static_cast(export_forms)); + for (int i = 0; i < export_forms; ++i) { + harmonic_fields.push_back(reconstruct_harmonic_ambient(mesh, evecs.col(i))); + } + write_harmonic_ambient_csv(cfg.output_dir + "/harmonic_ambient.csv", mesh, + harmonic_fields); + + write_circular_csv(cfg.output_dir + "/circular_coordinates.csv", mesh, theta_0, + theta_1); + write_circular_modes_csv(cfg.output_dir + "/circular_modes.csv", selected_eval_0, + selected_eval_1, cfg.circular_mode_0, + cfg.circular_mode_1, cfg.circular_lambda); + + if (export_ply) { std::vector field_0(static_cast(theta_0.size())); std::vector field_1(static_cast(theta_1.size())); for (int i = 0; i < theta_0.size(); ++i) { @@ -138,21 +441,15 @@ int main() { field_1[static_cast(i)] = theta_1[i]; } - io::export_ply_solid(mesh, field_0, "output_hodge/torus_angle_0.ply", 0.01); - io::export_ply_solid(mesh, field_1, "output_hodge/torus_angle_1.ply", 0.01); - } - - std::cout << "HODGE SPECTRUM (first 12 modes)\n"; - for (int i = 0; i < 12 && i < evals.size(); ++i) { - std::cout << "Mode " << i << ": lambda = " << evals[i] << "\n"; - } + io::export_ply_solid(mesh, field_0, + cfg.output_dir + "/torus_angle_0.ply", 0.01); + io::export_ply_solid(mesh, field_1, + cfg.output_dir + "/torus_angle_1.ply", 0.01); - if (!bench_mode) { - for (int i = 0; i < 3; ++i) { - const Eigen::VectorXf coeffs = evecs.col(i); - const auto vector_field = reconstruct_vector_field(mesh, coeffs, bandwidth); - const std::string fname = std::format("output_hodge/harmonic_form_{}.ply", i); - export_vector_field(fname, mesh, vector_field); + for (int i = 0; i < export_forms; ++i) { + const std::string fname = + std::format("{}/harmonic_form_{}.ply", cfg.output_dir, i); + export_vector_field(fname, mesh, harmonic_fields[static_cast(i)]); } } diff --git a/src/main_spectral.cpp b/src/main_spectral.cpp index 420921a..eac5ae3 100644 --- a/src/main_spectral.cpp +++ b/src/main_spectral.cpp @@ -30,7 +30,7 @@ int main(int argc, char **argv) { const float bandwidth = 0.005f; mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), bandwidth, 32}); + mesh.geometry.z_span(), 32}); const int n_basis = 16; ops::compute_eigenbasis(mesh, n_basis); diff --git a/tests/test_ops_hodge.cpp b/tests/test_ops_hodge.cpp index 7f97a48..bcb2424 100644 --- a/tests/test_ops_hodge.cpp +++ b/tests/test_ops_hodge.cpp @@ -25,7 +25,8 @@ make_torus_cloud(size_t n_points) { mesh.geometry.push_point({x, y, z}); } - mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), mesh.geometry.z_span(), 0.05f, 24}); + mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), + mesh.geometry.z_span(), 24}); return mesh; } @@ -57,11 +58,12 @@ TEST_CASE("Hodge operators produce finite spectrum and circular coordinates") { CHECK(std::isfinite(evals[i])); } - const auto theta = igneous::ops::compute_circular_coordinates(mesh, evecs.col(0), 0.05f); + const auto theta = + igneous::ops::compute_circular_coordinates(mesh, evecs.col(0), 0.0f); CHECK(theta.size() == static_cast(mesh.geometry.num_points())); for (int i = 0; i < theta.size(); ++i) { CHECK(std::isfinite(theta[i])); CHECK(theta[i] >= -0.01f); - CHECK(theta[i] <= 1.01f); + CHECK(theta[i] <= 6.30f); } } diff --git a/tests/test_ops_spectral_geometry.cpp b/tests/test_ops_spectral_geometry.cpp index 605c4ac..89f8fb4 100644 --- a/tests/test_ops_spectral_geometry.cpp +++ b/tests/test_ops_spectral_geometry.cpp @@ -18,7 +18,8 @@ make_diffusion_cloud(size_t n_points) { mesh.geometry.push_point({std::cos(t * 6.283185f), std::sin(t * 6.283185f), t}); } - mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), mesh.geometry.z_span(), 0.05f, 24}); + mesh.topology.build( + {mesh.geometry.x_span(), mesh.geometry.y_span(), mesh.geometry.z_span(), 24}); return mesh; } diff --git a/tests/test_topology_diffusion.cpp b/tests/test_topology_diffusion.cpp index 1b1d81e..4b539e3 100644 --- a/tests/test_topology_diffusion.cpp +++ b/tests/test_topology_diffusion.cpp @@ -41,7 +41,7 @@ TEST_CASE("DiffusionTopology produces stochastic Markov matrix and valid measure } igneous::data::DiffusionTopology topo; - topo.build({geometry.x_span(), geometry.y_span(), geometry.z_span(), 0.05f, 8}); + topo.build({geometry.x_span(), geometry.y_span(), geometry.z_span(), 8}); CHECK(topo.markov_row_offsets.size() == 17); CHECK(topo.markov_col_indices.size() == topo.markov_values.size()); @@ -85,7 +85,7 @@ TEST_CASE("Diffusion CSR markov step matches sparse matrix product") { } mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), 0.05f, 8}); + mesh.geometry.z_span(), 8}); const int n = static_cast(mesh.geometry.num_points()); Eigen::VectorXf u = Eigen::VectorXf::LinSpaced(n, -1.0f, 1.0f); @@ -113,7 +113,7 @@ TEST_CASE("Diffusion multi-step markov matches repeated single steps") { } mesh.topology.build({mesh.geometry.x_span(), mesh.geometry.y_span(), - mesh.geometry.z_span(), 0.05f, 8}); + mesh.geometry.z_span(), 8}); const int n = static_cast(mesh.geometry.num_points()); Eigen::VectorXf u0 = Eigen::VectorXf::LinSpaced(n, -1.0f, 1.0f); From 3c07519df86195e516050a793117ecea90b3c11d Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 09:32:42 -0700 Subject: [PATCH 3/9] Determinize reference harness and normalize spectral solve --- include/igneous/ops/spectral.hpp | 72 +++++++++++++++++++++++++--- notes/hodge/journal.md | 62 ++++++++++++++++++++++++ scripts/hodge/run_reference_hodge.py | 69 +++++++++++++++++++++++--- 3 files changed, 190 insertions(+), 13 deletions(-) diff --git a/include/igneous/ops/spectral.hpp b/include/igneous/ops/spectral.hpp index 0fd87d7..3f8002e 100644 --- a/include/igneous/ops/spectral.hpp +++ b/include/igneous/ops/spectral.hpp @@ -162,6 +162,61 @@ class SymmetricKernelCsrMatProd { int n_ = 0; }; +class NormalizedSymmetricKernelCsrMatProd { +public: + using Scalar = float; + + NormalizedSymmetricKernelCsrMatProd(const std::vector &row_offsets, + const std::vector &col_indices, + const std::vector &values, + const Eigen::VectorXf &inv_sqrt_rows, + int n) + : row_offsets_(row_offsets), col_indices_(col_indices), values_(values), + inv_sqrt_rows_(inv_sqrt_rows), n_(n) {} + + [[nodiscard]] int rows() const { return n_; } + [[nodiscard]] int cols() const { return n_; } + + void perform_op(const Scalar *x_in, Scalar *y_out) const { + const float *inv_sqrt_data = inv_sqrt_rows_.data(); + core::parallel_for_index( + 0, n_, + [&](int i) { + const int begin = row_offsets_[static_cast(i)]; + const int end = row_offsets_[static_cast(i) + 1]; + const int count = end - begin; + const int *cols = col_indices_.data() + begin; + const float *vals = values_.data() + begin; + + float acc = 0.0f; + int k = 0; + for (; k + 3 < count; k += 4) { + const int j0 = cols[k + 0]; + const int j1 = cols[k + 1]; + const int j2 = cols[k + 2]; + const int j3 = cols[k + 3]; + acc += vals[k + 0] * (x_in[j0] * inv_sqrt_data[j0]); + acc += vals[k + 1] * (x_in[j1] * inv_sqrt_data[j1]); + acc += vals[k + 2] * (x_in[j2] * inv_sqrt_data[j2]); + acc += vals[k + 3] * (x_in[j3] * inv_sqrt_data[j3]); + } + for (; k < count; ++k) { + const int j = cols[k]; + acc += vals[k] * (x_in[j] * inv_sqrt_data[j]); + } + y_out[i] = inv_sqrt_data[i] * acc; + }, + 8192); + } + +private: + const std::vector &row_offsets_; + const std::vector &col_indices_; + const std::vector &values_; + const Eigen::VectorXf &inv_sqrt_rows_; + int n_ = 0; +}; + // Computes the first k eigenvectors of the Markov Chain P. template void compute_eigenbasis(MeshT &mesh, int n_eigenvectors) { @@ -263,10 +318,17 @@ void compute_eigenbasis(MeshT &mesh, int n_eigenvectors) { const int compact_ncv = try_compact ? std::min(n, std::max(k_eval + 16, 20)) : full_ncv; - SymmetricKernelCsrMatProd op(sym_row_offsets, sym_col_indices, sym_values, - n); + Eigen::VectorXf inv_sqrt_rows(n); + for (int i = 0; i < n; ++i) { + inv_sqrt_rows[i] = 1.0f / std::sqrt(std::max(row_sums[i], 1e-12f)); + } + + NormalizedSymmetricKernelCsrMatProd op(sym_row_offsets, sym_col_indices, + sym_values, inv_sqrt_rows, n); const auto solve_symmetric = [&](int ncv) -> bool { - Spectra::SymEigsSolver eigs(op, k_eval, ncv); + Spectra::SymEigsSolver eigs(op, + k_eval, + ncv); eigs.init(); const int nconv = eigs.compute(Spectra::SortRule::LargestMagn); if (eigs.info() != Spectra::CompInfo::Successful || nconv <= 0) { @@ -274,10 +336,6 @@ void compute_eigenbasis(MeshT &mesh, int n_eigenvectors) { } Eigen::MatrixXf basis = eigs.eigenvectors(nconv); - Eigen::VectorXf inv_sqrt_rows(n); - for (int i = 0; i < n; ++i) { - inv_sqrt_rows[i] = 1.0f / std::sqrt(std::max(row_sums[i], 1e-12f)); - } basis = basis.array().colwise() * inv_sqrt_rows.array(); basis = basis.rowwise().reverse().eval(); if (std::abs(basis(0, 0)) > 1e-12f) { diff --git a/notes/hodge/journal.md b/notes/hodge/journal.md index 1802af7..4faa9ca 100644 --- a/notes/hodge/journal.md +++ b/notes/hodge/journal.md @@ -98,3 +98,65 @@ Use one entry per parity hypothesis. - Notes: - This is a structural parity fix in eigenmode selection policy, not a numerical micro-optimization. - Remaining gaps are now concentrated in harmonic-form thresholds and circular P95 tail. + +## 2026-02-14 Hypothesis H1: Deterministic Reference Basis in Harness +- Timestamp: 2026-02-14T16:26:56Z +- Commit: `388e076` + working tree (uncommitted) +- Hypothesis: Reference parity runs are unstable because the reference eigensolver initialisation is random; precomputing a deterministic function basis from the same reference kernel path should stabilise comparisons. +- Files touched: + - `scripts/hodge/run_reference_hodge.py` +- Commands: + - `scripts/hodge/run_parity_round.sh` (repeated) + - Determinism check: + - repeated `scripts/hodge/run_reference_hodge.py` runs and checksum comparison on `reference_harmonic_coeffs.csv` +- Numeric parity deltas: + - Deterministic baseline report: `notes/hodge/results/round_20260214-162656/report/parity_report.json` + - Composite score: `0.156485` + - Harmonic subspace max principal angle: `6.493845 deg` + - Harmonic Procrustes error: `0.098484` + - Circular correlation min: `0.714689` + - Circular wrapped P95 max: `1.572800 rad` +- Decision: `kept` +- Notes: + - Reference outputs are now repeatable across runs for identical inputs/parameters. + - This intentionally re-baselines parity metrics for reliable hypothesis gating. + +## 2026-02-14 Hypothesis A2: Normalized Symmetric-Kernel Eigensolve (Deterministic Baseline) +- Timestamp: 2026-02-14T16:28:53Z +- Commit: `388e076` + working tree (uncommitted) +- Hypothesis: Under deterministic reference outputs, solving the normalized symmetric operator `D^{-1/2} K D^{-1/2}` should materially reduce harmonic mismatch. +- Files touched: + - `include/igneous/ops/spectral.hpp` +- Commands: + - `cmake --build build -j8` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-162656/report/parity_report.json scripts/hodge/run_parity_round.sh` + - Stability check: repeated `scripts/hodge/run_parity_round.sh` runs +- Numeric parity deltas: + - Report: `notes/hodge/results/round_20260214-162853/report/parity_report.json` + - Composite: `0.066904` (`-57.25%` vs deterministic baseline) + - Harmonic subspace max principal angle: `6.376378 deg` + - Harmonic Procrustes error: `0.097143` + - Circular correlation min: `0.996571` (major improvement) + - Circular wrapped P95 max: `0.163400 rad` (major improvement) +- Decision: `kept` +- Notes: + - Cleared the 20% keep threshold with deterministic artifacts. + - Harmonic thresholds remain the dominant blocker. + +## 2026-02-14 Hypothesis H2: Double-Precision Hodge/Circular Solves +- Timestamp: 2026-02-14T16:30:47Z +- Commit: `388e076` + working tree (uncommitted) +- Hypothesis: Promoting Hodge spectral and circular generalized eigensolves to double precision will reduce harmonic/circular error tails. +- Files touched: + - `include/igneous/ops/hodge.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-162853/report/parity_report.json scripts/hodge/run_parity_round.sh` +- Numeric parity deltas: + - Report: `notes/hodge/results/round_20260214-163047/report/parity_report.json` + - Composite: `0.159232` (`+138.00%` regression) + - Circular correlation min dropped to `0.714585` + - Circular wrapped P95 max regressed to `1.566015 rad` +- Decision: `rejected` +- Notes: + - Severe regression; reverted immediately. diff --git a/scripts/hodge/run_reference_hodge.py b/scripts/hodge/run_reference_hodge.py index 04458c9..1df8249 100755 --- a/scripts/hodge/run_reference_hodge.py +++ b/scripts/hodge/run_reference_hodge.py @@ -6,8 +6,15 @@ from typing import List import numpy as np +from scipy.sparse import diags +from scipy.sparse.linalg import eigsh from diffusion_geometry import DiffusionGeometry +from diffusion_geometry.core import ( + build_symmetric_kernel_matrix, + knn_graph, + markov_chain, +) def load_points(path: pathlib.Path) -> np.ndarray: @@ -24,6 +31,40 @@ def write_table(path: pathlib.Path, header: List[str], rows: np.ndarray) -> None writer.writerow([float(v) for v in row]) +def compute_deterministic_basis( + kernel: np.ndarray, nbr_indices: np.ndarray, n_function_basis: int +) -> tuple[np.ndarray, np.ndarray]: + K_sym, row_sums = build_symmetric_kernel_matrix(kernel, nbr_indices) + row_sums = np.asarray(row_sums, dtype=float) + inv_sqrt = np.power(np.maximum(row_sums, 1e-12), -0.5) + K_normalized = diags(inv_sqrt) @ K_sym @ diags(inv_sqrt) + + n = K_normalized.shape[0] + n0 = int(min(max(1, n_function_basis), n)) + + if n0 < n: + _, eigenfunctions = eigsh( + K_normalized, + n0, + which="LM", + tol=1e-2, + v0=np.ones(n, dtype=float), + ) + else: + dense = K_normalized.toarray() + evals, evecs = np.linalg.eigh(dense) + top_idx = np.argsort(np.abs(evals))[-n0:] + eigenfunctions = evecs[:, top_idx[::-1]] + + basis = inv_sqrt[:, None] * eigenfunctions + basis = basis[:, ::-1] + if abs(basis[0, 0]) > 1e-12: + basis = basis / basis[0, 0] + + measure = row_sums / np.maximum(row_sums.sum(), 1e-12) + return basis, measure + + def circular_coordinates(form, lam: float, mode: int): dg = form.dg operator = form.sharp().operator - lam * dg.laplacian(0) @@ -63,16 +104,32 @@ def main() -> None: out_dir.mkdir(parents=True, exist_ok=True) points = load_points(pathlib.Path(args.input_csv)) + nbr_distances, nbr_indices = knn_graph(points, knn_kernel=args.k_neighbors) + kernel, bandwidths = markov_chain( + nbr_distances=nbr_distances, + nbr_indices=nbr_indices, + c=args.c, + bandwidth_variability=args.bandwidth_variability, + knn_bandwidth=args.knn_bandwidth, + ) + function_basis, measure = compute_deterministic_basis( + kernel=kernel, + nbr_indices=nbr_indices, + n_function_basis=args.n_function_basis, + ) - dg = DiffusionGeometry.from_point_cloud( - points, + dg = DiffusionGeometry.from_knn_kernel( + nbr_indices=nbr_indices, + kernel=kernel, + immersion_coords=points, + data_matrix=points, + bandwidths=bandwidths, n_function_basis=args.n_function_basis, n_coefficients=args.n_coefficients, - knn_kernel=args.k_neighbors, - knn_bandwidth=args.knn_bandwidth, - bandwidth_variability=args.bandwidth_variability, - c=args.c, + measure=measure, + function_basis=function_basis, use_mean_centres=True, + regularisation_method="diffusion", ) vals_1, forms_1 = dg.laplacian(1).spectrum() From 53d53c05dca9972b4c676f377c94de24f3ab371a Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 09:57:25 -0700 Subject: [PATCH 4/9] Align symmetric spectral eigenvector ordering with reference --- include/igneous/ops/spectral.hpp | 1 - notes/hodge/journal.md | 24 ++++++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/include/igneous/ops/spectral.hpp b/include/igneous/ops/spectral.hpp index 3f8002e..f6a0c12 100644 --- a/include/igneous/ops/spectral.hpp +++ b/include/igneous/ops/spectral.hpp @@ -337,7 +337,6 @@ void compute_eigenbasis(MeshT &mesh, int n_eigenvectors) { Eigen::MatrixXf basis = eigs.eigenvectors(nconv); basis = basis.array().colwise() * inv_sqrt_rows.array(); - basis = basis.rowwise().reverse().eval(); if (std::abs(basis(0, 0)) > 1e-12f) { basis /= basis(0, 0); } diff --git a/notes/hodge/journal.md b/notes/hodge/journal.md index 4faa9ca..d2b6d42 100644 --- a/notes/hodge/journal.md +++ b/notes/hodge/journal.md @@ -160,3 +160,27 @@ Use one entry per parity hypothesis. - Decision: `rejected` - Notes: - Severe regression; reverted immediately. + +## 2026-02-14 Hypothesis A3: Symmetric Eigenvector Ordering Parity +- Timestamp: 2026-02-14T16:57:05Z +- Commit: `3c07519` + working tree (uncommitted) +- Hypothesis: Reversing symmetric-eigensolver output columns after `LargestMagn` sorting is inconsistent with reference ordering and inflates parity error; removing the reversal should improve harmonic and circular alignment. +- Files touched: + - `include/igneous/ops/spectral.hpp` +- Commands: + - `cmake --build build -j8` + - `scripts/hodge/run_parity_round.sh` + - `scripts/hodge/run_parity_round.sh` (stability repeat) +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-162853/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-165639/report/parity_report.json` + - Stability report: `notes/hodge/results/round_20260214-165649/report/parity_report.json` + - Composite: `0.042955` (`-35.80%` vs baseline) + - Harmonic subspace max principal angle: `4.298653 deg` (improved from `6.376378 deg`) + - Harmonic Procrustes error: `0.066736` (improved from `0.097143`) + - Circular correlation min: `0.999483` (improved from `0.996571`) + - Circular wrapped P95 max: `0.059370 rad` (improved from `0.163400 rad`) +- Decision: `kept` +- Notes: + - Improvement is deterministic across repeated rounds with identical metrics. + - Final gates are not yet fully met; remaining blocker is harmonic-form error above target thresholds. From 35b254ff1b26343a430606816cecf68ccd33d6b3 Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 10:10:21 -0700 Subject: [PATCH 5/9] Add Hodge parity plot generator and hypothesis journal updates --- notes/hodge/journal.md | 52 ++++++ scripts/hodge/plot_hodge_outputs.py | 252 ++++++++++++++++++++++++++++ 2 files changed, 304 insertions(+) create mode 100755 scripts/hodge/plot_hodge_outputs.py diff --git a/notes/hodge/journal.md b/notes/hodge/journal.md index d2b6d42..917cfc5 100644 --- a/notes/hodge/journal.md +++ b/notes/hodge/journal.md @@ -184,3 +184,55 @@ Use one entry per parity hypothesis. - Notes: - Improvement is deterministic across repeated rounds with identical metrics. - Final gates are not yet fully met; remaining blocker is harmonic-form error above target thresholds. + +## 2026-02-14 Hypothesis A4: Down-Laplacian Composition via Function Gram Inverse +- Timestamp: 2026-02-14T17:02:39Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: Replacing `L_down = D D^T` with `L_down = D G0^{+} D^T` (spectral pseudo-inverse of function Gram) should match reference operator composition and reduce harmonic mismatch. +- Files touched: + - `include/igneous/ops/hodge.hpp` + - `src/main_hodge.cpp` + - `tests/test_ops_hodge.cpp` + - `benches/bench_pipelines.cpp` + - `benches/bench_dod.cpp` +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-165639/report/parity_report.json scripts/hodge/run_parity_round.sh` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-165639/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-170227/report/parity_report.json` + - Composite: `0.135324` (`+215.03%` regression) + - Harmonic subspace max principal angle: `4.298795 deg` (no meaningful change) + - Harmonic Procrustes error: `0.066738` (no meaningful change) + - Circular correlation min: `0.711944` (major regression) + - Circular wrapped P95 max: `1.541583 rad` (major regression) +- Decision: `rejected` +- Notes: + - Change severely degraded circular parity while not improving harmonic parity. + - Patch discarded and reverted. + +## 2026-02-14 Hypothesis A5: Full-NCV Symmetric Eigensolve +- Timestamp: 2026-02-14T17:06:35Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: Disabling the compact symmetric eigensolver path (`ncv = k + 16`) and always using full `ncv` parity settings should reduce harmonic mismatch. +- Files touched: + - `include/igneous/ops/spectral.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-165639/report/parity_report.json scripts/hodge/run_parity_round.sh` + - Circular mode sweep over `(mode_0, mode_1) in [0..4]x[0..4]` using `scripts/hodge/run_cpp_hodge.sh` + `scripts/hodge/compare_hodge_outputs.py` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-165639/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-170501/report/parity_report.json` + - Composite: `0.135310` (`+215.00%` regression) + - Harmonic subspace max principal angle: `4.297332 deg` (minor change) + - Harmonic Procrustes error: `0.066720` (minor change) + - Circular correlation min: `0.711942` (major regression) + - Circular wrapped P95 max: `1.541530 rad` (major regression) + - Best mode pair after sweep remained `(0,0)`; no mode remapping recovered baseline parity. +- Decision: `rejected` +- Notes: + - Regression is intrinsic to the full-`ncv` patch under this pipeline and not a mode-index artifact. + - Patch discarded and reverted. diff --git a/scripts/hodge/plot_hodge_outputs.py b/scripts/hodge/plot_hodge_outputs.py new file mode 100755 index 0000000..2773edd --- /dev/null +++ b/scripts/hodge/plot_hodge_outputs.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python3 +import argparse +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + + +def load_csv(path: Path): + if not path.exists(): + raise FileNotFoundError(f"Missing CSV: {path}") + return np.genfromtxt(path, delimiter=",", names=True) + + +def wrapped_angle_error(a: np.ndarray, b: np.ndarray) -> np.ndarray: + return np.abs(np.angle(np.exp(1j * (a - b)))) + + +def ambient_form_indices(dtype_names: tuple[str, ...]) -> list[int]: + out: list[int] = [] + i = 0 + while f"form{i}_x" in dtype_names: + if f"form{i}_y" in dtype_names and f"form{i}_z" in dtype_names: + out.append(i) + i += 1 + return out + + +def theta_indices(dtype_names: tuple[str, ...]) -> list[int]: + out: list[int] = [] + i = 0 + while f"theta_{i}" in dtype_names: + out.append(i) + i += 1 + return out + + +def sample_indices(n: int, max_points: int) -> np.ndarray: + if n <= max_points: + return np.arange(n, dtype=np.int64) + step = max(1, n // max_points) + return np.arange(0, n, step, dtype=np.int64)[:max_points] + + +def set_equal_axes(ax, x: np.ndarray, y: np.ndarray, z: np.ndarray): + mins = np.array([x.min(), y.min(), z.min()], dtype=np.float64) + maxs = np.array([x.max(), y.max(), z.max()], dtype=np.float64) + center = 0.5 * (mins + maxs) + radius = 0.5 * np.max(maxs - mins) + ax.set_xlim(center[0] - radius, center[0] + radius) + ax.set_ylim(center[1] - radius, center[1] + radius) + ax.set_zlim(center[2] - radius, center[2] + radius) + + +def plot_vector_pair( + out_path: Path, + points_ref: np.ndarray, + vec_ref: np.ndarray, + points_cpp: np.ndarray, + vec_cpp: np.ndarray, + title: str, +): + fig = plt.figure(figsize=(14, 6)) + ax_ref = fig.add_subplot(1, 2, 1, projection="3d") + ax_cpp = fig.add_subplot(1, 2, 2, projection="3d") + + for ax, pts, vecs, label in [ + (ax_ref, points_ref, vec_ref, "Reference"), + (ax_cpp, points_cpp, vec_cpp, "C++"), + ]: + ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], s=4, c="lightgray", alpha=0.35) + ax.quiver( + pts[:, 0], + pts[:, 1], + pts[:, 2], + vecs[:, 0], + vecs[:, 1], + vecs[:, 2], + length=0.18, + normalize=True, + linewidth=0.5, + color="#005f73", + ) + set_equal_axes(ax, pts[:, 0], pts[:, 1], pts[:, 2]) + ax.set_title(label) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_zticks([]) + + fig.suptitle(title) + fig.tight_layout() + fig.savefig(out_path, dpi=200) + plt.close(fig) + + +def plot_angle_pair( + out_path: Path, + points_ref: np.ndarray, + theta_ref: np.ndarray, + points_cpp: np.ndarray, + theta_cpp: np.ndarray, + title: str, +): + fig = plt.figure(figsize=(14, 6)) + ax_ref = fig.add_subplot(1, 2, 1, projection="3d") + ax_cpp = fig.add_subplot(1, 2, 2, projection="3d") + + for ax, pts, theta, label in [ + (ax_ref, points_ref, theta_ref, "Reference"), + (ax_cpp, points_cpp, theta_cpp, "C++"), + ]: + sc = ax.scatter( + pts[:, 0], + pts[:, 1], + pts[:, 2], + c=theta, + s=9, + cmap="twilight", + vmin=0.0, + vmax=2.0 * np.pi, + ) + set_equal_axes(ax, pts[:, 0], pts[:, 1], pts[:, 2]) + ax.set_title(label) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_zticks([]) + plt.colorbar(sc, ax=ax, shrink=0.65, pad=0.03) + + fig.suptitle(title) + fig.tight_layout() + fig.savefig(out_path, dpi=200) + plt.close(fig) + + +def plot_angle_error_hist( + out_path: Path, errors: np.ndarray, title: str, p95: float, mean: float +): + fig, ax = plt.subplots(figsize=(8, 5)) + ax.hist(errors, bins=48, color="#0a9396", alpha=0.9) + ax.axvline(p95, color="#ae2012", linestyle="--", linewidth=2, label=f"P95={p95:.3f}") + ax.axvline(mean, color="#005f73", linestyle="-.", linewidth=2, label=f"Mean={mean:.3f}") + ax.set_xlabel("Wrapped Angular Error (rad)") + ax.set_ylabel("Count") + ax.set_title(title) + ax.legend() + fig.tight_layout() + fig.savefig(out_path, dpi=200) + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser(description="Plot harmonic 1-forms and circular coordinates.") + parser.add_argument("--round-dir", required=True, help="Parity round directory") + parser.add_argument( + "--output-dir", + default=None, + help="Optional output directory (default: /report/plots)", + ) + parser.add_argument( + "--max-points", + type=int, + default=350, + help="Maximum points used for vector-field quiver plots", + ) + args = parser.parse_args() + + round_dir = Path(args.round_dir).resolve() + output_dir = ( + Path(args.output_dir).resolve() + if args.output_dir is not None + else round_dir / "report" / "plots" + ) + output_dir.mkdir(parents=True, exist_ok=True) + + ref_ambient = load_csv(round_dir / "reference" / "reference_harmonic_ambient.csv") + cpp_ambient = load_csv(round_dir / "cpp" / "harmonic_ambient.csv") + ref_circular = load_csv(round_dir / "reference" / "reference_circular_coordinates.csv") + cpp_circular = load_csv(round_dir / "cpp" / "circular_coordinates.csv") + + n = min(ref_ambient.shape[0], cpp_ambient.shape[0], ref_circular.shape[0], cpp_circular.shape[0]) + ref_ambient = ref_ambient[:n] + cpp_ambient = cpp_ambient[:n] + ref_circular = ref_circular[:n] + cpp_circular = cpp_circular[:n] + + points_ref = np.column_stack([ref_ambient["x"], ref_ambient["y"], ref_ambient["z"]]) + points_cpp = np.column_stack([cpp_ambient["x"], cpp_ambient["y"], cpp_ambient["z"]]) + + vec_idx = sample_indices(n, args.max_points) + + ref_forms = ambient_form_indices(ref_ambient.dtype.names) + cpp_forms = ambient_form_indices(cpp_ambient.dtype.names) + common_forms = sorted(set(ref_forms).intersection(cpp_forms)) + + generated_paths: list[Path] = [] + for form in common_forms[:2]: + ref_vec = np.column_stack( + [ref_ambient[f"form{form}_x"], ref_ambient[f"form{form}_y"], ref_ambient[f"form{form}_z"]] + ) + cpp_vec = np.column_stack( + [cpp_ambient[f"form{form}_x"], cpp_ambient[f"form{form}_y"], cpp_ambient[f"form{form}_z"]] + ) + out_path = output_dir / f"harmonic_form_{form}_vector_field.png" + plot_vector_pair( + out_path=out_path, + points_ref=points_ref[vec_idx], + vec_ref=ref_vec[vec_idx], + points_cpp=points_cpp[vec_idx], + vec_cpp=cpp_vec[vec_idx], + title=f"Harmonic 1-Form {form}: Vector Field", + ) + generated_paths.append(out_path) + + ref_thetas = theta_indices(ref_circular.dtype.names) + cpp_thetas = theta_indices(cpp_circular.dtype.names) + common_thetas = sorted(set(ref_thetas).intersection(cpp_thetas)) + + for theta_idx in common_thetas[:2]: + ref_theta = np.asarray(ref_circular[f"theta_{theta_idx}"], dtype=np.float64) + cpp_theta = np.asarray(cpp_circular[f"theta_{theta_idx}"], dtype=np.float64) + pair_path = output_dir / f"circular_theta_{theta_idx}_compare.png" + plot_angle_pair( + out_path=pair_path, + points_ref=points_ref, + theta_ref=ref_theta, + points_cpp=points_cpp, + theta_cpp=cpp_theta, + title=f"Circular Coordinate theta_{theta_idx}", + ) + generated_paths.append(pair_path) + + err = wrapped_angle_error(ref_theta, cpp_theta) + hist_path = output_dir / f"circular_theta_{theta_idx}_error_hist.png" + plot_angle_error_hist( + out_path=hist_path, + errors=err, + title=f"theta_{theta_idx}: Wrapped Angular Error", + p95=float(np.percentile(err, 95.0)), + mean=float(np.mean(err)), + ) + generated_paths.append(hist_path) + + print(f"plots_dir={output_dir}") + for p in generated_paths: + print(f"plot={p}") + + +if __name__ == "__main__": + main() From 2b1cf39145958a39f562afb77773b4958d9ba02c Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 10:10:45 -0700 Subject: [PATCH 6/9] Install matplotlib in Hodge reference venv setup --- scripts/hodge/setup_reference_env.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/hodge/setup_reference_env.sh b/scripts/hodge/setup_reference_env.sh index 541a451..26eef6d 100755 --- a/scripts/hodge/setup_reference_env.sh +++ b/scripts/hodge/setup_reference_env.sh @@ -24,7 +24,7 @@ fi source "${VENV_PATH}/bin/activate" python -m pip -q install --upgrade pip -python -m pip -q install numpy scipy scikit-learn opt_einsum +python -m pip -q install numpy scipy scikit-learn opt_einsum matplotlib python -m pip -q install -e "${ROOT_DIR}/DiffusionGeometry" echo "${VENV_PATH}" From 0c67e541dd55dd6ba525acb1b7bb929bfbba907c Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 10:42:13 -0700 Subject: [PATCH 7/9] Align reference harness immersion and default circular modes --- notes/hodge/journal.md | 153 +++++++++++++++++++++++++++ scripts/hodge/plot_hodge_outputs.py | 62 +++++++++++ scripts/hodge/run_cpp_hodge.sh | 2 +- scripts/hodge/run_parity_round.sh | 2 +- scripts/hodge/run_reference_hodge.py | 4 +- src/main_hodge.cpp | 2 +- 6 files changed, 220 insertions(+), 5 deletions(-) diff --git a/notes/hodge/journal.md b/notes/hodge/journal.md index 917cfc5..e88bae6 100644 --- a/notes/hodge/journal.md +++ b/notes/hodge/journal.md @@ -236,3 +236,156 @@ Use one entry per parity hypothesis. - Notes: - Regression is intrinsic to the full-`ncv` patch under this pipeline and not a mode-index artifact. - Patch discarded and reverted. + +## 2026-02-14 Hypothesis A6: NumPy-Exact Even-Length Median for Bandwidth Normalisation +- Timestamp: 2026-02-14T17:28:02Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: Matching NumPy median behavior for even-length arrays (`0.5*(x[n/2-1]+x[n/2])`) in `bandwidths_B` normalisation will reduce diffusion-kernel drift and improve harmonic parity. +- Files touched: + - `include/igneous/data/topology.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-170725/report/parity_report.json scripts/hodge/run_parity_round.sh` + - Circular mode sweep over `(mode_0, mode_1) in [0..3]x[0..3]` using `scripts/hodge/run_cpp_hodge.sh` + `scripts/hodge/compare_hodge_outputs.py` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-170725/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-172638/report/parity_report.json` + - Composite: `0.141622` (`+229.71%` regression) + - Harmonic subspace max principal angle: `4.297839 deg` (no meaningful change) + - Harmonic Procrustes error: `0.066728` (no meaningful change) + - Circular correlation min: `0.711971` (major regression) + - Circular wrapped P95 max: `1.541341 rad` (major regression) + - Mode sweep confirmed no recovery path; best remained `(0,0)`. +- Decision: `rejected` +- Notes: + - This median change materially perturbs circular operator spectrum without improving harmonic metrics. + - Patch discarded and reverted. + +## 2026-02-14 Hypothesis H3: Double Precision in Hodge Spectrum Solve Only +- Timestamp: 2026-02-14T17:29:37Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: Promoting only the generalized 1-form eigensolve (`compute_hodge_spectrum`) to `double` precision while keeping circular solve unchanged will improve harmonic parity without circular regression. +- Files touched: + - `include/igneous/ops/hodge.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-170725/report/parity_report.json scripts/hodge/run_parity_round.sh` + - Circular mode sweep over `(mode_0, mode_1) in [0..3]x[0..3]` using `scripts/hodge/run_cpp_hodge.sh` + `scripts/hodge/compare_hodge_outputs.py` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-170725/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-172859/report/parity_report.json` + - Composite: `0.076155` (`+77.29%` regression) + - Harmonic subspace max principal angle: `4.298632 deg` (no meaningful change) + - Harmonic Procrustes error: `0.066736` (no meaningful change) + - Circular correlation min: `0.915266` (regression) + - Circular wrapped P95 max: `0.622422 rad` (regression) + - Mode sweep confirmed no recovery path; best remained `(0,0)`. +- Decision: `rejected` +- Notes: + - Precision promotion in the Hodge spectrum alone destabilized circular parity with no harmonic gain. + - Patch discarded and reverted. + +## 2026-02-14 Hypothesis C1: Reference Harness Immersion Regularisation Parity +- Timestamp: 2026-02-14T17:37:38Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: The reference harness currently bypasses diffusion regularisation by forcing `immersion_coords=points`; switching to `immersion_coords=None` should align with reference constructor semantics and reduce harmonic mismatch. +- Files touched: + - `scripts/hodge/run_reference_hodge.py` +- Commands: + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-170725/report/parity_report.json scripts/hodge/run_parity_round.sh` + - `MODE_INDICES=0,1 scripts/hodge/run_parity_round.sh` + - Circular mode sweeps on corrected reference artifacts (`mode_0,mode_1` in `[0..5]` and wide `mode_1` sweep) +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-170725/report/parity_report.json` + - Corrected reference report (`MODE_INDICES=0,0`): `notes/hodge/results/round_20260214-173113/report/parity_report.json` + - Composite: `0.068078` (`+58.49%` vs baseline) + - Harmonic subspace max principal angle: `2.261020 deg` (improved from `4.298653`) + - Harmonic Procrustes error: `0.030792` (improved from `0.066736`) + - Circular correlation min: `0.855133` + - Circular wrapped P95 max: `0.788353 rad` + - Notebook-style mode reference (`MODE_INDICES=0,1`): `notes/hodge/results/round_20260214-173327/report/parity_report.json` + - Composite: `0.044602` + - Harmonic subspace max principal angle: `2.261020 deg` + - Harmonic Procrustes error: `0.030792` + - Circular correlation min: `0.966479` + - Circular wrapped P95 max: `0.471590 rad` +- Decision: `kept` +- Notes: + - This is a concrete reference-harness correctness fix: diffusion regularisation is now active for immersion coordinates. + - Harmonic parity now meets final gates; remaining blocker is circular wrapped-error tail. + +## 2026-02-14 Hypothesis D2: Circular Eigensolve via Gram-Orthonormal Restriction +- Timestamp: 2026-02-14T17:37:38Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: Replacing direct generalized eigensolve with the reference-style Gram-orthonormal restricted eigensolve (`rcond` cutoff + standard eig) will improve circular parity. +- Files touched: + - `include/igneous/ops/hodge.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `scripts/hodge/run_parity_round.sh` + - `MODE_INDICES=0,1 scripts/hodge/run_parity_round.sh` +- Numeric parity deltas: + - Corrected-reference report (`MODE_INDICES=0,0`): `notes/hodge/results/round_20260214-173647/report/parity_report.json` + - Composite: `0.068077` (no meaningful change) + - Circular correlation min: `0.855137` (no meaningful change) + - Circular wrapped P95 max: `0.788348 rad` (no meaningful change) + - Notebook-style mode report (`MODE_INDICES=0,1`): `notes/hodge/results/round_20260214-173656/report/parity_report.json` + - Composite: `0.044602` (no meaningful change) + - Circular correlation min: `0.966479` (no meaningful change) + - Circular wrapped P95 max: `0.471592 rad` (no meaningful change) +- Decision: `rejected` +- Notes: + - Functional behavior is effectively unchanged from the existing generalized formulation on this problem. + - Patch discarded and reverted. + +## 2026-02-14 Hypothesis D3: Double-Precision Circular Operator Eigensolve +- Timestamp: 2026-02-14T17:39:55Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: Promoting circular operator generalized eigensolve and angle reconstruction to `double` precision should reduce wrapped-angle tail error. +- Files touched: + - `include/igneous/ops/hodge.hpp` (reverted) +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-173327/report/parity_report.json MODE_INDICES=0,1 scripts/hodge/run_parity_round.sh` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-173327/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-173927/report/parity_report.json` + - Composite: `0.044602` (no meaningful change) + - Harmonic subspace max principal angle: `2.261020 deg` (unchanged) + - Harmonic Procrustes error: `0.030792` (unchanged) + - Circular correlation min: `0.966479` (unchanged) + - Circular wrapped P95 max: `0.471593 rad` (unchanged) +- Decision: `rejected` +- Notes: + - Precision upgrade did not materially move any parity metric. + - Patch discarded and reverted. + +## 2026-02-14 Hypothesis D4: Notebook-Style Circular Mode Defaults (0,1) +- Timestamp: 2026-02-14T17:41:20Z +- Commit: `53d53c0` + working tree (uncommitted) +- Hypothesis: After correcting immersion-regularisation parity, setting default circular mode indices to `(0,1)` (as in reference TDA torus workflow) should materially improve circular parity versus `(0,0)` while preserving harmonic gains. +- Files touched: + - `src/main_hodge.cpp` + - `scripts/hodge/run_parity_round.sh` + - `scripts/hodge/run_cpp_hodge.sh` + - `scripts/hodge/run_reference_hodge.py` +- Commands: + - `cmake --build build -j8` + - `ctest --test-dir build --output-on-failure` + - `PREVIOUS_REPORT_JSON=notes/hodge/results/round_20260214-173113/report/parity_report.json scripts/hodge/run_parity_round.sh` +- Numeric parity deltas: + - Baseline (corrected reference, `(0,0)`): `notes/hodge/results/round_20260214-173113/report/parity_report.json` + - Report (`(0,1)` defaults): `notes/hodge/results/round_20260214-174054/report/parity_report.json` + - Composite: `0.044602` (`-34.49%` improvement vs corrected `(0,0)` baseline) + - Harmonic subspace max principal angle: `2.261020 deg` (unchanged, passes target) + - Harmonic Procrustes error: `0.030792` (unchanged, passes target) + - Circular correlation min: `0.966479` (improved from `0.855133`) + - Circular wrapped P95 max: `0.471590 rad` (improved from `0.788353 rad`) +- Decision: `kept` +- Notes: + - This is a structural mode-selection alignment with the reference notebook workflow. + - Circular tail error remains above final gate threshold; further operator-level alignment is still required. diff --git a/scripts/hodge/plot_hodge_outputs.py b/scripts/hodge/plot_hodge_outputs.py index 2773edd..6bf7674 100755 --- a/scripts/hodge/plot_hodge_outputs.py +++ b/scripts/hodge/plot_hodge_outputs.py @@ -135,6 +135,53 @@ def plot_angle_pair( plt.close(fig) +def plot_circular_forms_grid( + out_path: Path, + points_ref: np.ndarray, + points_cpp: np.ndarray, + theta_ref0: np.ndarray, + theta_cpp0: np.ndarray, + theta_ref1: np.ndarray, + theta_cpp1: np.ndarray, +): + fig = plt.figure(figsize=(14, 12)) + axes = [ + fig.add_subplot(2, 2, 1, projection="3d"), + fig.add_subplot(2, 2, 2, projection="3d"), + fig.add_subplot(2, 2, 3, projection="3d"), + fig.add_subplot(2, 2, 4, projection="3d"), + ] + panels = [ + (axes[0], points_ref, theta_ref0, "Reference: Circular from Harmonic 1-Form 0"), + (axes[1], points_cpp, theta_cpp0, "C++: Circular from Harmonic 1-Form 0"), + (axes[2], points_ref, theta_ref1, "Reference: Circular from Harmonic 1-Form 1"), + (axes[3], points_cpp, theta_cpp1, "C++: Circular from Harmonic 1-Form 1"), + ] + + for ax, pts, theta, label in panels: + sc = ax.scatter( + pts[:, 0], + pts[:, 1], + pts[:, 2], + c=theta, + s=8, + cmap="twilight", + vmin=0.0, + vmax=2.0 * np.pi, + ) + set_equal_axes(ax, pts[:, 0], pts[:, 1], pts[:, 2]) + ax.set_title(label, fontsize=10) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_zticks([]) + plt.colorbar(sc, ax=ax, shrink=0.58, pad=0.02) + + fig.suptitle("Circular Coordinates by Source Harmonic 1-Form", fontsize=14) + fig.tight_layout() + fig.savefig(out_path, dpi=220) + plt.close(fig) + + def plot_angle_error_hist( out_path: Path, errors: np.ndarray, title: str, p95: float, mean: float ): @@ -243,6 +290,21 @@ def main(): ) generated_paths.append(hist_path) + if len(common_thetas) >= 2: + theta0 = common_thetas[0] + theta1 = common_thetas[1] + grid_path = output_dir / "circular_from_harmonic_forms_compare.png" + plot_circular_forms_grid( + out_path=grid_path, + points_ref=points_ref, + points_cpp=points_cpp, + theta_ref0=np.asarray(ref_circular[f"theta_{theta0}"], dtype=np.float64), + theta_cpp0=np.asarray(cpp_circular[f"theta_{theta0}"], dtype=np.float64), + theta_ref1=np.asarray(ref_circular[f"theta_{theta1}"], dtype=np.float64), + theta_cpp1=np.asarray(cpp_circular[f"theta_{theta1}"], dtype=np.float64), + ) + generated_paths.append(grid_path) + print(f"plots_dir={output_dir}") for p in generated_paths: print(f"plot={p}") diff --git a/scripts/hodge/run_cpp_hodge.sh b/scripts/hodge/run_cpp_hodge.sh index 7d30a14..cdba95f 100755 --- a/scripts/hodge/run_cpp_hodge.sh +++ b/scripts/hodge/run_cpp_hodge.sh @@ -17,7 +17,7 @@ BANDWIDTH_VARIABILITY="${BANDWIDTH_VARIABILITY:--0.5}" C_PARAM="${C_PARAM:-0.0}" CIRCULAR_LAMBDA="${CIRCULAR_LAMBDA:-1.0}" CIRCULAR_MODE_0="${CIRCULAR_MODE_0:-0}" -CIRCULAR_MODE_1="${CIRCULAR_MODE_1:-0}" +CIRCULAR_MODE_1="${CIRCULAR_MODE_1:-1}" if [[ ! -x "${EXE}" ]]; then echo "Executable not found: ${EXE}" >&2 diff --git a/scripts/hodge/run_parity_round.sh b/scripts/hodge/run_parity_round.sh index 247b7c9..b3338da 100755 --- a/scripts/hodge/run_parity_round.sh +++ b/scripts/hodge/run_parity_round.sh @@ -42,7 +42,7 @@ C_PARAM="${C_PARAM:-0.0}" CIRCULAR_LAMBDA="${CIRCULAR_LAMBDA:-1.0}" FORM_INDICES="${FORM_INDICES:-0,1}" -MODE_INDICES="${MODE_INDICES:-0,0}" +MODE_INDICES="${MODE_INDICES:-0,1}" INPUT_CSV="${INPUT_DIR}/torus.csv" "${PYTHON_BIN}" "${SCRIPT_DIR}/generate_input_torus.py" \ diff --git a/scripts/hodge/run_reference_hodge.py b/scripts/hodge/run_reference_hodge.py index 1df8249..9210a6e 100755 --- a/scripts/hodge/run_reference_hodge.py +++ b/scripts/hodge/run_reference_hodge.py @@ -97,7 +97,7 @@ def main() -> None: parser.add_argument("--c", type=float, default=0.0) parser.add_argument("--circular-lambda", type=float, default=1.0) parser.add_argument("--form-indices", default="0,1") - parser.add_argument("--mode-indices", default="0,0") + parser.add_argument("--mode-indices", default="0,1") args = parser.parse_args() out_dir = pathlib.Path(args.output_dir) @@ -121,7 +121,7 @@ def main() -> None: dg = DiffusionGeometry.from_knn_kernel( nbr_indices=nbr_indices, kernel=kernel, - immersion_coords=points, + immersion_coords=None, data_matrix=points, bandwidths=bandwidths, n_function_basis=args.n_function_basis, diff --git a/src/main_hodge.cpp b/src/main_hodge.cpp index 9509f4d..e16e717 100644 --- a/src/main_hodge.cpp +++ b/src/main_hodge.cpp @@ -36,7 +36,7 @@ struct Config { float circular_lambda = 1.0f; int circular_mode_0 = 0; - int circular_mode_1 = 0; + int circular_mode_1 = 1; bool export_ply = true; }; From 78ee203e5d4e58773d601bb25bd3177f8c14dc57 Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 11:31:45 -0700 Subject: [PATCH 8/9] Fix circular parity via harmonic form sign canonicalization --- include/igneous/ops/hodge.hpp | 21 +- notes/hodge/journal.md | 53 ++++ .../hodge/compare_circular_operator_debug.py | 243 ++++++++++++++++++ scripts/hodge/run_reference_hodge.py | 44 +++- src/main_hodge.cpp | 117 ++++++++- 5 files changed, 465 insertions(+), 13 deletions(-) create mode 100755 scripts/hodge/compare_circular_operator_debug.py diff --git a/include/igneous/ops/hodge.hpp b/include/igneous/ops/hodge.hpp index b125b2f..979566d 100644 --- a/include/igneous/ops/hodge.hpp +++ b/include/igneous/ops/hodge.hpp @@ -201,13 +201,23 @@ compute_hodge_spectrum(const Eigen::MatrixXf &laplacian, return std::make_pair(evals, evecs); } +template +struct CircularOperatorDiagnostics { + Eigen::MatrixXf x_weak; + Eigen::MatrixXf laplacian0_weak; + Eigen::MatrixXf function_gram; + Eigen::VectorXcf sorted_evals; + std::vector positive_imag_indices; +}; + template Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, const Eigen::VectorXf &alpha_coeffs, float bandwidth, float lambda = 1.0f, int positive_imag_mode = 0, - std::complex *selected_eval = nullptr) { + std::complex *selected_eval = nullptr, + CircularOperatorDiagnostics *diagnostics = nullptr) { const auto &U = mesh.topology.eigen_basis; const auto &mu = mesh.topology.mu; @@ -266,6 +276,11 @@ Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, Eigen::MatrixXf function_gram = U.transpose() * (U.array().colwise() * mu.array()).matrix(); Eigen::MatrixXf operator_weak = X_op - (lambda * laplacian0_weak); + if (diagnostics != nullptr) { + diagnostics->x_weak = X_op; + diagnostics->laplacian0_weak = laplacian0_weak; + diagnostics->function_gram = function_gram; + } Eigen::GeneralizedEigenSolver solver(operator_weak, function_gram); if (solver.info() != Eigen::Success) { @@ -307,6 +322,10 @@ Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, positive_imag_indices.push_back(i); } } + if (diagnostics != nullptr) { + diagnostics->sorted_evals = sorted_evals; + diagnostics->positive_imag_indices = positive_imag_indices; + } if (positive_imag_indices.empty()) { return Eigen::VectorXf::Zero(static_cast(n_verts)); diff --git a/notes/hodge/journal.md b/notes/hodge/journal.md index e88bae6..8348bc7 100644 --- a/notes/hodge/journal.md +++ b/notes/hodge/journal.md @@ -389,3 +389,56 @@ Use one entry per parity hypothesis. - Notes: - This is a structural mode-selection alignment with the reference notebook workflow. - Circular tail error remains above final gate threshold; further operator-level alignment is still required. + +## 2026-02-14 Hypothesis D5: Basis-Aligned Circular-Operator Diagnostics +- Timestamp: 2026-02-14T18:28:05Z +- Commit: `0c67e54` + working tree (uncommitted) +- Hypothesis: Circular matrix mismatches may be dominated by basis-coordinate gauge differences; compare operators after explicit function-basis alignment (`T^T W T`) to isolate true structural drift. +- Files touched: + - `include/igneous/ops/hodge.hpp` + - `src/main_hodge.cpp` + - `scripts/hodge/run_reference_hodge.py` + - `scripts/hodge/compare_circular_operator_debug.py` +- Commands: + - `cmake --build build --target igneous-hodge -j8` + - `scripts/hodge/run_parity_round.sh` + - `python3 scripts/hodge/compare_circular_operator_debug.py --reference-dir notes/hodge/results/round_20260214-182805/reference --cpp-dir notes/hodge/results/round_20260214-182805/cpp --output-markdown notes/hodge/results/round_20260214-182805/report/circular_operator_debug.md --output-json notes/hodge/results/round_20260214-182805/report/circular_operator_debug.json` +- Numeric parity deltas: + - Parity report: `notes/hodge/results/round_20260214-182805/report/parity_report.json` + - Composite: `0.044602` (unchanged) + - Circular correlation min: `0.966479` + - Circular wrapped P95 max: `0.471590 rad` + - Basis-aligned debug report: `notes/hodge/results/round_20260214-182805/report/circular_operator_debug.json` + - `form0_x_weak` aligned rel-Fro: `0.054460` + - `form1_x_weak` aligned rel-Fro: `1.995615` + - `form0_operator_weak` aligned rel-Fro: `0.036106` + - `form1_operator_weak` aligned rel-Fro: `0.469989` +- Decision: `kept` +- Notes: + - Diagnostic confirmed the remaining circular error is structural and form-specific, not a global basis-coordinate artifact. + - `form1_x_weak` being ~`2.0` after basis alignment indicates near sign-negation versus reference. + +## 2026-02-14 Hypothesis D6: Harmonic 1-Form Sign Canonicalization for Circular Coordinates +- Timestamp: 2026-02-14T18:30:22Z +- Commit: `0c67e54` + working tree (uncommitted) +- Hypothesis: Circular mismatch for harmonic form 1 is driven by eigenvector sign gauge; canonicalizing harmonic form signs via deterministic ambient orientation score before circular solve should collapse the circular tail error. +- Files touched: + - `src/main_hodge.cpp` +- Commands: + - `cmake --build build --target igneous-hodge -j8` + - `scripts/hodge/run_parity_round.sh` + - `python3 scripts/hodge/compare_circular_operator_debug.py --reference-dir notes/hodge/results/round_20260214-183022/reference --cpp-dir notes/hodge/results/round_20260214-183022/cpp --output-markdown notes/hodge/results/round_20260214-183022/report/circular_operator_debug.md --output-json notes/hodge/results/round_20260214-183022/report/circular_operator_debug.json` + - `notes/hodge/.venv_ref/bin/python scripts/hodge/plot_hodge_outputs.py --round-dir notes/hodge/results/round_20260214-183022 --output-dir notes/hodge/results/round_20260214-183022/report/plots` +- Numeric parity deltas: + - Baseline report: `notes/hodge/results/round_20260214-182805/report/parity_report.json` + - Report: `notes/hodge/results/round_20260214-183022/report/parity_report.json` + - Composite: `0.021085` (`-52.73%` improvement) + - Harmonic subspace max principal angle: `2.261020 deg` (unchanged, pass) + - Harmonic Procrustes error: `0.030792` (unchanged, pass) + - Circular correlation min: `0.999843` (improved from `0.966479`, pass) + - Circular wrapped P95 max: `0.031461 rad` (improved from `0.471590`, pass) + - Final pass gate: `true` +- Decision: `kept` +- Notes: + - Circular parity now passes all target thresholds on the canonical torus scenario. + - Basis-aligned operator debug also converged for form 1 (`form1_x_weak` aligned rel-Fro `0.062090`, `form1_operator_weak` `0.037623`). diff --git a/scripts/hodge/compare_circular_operator_debug.py b/scripts/hodge/compare_circular_operator_debug.py new file mode 100755 index 0000000..181da31 --- /dev/null +++ b/scripts/hodge/compare_circular_operator_debug.py @@ -0,0 +1,243 @@ +#!/usr/bin/env python3 +import argparse +import json +from pathlib import Path +from typing import Any + +import numpy as np + + +def load_columns(path: Path) -> dict[str, np.ndarray]: + arr = np.genfromtxt(path, delimiter=",", names=True) + if arr.shape == (): + arr = np.array([arr], dtype=arr.dtype) + return {name: np.asarray(arr[name]) for name in arr.dtype.names} + + +def load_matrix(path: Path) -> np.ndarray: + cols = load_columns(path) + keys = [k for k in cols.keys() if k.startswith("c")] + keys.sort(key=lambda x: int(x[1:])) + return np.column_stack([cols[k] for k in keys]).astype(np.float64) + + +def load_evals(path: Path) -> np.ndarray: + cols = load_columns(path) + return cols["real"].astype(np.float64) + 1j * cols["imag"].astype(np.float64) + + +def load_positive_evals(path: Path) -> np.ndarray: + cols = load_columns(path) + evals = cols["real"].astype(np.float64) + 1j * cols["imag"].astype(np.float64) + mask = cols["positive_imag"].astype(np.float64) > 0.5 + return evals[mask] + + +def rel_fro(a: np.ndarray, b: np.ndarray) -> float: + n = min(a.shape[0], b.shape[0]) + m = min(a.shape[1], b.shape[1]) + aa = a[:n, :m] + bb = b[:n, :m] + denom = max(np.linalg.norm(aa, ord="fro"), 1e-12) + return float(np.linalg.norm(aa - bb, ord="fro") / denom) + + +def align_basis( + reference_basis: np.ndarray, cpp_basis: np.ndarray +) -> tuple[np.ndarray, dict[str, float]]: + n = min(reference_basis.shape[0], cpp_basis.shape[0]) + k_ref = reference_basis.shape[1] + k_cpp = cpp_basis.shape[1] + + ref = reference_basis[:n, :k_ref] + cpp = cpp_basis[:n, :k_cpp] + + transform, _, rank, singular_values = np.linalg.lstsq(cpp, ref, rcond=None) + reconstructed = cpp @ transform + denom = max(np.linalg.norm(ref, ord="fro"), 1e-12) + recon_rel = float(np.linalg.norm(reconstructed - ref, ord="fro") / denom) + + cond = float("inf") + if singular_values.size > 0: + cond = float(singular_values[0] / max(singular_values[-1], 1e-12)) + + stats = { + "rows_used": float(n), + "reference_cols": float(k_ref), + "cpp_cols": float(k_cpp), + "transform_rank": float(rank), + "cpp_basis_condition_estimate": cond, + "reconstruction_rel_fro": recon_rel, + } + return transform, stats + + +def transform_weak_to_reference_basis(cpp_weak: np.ndarray, transform: np.ndarray) -> np.ndarray: + return transform.conj().T @ cpp_weak @ transform + + +def eval_alignment_error(ref: np.ndarray, cpp: np.ndarray) -> dict[str, float]: + k = min(ref.shape[0], cpp.shape[0]) + if k == 0: + return {"mean_abs_diff": float("nan"), "p95_abs_diff": float("nan")} + d = np.abs(ref[:k] - cpp[:k]) + return { + "mean_abs_diff": float(np.mean(d)), + "p95_abs_diff": float(np.percentile(d, 95.0)), + } + + +def render_markdown(payload: dict[str, Any]) -> str: + lines: list[str] = [] + lines.append("# Circular Operator Debug Comparison") + lines.append("") + lines.append(f"- Reference dir: `{payload['reference_dir']}`") + lines.append(f"- C++ dir: `{payload['cpp_dir']}`") + lines.append("") + lines.append("## Matrix Relative Frobenius Errors") + lines.append("") + for k, v in payload["matrix_rel_fro"].items(): + lines.append(f"- {k}: `{v:.6e}`") + lines.append("") + if payload.get("matrix_rel_fro_aligned"): + lines.append("## Matrix Relative Frobenius Errors (Basis Aligned)") + lines.append("") + for k, v in payload["matrix_rel_fro_aligned"].items(): + lines.append(f"- {k}: `{v:.6e}`") + lines.append("") + if payload.get("basis_alignment"): + lines.append("## Basis Alignment") + lines.append("") + stats = payload["basis_alignment"] + lines.append(f"- rows used: `{int(stats['rows_used'])}`") + lines.append(f"- reference columns: `{int(stats['reference_cols'])}`") + lines.append(f"- cpp columns: `{int(stats['cpp_cols'])}`") + lines.append(f"- transform rank: `{int(stats['transform_rank'])}`") + lines.append( + f"- cpp basis condition estimate: `{stats['cpp_basis_condition_estimate']:.6e}`" + ) + lines.append( + f"- reconstruction relative Frobenius error: `{stats['reconstruction_rel_fro']:.6e}`" + ) + lines.append("") + + lines.append("## Positive-Imag Eigenvalue Alignment") + lines.append("") + for form, stats in payload["eval_alignment"].items(): + lines.append( + f"- {form}: mean abs diff `{stats['mean_abs_diff']:.6e}`, " + f"P95 abs diff `{stats['p95_abs_diff']:.6e}`" + ) + lines.append("") + return "\n".join(lines) + + +def main() -> None: + parser = argparse.ArgumentParser( + description="Compare circular operator debug artifacts between reference and C++ outputs." + ) + parser.add_argument("--reference-dir", required=True) + parser.add_argument("--cpp-dir", required=True) + parser.add_argument("--output-markdown", required=True) + parser.add_argument("--output-json", required=True) + args = parser.parse_args() + + ref_dir = Path(args.reference_dir) + cpp_dir = Path(args.cpp_dir) + + matrix_pairs = { + "function_gram": ( + ref_dir / "reference_function_gram.csv", + cpp_dir / "function_gram.csv", + ), + "laplacian0_weak": ( + ref_dir / "reference_laplacian0_weak.csv", + cpp_dir / "laplacian0_weak.csv", + ), + "form0_x_weak": ( + ref_dir / "reference_circular_operator_form0_x_weak.csv", + cpp_dir / "circular_operator_form0_x_weak.csv", + ), + "form1_x_weak": ( + ref_dir / "reference_circular_operator_form1_x_weak.csv", + cpp_dir / "circular_operator_form1_x_weak.csv", + ), + "form0_operator_weak": ( + ref_dir / "reference_circular_operator_form0_operator_weak.csv", + cpp_dir / "circular_operator_form0_operator_weak.csv", + ), + "form1_operator_weak": ( + ref_dir / "reference_circular_operator_form1_operator_weak.csv", + cpp_dir / "circular_operator_form1_operator_weak.csv", + ), + } + + matrix_rel_fro: dict[str, float] = {} + for name, (rp, cp) in matrix_pairs.items(): + if not rp.exists() or not cp.exists(): + matrix_rel_fro[name] = float("nan") + continue + matrix_rel_fro[name] = rel_fro(load_matrix(rp), load_matrix(cp)) + + matrix_rel_fro_aligned: dict[str, float] = {} + basis_alignment: dict[str, float] = {} + ref_basis_path = ref_dir / "reference_function_basis.csv" + cpp_basis_path = cpp_dir / "function_basis.csv" + if ref_basis_path.exists() and cpp_basis_path.exists(): + ref_basis = load_matrix(ref_basis_path) + cpp_basis = load_matrix(cpp_basis_path) + transform, basis_alignment = align_basis(ref_basis, cpp_basis) + + for name, (rp, cp) in matrix_pairs.items(): + if not rp.exists() or not cp.exists(): + matrix_rel_fro_aligned[name] = float("nan") + continue + ref_m = load_matrix(rp) + cpp_m = load_matrix(cp) + rdim = min(ref_m.shape[0], ref_m.shape[1], transform.shape[1]) + cdim = min(cpp_m.shape[0], cpp_m.shape[1], transform.shape[0]) + if rdim <= 0 or cdim <= 0: + matrix_rel_fro_aligned[name] = float("nan") + continue + ref_crop = ref_m[:rdim, :rdim] + cpp_crop = cpp_m[:cdim, :cdim] + t_crop = transform[:cdim, :rdim] + cpp_aligned = transform_weak_to_reference_basis(cpp_crop, t_crop) + matrix_rel_fro_aligned[name] = rel_fro(ref_crop, cpp_aligned) + + eval_alignment: dict[str, dict[str, float]] = {} + for form in [0, 1]: + rp = ref_dir / f"reference_circular_operator_form{form}_evals.csv" + cp = cpp_dir / f"circular_operator_form{form}_evals.csv" + if not rp.exists() or not cp.exists(): + eval_alignment[f"form{form}"] = { + "mean_abs_diff": float("nan"), + "p95_abs_diff": float("nan"), + } + continue + eval_alignment[f"form{form}"] = eval_alignment_error( + load_positive_evals(rp), load_positive_evals(cp) + ) + + payload: dict[str, Any] = { + "reference_dir": str(ref_dir.resolve()), + "cpp_dir": str(cpp_dir.resolve()), + "matrix_rel_fro": matrix_rel_fro, + "matrix_rel_fro_aligned": matrix_rel_fro_aligned, + "basis_alignment": basis_alignment, + "eval_alignment": eval_alignment, + } + + out_md = Path(args.output_markdown) + out_json = Path(args.output_json) + out_md.parent.mkdir(parents=True, exist_ok=True) + out_json.parent.mkdir(parents=True, exist_ok=True) + out_md.write_text(render_markdown(payload)) + out_json.write_text(json.dumps(payload, indent=2)) + + print(f"Wrote debug markdown: {out_md}") + print(f"Wrote debug json: {out_json}") + + +if __name__ == "__main__": + main() diff --git a/scripts/hodge/run_reference_hodge.py b/scripts/hodge/run_reference_hodge.py index 9210a6e..ef7fd71 100755 --- a/scripts/hodge/run_reference_hodge.py +++ b/scripts/hodge/run_reference_hodge.py @@ -31,6 +31,28 @@ def write_table(path: pathlib.Path, header: List[str], rows: np.ndarray) -> None writer.writerow([float(v) for v in row]) +def write_matrix_csv(path: pathlib.Path, matrix: np.ndarray) -> None: + matrix = np.asarray(matrix, dtype=float) + rows, cols = matrix.shape + header = ["row"] + [f"c{i}" for i in range(cols)] + payload = np.column_stack([np.arange(rows, dtype=float), matrix]) + write_table(path, header, payload) + + +def write_complex_evals_csv(path: pathlib.Path, evals: np.ndarray) -> None: + evals = np.asarray(evals, dtype=complex) + rows = np.column_stack( + [ + np.arange(evals.shape[0], dtype=float), + np.real(evals), + np.imag(evals), + np.abs(evals), + (np.imag(evals) > 0).astype(float), + ] + ) + write_table(path, ["index", "real", "imag", "abs", "positive_imag"], rows) + + def compute_deterministic_basis( kernel: np.ndarray, nbr_indices: np.ndarray, n_function_basis: int ) -> tuple[np.ndarray, np.ndarray]: @@ -73,7 +95,7 @@ def circular_coordinates(form, lam: float, mode: int): circular_indices = np.where(evals.imag > 0)[0] if circular_indices.size == 0: zeros = np.zeros(dg.n, dtype=float) - return zeros, complex(0.0, 0.0), -1 + return zeros, complex(0.0, 0.0), -1, evals mode_idx = min(max(mode, 0), circular_indices.size - 1) chosen_global = int(circular_indices[mode_idx]) @@ -82,7 +104,7 @@ def circular_coordinates(form, lam: float, mode: int): circular_funcs = efunctions[circular_indices].to_ambient() chosen_func = circular_funcs[mode_idx] angles = np.mod(np.arctan2(chosen_func.imag, chosen_func.real), 2.0 * np.pi) - return angles, chosen_eval, chosen_global + return angles, chosen_eval, chosen_global, evals def main() -> None: @@ -143,6 +165,10 @@ def main() -> None: harmonic_ambient = [] circular_thetas = [] circular_meta = [] + laplacian0_weak = np.asarray(dg.laplacian(0).weak) + function_gram = np.asarray(dg.function_space.gram) + write_matrix_csv(out_dir / "reference_function_gram.csv", function_gram) + write_matrix_csv(out_dir / "reference_laplacian0_weak.csv", laplacian0_weak) for idx, form_idx in enumerate(form_indices): form = forms_1[form_idx] @@ -150,7 +176,7 @@ def main() -> None: ambient = np.asarray(form.to_ambient(), dtype=float) harmonic_ambient.append(ambient) - theta, eigval, eig_idx = circular_coordinates( + theta, eigval, eig_idx, evals = circular_coordinates( form, args.circular_lambda, mode_indices[idx] ) circular_thetas.append(theta) @@ -166,11 +192,23 @@ def main() -> None: } ) + x_weak = np.asarray(form.sharp().operator.weak) + operator_weak = x_weak - args.circular_lambda * laplacian0_weak + write_matrix_csv(out_dir / f"reference_circular_operator_form{idx}_x_weak.csv", x_weak) + write_matrix_csv( + out_dir / f"reference_circular_operator_form{idx}_operator_weak.csv", + operator_weak, + ) + write_complex_evals_csv( + out_dir / f"reference_circular_operator_form{idx}_evals.csv", evals + ) + write_table( out_dir / "reference_points.csv", ["x", "y", "z"], points, ) + write_matrix_csv(out_dir / "reference_function_basis.csv", function_basis) spectrum_rows = np.column_stack([np.arange(vals_1.shape[0]), np.asarray(vals_1, dtype=float)]) write_table( diff --git a/src/main_hodge.cpp b/src/main_hodge.cpp index e16e717..bbebed3 100644 --- a/src/main_hodge.cpp +++ b/src/main_hodge.cpp @@ -281,6 +281,42 @@ reconstruct_harmonic_ambient(const DiffusionMesh &mesh, return field; } +static double orientation_score(const DiffusionMesh &mesh, + const std::vector &field) { + const size_t n_verts = mesh.geometry.num_points(); + if (n_verts == 0) { + return 0.0; + } + + double accum = 0.0; + for (size_t i = 0; i < n_verts; ++i) { + const auto p = mesh.geometry.get_vec3(i); + const auto v = field[i]; + accum += static_cast(p.x) * static_cast(v.x) + + static_cast(p.y) * static_cast(v.y) + + static_cast(p.z) * static_cast(v.z); + } + + return accum / static_cast(n_verts); +} + +static bool canonicalize_form_sign(const DiffusionMesh &mesh, + Eigen::Ref coeffs, + std::vector &ambient_field) { + ambient_field = reconstruct_harmonic_ambient(mesh, coeffs); + if (orientation_score(mesh, ambient_field) >= 0.0) { + return false; + } + + coeffs *= -1.0f; + for (auto &v : ambient_field) { + v.x *= -1.0f; + v.y *= -1.0f; + v.z *= -1.0f; + } + return true; +} + static void export_vector_field(const std::string &filename, const DiffusionMesh &mesh, const std::vector &vectors) { @@ -353,6 +389,43 @@ static void write_circular_modes_csv(const std::string &filename, << "," << eval_1.imag() << "\n"; } +static void write_matrix_csv(const std::string &filename, + const Eigen::MatrixXf &matrix) { + std::ofstream file(filename); + file << "row"; + for (int j = 0; j < matrix.cols(); ++j) { + file << ",c" << j; + } + file << "\n"; + + for (int i = 0; i < matrix.rows(); ++i) { + file << i; + for (int j = 0; j < matrix.cols(); ++j) { + file << "," << matrix(i, j); + } + file << "\n"; + } +} + +static void write_complex_evals_csv( + const std::string &filename, const Eigen::VectorXcf &evals, + const std::vector &positive_imag_indices) { + std::vector is_positive(static_cast(evals.size()), 0); + for (int idx : positive_imag_indices) { + if (idx >= 0 && idx < evals.size()) { + is_positive[static_cast(idx)] = 1; + } + } + + std::ofstream file(filename); + file << "index,real,imag,abs,positive_imag\n"; + for (int i = 0; i < evals.size(); ++i) { + const auto v = evals[i]; + file << i << "," << v.real() << "," << v.imag() << "," << std::abs(v) << "," + << is_positive[static_cast(i)] << "\n"; + } +} + int main(int argc, char **argv) { Config cfg; if (!parse_args(argc, argv, cfg)) { @@ -398,14 +471,25 @@ int main(int argc, char **argv) { return 1; } + const int export_forms = std::min(3, static_cast(evecs.cols())); + std::vector> harmonic_fields; + harmonic_fields.reserve(static_cast(export_forms)); + for (int i = 0; i < export_forms; ++i) { + std::vector field; + canonicalize_form_sign(mesh, evecs.col(i), field); + harmonic_fields.push_back(std::move(field)); + } + std::complex selected_eval_0(0.0f, 0.0f); std::complex selected_eval_1(0.0f, 0.0f); + ops::CircularOperatorDiagnostics circular_diag_0; + ops::CircularOperatorDiagnostics circular_diag_1; const auto theta_0 = ops::compute_circular_coordinates( mesh, evecs.col(0), 0.0f, cfg.circular_lambda, cfg.circular_mode_0, - &selected_eval_0); + &selected_eval_0, &circular_diag_0); const auto theta_1 = ops::compute_circular_coordinates( mesh, evecs.col(1), 0.0f, cfg.circular_lambda, cfg.circular_mode_1, - &selected_eval_1); + &selected_eval_1, &circular_diag_1); std::cout << "HODGE SPECTRUM (first 12 modes)\n"; for (int i = 0; i < 12 && i < evals.size(); ++i) { @@ -415,15 +499,8 @@ int main(int argc, char **argv) { write_points_csv(cfg.output_dir + "/points.csv", mesh); write_spectrum_csv(cfg.output_dir + "/hodge_spectrum.csv", evals); - const int export_forms = std::min(3, static_cast(evecs.cols())); write_harmonic_coeffs_csv(cfg.output_dir + "/harmonic_coeffs.csv", evecs, export_forms); - - std::vector> harmonic_fields; - harmonic_fields.reserve(static_cast(export_forms)); - for (int i = 0; i < export_forms; ++i) { - harmonic_fields.push_back(reconstruct_harmonic_ambient(mesh, evecs.col(i))); - } write_harmonic_ambient_csv(cfg.output_dir + "/harmonic_ambient.csv", mesh, harmonic_fields); @@ -432,6 +509,28 @@ int main(int argc, char **argv) { write_circular_modes_csv(cfg.output_dir + "/circular_modes.csv", selected_eval_0, selected_eval_1, cfg.circular_mode_0, cfg.circular_mode_1, cfg.circular_lambda); + write_matrix_csv(cfg.output_dir + "/function_gram.csv", + circular_diag_0.function_gram); + write_matrix_csv(cfg.output_dir + "/laplacian0_weak.csv", + circular_diag_0.laplacian0_weak); + write_matrix_csv(cfg.output_dir + "/circular_operator_form0_x_weak.csv", + circular_diag_0.x_weak); + write_matrix_csv(cfg.output_dir + "/circular_operator_form1_x_weak.csv", + circular_diag_1.x_weak); + write_matrix_csv(cfg.output_dir + "/circular_operator_form0_operator_weak.csv", + circular_diag_0.x_weak - + (cfg.circular_lambda * circular_diag_0.laplacian0_weak)); + write_matrix_csv(cfg.output_dir + "/circular_operator_form1_operator_weak.csv", + circular_diag_1.x_weak - + (cfg.circular_lambda * circular_diag_1.laplacian0_weak)); + write_complex_evals_csv(cfg.output_dir + "/circular_operator_form0_evals.csv", + circular_diag_0.sorted_evals, + circular_diag_0.positive_imag_indices); + write_complex_evals_csv(cfg.output_dir + "/circular_operator_form1_evals.csv", + circular_diag_1.sorted_evals, + circular_diag_1.positive_imag_indices); + write_matrix_csv(cfg.output_dir + "/function_basis.csv", + mesh.topology.eigen_basis); if (export_ply) { std::vector field_0(static_cast(theta_0.size())); From a32521a0bab2de2a7fc1d9a6bc8dc476a89ad28b Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Sat, 14 Feb 2026 12:40:56 -0700 Subject: [PATCH 9/9] Clean hodge runtime outputs and harden parity regressions --- .github/workflows/ci.yml | 2 + CMakeLists.txt | 11 + README.md | 26 ++ include/igneous/ops/hodge.hpp | 23 +- .../hodge/compare_circular_operator_debug.py | 243 ------------------ .../{ => diagnostics}/plot_hodge_outputs.py | 0 scripts/hodge/run_reference_hodge.py | 41 +-- src/main_hodge.cpp | 65 +---- tests/test_hodge_cli_outputs.sh | 67 +++++ tests/test_hodge_parity_optional.sh | 52 ++++ tests/test_ops_hodge.cpp | 42 ++- 11 files changed, 219 insertions(+), 353 deletions(-) delete mode 100755 scripts/hodge/compare_circular_operator_debug.py rename scripts/hodge/{ => diagnostics}/plot_hodge_outputs.py (100%) create mode 100755 tests/test_hodge_cli_outputs.sh create mode 100755 tests/test_hodge_parity_optional.sh diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e28d7cc..02d1c67 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,6 +70,7 @@ jobs: test_ops_spectral_geometry test_ops_hodge test_io_meshes + igneous-hodge - name: Run tests shell: bash @@ -145,6 +146,7 @@ jobs: test_ops_spectral_geometry test_ops_hodge test_io_meshes + igneous-hodge - name: Run tests (Sanitizers) shell: bash diff --git a/CMakeLists.txt b/CMakeLists.txt index e3b9189..d97f674 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -85,3 +85,14 @@ target_link_libraries(igneous-spectral PRIVATE igneous fmt::fmt) add_executable(igneous-hodge src/main_hodge.cpp) target_link_libraries(igneous-hodge PRIVATE igneous fmt::fmt) + +add_test( + NAME test_hodge_cli_outputs + COMMAND bash ${CMAKE_SOURCE_DIR}/tests/test_hodge_cli_outputs.sh + $) +set_tests_properties(test_hodge_cli_outputs PROPERTIES LABELS "hodge;cli") + +add_test( + NAME test_hodge_parity_optional + COMMAND bash ${CMAKE_SOURCE_DIR}/tests/test_hodge_parity_optional.sh) +set_tests_properties(test_hodge_parity_optional PROPERTIES LABELS "hodge;parity") diff --git a/README.md b/README.md index 2140a0f..f027262 100644 --- a/README.md +++ b/README.md @@ -46,6 +46,30 @@ cmake --build build -j8 ./build/igneous-hodge ``` +## Reference Implementation + +The diffusion/Hodge parity work in this repository is aligned against the Python +reference implementation: + +- [DiffusionGeometry](https://github.com/Iolo-Jones/DiffusionGeometry) + +## Hodge Parity Workflow + +Standard parity round: + +```bash +./scripts/hodge/run_parity_round.sh +``` + +Optional diagnostic plots for a parity round: + +```bash +ROUND_DIR="$(ls -1dt notes/hodge/results/round_* | head -n1)" +notes/hodge/.venv_ref/bin/python \ + ./scripts/hodge/diagnostics/plot_hodge_outputs.py \ + --round-dir "${ROUND_DIR}" +``` + Set `IGNEOUS_BENCH_MODE=1` to disable heavy export paths in runtime apps. Runtime backend controls: @@ -74,6 +98,8 @@ Current suites: - `test_ops_spectral_geometry` - `test_ops_hodge` - `test_io_meshes` +- `test_hodge_cli_outputs` +- `test_hodge_parity_optional` (skips unless `DiffusionGeometry/` is available, unless forced) ## Benchmarks diff --git a/include/igneous/ops/hodge.hpp b/include/igneous/ops/hodge.hpp index 979566d..0c3c719 100644 --- a/include/igneous/ops/hodge.hpp +++ b/include/igneous/ops/hodge.hpp @@ -201,23 +201,13 @@ compute_hodge_spectrum(const Eigen::MatrixXf &laplacian, return std::make_pair(evals, evecs); } -template -struct CircularOperatorDiagnostics { - Eigen::MatrixXf x_weak; - Eigen::MatrixXf laplacian0_weak; - Eigen::MatrixXf function_gram; - Eigen::VectorXcf sorted_evals; - std::vector positive_imag_indices; -}; - template Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, const Eigen::VectorXf &alpha_coeffs, float bandwidth, float lambda = 1.0f, int positive_imag_mode = 0, - std::complex *selected_eval = nullptr, - CircularOperatorDiagnostics *diagnostics = nullptr) { + std::complex *selected_eval = nullptr) { const auto &U = mesh.topology.eigen_basis; const auto &mu = mesh.topology.mu; @@ -276,12 +266,6 @@ Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, Eigen::MatrixXf function_gram = U.transpose() * (U.array().colwise() * mu.array()).matrix(); Eigen::MatrixXf operator_weak = X_op - (lambda * laplacian0_weak); - if (diagnostics != nullptr) { - diagnostics->x_weak = X_op; - diagnostics->laplacian0_weak = laplacian0_weak; - diagnostics->function_gram = function_gram; - } - Eigen::GeneralizedEigenSolver solver(operator_weak, function_gram); if (solver.info() != Eigen::Success) { return Eigen::VectorXf::Zero(static_cast(n_verts)); @@ -322,11 +306,6 @@ Eigen::VectorXf compute_circular_coordinates(const MeshT &mesh, positive_imag_indices.push_back(i); } } - if (diagnostics != nullptr) { - diagnostics->sorted_evals = sorted_evals; - diagnostics->positive_imag_indices = positive_imag_indices; - } - if (positive_imag_indices.empty()) { return Eigen::VectorXf::Zero(static_cast(n_verts)); } diff --git a/scripts/hodge/compare_circular_operator_debug.py b/scripts/hodge/compare_circular_operator_debug.py deleted file mode 100755 index 181da31..0000000 --- a/scripts/hodge/compare_circular_operator_debug.py +++ /dev/null @@ -1,243 +0,0 @@ -#!/usr/bin/env python3 -import argparse -import json -from pathlib import Path -from typing import Any - -import numpy as np - - -def load_columns(path: Path) -> dict[str, np.ndarray]: - arr = np.genfromtxt(path, delimiter=",", names=True) - if arr.shape == (): - arr = np.array([arr], dtype=arr.dtype) - return {name: np.asarray(arr[name]) for name in arr.dtype.names} - - -def load_matrix(path: Path) -> np.ndarray: - cols = load_columns(path) - keys = [k for k in cols.keys() if k.startswith("c")] - keys.sort(key=lambda x: int(x[1:])) - return np.column_stack([cols[k] for k in keys]).astype(np.float64) - - -def load_evals(path: Path) -> np.ndarray: - cols = load_columns(path) - return cols["real"].astype(np.float64) + 1j * cols["imag"].astype(np.float64) - - -def load_positive_evals(path: Path) -> np.ndarray: - cols = load_columns(path) - evals = cols["real"].astype(np.float64) + 1j * cols["imag"].astype(np.float64) - mask = cols["positive_imag"].astype(np.float64) > 0.5 - return evals[mask] - - -def rel_fro(a: np.ndarray, b: np.ndarray) -> float: - n = min(a.shape[0], b.shape[0]) - m = min(a.shape[1], b.shape[1]) - aa = a[:n, :m] - bb = b[:n, :m] - denom = max(np.linalg.norm(aa, ord="fro"), 1e-12) - return float(np.linalg.norm(aa - bb, ord="fro") / denom) - - -def align_basis( - reference_basis: np.ndarray, cpp_basis: np.ndarray -) -> tuple[np.ndarray, dict[str, float]]: - n = min(reference_basis.shape[0], cpp_basis.shape[0]) - k_ref = reference_basis.shape[1] - k_cpp = cpp_basis.shape[1] - - ref = reference_basis[:n, :k_ref] - cpp = cpp_basis[:n, :k_cpp] - - transform, _, rank, singular_values = np.linalg.lstsq(cpp, ref, rcond=None) - reconstructed = cpp @ transform - denom = max(np.linalg.norm(ref, ord="fro"), 1e-12) - recon_rel = float(np.linalg.norm(reconstructed - ref, ord="fro") / denom) - - cond = float("inf") - if singular_values.size > 0: - cond = float(singular_values[0] / max(singular_values[-1], 1e-12)) - - stats = { - "rows_used": float(n), - "reference_cols": float(k_ref), - "cpp_cols": float(k_cpp), - "transform_rank": float(rank), - "cpp_basis_condition_estimate": cond, - "reconstruction_rel_fro": recon_rel, - } - return transform, stats - - -def transform_weak_to_reference_basis(cpp_weak: np.ndarray, transform: np.ndarray) -> np.ndarray: - return transform.conj().T @ cpp_weak @ transform - - -def eval_alignment_error(ref: np.ndarray, cpp: np.ndarray) -> dict[str, float]: - k = min(ref.shape[0], cpp.shape[0]) - if k == 0: - return {"mean_abs_diff": float("nan"), "p95_abs_diff": float("nan")} - d = np.abs(ref[:k] - cpp[:k]) - return { - "mean_abs_diff": float(np.mean(d)), - "p95_abs_diff": float(np.percentile(d, 95.0)), - } - - -def render_markdown(payload: dict[str, Any]) -> str: - lines: list[str] = [] - lines.append("# Circular Operator Debug Comparison") - lines.append("") - lines.append(f"- Reference dir: `{payload['reference_dir']}`") - lines.append(f"- C++ dir: `{payload['cpp_dir']}`") - lines.append("") - lines.append("## Matrix Relative Frobenius Errors") - lines.append("") - for k, v in payload["matrix_rel_fro"].items(): - lines.append(f"- {k}: `{v:.6e}`") - lines.append("") - if payload.get("matrix_rel_fro_aligned"): - lines.append("## Matrix Relative Frobenius Errors (Basis Aligned)") - lines.append("") - for k, v in payload["matrix_rel_fro_aligned"].items(): - lines.append(f"- {k}: `{v:.6e}`") - lines.append("") - if payload.get("basis_alignment"): - lines.append("## Basis Alignment") - lines.append("") - stats = payload["basis_alignment"] - lines.append(f"- rows used: `{int(stats['rows_used'])}`") - lines.append(f"- reference columns: `{int(stats['reference_cols'])}`") - lines.append(f"- cpp columns: `{int(stats['cpp_cols'])}`") - lines.append(f"- transform rank: `{int(stats['transform_rank'])}`") - lines.append( - f"- cpp basis condition estimate: `{stats['cpp_basis_condition_estimate']:.6e}`" - ) - lines.append( - f"- reconstruction relative Frobenius error: `{stats['reconstruction_rel_fro']:.6e}`" - ) - lines.append("") - - lines.append("## Positive-Imag Eigenvalue Alignment") - lines.append("") - for form, stats in payload["eval_alignment"].items(): - lines.append( - f"- {form}: mean abs diff `{stats['mean_abs_diff']:.6e}`, " - f"P95 abs diff `{stats['p95_abs_diff']:.6e}`" - ) - lines.append("") - return "\n".join(lines) - - -def main() -> None: - parser = argparse.ArgumentParser( - description="Compare circular operator debug artifacts between reference and C++ outputs." - ) - parser.add_argument("--reference-dir", required=True) - parser.add_argument("--cpp-dir", required=True) - parser.add_argument("--output-markdown", required=True) - parser.add_argument("--output-json", required=True) - args = parser.parse_args() - - ref_dir = Path(args.reference_dir) - cpp_dir = Path(args.cpp_dir) - - matrix_pairs = { - "function_gram": ( - ref_dir / "reference_function_gram.csv", - cpp_dir / "function_gram.csv", - ), - "laplacian0_weak": ( - ref_dir / "reference_laplacian0_weak.csv", - cpp_dir / "laplacian0_weak.csv", - ), - "form0_x_weak": ( - ref_dir / "reference_circular_operator_form0_x_weak.csv", - cpp_dir / "circular_operator_form0_x_weak.csv", - ), - "form1_x_weak": ( - ref_dir / "reference_circular_operator_form1_x_weak.csv", - cpp_dir / "circular_operator_form1_x_weak.csv", - ), - "form0_operator_weak": ( - ref_dir / "reference_circular_operator_form0_operator_weak.csv", - cpp_dir / "circular_operator_form0_operator_weak.csv", - ), - "form1_operator_weak": ( - ref_dir / "reference_circular_operator_form1_operator_weak.csv", - cpp_dir / "circular_operator_form1_operator_weak.csv", - ), - } - - matrix_rel_fro: dict[str, float] = {} - for name, (rp, cp) in matrix_pairs.items(): - if not rp.exists() or not cp.exists(): - matrix_rel_fro[name] = float("nan") - continue - matrix_rel_fro[name] = rel_fro(load_matrix(rp), load_matrix(cp)) - - matrix_rel_fro_aligned: dict[str, float] = {} - basis_alignment: dict[str, float] = {} - ref_basis_path = ref_dir / "reference_function_basis.csv" - cpp_basis_path = cpp_dir / "function_basis.csv" - if ref_basis_path.exists() and cpp_basis_path.exists(): - ref_basis = load_matrix(ref_basis_path) - cpp_basis = load_matrix(cpp_basis_path) - transform, basis_alignment = align_basis(ref_basis, cpp_basis) - - for name, (rp, cp) in matrix_pairs.items(): - if not rp.exists() or not cp.exists(): - matrix_rel_fro_aligned[name] = float("nan") - continue - ref_m = load_matrix(rp) - cpp_m = load_matrix(cp) - rdim = min(ref_m.shape[0], ref_m.shape[1], transform.shape[1]) - cdim = min(cpp_m.shape[0], cpp_m.shape[1], transform.shape[0]) - if rdim <= 0 or cdim <= 0: - matrix_rel_fro_aligned[name] = float("nan") - continue - ref_crop = ref_m[:rdim, :rdim] - cpp_crop = cpp_m[:cdim, :cdim] - t_crop = transform[:cdim, :rdim] - cpp_aligned = transform_weak_to_reference_basis(cpp_crop, t_crop) - matrix_rel_fro_aligned[name] = rel_fro(ref_crop, cpp_aligned) - - eval_alignment: dict[str, dict[str, float]] = {} - for form in [0, 1]: - rp = ref_dir / f"reference_circular_operator_form{form}_evals.csv" - cp = cpp_dir / f"circular_operator_form{form}_evals.csv" - if not rp.exists() or not cp.exists(): - eval_alignment[f"form{form}"] = { - "mean_abs_diff": float("nan"), - "p95_abs_diff": float("nan"), - } - continue - eval_alignment[f"form{form}"] = eval_alignment_error( - load_positive_evals(rp), load_positive_evals(cp) - ) - - payload: dict[str, Any] = { - "reference_dir": str(ref_dir.resolve()), - "cpp_dir": str(cpp_dir.resolve()), - "matrix_rel_fro": matrix_rel_fro, - "matrix_rel_fro_aligned": matrix_rel_fro_aligned, - "basis_alignment": basis_alignment, - "eval_alignment": eval_alignment, - } - - out_md = Path(args.output_markdown) - out_json = Path(args.output_json) - out_md.parent.mkdir(parents=True, exist_ok=True) - out_json.parent.mkdir(parents=True, exist_ok=True) - out_md.write_text(render_markdown(payload)) - out_json.write_text(json.dumps(payload, indent=2)) - - print(f"Wrote debug markdown: {out_md}") - print(f"Wrote debug json: {out_json}") - - -if __name__ == "__main__": - main() diff --git a/scripts/hodge/plot_hodge_outputs.py b/scripts/hodge/diagnostics/plot_hodge_outputs.py similarity index 100% rename from scripts/hodge/plot_hodge_outputs.py rename to scripts/hodge/diagnostics/plot_hodge_outputs.py diff --git a/scripts/hodge/run_reference_hodge.py b/scripts/hodge/run_reference_hodge.py index ef7fd71..0a2cb87 100755 --- a/scripts/hodge/run_reference_hodge.py +++ b/scripts/hodge/run_reference_hodge.py @@ -120,6 +120,11 @@ def main() -> None: parser.add_argument("--circular-lambda", type=float, default=1.0) parser.add_argument("--form-indices", default="0,1") parser.add_argument("--mode-indices", default="0,1") + parser.add_argument( + "--emit-operator-debug", + action="store_true", + help="Emit operator-level debug artifacts for diagnostics", + ) args = parser.parse_args() out_dir = pathlib.Path(args.output_dir) @@ -165,10 +170,12 @@ def main() -> None: harmonic_ambient = [] circular_thetas = [] circular_meta = [] - laplacian0_weak = np.asarray(dg.laplacian(0).weak) - function_gram = np.asarray(dg.function_space.gram) - write_matrix_csv(out_dir / "reference_function_gram.csv", function_gram) - write_matrix_csv(out_dir / "reference_laplacian0_weak.csv", laplacian0_weak) + laplacian0_weak = None + if args.emit_operator_debug: + laplacian0_weak = np.asarray(dg.laplacian(0).weak) + function_gram = np.asarray(dg.function_space.gram) + write_matrix_csv(out_dir / "reference_function_gram.csv", function_gram) + write_matrix_csv(out_dir / "reference_laplacian0_weak.csv", laplacian0_weak) for idx, form_idx in enumerate(form_indices): form = forms_1[form_idx] @@ -192,23 +199,27 @@ def main() -> None: } ) - x_weak = np.asarray(form.sharp().operator.weak) - operator_weak = x_weak - args.circular_lambda * laplacian0_weak - write_matrix_csv(out_dir / f"reference_circular_operator_form{idx}_x_weak.csv", x_weak) - write_matrix_csv( - out_dir / f"reference_circular_operator_form{idx}_operator_weak.csv", - operator_weak, - ) - write_complex_evals_csv( - out_dir / f"reference_circular_operator_form{idx}_evals.csv", evals - ) + if args.emit_operator_debug: + x_weak = np.asarray(form.sharp().operator.weak) + operator_weak = x_weak - args.circular_lambda * laplacian0_weak + write_matrix_csv( + out_dir / f"reference_circular_operator_form{idx}_x_weak.csv", x_weak + ) + write_matrix_csv( + out_dir / f"reference_circular_operator_form{idx}_operator_weak.csv", + operator_weak, + ) + write_complex_evals_csv( + out_dir / f"reference_circular_operator_form{idx}_evals.csv", evals + ) write_table( out_dir / "reference_points.csv", ["x", "y", "z"], points, ) - write_matrix_csv(out_dir / "reference_function_basis.csv", function_basis) + if args.emit_operator_debug: + write_matrix_csv(out_dir / "reference_function_basis.csv", function_basis) spectrum_rows = np.column_stack([np.arange(vals_1.shape[0]), np.asarray(vals_1, dtype=float)]) write_table( diff --git a/src/main_hodge.cpp b/src/main_hodge.cpp index bbebed3..c2675fd 100644 --- a/src/main_hodge.cpp +++ b/src/main_hodge.cpp @@ -389,43 +389,6 @@ static void write_circular_modes_csv(const std::string &filename, << "," << eval_1.imag() << "\n"; } -static void write_matrix_csv(const std::string &filename, - const Eigen::MatrixXf &matrix) { - std::ofstream file(filename); - file << "row"; - for (int j = 0; j < matrix.cols(); ++j) { - file << ",c" << j; - } - file << "\n"; - - for (int i = 0; i < matrix.rows(); ++i) { - file << i; - for (int j = 0; j < matrix.cols(); ++j) { - file << "," << matrix(i, j); - } - file << "\n"; - } -} - -static void write_complex_evals_csv( - const std::string &filename, const Eigen::VectorXcf &evals, - const std::vector &positive_imag_indices) { - std::vector is_positive(static_cast(evals.size()), 0); - for (int idx : positive_imag_indices) { - if (idx >= 0 && idx < evals.size()) { - is_positive[static_cast(idx)] = 1; - } - } - - std::ofstream file(filename); - file << "index,real,imag,abs,positive_imag\n"; - for (int i = 0; i < evals.size(); ++i) { - const auto v = evals[i]; - file << i << "," << v.real() << "," << v.imag() << "," << std::abs(v) << "," - << is_positive[static_cast(i)] << "\n"; - } -} - int main(int argc, char **argv) { Config cfg; if (!parse_args(argc, argv, cfg)) { @@ -482,14 +445,12 @@ int main(int argc, char **argv) { std::complex selected_eval_0(0.0f, 0.0f); std::complex selected_eval_1(0.0f, 0.0f); - ops::CircularOperatorDiagnostics circular_diag_0; - ops::CircularOperatorDiagnostics circular_diag_1; const auto theta_0 = ops::compute_circular_coordinates( mesh, evecs.col(0), 0.0f, cfg.circular_lambda, cfg.circular_mode_0, - &selected_eval_0, &circular_diag_0); + &selected_eval_0); const auto theta_1 = ops::compute_circular_coordinates( mesh, evecs.col(1), 0.0f, cfg.circular_lambda, cfg.circular_mode_1, - &selected_eval_1, &circular_diag_1); + &selected_eval_1); std::cout << "HODGE SPECTRUM (first 12 modes)\n"; for (int i = 0; i < 12 && i < evals.size(); ++i) { @@ -509,28 +470,6 @@ int main(int argc, char **argv) { write_circular_modes_csv(cfg.output_dir + "/circular_modes.csv", selected_eval_0, selected_eval_1, cfg.circular_mode_0, cfg.circular_mode_1, cfg.circular_lambda); - write_matrix_csv(cfg.output_dir + "/function_gram.csv", - circular_diag_0.function_gram); - write_matrix_csv(cfg.output_dir + "/laplacian0_weak.csv", - circular_diag_0.laplacian0_weak); - write_matrix_csv(cfg.output_dir + "/circular_operator_form0_x_weak.csv", - circular_diag_0.x_weak); - write_matrix_csv(cfg.output_dir + "/circular_operator_form1_x_weak.csv", - circular_diag_1.x_weak); - write_matrix_csv(cfg.output_dir + "/circular_operator_form0_operator_weak.csv", - circular_diag_0.x_weak - - (cfg.circular_lambda * circular_diag_0.laplacian0_weak)); - write_matrix_csv(cfg.output_dir + "/circular_operator_form1_operator_weak.csv", - circular_diag_1.x_weak - - (cfg.circular_lambda * circular_diag_1.laplacian0_weak)); - write_complex_evals_csv(cfg.output_dir + "/circular_operator_form0_evals.csv", - circular_diag_0.sorted_evals, - circular_diag_0.positive_imag_indices); - write_complex_evals_csv(cfg.output_dir + "/circular_operator_form1_evals.csv", - circular_diag_1.sorted_evals, - circular_diag_1.positive_imag_indices); - write_matrix_csv(cfg.output_dir + "/function_basis.csv", - mesh.topology.eigen_basis); if (export_ply) { std::vector field_0(static_cast(theta_0.size())); diff --git a/tests/test_hodge_cli_outputs.sh b/tests/test_hodge_cli_outputs.sh new file mode 100755 index 0000000..ed1577e --- /dev/null +++ b/tests/test_hodge_cli_outputs.sh @@ -0,0 +1,67 @@ +#!/usr/bin/env bash +set -euo pipefail + +if [[ $# -lt 1 ]]; then + echo "usage: $0 " >&2 + exit 2 +fi + +HODGE_EXE="$1" +if [[ ! -x "${HODGE_EXE}" ]]; then + echo "missing executable: ${HODGE_EXE}" >&2 + exit 2 +fi + +TMP_DIR="$(mktemp -d)" +trap 'rm -rf "${TMP_DIR}"' EXIT +OUT_DIR="${TMP_DIR}/out" + +"${HODGE_EXE}" \ + --output-dir "${OUT_DIR}" \ + --n-points 160 \ + --n-basis 24 \ + --k-neighbors 24 \ + --knn-bandwidth 8 \ + --bandwidth-variability -0.5 \ + --c 0.0 \ + --circular-lambda 1.0 \ + --circular-mode-0 0 \ + --circular-mode-1 1 \ + --no-ply >/dev/null + +required_files=( + "points.csv" + "hodge_spectrum.csv" + "harmonic_coeffs.csv" + "harmonic_ambient.csv" + "circular_coordinates.csv" + "circular_modes.csv" +) + +for rel in "${required_files[@]}"; do + if [[ ! -f "${OUT_DIR}/${rel}" ]]; then + echo "missing required output file: ${rel}" >&2 + exit 1 + fi +done + +forbidden_files=( + "function_gram.csv" + "laplacian0_weak.csv" + "function_basis.csv" + "circular_operator_form0_x_weak.csv" + "circular_operator_form1_x_weak.csv" + "circular_operator_form0_operator_weak.csv" + "circular_operator_form1_operator_weak.csv" + "circular_operator_form0_evals.csv" + "circular_operator_form1_evals.csv" +) + +for rel in "${forbidden_files[@]}"; do + if [[ -f "${OUT_DIR}/${rel}" ]]; then + echo "found unexpected diagnostic output file: ${rel}" >&2 + exit 1 + fi +done + +echo "hodge CLI output contract check passed" diff --git a/tests/test_hodge_parity_optional.sh b/tests/test_hodge_parity_optional.sh new file mode 100755 index 0000000..462ed55 --- /dev/null +++ b/tests/test_hodge_parity_optional.sh @@ -0,0 +1,52 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +REQUIRE_PARITY="${IGNEOUS_REQUIRE_PARITY:-0}" +REF_DIR="${ROOT_DIR}/DiffusionGeometry" + +if [[ ! -d "${REF_DIR}" && "${REQUIRE_PARITY}" != "1" ]]; then + echo "optional parity test skipped: ${REF_DIR} not found" + exit 0 +fi + +if [[ ! -d "${REF_DIR}" && "${REQUIRE_PARITY}" == "1" ]]; then + echo "parity required but reference repo is missing: ${REF_DIR}" >&2 + exit 1 +fi + +run_output="$("${ROOT_DIR}/scripts/hodge/run_parity_round.sh")" +echo "${run_output}" + +round_dir="$(echo "${run_output}" | awk -F= '/^round_dir=/{print $2}' | tail -n1)" +if [[ -z "${round_dir}" ]]; then + echo "failed to parse round_dir from parity output" >&2 + exit 1 +fi + +report_json="${round_dir}/report/parity_report.json" +if [[ ! -f "${report_json}" ]]; then + echo "missing parity report: ${report_json}" >&2 + exit 1 +fi + +python3 - "${report_json}" <<'PY' +import json +import pathlib +import sys + +report_path = pathlib.Path(sys.argv[1]) +payload = json.loads(report_path.read_text()) +gates = payload.get("gates", {}) +summary = payload.get("summary", {}) + +if not gates.get("final_pass", False): + score = summary.get("composite_score", "nan") + corr = summary.get("circular_complex_correlation_min", "nan") + p95 = summary.get("circular_wrapped_p95_rad_max", "nan") + raise SystemExit( + f"parity final gate failed: composite={score}, corr_min={corr}, p95={p95}" + ) + +print("optional parity regression check passed") +PY diff --git a/tests/test_ops_hodge.cpp b/tests/test_ops_hodge.cpp index bcb2424..44966df 100644 --- a/tests/test_ops_hodge.cpp +++ b/tests/test_ops_hodge.cpp @@ -31,8 +31,8 @@ make_torus_cloud(size_t n_points) { } TEST_CASE("Hodge operators produce finite spectrum and circular coordinates") { - auto mesh = make_torus_cloud(320); - igneous::ops::compute_eigenbasis(mesh, 10); + auto mesh = make_torus_cloud(1000); + igneous::ops::compute_eigenbasis(mesh, 50); igneous::ops::GeometryWorkspace geom_ws; const auto G = igneous::ops::compute_1form_gram_matrix(mesh, 0.05f, geom_ws); @@ -54,16 +54,38 @@ TEST_CASE("Hodge operators produce finite spectrum and circular coordinates") { auto [evals, evecs] = igneous::ops::compute_hodge_spectrum(L, G); CHECK(evals.size() > 0); + CHECK(evecs.cols() >= 2); for (int i = 0; i < evals.size(); ++i) { CHECK(std::isfinite(evals[i])); } - const auto theta = - igneous::ops::compute_circular_coordinates(mesh, evecs.col(0), 0.0f); - CHECK(theta.size() == static_cast(mesh.geometry.num_points())); - for (int i = 0; i < theta.size(); ++i) { - CHECK(std::isfinite(theta[i])); - CHECK(theta[i] >= -0.01f); - CHECK(theta[i] <= 6.30f); - } + const auto theta_0 = igneous::ops::compute_circular_coordinates( + mesh, evecs.col(0), 0.0f, 1.0f, 0); + const auto theta_1 = igneous::ops::compute_circular_coordinates( + mesh, evecs.col(1), 0.0f, 1.0f, 1); + + const auto check_theta = [&](const Eigen::VectorXf &theta) { + CHECK(theta.size() == static_cast(mesh.geometry.num_points())); + + float mean = 0.0f; + for (int i = 0; i < theta.size(); ++i) { + CHECK(std::isfinite(theta[i])); + CHECK(theta[i] >= -0.01f); + CHECK(theta[i] <= 6.30f); + mean += theta[i]; + } + mean /= static_cast(std::max(1, theta.size())); + + float var = 0.0f; + for (int i = 0; i < theta.size(); ++i) { + const float d = theta[i] - mean; + var += d * d; + } + var /= static_cast(std::max(1, theta.size())); + const float stddev = std::sqrt(var); + CHECK(stddev > 0.03f); + }; + + check_theta(theta_0); + check_theta(theta_1); }