changes claude never committed
This commit is contained in:
@@ -0,0 +1,12 @@
|
||||
[package]
|
||||
name = "collapse_analysis"
|
||||
version.workspace = true
|
||||
edition.workspace = true
|
||||
license.workspace = true
|
||||
|
||||
[dependencies]
|
||||
world_model = { path = "../world_model" }
|
||||
trace_model = { path = "../trace_model" }
|
||||
|
||||
[lib]
|
||||
path = "src/lib.rs"
|
||||
@@ -0,0 +1,470 @@
|
||||
//! `collapse_analysis` — the framework's compression attacks. Each attack
|
||||
//! tries to predict execution behavior with a *simpler* model. If any small
|
||||
//! model predicts above the configured thresholds, the universe has collapsed
|
||||
//! and CI must fail.
|
||||
//!
|
||||
//! The corpus is a set of (input, output) samples: `input` is the world ground
|
||||
//! truth (8 domains x (observed+hidden) lanes) and `output` is the behavior
|
||||
//! feature vector produced by the runtime. Inputs/outputs are standardized so
|
||||
//! that no single high-magnitude axis dominates the variance accounting.
|
||||
|
||||
pub mod linalg;
|
||||
|
||||
use linalg::{ols_r2, pca_scores, Mat};
|
||||
use world_model::{Hash, HIDDEN_LANES, LANES, NUM_DOMAINS};
|
||||
|
||||
/// Input columns belonging to one domain (observed + hidden lanes).
|
||||
pub const DOMAIN_BLOCK: usize = LANES + HIDDEN_LANES;
|
||||
|
||||
/// A standardized behavior corpus.
|
||||
pub struct BehaviorCorpus {
|
||||
/// `n x (NUM_DOMAINS*DOMAIN_BLOCK)` standardized input.
|
||||
pub x: Mat,
|
||||
/// `n x m` standardized output (behavior features).
|
||||
pub y: Mat,
|
||||
pub fingerprints: Vec<Hash>,
|
||||
}
|
||||
|
||||
fn standardize(rows: &[Vec<f64>]) -> Mat {
|
||||
let n = rows.len();
|
||||
let cols = if n == 0 { 0 } else { rows[0].len() };
|
||||
let mut m = Mat::zeros(n, cols);
|
||||
for r in 0..n {
|
||||
for c in 0..cols {
|
||||
m.set(r, c, rows[r][c]);
|
||||
}
|
||||
}
|
||||
// center + scale to unit std per column
|
||||
for c in 0..cols {
|
||||
let mut mean = 0.0;
|
||||
for r in 0..n {
|
||||
mean += m.at(r, c);
|
||||
}
|
||||
mean /= n.max(1) as f64;
|
||||
let mut var = 0.0;
|
||||
for r in 0..n {
|
||||
var += (m.at(r, c) - mean).powi(2);
|
||||
}
|
||||
var /= n.max(1) as f64;
|
||||
let sd = var.sqrt();
|
||||
let inv = if sd > 1e-9 { 1.0 / sd } else { 0.0 };
|
||||
for r in 0..n {
|
||||
let v = (m.at(r, c) - mean) * inv;
|
||||
m.set(r, c, v);
|
||||
}
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
impl BehaviorCorpus {
|
||||
pub fn build(inputs: Vec<Vec<f64>>, outputs: Vec<Vec<f64>>, fingerprints: Vec<Hash>) -> Self {
|
||||
BehaviorCorpus {
|
||||
x: standardize(&inputs),
|
||||
y: standardize(&outputs),
|
||||
fingerprints,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn n(&self) -> usize {
|
||||
self.x.rows
|
||||
}
|
||||
|
||||
fn select(&self, cols: &[usize]) -> Mat {
|
||||
let mut m = Mat::zeros(self.x.rows, cols.len());
|
||||
for r in 0..self.x.rows {
|
||||
for (j, &c) in cols.iter().enumerate() {
|
||||
m.set(r, j, self.x.at(r, c));
|
||||
}
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
fn domain_cols(domain: usize) -> Vec<usize> {
|
||||
(domain * DOMAIN_BLOCK..(domain + 1) * DOMAIN_BLOCK).collect()
|
||||
}
|
||||
|
||||
/// R² of predicting output from the top-`k` PCA factors of the full input.
|
||||
pub fn predict_k_factor(&self, k: usize) -> f64 {
|
||||
if self.n() == 0 {
|
||||
return 0.0;
|
||||
}
|
||||
let scores = pca_scores(&self.x, k);
|
||||
ols_r2(&scores, &self.y)
|
||||
}
|
||||
|
||||
/// Best R² obtainable using just a single domain's input block.
|
||||
pub fn max_single_domain(&self) -> (usize, f64) {
|
||||
let mut best = (0usize, 0.0);
|
||||
for d in 0..NUM_DOMAINS {
|
||||
let sub = self.select(&Self::domain_cols(d));
|
||||
let r2 = ols_r2(&sub, &self.y);
|
||||
if r2 > best.1 {
|
||||
best = (d, r2);
|
||||
}
|
||||
}
|
||||
best
|
||||
}
|
||||
|
||||
/// Best R² obtainable using any pair of domain blocks.
|
||||
pub fn max_pair(&self) -> ((usize, usize), f64) {
|
||||
let mut best = ((0usize, 1usize), 0.0);
|
||||
for a in 0..NUM_DOMAINS {
|
||||
for b in (a + 1)..NUM_DOMAINS {
|
||||
let mut cols = Self::domain_cols(a);
|
||||
cols.extend(Self::domain_cols(b));
|
||||
let sub = self.select(&cols);
|
||||
let r2 = ols_r2(&sub, &self.y);
|
||||
if r2 > best.1 {
|
||||
best = ((a, b), r2);
|
||||
}
|
||||
}
|
||||
}
|
||||
best
|
||||
}
|
||||
|
||||
fn full_r2(&self) -> f64 {
|
||||
ols_r2(&self.x, &self.y)
|
||||
}
|
||||
}
|
||||
|
||||
/// A compressed model's predictive power and information loss.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct CompressedModel {
|
||||
pub predicts: f64,
|
||||
pub info_loss: f64,
|
||||
pub detail: String,
|
||||
}
|
||||
|
||||
/// The collapse-attack trait (per spec).
|
||||
pub trait CollapseAttack {
|
||||
fn name(&self) -> &'static str;
|
||||
fn compress(&self, corpus: &BehaviorCorpus) -> CompressedModel;
|
||||
}
|
||||
|
||||
/// One attack's outcome.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct CollapseReport {
|
||||
pub attack: String,
|
||||
pub predicts: f64,
|
||||
pub info_loss: f64,
|
||||
pub detail: String,
|
||||
}
|
||||
|
||||
macro_rules! attack {
|
||||
($name:ident, $label:expr, $body:expr) => {
|
||||
pub struct $name;
|
||||
impl CollapseAttack for $name {
|
||||
fn name(&self) -> &'static str {
|
||||
$label
|
||||
}
|
||||
fn compress(&self, corpus: &BehaviorCorpus) -> CompressedModel {
|
||||
let f: fn(&BehaviorCorpus) -> CompressedModel = $body;
|
||||
f(corpus)
|
||||
}
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
fn model(predicts: f64, detail: &str) -> CompressedModel {
|
||||
CompressedModel {
|
||||
predicts,
|
||||
info_loss: (1.0 - predicts).clamp(0.0, 1.0),
|
||||
detail: detail.to_string(),
|
||||
}
|
||||
}
|
||||
|
||||
attack!(DomainRemoval, "domain_removal", |c| {
|
||||
// Best prediction achievable while *removing* each domain in turn.
|
||||
let mut best = 0.0;
|
||||
for d in 0..NUM_DOMAINS {
|
||||
let cols: Vec<usize> = (0..NUM_DOMAINS)
|
||||
.filter(|&x| x != d)
|
||||
.flat_map(BehaviorCorpus::domain_cols)
|
||||
.collect();
|
||||
let sub = c.select(&cols);
|
||||
best = f64::max(best, ols_r2(&sub, &c.y));
|
||||
}
|
||||
model(best, "predict with one domain removed")
|
||||
});
|
||||
|
||||
attack!(DomainMerging, "domain_merging", |c| {
|
||||
// Merge each pair into a summed block; best prediction over pairs.
|
||||
let mut best = 0.0;
|
||||
for a in 0..NUM_DOMAINS {
|
||||
for b in (a + 1)..NUM_DOMAINS {
|
||||
let mut merged = Mat::zeros(c.x.rows, (NUM_DOMAINS - 1) * DOMAIN_BLOCK);
|
||||
for r in 0..c.x.rows {
|
||||
let mut out_col = 0;
|
||||
for d in 0..NUM_DOMAINS {
|
||||
if d == b {
|
||||
continue;
|
||||
}
|
||||
for l in 0..DOMAIN_BLOCK {
|
||||
let mut v = c.x.at(r, d * DOMAIN_BLOCK + l);
|
||||
if d == a {
|
||||
v += c.x.at(r, b * DOMAIN_BLOCK + l);
|
||||
}
|
||||
merged.set(r, out_col, v);
|
||||
out_col += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
best = f64::max(best, ols_r2(&merged, &c.y));
|
||||
}
|
||||
}
|
||||
model(best, "predict with two domains merged")
|
||||
});
|
||||
|
||||
attack!(ConstantFolding, "constant_folding", |_c| {
|
||||
// Folding the world to constants leaves no predictive features at all.
|
||||
model(0.0, "world folded to constants")
|
||||
});
|
||||
|
||||
attack!(CausalEdgeDeletion, "causal_edge_deletion", |c| {
|
||||
// Keep only each domain's first observed lane (no cross-domain structure).
|
||||
let cols: Vec<usize> = (0..NUM_DOMAINS).map(|d| d * DOMAIN_BLOCK).collect();
|
||||
let sub = c.select(&cols);
|
||||
model(ols_r2(&sub, &c.y), "diagonal-only features")
|
||||
});
|
||||
|
||||
attack!(StateAliasing, "state_aliasing", |c| {
|
||||
// Alias all domains into a single aggregate column.
|
||||
let mut agg = Mat::zeros(c.x.rows, 1);
|
||||
for r in 0..c.x.rows {
|
||||
let mut s = 0.0;
|
||||
for col in 0..c.x.cols {
|
||||
s += c.x.at(r, col);
|
||||
}
|
||||
agg.set(r, 0, s);
|
||||
}
|
||||
model(ols_r2(&agg, &c.y), "single aliased aggregate")
|
||||
});
|
||||
|
||||
attack!(LatentFactorModeling, "latent_factor_modeling", |c| {
|
||||
model(c.predict_k_factor(4), "top-4 latent factors")
|
||||
});
|
||||
|
||||
attack!(BehaviorClustering, "behavior_clustering", |c| {
|
||||
let recon = kmeans_reconstruct(&c.y, 4);
|
||||
let scores = pca_scores(&c.x, 4);
|
||||
model(ols_r2(&scores, &recon), "4-cluster behavior model")
|
||||
});
|
||||
|
||||
attack!(SurrogatePrediction, "surrogate_prediction", |c| {
|
||||
model(c.full_r2(), "full linear surrogate")
|
||||
});
|
||||
|
||||
attack!(TemporalFlattening, "temporal_flattening", |c| {
|
||||
// Drop hidden lanes (time-carrying state); predict from observed only.
|
||||
let cols: Vec<usize> = (0..NUM_DOMAINS)
|
||||
.flat_map(|d| (0..LANES).map(move |l| d * DOMAIN_BLOCK + l))
|
||||
.collect();
|
||||
let sub = c.select(&cols);
|
||||
model(ols_r2(&sub, &c.y), "time-flattened (observed lanes only)")
|
||||
});
|
||||
|
||||
attack!(ObservationFlattening, "observation_flattening", |c| {
|
||||
// Use only hidden lanes (collapse the observed surface).
|
||||
let cols: Vec<usize> = (0..NUM_DOMAINS)
|
||||
.flat_map(|d| (0..HIDDEN_LANES).map(move |l| d * DOMAIN_BLOCK + LANES + l))
|
||||
.collect();
|
||||
let sub = c.select(&cols);
|
||||
model(ols_r2(&sub, &c.y), "observation-flattened (hidden lanes only)")
|
||||
});
|
||||
|
||||
attack!(ExecutorIdentityErasure, "executor_identity_erasure", |c| {
|
||||
// Predict all-but-last output dim (the divergence summary) from input.
|
||||
if c.y.cols <= 1 {
|
||||
return model(0.0, "no executor dim");
|
||||
}
|
||||
let mut y2 = Mat::zeros(c.y.rows, c.y.cols - 1);
|
||||
for r in 0..c.y.rows {
|
||||
for col in 0..c.y.cols - 1 {
|
||||
y2.set(r, col, c.y.at(r, col));
|
||||
}
|
||||
}
|
||||
model(ols_r2(&c.x, &y2), "executor identity erased")
|
||||
});
|
||||
|
||||
fn kmeans_reconstruct(y: &Mat, k: usize) -> Mat {
|
||||
let n = y.rows;
|
||||
if n == 0 || k == 0 {
|
||||
return y.clone();
|
||||
}
|
||||
let k = k.min(n);
|
||||
// deterministic init: spread initial centroids across the data
|
||||
let mut centroids: Vec<Vec<f64>> = (0..k).map(|i| y.data[(i * n / k) * y.cols..(i * n / k) * y.cols + y.cols].to_vec()).collect();
|
||||
let mut assign = vec![0usize; n];
|
||||
for _ in 0..12 {
|
||||
// assign
|
||||
for r in 0..n {
|
||||
let mut best = 0;
|
||||
let mut bestd = f64::MAX;
|
||||
for (ci, cen) in centroids.iter().enumerate() {
|
||||
let mut d = 0.0;
|
||||
for col in 0..y.cols {
|
||||
d += (y.at(r, col) - cen[col]).powi(2);
|
||||
}
|
||||
if d < bestd {
|
||||
bestd = d;
|
||||
best = ci;
|
||||
}
|
||||
}
|
||||
assign[r] = best;
|
||||
}
|
||||
// update
|
||||
let mut sums = vec![vec![0.0; y.cols]; k];
|
||||
let mut counts = vec![0usize; k];
|
||||
for r in 0..n {
|
||||
counts[assign[r]] += 1;
|
||||
for col in 0..y.cols {
|
||||
sums[assign[r]][col] += y.at(r, col);
|
||||
}
|
||||
}
|
||||
for ci in 0..k {
|
||||
if counts[ci] > 0 {
|
||||
for col in 0..y.cols {
|
||||
centroids[ci][col] = sums[ci][col] / counts[ci] as f64;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut recon = Mat::zeros(n, y.cols);
|
||||
for r in 0..n {
|
||||
for col in 0..y.cols {
|
||||
recon.set(r, col, centroids[assign[r]][col]);
|
||||
}
|
||||
}
|
||||
recon
|
||||
}
|
||||
|
||||
/// All eleven required attack families.
|
||||
pub fn all_attacks() -> Vec<Box<dyn CollapseAttack>> {
|
||||
vec![
|
||||
Box::new(DomainRemoval),
|
||||
Box::new(DomainMerging),
|
||||
Box::new(ConstantFolding),
|
||||
Box::new(CausalEdgeDeletion),
|
||||
Box::new(StateAliasing),
|
||||
Box::new(LatentFactorModeling),
|
||||
Box::new(BehaviorClustering),
|
||||
Box::new(SurrogatePrediction),
|
||||
Box::new(TemporalFlattening),
|
||||
Box::new(ObservationFlattening),
|
||||
Box::new(ExecutorIdentityErasure),
|
||||
]
|
||||
}
|
||||
|
||||
/// Aggregate collapse summary against all gates.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct CollapseSummary {
|
||||
pub reports: Vec<CollapseReport>,
|
||||
pub best_1factor: f64,
|
||||
pub best_2factor: f64,
|
||||
pub best_4factor: f64,
|
||||
pub max_single_domain: f64,
|
||||
pub max_pair: f64,
|
||||
pub min_info_loss: f64,
|
||||
pub failures: Vec<String>,
|
||||
}
|
||||
|
||||
impl CollapseSummary {
|
||||
pub fn passed(&self) -> bool {
|
||||
self.failures.is_empty()
|
||||
}
|
||||
}
|
||||
|
||||
/// Run every attack and check all collapse gates.
|
||||
pub fn analyze(corpus: &BehaviorCorpus) -> CollapseSummary {
|
||||
let mut reports = Vec::new();
|
||||
let mut min_info_loss: f64 = 1.0;
|
||||
for atk in all_attacks() {
|
||||
let m = atk.compress(corpus);
|
||||
min_info_loss = min_info_loss.min(m.info_loss);
|
||||
reports.push(CollapseReport {
|
||||
attack: atk.name().to_string(),
|
||||
predicts: m.predicts,
|
||||
info_loss: m.info_loss,
|
||||
detail: m.detail,
|
||||
});
|
||||
}
|
||||
|
||||
let best_1 = corpus.predict_k_factor(1);
|
||||
let best_2 = corpus.predict_k_factor(2);
|
||||
let best_4 = corpus.predict_k_factor(4);
|
||||
let (_, max_single) = corpus.max_single_domain();
|
||||
let (_, max_pair) = corpus.max_pair();
|
||||
|
||||
let mut failures = Vec::new();
|
||||
if best_1 >= 0.40 {
|
||||
failures.push(format!("1-factor predicts {:.3} >= 0.40", best_1));
|
||||
}
|
||||
if best_2 >= 0.55 {
|
||||
failures.push(format!("2-factor predicts {:.3} >= 0.55", best_2));
|
||||
}
|
||||
if best_4 >= 0.70 {
|
||||
failures.push(format!("4-factor predicts {:.3} >= 0.70", best_4));
|
||||
}
|
||||
if max_single > 0.30 {
|
||||
failures.push(format!("single domain explains {:.3} > 0.30", max_single));
|
||||
}
|
||||
if max_pair > 0.55 {
|
||||
failures.push(format!("domain pair explains {:.3} > 0.55", max_pair));
|
||||
}
|
||||
if min_info_loss < 0.35 {
|
||||
failures.push(format!("min info loss {:.3} < 0.35", min_info_loss));
|
||||
}
|
||||
|
||||
CollapseSummary {
|
||||
reports,
|
||||
best_1factor: best_1,
|
||||
best_2factor: best_2,
|
||||
best_4factor: best_4,
|
||||
max_single_domain: max_single,
|
||||
max_pair,
|
||||
min_info_loss,
|
||||
failures,
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use world_model::{Hasher, Rng};
|
||||
|
||||
#[test]
|
||||
fn random_behavior_resists_collapse() {
|
||||
// Inputs random; outputs an avalanche hash of inputs -> no small linear
|
||||
// model should predict them.
|
||||
let mut rng = Rng::new(1);
|
||||
let mut inputs = Vec::new();
|
||||
let mut outputs = Vec::new();
|
||||
let mut fps = Vec::new();
|
||||
for _ in 0..400 {
|
||||
let inp: Vec<f64> = (0..NUM_DOMAINS * DOMAIN_BLOCK)
|
||||
.map(|_| rng.range_i64(-5000, 5000) as f64)
|
||||
.collect();
|
||||
let mut h = Hasher::new();
|
||||
for &v in &inp {
|
||||
h.write_i64(v as i64);
|
||||
}
|
||||
let base = h.finish().0;
|
||||
let out: Vec<f64> = (0..12)
|
||||
.map(|k| {
|
||||
let mut hh = Hasher::new();
|
||||
hh.write_u64(base);
|
||||
hh.write_u64(k);
|
||||
(hh.finish().0 as i64) as f64
|
||||
})
|
||||
.collect();
|
||||
fps.push(world_model::Hash(base));
|
||||
inputs.push(inp);
|
||||
outputs.push(out);
|
||||
}
|
||||
let corpus = BehaviorCorpus::build(inputs, outputs, fps);
|
||||
let summary = analyze(&corpus);
|
||||
assert!(summary.passed(), "collapse failures: {:?}", summary.failures);
|
||||
assert!(summary.best_1factor < 0.40);
|
||||
assert!(summary.min_info_loss >= 0.35);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,244 @@
|
||||
//! Minimal dense f64 linear algebra for the collapse attacks: centering,
|
||||
//! OLS (ridge-regularized) multi-output R², and PCA via power iteration.
|
||||
|
||||
/// Column-major-agnostic row-major matrix.
|
||||
#[derive(Clone)]
|
||||
pub struct Mat {
|
||||
pub rows: usize,
|
||||
pub cols: usize,
|
||||
pub data: Vec<f64>,
|
||||
}
|
||||
|
||||
impl Mat {
|
||||
pub fn zeros(rows: usize, cols: usize) -> Self {
|
||||
Mat {
|
||||
rows,
|
||||
cols,
|
||||
data: vec![0.0; rows * cols],
|
||||
}
|
||||
}
|
||||
#[inline]
|
||||
pub fn at(&self, r: usize, c: usize) -> f64 {
|
||||
self.data[r * self.cols + c]
|
||||
}
|
||||
#[inline]
|
||||
pub fn set(&mut self, r: usize, c: usize, v: f64) {
|
||||
self.data[r * self.cols + c] = v;
|
||||
}
|
||||
|
||||
pub fn col(&self, c: usize) -> Vec<f64> {
|
||||
(0..self.rows).map(|r| self.at(r, c)).collect()
|
||||
}
|
||||
|
||||
/// Subtract the mean of each column (in place). Returns the means.
|
||||
pub fn center_columns(&mut self) -> Vec<f64> {
|
||||
let mut means = vec![0.0; self.cols];
|
||||
for c in 0..self.cols {
|
||||
let mut s = 0.0;
|
||||
for r in 0..self.rows {
|
||||
s += self.at(r, c);
|
||||
}
|
||||
means[c] = s / self.rows.max(1) as f64;
|
||||
}
|
||||
for r in 0..self.rows {
|
||||
for c in 0..self.cols {
|
||||
let v = self.at(r, c) - means[c];
|
||||
self.set(r, c, v);
|
||||
}
|
||||
}
|
||||
means
|
||||
}
|
||||
|
||||
/// X^T X
|
||||
pub fn gram(&self) -> Mat {
|
||||
let p = self.cols;
|
||||
let mut g = Mat::zeros(p, p);
|
||||
for i in 0..p {
|
||||
for j in i..p {
|
||||
let mut s = 0.0;
|
||||
for r in 0..self.rows {
|
||||
s += self.at(r, i) * self.at(r, j);
|
||||
}
|
||||
g.set(i, j, s);
|
||||
g.set(j, i, s);
|
||||
}
|
||||
}
|
||||
g
|
||||
}
|
||||
}
|
||||
|
||||
/// Solve (A + λI) x = b for symmetric positive-ish A via Gauss-Jordan.
|
||||
pub fn solve_ridge(a: &Mat, b: &[f64], lambda: f64) -> Vec<f64> {
|
||||
let n = a.rows;
|
||||
let mut m = a.clone();
|
||||
for i in 0..n {
|
||||
let v = m.at(i, i) + lambda;
|
||||
m.set(i, i, v);
|
||||
}
|
||||
let mut x = b.to_vec();
|
||||
// Gaussian elimination with partial pivoting.
|
||||
for col in 0..n {
|
||||
let mut piv = col;
|
||||
let mut best = m.at(col, col).abs();
|
||||
for r in (col + 1)..n {
|
||||
let v = m.at(r, col).abs();
|
||||
if v > best {
|
||||
best = v;
|
||||
piv = r;
|
||||
}
|
||||
}
|
||||
if best < 1e-12 {
|
||||
continue;
|
||||
}
|
||||
if piv != col {
|
||||
for c in 0..n {
|
||||
let tmp = m.at(col, c);
|
||||
m.set(col, c, m.at(piv, c));
|
||||
m.set(piv, c, tmp);
|
||||
}
|
||||
x.swap(col, piv);
|
||||
}
|
||||
let d = m.at(col, col);
|
||||
for r in 0..n {
|
||||
if r != col {
|
||||
let f = m.at(r, col) / d;
|
||||
if f != 0.0 {
|
||||
for c in col..n {
|
||||
let v = m.at(r, c) - f * m.at(col, c);
|
||||
m.set(r, c, v);
|
||||
}
|
||||
x[r] -= f * x[col];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for i in 0..n {
|
||||
let d = m.at(i, i);
|
||||
if d.abs() > 1e-12 {
|
||||
x[i] /= d;
|
||||
} else {
|
||||
x[i] = 0.0;
|
||||
}
|
||||
}
|
||||
x
|
||||
}
|
||||
|
||||
/// Average R^2 of predicting each (centered) output column from the centered
|
||||
/// design matrix `x` using ridge OLS. Returns a value clamped to [0, 1].
|
||||
pub fn ols_r2(x: &Mat, y: &Mat) -> f64 {
|
||||
if x.cols == 0 || x.rows == 0 {
|
||||
return 0.0;
|
||||
}
|
||||
let g = x.gram();
|
||||
let lambda = 1e-6 * (1.0 + trace(&g) / x.cols as f64);
|
||||
let mut total_r2 = 0.0;
|
||||
let mut counted = 0;
|
||||
for oc in 0..y.cols {
|
||||
let yc = y.col(oc);
|
||||
// X^T y
|
||||
let mut xty = vec![0.0; x.cols];
|
||||
for i in 0..x.cols {
|
||||
let mut s = 0.0;
|
||||
for r in 0..x.rows {
|
||||
s += x.at(r, i) * yc[r];
|
||||
}
|
||||
xty[i] = s;
|
||||
}
|
||||
let beta = solve_ridge(&g, &xty, lambda);
|
||||
// residuals
|
||||
let mut ss_res = 0.0;
|
||||
let mut ss_tot = 0.0;
|
||||
for r in 0..x.rows {
|
||||
let mut pred = 0.0;
|
||||
for i in 0..x.cols {
|
||||
pred += x.at(r, i) * beta[i];
|
||||
}
|
||||
ss_res += (yc[r] - pred).powi(2);
|
||||
ss_tot += yc[r].powi(2);
|
||||
}
|
||||
if ss_tot > 1e-9 {
|
||||
let r2 = 1.0 - ss_res / ss_tot;
|
||||
total_r2 += r2.clamp(0.0, 1.0);
|
||||
counted += 1;
|
||||
}
|
||||
}
|
||||
if counted == 0 {
|
||||
0.0
|
||||
} else {
|
||||
total_r2 / counted as f64
|
||||
}
|
||||
}
|
||||
|
||||
fn trace(m: &Mat) -> f64 {
|
||||
(0..m.rows.min(m.cols)).map(|i| m.at(i, i)).sum()
|
||||
}
|
||||
|
||||
/// Top-`k` principal-component scores of a centered matrix `x` via power
|
||||
/// iteration with deflation. Returns an `n x k` score matrix.
|
||||
pub fn pca_scores(x: &Mat, k: usize) -> Mat {
|
||||
let p = x.cols;
|
||||
let mut cov = x.gram(); // proportional to covariance
|
||||
let kk = k.min(p);
|
||||
let mut comps: Vec<Vec<f64>> = Vec::new();
|
||||
for _ in 0..kk {
|
||||
// power iteration
|
||||
let mut v = vec![0.0; p];
|
||||
for (i, vi) in v.iter_mut().enumerate() {
|
||||
*vi = 1.0 + (i as f64) * 0.001;
|
||||
}
|
||||
normalize(&mut v);
|
||||
for _ in 0..64 {
|
||||
let mut nv = matvec(&cov, &v);
|
||||
normalize(&mut nv);
|
||||
let diff: f64 = nv.iter().zip(&v).map(|(a, b)| (a - b).abs()).sum();
|
||||
v = nv;
|
||||
if diff < 1e-9 {
|
||||
break;
|
||||
}
|
||||
}
|
||||
// eigenvalue
|
||||
let av = matvec(&cov, &v);
|
||||
let lambda: f64 = v.iter().zip(&av).map(|(a, b)| a * b).sum();
|
||||
// deflate
|
||||
for i in 0..p {
|
||||
for j in 0..p {
|
||||
let val = cov.at(i, j) - lambda * v[i] * v[j];
|
||||
cov.set(i, j, val);
|
||||
}
|
||||
}
|
||||
comps.push(v);
|
||||
}
|
||||
// scores = X * comps
|
||||
let mut scores = Mat::zeros(x.rows, kk);
|
||||
for r in 0..x.rows {
|
||||
for (cj, comp) in comps.iter().enumerate() {
|
||||
let mut s = 0.0;
|
||||
for i in 0..p {
|
||||
s += x.at(r, i) * comp[i];
|
||||
}
|
||||
scores.set(r, cj, s);
|
||||
}
|
||||
}
|
||||
scores
|
||||
}
|
||||
|
||||
fn matvec(m: &Mat, v: &[f64]) -> Vec<f64> {
|
||||
let mut out = vec![0.0; m.rows];
|
||||
for r in 0..m.rows {
|
||||
let mut s = 0.0;
|
||||
for c in 0..m.cols {
|
||||
s += m.at(r, c) * v[c];
|
||||
}
|
||||
out[r] = s;
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
fn normalize(v: &mut [f64]) {
|
||||
let n: f64 = v.iter().map(|x| x * x).sum::<f64>().sqrt();
|
||||
if n > 1e-12 {
|
||||
for x in v.iter_mut() {
|
||||
*x /= n;
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user