use std::sync::Arc;
use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
use crate::atom_codes::SparseAtomCodes;
use crate::basis::{
CylinderHarmonicEvaluator, DuchonCoordinateEvaluator, EuclideanPatchEvaluator,
MobiusHarmonicEvaluator, PeriodicHarmonicEvaluator, SaeBasisSecondJet, SphereChartEvaluator,
TorusHarmonicEvaluator,
};
use crate::description_length::{BirthMdlPrescreen, predicted_birth_dl_bits};
use crate::frames::GrassmannFrame;
use crate::manifold::{
AssignmentMode, AtlasSeamKind, GraphStructureSelection, LearnedGraphAtom, OccupancyLaw,
SaeAtomBasisKind, SaeManifoldAtom, SaeManifoldRho, SaeManifoldTerm, SphereChartTransition,
UnitSpeedChartTransition, amplitude_concentration_certificate, classify_occupancy_interval,
};
use crate::migration_ledger::SaeMigrationLedger;
use crate::null_sampler::{NULL_REPLICATES, coactivation_exceedance_for_pairs};
use gam_linalg::faer_ndarray::FaerSvd;
use gam_runtime::warm_start::Fingerprinter;
use gam_solve::gaussian_reml::gaussian_reml_multi_closed_form;
use gam_solve::inference::residual_factor::{ResidualFactorInput, StructuredResidualModel};
use gam_solve::structure_search::{
ChartGlueOutcome, CollapseAction, MoveBudget, MoveProposal, SearchLedger, SearchOutcome,
StructureMove, search,
};
use gam_solve::{
AutoTopologyKind, TopologyAutoFitEvidence, TopologyAutoSelector, TopologyScoreScale,
select_topology_with_fit,
};
use gam_terms::inference::structure_evidence::{ClaimKind, StructureLedger};
use gam_terms::latent::{LatentIdMode, LatentManifold};
use gam_terms::structure::anova_atom::{
CarveReport, FissionDecision, carve, carve_input_from_fitted_atom, fission_decision,
};
use std::sync::atomic::{AtomicBool, Ordering};
const ACTIVE_SUPPORT_REL_FLOOR: f64 = 0.5;
const ARD_DIVERGENCE_LOG_PRECISION: f64 = 12.0;
const FUSION_DEPENDENCE_FLOOR: f64 = 0.6;
const NULL_EXCEEDANCE_ALPHA: f64 = 0.05;
fn null_exceedance_z_floor() -> f64 {
use statrs::distribution::{ContinuousCDF, Normal};
Normal::new(0.0, 1.0)
.expect("standard normal is well-defined")
.inverse_cdf(1.0 - NULL_EXCEEDANCE_ALPHA)
}
const ABSORPTION_ASYMMETRY_FLOOR: f64 = 0.5;
const FISSION_SYMMETRY_BREAK_EPS: f64 = 0.05;
const WITHIN_ATOM_CARVE_ALPHA: f64 = 0.05;
#[derive(Clone, Copy, Debug)]
pub struct HarvestParams {
pub max_fusions: usize,
pub max_fissions: usize,
pub max_births: usize,
}
impl Default for HarvestParams {
fn default() -> Self {
Self {
max_fusions: 4,
max_fissions: 4,
max_births: 4,
}
}
}
pub fn sparse_codes_from_term(term: &SaeManifoldTerm) -> SparseAtomCodes {
let assignments = term.assignment.assignments();
let n = assignments.nrows();
let k = assignments.ncols();
let floor = if k == 0 {
0.0
} else {
ACTIVE_SUPPORT_REL_FLOOR / k as f64
};
let mut codes = SparseAtomCodes::empty(n, k);
for row in 0..n {
for atom in 0..k {
let mass = assignments[[row, atom]];
if mass > floor {
codes.row_mut(row).assign(atom, mass);
}
}
}
codes
}
fn per_atom_max_mass(term: &SaeManifoldTerm) -> Array1<f64> {
let assignments = term.assignment.assignments();
let k = assignments.ncols();
let mut out = Array1::<f64>::zeros(k);
for atom in 0..k {
let mut max = 0.0_f64;
for &m in assignments.column(atom).iter() {
if m > max {
max = m;
}
}
out[atom] = max;
}
out
}
fn participation_ratio(spectrum: &[f64]) -> f64 {
let sum: f64 = spectrum.iter().map(|&e| e.max(0.0)).sum();
let sum_sq: f64 = spectrum.iter().map(|&e| e.max(0.0) * e.max(0.0)).sum();
if sum_sq > 0.0 {
(sum * sum) / sum_sq
} else {
0.0
}
}
fn curved_topology_for_span(span: f64) -> (usize, usize) {
match span.round().max(1.0) as usize {
0 | 1 | 2 => (1, 3), 3 => (2, 7), _ => (2, 25), }
}
fn mean_active_atoms(assignments: ArrayView2<'_, f64>) -> f64 {
let n = assignments.nrows();
let k = assignments.ncols();
if n == 0 || k == 0 {
return 1.0;
}
let floor = ACTIVE_SUPPORT_REL_FLOOR / k as f64;
let mut total_active = 0usize;
for row in 0..n {
for atom in 0..k {
if assignments[[row, atom]] > floor {
total_active += 1;
}
}
}
(total_active as f64 / n as f64).max(1.0)
}
fn per_atom_ard_divergence(rho: &SaeManifoldRho, atom: usize) -> f64 {
rho.log_ard
.get(atom)
.and_then(|axes| axes.iter().copied().reduce(f64::max))
.unwrap_or(f64::NEG_INFINITY)
}
fn post_move_structure_hash(term: &SaeManifoldTerm, mv: &StructureMove) -> u64 {
let mut fp = Fingerprinter::new();
fp.write_str("sae_structure_move");
match mv {
StructureMove::Birth { candidate } => {
fp.write_str("birth");
fp.write_usize(*candidate);
}
StructureMove::Death { atom } => {
fp.write_str("death");
fp.write_usize(*atom);
}
StructureMove::Fission { atom } => {
fp.write_str("fission");
fp.write_usize(*atom);
}
StructureMove::Fusion { a, b } => {
fp.write_str("fusion");
fp.write_usize((*a).min(*b));
fp.write_usize((*a).max(*b));
}
StructureMove::Glue { a, b, outcome } => {
fp.write_str(match outcome {
ChartGlueOutcome::Fuse => "fusion",
ChartGlueOutcome::RegisterAtlas => "atlas_register",
});
fp.write_usize((*a).min(*b));
fp.write_usize((*a).max(*b));
}
}
fp.write_usize(term.atoms.len());
for atom in &term.atoms {
fp.write_str(basis_kind_tag(&atom.basis_kind));
fp.write_usize(atom.latent_dim);
}
let digest = fp.finalize();
let bytes = digest.as_bytes();
u64::from_le_bytes([
bytes[0], bytes[1], bytes[2], bytes[3], bytes[4], bytes[5], bytes[6], bytes[7],
])
}
fn dedup_most_suspect_per_parent(candidates: Vec<(usize, f64)>) -> Vec<(usize, f64)> {
let mut best_per_parent: std::collections::HashMap<usize, f64> =
std::collections::HashMap::new();
for (atom, significance) in candidates {
best_per_parent
.entry(atom)
.and_modify(|s| {
if significance < *s {
*s = significance;
}
})
.or_insert(significance);
}
let mut out: Vec<(usize, f64)> = best_per_parent.into_iter().collect();
out.sort_by(|x, y| x.1.total_cmp(&y.1).then(x.0.cmp(&y.0)));
out
}
fn basis_kind_tag(kind: &SaeAtomBasisKind) -> &str {
match kind {
SaeAtomBasisKind::Duchon => "duchon",
SaeAtomBasisKind::Periodic => "periodic",
SaeAtomBasisKind::Sphere => "sphere",
SaeAtomBasisKind::Torus => "torus",
SaeAtomBasisKind::Linear => "linear",
SaeAtomBasisKind::EuclideanPatch => "euclidean_patch",
SaeAtomBasisKind::Poincare => "poincare",
SaeAtomBasisKind::Cylinder => "cylinder",
SaeAtomBasisKind::Mobius => "mobius",
SaeAtomBasisKind::FiniteSet => "finite_set",
SaeAtomBasisKind::Precomputed(_) => "precomputed",
}
}
fn proposal(term: &SaeManifoldTerm, mv: StructureMove, trigger: f64) -> MoveProposal {
let structure_hash = post_move_structure_hash(term, &mv);
let claim = match &mv {
StructureMove::Birth { candidate } => ClaimKind::AtomExists {
atom: term.k_atoms() + *candidate,
},
StructureMove::Death { atom } => ClaimKind::AtomExists { atom: *atom },
StructureMove::Fusion { a, b } => ClaimKind::BindingEdge { a: *a, b: *b },
StructureMove::Fission { atom } => ClaimKind::Custom {
label: format!("fission:{atom}"),
},
StructureMove::Glue { a, b, .. } => ClaimKind::Custom {
label: format!(
"seam_glue:{structure_hash:016x}:{}:{}",
(*a).min(*b),
(*a).max(*b)
),
},
};
MoveProposal {
mv,
trigger,
structure_hash,
claim,
}
}
pub fn harvest_move_proposals(
term: &SaeManifoldTerm,
rho: &SaeManifoldRho,
residuals: ArrayView2<'_, f64>,
params: &HarvestParams,
) -> Result<HarvestReport, String> {
let k = term.k_atoms();
let mut proposals: Vec<MoveProposal> = Vec::new();
let max_mass = per_atom_max_mass(term);
let terminal: std::collections::HashSet<usize> = term
.collapse_events()
.iter()
.filter(|e| matches!(e.action, CollapseAction::Terminal))
.map(|e| e.atom)
.collect();
for atom in 0..k {
let ard = per_atom_ard_divergence(rho, atom);
let diverged = ard >= ARD_DIVERGENCE_LOG_PRECISION;
let collapsed = terminal.contains(&atom);
if diverged || collapsed {
let trigger = if collapsed { f64::MAX / 2.0 } else { ard };
let trigger = trigger - max_mass[atom].min(1.0) * 1e-9;
proposals.push(proposal(term, StructureMove::Death { atom }, trigger));
}
}
let codes = sparse_codes_from_term(term);
let want_coactivation = params.max_fusions > 0 || params.max_fissions > 0;
let coactive_pairs = if want_coactivation {
codes.coactive_pair_stats()
} else {
Vec::new()
};
let coactive_pair_keys: Vec<(usize, usize)> =
coactive_pairs.iter().map(|(a, b, _)| (*a, *b)).collect();
let exceedance_z = if want_coactivation {
coactivation_exceedance_for_pairs(&codes, &coactive_pair_keys, NULL_REPLICATES)
} else {
Vec::new()
};
let z_floor = null_exceedance_z_floor();
let mut fusion_pairs: Vec<(usize, usize, f64)> = Vec::new();
for (pair_idx, &(a, b, stats)) in coactive_pairs.iter().enumerate() {
let dep = stats.dependence();
if dep < FUSION_DEPENDENCE_FLOOR {
continue;
}
let z = exceedance_z[pair_idx];
if z >= z_floor {
fusion_pairs.push((a, b, z));
}
}
fusion_pairs.sort_by(|x, y| y.2.total_cmp(&x.2).then(x.0.cmp(&y.0)).then(x.1.cmp(&y.1)));
for &(a, b, z) in fusion_pairs.iter().take(params.max_fusions) {
proposals.push(proposal(term, StructureMove::Fusion { a, b }, z));
}
let mut certified_glues = Vec::new();
let (glues_proposed, glue_candidates_screened) = harvest_glue_proposals(
term,
residuals,
params.max_fusions,
&mut proposals,
&mut certified_glues,
);
let mut fission_atoms: Vec<(usize, f64)> = Vec::new();
for (pair_idx, &(a, b, stats)) in coactive_pairs.iter().enumerate() {
let asym = stats.absorption_asymmetry();
if asym < ABSORPTION_ASYMMETRY_FLOOR {
continue;
}
let z = exceedance_z[pair_idx];
if z < z_floor {
continue;
}
let parent = if stats.p_a_given_b >= stats.p_b_given_a {
a
} else {
b
};
let significance = (1.0 - asym).max(0.0);
fission_atoms.push((parent, significance));
}
let fission_atoms = dedup_most_suspect_per_parent(fission_atoms);
let mut carve_results: Vec<FissionCarveResult> = Vec::new();
let mut fission_carve_ran_count = 0usize;
let mut fission_carve_unavailable_count = 0usize;
let mut fission_carve_blocked_count = 0usize;
let mut gated_fissions: Vec<(usize, f64)> = Vec::new();
for &(atom, significance) in fission_atoms.iter().take(params.max_fissions) {
match run_within_atom_carve(term, atom) {
Some(Ok(report)) => {
fission_carve_ran_count += 1;
let decision = fission_decision(&report, None);
let edge_p = report.edge_p_value;
let interaction = report.interaction_fraction;
carve_results.push(FissionCarveResult {
atom,
edge_p_value: edge_p,
interaction_fraction: interaction,
decision,
});
match decision {
FissionDecision::Keep => {
fission_carve_blocked_count += 1;
log::debug!(
"[structure-harvest] #993 carve KEEPS atom {atom}: binding proven \
(edge_p={edge_p:?}, interaction_fraction={interaction:.3e}); no fission proposed",
);
}
FissionDecision::SplitReconstructionOnly
| FissionDecision::SplitCertifiedJoint => {
gated_fissions.push((atom, interaction));
}
}
}
Some(Err(err)) => {
fission_carve_unavailable_count += 1;
log::debug!(
"[structure-harvest] #993 carve could not run on atom {atom}: {err}; \
fission audit rides on co-activation significance, e-gate owns acceptance",
);
gated_fissions.push((atom, significance));
}
None => {
fission_carve_unavailable_count += 1;
gated_fissions.push((atom, significance));
}
}
}
for &(atom, trigger) in &gated_fissions {
proposals.push(proposal(term, StructureMove::Fission { atom }, trigger));
}
let n = residuals.nrows();
let assignments = term.assignment.assignments();
let activity: Array1<f64> = (0..n).map(|r| assignments.row(r).sum()).collect();
let mut births_proposed = 0usize;
let mut birth_predictions: Vec<(usize, f64)> = Vec::new();
let mut births_deferred = 0usize;
let mut deferred_predicted_bits = 0.0_f64;
let mut birth_skipped_reason: Option<String> = None;
if params.max_births > 0 && n > 0 && residuals.ncols() > 0 {
let p = residuals.ncols();
let max_rank = params.max_births.min(p.saturating_sub(1));
match StructuredResidualModel::fit(ResidualFactorInput {
residuals,
activity: activity.view(),
max_factor_rank: max_rank,
}) {
Ok(model) => {
let factor = model.factor();
let diagonal = model.diagonal();
let r = model.factor_rank();
let energies: Vec<f64> = (0..r)
.map(|j| factor.column(j).iter().map(|v| v * v).sum::<f64>())
.collect();
let span = participation_ratio(&energies);
let (intrinsic_dim, basis_size) = curved_topology_for_span(span);
let g_dict = term.k_atoms();
let l0 = mean_active_atoms(assignments.view());
let n_tokens = n as f64;
let mut scored: Vec<(usize, f64)> = Vec::with_capacity(r);
for j in 0..r {
let col = factor.column(j);
let energy = energies[j];
if !(energy > 0.0) {
births_deferred += 1;
continue;
}
let norm = energy.sqrt();
let mut noise_floor = 0.0_f64;
for out in 0..p {
let u = col[out] / norm;
noise_floor += u * u * diagonal[out];
}
let mut active = 0usize;
for row in 0..n {
let res_row = residuals.row(row);
let mut proj = 0.0_f64;
for out in 0..p {
proj += res_row[out] * col[out];
}
proj /= norm;
if proj * proj > noise_floor {
active += 1;
}
}
let rho = active as f64 / n_tokens;
let predicted = predicted_birth_dl_bits(&BirthMdlPrescreen {
rho,
span,
intrinsic_dim,
basis_size,
signal_var: energy,
noise_floor,
n_tokens,
p_out: p,
g_dict,
l0,
});
if predicted.is_finite() && predicted > 0.0 {
scored.push((j, predicted));
} else {
births_deferred += 1;
if predicted.is_finite() {
deferred_predicted_bits += predicted;
}
}
}
scored.sort_by(|a, b| b.1.total_cmp(&a.1).then(a.0.cmp(&b.0)));
for &(candidate, predicted) in scored.iter().take(params.max_births) {
proposals.push(proposal(
term,
StructureMove::Birth { candidate },
predicted,
));
birth_predictions.push((candidate, predicted));
births_proposed += 1;
}
for &(_, predicted) in scored.iter().skip(params.max_births) {
births_deferred += 1;
deferred_predicted_bits += predicted;
}
if births_deferred > 0 {
log::debug!(
"[structure-harvest] #2233 MDL pre-screen deferred {births_deferred} \
birth(s) (total predicted ΔMDL {deferred_predicted_bits:.1} bits; span \
ŝ={span:.2}, d={intrinsic_dim}, m={basis_size}); proposed {births_proposed} \
ordered by predicted ΔMDL",
);
}
}
Err(e) => {
birth_skipped_reason = Some(e);
}
}
} else if params.max_births > 0 {
birth_skipped_reason =
Some("residuals empty or single-channel; no factor subspace to mine".to_string());
}
Ok(HarvestReport {
proposals,
fission_carve_results: carve_results,
fission_carve_ran_count,
fission_carve_unavailable_count,
fission_carve_blocked_count,
births_proposed,
birth_predictions,
births_deferred,
deferred_predicted_bits,
birth_skipped_reason,
glues_proposed,
glue_candidates_screened,
certified_glues,
})
}
#[derive(Clone, Debug)]
pub struct FissionCarveResult {
pub atom: usize,
pub edge_p_value: Option<f64>,
pub interaction_fraction: f64,
pub decision: FissionDecision,
}
fn run_within_atom_carve(
term: &SaeManifoldTerm,
atom: usize,
) -> Option<Result<CarveReport, String>> {
let a = &term.atoms[atom];
if a.latent_dim != 2 {
return None;
}
let evaluator = a.basis_evaluator.as_ref()?;
let (m_a, m_b) = evaluator.factor_basis_sizes()?;
let build = carve_input_from_fitted_atom(
a.basis_values.view(),
a.decoder_coefficients.view(),
m_a,
m_b,
);
let bundle = match build {
Ok(b) => b,
Err(e) => return Some(Err(e)),
};
let input = bundle.representational_carve_input();
Some(carve(&input, WITHIN_ATOM_CARVE_ALPHA))
}
#[derive(Clone, Debug)]
enum CertifiedGlueTransition {
UnitSpeed {
transition: UnitSpeedChartTransition,
rows_b: Vec<usize>,
},
Sphere(SphereChartTransition),
}
#[derive(Clone, Debug)]
struct CertifiedGlue {
a: usize,
b: usize,
outcome: ChartGlueOutcome,
transition: CertifiedGlueTransition,
}
#[derive(Clone, Debug)]
pub struct HarvestReport {
pub proposals: Vec<MoveProposal>,
pub fission_carve_results: Vec<FissionCarveResult>,
pub fission_carve_ran_count: usize,
pub fission_carve_unavailable_count: usize,
pub fission_carve_blocked_count: usize,
pub births_proposed: usize,
pub birth_predictions: Vec<(usize, f64)>,
pub births_deferred: usize,
pub deferred_predicted_bits: f64,
pub birth_skipped_reason: Option<String>,
pub glues_proposed: usize,
pub glue_candidates_screened: usize,
certified_glues: Vec<CertifiedGlue>,
}
const GLUE_DEFAULT_PERIOD: f64 = 1.0;
const GLUE_LOG_E_CLAMP: f64 = 50.0;
#[derive(Clone, Copy, Debug)]
pub struct ChartTransition {
pub sign: i8,
pub offset: f64,
pub log_e_value: f64,
}
struct SeamTransition {
sign: f64,
offset: f64,
period: f64,
rows_a: Vec<usize>,
points_a: Array2<f64>,
rows_b: Vec<usize>,
points_b: Array2<f64>,
mapped_b_to_a: Array2<f64>,
mapped_a_to_b: Array2<f64>,
}
fn atom_axis_period(term: &SaeManifoldTerm, atom: usize) -> f64 {
let coords = &term.assignment.coords;
if atom < coords.len() {
if let Some(Some(p)) = coords[atom].effective_axis_periods().first().copied() {
if p.is_finite() && p > 0.0 {
return p;
}
}
}
GLUE_DEFAULT_PERIOD
}
fn atom_active_rows(term: &SaeManifoldTerm, atom: usize) -> Vec<usize> {
let assignments = term.assignment.assignments();
let k = assignments.ncols();
let floor = if k == 0 {
0.0
} else {
ACTIVE_SUPPORT_REL_FLOOR / k as f64
};
(0..assignments.nrows())
.filter(|&r| assignments[[r, atom]] > floor)
.collect()
}
fn decoded_points_at(atom: &SaeManifoldAtom, rows: &[usize]) -> Array2<f64> {
let phi_sub = atom.basis_values.select(Axis(0), rows);
phi_sub.dot(&atom.decoder_coefficients)
}
fn periodic_decoded_points(
decoder: ArrayView2<'_, f64>,
coordinates: &[f64],
) -> Option<Array2<f64>> {
let m = decoder.nrows();
if m == 0 || m % 2 == 0 {
return None;
}
let p = decoder.ncols();
let harmonics = (m - 1) / 2;
let mut points = Array2::<f64>::zeros((coordinates.len(), p));
for (row, &coordinate) in coordinates.iter().enumerate() {
for output in 0..p {
let mut value = decoder[[0, output]];
for harmonic in 1..=harmonics {
let angle = std::f64::consts::TAU * harmonic as f64 * coordinate;
value += angle.sin() * decoder[[2 * harmonic - 1, output]]
+ angle.cos() * decoder[[2 * harmonic, output]];
}
points[[row, output]] = value;
}
}
Some(points)
}
fn periodic_decoder_under_transition(
decoder_a: ArrayView2<'_, f64>,
sign: i8,
offset: f64,
) -> Option<Array2<f64>> {
if !matches!(sign, -1 | 1) || decoder_a.nrows() == 0 || decoder_a.nrows() % 2 == 0 {
return None;
}
let mut mapped = decoder_a.to_owned();
let harmonics = (decoder_a.nrows() - 1) / 2;
for harmonic in 1..=harmonics {
let angle = std::f64::consts::TAU * harmonic as f64 * offset;
let (cosine, sine) = (angle.cos(), angle.sin());
for output in 0..decoder_a.ncols() {
let a_sin = decoder_a[[2 * harmonic - 1, output]];
let a_cos = decoder_a[[2 * harmonic, output]];
if sign == 1 {
mapped[[2 * harmonic - 1, output]] = cosine * a_sin - sine * a_cos;
mapped[[2 * harmonic, output]] = sine * a_sin + cosine * a_cos;
} else {
mapped[[2 * harmonic - 1, output]] = -cosine * a_sin + sine * a_cos;
mapped[[2 * harmonic, output]] = sine * a_sin + cosine * a_cos;
}
}
}
Some(mapped)
}
fn fit_periodic_transition_from_decoders(
decoder_a: ArrayView2<'_, f64>,
decoder_b: ArrayView2<'_, f64>,
) -> Option<(i8, f64)> {
if decoder_a.dim() != decoder_b.dim() || decoder_a.nrows() < 3 || decoder_a.nrows() % 2 == 0 {
return None;
}
let dot = |left_row: usize, right_row: usize| -> f64 {
(0..decoder_a.ncols())
.map(|output| decoder_b[[left_row, output]] * decoder_a[[right_row, output]])
.sum()
};
let harmonics = (decoder_a.nrows() - 1) / 2;
let mut candidates = Vec::new();
for sign in [1_i8, -1_i8] {
for harmonic in 1..=harmonics {
let sin_row = 2 * harmonic - 1;
let cos_row = 2 * harmonic;
let (cos_score, sin_score) = if sign == 1 {
(
dot(sin_row, sin_row) + dot(cos_row, cos_row),
-dot(sin_row, cos_row) + dot(cos_row, sin_row),
)
} else {
(
-dot(sin_row, sin_row) + dot(cos_row, cos_row),
dot(sin_row, cos_row) + dot(cos_row, sin_row),
)
};
if cos_score.hypot(sin_score) > 0.0 {
let harmonic_phase = sin_score.atan2(cos_score).rem_euclid(std::f64::consts::TAU);
for branch in 0..harmonic {
let phase =
(harmonic_phase + std::f64::consts::TAU * branch as f64) / harmonic as f64;
candidates.push((sign, phase));
}
break;
}
}
}
candidates
.into_iter()
.filter_map(|(sign, angle)| {
let offset = angle / std::f64::consts::TAU;
let mapped = periodic_decoder_under_transition(decoder_a, sign, offset)?;
let residual = mapped
.iter()
.zip(decoder_b.iter())
.map(|(predicted, observed)| (predicted - observed).powi(2))
.sum::<f64>();
residual.is_finite().then_some((sign, offset, residual))
})
.min_by(|left, right| {
left.2.total_cmp(&right.2).then_with(|| {
right.0.cmp(&left.0)
})
})
.map(|(sign, offset, _)| (sign, offset))
}
fn fit_seam_transition(term: &SaeManifoldTerm, a: usize, b: usize) -> Option<SeamTransition> {
let k = term.k_atoms();
if a >= k || b >= k || a == b {
return None;
}
let atom_a = &term.atoms[a];
let atom_b = &term.atoms[b];
if atom_a.latent_dim != 1
|| atom_b.latent_dim != 1
|| !matches!(atom_a.basis_kind, SaeAtomBasisKind::Periodic)
|| !matches!(atom_b.basis_kind, SaeAtomBasisKind::Periodic)
{
return None;
}
let p = atom_a.decoder_coefficients.ncols();
if p == 0 || p != atom_b.decoder_coefficients.ncols() {
return None;
}
let rows_a = atom_active_rows(term, a);
let rows_b = atom_active_rows(term, b);
if rows_a.is_empty() || rows_b.is_empty() {
return None;
}
let coords = &term.assignment.coords;
if b >= coords.len() || coords[b].latent_dim() < 1 {
return None;
}
let period_a = atom_axis_period(term, a);
let period_b = atom_axis_period(term, b);
if period_a.to_bits() != 1.0_f64.to_bits() || period_b.to_bits() != period_a.to_bits() {
return None;
}
let decoder_a = atom_a.full_width_decoder();
let decoder_b = atom_b.full_width_decoder();
let (sign, offset) = fit_periodic_transition_from_decoders(decoder_a.view(), decoder_b.view())?;
let points_a = decoded_points_at(atom_a, &rows_a);
let points_b = decoded_points_at(atom_b, &rows_b);
let mapped_b_coords: Vec<f64> = rows_b
.iter()
.map(|&row| (sign as f64 * coords[b].row(row)[0] + offset).rem_euclid(period_a))
.collect();
let mapped_a_coords: Vec<f64> = rows_a
.iter()
.map(|&row| {
(sign as f64 * (coords[a].row(row)[0] - offset)).rem_euclid(period_a)
})
.collect();
let mapped_b_to_a = periodic_decoded_points(decoder_a.view(), &mapped_b_coords)?;
let mapped_a_to_b = periodic_decoded_points(decoder_b.view(), &mapped_a_coords)?;
Some(SeamTransition {
sign: sign as f64,
offset,
period: period_a,
rows_a,
points_a,
rows_b,
points_b,
mapped_b_to_a,
mapped_a_to_b,
})
}
fn unit_speed_glue_certificate(
term: &SaeManifoldTerm,
residuals: ArrayView2<'_, f64>,
a: usize,
b: usize,
) -> Option<(ChartTransition, CertifiedGlue)> {
let seam = fit_seam_transition(term, a, b)?;
let log_e = seam_equivalence_log_e(
residuals,
&seam.rows_a,
&seam.points_a,
&seam.mapped_a_to_b,
&seam.rows_b,
&seam.points_b,
&seam.mapped_b_to_a,
)?;
let chart_transition = ChartTransition {
sign: seam.sign as i8,
offset: seam.offset,
log_e_value: log_e,
};
let outcome = if chart_transition.sign == 1 {
ChartGlueOutcome::Fuse
} else {
ChartGlueOutcome::RegisterAtlas
};
let transition = UnitSpeedChartTransition::new(
b,
a,
chart_transition.sign,
chart_transition.offset,
seam.period,
AtlasSeamKind::Regular,
)
.ok()?;
Some((
chart_transition,
CertifiedGlue {
a,
b,
outcome,
transition: CertifiedGlueTransition::UnitSpeed {
transition,
rows_b: seam.rows_b,
},
},
))
}
fn seam_equivalence_log_e(
residuals: ArrayView2<'_, f64>,
rows_a: &[usize],
points_a: &Array2<f64>,
mapped_a_to_b: &Array2<f64>,
rows_b: &[usize],
points_b: &Array2<f64>,
mapped_b_to_a: &Array2<f64>,
) -> Option<f64> {
let p = points_a.ncols();
if p == 0 || residuals.ncols() != p || points_b.ncols() != p {
return None;
}
let na = points_a.nrows();
let nb = points_b.nrows();
if na != rows_a.len() || nb != rows_b.len() {
return None;
}
if mapped_a_to_b.dim() != (na, p) || mapped_b_to_a.dim() != (nb, p) {
return None;
}
if na < 2 || nb < 2 {
return None;
}
let a_est: Vec<usize> = (0..na).filter(|i| i % 2 == 0).collect();
let a_eval: Vec<usize> = (0..na).filter(|i| i % 2 == 1).collect();
let b_est: Vec<usize> = (0..nb).filter(|i| i % 2 == 0).collect();
let b_eval: Vec<usize> = (0..nb).filter(|i| i % 2 == 1).collect();
let n_est = a_est.len() + b_est.len();
let n_eval = a_eval.len() + b_eval.len();
if n_est == 0 || n_eval == 0 {
return None;
}
let mut mu = vec![0.0_f64; p];
for &i in &a_est {
for c in 0..p {
mu[c] += points_a[[i, c]];
}
}
for &i in &b_est {
for c in 0..p {
mu[c] += points_b[[i, c]];
}
}
for c in 0..p {
mu[c] /= n_est as f64;
}
let point_null_sq =
|pt: ArrayView1<'_, f64>| -> f64 { (0..p).map(|c| (pt[c] - mu[c]).powi(2)).sum::<f64>() };
let mut pool_acc = 0.0_f64;
for &i in &a_est {
pool_acc += point_null_sq(points_a.row(i));
}
for &i in &b_est {
pool_acc += point_null_sq(points_b.row(i));
}
let pool_sq = pool_acc / (n_est as f64 * p as f64);
if !(pool_sq.is_finite() && pool_sq > 0.0) {
return None;
}
let mut band_acc = 0.0_f64;
let mut band_rows = 0usize;
for &i in &a_est {
let r = rows_a[i];
for c in 0..p {
band_acc += residuals[[r, c]].powi(2);
}
band_rows += 1;
}
for &i in &b_est {
let r = rows_b[i];
for c in 0..p {
band_acc += residuals[[r, c]].powi(2);
}
band_rows += 1;
}
let band_raw = if band_rows == 0 {
0.0
} else {
band_acc / (band_rows as f64 * p as f64)
};
let band_sq = band_raw.max(pool_sq * f64::EPSILON);
if !(pool_sq > band_sq) {
return None;
}
let norm_term = (p as f64 / 2.0) * (pool_sq / band_sq).ln();
let mut log_e = 0.0_f64;
for &i in &b_eval {
let e_glue: f64 = (0..p)
.map(|c| (points_b[[i, c]] - mapped_b_to_a[[i, c]]).powi(2))
.sum();
let e_null = point_null_sq(points_b.row(i));
log_e += norm_term - e_glue / (2.0 * band_sq) + e_null / (2.0 * pool_sq);
}
for &i in &a_eval {
let e_glue: f64 = (0..p)
.map(|c| (points_a[[i, c]] - mapped_a_to_b[[i, c]]).powi(2))
.sum();
let e_null = point_null_sq(points_a.row(i));
log_e += norm_term - e_glue / (2.0 * band_sq) + e_null / (2.0 * pool_sq);
}
if !log_e.is_finite() {
log_e = if log_e < 0.0 {
-GLUE_LOG_E_CLAMP
} else {
GLUE_LOG_E_CLAMP
};
}
Some(log_e.clamp(-GLUE_LOG_E_CLAMP, GLUE_LOG_E_CLAMP))
}
fn sphere_linear_block(atom: &SaeManifoldAtom) -> Option<Array2<f64>> {
let decoder = atom.full_width_decoder();
if decoder.nrows() != 7 {
return None;
}
Some(decoder.slice(ndarray::s![1..4, ..]).to_owned())
}
fn inverse_3x3(m: &[[f64; 3]; 3]) -> Option<[[f64; 3]; 3]> {
let det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
- m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
+ m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
if !det.is_finite() || det.abs() < 1e-12 {
return None;
}
let inv_det = 1.0 / det;
let mut inv = [[0.0; 3]; 3];
inv[0][0] = (m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inv_det;
inv[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * inv_det;
inv[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inv_det;
inv[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * inv_det;
inv[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inv_det;
inv[1][2] = (m[0][2] * m[1][0] - m[0][0] * m[1][2]) * inv_det;
inv[2][0] = (m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inv_det;
inv[2][1] = (m[0][1] * m[2][0] - m[0][0] * m[2][1]) * inv_det;
inv[2][2] = (m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inv_det;
Some(inv)
}
fn nearest_orthogonal_3x3(m: [[f64; 3]; 3]) -> Option<[[f64; 3]; 3]> {
let mut q = m;
for _ in 0..128 {
let inv = inverse_3x3(&q)?;
let mut next = [[0.0; 3]; 3];
let mut diff = 0.0;
for i in 0..3 {
for j in 0..3 {
next[i][j] = 0.5 * (q[i][j] + inv[j][i]);
diff += (next[i][j] - q[i][j]).abs();
}
}
q = next;
if diff < 1e-15 {
break;
}
}
if q.iter().flatten().any(|x| !x.is_finite()) {
return None;
}
Some(q)
}
fn sphere_decoded_points_at_units(
decoder: ArrayView2<'_, f64>,
units: &[[f64; 3]],
) -> Option<Array2<f64>> {
if decoder.nrows() != 7 {
return None;
}
let p = decoder.ncols();
let mut points = Array2::<f64>::zeros((units.len(), p));
for (row, &[x, y, z]) in units.iter().enumerate() {
let phi = [1.0, x, y, z, x * y, y * z, x * z];
for output in 0..p {
let mut value = 0.0;
for (basis, &phi_b) in phi.iter().enumerate() {
value += phi_b * decoder[[basis, output]];
}
points[[row, output]] = value;
}
}
Some(points)
}
struct SphereSeamTransition {
rotation: [[f64; 3]; 3],
seam_kind: AtlasSeamKind,
rows_a: Vec<usize>,
points_a: Array2<f64>,
rows_b: Vec<usize>,
points_b: Array2<f64>,
mapped_b_to_a: Array2<f64>,
mapped_a_to_b: Array2<f64>,
}
fn sphere_row_unit(atom: &SaeManifoldAtom, row: usize) -> [f64; 3] {
[
atom.basis_values[[row, 1]],
atom.basis_values[[row, 2]],
atom.basis_values[[row, 3]],
]
}
fn rotate_unit(r: &[[f64; 3]; 3], u: [f64; 3]) -> [f64; 3] {
[
r[0][0] * u[0] + r[0][1] * u[1] + r[0][2] * u[2],
r[1][0] * u[0] + r[1][1] * u[1] + r[1][2] * u[2],
r[2][0] * u[0] + r[2][1] * u[1] + r[2][2] * u[2],
]
}
fn is_sphere_pair(term: &SaeManifoldTerm, a: usize, b: usize) -> bool {
let k = term.k_atoms();
if a >= k || b >= k || a == b {
return false;
}
let sa = &term.atoms[a];
let sb = &term.atoms[b];
sa.latent_dim == 2
&& sb.latent_dim == 2
&& matches!(sa.basis_kind, SaeAtomBasisKind::Sphere)
&& matches!(sb.basis_kind, SaeAtomBasisKind::Sphere)
}
fn fit_sphere_seam_transition(
term: &SaeManifoldTerm,
a: usize,
b: usize,
) -> Option<SphereSeamTransition> {
if !is_sphere_pair(term, a, b) {
return None;
}
let atom_a = &term.atoms[a];
let atom_b = &term.atoms[b];
let p = atom_a.decoder_coefficients.ncols();
if p == 0 || p != atom_b.decoder_coefficients.ncols() {
return None;
}
let l_a = sphere_linear_block(atom_a)?;
let l_b = sphere_linear_block(atom_b)?;
let m_prod = l_a.dot(&l_b.t());
let mut m = [[0.0_f64; 3]; 3];
for i in 0..3 {
for j in 0..3 {
m[i][j] = m_prod[[i, j]];
}
}
let rotation = nearest_orthogonal_3x3(m)?;
let rows_a = atom_active_rows(term, a);
let rows_b = atom_active_rows(term, b);
if rows_a.is_empty() || rows_b.is_empty() {
return None;
}
let lat_of = |u: [f64; 3]| -> f64 { u[2].clamp(-1.0, 1.0).asin() };
let (mut a_lat_lo, mut a_lat_hi) = (f64::INFINITY, f64::NEG_INFINITY);
for &r in &rows_a {
let lat = lat_of(sphere_row_unit(atom_a, r));
a_lat_lo = a_lat_lo.min(lat);
a_lat_hi = a_lat_hi.max(lat);
}
let (mut b_lat_lo, mut b_lat_hi) = (f64::INFINITY, f64::NEG_INFINITY);
for &r in &rows_b {
let lat = lat_of(sphere_row_unit(atom_b, r));
b_lat_lo = b_lat_lo.min(lat);
b_lat_hi = b_lat_hi.max(lat);
}
let b_pole_in_a = lat_of(rotate_unit(&rotation, [0.0, 0.0, 1.0]));
let inv_rotation = inverse_3x3(&rotation)?;
let a_pole_in_b = lat_of(rotate_unit(&inv_rotation, [0.0, 0.0, 1.0]));
let b_pole_interior_to_a = b_pole_in_a > a_lat_lo && b_pole_in_a < a_lat_hi;
let a_pole_interior_to_b = a_pole_in_b > b_lat_lo && a_pole_in_b < b_lat_hi;
let seam_kind = if b_pole_interior_to_a && a_pole_interior_to_b {
AtlasSeamKind::Pole
} else {
AtlasSeamKind::Regular
};
let units_a: Vec<[f64; 3]> = rows_a.iter().map(|&r| sphere_row_unit(atom_a, r)).collect();
let units_b: Vec<[f64; 3]> = rows_b.iter().map(|&r| sphere_row_unit(atom_b, r)).collect();
let points_a = sphere_decoded_points_at_units(atom_a.full_width_decoder().view(), &units_a)?;
let points_b = sphere_decoded_points_at_units(atom_b.full_width_decoder().view(), &units_b)?;
let mapped_b_units: Vec<[f64; 3]> =
units_b.iter().map(|&u| rotate_unit(&rotation, u)).collect();
let mapped_b_to_a =
sphere_decoded_points_at_units(atom_a.full_width_decoder().view(), &mapped_b_units)?;
let mapped_a_units: Vec<[f64; 3]> = units_a
.iter()
.map(|&u| rotate_unit(&inv_rotation, u))
.collect();
let mapped_a_to_b =
sphere_decoded_points_at_units(atom_b.full_width_decoder().view(), &mapped_a_units)?;
Some(SphereSeamTransition {
rotation,
seam_kind,
rows_a,
points_a,
rows_b,
points_b,
mapped_b_to_a,
mapped_a_to_b,
})
}
fn sphere_glue_pair_evalue(
term: &SaeManifoldTerm,
residuals: ArrayView2<'_, f64>,
a: usize,
b: usize,
) -> Option<(SphereChartTransition, f64)> {
let seam = fit_sphere_seam_transition(term, a, b)?;
if !matches!(seam.seam_kind, AtlasSeamKind::Pole) {
return None;
}
let log_e = seam_equivalence_log_e(
residuals,
&seam.rows_a,
&seam.points_a,
&seam.mapped_a_to_b,
&seam.rows_b,
&seam.points_b,
&seam.mapped_b_to_a,
)?;
let transition = SphereChartTransition::new(b, a, seam.rotation, AtlasSeamKind::Pole).ok()?;
Some((transition, log_e))
}
fn harvest_glue_proposals(
term: &SaeManifoldTerm,
residuals: ArrayView2<'_, f64>,
budget: usize,
proposals: &mut Vec<MoveProposal>,
certified_glues: &mut Vec<CertifiedGlue>,
) -> (usize, usize) {
let k = term.k_atoms();
if k < 2 || budget == 0 {
return (0, 0);
}
let assignments = term.assignment.assignments();
let n_rows = assignments.nrows();
if n_rows == 0 {
return (0, 0);
}
let floor = ACTIVE_SUPPORT_REL_FLOOR / k as f64;
let supports: Vec<Vec<bool>> = (0..k)
.map(|atom| {
(0..n_rows)
.map(|r| assignments[[r, atom]] > floor)
.collect()
})
.collect();
let support_sizes: Vec<usize> = supports
.iter()
.map(|s| s.iter().filter(|&&x| x).count())
.collect();
let frames: Vec<Option<GrassmannFrame>> = (0..k)
.map(|atom| {
let at = &term.atoms[atom];
if at.latent_dim == 1 && matches!(at.basis_kind, SaeAtomBasisKind::Periodic) {
GrassmannFrame::from_decoder_row_space(at.decoder_coefficients.view())
} else {
None
}
})
.collect();
let mut screened = 0usize;
let mut candidates: Vec<(usize, usize, f64)> = Vec::new();
for a in 0..k {
let fa = match &frames[a] {
Some(f) => f,
None => continue,
};
if support_sizes[a] == 0 {
continue;
}
for b in (a + 1)..k {
if term.charts_share_atlas(a, b) {
continue;
}
let fb = match &frames[b] {
Some(f) => f,
None => continue,
};
if support_sizes[b] == 0 {
continue;
}
let inter = (0..n_rows)
.filter(|&r| supports[a][r] && supports[b][r])
.count();
let expected = support_sizes[a] as f64 * support_sizes[b] as f64 / n_rows as f64;
if inter as f64 > expected {
continue;
}
let alignment = match fa.max_principal_angle(fb.frame()) {
Ok(theta) => theta.cos(),
Err(_) => continue,
};
if !alignment.is_finite() {
continue;
}
screened += 1;
candidates.push((a, b, alignment));
}
}
candidates.sort_by(|x, y| y.2.total_cmp(&x.2).then(x.0.cmp(&y.0)).then(x.1.cmp(&y.1)));
let mut proposed = 0usize;
for &(a, b, _score) in candidates.iter().take(budget) {
if let Some((tr, certificate)) = unit_speed_glue_certificate(term, residuals, a, b) {
proposals.push(proposal(
term,
StructureMove::Glue {
a,
b,
outcome: certificate.outcome,
},
tr.log_e_value,
));
certified_glues.push(certificate);
proposed += 1;
}
}
let mut sphere_candidates: Vec<(usize, usize)> = Vec::new();
for a in 0..k {
if support_sizes[a] == 0 || !matches!(term.atoms[a].basis_kind, SaeAtomBasisKind::Sphere) {
continue;
}
for b in (a + 1)..k {
if support_sizes[b] == 0
|| !matches!(term.atoms[b].basis_kind, SaeAtomBasisKind::Sphere)
|| term.charts_share_atlas(a, b)
{
continue;
}
let inter = (0..n_rows)
.filter(|&r| supports[a][r] && supports[b][r])
.count();
let expected = support_sizes[a] as f64 * support_sizes[b] as f64 / n_rows as f64;
if inter as f64 > expected {
continue;
}
sphere_candidates.push((a, b));
}
}
for &(a, b) in sphere_candidates.iter().take(budget) {
screened += 1;
if let Some((transition, log_e)) = sphere_glue_pair_evalue(term, residuals, a, b) {
proposals.push(proposal(
term,
StructureMove::Glue {
a,
b,
outcome: ChartGlueOutcome::RegisterAtlas,
},
log_e,
));
certified_glues.push(CertifiedGlue {
a,
b,
outcome: ChartGlueOutcome::RegisterAtlas,
transition: CertifiedGlueTransition::Sphere(transition),
});
proposed += 1;
}
}
(proposed, screened)
}
fn transplant_glued_coords(
term: &mut SaeManifoldTerm,
a: usize,
b: usize,
transition: &UnitSpeedChartTransition,
rows_b: &[usize],
) -> Result<(), String> {
if transition.from_chart != b || transition.to_chart != a {
return Err(format!(
"transplant_glued_coords: transition {}->{} does not match glue ({a},{b})",
transition.from_chart, transition.to_chart
));
}
let coords = &mut term.assignment.coords;
if a >= coords.len() || b >= coords.len() {
return Err(format!(
"transplant_glued_coords: glue ({a},{b}) outside {} coordinate blocks",
coords.len()
));
}
let da = coords[a].latent_dim();
let db = coords[b].latent_dim();
if da < 1 || db < 1 || coords[a].n_obs() != coords[b].n_obs() {
return Err(format!(
"transplant_glued_coords: incompatible coordinate blocks for glue ({a},{b})"
));
}
let flat_b = coords[b].as_flat().to_owned();
let mut flat_a = coords[a].as_flat().to_owned();
let n = coords[b].n_obs();
for &r in rows_b {
if r >= n {
return Err(format!(
"transplant_glued_coords: certified row {r} outside n={n} for glue ({a},{b})"
));
}
let t_b = flat_b[r * db];
flat_a[r * da] = transition.apply(t_b);
}
coords[a].set_flat(flat_a.view());
Ok(())
}
pub fn apply_structure_move(
term: &SaeManifoldTerm,
rho: &SaeManifoldRho,
mv: &StructureMove,
birth_decoders: &[Array2<f64>],
) -> Result<(SaeManifoldTerm, SaeManifoldRho), String> {
match mv {
StructureMove::Death { atom } => {
let mut child = term.clone();
demote_atom(&mut child, *atom)?;
Ok((child, rho.clone()))
}
StructureMove::Fusion { a, b } => {
let mut child = term.clone();
fold_atom_into(&mut child, *a, *b)?;
Ok((child, rho.clone()))
}
StructureMove::Glue { a, b, outcome } => {
if is_sphere_pair(term, *a, *b) {
let mut child = term.clone();
match outcome {
ChartGlueOutcome::RegisterAtlas => {
let seam = fit_sphere_seam_transition(term, *a, *b).ok_or_else(|| {
format!(
"apply_structure_move: sphere seam ({a},{b}) is no longer identifiable"
)
})?;
if !matches!(seam.seam_kind, AtlasSeamKind::Pole) {
return Err(format!(
"apply_structure_move: sphere seam ({a},{b}) is not a pole seam"
));
}
child.register_sphere_chart_transition(SphereChartTransition::new(
*b,
*a,
seam.rotation,
AtlasSeamKind::Pole,
)?)?;
}
ChartGlueOutcome::Fuse => {
return Err(format!(
"apply_structure_move: sphere pole seam ({a},{b}) cannot be destructively fused"
));
}
}
return Ok((child, rho.clone()));
}
let seam = fit_seam_transition(term, *a, *b).ok_or_else(|| {
format!("apply_structure_move: chart seam ({a},{b}) is no longer identifiable")
})?;
let mut child = term.clone();
match outcome {
ChartGlueOutcome::Fuse => {
if seam.sign != 1.0 {
return Err(format!(
"apply_structure_move: refusing to fuse orientation-reversing seam ({a},{b})"
));
}
let transition = UnitSpeedChartTransition::new(
*b,
*a,
1,
seam.offset,
seam.period,
AtlasSeamKind::Regular,
)?;
fold_atom_into(&mut child, *a, *b)?;
transplant_glued_coords(&mut child, *a, *b, &transition, &seam.rows_b)?;
}
ChartGlueOutcome::RegisterAtlas => {
if seam.sign != -1.0 {
return Err(format!(
"apply_structure_move: atlas registration requires an orientation-reversing seam, got sign {}",
seam.sign
));
}
child.register_chart_transition(UnitSpeedChartTransition::new(
*b,
*a,
-1,
seam.offset,
seam.period,
AtlasSeamKind::Regular,
)?)?;
}
}
Ok((child, rho.clone()))
}
StructureMove::Fission { atom } => {
let (child, child_rho) = duplicate_atom(term, rho, *atom)?;
Ok((child, child_rho))
}
StructureMove::Birth { candidate } => {
let decoder = birth_decoders.get(*candidate).ok_or_else(|| {
format!(
"apply_structure_move: birth candidate {candidate} out of range \
({} residual-factor decoders)",
birth_decoders.len()
)
})?;
born_atom(term, rho, decoder.view())
}
}
}
#[derive(Clone, Debug)]
pub enum BirthSeed {
ResidualFactor(Array2<f64>),
Circle {
decoder: Array2<f64>,
phase_coords: Array2<f64>,
gate: Vec<f64>,
},
}
pub fn apply_structure_move_seeded(
term: &SaeManifoldTerm,
rho: &SaeManifoldRho,
mv: &StructureMove,
birth_seeds: &[BirthSeed],
) -> Result<(SaeManifoldTerm, SaeManifoldRho), String> {
match mv {
StructureMove::Birth { candidate } => {
let seed = birth_seeds.get(*candidate).ok_or_else(|| {
format!(
"apply_structure_move_seeded: birth candidate {candidate} out of range \
({} birth seeds)",
birth_seeds.len()
)
})?;
match seed {
BirthSeed::ResidualFactor(decoder) => born_atom(term, rho, decoder.view()),
BirthSeed::Circle {
decoder,
phase_coords,
gate,
} => born_circle_atom(
term,
rho,
decoder.clone(),
phase_coords.clone(),
gate.clone(),
),
}
}
other => apply_structure_move(term, rho, other, &[]),
}
}
const DEMOTE_LOGIT: f64 = -40.0;
fn demote_atom(term: &mut SaeManifoldTerm, atom: usize) -> Result<(), String> {
let k = term.k_atoms();
if atom >= k {
return Err(format!("demote_atom: atom {atom} out of range (K={k})"));
}
for row in 0..term.assignment.logits.nrows() {
term.assignment.logits[[row, atom]] = DEMOTE_LOGIT;
}
Ok(())
}
fn fold_atom_into(term: &mut SaeManifoldTerm, a: usize, b: usize) -> Result<(), String> {
let k = term.k_atoms();
if a >= k || b >= k {
return Err(format!(
"fold_atom_into: atoms ({a},{b}) out of range (K={k})"
));
}
if a == b {
return Err("fold_atom_into: cannot fuse an atom with itself".to_string());
}
let softmax_routing = matches!(term.assignment.mode, AssignmentMode::Softmax { .. });
for row in 0..term.assignment.logits.nrows() {
let la = term.assignment.logits[[row, a]];
let lb = term.assignment.logits[[row, b]];
term.assignment.logits[[row, a]] = if softmax_routing {
let m = la.max(lb);
if m == f64::NEG_INFINITY {
f64::NEG_INFINITY
} else {
m + ((la - m).exp() + (lb - m).exp()).ln()
}
} else {
la.max(lb)
};
}
demote_atom(term, b)?;
Ok(())
}
fn remove_atoms(
term: &mut SaeManifoldTerm,
rho: &mut SaeManifoldRho,
remove: &std::collections::BTreeSet<usize>,
) -> Result<(), String> {
let k = term.k_atoms();
if let Some(&bad) = remove.iter().find(|&&j| j >= k) {
return Err(format!("remove_atoms: atom {bad} out of range (K={k})"));
}
if remove.len() >= k {
return Err("remove_atoms: cannot remove every atom".to_string());
}
if remove.is_empty() {
return Ok(());
}
let n = term.assignment.logits.nrows();
if term.assignment.logits.ncols() != k
|| term.assignment.coords.len() != k
|| term.assignment.ungated.len() != k
{
return Err(format!(
"remove_atoms: atom-indexed assignment shape mismatch: atoms={k}, \
logits={:?}, coords={}, ungated={}",
term.assignment.logits.dim(),
term.assignment.coords.len(),
term.assignment.ungated.len()
));
}
if let Some(frozen) = term.assignment.frozen_logits.as_ref() {
if frozen.dim() != (n, k) {
return Err(format!(
"remove_atoms: frozen logits shape {:?} must equal ({n}, {k})",
frozen.dim()
));
}
}
if rho.log_lambda_smooth.len() != k || rho.log_ard.len() != k {
return Err(format!(
"remove_atoms: rho per-atom lengths (smooth {}, ard {}) must equal K={k}",
rho.log_lambda_smooth.len(),
rho.log_ard.len()
));
}
let keep: Vec<usize> = (0..k).filter(|j| !remove.contains(j)).collect();
let mut old_to_new = vec![None; k];
for (new, &old) in keep.iter().enumerate() {
old_to_new[old] = Some(new);
}
term.remap_chart_atlases(&old_to_new)?;
term.atoms = keep.iter().map(|&j| term.atoms[j].clone()).collect();
let compacted_logits = Array2::from_shape_fn((n, keep.len()), |(row, new_atom)| {
term.assignment.logits[[row, keep[new_atom]]]
});
term.assignment.logits = compacted_logits;
term.assignment.coords = keep
.iter()
.map(|&j| term.assignment.coords[j].clone())
.collect();
term.assignment.ungated = keep.iter().map(|&j| term.assignment.ungated[j]).collect();
term.assignment.frozen_logits = None;
rho.log_lambda_smooth = keep.iter().map(|&j| rho.log_lambda_smooth[j]).collect();
rho.log_ard = keep.iter().map(|&j| rho.log_ard[j].clone()).collect();
term.collapse_events.clear();
term.last_row_layout = None;
term.last_frames_active = false;
term.fixed_decoder_assembly = false;
term.border_hbb_workspace = Array2::<f64>::zeros((0, 0));
term.decoder_repulsion_gate = None;
term.barrier_coactivation_gate = None;
term.streaming_gates_frozen = false;
term.curvature_walk_report = None;
term.expected_evidence_gauge_deflated_directions = None;
term.evidence_gauge_deflation_reanchors = 0;
term.evidence_gauge_deflation_last_delta_sign = 0;
term.dictionary_cocollapse_reseeds = 0;
term.structural_cocollapse_reseeds = 0;
term.atom_inner_fits = None;
term.oos_linear_images = None;
term.hybrid_split_report = None;
term.best_cocollapse_incumbent = None;
term.best_fit_incumbent = None;
let softmax_active_cap = term.softmax_active_cap;
term.set_softmax_active_cap(softmax_active_cap);
Ok(())
}
fn compact_glued_atoms(
term: &mut SaeManifoldTerm,
rho: &mut SaeManifoldRho,
round_ledger: &SearchLedger,
certified_glues: &[CertifiedGlue],
) -> Result<usize, String> {
use gam_solve::structure_search::MoveVerdict;
let accepted_glues: Vec<(usize, usize, ChartGlueOutcome)> = round_ledger
.moves
.iter()
.filter_map(|rec| {
if let (StructureMove::Glue { a, b, outcome }, MoveVerdict::Accepted { .. }) =
(&rec.mv, &rec.verdict)
{
Some((*a, *b, *outcome))
} else {
None
}
})
.collect();
if accepted_glues.is_empty() {
return Ok(0);
}
let k = term.k_atoms();
let mut touched = std::collections::BTreeSet::new();
for &(a, b, _) in &accepted_glues {
if a >= k || b >= k || a == b {
return Err(format!(
"compact_glued_atoms: accepted glue ({a},{b}) out of range or self-gluing (K={k})"
));
}
if !touched.insert(a) || !touched.insert(b) {
return Err(format!(
"compact_glued_atoms: accepted glues are not an atom-disjoint matching; \
atom reused by ({a},{b})"
));
}
}
let mut adopted: Vec<CertifiedGlue> = Vec::with_capacity(accepted_glues.len());
for &(a, b, outcome) in &accepted_glues {
let mut matches = certified_glues
.iter()
.filter(|certificate| certificate.a == a && certificate.b == b);
let certificate = matches.next().ok_or_else(|| {
format!("compact_glued_atoms: accepted glue ({a},{b}) has no harvest-time certificate")
})?;
if matches.next().is_some() {
return Err(format!(
"compact_glued_atoms: accepted glue ({a},{b}) has duplicate harvest-time certificates"
));
}
if certificate.outcome != outcome {
return Err(format!(
"compact_glued_atoms: accepted glue ({a},{b}) outcome {outcome:?} does not match certified {:?}",
certificate.outcome
));
}
match (&certificate.transition, outcome) {
(CertifiedGlueTransition::UnitSpeed { transition, .. }, ChartGlueOutcome::Fuse)
if transition.from_chart == b
&& transition.to_chart == a
&& transition.sign == 1
&& matches!(transition.seam_kind, AtlasSeamKind::Regular) => {}
(
CertifiedGlueTransition::UnitSpeed { transition, .. },
ChartGlueOutcome::RegisterAtlas,
) if transition.from_chart == b
&& transition.to_chart == a
&& transition.sign == -1
&& matches!(transition.seam_kind, AtlasSeamKind::Regular) => {}
(CertifiedGlueTransition::Sphere(transition), ChartGlueOutcome::RegisterAtlas)
if transition.from_chart == b
&& transition.to_chart == a
&& matches!(transition.seam_kind, AtlasSeamKind::Pole) => {}
_ => {
return Err(format!(
"compact_glued_atoms: accepted glue ({a},{b}) is incompatible with its certified transition"
));
}
}
adopted.push(certificate.clone());
}
let mut child_term = term.clone();
let mut child_rho = rho.clone();
let mut to_remove: std::collections::BTreeSet<usize> = std::collections::BTreeSet::new();
for certificate in adopted {
let (a, b) = (certificate.a, certificate.b);
match certificate.transition {
CertifiedGlueTransition::UnitSpeed { transition, rows_b }
if matches!(certificate.outcome, ChartGlueOutcome::Fuse) =>
{
fold_atom_into(&mut child_term, a, b)?;
transplant_glued_coords(&mut child_term, a, b, &transition, &rows_b)?;
to_remove.insert(b);
}
CertifiedGlueTransition::UnitSpeed { transition, .. } => {
child_term.register_chart_transition(transition)?;
}
CertifiedGlueTransition::Sphere(transition) => {
child_term.register_sphere_chart_transition(transition)?;
}
}
}
remove_atoms(&mut child_term, &mut child_rho, &to_remove)?;
*term = child_term;
*rho = child_rho;
Ok(to_remove.len())
}
fn refresh_registered_atlas_transitions(term: &mut SaeManifoldTerm) -> Result<(), String> {
let registered: Vec<UnitSpeedChartTransition> = term
.chart_atlases()
.iter()
.flat_map(|atlas| atlas.transitions().iter().copied())
.filter(|transition| matches!(transition.seam_kind, AtlasSeamKind::Regular))
.collect();
for transition in registered {
let seam = fit_seam_transition(term, transition.to_chart, transition.from_chart)
.ok_or_else(|| {
format!(
"terminal atlas seam {}->{} is no longer identifiable",
transition.from_chart, transition.to_chart
)
})?;
if seam.sign as i8 != transition.sign {
return Err(format!(
"terminal atlas seam {}->{} changed orientation ({} -> {})",
transition.from_chart, transition.to_chart, transition.sign, seam.sign as i8
));
}
term.refresh_chart_transition(UnitSpeedChartTransition::new(
transition.from_chart,
transition.to_chart,
transition.sign,
seam.offset,
seam.period,
transition.seam_kind,
)?)?;
}
let registered_spheres: Vec<SphereChartTransition> = term
.chart_atlases()
.iter()
.flat_map(|atlas| atlas.sphere_transitions().iter().copied())
.collect();
for transition in registered_spheres {
let seam = fit_sphere_seam_transition(term, transition.to_chart, transition.from_chart)
.ok_or_else(|| {
format!(
"terminal sphere atlas seam {}->{} is no longer identifiable",
transition.from_chart, transition.to_chart
)
})?;
if seam.seam_kind != transition.seam_kind {
return Err(format!(
"terminal sphere atlas seam {}->{} changed kind ({:?} -> {:?})",
transition.from_chart, transition.to_chart, transition.seam_kind, seam.seam_kind
));
}
term.refresh_sphere_chart_transition(SphereChartTransition::new(
transition.from_chart,
transition.to_chart,
seam.rotation,
transition.seam_kind,
)?)?;
}
Ok(())
}
fn duplicate_atom(
term: &SaeManifoldTerm,
rho: &SaeManifoldRho,
parent: usize,
) -> Result<(SaeManifoldTerm, SaeManifoldRho), String> {
let k = term.k_atoms();
if parent >= k {
return Err(format!(
"duplicate_atom: parent {parent} out of range (K={k})"
));
}
let mut atoms = term.atoms.clone();
let mut child_atom = term.atoms[parent].clone();
{
let (m, p) = atoms[parent].decoder_coefficients.dim();
let s = |i: usize, j: usize| -> f64 {
let raw = ((i * 7 + j * 13) % 11) as f64 / 5.0 - 1.0;
if raw.abs() < 0.2 { 0.3 } else { raw }
};
for i in 0..m {
for j in 0..p {
let f = FISSION_SYMMETRY_BREAK_EPS * s(i, j);
atoms[parent].decoder_coefficients[[i, j]] *= 1.0 - f;
child_atom.decoder_coefficients[[i, j]] *= 1.0 + f;
}
}
atoms[parent].decoder_frame = None;
child_atom.decoder_frame = None;
}
atoms.push(child_atom);
let n = term.assignment.logits.nrows();
let mut logits = Array2::<f64>::zeros((n, k + 1));
let split = std::f64::consts::LN_2;
for row in 0..n {
for col in 0..k {
let mut v = term.assignment.logits[[row, col]];
if col == parent {
v -= split;
}
logits[[row, col]] = v;
}
logits[[row, k]] = term.assignment.logits[[row, parent]] - split;
}
let mut coords = term.assignment.coords.clone();
coords.push(term.assignment.coords[parent].clone());
let assignment =
crate::manifold::SaeAssignment::with_mode(logits, coords, term.assignment.mode)?;
let child = SaeManifoldTerm::new(atoms, assignment)?;
let mut child_rho = rho.clone();
if parent < child_rho.log_ard.len() {
let inherited = child_rho.log_ard[parent].clone();
child_rho.log_ard.push(inherited);
} else {
child_rho.log_ard.push(Array1::<f64>::zeros(0));
}
let inherited_smooth = child_rho
.log_lambda_smooth
.get(parent)
.or_else(|| child_rho.log_lambda_smooth.first())
.copied()
.unwrap_or(0.0);
child_rho.log_lambda_smooth.push(inherited_smooth);
Ok((child, child_rho))
}
#[derive(Clone)]
struct TopologyRaceFit {
evaluator: Arc<dyn SaeBasisSecondJet>,
basis_kind: SaeAtomBasisKind,
manifold: LatentManifold,
latent_dim: usize,
coords: Array2<f64>,
phi: Array2<f64>,
jet: ndarray::Array3<f64>,
decoder: Array2<f64>,
penalty: Array2<f64>,
}
struct TopologyCandidateSpec {
kind: AutoTopologyKind,
basis_kind: SaeAtomBasisKind,
manifold: LatentManifold,
latent_dim: usize,
evaluator: Arc<dyn SaeBasisSecondJet>,
coords: Array2<f64>,
}
fn topology_candidates_for_dim(
coords: ArrayView2<'_, f64>,
d_k: usize,
) -> Result<Vec<TopologyCandidateSpec>, String> {
let n = coords.nrows();
let d_seed = coords.ncols();
if d_k == 0 {
return Ok(Vec::new());
}
let coords_d = |d: usize| -> Array2<f64> {
let mut out = Array2::<f64>::zeros((n, d));
for row in 0..n {
for col in 0..d {
let src = col.min(d_seed.saturating_sub(1));
out[[row, col]] = coords[[row, src]];
}
}
out
};
let mut specs: Vec<TopologyCandidateSpec> = Vec::new();
match d_k {
1 => {
let n_harmonics = (2 * d_k + 1).max(3) | 1; specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Circle,
basis_kind: SaeAtomBasisKind::Periodic,
manifold: LatentManifold::Circle { period: 1.0 },
latent_dim: 1,
evaluator: Arc::new(PeriodicHarmonicEvaluator::new(n_harmonics)?),
coords: coords_d(1),
});
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Euclidean,
basis_kind: SaeAtomBasisKind::EuclideanPatch,
manifold: LatentManifold::Euclidean,
latent_dim: 1,
evaluator: Arc::new(EuclideanPatchEvaluator::new(1, 3)?),
coords: coords_d(1),
});
}
2 => {
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Torus,
basis_kind: SaeAtomBasisKind::Torus,
manifold: LatentManifold::Product(vec![
LatentManifold::Circle { period: 1.0 },
LatentManifold::Circle { period: 1.0 },
]),
latent_dim: 2,
evaluator: Arc::new(TorusHarmonicEvaluator::new(2, 2)?),
coords: coords_d(2),
});
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Sphere,
basis_kind: SaeAtomBasisKind::Sphere,
manifold: LatentManifold::Product(vec![
LatentManifold::Interval {
lo: -std::f64::consts::FRAC_PI_2,
hi: std::f64::consts::FRAC_PI_2,
},
LatentManifold::Circle {
period: std::f64::consts::TAU,
},
]),
latent_dim: 2,
evaluator: Arc::new(SphereChartEvaluator),
coords: coords_d(2),
});
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Euclidean,
basis_kind: SaeAtomBasisKind::EuclideanPatch,
manifold: LatentManifold::Euclidean,
latent_dim: 2,
evaluator: Arc::new(EuclideanPatchEvaluator::new(2, 2)?),
coords: coords_d(2),
});
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Cylinder,
basis_kind: SaeAtomBasisKind::Cylinder,
manifold: LatentManifold::Product(vec![
LatentManifold::Circle { period: 1.0 },
LatentManifold::Euclidean,
]),
latent_dim: 2,
evaluator: Arc::new(CylinderHarmonicEvaluator::new(2, 2)?),
coords: coords_d(2),
});
}
_ => {
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Euclidean,
basis_kind: SaeAtomBasisKind::EuclideanPatch,
manifold: LatentManifold::Euclidean,
latent_dim: d_k,
evaluator: Arc::new(EuclideanPatchEvaluator::new(d_k, 2)?),
coords: coords_d(d_k),
});
}
}
Ok(specs)
}
fn fit_topology_candidate(
spec: &TopologyCandidateSpec,
target: ArrayView2<'_, f64>,
weights: ArrayView1<'_, f64>,
) -> Result<TopologyAutoFitEvidence<TopologyRaceFit>, String> {
let n = target.nrows();
let (phi, jet) = spec.evaluator.evaluate(spec.coords.view())?;
let m = phi.ncols();
if phi.nrows() != n {
return Err(format!(
"fit_topology_candidate: basis rows {} != target rows {n}",
phi.nrows()
));
}
if weights.len() != n {
return Err(format!(
"fit_topology_candidate: weights length {} != target rows {n}",
weights.len()
));
}
let mut w_sum = 0.0_f64;
for row in 0..n {
let w = weights[row];
if !(w.is_finite() && w >= 0.0) {
return Err("fit_topology_candidate: weights must be finite and non-negative".into());
}
w_sum += w;
}
if !(w_sum > 0.0 && w_sum.is_finite()) {
return Err("fit_topology_candidate: degenerate (zero-mass) birth target".into());
}
let second_jet = spec.evaluator.second_jet(spec.coords.view())?; let d = spec.latent_dim;
let mut s_raw = Array2::<f64>::zeros((m, m));
for row in 0..n {
for a in 0..d {
for c in 0..d {
for mu in 0..m {
let hmu = second_jet[[row, mu, a, c]];
if hmu == 0.0 {
continue;
}
for nu in mu..m {
s_raw[[mu, nu]] += hmu * second_jet[[row, nu, a, c]];
}
}
}
}
}
for mu in 0..m {
for nu in (mu + 1)..m {
s_raw[[nu, mu]] = s_raw[[mu, nu]];
}
}
let reml_fit =
gaussian_reml_multi_closed_form(phi.view(), target, s_raw.view(), Some(weights), None)
.map_err(|e| format!("fit_topology_candidate: REML evidence: {e:?}"))?;
let lambda = reml_fit.lambda;
if !(lambda.is_finite() && lambda >= 0.0) {
return Err(format!(
"fit_topology_candidate: REML returned a non-finite/negative λ ({lambda})"
));
}
let raw_reml = reml_fit.reml_score;
if !raw_reml.is_finite() {
return Err("fit_topology_candidate: non-finite REML score".into());
}
let decoder = reml_fit.coefficients.clone(); let mut effective_dim = reml_fit.edf;
if !(effective_dim.is_finite() && effective_dim > 0.0) {
effective_dim = 1.0;
}
let penalty = s_raw.clone();
Ok(TopologyAutoFitEvidence {
topology_name: spec.kind.as_str().to_string(),
raw_reml,
null_dim: 0.0,
null_space_logdet: None,
effective_dim,
n_obs: n,
fit_handle: TopologyRaceFit {
evaluator: spec.evaluator.clone(),
basis_kind: spec.basis_kind.clone(),
manifold: spec.manifold.clone(),
latent_dim: spec.latent_dim,
coords: spec.coords.clone(),
phi,
jet,
decoder,
penalty,
},
})
}
fn birth_row_amplitudes(target: ArrayView2<'_, f64>) -> Array1<f64> {
let n = target.nrows();
let mut amps = Array1::<f64>::zeros(n);
for i in 0..n {
let mut ss = 0.0_f64;
for &v in target.row(i).iter() {
ss += v * v;
}
amps[i] = ss.sqrt();
}
amps
}
fn standardized_log_birth_amplitudes(amps: ArrayView1<'_, f64>) -> Option<Array1<f64>> {
let n = amps.len();
if n == 0 {
return None;
}
let mut logs = Array1::<f64>::zeros(n);
for (i, &) in amps.iter().enumerate() {
if !amp.is_finite() || amp < 0.0 {
return None;
}
logs[i] = amp.max(f64::MIN_POSITIVE).ln();
}
let mean = logs.sum() / n as f64;
let mut var = 0.0_f64;
for &value in logs.iter() {
let centered = value - mean;
var += centered * centered;
}
let std = (var / n as f64).sqrt();
if !std.is_finite() || std <= 0.0 {
return None;
}
for value in logs.iter_mut() {
*value = (*value - mean) / std;
}
Some(logs)
}
fn radial_promoted_specs(
coords: ArrayView2<'_, f64>,
target: ArrayView2<'_, f64>,
d_k: usize,
) -> Result<Option<Vec<TopologyCandidateSpec>>, String> {
if d_k != 1 {
return Ok(None);
}
let amps = birth_row_amplitudes(target);
let cert = amplitude_concentration_certificate(amps.view());
if !cert.recommends_radial_axis() {
return Ok(None);
}
let log_amp_coord = standardized_log_birth_amplitudes(amps.view())
.ok_or_else(|| "radial_promoted_specs: degenerate log-amplitude spread".to_string())?;
let mut promoted: Vec<TopologyCandidateSpec> = Vec::with_capacity(3);
for spec in topology_candidates_for_dim(coords, 1)? {
if spec.kind == AutoTopologyKind::Circle {
promoted.push(spec);
}
}
for mut spec in topology_candidates_for_dim(coords, 2)? {
if matches!(
spec.kind,
AutoTopologyKind::Cylinder | AutoTopologyKind::Euclidean
) {
for row in 0..spec.coords.nrows() {
spec.coords[[row, 1]] = log_amp_coord[row];
}
promoted.push(spec);
}
}
if promoted.len() < 2 {
return Ok(None);
}
Ok(Some(promoted))
}
static FINITE_SET_RACE_ENROLLED: AtomicBool = AtomicBool::new(false);
pub fn finite_set_race_enrolled() -> bool {
FINITE_SET_RACE_ENROLLED.load(Ordering::Relaxed)
}
pub fn set_finite_set_race_enrolled(enrolled: bool) {
FINITE_SET_RACE_ENROLLED.store(enrolled, Ordering::Relaxed);
}
pub fn finite_set_candidate_for_birth(coords: ArrayView2<'_, f64>) -> Option<(usize, Array2<f64>)> {
if coords.ncols() != 1 {
return None;
}
let n = coords.nrows();
if n < 4 {
return None;
}
let col = coords.column(0);
let (mut lo, mut hi) = (f64::INFINITY, f64::NEG_INFINITY);
for &t in col.iter() {
if !t.is_finite() {
return None;
}
lo = lo.min(t);
hi = hi.max(t);
}
let span = hi - lo;
if !(span > 0.0) {
return None;
}
let r: Vec<f64> = col
.iter()
.map(|&t| ((t - lo) / span).clamp(0.0, 1.0))
.collect();
match classify_occupancy_interval(&r) {
OccupancyLaw::Discrete { anchors } if anchors >= 2 => {
let mut idx = Array2::<f64>::zeros((n, 1));
for i in 0..n {
let bin = (r[i] * anchors as f64).floor();
idx[[i, 0]] = bin.clamp(0.0, (anchors - 1) as f64);
}
Some((anchors, idx))
}
_ => None,
}
}
#[derive(Clone, Debug)]
pub struct GraphBirthCandidate {
pub atom: LearnedGraphAtom,
pub selection: GraphStructureSelection,
}
pub fn graph_birth_candidate_for_structure_search(
anchor_embeddings: ArrayView2<'_, f64>,
row_coordinates: &[f64],
n_eff: f64,
edge_precisions: &[f64],
edge_delta_loss: &[f64],
) -> Result<GraphBirthCandidate, String> {
let atom = LearnedGraphAtom::from_reml_knn_edges(
anchor_embeddings,
row_coordinates,
n_eff,
edge_precisions,
edge_delta_loss,
)?;
let selection = atom.structure_selection();
Ok(GraphBirthCandidate { atom, selection })
}
fn race_birth_topology(
coords: ArrayView2<'_, f64>,
target: ArrayView2<'_, f64>,
weights: ArrayView1<'_, f64>,
d_k: usize,
) -> Result<Option<TopologyRaceFit>, String> {
let base_specs = topology_candidates_for_dim(coords, d_k)?;
if base_specs.is_empty() {
return Ok(None);
}
if let Ok(Some(promoted)) = radial_promoted_specs(coords, target, d_k) {
if !promoted.is_empty() {
if let Ok(Some(fit)) = race_spec_set(promoted, target, weights) {
return Ok(Some(fit));
}
}
}
race_spec_set(base_specs, target, weights)
}
fn race_spec_set(
specs: Vec<TopologyCandidateSpec>,
target: ArrayView2<'_, f64>,
weights: ArrayView1<'_, f64>,
) -> Result<Option<TopologyRaceFit>, String> {
if specs.is_empty() {
return Ok(None);
}
let selector = TopologyAutoSelector {
candidates: specs.iter().map(|s| s.kind).collect(),
score_scale: TopologyScoreScale::PerObservation,
latent: None,
};
let mut by_kind: std::collections::HashMap<AutoTopologyKind, &TopologyCandidateSpec> =
std::collections::HashMap::with_capacity(specs.len() + 1);
for spec in &specs {
by_kind.insert(spec.kind, spec);
}
if !by_kind.contains_key(&AutoTopologyKind::ConstantCurvature) {
if let Some(sphere) = specs.iter().find(|s| s.kind == AutoTopologyKind::Sphere) {
by_kind.insert(AutoTopologyKind::ConstantCurvature, sphere);
} else if let Some(euclid) = specs.iter().find(|s| s.kind == AutoTopologyKind::Euclidean) {
by_kind.insert(AutoTopologyKind::ConstantCurvature, euclid);
}
}
let ranked = select_topology_with_fit(&selector, |kind| {
let spec = by_kind.get(&kind).ok_or_else(|| {
format!(
"race_birth_topology: no realized candidate for fused topology {:?}",
kind.as_str()
)
})?;
fit_topology_candidate(spec, target, weights)
})?;
let winner = ranked
.winner()
.ok_or_else(|| "race_birth_topology: empty ranking".to_string())?;
Ok(Some(winner.fit_handle.clone()))
}
pub struct PrimaryTopologyChoice {
pub basis_kind: SaeAtomBasisKind,
pub latent_dim: usize,
pub n_harmonics: Option<usize>,
pub n_duchon_centers: Option<usize>,
}
pub fn discover_primary_atom_topologies(
target: ArrayView2<'_, f64>,
labels: &[usize],
k_atoms: usize,
max_dims: &[usize],
) -> Result<Vec<PrimaryTopologyChoice>, String> {
let n_obs = target.nrows();
let p_out = target.ncols();
if labels.len() != n_obs {
return Err(format!(
"discover_primary_atom_topologies: labels must have N={n_obs} entries; got {}",
labels.len()
));
}
if max_dims.len() != k_atoms {
return Err(format!(
"discover_primary_atom_topologies: max_dims must have K={k_atoms} entries; got {}",
max_dims.len()
));
}
if p_out < 2 {
return Err(format!(
"discover_primary_atom_topologies: evidence racing needs at least two output dimensions; got P={p_out}"
));
}
(0..k_atoms)
.map(|atom_idx| -> Result<PrimaryTopologyChoice, String> {
let rows: Vec<usize> =
(0..n_obs).filter(|&row| labels[row] == atom_idx).collect();
if rows.len() < 16 {
return Err(format!(
"discover_primary_atom_topologies: auto atom {atom_idx} has only {} seed-cluster rows; at least 16 are required for an evidence race (name an explicit topology when discovery is not identifiable)",
rows.len()
));
}
let mut mean = vec![0.0_f64; p_out];
for &row in &rows {
for col in 0..p_out {
mean[col] += target[[row, col]];
}
}
let inv_count = 1.0 / rows.len() as f64;
for value in &mut mean {
*value *= inv_count;
}
let mut local = Array2::<f64>::zeros((rows.len(), p_out));
for (out_row, &src_row) in rows.iter().enumerate() {
for col in 0..p_out {
local[[out_row, col]] = target[[src_row, col]] - mean[col];
}
}
let (_u, _s, vt_opt) = local.svd(false, true).map_err(|error| {
format!(
"discover_primary_atom_topologies: SVD failed for auto atom {atom_idx}: {error}"
)
})?;
let vt = vt_opt.ok_or_else(|| {
format!(
"discover_primary_atom_topologies: SVD returned no right-singular frame for auto atom {atom_idx}"
)
})?;
let n_pcs = vt.nrows().min(4);
if n_pcs < 2 {
return Err(format!(
"discover_primary_atom_topologies: auto atom {atom_idx} has principal rank {n_pcs}; at least two directions are required"
));
}
let mut proj = Array2::<f64>::zeros((n_obs, n_pcs));
for row in 0..n_obs {
for pc in 0..n_pcs {
let mut acc = 0.0_f64;
for col in 0..p_out {
acc += (target[[row, col]] - mean[col]) * vt[[pc, col]];
}
proj[[row, pc]] = acc;
}
}
let cluster_sd = |pc: usize| -> f64 {
let mut acc = 0.0_f64;
for &row in &rows {
acc += proj[[row, pc]] * proj[[row, pc]];
}
(acc * inv_count).sqrt().max(1e-12)
};
let phase = |a: f64, b: f64| -> f64 {
let frac = b.atan2(a) / std::f64::consts::TAU;
frac - frac.floor()
};
let mut specs: Vec<TopologyCandidateSpec> = Vec::with_capacity(4);
let circle_coords = {
let mut coords = Array2::<f64>::zeros((n_obs, 1));
for row in 0..n_obs {
coords[[row, 0]] = phase(proj[[row, 0]], proj[[row, 1]]);
}
let n_harmonics = 3;
let evaluator = PeriodicHarmonicEvaluator::new(n_harmonics).map_err(|error| {
format!(
"discover_primary_atom_topologies: circle evaluator failed for auto atom {atom_idx}: {error}"
)
})?;
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Circle,
basis_kind: SaeAtomBasisKind::Periodic,
manifold: LatentManifold::Circle { period: 1.0 },
latent_dim: 1,
evaluator: Arc::new(evaluator),
coords: coords.clone(),
});
coords
};
let mut sheet_coords: Option<Array2<f64>> = None;
if max_dims[atom_idx] >= 2 {
let (sd0, sd1) = (cluster_sd(0), cluster_sd(1));
let mut coords = Array2::<f64>::zeros((n_obs, 2));
for row in 0..n_obs {
coords[[row, 0]] = proj[[row, 0]] / sd0;
coords[[row, 1]] = proj[[row, 1]] / sd1;
}
let evaluator = EuclideanPatchEvaluator::new(2, 2).map_err(|error| {
format!(
"discover_primary_atom_topologies: flat evaluator failed for auto atom {atom_idx}: {error}"
)
})?;
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Euclidean,
basis_kind: SaeAtomBasisKind::EuclideanPatch,
manifold: LatentManifold::Euclidean,
latent_dim: 2,
evaluator: Arc::new(evaluator),
coords: coords.clone(),
});
if let Some(centers) =
duchon_sheet_centers(&coords, &rows, duchon_sheet_race_center_budget(rows.len()))
{
let evaluator = DuchonCoordinateEvaluator::new(centers, DUCHON_SHEET_M)
.map_err(|error| {
format!(
"discover_primary_atom_topologies: duchon-sheet evaluator failed for auto atom {atom_idx}: {error}"
)
})?;
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::DuchonSheet,
basis_kind: SaeAtomBasisKind::Duchon,
manifold: LatentManifold::Euclidean,
latent_dim: 2,
evaluator: Arc::new(evaluator),
coords: coords.clone(),
});
}
sheet_coords = Some(coords);
if n_pcs >= 3 {
let mut coords = Array2::<f64>::zeros((n_obs, 2));
for row in 0..n_obs {
let (x, y, z) = (proj[[row, 0]], proj[[row, 1]], proj[[row, 2]]);
let norm = (x * x + y * y + z * z).sqrt().max(1e-12);
coords[[row, 0]] = (z / norm).clamp(-1.0, 1.0).asin();
coords[[row, 1]] = y.atan2(x);
}
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Sphere,
basis_kind: SaeAtomBasisKind::Sphere,
manifold: LatentManifold::Product(vec![
LatentManifold::Interval {
lo: -std::f64::consts::FRAC_PI_2,
hi: std::f64::consts::FRAC_PI_2,
},
LatentManifold::Circle {
period: std::f64::consts::TAU,
},
]),
latent_dim: 2,
evaluator: Arc::new(SphereChartEvaluator),
coords,
});
}
if n_pcs >= 3 {
if let (Ok(coords), Ok(evaluator)) = (
crate::manifold::mobius_double_cover_coords_from_projection(
proj.view(),
&rows,
),
MobiusHarmonicEvaluator::new(3, 2),
) {
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Mobius,
basis_kind: SaeAtomBasisKind::Mobius,
manifold: LatentManifold::Product(vec![
LatentManifold::Circle { period: 2.0 },
LatentManifold::Interval { lo: -1.0, hi: 1.0 },
]),
latent_dim: 2,
evaluator: Arc::new(evaluator),
coords,
});
}
}
if n_pcs >= 4 {
let mut coords = Array2::<f64>::zeros((n_obs, 2));
for row in 0..n_obs {
coords[[row, 0]] = phase(proj[[row, 0]], proj[[row, 1]]);
coords[[row, 1]] = phase(proj[[row, 2]], proj[[row, 3]]);
}
let evaluator = TorusHarmonicEvaluator::new(2, 2).map_err(|error| {
format!(
"discover_primary_atom_topologies: torus evaluator failed for auto atom {atom_idx}: {error}"
)
})?;
specs.push(TopologyCandidateSpec {
kind: AutoTopologyKind::Torus,
basis_kind: SaeAtomBasisKind::Torus,
manifold: LatentManifold::Product(vec![
LatentManifold::Circle { period: 1.0 },
LatentManifold::Circle { period: 1.0 },
]),
latent_dim: 2,
evaluator: Arc::new(evaluator),
coords,
});
}
}
if specs.is_empty() {
return Err(format!(
"discover_primary_atom_topologies: auto atom {atom_idx} produced no realizable candidates"
));
}
let mut weights = Array1::<f64>::zeros(n_obs);
for &row in &rows {
weights[row] = 1.0;
}
let fit = race_spec_set(specs, target, weights.view()).map_err(|error| {
format!(
"discover_primary_atom_topologies: evidence race failed for auto atom {atom_idx}: {error}"
)
})?;
let fit = fit.ok_or_else(|| {
format!(
"discover_primary_atom_topologies: evidence race returned no winner for auto atom {atom_idx}"
)
})?;
let n_harmonics = if fit.basis_kind == SaeAtomBasisKind::Periodic {
Some(select_periodic_resolution(
circle_coords.view(),
target,
weights.view(),
rows.len(),
)?)
} else {
None
};
let n_duchon_centers = if fit.basis_kind == SaeAtomBasisKind::Duchon {
let coords = sheet_coords.as_ref().ok_or_else(|| {
format!(
"discover_primary_atom_topologies: duchon-sheet winner without a 2-D chart for auto atom {atom_idx}"
)
})?;
Some(select_duchon_sheet_resolution(
coords,
target,
weights.view(),
&rows,
)?)
} else {
None
};
Ok(PrimaryTopologyChoice {
basis_kind: fit.basis_kind,
latent_dim: fit.latent_dim,
n_harmonics,
n_duchon_centers,
})
})
.collect()
}
const DUCHON_SHEET_M: usize = 3;
const DUCHON_SHEET_NULLSPACE_DIM: usize = 6;
fn duchon_sheet_race_center_budget(n_cluster: usize) -> usize {
let floor = DUCHON_SHEET_NULLSPACE_DIM + 2 + 1;
if n_cluster <= floor {
return 0;
}
n_cluster.min(32).max(floor)
}
fn duchon_sheet_centers(
coords: &Array2<f64>,
rows: &[usize],
n_centers: usize,
) -> Option<Array2<f64>> {
if n_centers == 0 || rows.len() < n_centers {
return None;
}
let mut centers = Array2::<f64>::zeros((n_centers, 2));
for i in 0..n_centers {
let row = rows[i * rows.len() / n_centers];
centers[[i, 0]] = coords[[row, 0]];
centers[[i, 1]] = coords[[row, 1]];
}
Some(centers)
}
fn select_duchon_sheet_resolution(
sheet_coords: &Array2<f64>,
target: ArrayView2<'_, f64>,
weights: ArrayView1<'_, f64>,
rows: &[usize],
) -> Result<usize, String> {
let floor = duchon_sheet_race_center_budget(rows.len());
if floor == 0 {
return Err(
"select_duchon_sheet_resolution: cluster too small to identify the thin-plate nullspace"
.to_string(),
);
}
let ceiling = rows.len().saturating_sub(1).max(floor);
let mut ladder: Vec<usize> = Vec::new();
let mut c = floor;
while c < ceiling {
ladder.push(c);
c = c.saturating_mul(2);
}
ladder.push(ceiling);
let mut best_c = 0usize;
let mut best_score = f64::INFINITY;
for &n_centers in &ladder {
let Some(centers) = duchon_sheet_centers(sheet_coords, rows, n_centers) else {
continue;
};
let evaluator = match DuchonCoordinateEvaluator::new(centers, DUCHON_SHEET_M) {
Ok(evaluator) => evaluator,
Err(_) => continue,
};
let spec = TopologyCandidateSpec {
kind: AutoTopologyKind::DuchonSheet,
basis_kind: SaeAtomBasisKind::Duchon,
manifold: LatentManifold::Euclidean,
latent_dim: 2,
evaluator: Arc::new(evaluator),
coords: sheet_coords.clone(),
};
let score = match fit_topology_candidate(&spec, target, weights) {
Ok(evidence) => evidence.raw_reml,
Err(_) => continue,
};
if score.is_finite() && score < best_score {
best_score = score;
best_c = n_centers;
}
}
if best_c == 0 {
return Err(
"select_duchon_sheet_resolution: no fittable center count for the duchon-sheet winner"
.to_string(),
);
}
Ok(best_c)
}
fn select_periodic_resolution(
circle_coords: ArrayView2<'_, f64>,
target: ArrayView2<'_, f64>,
weights: ArrayView1<'_, f64>,
n_cluster: usize,
) -> Result<usize, String> {
let n_obs = target.nrows();
let p_out = target.ncols();
let ident_ceiling = (n_cluster.saturating_sub(2) / 2).max(1);
let mut peak_energy = 0.0_f64;
let mut energies = Vec::with_capacity(ident_ceiling);
for h in 1..=ident_ceiling {
let mut energy = 0.0_f64;
for col in 0..p_out {
let (mut re, mut im) = (0.0_f64, 0.0_f64);
for row in 0..n_obs {
let w = weights[row];
if w == 0.0 {
continue;
}
let angle = std::f64::consts::TAU * h as f64 * circle_coords[[row, 0]];
re += w * target[[row, col]] * angle.cos();
im += w * target[[row, col]] * angle.sin();
}
energy += re * re + im * im;
}
peak_energy = peak_energy.max(energy);
energies.push(energy);
}
if !(peak_energy > 0.0) {
return Err(
"select_periodic_resolution: the circle winner carries no angular energy".to_string(),
);
}
let energy_floor = peak_energy * 1e-12;
let bandwidth = energies
.iter()
.rposition(|&energy| energy > energy_floor)
.map(|idx| idx + 1)
.unwrap_or(1);
let ceiling = bandwidth.min(ident_ceiling).max(1);
let mut best_h = 0usize;
let mut best_score = f64::INFINITY;
for h in 1..=ceiling {
let evaluator = match PeriodicHarmonicEvaluator::new(h) {
Ok(evaluator) => evaluator,
Err(_) => continue,
};
let spec = TopologyCandidateSpec {
kind: AutoTopologyKind::Circle,
basis_kind: SaeAtomBasisKind::Periodic,
manifold: LatentManifold::Circle { period: 1.0 },
latent_dim: 1,
evaluator: Arc::new(evaluator),
coords: circle_coords.to_owned(),
};
let score = match fit_topology_candidate(&spec, target, weights) {
Ok(evidence) => evidence.raw_reml,
Err(_) => continue,
};
if score.is_finite() && score < best_score {
best_score = score;
best_h = h;
}
}
if best_h == 0 {
return Err(
"select_periodic_resolution: no fittable harmonic resolution for the circle winner"
.to_string(),
);
}
Ok(best_h)
}
pub fn resolve_auto_primary_atoms(
target: ArrayView2<'_, f64>,
labels: &[usize],
atom_basis: &mut [String],
atom_dim: &mut [usize],
) -> Result<Vec<Option<usize>>, String> {
let k_atoms = atom_basis.len();
if atom_dim.len() != k_atoms {
return Err(format!(
"resolve_auto_primary_atoms: atom_basis and atom_dim must both have K={k_atoms} entries; atom_dim has {}",
atom_dim.len()
));
}
let mut duchon_center_overrides: Vec<Option<usize>> = vec![None; k_atoms];
if !atom_basis.iter().any(|basis| basis == "auto") {
return Ok(duchon_center_overrides);
}
let choices = discover_primary_atom_topologies(target, labels, k_atoms, atom_dim)?;
for atom_idx in 0..k_atoms {
if atom_basis[atom_idx] != "auto" {
continue;
}
let choice = &choices[atom_idx];
match choice.basis_kind {
SaeAtomBasisKind::Torus => {
atom_basis[atom_idx] = "torus".to_string();
atom_dim[atom_idx] = choice.latent_dim;
}
SaeAtomBasisKind::Sphere => {
atom_basis[atom_idx] = "sphere".to_string();
atom_dim[atom_idx] = choice.latent_dim;
}
SaeAtomBasisKind::Mobius => {
atom_basis[atom_idx] = "mobius".to_string();
atom_dim[atom_idx] = choice.latent_dim;
}
SaeAtomBasisKind::EuclideanPatch => {
atom_basis[atom_idx] = "duchon".to_string();
atom_dim[atom_idx] = choice.latent_dim;
}
SaeAtomBasisKind::Duchon => {
atom_basis[atom_idx] = "duchon".to_string();
atom_dim[atom_idx] = choice.latent_dim;
duchon_center_overrides[atom_idx] = choice.n_duchon_centers;
}
SaeAtomBasisKind::Periodic => {
atom_basis[atom_idx] = "periodic".to_string();
if let Some(n_harmonics) = choice.n_harmonics {
atom_dim[atom_idx] = n_harmonics;
}
}
ref unexpected => {
return Err(format!(
"resolve_auto_primary_atoms: evidence race selected unsupported primary basis {unexpected:?} for auto atom {atom_idx}"
));
}
}
}
Ok(duchon_center_overrides)
}
const BIRTH_SEED_LOGIT: f64 = -4.0;
fn born_atom(
term: &SaeManifoldTerm,
rho: &SaeManifoldRho,
factor_dir: ArrayView2<'_, f64>,
) -> Result<(SaeManifoldTerm, SaeManifoldRho), String> {
let k = term.k_atoms();
if term.atoms.is_empty() {
return Err(
"born_atom: cannot birth from an empty dictionary (no template atom to seed the \
coordinate block / basis from)"
.to_string(),
);
}
let template = &term.atoms[0];
let m = template.basis_size();
let p = term.output_dim();
if factor_dir.dim() != (m, p) {
return Err(format!(
"born_atom: residual-factor decoder must be ({m}, {p}); got {:?}",
factor_dir.dim()
));
}
let mut atoms = term.atoms.clone();
let template_coords = term.assignment.coords[0].as_matrix();
let birth_target = template.basis_values.dot(&factor_dir); let weights = Array1::<f64>::ones(birth_target.nrows());
let raced = race_birth_topology(
template_coords.view(),
birth_target.view(),
weights.view(),
template.latent_dim,
)?;
let (born, born_coord_block) = match raced {
Some(fit) => {
let mut atom = SaeManifoldAtom::new(
format!("atom_born_{k}"),
fit.basis_kind.clone(),
fit.latent_dim,
fit.phi.clone(),
fit.jet.clone(),
fit.decoder.clone(),
fit.penalty.clone(),
)?
.with_basis_second_jet(fit.evaluator.clone());
atom.refresh_intrinsic_smooth_penalty();
let coord_block = gam_terms::latent::LatentCoordValues::from_matrix_with_manifold(
fit.coords.view(),
LatentIdMode::None,
fit.manifold.clone(),
);
(atom, coord_block)
}
None => {
let mut atom = template.clone();
atom.decoder_coefficients = factor_dir.to_owned();
atom.refresh_intrinsic_smooth_penalty();
(atom, term.assignment.coords[0].clone())
}
};
atoms.push(born);
let n = term.assignment.logits.nrows();
let mut logits = Array2::<f64>::zeros((n, k + 1));
for row in 0..n {
for col in 0..k {
logits[[row, col]] = term.assignment.logits[[row, col]];
}
logits[[row, k]] = BIRTH_SEED_LOGIT;
}
let mut coords = term.assignment.coords.clone();
coords.push(born_coord_block);
let assignment =
crate::manifold::SaeAssignment::with_mode(logits, coords, term.assignment.mode)?;
let child = SaeManifoldTerm::new(atoms, assignment)?;
let mut child_rho = rho.clone();
let inherited = child_rho
.log_ard
.first()
.cloned()
.unwrap_or_else(|| Array1::<f64>::zeros(0));
child_rho.log_ard.push(inherited);
let inherited_smooth = child_rho.log_lambda_smooth.first().copied().unwrap_or(0.0);
child_rho.log_lambda_smooth.push(inherited_smooth);
Ok((child, child_rho))
}
pub(crate) fn born_circle_atom(
term: &SaeManifoldTerm,
rho: &SaeManifoldRho,
harmonic_decoder: Array2<f64>,
phase_coords: Array2<f64>,
circle_gate: Vec<f64>,
) -> Result<(SaeManifoldTerm, SaeManifoldRho), String> {
let k = term.k_atoms();
if term.atoms.is_empty() {
return Err("born_circle_atom: cannot birth from an empty dictionary".to_string());
}
let m = harmonic_decoder.nrows();
let p = term.output_dim();
if m % 2 != 1 || m < 3 {
return Err(format!(
"born_circle_atom: harmonic decoder must have odd height >= 3 (constant + \
>= 1 sin/cos harmonic pair); got height {m}"
));
}
if harmonic_decoder.ncols() != p {
return Err(format!(
"born_circle_atom: harmonic decoder must have {p} columns (output dim); got {}",
harmonic_decoder.ncols()
));
}
let n = term.assignment.logits.nrows();
if phase_coords.dim() != (n, 1) {
return Err(format!(
"born_circle_atom: phase coords must be ({n}, 1); got {:?}",
phase_coords.dim()
));
}
let evaluator = std::sync::Arc::new(crate::manifold::PeriodicHarmonicEvaluator::new(m)?);
let (phi, jet) = {
use crate::manifold::SaeBasisEvaluator;
evaluator.evaluate(phase_coords.view())?
};
let mut born = SaeManifoldAtom::new(
format!("atom_born_{k}"),
SaeAtomBasisKind::Periodic,
1,
phi,
jet,
harmonic_decoder,
Array2::<f64>::eye(m),
)?
.with_basis_second_jet(evaluator.clone());
born.refresh_intrinsic_smooth_penalty();
let born_coord_block = gam_terms::latent::LatentCoordValues::from_matrix_with_manifold(
phase_coords.view(),
LatentIdMode::None,
LatentManifold::Circle { period: 1.0 },
);
let mut atoms = term.atoms.clone();
atoms.push(born);
let mut logits = Array2::<f64>::zeros((n, k + 1));
for row in 0..n {
for col in 0..k {
logits[[row, col]] = term.assignment.logits[[row, col]];
}
let own_gate = circle_gate.get(row).copied().unwrap_or(f64::NEG_INFINITY);
let inc_max = (0..k)
.map(|c| term.assignment.logits[[row, c]])
.fold(f64::NEG_INFINITY, f64::max);
logits[[row, k]] = if own_gate.is_finite() {
if inc_max.is_finite() {
inc_max.max(own_gate)
} else {
own_gate
}
} else {
BIRTH_SEED_LOGIT
};
}
let mut coords = term.assignment.coords.clone();
coords.push(born_coord_block);
let assignment =
crate::manifold::SaeAssignment::with_mode(logits, coords, term.assignment.mode)?;
let child = SaeManifoldTerm::new(atoms, assignment)?;
let mut child_rho = rho.clone();
let inherited = child_rho
.log_ard
.first()
.cloned()
.unwrap_or_else(|| Array1::<f64>::zeros(0));
child_rho.log_ard.push(inherited);
let inherited_smooth = child_rho.log_lambda_smooth.first().copied().unwrap_or(0.0);
child_rho.log_lambda_smooth.push(inherited_smooth);
Ok((child, child_rho))
}
#[derive(Clone, Debug)]
pub struct RowBlockShard {
pub target: std::sync::Arc<Array2<f64>>,
pub rows: Vec<usize>,
}
#[derive(Clone, Debug)]
pub struct EstimationEvalSplit {
pub estimation_rows: Vec<usize>,
pub shards: Vec<RowBlockShard>,
}
const ESTIMATION_FRACTION: f64 = 0.6;
pub fn estimation_eval_split(target: ArrayView2<'_, f64>, n_shards: usize) -> EstimationEvalSplit {
let n = target.nrows();
if n == 0 {
return EstimationEvalSplit {
estimation_rows: Vec::new(),
shards: Vec::new(),
};
}
let shared = std::sync::Arc::new(target.to_owned());
let n_est =
((n as f64 * ESTIMATION_FRACTION).round() as usize).clamp(1, n.saturating_sub(1).max(1));
let estimation_rows: Vec<usize> = (0..n_est).collect();
let eval_rows: Vec<usize> = (n_est..n).collect();
let n_eval = eval_rows.len();
let n_shards = n_shards.min(n_eval).max(usize::from(n_eval > 0));
let mut shards = Vec::new();
if n_eval > 0 && n_shards > 0 {
let base = n_eval / n_shards;
let rem = n_eval % n_shards;
let mut cursor = 0usize;
for s in 0..n_shards {
let len = base + usize::from(s < rem);
let rows: Vec<usize> = eval_rows[cursor..cursor + len].to_vec();
shards.push(RowBlockShard {
target: shared.clone(),
rows,
});
cursor += len;
}
}
EstimationEvalSplit {
estimation_rows,
shards,
}
}
pub struct StructureSearchResult {
pub term: SaeManifoldTerm,
pub rho: SaeManifoldRho,
pub rounds: Vec<SearchLedger>,
pub migration: SaeMigrationLedger,
}
impl StructureSearchResult {
#[must_use]
pub fn from_rounds(
term: SaeManifoldTerm,
rho: SaeManifoldRho,
rounds: Vec<SearchLedger>,
) -> Self {
Self::from_rounds_with_predictions(term, rho, rounds, &[])
}
#[must_use]
pub fn from_rounds_with_predictions(
term: SaeManifoldTerm,
rho: SaeManifoldRho,
rounds: Vec<SearchLedger>,
birth_predictions: &[std::collections::HashMap<usize, f64>],
) -> Self {
let mut migration = SaeMigrationLedger::new();
let empty = std::collections::HashMap::new();
for (round_idx, round_ledger) in rounds.iter().enumerate() {
let preds = birth_predictions.get(round_idx).unwrap_or(&empty);
migration.record_search_round(round_idx, round_ledger, preds);
}
Self {
term,
rho,
rounds,
migration,
}
}
#[must_use]
pub fn structure_changed(&self) -> bool {
use gam_solve::structure_search::MoveVerdict;
self.rounds.iter().any(|round| {
round.moves.iter().any(|record| {
matches!(
record.verdict,
MoveVerdict::Accepted { .. } | MoveVerdict::Demoted { .. }
)
})
})
}
}
#[derive(Clone, Copy, Debug)]
pub struct RoundDriverConfig {
pub n_shards: usize,
pub budget: MoveBudget,
pub max_rounds: usize,
pub harvest_params: HarvestParams,
pub curl: Option<CurlConfig>,
}
#[derive(Clone, Copy, Debug)]
pub struct CurlConfig {
pub coalesce_cos_threshold: f64,
pub coalesce_max_overlap: f64,
pub min_cooccurrence: usize,
pub subsample_rows: usize,
pub harmonics: usize,
pub max_curls: usize,
pub flatten: bool,
pub cooldown_rounds: usize,
}
impl Default for CurlConfig {
fn default() -> Self {
Self {
coalesce_cos_threshold: -0.85,
coalesce_max_overlap: 0.15,
min_cooccurrence: 8,
subsample_rows: 4096,
harmonics: 1,
max_curls: 4,
flatten: true,
cooldown_rounds: 2,
}
}
}
pub fn run_structure_search_rounds(
mut term: SaeManifoldTerm,
mut rho: SaeManifoldRho,
target: ArrayView2<'_, f64>,
config: RoundDriverConfig,
ledger: &mut StructureLedger,
mut candidate_fit: impl FnMut(
SaeManifoldTerm,
SaeManifoldRho,
&[usize],
) -> (SaeManifoldTerm, SaeManifoldRho),
mut finalize_round: impl FnMut(
SaeManifoldTerm,
SaeManifoldRho,
&[usize],
) -> (SaeManifoldTerm, SaeManifoldRho),
) -> Result<StructureSearchResult, String> {
let RoundDriverConfig {
n_shards,
budget,
max_rounds,
harvest_params,
curl,
} = config;
let split = estimation_eval_split(target, n_shards);
let mut rounds: Vec<SearchLedger> = Vec::new();
let mut round_predictions: Vec<std::collections::HashMap<usize, f64>> = Vec::new();
let mut cooldown = crate::manifold::CurlCooldownLedger::new();
for _ in 0..max_rounds {
let fitted = term.try_fitted_target_aware(target, None)?;
let residuals = &target.to_owned() - &fitted;
let mut report = harvest_move_proposals(&term, &rho, residuals.view(), &harvest_params)?;
let birth_predictions: std::collections::HashMap<usize, f64> =
report.birth_predictions.iter().copied().collect();
let residual_decoders = build_birth_decoders(&term, residuals.view(), &harvest_params)?;
let mut birth_seeds: Vec<BirthSeed> = residual_decoders
.into_iter()
.map(BirthSeed::ResidualFactor)
.collect();
let mut curl_atoms: std::collections::HashMap<usize, Vec<usize>> =
std::collections::HashMap::new();
let mut flatten_atoms: std::collections::HashSet<usize> = std::collections::HashSet::new();
if let Some(cfg) = curl {
for cand in curl_candidates(&term, residuals.view(), &cfg) {
if cooldown.blocked(&cand.members) {
continue;
}
let candidate = birth_seeds.len();
birth_seeds.push(cand.seed);
curl_atoms.insert(candidate, cand.members.clone());
report.proposals.push(proposal(
&term,
StructureMove::Birth { candidate },
cand.net_evidence,
));
}
if cfg.flatten {
for atom in flatten_candidates(&term) {
if cooldown.blocked(&[atom]) {
continue;
}
flatten_atoms.insert(atom);
report.proposals.push(proposal(
&term,
StructureMove::Death { atom },
f64::MAX / 4.0,
));
}
}
}
if report.proposals.is_empty() || split.shards.is_empty() {
rounds.push(SearchLedger {
alpha: budget.alpha,
moves: Vec::new(),
collapse_events: term.collapse_events().to_vec(),
});
round_predictions.push(birth_predictions);
break;
}
type State = (SaeManifoldTerm, SaeManifoldRho);
let collapse_events = term.collapse_events().to_vec();
let decoders = birth_seeds;
let estimation_rows = split.estimation_rows.clone();
let certified_glues = std::mem::take(&mut report.certified_glues);
let proposals = std::mem::take(&mut report.proposals);
let outcome: SearchOutcome<State> = search(
(term, rho),
proposals,
&split.shards,
&budget,
ledger,
|state: &State, mv: &StructureMove| {
if matches!(mv, StructureMove::Glue { .. }) {
return Ok(state.clone());
}
let (cand_term, cand_rho) =
apply_structure_move_seeded(&state.0, &state.1, mv, &decoders)?;
Ok(candidate_fit(cand_term, cand_rho, &estimation_rows))
},
|state: &State, shard: &RowBlockShard| eval_log_lik(&state.0, shard),
|state: &State, shard: &RowBlockShard| eval_log_lik(&state.0, shard),
|state: State, _: &RowBlockShard| state,
)?;
let (next_term, next_rho) = outcome.state;
let mut round_ledger = outcome.ledger;
round_ledger.collapse_events = collapse_events;
let applied = round_ledger.moves.iter().any(|m| {
matches!(
m.verdict,
gam_solve::structure_search::MoveVerdict::Accepted { .. }
| gam_solve::structure_search::MoveVerdict::Demoted { .. }
)
});
let requires_polish = round_ledger.moves.iter().any(|record| {
let fired = matches!(
record.verdict,
gam_solve::structure_search::MoveVerdict::Accepted { .. }
| gam_solve::structure_search::MoveVerdict::Demoted { .. }
);
fired && !matches!(record.mv, StructureMove::Glue { .. })
});
if let Some(cfg) = curl {
use gam_solve::structure_search::MoveVerdict;
for rec in &round_ledger.moves {
let fired = matches!(
rec.verdict,
MoveVerdict::Accepted { .. } | MoveVerdict::Demoted { .. }
);
if !fired {
continue;
}
match &rec.mv {
StructureMove::Birth { candidate } => {
if let Some(members) = curl_atoms.get(candidate) {
cooldown.record(members, cfg.cooldown_rounds);
}
}
StructureMove::Death { atom } if flatten_atoms.contains(atom) => {
cooldown.record(&[*atom], cfg.cooldown_rounds);
}
_ => {}
}
}
cooldown.tick();
}
rounds.push(round_ledger);
round_predictions.push(birth_predictions);
if applied {
let (mut next_term, mut next_rho) = (next_term, next_rho);
compact_glued_atoms(
&mut next_term,
&mut next_rho,
rounds.last().expect("round ledger pushed above"),
&certified_glues,
)?;
if requires_polish {
let (mut polished_term, polished_rho) =
finalize_round(next_term, next_rho, &split.estimation_rows);
refresh_registered_atlas_transitions(&mut polished_term)?;
term = polished_term;
rho = polished_rho;
} else {
term = next_term;
rho = next_rho;
}
} else {
term = next_term;
rho = next_rho;
break;
}
}
Ok(StructureSearchResult::from_rounds_with_predictions(
term,
rho,
rounds,
&round_predictions,
))
}
fn build_birth_decoders(
term: &SaeManifoldTerm,
residuals: ArrayView2<'_, f64>,
params: &HarvestParams,
) -> Result<Vec<Array2<f64>>, String> {
let n = residuals.nrows();
let p = residuals.ncols();
if params.max_births == 0 || n == 0 || p == 0 {
return Ok(Vec::new());
}
let assignments = term.assignment.assignments();
let activity: Array1<f64> = (0..n).map(|r| assignments.row(r).sum()).collect();
let max_rank = params.max_births.min(p.saturating_sub(1));
let model = StructuredResidualModel::fit(ResidualFactorInput {
residuals,
activity: activity.view(),
max_factor_rank: max_rank,
})
.map_err(|e| format!("build_birth_decoders: structured-residual fit failed: {e}"))?;
let factor = model.factor();
let r = factor.ncols();
let m = term.atoms[0].basis_size();
let mut decoders = Vec::with_capacity(r);
for j in 0..r {
let mut decoder = Array2::<f64>::zeros((m, p));
for out in 0..p {
decoder[[0, out]] = factor[[out, j]];
}
decoders.push(decoder);
}
Ok(decoders)
}
struct CurlCandidate {
members: Vec<usize>,
seed: BirthSeed,
net_evidence: f64,
}
fn is_linear_like(kind: &SaeAtomBasisKind) -> bool {
matches!(
kind,
SaeAtomBasisKind::Linear | SaeAtomBasisKind::EuclideanPatch
)
}
fn atom_ambient_image(atom: &SaeManifoldAtom) -> Array2<f64> {
atom.basis_values.dot(&atom.decoder_coefficients)
}
fn power_iter_top_dir(
img: ArrayView2<'_, f64>,
center: &Array1<f64>,
active: &[usize],
) -> Array1<f64> {
let p = img.ncols();
let mut v = Array1::<f64>::zeros(p);
let mut best_norm = 0.0_f64;
for &r in active {
let mut nrm = 0.0_f64;
for j in 0..p {
let d = img[[r, j]] - center[j];
nrm += d * d;
}
if nrm > best_norm {
best_norm = nrm;
for j in 0..p {
v[j] = img[[r, j]] - center[j];
}
}
}
let vn = v.dot(&v).sqrt();
if vn <= 0.0 {
return v;
}
v.mapv_inplace(|x| x / vn);
for _ in 0..5 {
let mut w = Array1::<f64>::zeros(p);
for &r in active {
let mut dot = 0.0_f64;
for j in 0..p {
dot += (img[[r, j]] - center[j]) * v[j];
}
for j in 0..p {
w[j] += (img[[r, j]] - center[j]) * dot;
}
}
let wn = w.dot(&w).sqrt();
if wn <= 0.0 {
break;
}
w.mapv_inplace(|x| x / wn);
v = w;
}
v
}
fn linear_atom_frames(term: &SaeManifoldTerm) -> Vec<(usize, Array1<f64>, Vec<bool>, Array2<f64>)> {
let assignments = term.assignment.assignments();
let n = assignments.nrows();
let k = assignments.ncols();
let floor = if k == 0 {
0.0
} else {
ACTIVE_SUPPORT_REL_FLOOR / k as f64
};
let mut out = Vec::new();
for (a, atom) in term.atoms.iter().enumerate() {
if !is_linear_like(&atom.basis_kind) {
continue;
}
let active_mask: Vec<bool> = (0..n).map(|r| assignments[[r, a]] > floor).collect();
let active_idx: Vec<usize> = (0..n).filter(|&r| active_mask[r]).collect();
if active_idx.len() < 2 {
continue;
}
let img = atom_ambient_image(atom);
if img.ncols() == 0 {
continue;
}
let p = img.ncols();
let mut center = Array1::<f64>::zeros(p);
for &r in &active_idx {
for j in 0..p {
center[j] += img[[r, j]];
}
}
center.mapv_inplace(|x| x / active_idx.len() as f64);
let dir = power_iter_top_dir(img.view(), ¢er, &active_idx);
if dir.dot(&dir).sqrt() <= 0.0 {
continue;
}
out.push((a, dir, active_mask, img));
}
out
}
fn curl_candidates(
term: &SaeManifoldTerm,
residuals: ArrayView2<'_, f64>,
cfg: &CurlConfig,
) -> Vec<CurlCandidate> {
let frames = linear_atom_frames(term);
if frames.len() < 2 {
return Vec::new();
}
let n = term.assignment.logits.nrows();
let p = term.output_dim();
let mut sse = 0.0_f64;
let mut cnt = 0usize;
for r in 0..residuals.nrows() {
for j in 0..residuals.ncols() {
sse += residuals[[r, j]] * residuals[[r, j]];
cnt += 1;
}
}
let sigma = if cnt > 0 {
(sse / cnt as f64).sqrt().max(1e-9)
} else {
1e-9
};
let dirs: Vec<ArrayView1<f64>> = frames.iter().map(|(_, d, _, _)| d.view()).collect();
let actives: Vec<Vec<bool>> = frames.iter().map(|(_, _, m, _)| m.clone()).collect();
let ids: Vec<usize> = frames.iter().map(|(a, _, _, _)| *a).collect();
let signed = crate::manifold::coalesce_antipodal(
&dirs,
&actives,
&ids,
cfg.coalesce_cos_threshold,
cfg.coalesce_max_overlap,
);
if signed.len() < 2 {
return Vec::new();
}
let frame_of: std::collections::HashMap<usize, usize> =
ids.iter().enumerate().map(|(i, a)| (*a, i)).collect();
let signed_active: Vec<Vec<bool>> = signed.iter().map(|s| s.active.clone()).collect();
let rows: Vec<usize> = if n <= cfg.subsample_rows {
(0..n).collect()
} else {
let stride = n / cfg.subsample_rows;
(0..n).step_by(stride.max(1)).collect()
};
let pairs = crate::manifold::cooccurrence_pairs(&signed_active, &rows, cfg.min_cooccurrence);
let mut cands: Vec<CurlCandidate> = Vec::new();
for (si, sj, _count) in pairs {
let di = &signed[si];
let dj = &signed[sj];
let mut co_fire: Vec<usize> = (0..n)
.filter(|&r| {
di.active.get(r).copied().unwrap_or(false)
&& dj.active.get(r).copied().unwrap_or(false)
})
.collect();
if co_fire.len() < cfg.min_cooccurrence.max(2) {
continue;
}
if co_fire.len() > cfg.subsample_rows {
let stride = (co_fire.len() / cfg.subsample_rows).max(1);
co_fire = co_fire.iter().copied().step_by(stride).collect();
}
let members: Vec<usize> = di
.members
.iter()
.chain(dj.members.iter())
.copied()
.collect();
let mut x = Array2::<f64>::zeros((co_fire.len(), p));
for (row_out, &r) in co_fire.iter().enumerate() {
for &atom in &members {
if let Some(&fi) = frame_of.get(&atom) {
let img = &frames[fi].3;
for j in 0..p {
x[[row_out, j]] += img[[r, j]];
}
}
}
}
let mut center = Array1::<f64>::zeros(p);
for row_out in 0..co_fire.len() {
for j in 0..p {
center[j] += x[[row_out, j]];
}
}
center.mapv_inplace(|v| v / co_fire.len() as f64);
let (alpha, beta, e1, e2) = match crate::manifold::orthonormal_pair_coords(
x.view(),
di.dir.view(),
dj.dir.view(),
center.view(),
) {
Ok(t) => t,
Err(_) => continue,
};
let n_eff = co_fire.len() as f64;
let m_circle = (2 * cfg.harmonics + 1) as f64;
let delta_charge = 0.5 * m_circle * n_eff.max(2.0).ln();
let verdict = match crate::manifold::curl_verdict(
alpha.view(),
beta.view(),
sigma,
n_eff,
delta_charge,
) {
Ok(v) => v,
Err(_) => continue,
};
if !verdict.recommend_curl {
continue;
}
let seed_circle = match crate::manifold::curl_seed(
e1.view(),
e2.view(),
alpha.view(),
beta.view(),
cfg.harmonics,
center.view(),
) {
Ok(s) => s,
Err(_) => continue,
};
let mut phase_coords = Array2::<f64>::zeros((n, 1));
let mut gate = vec![f64::NEG_INFINITY; n];
let own = verdict.gain_nats_per_row.max(0.5);
for (idx, &r) in co_fire.iter().enumerate() {
phase_coords[[r, 0]] = seed_circle.theta_turns[idx];
gate[r] = own;
}
cands.push(CurlCandidate {
members,
seed: BirthSeed::Circle {
decoder: seed_circle.decoder,
phase_coords,
gate,
},
net_evidence: verdict.net_evidence_nats,
});
}
cands.sort_by(|a, b| b.net_evidence.total_cmp(&a.net_evidence));
let mut claimed: std::collections::HashSet<usize> = std::collections::HashSet::new();
let mut out = Vec::new();
for c in cands {
if c.members.iter().any(|a| claimed.contains(a)) {
continue;
}
for a in &c.members {
claimed.insert(*a);
}
out.push(c);
if out.len() >= cfg.max_curls {
break;
}
}
out
}
fn flatten_candidates(term: &SaeManifoldTerm) -> Vec<usize> {
let assignments = term.assignment.assignments();
let n = assignments.nrows();
let k = assignments.ncols();
let floor = if k == 0 {
0.0
} else {
ACTIVE_SUPPORT_REL_FLOOR / k as f64
};
let mut out = Vec::new();
for (a, atom) in term.atoms.iter().enumerate() {
if !matches!(atom.basis_kind, SaeAtomBasisKind::Periodic) || atom.latent_dim != 1 {
continue;
}
let active_idx: Vec<usize> = (0..n).filter(|&r| assignments[[r, a]] > floor).collect();
if active_idx.len() < 8 {
continue;
}
let img = atom_ambient_image(atom);
let p = img.ncols();
let mut center = Array1::<f64>::zeros(p);
for &r in &active_idx {
for j in 0..p {
center[j] += img[[r, j]];
}
}
center.mapv_inplace(|x| x / active_idx.len() as f64);
let coords = term.assignment.coords[a].as_matrix();
if coords.ncols() == 0 {
continue;
}
let mut radii = Array1::<f64>::zeros(active_idx.len());
let mut angles = Array1::<f64>::zeros(active_idx.len());
for (i, &r) in active_idx.iter().enumerate() {
let mut rr = 0.0_f64;
for j in 0..p {
let d = img[[r, j]] - center[j];
rr += d * d;
}
radii[i] = rr.sqrt();
angles[i] = std::f64::consts::TAU * coords[[r, 0]];
}
if let Ok(v) = crate::manifold::flatten_verdict(radii.view(), angles.view()) {
if v.recommend_flatten {
out.push(a);
}
}
}
out
}
fn eval_log_lik(term: &SaeManifoldTerm, shard: &RowBlockShard) -> f64 {
let fitted = match term.try_fitted_target_aware(shard.target.view(), None) {
Ok(f) => f,
Err(_) => return f64::NEG_INFINITY,
};
let n_full = fitted.nrows();
let p = fitted.ncols();
if p != shard.target.ncols() || n_full != shard.target.nrows() {
return f64::NEG_INFINITY;
}
let mut sse = 0.0_f64;
let mut count = 0usize;
for &row in &shard.rows {
if row >= n_full {
continue;
}
for out in 0..p {
let d = fitted[[row, out]] - shard.target[[row, out]];
sse_accumulate(&mut sse, d);
}
count += p;
}
if count == 0 {
return f64::NEG_INFINITY;
}
let reconstruction = -0.5 * sse;
let gate_evidence = gate_block_log_evidence(term, shard);
reconstruction + gate_evidence
}
fn gate_block_log_evidence(term: &SaeManifoldTerm, shard: &RowBlockShard) -> f64 {
use gam_solve::inference::pg_gate_evidence::{GateBlock, pg_gate_evidence};
let logits = &term.assignment.logits;
let n_full = logits.nrows();
let k = logits.ncols();
if k == 0 {
return 0.0;
}
let rows: Vec<usize> = shard.rows.iter().copied().filter(|&r| r < n_full).collect();
let m = rows.len();
if m == 0 {
return 0.0;
}
let design = Array2::<f64>::ones((m, 1));
let b = Array1::<f64>::ones(m);
let penalty = Array2::<f64>::eye(1);
let mut total = 0.0_f64;
for atom in 0..k {
let mut psi = Array1::<f64>::zeros(m);
let mut y = Array1::<f64>::zeros(m);
for (i, &row) in rows.iter().enumerate() {
let logit = logits[[row, atom]];
if !logit.is_finite() {
return 0.0;
}
psi[i] = logit;
y[i] = if logit > 0.0 { 1.0 } else { 0.0 };
}
let block = GateBlock {
design: design.view(),
y: y.view(),
b: b.view(),
offset: None,
psi_hat: Some(psi.view()),
penalty: Some(penalty.view()),
hess_rest: None,
h_rest: None,
};
match pg_gate_evidence(&block) {
Ok(ev) => total -= ev.neg_log_evidence,
Err(_) => return 0.0,
}
}
total
}
#[inline]
fn sse_accumulate(sse: &mut f64, d: f64) {
*sse += d * d;
}
#[derive(Clone, Copy, Debug)]
pub struct ProductionRefitParams {
pub inner_max_iter: usize,
pub scoring_inner_max_iter: usize,
pub learning_rate: f64,
pub ridge_ext_coord: f64,
pub ridge_beta: f64,
}
pub fn run_production_structure_search(
term: SaeManifoldTerm,
rho: SaeManifoldRho,
target: ArrayView2<'_, f64>,
config: RoundDriverConfig,
refit_params: ProductionRefitParams,
ledger: &mut StructureLedger,
) -> Result<StructureSearchResult, String> {
let n = target.nrows();
let refit_at = |full_target: ArrayView2<'_, f64>,
mut cand_term: SaeManifoldTerm,
mut cand_rho: SaeManifoldRho,
estimation_rows: &[usize],
inner_max_iter: usize|
-> (SaeManifoldTerm, SaeManifoldRho) {
const HELD_OUT_WEIGHT: f64 = 1e-12;
let mut weights = vec![HELD_OUT_WEIGHT; n];
for &r in estimation_rows {
if r < n {
weights[r] = 1.0;
}
}
if cand_term.set_row_loss_weights(weights).is_err() {
return (cand_term, cand_rho);
}
if cand_term
.run_joint_fit_arrow_schur(
full_target,
&mut cand_rho,
None,
inner_max_iter,
refit_params.learning_rate,
refit_params.ridge_ext_coord,
refit_params.ridge_beta,
)
.is_err()
{
return (cand_term, cand_rho);
}
(cand_term, cand_rho)
};
let scoring_iters = refit_params.scoring_inner_max_iter;
let full_iters = refit_params.inner_max_iter;
let full_target_score = target.to_owned();
let full_target_polish = target.to_owned();
run_structure_search_rounds(
term,
rho,
target,
config,
ledger,
move |cand_term, cand_rho, estimation_rows| {
refit_at(
full_target_score.view(),
cand_term,
cand_rho,
estimation_rows,
scoring_iters,
)
},
move |adopted_term, adopted_rho, _estimation_rows| {
let all_rows: Vec<usize> = (0..n).collect();
refit_at(
full_target_polish.view(),
adopted_term,
adopted_rho,
&all_rows,
full_iters,
)
},
)
}
pub fn rounds_to_json(rounds: &[SearchLedger]) -> Result<String, String> {
serde_json::to_string(rounds)
.map_err(|e| format!("rounds_to_json: serialize search ledger: {e}"))
}
#[cfg(test)]
mod tests {
use super::*;
use crate::manifold::{
AssignmentMode, PeriodicHarmonicEvaluator, SaeAssignment, SaeAtomBasisKind,
SaeBasisEvaluator, SaeManifoldAtom,
};
use gam_solve::structure_search::{CollapseAction, CollapseEvent};
use gam_terms::latent::LatentManifold;
use ndarray::Array2;
use std::sync::Arc;
#[test]
fn dedup_most_suspect_keeps_one_per_parent() {
let raw = vec![
(2usize, 0.4_f64),
(5, 0.9),
(1, 0.6),
(2, 0.1),
(5, 0.3),
(2, 0.7),
];
let out = dedup_most_suspect_per_parent(raw);
assert_eq!(out.len(), 3, "one entry per distinct parent: {out:?}");
let mut atoms: Vec<usize> = out.iter().map(|(a, _)| *a).collect();
atoms.sort_unstable();
assert_eq!(atoms, vec![1, 2, 5], "all distinct parents kept");
let sig = |atom: usize| out.iter().find(|(a, _)| *a == atom).unwrap().1;
assert_eq!(sig(2), 0.1, "atom 2 keeps its most-suspect nomination");
assert_eq!(sig(5), 0.3, "atom 5 keeps its most-suspect nomination");
assert_eq!(sig(1), 0.6, "the singly-nominated atom is unchanged");
assert_eq!(
out,
vec![(2, 0.1), (5, 0.3), (1, 0.6)],
"deterministic most-suspect-first order"
);
}
const ON: f64 = 6.0;
const OFF: f64 = -6.0;
#[test]
fn auto_primary_topology_never_falls_back_on_race_failure_2238_2239() {
let target = Array2::<f64>::zeros((15, 2));
let labels = vec![0usize; 15];
let mut basis = vec!["auto".to_string()];
let mut dims = vec![2usize];
let error = resolve_auto_primary_atoms(target.view(), &labels, &mut basis, &mut dims)
.expect_err("an undersupported automatic race must be rejected");
assert!(error.contains("auto atom 0"), "unexpected error: {error}");
assert!(error.contains("at least 16"), "unexpected error: {error}");
assert_eq!(basis, vec!["auto"], "failure must not install a fallback");
assert_eq!(dims, vec![2], "failure must not rewrite latent dimension");
}
#[test]
fn auto_primary_topology_selects_two_dimensional_factor_2238() {
let side = 8usize;
let target = Array2::<f64>::from_shape_fn((side * side, 2), |(row, col)| {
let i = row / side;
let j = row % side;
if col == 0 {
i as f64 - 0.5 * (side - 1) as f64
} else {
j as f64 - 0.5 * (side - 1) as f64
}
});
let labels = vec![0usize; target.nrows()];
let choices = discover_primary_atom_topologies(target.view(), &labels, 1, &[2])
.expect("the supported planar race must produce a winner");
assert_eq!(choices.len(), 1);
assert_eq!(choices[0].latent_dim, 2);
assert_eq!(choices[0].basis_kind, SaeAtomBasisKind::EuclideanPatch);
}
#[test]
fn auto_primary_topology_selects_curved_sphere_and_beats_circle_2238_2239() {
use crate::basis::SphereChartEvaluator;
use gam_solve::AutoTopologyKind;
use ndarray::Array1;
let (n_lat, n_lon) = (12usize, 14usize);
let n = n_lat * n_lon;
let mut lat = Vec::with_capacity(n);
let mut lon = Vec::with_capacity(n);
let mut target = Array2::<f64>::zeros((n, 3));
for i in 0..n_lat {
let theta = -std::f64::consts::FRAC_PI_2
+ std::f64::consts::PI * (i as f64 + 1.0) / (n_lat as f64 + 1.0);
for j in 0..n_lon {
let phi = std::f64::consts::TAU * j as f64 / n_lon as f64;
let row = i * n_lon + j;
target[[row, 0]] = theta.cos() * phi.cos();
target[[row, 1]] = theta.cos() * phi.sin();
target[[row, 2]] = theta.sin();
lat.push(theta);
lon.push(phi);
}
}
let labels = vec![0usize; n];
let choices = discover_primary_atom_topologies(target.view(), &labels, 1, &[2])
.expect("the supported sphere race must produce a winner");
assert_eq!(choices.len(), 1);
assert_eq!(
choices[0].basis_kind,
SaeAtomBasisKind::Sphere,
"the curved 2-sphere factor must be discovered as a sphere chart, not a circle/patch"
);
assert_eq!(choices[0].latent_dim, 2, "a sphere is intrinsically 2-D");
let weights = Array1::<f64>::ones(n);
let recon_r2 = |spec: &TopologyCandidateSpec| -> f64 {
let fit = fit_topology_candidate(spec, target.view(), weights.view())
.expect("candidate fit")
.fit_handle;
let recon = fit.phi.dot(&fit.decoder);
let mut means = [0.0_f64; 3];
for col in 0..3 {
let mut acc = 0.0;
for row in 0..n {
acc += target[[row, col]];
}
means[col] = acc / n as f64;
}
let (mut ss_res, mut ss_tot) = (0.0_f64, 0.0_f64);
for row in 0..n {
for col in 0..3 {
let r = target[[row, col]] - recon[[row, col]];
ss_res += r * r;
let c = target[[row, col]] - means[col];
ss_tot += c * c;
}
}
1.0 - ss_res / ss_tot.max(1e-12)
};
let mut circle_coords = Array2::<f64>::zeros((n, 1));
for row in 0..n {
circle_coords[[row, 0]] = lon[row] / std::f64::consts::TAU;
}
let circle_spec = TopologyCandidateSpec {
kind: AutoTopologyKind::Circle,
basis_kind: SaeAtomBasisKind::Periodic,
manifold: LatentManifold::Circle { period: 1.0 },
latent_dim: 1,
evaluator: Arc::new(PeriodicHarmonicEvaluator::new(3).expect("periodic evaluator")),
coords: circle_coords,
};
let mut sphere_coords = Array2::<f64>::zeros((n, 2));
for row in 0..n {
sphere_coords[[row, 0]] = lat[row];
sphere_coords[[row, 1]] = lon[row];
}
let sphere_spec = TopologyCandidateSpec {
kind: AutoTopologyKind::Sphere,
basis_kind: SaeAtomBasisKind::Sphere,
manifold: LatentManifold::Product(vec![
LatentManifold::Interval {
lo: -std::f64::consts::FRAC_PI_2,
hi: std::f64::consts::FRAC_PI_2,
},
LatentManifold::Circle {
period: std::f64::consts::TAU,
},
]),
latent_dim: 2,
evaluator: Arc::new(SphereChartEvaluator),
coords: sphere_coords,
};
let circle_r2 = recon_r2(&circle_spec);
let sphere_r2 = recon_r2(&sphere_spec);
eprintln!(
"[topology-2238] planted 2-sphere R²: circle-pinned={circle_r2:.4} sphere-chart={sphere_r2:.4}"
);
assert!(
sphere_r2 > 0.9,
"the discovered sphere chart must recover the planted 2-sphere (R²={sphere_r2:.4})"
);
assert!(
circle_r2 < 0.75,
"the 1-D circle default structurally caps the 2-sphere recovery (R²={circle_r2:.4})"
);
assert!(
sphere_r2 > circle_r2 + 0.2,
"discovery must strictly beat the circle-pinned recovery (sphere={sphere_r2:.4} vs circle={circle_r2:.4})"
);
}
#[test]
fn select_periodic_resolution_grows_past_default_over_harmonic_gap_2243() {
use gam_solve::AutoTopologyKind;
use ndarray::Array1;
let n = 240usize;
let mut coords = Array2::<f64>::zeros((n, 1));
let mut target = Array2::<f64>::zeros((n, 4));
for row in 0..n {
let t = row as f64 / n as f64;
let angle = std::f64::consts::TAU * t;
coords[[row, 0]] = t;
target[[row, 0]] = angle.cos();
target[[row, 1]] = angle.sin();
target[[row, 2]] = (4.0 * angle).cos();
target[[row, 3]] = (4.0 * angle).sin();
}
let weights = Array1::<f64>::ones(n);
let selected = select_periodic_resolution(coords.view(), target.view(), weights.view(), n)
.expect("resolution selection must succeed on a supported periodic signal");
assert!(
selected >= 4,
"the 4th-harmonic content (past a gap) requires at least 4 harmonics; the fixed \
2-harmonic default under-resolves it (selected={selected})"
);
let circle_r2 = |h: usize| -> f64 {
let spec = TopologyCandidateSpec {
kind: AutoTopologyKind::Circle,
basis_kind: SaeAtomBasisKind::Periodic,
manifold: LatentManifold::Circle { period: 1.0 },
latent_dim: 1,
evaluator: Arc::new(PeriodicHarmonicEvaluator::new(h).expect("periodic evaluator")),
coords: coords.clone(),
};
let fit = fit_topology_candidate(&spec, target.view(), weights.view())
.expect("candidate fit")
.fit_handle;
let recon = fit.phi.dot(&fit.decoder);
let (mut ss_res, mut ss_tot) = (0.0_f64, 0.0_f64);
for col in 0..4 {
let mut mean = 0.0;
for row in 0..n {
mean += target[[row, col]];
}
mean /= n as f64;
for row in 0..n {
let r = target[[row, col]] - recon[[row, col]];
ss_res += r * r;
let c = target[[row, col]] - mean;
ss_tot += c * c;
}
}
1.0 - ss_res / ss_tot.max(1e-12)
};
let default_r2 = circle_r2(2);
let selected_r2 = circle_r2(selected);
eprintln!(
"[resolution-2243] circle R²: default(2 harmonics)={default_r2:.4} selected({selected})={selected_r2:.4}"
);
assert!(
selected_r2 > 0.99,
"the evidence-selected resolution must recover the signal (R²={selected_r2:.4})"
);
assert!(
default_r2 < 0.75,
"the 2-harmonic default cannot represent the 4th-harmonic half of the energy (R²={default_r2:.4})"
);
}
fn vdc(n: usize) -> Vec<f64> {
(0..n)
.map(|i| {
let (mut x, mut denom, mut k) = (0.0_f64, 2.0_f64, i + 1);
while k > 0 {
x += (k & 1) as f64 / denom;
denom *= 2.0;
k >>= 1;
}
x
})
.collect()
}
#[test]
fn radial_promotion_fires_only_on_continuous_amplitude() {
let n = 400;
let coords = Array2::from_shape_fn((n, 1), |(i, _)| i as f64 / n as f64);
let u = vdc(n);
let disk = Array2::from_shape_fn((n, 2), |(i, j)| {
let r = u[i].sqrt();
let theta = std::f64::consts::TAU * (i as f64 / n as f64);
if j == 0 {
r * theta.cos()
} else {
r * theta.sin()
}
});
let promoted = radial_promoted_specs(coords.view(), disk.view(), 1)
.expect("promotion decision")
.expect("disk amplitude is continuous ⇒ promotion fires");
let kinds: std::collections::HashSet<_> = promoted.iter().map(|s| s.kind).collect();
assert!(kinds.contains(&AutoTopologyKind::Circle), "{kinds:?}");
assert!(kinds.contains(&AutoTopologyKind::Cylinder), "{kinds:?}");
assert!(kinds.contains(&AutoTopologyKind::Euclidean), "{kinds:?}");
assert_eq!(kinds.len(), promoted.len());
let expected_radial =
standardized_log_birth_amplitudes(birth_row_amplitudes(disk.view()).view())
.expect("disk log-amplitude spread");
for spec in promoted.iter().filter(|spec| {
matches!(
spec.kind,
AutoTopologyKind::Cylinder | AutoTopologyKind::Euclidean
)
}) {
for row in 0..n {
assert!(
(spec.coords[[row, 1]] - expected_radial[row]).abs() < 1.0e-12,
"promoted {:?} row {row} axis 1 must be standardized log-amplitude",
spec.kind
);
}
}
let ring = Array2::from_shape_fn((n, 2), |(i, j)| {
if i % 2 == 0 {
0.0
} else {
let theta = std::f64::consts::TAU * (i as f64 / n as f64);
if j == 0 { theta.cos() } else { theta.sin() }
}
});
assert!(
radial_promoted_specs(coords.view(), ring.view(), 1)
.expect("promotion decision")
.is_none(),
"present/absent birth must not promote"
);
assert!(
radial_promoted_specs(coords.view(), disk.view(), 2)
.expect("promotion decision")
.is_none()
);
}
fn topology_fit_sse(fit: &TopologyRaceFit, target: ArrayView2<'_, f64>) -> f64 {
let fitted = fit.phi.dot(&fit.decoder);
let mut sse = 0.0_f64;
for row in 0..target.nrows() {
for col in 0..target.ncols() {
let err = target[[row, col]] - fitted[[row, col]];
sse += err * err;
}
}
sse
}
#[test]
fn radial_promotion_seed_coordinate_expresses_annulus_radius() {
let n_angles = 16;
let n_radii = 25;
let n = n_angles * n_radii;
let radial = vdc(n_radii);
let coords = Array2::from_shape_fn((n, 1), |(row, _)| {
let angle_idx = row / n_radii;
angle_idx as f64 / n_angles as f64
});
let annulus = Array2::from_shape_fn((n, 2), |(row, col)| {
let angle_idx = row / n_radii;
let radius_idx = row % n_radii;
let theta = std::f64::consts::TAU * (angle_idx as f64 / n_angles as f64);
let radius = 0.3 + 0.7 * radial[radius_idx].sqrt();
if col == 0 {
radius * theta.cos()
} else {
radius * theta.sin()
}
});
let promoted = radial_promoted_specs(coords.view(), annulus.view(), 1)
.expect("promotion decision")
.expect("annulus radius spread promotes a radial axis");
let circle = promoted
.iter()
.find(|spec| spec.kind == AutoTopologyKind::Circle)
.expect("promoted race includes the circle alternative");
let cylinder = promoted
.iter()
.find(|spec| spec.kind == AutoTopologyKind::Cylinder)
.expect("promoted race includes the cylinder alternative");
let weights = Array1::<f64>::ones(n);
let circle_fit =
fit_topology_candidate(circle, annulus.view(), weights.view()).expect("circle fit");
let cylinder_fit =
fit_topology_candidate(cylinder, annulus.view(), weights.view()).expect("cylinder fit");
let circle_sse = topology_fit_sse(&circle_fit.fit_handle, annulus.view());
let cylinder_sse = topology_fit_sse(&cylinder_fit.fit_handle, annulus.view());
assert!(
cylinder_sse < 0.75 * circle_sse,
"radial seed should let cylinder express radius variation: cylinder_sse={cylinder_sse}, circle_sse={circle_sse}"
);
}
#[test]
fn birth_row_amplitudes_are_row_norms() {
let y = Array2::from_shape_vec((2, 2), vec![3.0, 4.0, 0.0, 0.0]).unwrap();
let a = birth_row_amplitudes(y.view());
assert!((a[0] - 5.0).abs() < 1e-12);
assert!((a[1]).abs() < 1e-12);
}
#[test]
fn finite_set_race_is_not_enrolled_by_default() {
assert!(!finite_set_race_enrolled());
set_finite_set_race_enrolled(true);
assert!(finite_set_race_enrolled());
set_finite_set_race_enrolled(false);
assert!(!finite_set_race_enrolled());
}
#[test]
fn finite_set_candidate_fires_on_discrete_occupancy() {
let per = 100;
let mut rows = Vec::new();
for i in 0..(7 * per) {
rows.push((i % 7) as f64 + 0.001 * ((i as f64).sin()));
}
let coords = Array2::from_shape_vec((7 * per, 1), rows).unwrap();
let (anchors, idx) =
finite_set_candidate_for_birth(coords.view()).expect("discrete ⇒ finite-set candidate");
assert_eq!(anchors, 7, "anchors");
assert_eq!(crate::manifold::finite_set_rank_charge(anchors), 6);
assert!(
idx.iter()
.all(|&v| (0.0..=6.0).contains(&v) && v.fract() == 0.0)
);
let n = 400;
let uni = Array2::from_shape_fn((n, 1), |(i, _)| i as f64 / n as f64);
assert!(finite_set_candidate_for_birth(uni.view()).is_none());
}
#[test]
fn anchor_indicator_evaluator_is_one_hot_with_zero_jets() {
use crate::basis::{AnchorIndicatorEvaluator, SaeBasisEvaluator, SaeBasisSecondJet};
let ev = AnchorIndicatorEvaluator::new(3).unwrap();
let coords = Array2::from_shape_vec((4, 1), vec![0.0, 1.0, 2.0, 1.4]).unwrap();
let (phi, jet) = ev.evaluate(coords.view()).unwrap();
assert_eq!(phi.dim(), (4, 3));
for r in 0..4 {
assert!((phi.row(r).sum() - 1.0).abs() < 1e-12);
}
assert!((phi[[0, 0]] - 1.0).abs() < 1e-12);
assert!((phi[[1, 1]] - 1.0).abs() < 1e-12);
assert!((phi[[2, 2]] - 1.0).abs() < 1e-12);
assert!((phi[[3, 1]] - 1.0).abs() < 1e-12); assert!(jet.iter().all(|&v| v == 0.0));
let h = ev.second_jet(coords.view()).unwrap();
assert!(h.iter().all(|&v| v == 0.0));
}
fn planted_term(active: &[Vec<bool>]) -> (SaeManifoldTerm, SaeManifoldRho) {
let n = active.len();
let k = active[0].len();
let p = 4usize;
let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(3).unwrap());
let coords = Array2::<f64>::from_shape_fn((n, 1), |(row, _)| row as f64 / n as f64);
let (phi, jet) = evaluator.evaluate(coords.view()).unwrap();
let mut atoms = Vec::with_capacity(k);
let mut coord_blocks = Vec::with_capacity(k);
for atom_idx in 0..k {
let mut decoder = Array2::<f64>::zeros((3, p));
decoder[[1, atom_idx % p]] = 1.0;
decoder[[2, (atom_idx + 1) % p]] = 1.0;
let atom = SaeManifoldAtom::new(
format!("atom_{atom_idx}"),
SaeAtomBasisKind::Periodic,
1,
phi.clone(),
jet.clone(),
decoder,
Array2::<f64>::eye(3),
)
.unwrap()
.with_basis_second_jet(evaluator.clone());
atoms.push(atom);
coord_blocks.push(coords.clone());
}
let mut logits = Array2::<f64>::zeros((n, k));
for (row, atom_active) in active.iter().enumerate() {
for (atom, &on) in atom_active.iter().enumerate() {
logits[[row, atom]] = if on { ON } else { OFF };
}
}
let assignment = SaeAssignment::from_blocks_with_mode_and_manifolds(
logits,
coord_blocks,
vec![LatentManifold::Circle { period: 1.0 }; k],
AssignmentMode::softmax(1.0),
)
.unwrap();
let term = SaeManifoldTerm::new(atoms, assignment).unwrap();
let rho = SaeManifoldRho::new(0.0, 0.0, vec![Array1::<f64>::zeros(1); k]);
(term, rho)
}
fn residuals_of(term: &SaeManifoldTerm) -> Array2<f64> {
let fitted = term.try_fitted().unwrap();
-&fitted
}
#[test]
fn structure_changed_is_true_only_when_a_move_lands() {
use gam_solve::structure_search::{MoveRecord, MoveVerdict};
fn ledger_with(verdicts: Vec<MoveVerdict>) -> SearchLedger {
SearchLedger {
alpha: 0.05,
moves: verdicts
.into_iter()
.enumerate()
.map(|(i, verdict)| MoveRecord {
mv: StructureMove::Death { atom: i },
trigger: 0.0,
structure_hash: i as u64,
claim: ClaimKind::AtomExists { atom: i },
verdict,
})
.collect(),
collapse_events: Vec::new(),
}
}
let (term0, rho0) = planted_term(&[vec![true], vec![true]]);
let empty = StructureSearchResult::from_rounds(term0.clone(), rho0.clone(), Vec::new());
assert!(
!empty.structure_changed(),
"no rounds ⇒ the term/rho are the pre-search fit ⇒ structure_changed() must be false"
);
let no_landed = StructureSearchResult::from_rounds(
term0.clone(),
rho0.clone(),
vec![ledger_with(vec![
MoveVerdict::Contested { log_e: -1.0 },
MoveVerdict::Vetoed { log_e: -2.0 },
])],
);
assert!(
!no_landed.structure_changed(),
"all-contested/vetoed rounds leave the model unchanged ⇒ structure_changed() must be false"
);
let accepted = StructureSearchResult::from_rounds(
term0.clone(),
rho0.clone(),
vec![ledger_with(vec![
MoveVerdict::Contested { log_e: -1.0 },
MoveVerdict::Accepted { log_e: 3.0 },
])],
);
assert!(
accepted.structure_changed(),
"a landed Accepted move mutates term/rho ⇒ structure_changed() must be true (recompute bands)"
);
let demoted = StructureSearchResult::from_rounds(
term0.clone(),
rho0.clone(),
vec![ledger_with(vec![MoveVerdict::Demoted { log_e: -1.0 }])],
);
assert!(
demoted.structure_changed(),
"a landed Demoted death folds an atom to ~0 routing ⇒ structure_changed() must be true"
);
}
#[test]
fn residual_bearing_fit_harvests_birth_proposal() {
let n = 40usize;
let active: Vec<Vec<bool>> = (0..n).map(|_| vec![true]).collect();
let (term, rho) = planted_term(&active);
let p = term.output_dim();
let mut residuals = Array2::<f64>::zeros((n, p));
let u = [0.6_f64, -0.4, 0.5, -0.3];
for row in 0..n {
let amp = 1.0 + (row as f64) / (n as f64);
for c in 0..p {
residuals[[row, c]] = amp * u[c % u.len()];
}
}
let params = HarvestParams {
max_fusions: 0,
max_fissions: 0,
max_births: 2,
};
let report = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
let births: usize = report
.proposals
.iter()
.filter(|p| matches!(p.mv, StructureMove::Birth { .. }))
.count();
assert!(
births >= 1,
"a residual-bearing fit with births enabled must harvest at least \
one birth proposal (so K can be discovered); got {:?}",
report.proposals.iter().map(|p| &p.mv).collect::<Vec<_>>()
);
assert!(
report.births_proposed >= 1,
"births_proposed must count the harvested births; got {}",
report.births_proposed
);
assert!(
report.birth_skipped_reason.is_none(),
"the birth channel must run (no skip) on a non-degenerate residual; got {:?}",
report.birth_skipped_reason
);
}
#[test]
fn fully_reconstructed_null_harvests_no_birth() {
let n = 40usize;
let active: Vec<Vec<bool>> = (0..n).map(|_| vec![true]).collect();
let (term, rho) = planted_term(&active);
let p = term.output_dim();
let zero_residual = Array2::<f64>::zeros((n, p));
let params = HarvestParams {
max_fusions: 0,
max_fissions: 0,
max_births: 2,
};
let report = harvest_move_proposals(&term, &rho, zero_residual.view(), ¶ms).unwrap();
let births: usize = report
.proposals
.iter()
.filter(|p| matches!(p.mv, StructureMove::Birth { .. }))
.count();
assert_eq!(
births, 0,
"a fully-reconstructed (zero-residual) null must harvest no birth \
proposal; got {births} births"
);
}
#[test]
fn planted_shatter_harvests_fusion_not_fission() {
let n = 30usize;
let active: Vec<Vec<bool>> = (0..n)
.map(|row| {
let dup = row % 3 == 0;
vec![dup, dup, row % 2 == 0]
})
.collect();
let (term, rho) = planted_term(&active);
let residuals = residuals_of(&term);
let params = HarvestParams {
max_fusions: 4,
max_fissions: 4,
max_births: 0,
};
let report = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
let has_fusion_01 = report.proposals.iter().any(|p| {
matches!(p.mv, StructureMove::Fusion { a, b } if (a, b) == (0, 1) || (a, b) == (1, 0))
});
assert!(
has_fusion_01,
"shattered duplicate pair (0,1) must yield a fusion proposal; got {:?}",
report.proposals.iter().map(|p| &p.mv).collect::<Vec<_>>()
);
let has_fission = report
.proposals
.iter()
.any(|p| matches!(p.mv, StructureMove::Fission { .. }));
assert!(
!has_fission,
"symmetric duplicate supports must not trigger an absorption fission audit"
);
}
#[test]
fn planted_absorption_harvests_fission_audit_with_loud_carve_skip() {
let n = 40usize;
let active: Vec<Vec<bool>> = (0..n)
.map(|row| {
let child = row % 4 == 0;
let parent = row % 2 == 0 || row % 4 == 1;
vec![parent, child, row % 5 == 0]
})
.collect();
let (term, rho) = planted_term(&active);
let residuals = residuals_of(&term);
let params = HarvestParams {
max_fusions: 4,
max_fissions: 4,
max_births: 0,
};
let report = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
let fissioned_parent = report
.proposals
.iter()
.any(|p| matches!(p.mv, StructureMove::Fission { atom: 0 }));
assert!(
fissioned_parent,
"nested-support parent (atom 0) must be flagged for a fission audit; got {:?}",
report.proposals.iter().map(|p| &p.mv).collect::<Vec<_>>()
);
assert_eq!(
report.fission_carve_ran_count, 0,
"1-D periodic atoms are not a product manifold; the within-atom carve cannot run"
);
assert!(
report.fission_carve_unavailable_count >= 1,
"the non-product fission candidate must be recorded as carve-unavailable, not silent"
);
assert!(
report.fission_carve_results.is_empty(),
"no carve ran, so there are no carve results to report"
);
}
#[test]
fn independent_atoms_harvest_no_fusion() {
let n = 60usize;
let active: Vec<Vec<bool>> = (0..n)
.map(|row| vec![row % 2 == 0, row % 3 == 0, row % 5 == 0])
.collect();
let (term, rho) = planted_term(&active);
let residuals = residuals_of(&term);
let params = HarvestParams {
max_fusions: 4,
max_fissions: 4,
max_births: 0,
};
let report = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
let has_fusion = report
.proposals
.iter()
.any(|p| matches!(p.mv, StructureMove::Fusion { .. }));
assert!(
!has_fusion,
"independent atom supports must not produce fusion proposals; got {:?}",
report.proposals.iter().map(|p| &p.mv).collect::<Vec<_>>()
);
}
fn tiled_circle_term(
n: usize,
k: usize,
decoder_scale: &[f64],
) -> (SaeManifoldTerm, SaeManifoldRho) {
assert_eq!(decoder_scale.len(), k, "one decoder scale per arc atom");
let p = 4usize;
let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(3).unwrap());
let coords = Array2::<f64>::from_shape_fn((n, 1), |(row, _)| row as f64 / n as f64);
let (phi, jet) = evaluator.evaluate(coords.view()).unwrap();
let m = phi.ncols();
let mut atoms = Vec::with_capacity(k);
let mut coord_blocks = Vec::with_capacity(k);
for (j, &scale) in decoder_scale.iter().enumerate() {
let mut decoder = Array2::<f64>::zeros((m, p));
decoder[[1, 0]] = scale;
decoder[[2, 1]] = scale;
let atom = SaeManifoldAtom::new(
format!("arc_{j}"),
SaeAtomBasisKind::Periodic,
1,
phi.clone(),
jet.clone(),
decoder,
Array2::<f64>::eye(m),
)
.unwrap()
.with_basis_second_jet(evaluator.clone());
atoms.push(atom);
coord_blocks.push(coords.clone());
}
let mut logits = Array2::<f64>::zeros((n, k));
for row in 0..n {
let owner = (row * k) / n; for j in 0..k {
logits[[row, j]] = if j == owner { ON } else { OFF };
}
}
let assignment = SaeAssignment::from_blocks_with_mode_and_manifolds(
logits,
coord_blocks,
vec![LatentManifold::Circle { period: 1.0 }; k],
AssignmentMode::softmax(1.0),
)
.unwrap();
let term = SaeManifoldTerm::new(atoms, assignment).unwrap();
let rho = SaeManifoldRho::new(0.0, 0.0, vec![Array1::<f64>::zeros(1); k]);
(term, rho)
}
#[test]
fn distinct_concentric_circles_do_not_glue() {
let n = 40usize;
let (term, rho) = tiled_circle_term(n, 2, &[1.0, 2.0]); let residuals = Array2::<f64>::zeros((n, 4));
let params = HarvestParams {
max_fusions: 4,
max_fissions: 4,
max_births: 0,
};
let report = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
assert!(
!report
.proposals
.iter()
.any(|p| matches!(p.mv, StructureMove::Fusion { .. })),
"disjoint supports must yield no co-activation fusion"
);
assert!(
report.glue_candidates_screened >= 1,
"the shared-plane pair must be geometrically screened"
);
let (e_distinct, _) = unit_speed_glue_certificate(&term, residuals.view(), 0, 1)
.expect("the aligned pair yields a seam e-value and transition certificate");
assert!(
e_distinct.log_e_value < 0.0,
"distinct concentric circles must NOT glue (negative log-e), got {}",
e_distinct.log_e_value
);
for p in &report.proposals {
if matches!(p.mv, StructureMove::Glue { .. }) {
assert!(
p.trigger < 0.0,
"the distinct-circle glue must carry negative evidence, got {}",
p.trigger
);
}
}
}
#[test]
fn closed_form_transition_uses_first_nonzero_harmonic_without_scanning() {
let mut decoder_a = Array2::<f64>::zeros((5, 3));
decoder_a[[3, 0]] = 2.0;
decoder_a[[4, 1]] = 1.0;
decoder_a[[0, 2]] = 0.25;
let decoder_b = periodic_decoder_under_transition(decoder_a.view(), -1, 0.125).unwrap();
let (sign, offset) =
fit_periodic_transition_from_decoders(decoder_a.view(), decoder_b.view()).unwrap();
assert_eq!(sign, -1);
let recovered = periodic_decoder_under_transition(decoder_a.view(), sign, offset).unwrap();
for (actual, expected) in recovered.iter().zip(decoder_b.iter()) {
assert!((actual - expected).abs() < 32.0 * f64::EPSILON);
}
}
#[test]
fn orientation_reversing_seam_registers_atlas_without_destructive_fusion() {
let n = 40usize;
let (mut term, rho) = tiled_circle_term(n, 2, &[1.0, 1.0]);
term.atoms[1].decoder_coefficients[[1, 0]] = -1.0;
let residuals = Array2::<f64>::zeros((n, 4));
let (transition, _) = unit_speed_glue_certificate(&term, residuals.view(), 0, 1)
.expect("reflected charts have an exact certified seam");
assert_eq!(transition.sign, -1);
assert!(transition.log_e_value > 5.0);
let report = harvest_move_proposals(
&term,
&rho,
residuals.view(),
&HarvestParams {
max_fusions: 4,
max_fissions: 0,
max_births: 0,
},
)
.unwrap();
let mv = report
.proposals
.iter()
.find_map(|proposal| match proposal.mv {
StructureMove::Glue {
a,
b,
outcome: ChartGlueOutcome::RegisterAtlas,
} => Some(StructureMove::Glue {
a,
b,
outcome: ChartGlueOutcome::RegisterAtlas,
}),
_ => None,
})
.expect("negative seam must propose atlas registration");
let fitted_before = term.try_fitted().unwrap();
let (registered, _) = apply_structure_move(&term, &rho, &mv, &[]).unwrap();
assert_eq!(registered.k_atoms(), 2, "both local charts must survive");
assert_eq!(registered.semantic_atom_count(), 1);
assert_eq!(registered.chart_atlases().len(), 1);
assert_eq!(registered.chart_atlases()[0].transitions()[0].sign, -1);
assert_eq!(
registered.try_fitted().unwrap(),
fitted_before,
"atlas registration is an image-exact quotient"
);
let assignments = registered.assignment.assignments();
for row in 0..n {
let (activation, partition) = registered
.atlas_partition_of_unity(0, assignments.row(row))
.unwrap();
assert!((partition.sum() - 1.0).abs() < 8.0 * f64::EPSILON);
for (slot, &chart) in registered.chart_atlases()[0].charts().iter().enumerate() {
assert!(
(activation * partition[slot] - assignments[[row, chart]]).abs()
< 8.0 * f64::EPSILON
);
}
}
}
#[test]
fn over_tiling_physical_excision_reduces_k_toward_one() {
use gam_solve::structure_search::{MoveRecord, MoveVerdict};
let n = 32usize;
let (mut term, mut rho) = tiled_circle_term(n, 4, &[1.0; 4]);
assert_eq!(term.k_atoms(), 4);
let residuals0 = Array2::<f64>::zeros((n, 4));
let (e_arc, _) = unit_speed_glue_certificate(&term, residuals0.view(), 0, 2)
.expect("a d=1 aligned disjoint pair yields a certified seam e-value");
assert!(
e_arc.log_e_value > 5.0,
"e_glue must certify two arcs of one circle, got {}",
e_arc.log_e_value
);
let params0 = HarvestParams {
max_fusions: 16,
max_fissions: 0,
max_births: 0,
};
let report0 = harvest_move_proposals(&term, &rho, residuals0.view(), ¶ms0).unwrap();
assert!(
!report0
.proposals
.iter()
.any(|p| matches!(p.mv, StructureMove::Fusion { .. })),
"disjoint tiling must yield no co-activation fusion"
);
assert!(
report0.glues_proposed >= 3,
"a spanning set (≥ k−1 = 3 edges) must reassemble the 4 arcs, got {}",
report0.glues_proposed
);
let first_epoch_glue_claims: Vec<ClaimKind> = report0
.proposals
.iter()
.filter(|proposal| matches!(proposal.mv, StructureMove::Glue { .. }))
.map(|proposal| proposal.claim.clone())
.collect();
assert!(!first_epoch_glue_claims.is_empty());
term.assignment.frozen_logits = Some(term.assignment.logits.clone());
term.last_frames_active = true;
term.fixed_decoder_assembly = true;
term.border_hbb_workspace = Array2::<f64>::ones((3, 3));
term.decoder_repulsion_gate = Some(vec![(0, 1, 1.0)]);
term.streaming_gates_frozen = true;
term.expected_evidence_gauge_deflated_directions = Some(7);
term.evidence_gauge_deflation_reanchors = 2;
term.evidence_gauge_deflation_last_delta_sign = -1;
term.dictionary_cocollapse_reseeds = 3;
term.structural_cocollapse_reseeds = 4;
term.softmax_active_cap = Some(3);
let accepted_glue = |a: usize, b: usize| MoveRecord {
mv: StructureMove::Glue {
a,
b,
outcome: ChartGlueOutcome::Fuse,
},
trigger: 40.0,
structure_hash: 0,
claim: ClaimKind::Custom {
label: format!("seam_glue:{a}:{b}"),
},
verdict: MoveVerdict::Accepted { log_e: 40.0 },
};
let ledger = SearchLedger {
alpha: 0.05,
moves: vec![accepted_glue(0, 1), accepted_glue(2, 3)],
collapse_events: Vec::new(),
};
let removed =
compact_glued_atoms(&mut term, &mut rho, &ledger, &report0.certified_glues).unwrap();
assert_eq!(
removed, 2,
"both folded partners must be physically excised"
);
assert_eq!(
term.k_atoms(),
2,
"physical excision must reduce K from 4 to 2 (no active-mass resurrection)"
);
assert!(
term.assignment
.logits
.rows()
.into_iter()
.all(|row| row.as_slice().is_some()),
"compaction must materialize a row-contiguous router for the polish refit"
);
assert_eq!(
rho.log_ard.len(),
2,
"ρ ARD blocks must fall in lock-step with the atoms"
);
assert_eq!(
rho.log_lambda_smooth.len(),
2,
"ρ smoothness blocks must fall in lock-step with the atoms"
);
assert!(
term.assignment.frozen_logits.is_none(),
"an old-K frozen router cannot survive compaction"
);
assert!(!term.last_frames_active);
assert!(!term.fixed_decoder_assembly);
assert_eq!(term.border_hbb_workspace.dim(), (0, 0));
assert!(term.decoder_repulsion_gate.is_none());
assert!(!term.streaming_gates_frozen);
assert_eq!(term.expected_evidence_gauge_deflated_directions, None);
assert_eq!(term.evidence_gauge_deflation_reanchors, 0);
assert_eq!(term.evidence_gauge_deflation_last_delta_sign, 0);
assert_eq!(term.dictionary_cocollapse_reseeds, 0);
assert_eq!(term.structural_cocollapse_reseeds, 0);
assert_eq!(
term.softmax_active_cap, None,
"a width-3 cap is inert after compaction to K=2"
);
let residuals = Array2::<f64>::zeros((n, 4));
let params = HarvestParams {
max_fusions: 4,
max_fissions: 0,
max_births: 0,
};
let report2 = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
assert!(
report2.glues_proposed >= 1,
"the two reassembled half-circle survivors must still glue toward K=1"
);
for proposal in report2
.proposals
.iter()
.filter(|proposal| matches!(proposal.mv, StructureMove::Glue { .. }))
{
assert!(
!first_epoch_glue_claims.contains(&proposal.claim),
"a reduced dictionary must not reuse old atom-index evidence: {:?}",
proposal.claim
);
}
}
#[test]
fn physical_excision_transplants_coords_from_the_live_seam() {
use gam_solve::structure_search::{MoveRecord, MoveVerdict};
let (mut term, mut rho) = tiled_circle_term(32, 2, &[1.0; 2]);
assert!(term.assignment.frozen_logits.is_none());
let seam = fit_seam_transition(&term, 0, 1).expect("live pair has a seam");
let residuals = Array2::<f64>::zeros((32, 4));
let (_, certificate) = unit_speed_glue_certificate(&term, residuals.view(), 0, 1)
.expect("live pair has a harvest-time glue certificate");
let da = term.assignment.coords[0].latent_dim();
let db = term.assignment.coords[1].latent_dim();
let flat_b = term.assignment.coords[1].as_flat().to_owned();
let expected: Vec<(usize, f64)> = seam
.rows_b
.iter()
.map(|&row| {
let mapped = (seam.sign * flat_b[row * db] + seam.offset).rem_euclid(seam.period);
(row, mapped)
})
.collect();
let mut flat_a = term.assignment.coords[0].as_flat().to_owned();
for &(row, mapped) in &expected {
flat_a[row * da] = (mapped + 0.37).rem_euclid(seam.period);
}
term.assignment.coords[0].set_flat(flat_a.view());
term.atoms[1].decoder_coefficients.fill(0.0);
assert!(
fit_seam_transition(&term, 0, 1).is_none(),
"post-harvest fixture must make seam re-fitting impossible"
);
let ledger = SearchLedger {
alpha: 0.05,
moves: vec![MoveRecord {
mv: StructureMove::Glue {
a: 0,
b: 1,
outcome: ChartGlueOutcome::Fuse,
},
trigger: 40.0,
structure_hash: 0,
claim: ClaimKind::Custom {
label: "test-live-seam".to_string(),
},
verdict: MoveVerdict::Accepted { log_e: 40.0 },
}],
collapse_events: Vec::new(),
};
compact_glued_atoms(&mut term, &mut rho, &ledger, &[certificate]).unwrap();
assert_eq!(term.k_atoms(), 1);
let survivor = term.assignment.coords[0].as_flat();
for (row, mapped) in expected {
assert!(
(survivor[row * da] - mapped).abs() < 1.0e-12,
"row {row}: survivor coordinate {} != seam-mapped {mapped}",
survivor[row * da]
);
}
}
#[test]
fn physical_excision_validation_is_transactional() {
use gam_solve::structure_search::{MoveRecord, MoveVerdict};
let (mut term, mut rho) = tiled_circle_term(16, 3, &[1.0; 3]);
let atoms_before = term.k_atoms();
let logits_before = term.assignment.logits.clone();
rho.log_ard.pop();
let remove = std::collections::BTreeSet::from([1usize]);
let err = remove_atoms(&mut term, &mut rho, &remove).unwrap_err();
assert!(
err.contains("rho per-atom lengths"),
"unexpected error: {err}"
);
assert_eq!(term.k_atoms(), atoms_before);
assert_eq!(term.assignment.logits, logits_before);
rho.log_ard.push(Array1::zeros(1));
let accepted = |a: usize, b: usize| MoveRecord {
mv: StructureMove::Glue {
a,
b,
outcome: ChartGlueOutcome::Fuse,
},
trigger: 40.0,
structure_hash: 0,
claim: ClaimKind::Custom {
label: format!("test-glue:{a}:{b}"),
},
verdict: MoveVerdict::Accepted { log_e: 40.0 },
};
let overlapping = SearchLedger {
alpha: 0.05,
moves: vec![accepted(0, 1), accepted(1, 2)],
collapse_events: Vec::new(),
};
let err = compact_glued_atoms(&mut term, &mut rho, &overlapping, &[]).unwrap_err();
assert!(
err.contains("not an atom-disjoint matching"),
"unexpected error: {err}"
);
assert_eq!(term.k_atoms(), atoms_before);
assert_eq!(term.assignment.logits, logits_before);
let missing_certificate = SearchLedger {
alpha: 0.05,
moves: vec![accepted(0, 1)],
collapse_events: Vec::new(),
};
let err = compact_glued_atoms(&mut term, &mut rho, &missing_certificate, &[]).unwrap_err();
assert!(
err.contains("no harvest-time certificate"),
"unexpected error: {err}"
);
assert_eq!(term.k_atoms(), atoms_before);
assert_eq!(term.assignment.logits, logits_before);
}
#[test]
fn diverged_ard_and_terminal_collapse_harvest_deaths() {
let n = 20usize;
let active: Vec<Vec<bool>> = (0..n).map(|row| vec![true, row % 2 == 0, false]).collect();
let (mut term, mut rho) = planted_term(&active);
rho.log_ard[2] = Array1::from_elem(1, ARD_DIVERGENCE_LOG_PRECISION + 5.0);
term.record_collapse_event(CollapseEvent {
iteration: 3,
atom: 1,
max_active_mass: 1e-6,
floor: 1e-3,
action: CollapseAction::Terminal,
});
let residuals = residuals_of(&term);
let params = HarvestParams {
max_fusions: 0,
max_fissions: 0,
max_births: 0,
};
let report = harvest_move_proposals(&term, &rho, residuals.view(), ¶ms).unwrap();
let death_atoms: Vec<usize> = report
.proposals
.iter()
.filter_map(|p| match p.mv {
StructureMove::Death { atom } => Some(atom),
_ => None,
})
.collect();
assert!(
death_atoms.contains(&2),
"diverged ARD on atom 2 must yield a death proposal; got {death_atoms:?}"
);
assert!(
death_atoms.contains(&1),
"terminal collapse on atom 1 must yield a death proposal; got {death_atoms:?}"
);
}
#[test]
fn apply_move_restructures_warm() {
let n = 12usize;
let active: Vec<Vec<bool>> = (0..n).map(|row| vec![true, row % 2 == 0]).collect();
let (term, rho) = planted_term(&active);
let k0 = term.k_atoms();
let (fissioned, fissioned_rho) =
apply_structure_move(&term, &rho, &StructureMove::Fission { atom: 0 }, &[]).unwrap();
assert_eq!(fissioned.k_atoms(), k0 + 1);
assert_eq!(fissioned_rho.log_ard.len(), k0 + 1);
assert_eq!(
fissioned_rho.log_lambda_smooth.len(),
fissioned.k_atoms(),
"fission must grow per-atom log_lambda_smooth in lockstep with K"
);
let (fused, _) =
apply_structure_move(&term, &rho, &StructureMove::Fusion { a: 0, b: 1 }, &[]).unwrap();
assert_eq!(fused.k_atoms(), k0);
let fused_assign = fused.assignment.assignments();
assert!(
fused_assign.column(1).iter().all(|&m| m < 1e-6),
"fused-away atom 1 must route to ~0 mass"
);
let (dead, _) =
apply_structure_move(&term, &rho, &StructureMove::Death { atom: 1 }, &[]).unwrap();
assert_eq!(dead.k_atoms(), k0);
let dead_assign = dead.assignment.assignments();
assert!(dead_assign.column(1).iter().all(|&m| m < 1e-6));
let p = term.output_dim();
let m = term.atoms[0].basis_size();
let mut decoder = Array2::<f64>::zeros((m, p));
decoder[[0, 0]] = 0.7;
let birth_target = term.atoms[0].basis_values.dot(&decoder); let (born, born_rho) = apply_structure_move(
&term,
&rho,
&StructureMove::Birth { candidate: 0 },
&[decoder],
)
.unwrap();
assert_eq!(born.k_atoms(), k0 + 1);
assert_eq!(born_rho.log_ard.len(), k0 + 1);
assert_eq!(born_rho.log_lambda_smooth.len(), k0 + 1);
let born_atom = &born.atoms[k0];
let born_image = born_atom.basis_values.dot(&born_atom.decoder_coefficients);
assert_eq!(born_image.dim(), birth_target.dim());
let mut max_recon_err = 0.0_f64;
for (a, b) in born_image.iter().zip(birth_target.iter()) {
max_recon_err = max_recon_err.max((a - b).abs());
}
assert!(
max_recon_err < 1e-3,
"born atom must reconstruct the residual-factor image (penalized fit); \
max |Φ_born·B_born − Φ_template·factor_dir| = {max_recon_err:.3e} (> 1e-3)"
);
}
#[test]
fn grown_atom_count_assembles_without_lambda_smooth_oob_357() {
let n = 16usize;
let active: Vec<Vec<bool>> = (0..n).map(|row| vec![true, row % 2 == 0]).collect();
let (term, rho) = planted_term(&active);
let target = Array2::<f64>::from_shape_fn((n, term.output_dim()), |(row, col)| {
0.1 * (row as f64) - 0.05 * (col as f64)
});
let (fissioned, fissioned_rho) =
apply_structure_move(&term, &rho, &StructureMove::Fission { atom: 0 }, &[]).unwrap();
assert_eq!(fissioned_rho.log_lambda_smooth.len(), fissioned.k_atoms());
let mut fissioned = fissioned;
fissioned
.assemble_arrow_schur_scaled(target.view(), &fissioned_rho, None, 1.0)
.expect("post-fission assembly must not panic or error on the grown atom set");
let p = term.output_dim();
let m = term.atoms[0].basis_size();
let mut decoder = Array2::<f64>::zeros((m, p));
decoder[[0, 0]] = 0.5;
let (born, born_rho) = apply_structure_move(
&term,
&rho,
&StructureMove::Birth { candidate: 0 },
&[decoder],
)
.unwrap();
assert_eq!(born_rho.log_lambda_smooth.len(), born.k_atoms());
let mut born = born;
born.assemble_arrow_schur_scaled(target.view(), &born_rho, None, 1.0)
.expect("post-birth assembly must not panic or error on the grown atom set");
}
#[test]
fn round_driver_ledger_is_byte_deterministic() {
let n = 24usize;
let active: Vec<Vec<bool>> = (0..n)
.map(|row| {
let dup = row % 3 == 0;
vec![dup, dup, row % 2 == 0]
})
.collect();
let run = || {
let (term, rho) = planted_term(&active);
let target = Array2::<f64>::zeros((n, term.output_dim()));
let mut ledger = gam_terms::inference::structure_evidence::StructureLedger::new();
let budget = MoveBudget {
max_moves: 4,
alpha: 0.05,
};
let params = HarvestParams {
max_fusions: 4,
max_fissions: 0,
max_births: 0,
};
let config = RoundDriverConfig {
n_shards: 3,
budget,
max_rounds: 2,
harvest_params: params,
curl: None,
};
run_structure_search_rounds(
term,
rho,
target.view(),
config,
&mut ledger,
|t, r, _| (t, r),
|t, r, _| (t, r),
)
.unwrap()
};
let a = run();
let b = run();
let sa = serde_json::to_string(&a.rounds).unwrap();
let sb = serde_json::to_string(&b.rounds).unwrap();
assert_eq!(
sa, sb,
"identical inputs must produce a byte-identical ledger"
);
assert_eq!(a.term.k_atoms(), b.term.k_atoms());
}
#[test]
fn scoring_iter_cap_preserves_moves_and_adopted_fit() {
let n = 40usize;
let active: Vec<Vec<bool>> = (0..n).map(|_| vec![true]).collect();
let p = 4usize;
let u = [0.6_f64, -0.4, 0.5, -0.3];
let mut target = Array2::<f64>::zeros((n, p));
for row in 0..n {
let amp = 1.0 + (row as f64) / (n as f64);
for c in 0..p {
target[[row, c]] = amp * u[c % u.len()];
}
}
let config = RoundDriverConfig {
n_shards: 4,
budget: MoveBudget {
max_moves: 4,
alpha: 0.05,
},
max_rounds: 2,
harvest_params: HarvestParams {
max_fusions: 2,
max_fissions: 2,
max_births: 2,
},
curl: None,
};
let full_iters = 24usize;
let run = |scoring_inner_max_iter: usize| {
let (term, rho) = planted_term(&active);
let mut ledger = StructureLedger::new();
let refit_params = ProductionRefitParams {
inner_max_iter: full_iters,
scoring_inner_max_iter,
learning_rate: 1.0,
ridge_ext_coord: 1e-6,
ridge_beta: 1e-6,
};
let result = run_production_structure_search(
term,
rho,
target.view(),
config,
refit_params,
&mut ledger,
)
.unwrap();
let fitted = result.term.try_fitted().unwrap();
(result, fitted)
};
let (reference, ref_fitted) = run(full_iters);
let (capped, cap_fitted) = run(4);
use gam_solve::structure_search::MoveVerdict;
let verdict_kind = |v: &MoveVerdict| -> &'static str {
match v {
MoveVerdict::Accepted { .. } => "Accepted",
MoveVerdict::Contested { .. } => "Contested",
MoveVerdict::Demoted { .. } => "Demoted",
MoveVerdict::Vetoed { .. } => "Vetoed",
MoveVerdict::Deduplicated => "Deduplicated",
MoveVerdict::Stale => "Stale",
MoveVerdict::Deferred => "Deferred",
}
};
let round_moves = |rounds: &[SearchLedger]| -> String {
serde_json::to_string(
&rounds
.iter()
.map(|r| {
r.moves
.iter()
.map(|m| {
(
serde_json::to_string(&m.mv).unwrap(),
m.structure_hash,
serde_json::to_string(&m.claim).unwrap(),
verdict_kind(&m.verdict),
)
})
.collect::<Vec<_>>()
})
.collect::<Vec<_>>(),
)
.unwrap()
};
assert_eq!(
round_moves(&reference.rounds),
round_moves(&capped.rounds),
"scoring-iteration cap changed the accepted-move trajectory — the e-gate \
decisions are NOT cap-invariant (the #1026 economy is unsound)"
);
assert_eq!(
reference.term.k_atoms(),
capped.term.k_atoms(),
"scoring cap changed the discovered dictionary size"
);
assert_eq!(ref_fitted.dim(), cap_fitted.dim());
let mut max_abs = 0.0_f64;
for (a, b) in ref_fitted.iter().zip(cap_fitted.iter()) {
max_abs = max_abs.max((a - b).abs());
}
assert!(
max_abs < 1e-6,
"capped-scoring adopted fit diverged from the full-iter reference by \
{max_abs:.3e} (> 1e-6); the polish did not reach the same optimum"
);
}
#[test]
fn estimation_eval_split_is_disjoint() {
let target = Array2::<f64>::zeros((20, 3));
let split = estimation_eval_split(target.view(), 4);
assert!(!split.estimation_rows.is_empty());
assert!(!split.shards.is_empty());
let est: std::collections::HashSet<usize> = split.estimation_rows.iter().copied().collect();
for shard in &split.shards {
for &row in &shard.rows {
assert!(
!est.contains(&row),
"eval shard row {row} must not be in the estimation set"
);
}
}
}
#[test]
fn birth_topology_race_assigns_circle_vs_line_by_evidence() {
use std::f64::consts::TAU;
let n = 80usize;
let coords = Array2::<f64>::from_shape_fn((n, 1), |(row, _)| row as f64 / n as f64);
let p = 4usize;
let mut circle_target = Array2::<f64>::zeros((n, p));
for row in 0..n {
let t = coords[[row, 0]];
circle_target[[row, 0]] = (TAU * t).cos();
circle_target[[row, 1]] = (TAU * t).sin();
}
let mut line_target = Array2::<f64>::zeros((n, p));
let u = [0.7_f64, -0.4, 0.5, -0.2];
for row in 0..n {
let t = coords[[row, 0]];
for c in 0..p {
line_target[[row, c]] = t * u[c];
}
}
let weights = Array1::<f64>::ones(n);
let circle_fit =
race_birth_topology(coords.view(), circle_target.view(), weights.view(), 1)
.expect("circle race runs")
.expect("circle race has a realizable candidate");
let line_fit = race_birth_topology(coords.view(), line_target.view(), weights.view(), 1)
.expect("line race runs")
.expect("line race has a realizable candidate");
assert_eq!(
circle_fit.basis_kind,
SaeAtomBasisKind::Periodic,
"a circular birth residual must win the circle (Periodic) topology"
);
assert_eq!(
line_fit.basis_kind,
SaeAtomBasisKind::EuclideanPatch,
"a straight birth residual must win the line (EuclideanPatch) topology"
);
assert_ne!(
circle_fit.basis_kind, line_fit.basis_kind,
"the discovery must assign DIFFERENT topologies to the circle and line \
atoms (evidence-chosen, not inherited)"
);
}
#[test]
fn birth_topology_race_d2_includes_and_selects_cylinder() {
use std::f64::consts::TAU;
let n = 120usize;
let coords = Array2::<f64>::from_shape_fn((n, 2), |(row, axis)| {
if axis == 0 {
(row as f64 / n as f64) * 2.0
} else {
(row as f64 / n as f64) * 3.0 - 1.5
}
});
let specs = topology_candidates_for_dim(coords.view(), 2).expect("d=2 candidates build");
let has_cylinder = specs
.iter()
.any(|s| s.basis_kind == SaeAtomBasisKind::Cylinder);
assert!(
has_cylinder,
"the d=2 topology-race candidate set MUST include the Cylinder kind; got {:?}",
specs.iter().map(|s| &s.basis_kind).collect::<Vec<_>>()
);
let has_torus = specs
.iter()
.any(|s| s.basis_kind == SaeAtomBasisKind::Torus);
let has_sphere = specs
.iter()
.any(|s| s.basis_kind == SaeAtomBasisKind::Sphere);
let has_patch = specs
.iter()
.any(|s| s.basis_kind == SaeAtomBasisKind::EuclideanPatch);
assert!(
has_torus && has_sphere && has_patch,
"the d=2 race must be COMPLETE (torus + sphere + euclidean + cylinder)"
);
let p = 4usize;
let mut cyl_target = Array2::<f64>::zeros((n, p));
for row in 0..n {
let phase = coords[[row, 0]];
let mag = coords[[row, 1]];
cyl_target[[row, 0]] = (TAU * phase).cos();
cyl_target[[row, 1]] = (TAU * phase).sin();
cyl_target[[row, 2]] = mag;
}
let weights = Array1::<f64>::ones(n);
let cyl_fit = race_birth_topology(coords.view(), cyl_target.view(), weights.view(), 2)
.expect("cylinder race runs")
.expect("cylinder race has a realizable candidate");
assert_eq!(
cyl_fit.basis_kind,
SaeAtomBasisKind::Cylinder,
"a cylindrical birth residual (periodic along one axis, linear along the \
other) must win the Cylinder topology by evidence; got {:?}",
cyl_fit.basis_kind
);
}
#[test]
fn born_atom_reports_finite_uncertainty_band() {
let n = 48usize;
let active: Vec<Vec<bool>> = (0..n).map(|_| vec![true]).collect();
let (term, rho) = planted_term(&active);
let k_seed = term.k_atoms();
let p = term.output_dim();
let m = term.atoms[0].basis_size();
let mut decoder = Array2::<f64>::zeros((m, p));
decoder[[1, 0]] = 0.9;
decoder[[2, 1]] = -0.6;
let (mut born, _born_rho) = apply_structure_move(
&term,
&rho,
&StructureMove::Birth { candidate: 0 },
&[decoder],
)
.expect("birth applies");
assert_eq!(born.k_atoms(), k_seed + 1, "the birth grows K by one");
let target = born.try_fitted().expect("born term reconstructs");
let dispersion = 1.0e-2_f64;
born.set_atom_inner_fits(target.view(), dispersion)
.expect("inner fits build");
let mut unc = born.shape_uncertainty_without_decoder_covariance(dispersion);
unc.atoms.truncate(k_seed);
assert_eq!(
unc.atoms.len(),
k_seed,
"seed-K Schur band omits the born atom"
);
born.complete_born_atom_shape_bands(&mut unc)
.expect("born-atom band completes");
assert_eq!(
unc.atoms.len(),
born.k_atoms(),
"completion must grow the band list to the post-search atom count"
);
let born_band = &unc.atoms[k_seed];
assert!(
born_band.band_sd.nrows() > 0 && born_band.band_sd.ncols() == p,
"the born atom's band must be shaped (G>0, p)"
);
let mut any_positive = false;
for &sd in born_band.band_sd.iter() {
assert!(
sd.is_finite() && sd >= 0.0,
"born-atom band sd must be finite and non-negative; got {sd}"
);
if sd > 0.0 {
any_positive = true;
}
}
assert!(
any_positive,
"a born atom with a non-degenerate inner Hessian must report a strictly \
positive uncertainty somewhere (a finite band, never all-zero / missing)"
);
}
#[test]
fn production_gate_consumes_corrected_pg_normalizer() {
let n = 32usize;
let null_active: Vec<Vec<bool>> = (0..n).map(|_| vec![true, true]).collect();
let cand_active: Vec<Vec<bool>> = (0..n).map(|_| vec![true, true, true]).collect();
let (null_term, _) = planted_term(&null_active);
let (cand_term, _) = planted_term(&cand_active);
assert_eq!(null_term.k_atoms(), 2);
assert_eq!(cand_term.k_atoms(), 3, "candidate grows K by one atom");
let p = null_term.output_dim();
let target = Arc::new(Array2::<f64>::zeros((n, p)));
let shard = RowBlockShard {
target: target.clone(),
rows: (0..n).collect(),
};
let null_gate = gate_block_log_evidence(&null_term, &shard);
let cand_gate = gate_block_log_evidence(&cand_term, &shard);
assert!(
null_gate.is_finite() && cand_gate.is_finite(),
"gate-block evidence must be finite on a well-posed gate block"
);
let log_2pi = (2.0 * std::f64::consts::PI).ln();
let gate_delta = cand_gate - null_gate;
let per_atom_no_norm = |term: &SaeManifoldTerm| -> f64 {
let dg = term.k_atoms() as f64; gate_block_log_evidence(term, &shard) + 0.5 * dg * log_2pi
};
let no_norm_delta = per_atom_no_norm(&cand_term) - per_atom_no_norm(&null_term);
let normalizer_in_delta = gate_delta - no_norm_delta;
assert!(
(normalizer_in_delta + 0.5 * log_2pi).abs() < 1e-9,
"the gate-block normalizer in the K→K+1 difference must be the \
corrected −½·log(2π) Occam penalty, got {normalizer_in_delta} \
(buggy +½·log(2π) = {})",
0.5 * log_2pi
);
let full = eval_log_lik(&cand_term, &shard);
let recon_only = {
let fitted = cand_term.try_fitted().unwrap();
let mut sse = 0.0;
for &row in &shard.rows {
for out in 0..p {
let d = fitted[[row, out]] - shard.target[[row, out]];
sse += d * d;
}
}
-0.5 * sse
};
assert!(
(full - (recon_only + cand_gate)).abs() < 1e-9,
"the live per-shard likelihood must equal reconstruction + the \
PG gate-block evidence (so the corrected normalizer reaches the gate)"
);
}
#[test]
fn fission_breaks_symmetry_so_children_can_separate() {
let (term, rho) = planted_term(&vec![vec![true]; 8]);
assert_eq!(term.k_atoms(), 1);
let orig = term.atoms[0].decoder_coefficients.clone();
let (child, _child_rho) =
apply_structure_move(&term, &rho, &StructureMove::Fission { atom: 0 }, &[]).unwrap();
assert_eq!(child.k_atoms(), 2, "fission must add one atom");
let d0 = &child.atoms[0].decoder_coefficients;
let d1 = &child.atoms[1].decoder_coefficients;
let sep = (d0 - d1).iter().map(|x| x * x).sum::<f64>().sqrt();
let scale = orig.iter().map(|x| x * x).sum::<f64>().sqrt().max(1e-12);
assert!(
sep / scale > 1.0e-3,
"fission children must NOT be identical (symmetric saddle); rel sep = {}",
sep / scale
);
let combined = (d0 + d1).mapv(|x| 0.5 * x);
let warm_err = (&combined - &orig)
.iter()
.map(|x| x * x)
.sum::<f64>()
.sqrt();
assert!(
warm_err < 1.0e-12,
"mass-split combined decoder must equal the original; err = {warm_err}"
);
for row in 0..child.assignment.logits.nrows() {
assert!(
(child.assignment.logits[[row, 0]] - child.assignment.logits[[row, 1]]).abs()
< 1e-12,
"fission must split routing mass 50/50 (equal child logits)"
);
}
}
#[test]
fn fusion_preserves_combined_softmax_mass() {
let (term, rho) = planted_term(&vec![vec![true, true, true]; 6]);
let combined: Vec<f64> = (0..6)
.map(|r| {
let a = term.assignment.try_assignments_row(r).unwrap();
a[0] + a[1]
})
.collect();
let (fused, _) =
apply_structure_move(&term, &rho, &StructureMove::Fusion { a: 0, b: 1 }, &[]).unwrap();
for r in 0..6 {
let a = fused.assignment.try_assignments_row(r).unwrap();
assert!(
(a[0] - combined[r]).abs() < 1e-6,
"fused atom must carry the COMBINED softmax mass (logsumexp, not \
max): got {}, want {} (row {r})",
a[0],
combined[r]
);
assert!(
combined[r] > 0.6,
"fixture must exercise a co-active pair (combined mass {} should be ~⅔)",
combined[r]
);
}
}
#[test]
fn fusion_of_zero_mass_pair_yields_neg_inf_not_nan() {
let (mut term, rho) = planted_term(&vec![vec![true, true, true]; 6]);
assert!(
matches!(term.assignment.mode, AssignmentMode::Softmax { .. }),
"fixture must be softmax-routed to exercise the logsumexp combine"
);
term.assignment.logits[[0, 0]] = f64::NEG_INFINITY;
term.assignment.logits[[0, 1]] = f64::NEG_INFINITY;
let (fused, _) =
apply_structure_move(&term, &rho, &StructureMove::Fusion { a: 0, b: 1 }, &[]).unwrap();
let folded = fused.assignment.logits[[0, 0]];
assert!(
!folded.is_nan(),
"fused zero-mass logit must not be NaN (got {folded})"
);
assert_eq!(
folded,
f64::NEG_INFINITY,
"combined mass of two zero-mass atoms is zero → logit -∞"
);
for c in 0..fused.assignment.logits.ncols() {
assert!(
!fused.assignment.logits[[0, c]].is_nan(),
"row 0 col {c} must not be NaN after the fold"
);
}
}
fn linear_line_atom(name: &str, coord: &Array1<f64>, dir: &Array1<f64>) -> SaeManifoldAtom {
let n = coord.len();
let p = dir.len();
let mut phi = Array2::<f64>::zeros((n, 2));
let mut jet = ndarray::Array3::<f64>::zeros((n, 2, 1));
for r in 0..n {
phi[[r, 0]] = 1.0;
phi[[r, 1]] = coord[r];
jet[[r, 0, 0]] = 0.0;
jet[[r, 1, 0]] = 1.0;
}
let mut decoder = Array2::<f64>::zeros((2, p));
for j in 0..p {
decoder[[1, j]] = dir[j];
}
SaeManifoldAtom::new(
name.to_string(),
SaeAtomBasisKind::Linear,
1,
phi,
jet,
decoder,
Array2::<f64>::eye(2),
)
.unwrap()
}
fn shattered_plane_term(gaussian: bool) -> (SaeManifoldTerm, SaeManifoldRho) {
let n = 600usize;
let radius = 3.0_f64;
let u = Array1::from_vec(vec![1.0, 0.0, 0.0, 0.0]);
let v = Array1::from_vec(vec![0.0, 1.0, 0.0, 0.0]);
let neg_u = u.mapv(|x| -x);
let neg_v = v.mapv(|x| -x);
let mut s = 0xC0FFEE_u64;
let lcg = |st: &mut u64| -> f64 {
*st = st
.wrapping_mul(6364136223846793005)
.wrapping_add(1442695040888963407);
((*st >> 11) as f64) / ((1u64 << 53) as f64)
};
let mut xs = Array1::<f64>::zeros(n);
let mut ys = Array1::<f64>::zeros(n);
for r in 0..n {
if gaussian {
let u1 = lcg(&mut s).max(1e-12);
let u2 = lcg(&mut s);
let g0 = (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos();
let g1 = (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).sin();
xs[r] = radius * g0;
ys[r] = radius * g1;
} else {
let th = std::f64::consts::TAU * (r as f64 + 0.5) / n as f64;
xs[r] = radius * th.cos();
ys[r] = radius * th.sin();
}
}
let cu: Array1<f64> = xs.mapv(|x| x.max(0.0));
let cnu: Array1<f64> = xs.mapv(|x| (-x).max(0.0));
let cv: Array1<f64> = ys.mapv(|y| y.max(0.0));
let cnv: Array1<f64> = ys.mapv(|y| (-y).max(0.0));
let atoms = vec![
linear_line_atom("half_+u", &cu, &u),
linear_line_atom("half_-u", &cnu, &neg_u),
linear_line_atom("half_+v", &cv, &v),
linear_line_atom("half_-v", &cnv, &neg_v),
];
let coord_blocks = vec![
cu.clone().insert_axis(ndarray::Axis(1)),
cnu.clone().insert_axis(ndarray::Axis(1)),
cv.clone().insert_axis(ndarray::Axis(1)),
cnv.clone().insert_axis(ndarray::Axis(1)),
];
let k = atoms.len();
let lobes = [&cu, &cnu, &cv, &cnv];
let mut logits = Array2::<f64>::zeros((n, k));
for r in 0..n {
for (a, lobe) in lobes.iter().enumerate() {
logits[[r, a]] = if lobe[r] > 1e-9 { ON } else { OFF };
}
}
let assignment = SaeAssignment::from_blocks_with_mode_and_manifolds(
logits,
coord_blocks,
vec![LatentManifold::Euclidean; k],
AssignmentMode::softmax(1.0),
)
.unwrap();
let term = SaeManifoldTerm::new(atoms, assignment).unwrap();
let rho = SaeManifoldRho::new(0.0, 0.0, vec![Array1::<f64>::zeros(1); k]);
(term, rho)
}
#[test]
fn curl_recovers_shattered_centered_circle() {
let (term, rho) = shattered_plane_term(false);
let residuals = residuals_of(&term);
let cfg = CurlConfig::default();
let cands = curl_candidates(&term, residuals.view(), &cfg);
assert!(
!cands.is_empty(),
"curl must recover the shattered circle (got no candidate)"
);
let cand = &cands[0];
let mut members = cand.members.clone();
members.sort_unstable();
members.dedup();
assert_eq!(
members,
vec![0, 1, 2, 3],
"the circle's donor set is all four rectified halves"
);
assert!(
cand.net_evidence > 0.0,
"net evidence must favour the circle"
);
let mv = StructureMove::Birth { candidate: 0 };
let seeds = vec![cand.seed.clone()];
let (born, _born_rho) = apply_structure_move_seeded(&term, &rho, &mv, &seeds).unwrap();
let circle = born.k_atoms() - 1;
assert_eq!(
born.atoms[circle].basis_kind,
SaeAtomBasisKind::Periodic,
"curl births a Periodic (circle) atom"
);
let img = atom_ambient_image(&born.atoms[circle]);
let ncols = img.ncols();
let mut center = Array1::<f64>::zeros(ncols);
for r in 0..img.nrows() {
for j in 0..ncols {
center[j] += img[[r, j]];
}
}
center.mapv_inplace(|x| x / img.nrows() as f64);
let mut min_r = f64::INFINITY;
let mut max_r = 0.0_f64;
for r in 0..img.nrows() {
let mut rr = 0.0_f64;
for j in 0..ncols {
let d = img[[r, j]] - center[j];
rr += d * d;
}
let rr = rr.sqrt();
min_r = min_r.min(rr);
max_r = max_r.max(rr);
}
assert!(
max_r > 0.0 && (max_r - min_r) / max_r < 0.1,
"born circle must trace a constant-radius ring (min={min_r:.3}, max={max_r:.3})"
);
}
#[test]
fn curl_rejects_gaussian_fill_plane() {
let (term, _rho) = shattered_plane_term(true);
let residuals = residuals_of(&term);
let cfg = CurlConfig::default();
let cands = curl_candidates(&term, residuals.view(), &cfg);
assert!(
cands.is_empty(),
"a Gaussian-fill plane must not be curled (κ ≈ 2)"
);
}
fn single_circle_term(phase_turns: &Array1<f64>) -> (SaeManifoldTerm, SaeManifoldRho) {
let n = phase_turns.len();
let p = 4usize;
let radius = 3.0_f64;
let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(3).unwrap());
let coords = phase_turns.clone().insert_axis(ndarray::Axis(1));
let (phi, jet) = evaluator.evaluate(coords.view()).unwrap();
let mut decoder = Array2::<f64>::zeros((3, p));
decoder[[2, 0]] = radius; decoder[[1, 1]] = radius; let atom = SaeManifoldAtom::new(
"circle".to_string(),
SaeAtomBasisKind::Periodic,
1,
phi,
jet,
decoder,
Array2::<f64>::eye(3),
)
.unwrap()
.with_basis_second_jet(evaluator.clone());
let logits = Array2::<f64>::from_elem((n, 1), ON);
let assignment = SaeAssignment::from_blocks_with_mode_and_manifolds(
logits,
vec![coords],
vec![LatentManifold::Circle { period: 1.0 }],
AssignmentMode::softmax(1.0),
)
.unwrap();
let term = SaeManifoldTerm::new(vec![atom], assignment).unwrap();
let rho = SaeManifoldRho::new(0.0, 0.0, vec![Array1::<f64>::zeros(1)]);
(term, rho)
}
#[test]
fn flatten_flags_diameter_and_spares_healthy_ring() {
let n = 400usize;
let diameter_phases = Array1::from_shape_fn(n, |r| if r % 2 == 0 { 0.0 } else { 0.5 });
let (diam_term, _) = single_circle_term(&diameter_phases);
let flagged = flatten_candidates(&diam_term);
assert_eq!(flagged, vec![0], "a diameter-collapsed circle must flatten");
let ring_phases = Array1::from_shape_fn(n, |r| r as f64 / n as f64);
let (ring_term, _) = single_circle_term(&ring_phases);
let flagged = flatten_candidates(&ring_term);
assert!(
flagged.is_empty(),
"a healthy full-coverage ring must NOT be flattened"
);
}
#[test]
fn curl_killer_demo_planted_circle_wins_race() {
let (term, _rho) = shattered_plane_term(false);
let residuals = residuals_of(&term);
let cands = curl_candidates(&term, residuals.view(), &CurlConfig::default());
assert!(
!cands.is_empty(),
"curl must recover the shattered circle before the race"
);
let mut members = cands[0].members.clone();
members.sort_unstable();
members.dedup();
assert_eq!(
members,
vec![0, 1, 2, 3],
"the recovered circle must claim all four rectified halves"
);
let budget = MoveBudget {
max_moves: 4,
alpha: 0.05,
};
let harvest_params = HarvestParams {
max_fusions: 0,
max_fissions: 0,
max_births: 0,
};
let run = |curl: Option<CurlConfig>| -> StructureSearchResult {
let (term, rho) = shattered_plane_term(false);
let target = 2.0 * term.try_fitted().unwrap();
let mut ledger = StructureLedger::new();
let config = RoundDriverConfig {
n_shards: 3,
budget,
max_rounds: 1,
harvest_params,
curl,
};
run_structure_search_rounds(
term,
rho,
target.view(),
config,
&mut ledger,
|t: SaeManifoldTerm, r: SaeManifoldRho, _rows: &[usize]| (t, r),
|t: SaeManifoldTerm, r: SaeManifoldRho, _rows: &[usize]| (t, r),
)
.unwrap()
};
let off = run(None);
let off_births = off
.rounds
.iter()
.flat_map(|r| r.moves.iter())
.filter(|m| matches!(m.mv, StructureMove::Birth { .. }))
.count();
assert_eq!(off_births, 0, "curl OFF (default) must inject no births");
let on = run(Some(CurlConfig::default()));
let accepted_curl_births = on
.rounds
.iter()
.flat_map(|r| r.moves.iter())
.filter(|m| {
matches!(m.mv, StructureMove::Birth { .. })
&& matches!(
m.verdict,
gam_solve::structure_search::MoveVerdict::Accepted { .. }
)
})
.count();
assert_eq!(
accepted_curl_births, 1,
"curl ON must certify exactly one circle Birth winner"
);
assert_eq!(
on.term.atoms.last().map(|a| &a.basis_kind),
Some(&SaeAtomBasisKind::Periodic),
"the accepted curl winner must be the recovered circle atom"
);
assert!(
on.structure_changed(),
"accepted curl winner must mutate the returned term"
);
let (gauss_term, _) = shattered_plane_term(true);
let gauss_residuals = residuals_of(&gauss_term);
let gauss_cands =
curl_candidates(&gauss_term, gauss_residuals.view(), &CurlConfig::default());
assert!(
gauss_cands.is_empty(),
"a Gaussian-fill plane must not be curled"
);
let n = 400usize;
let radii = Array1::<f64>::from_elem(n, 3.0);
let angles = Array1::<f64>::from_shape_fn(n, |r| {
if r % 2 == 0 {
0.0
} else {
std::f64::consts::PI
}
});
let flatten = crate::manifold::flatten_verdict(radii.view(), angles.view()).unwrap();
assert!(flatten.recommend_flatten, "diameter must flatten");
assert_eq!(flatten.residual_rank, 1, "diameter must flatten to rank 1");
}
}