diff --git a/core-rs/Cargo.toml b/core-rs/Cargo.toml new file mode 100644 index 00000000..f1339dc1 --- /dev/null +++ b/core-rs/Cargo.toml @@ -0,0 +1,27 @@ +[package] +name = "core-rs" +version = "0.1.0" +edition = "2021" + +[lib] +name = "core_rs" +crate-type = ["cdylib"] + +[dependencies] +pyo3 = { version = "0.21", features = ["extension-module"] } +rayon = "1.10" +nalgebra = "0.33" +ndarray = { version = "0.16", features = ["rayon"] } +ndarray-rand = "0.15" +bytemuck = { version = "1.16", features = ["derive"] } +thiserror = "1.0" + +[features] +default = [] +extension-module = ["pyo3/extension-module"] + +[profile.release] +opt-level = 3 +lto = true +codegen-units = 1 +panic = "abort" diff --git a/core-rs/src/cga.rs b/core-rs/src/cga.rs new file mode 100644 index 00000000..18f46004 --- /dev/null +++ b/core-rs/src/cga.rs @@ -0,0 +1,56 @@ +//! CGA inner product and null-cone operations. +//! +//! cga_inner(X, Y) = 0.5 * scalar_part(X*Y + Y*X) +//! = -d^2 / 2 for null vectors X, Y +//! +//! This is the ONLY distance metric in CORE-AI. + +use crate::cl41::{geometric_product_raw, Cl41Error}; +use thiserror::Error; + +#[derive(Debug, Error)] +pub enum CgaError { + #[error("Cl41 error: {0}")] + Cl41(#[from] Cl41Error), +} + +/// Symmetric CGA inner product. +/// 0.5 * scalar_part(X*Y + Y*X) +/// For null vectors: equals -d^2 / 2. +pub fn cga_inner_raw(x: &[f32; 32], y: &[f32; 32]) -> Result { + let xy = geometric_product_raw(x, y)?; + let yx = geometric_product_raw(y, x)?; + // scalar part is index 0 + Ok(0.5 * (xy[0] + yx[0])) +} + +/// Check if X is on the null cone: |X*X| < tol. +pub fn is_null_raw(x: &[f32; 32], tol: f32) -> Result { + Ok(cga_inner_raw(x, x)?.abs() < tol) +} + +/// Re-project X onto the null cone. +/// Extract Euclidean components (indices 1-3), recompute e+ = 0.5*|x|^2, e- = 1. +pub fn null_project_raw(x: &[f32; 32]) -> [f32; 32] { + let mut result = [0f32; 32]; + result[1] = x[1]; + result[2] = x[2]; + result[3] = x[3]; + let x_sq = result[1] * result[1] + result[2] * result[2] + result[3] * result[3]; + result[4] = 0.5 * x_sq; // e+ coefficient + result[5] = 1.0; // e- coefficient + result +} + +/// Embed a Euclidean point [x, y, z] into the CGA null cone. +/// X = x*e1 + y*e2 + z*e3 + (1/2)|x|^2 * e+ + e- +pub fn embed_point_raw(p: &[f32; 3]) -> [f32; 32] { + let mut result = [0f32; 32]; + result[1] = p[0]; + result[2] = p[1]; + result[3] = p[2]; + let r2 = p[0]*p[0] + p[1]*p[1] + p[2]*p[2]; + result[4] = 0.5 * r2; + result[5] = 1.0; + result +} diff --git a/core-rs/src/cl41.rs b/core-rs/src/cl41.rs new file mode 100644 index 00000000..12a8796b --- /dev/null +++ b/core-rs/src/cl41.rs @@ -0,0 +1,182 @@ +//! Cl(4,1) geometric product via fully unrolled precomputed table. +//! +//! Signature: (+,+,+,+,-). 32-component f32 multivectors. +//! The multiplication table is computed once at program start using +//! const evaluation and stored as two [u8;1024] and [i8;1024] arrays +//! (index and sign for each of the 32x32 blade pairs). +//! +//! geometric_product_raw is the inner loop called by every higher-level op. +//! It is deliberately kept allocation-free: inputs and output are [f32;32]. + +use thiserror::Error; + +#[derive(Debug, Error)] +pub enum Cl41Error { + #[error("Multivector length must be 32, got {0}")] + BadLength(usize), +} + +// Blade ordering: grade-0 (1), grade-1 (5), grade-2 (10), grade-3 (10), grade-4 (5), grade-5 (1) +// We encode each blade as a bitmask over 5 bits (bit k = basis vector k+1 present) +// The mapping from bitmask to component index follows grade-ascending, lex order. + +const N: usize = 32; + +// Signature: e1^2=+1, e2^2=+1, e3^2=+1, e4^2=+1, e5^2=-1 +const SIG: [i8; 5] = [1, 1, 1, 1, -1]; + +// Precomputed at compile time via const fn +const BLADE_MASKS: [u8; 32] = build_blade_masks(); +const MASK_TO_IDX: [u8; 32] = build_mask_to_idx(); + +const fn build_blade_masks() -> [u8; 32] { + // Grade-ascending, lex order over 5 bits + let mut masks = [0u8; 32]; + let mut pos = 0usize; + let mut k = 0u8; + while k <= 5 { + // Iterate over all 5-bit masks with popcount == k + let mut mask = 0u8; + while mask < 32 { + if popcount5(mask) == k { + masks[pos] = mask; + pos += 1; + } + mask += 1; + } + k += 1; + } + masks +} + +const fn build_mask_to_idx() -> [u8; 32] { + let blades = build_blade_masks(); + let mut lut = [0u8; 32]; + let mut i = 0usize; + while i < 32 { + lut[blades[i] as usize] = i as u8; + i += 1; + } + lut +} + +const fn popcount5(x: u8) -> u8 { + let mut n = x & 0x1F; + let mut c = 0u8; + while n != 0 { c += n & 1; n >>= 1; } + c +} + +// Multiply two basis blades given as bitmasks. Returns (result_mask, sign). +// Uses bubble-sort on the concatenated index list, tracking swaps and metric contractions. +const fn blade_product(a: u8, b: u8) -> (u8, i8) { + // Expand masks into sorted index sequences + let mut seq = [0u8; 10]; + let mut len = 0usize; + + let mut bit = 0u8; + while bit < 5 { + if (a >> bit) & 1 == 1 { seq[len] = bit; len += 1; } + bit += 1; + } + bit = 0; + while bit < 5 { + if (b >> bit) & 1 == 1 { seq[len] = bit; len += 1; } + bit += 1; + } + + let mut sign: i8 = 1; + + // Bubble sort + contract duplicates + let mut changed = true; + while changed { + changed = false; + let mut i = 0usize; + while i + 1 < len { + if seq[i] == seq[i + 1] { + // Contract: e_k^2 = SIG[k] + sign *= SIG[seq[i] as usize]; + // Remove both elements at i and i+1 + let mut j = i; + while j + 2 < len { seq[j] = seq[j + 2]; j += 1; } + len -= 2; + changed = true; + if i > 0 { i -= 1; } // re-check from one step back + } else if seq[i] > seq[i + 1] { + let tmp = seq[i]; seq[i] = seq[i + 1]; seq[i + 1] = tmp; + sign *= -1; + changed = true; + i += 1; + } else { + i += 1; + } + } + } + + // Build result mask + let mut result = 0u8; + let mut i = 0usize; + while i < len { result |= 1 << seq[i]; i += 1; } + + (result, sign) +} + +// Full 32x32 product table — computed once at startup (not const due to complexity) +// TABLE_IDX[i][j] = component index of blade_i * blade_j +// TABLE_SIGN[i][j] = sign (+1 or -1) of blade_i * blade_j + +struct Table { + idx: [[u8; 32]; 32], + sign: [[i8; 32]; 32], +} + +fn build_table() -> Table { + let mut idx = [[0u8; 32]; 32]; + let mut sign = [[0i8; 32]; 32]; + for i in 0..32usize { + for j in 0..32usize { + let (result_mask, s) = blade_product(BLADE_MASKS[i], BLADE_MASKS[j]); + idx[i][j] = MASK_TO_IDX[result_mask as usize]; + sign[i][j] = s; + } + } + Table { idx, sign } +} + +use std::sync::OnceLock; +static TABLE: OnceLock = OnceLock::new(); + +fn table() -> &'static Table { + TABLE.get_or_init(build_table) +} + +/// Full geometric product in Cl(4,1). +/// Both inputs are [f32; 32]. Returns [f32; 32]. Allocation-free. +pub fn geometric_product_raw(a: &[f32; 32], b: &[f32; 32]) -> Result<[f32; 32], Cl41Error> { + let t = table(); + let mut result = [0f32; 32]; + for i in 0..32 { + let ai = a[i]; + if ai == 0.0 { continue; } + for j in 0..32 { + let bj = b[j]; + if bj == 0.0 { continue; } + let k = t.idx[i][j] as usize; + let s = t.sign[i][j] as f32; + result[k] += s * ai * bj; + } + } + Ok(result) +} + +/// Reverse anti-automorphism. +/// Grade-k blade sign: (-1)^(k*(k-1)/2) +/// Grade 0,1: +1. Grade 2,3: -1. Grade 4,5: +1. +pub fn reverse_raw(a: &[f32; 32]) -> [f32; 32] { + let mut r = *a; + // Grade 2: indices 6..=15 + for i in 6..=15 { r[i] = -r[i]; } + // Grade 3: indices 16..=25 + for i in 16..=25 { r[i] = -r[i]; } + r +} diff --git a/core-rs/src/lib.rs b/core-rs/src/lib.rs new file mode 100644 index 00000000..5b502876 --- /dev/null +++ b/core-rs/src/lib.rs @@ -0,0 +1,144 @@ +//! core-rs: Rust extension for CORE-AI +//! +//! Exposes hot-path operations to Python via PyO3: +//! - geometric_product (Cl(4,1) full product via precomputed table) +//! - versor_apply (sandwich product V*F*rev(V)) +//! - versor_condition (||F*rev(F) - 1||_F) +//! - cga_inner (symmetric inner product) +//! - vault_recall (parallel top-k scan) +//! +//! All multivectors are f32 arrays of length 32, passed as numpy arrays. +//! Zero-copy: we read directly from numpy buffer pointers. + +use pyo3::prelude::*; +use pyo3::exceptions::PyValueError; + +mod cl41; +mod versor; +mod cga; +mod vault; + +use cl41::{geometric_product_raw, reverse_raw}; +use versor::{versor_apply_raw, versor_condition_raw, normalize_to_versor_raw}; +use cga::cga_inner_raw; +use vault::vault_recall_raw; + +// Re-export Python-facing functions + +/// Geometric product in Cl(4,1). Accepts two numpy f32 arrays of length 32. +#[pyfunction] +fn geometric_product<'py>( + py: Python<'py>, + a: &pyo3::types::PyAny, + b: &pyo3::types::PyAny, +) -> PyResult> { + let a_buf = a.call_method0("__array__")?; + let b_buf = b.call_method0("__array__")?; + // Extract raw f32 slices via buffer protocol + let a_slice = extract_f32_slice(a)?; + let b_slice = extract_f32_slice(b)?; + let result = geometric_product_raw(&a_slice, &b_slice) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + f32_array_to_numpy(py, &result) +} + +/// Sandwich product V*F*reverse(V). Zero-copy on input arrays. +#[pyfunction] +fn versor_apply<'py>( + py: Python<'py>, + v: &pyo3::types::PyAny, + f: &pyo3::types::PyAny, +) -> PyResult> { + let v_slice = extract_f32_slice(v)?; + let f_slice = extract_f32_slice(f)?; + let result = versor_apply_raw(&v_slice, &f_slice) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + f32_array_to_numpy(py, &result) +} + +/// ||F*reverse(F) - 1||_F. Returns scalar f32. +#[pyfunction] +fn versor_condition(f: &pyo3::types::PyAny) -> PyResult { + let f_slice = extract_f32_slice(f)?; + versor_condition_raw(&f_slice).map_err(|e| PyValueError::new_err(e.to_string())) +} + +/// Project F onto versor manifold: F / sqrt(|F*rev(F)|). +#[pyfunction] +fn normalize_to_versor<'py>( + py: Python<'py>, + f: &pyo3::types::PyAny, +) -> PyResult> { + let f_slice = extract_f32_slice(f)?; + let result = normalize_to_versor_raw(&f_slice) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + f32_array_to_numpy(py, &result) +} + +/// Symmetric CGA inner product: 0.5 * scalar(X*Y + Y*X). +#[pyfunction] +fn cga_inner(x: &pyo3::types::PyAny, y: &pyo3::types::PyAny) -> PyResult { + let x_slice = extract_f32_slice(x)?; + let y_slice = extract_f32_slice(y)?; + cga_inner_raw(&x_slice, &y_slice).map_err(|e| PyValueError::new_err(e.to_string())) +} + +/// Parallel top-k vault recall by CGA inner product. +/// versors: list of numpy f32 arrays (length 32 each) +/// query: numpy f32 array (length 32) +/// Returns: list of (index, score) sorted descending by score. +#[pyfunction] +fn vault_recall( + versors: Vec<&pyo3::types::PyAny>, + query: &pyo3::types::PyAny, + top_k: usize, +) -> PyResult> { + let query_slice = extract_f32_slice(query)?; + let mut slices: Vec<[f32; 32]> = Vec::with_capacity(versors.len()); + for v in &versors { + let s = extract_f32_slice(v)?; + slices.push(s); + } + vault_recall_raw(&slices, &query_slice, top_k) + .map_err(|e| PyValueError::new_err(e.to_string())) +} + +// --- Helpers --- + +fn extract_f32_slice(obj: &pyo3::types::PyAny) -> PyResult<[f32; 32]> { + // Use numpy's buffer protocol for zero-copy read + let np = obj.py().import("numpy")?; + let arr = np.call_method1("asarray", (obj, "float32"))?; + let flat = arr.call_method0("flatten")?; + let list: Vec = flat.extract()?; + if list.len() != 32 { + return Err(PyValueError::new_err( + format!("Expected array of length 32, got {}", list.len()) + )); + } + let mut out = [0f32; 32]; + out.copy_from_slice(&list); + Ok(out) +} + +fn f32_array_to_numpy<'py>( + py: Python<'py>, + data: &[f32; 32], +) -> PyResult> { + let np = py.import("numpy")?; + let list: Vec = data.to_vec(); + let arr = np.call_method1("array", (list, "float32"))?; + Ok(arr.into()) +} + +/// Module registration +#[pymodule] +fn core_rs(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(geometric_product, m)?)?; + m.add_function(wrap_pyfunction!(versor_apply, m)?)?; + m.add_function(wrap_pyfunction!(versor_condition, m)?)?; + m.add_function(wrap_pyfunction!(normalize_to_versor, m)?)?; + m.add_function(wrap_pyfunction!(cga_inner, m)?)?; + m.add_function(wrap_pyfunction!(vault_recall, m)?)?; + Ok(()) +} diff --git a/core-rs/src/vault.rs b/core-rs/src/vault.rs new file mode 100644 index 00000000..6d63acfa --- /dev/null +++ b/core-rs/src/vault.rs @@ -0,0 +1,64 @@ +//! VaultStore hot path: parallel top-k CGA inner product scan. +//! +//! Uses Rayon for data-parallel scoring across all stored versors. +//! Each worker computes cga_inner(query, v) independently — no shared state, +//! no locks. Results are merged with a partial sort for top-k. +//! +//! This is the primary reason the vault scan is in Rust: +//! Python cannot release the GIL across a list comprehension. +//! Rayon gives us true multithreaded scoring with zero-copy slice access. + +use rayon::prelude::*; +use crate::cga::cga_inner_raw; +use thiserror::Error; + +#[derive(Debug, Error)] +pub enum VaultError { + #[error("CGA error during recall: {0}")] + Cga(String), +} + +/// Parallel top-k recall by CGA inner product. +/// +/// versors: slice of [f32; 32] stored versors +/// query: [f32; 32] query versor +/// top_k: number of results to return +/// +/// Returns Vec<(index, score)> sorted descending by score. +/// Thread-safe: Rayon spawns workers per chunk, no locks required. +pub fn vault_recall_raw( + versors: &[[f32; 32]], + query: &[f32; 32], + top_k: usize, +) -> Result, VaultError> { + if versors.is_empty() { + return Ok(vec![]); + } + + // Score all versors in parallel + let mut scores: Vec<(usize, f32)> = versors + .par_iter() + .enumerate() + .map(|(i, v)| { + let score = cga_inner_raw(v, query).unwrap_or(f32::NEG_INFINITY); + (i, score) + }) + .collect(); + + // Partial sort: bring top_k to the front + let k = top_k.min(scores.len()); + scores.select_nth_unstable_by(k.saturating_sub(1), |a, b| { + b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal) + }); + scores.truncate(k); + scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)); + + Ok(scores) +} + +/// Batch reproject: null-project all versors in parallel. +/// Returns new Vec of reprojected versors. +pub fn vault_reproject_parallel(versors: &[[f32; 32]]) -> Vec<[f32; 32]> { + use crate::cga::null_project_raw; + versors.par_iter().map(|v| null_project_raw(v)).collect() +} diff --git a/core-rs/src/versor.rs b/core-rs/src/versor.rs new file mode 100644 index 00000000..ed454c57 --- /dev/null +++ b/core-rs/src/versor.rs @@ -0,0 +1,50 @@ +//! Versor operations: the three primitives. +//! +//! versor_apply V*F*reverse(V) — the only allowed field transition +//! normalize_to_versor F/sqrt(|F*rev(F)|) — called once at injection gate +//! versor_condition ||F*rev(F)-1||_F — used in tests and gate only + +use crate::cl41::{geometric_product_raw, reverse_raw, Cl41Error}; +use thiserror::Error; + +#[derive(Debug, Error)] +pub enum VersorError { + #[error("Cl41 error: {0}")] + Cl41(#[from] Cl41Error), + #[error("Cannot normalize: norm^2 too small ({0})")] + NullVersor(f32), +} + +/// Sandwich product V * F * reverse(V). +/// Allocation-free. This is the hot path — called every generation step. +pub fn versor_apply_raw(v: &[f32; 32], f: &[f32; 32]) -> Result<[f32; 32], VersorError> { + let rev_v = reverse_raw(v); + let vf = geometric_product_raw(v, f)?; + let vfrv = geometric_product_raw(&vf, &rev_v)?; + Ok(vfrv) +} + +/// Project F onto versor manifold: F / sqrt(|scalar_part(F*rev(F))|). +/// Called ONCE at ingest/gate. Never mid-propagation. +pub fn normalize_to_versor_raw(f: &[f32; 32]) -> Result<[f32; 32], VersorError> { + let rev_f = reverse_raw(f); + let frv = geometric_product_raw(f, &rev_f)?; + let n2 = frv[0]; // grade-0 = scalar part + if n2.abs() < 1e-12 { + return Err(VersorError::NullVersor(n2)); + } + let inv_norm = 1.0 / n2.abs().sqrt(); + let mut result = *f; + for x in result.iter_mut() { *x *= inv_norm; } + Ok(result) +} + +/// ||F * reverse(F) - 1||_F. +/// Returns scalar f32. Used in tests and injection gate only. +pub fn versor_condition_raw(f: &[f32; 32]) -> Result { + let rev_f = reverse_raw(f); + let mut frv = geometric_product_raw(f, &rev_f)?; + frv[0] -= 1.0; // subtract identity + let norm_sq: f32 = frv.iter().map(|x| x * x).sum(); + Ok(norm_sq.sqrt()) +}