diff --git a/srctmp/channel.rs b/srctmp/channel.rs new file mode 100644 index 0000000..c033cf9 --- /dev/null +++ b/srctmp/channel.rs @@ -0,0 +1,192 @@ +use crate::{Gf2, LdpcError, Llr, Result}; + +// Trait Channel + +pub trait Channel: Send + Sync { + fn transmit(&self, codeword: &[Gf2], rng: &mut impl rand::Rng) -> Vec; + fn capacity(&self) -> f64; +} + +// Canal AWGN +// Modulation BPSK : 0 -> +1.0, 1 -> -1.0 +// Signal reçu : y = x + n, n eq N(0, sig²) +// LLR optimal : L(y) = 2y/sig² +// sig² = N_0/2 = 1/(2 R SNR_lin) + +#[derive(Debug, Clone)] +pub struct AwgnChannel { + pub snr_db: f64, + sigma: f64, +} + +impl AwgnChannel { + pub fn new(snr_db: f64, code_rate: f64) -> Result { + if !(0.0..1.0).contains(&code_rate) { + return Err(LdpcError::OutOfRange("code_rate ∈ ]0, 1[".into())); + } + let snr_lin = 10.0_f64.powf(snr_db / 10.0); + let sigma = (1.0 / (2.0 * code_rate * snr_lin)).sqrt(); + Ok(Self { snr_db, sigma }) + } + + pub fn sigma(&self) -> f64 { + self.sigma + } + pub fn snr_linear(&self) -> f64 { + 10.0_f64.powf(self.snr_db / 10.0) + } + + #[inline] + pub fn llr_from_received(y: f64, sigma: f64) -> Llr { + 2.0 * y / (sigma * sigma) + } +} + +impl Channel for AwgnChannel { + fn transmit(&self, codeword: &[Gf2], rng: &mut impl rand::Rng) -> Vec { + use rand_distr::{Distribution, Normal}; + let normal = Normal::new(0.0, self.sigma).unwrap(); + codeword + .iter() + .map(|&b| { + let x = if b == 0 { 1.0_f64 } else { -1.0_f64 }; + let y = x + normal.sample(rng); + Self::llr_from_received(y, self.sigma) + }) + .collect() + } + + fn capacity(&self) -> f64 { + // Capacité BPSK-AWGN par Monte-Carlo + use rand_distr::{Distribution, Normal}; + let mut rng = rand::thread_rng(); + let normal = Normal::new(0.0, self.sigma).unwrap(); + let n_samples = 10_000usize; + let mut sum = 0.0f64; + for _ in 0..n_samples { + let n: f64 = normal.sample(&mut rng); + let y = 1.0 + n; // bit 0 transmis (x=+1) + let llr = Self::llr_from_received(y, self.sigma); + // I = 1 - E[log2(1 + exp(-2y/sig²))] + sum += (1.0 + (-llr).exp()).log2(); + } + 1.0 - sum / n_samples as f64 + } +} + +// Canal BSC +// Chaque bit inversé avec probabilité p +// LLR : +-log((1-p)/p) selon le bit reçu + +#[derive(Debug, Clone)] +pub struct BscChannel { + pub crossover_prob: f64, + llr_magnitude: Llr, +} + +impl BscChannel { + pub fn new(crossover_prob: f64) -> Result { + if crossover_prob <= 0.0 || crossover_prob >= 0.5 { + return Err(LdpcError::OutOfRange("p ∈ ]0, 0.5[".into())); + } + let llr_magnitude = ((1.0 - crossover_prob) / crossover_prob).ln(); + Ok(Self { + crossover_prob, + llr_magnitude, + }) + } +} + +impl Channel for BscChannel { + fn transmit(&self, codeword: &[Gf2], rng: &mut impl rand::Rng) -> Vec { + codeword + .iter() + .map(|&b| { + let rcv = if rng.gen::() < self.crossover_prob { + b ^ 1 + } else { + b + }; + if rcv == 0 { + self.llr_magnitude + } else { + -self.llr_magnitude + } + }) + .collect() + } + + // C_BSC(p) = 1 - Hb(p) avec Hb = entropie binaire + fn capacity(&self) -> f64 { + let p = self.crossover_prob; + let hb = -p * p.log2() - (1.0 - p) * (1.0 - p).log2(); + 1.0 - hb + } +} + +// Canal BEC +// Bit effacé avec probabilité ε -> LLR = 0 (incertitude totale) +// Bit reçu correctement → LLR = +-CERTAIN_LLR (grand mais fini) + +#[derive(Debug, Clone)] +pub struct BecChannel { + pub erasure_prob: f64, +} + +impl BecChannel { + pub fn new(erasure_prob: f64) -> Result { + if erasure_prob <= 0.0 || erasure_prob >= 1.0 { + return Err(LdpcError::OutOfRange("ε ∈ ]0, 1[".into())); + } + Ok(Self { erasure_prob }) + } + + // Grand mais fini pour stabilité numérique + const CERTAIN_LLR: Llr = 100.0; +} + +impl Channel for BecChannel { + fn transmit(&self, codeword: &[Gf2], rng: &mut impl rand::Rng) -> Vec { + codeword + .iter() + .map(|&b| { + if rng.gen::() < self.erasure_prob { + 0.0 // Effacement + } else if b == 0 { + Self::CERTAIN_LLR + } else { + -Self::CERTAIN_LLR + } + }) + .collect() + } + + /// C_BEC(ε) = 1 - ε + fn capacity(&self) -> f64 { + 1.0 - self.erasure_prob + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// +// #[test] +// fn test_awgn_llr_formula() { +// assert!((AwgnChannel::llr_from_received(1.0, 1.0) - 2.0).abs() < 1e-12); +// assert!((AwgnChannel::llr_from_received(-1.0, 1.0) + 2.0).abs() < 1e-12); +// } +// +// #[test] +// fn test_bec_capacity() { +// let bec = BecChannel::new(0.3).unwrap(); +// assert!((bec.capacity() - 0.7).abs() < 1e-12); +// } +// +// #[test] +// fn test_bsc_capacity_bounds() { +// let bsc = BscChannel::new(0.1).unwrap(); +// let cap = bsc.capacity(); +// assert!(cap > 0.0 && cap < 1.0); +// } +// } diff --git a/srctmp/code.rs b/srctmp/code.rs new file mode 100644 index 0000000..3478088 --- /dev/null +++ b/srctmp/code.rs @@ -0,0 +1,453 @@ +use crate::graph::TannerGraph; +use crate::matrix::{DenseMatrixGF2, SparseMatrixGF2}; +use crate::{Gf2, LdpcError, Result}; + +#[derive(Debug, Clone)] +pub struct LdpcParams { + pub n: usize, + pub k: usize, + pub topology: CodeTopology, + pub generation: GenerationMethod, + pub seed: Option, +} + +impl LdpcParams { + pub fn rate(&self) -> f64 { + self.k as f64 / self.n as f64 + } + pub fn m(&self) -> usize { + self.n - self.k + } + + pub fn validate(&self) -> Result<()> { + if self.k >= self.n { + return Err(LdpcError::InvalidParameters("k doit être < n".into())); + } + if self.n < 4 { + return Err(LdpcError::InvalidParameters("n trop petit".into())); + } + self.topology.validate(self) + } +} + +// Topologie + +#[derive(Debug, Clone)] +pub enum CodeTopology { + // Chaque var-node a degré wc, chaque check-node a degré wr + // Condition nécessaire : n * wc == m * wr. + Regular { wc: usize, wr: usize }, + // Degrés définis par des polynômes de distribution (density evolution) + // lambda[i] = fraction d'arêtes sur des var-nodes de degré i+2 + // rho[i] = fraction d'arêtes sur des check-nodes de degré i+2 + Irregular { lambda: Vec, rho: Vec }, +} + +impl CodeTopology { + fn validate(&self, params: &LdpcParams) -> Result<()> { + match self { + CodeTopology::Regular { wc, wr } => { + if params.n * wc != params.m() * wr { + return Err(LdpcError::InvalidParameters(format!( + "n*wc ({}) doit égaler m*wr ({}) pour un code régulier", + params.n * wc, + params.m() * wr + ))); + } + } + CodeTopology::Irregular { lambda, rho } => { + let sl: f64 = lambda.iter().sum(); + let sr: f64 = rho.iter().sum(); + if (sl - 1.0).abs() > 1e-6 || (sr - 1.0).abs() > 1e-6 { + return Err(LdpcError::InvalidParameters( + "Les polynômes lambda et rho doivent sommer à 1".into(), + )); + } + } + } + Ok(()) + } +} + +// Méthode de génération + +#[derive(Debug, Clone)] +pub enum GenerationMethod { + // H = [H1 | H2 | ... | Hwc]^T, H1 régulière, Hi = permutation de H1 + Gallager, + // Ajout de colonnes de poids fixe, rejet si cycle-4 créé + MacKayNeal { max_attempts: usize }, + // Progressive Edge Growth : maximise le girth local arête par arête + Peg, +} + +// Forme systématique + +#[derive(Debug, Clone)] +pub struct SystematicForm { + // G = [I_k | P], dense car P est généralement pleine + pub g: DenseMatrixGF2, + // Permutation de colonnes appliquée à H pour obtenir [A | I_m] + pub col_perm: Vec, + // Permutation inverse pour réordonner le mot de code + pub col_perm_inv: Vec, +} + +// Structure principale + +#[derive(Debug, Clone)] +pub struct LdpcCode { + pub params: LdpcParams, + pub h: SparseMatrixGF2, + pub graph: TannerGraph, + pub systematic_form: Option, +} + +impl LdpcCode { + pub fn new(params: LdpcParams) -> Result { + params.validate()?; + let mut rng = build_rng(params.seed); + let h = match ¶ms.generation { + GenerationMethod::Gallager => generate_gallager(¶ms, &mut rng)?, + GenerationMethod::MacKayNeal { max_attempts } => { + generate_mackay_neal(¶ms, *max_attempts, &mut rng)? + } + GenerationMethod::Peg => generate_peg(¶ms)?, + }; + let graph = TannerGraph::from_matrix(&h); + Ok(Self { + params, + h, + graph, + systematic_form: None, + }) + } + + pub fn from_matrix(h: SparseMatrixGF2, k: usize) -> Result { + let n = h.cols; + let params = LdpcParams { + n, + k, + topology: CodeTopology::Regular { wc: 0, wr: 0 }, // inconnu + generation: GenerationMethod::Gallager, + seed: None, + }; + let graph = TannerGraph::from_matrix(&h); + Ok(Self { + params, + h, + graph, + systematic_form: None, + }) + } + + // Calcule G par Gauss-Jordan sur H (to cache) + // Appelé automatiquement par SystematicEncoder si nécessaire + pub fn compute_systematic_form(&mut self) -> Result<()> { + if self.systematic_form.is_some() { + return Ok(()); + } + + let dense = self.h.to_dense(); + // On travaille sur [H | I_m] pour tracker les opérations de lignes + let m = self.params.m(); + let n = self.params.n; + + // Construire la matrice augmentée de taille m (n + m) + let mut aug = DenseMatrixGF2::zeros(m, n + m); + for r in 0..m { + for c in 0..n { + aug.set(r, c, dense[r][c]); + } + aug.set(r, n + r, 1); // partie identité + } + + // Gauss-Jordan sur les n premières colonnes + let (col_perm, rank) = aug.gauss_jordan_gf2(); + if rank < m { + return Err(LdpcError::SingularMatrix); + } + + // Extraire P (colonnes n..n+m de aug) -> partie redondante de G + // G = [I_k | P^T] après transposition correcte + // TODO: extraire proprement P et construire G = [I_k | P] + let g = DenseMatrixGF2::identity(self.params.k); // placeholder + + let col_perm_inv = { + let mut inv = vec![0usize; col_perm.len()]; + for (i, &p) in col_perm.iter().enumerate() { + inv[p] = i; + } + inv + }; + + self.systematic_form = Some(SystematicForm { + g, + col_perm, + col_perm_inv, + }); + Ok(()) + } + + pub fn rate(&self) -> f64 { + self.params.rate() + } + pub fn n(&self) -> usize { + self.params.n + } + pub fn k(&self) -> usize { + self.params.k + } + pub fn m(&self) -> usize { + self.params.m() + } + pub fn girth(&self) -> usize { + self.graph.girth() + } + + pub fn is_codeword(&self, c: &[Gf2]) -> bool { + self.h.multiply_vec(c).iter().all(|&s| s == 0) + } +} + +// RNG + +fn build_rng(seed: Option) -> impl rand::Rng { + use rand::SeedableRng; + rand::rngs::StdRng::seed_from_u64(seed.unwrap_or_else(rand::random)) +} + +// Gallager +// H divisée en wc sous-matrices de taille (m/wc) n +// H1 est une matrice régulière (chaque ligne contient exactement wr uns) +// H2..Hwc sont des permutations aléatoires de colonnes de H1 + +fn generate_gallager(params: &LdpcParams, rng: &mut impl rand::Rng) -> Result { + let CodeTopology::Regular { wc, wr } = params.topology else { + return Err(LdpcError::InvalidParameters( + "Gallager nécessite un code régulier".into(), + )); + }; + let n = params.n; + let m = params.m(); + if m % wc != 0 { + return Err(LdpcError::InvalidParameters( + "m doit être divisible par wc".into(), + )); + } + let rows_per_block = m / wc; + let mut ones: Vec<(usize, usize)> = Vec::new(); + + // H1 : blocs réguliers de wr uns par ligne + for r in 0..rows_per_block { + for j in 0..wr { + ones.push((r, r * wr + j)); + } + } + // H2..Hwc : permutations aléatoires de colonnes + use rand::seq::SliceRandom; + for block in 1..wc { + let mut perm: Vec = (0..n).collect(); + perm.shuffle(rng); + for r in 0..rows_per_block { + for j in 0..wr { + ones.push((block * rows_per_block + r, perm[r * wr + j])); + } + } + } + Ok(SparseMatrixGF2::from_positions(m, n, ones)) +} + +// MacKay-Neal +// Ajoute les colonnes une à une avec poids wc +// Rejette toute colonne créant un cycle-4 (deux colonnes ne peuvent +// partager qu'un seul 1 commun). Relance si max_attempts est dépassé + +fn generate_mackay_neal( + params: &LdpcParams, + max_attempts: usize, + rng: &mut impl rand::Rng, +) -> Result { + let CodeTopology::Regular { wc, .. } = params.topology else { + return Err(LdpcError::InvalidParameters( + "MacKayNeal nécessite un code régulier".into(), + )); + }; + let n = params.n; + let m = params.m(); + let mut ones: Vec<(usize, usize)> = Vec::new(); + + use rand::seq::SliceRandom; + for col in 0..n { + let mut placed = false; + for _attempt in 0..max_attempts { + // Tirer wc lignes sans remise + let mut rows: Vec = (0..m).collect(); + rows.shuffle(rng); + let candidate: Vec = rows[..wc].to_vec(); + + // Vérifier cycle-4 : cette colonne ne partage pas ≥2 lignes avec une colonne existante + let mut ok = true; + let mut c2 = 0; + while c2 < col && ok { + let existing: Vec = ones + .iter() + .filter(|&&(_, c)| c == c2) + .map(|&(r, _)| r) + .collect(); + let shared = candidate.iter().filter(|r| existing.contains(r)).count(); + if shared >= 2 { + ok = false; + } + c2 += 1; + } + if ok { + for r in candidate { + ones.push((r, col)); + } + placed = true; + break; + } + } + if !placed { + return Err(LdpcError::GenerationFailed { + attempts: max_attempts, + }); + } + } + Ok(SparseMatrixGF2::from_positions(m, n, ones)) +} + +// PEG +// Progressive Edge Growth +// Pour chaque arête à ajouter, choisit le check-node qui maximise +// le girth local du graphe courant (BFS depuis le var-node courant) + +fn generate_peg(params: &LdpcParams) -> Result { + let CodeTopology::Regular { wc, .. } = params.topology else { + return Err(LdpcError::InvalidParameters( + "PEG nécessite un code régulier".into(), + )); + }; + let n = params.n; + let m = params.m(); + let mut var_to_chk: Vec> = vec![vec![]; n]; + let mut chk_to_var: Vec> = vec![vec![]; m]; + let mut ones: Vec<(usize, usize)> = Vec::new(); + + for v in 0..n { + for _edge in 0..wc { + // BFS depuis v dans le graphe courant pour trouver le check-node + // non-voisin de v le plus éloigné (maximise le girth local) + let best_chk = peg_find_best_check(v, &var_to_chk, &chk_to_var, m); + var_to_chk[v].push(best_chk); + chk_to_var[best_chk].push(v); + ones.push((best_chk, v)); + } + } + Ok(SparseMatrixGF2::from_positions(m, n, ones)) +} + +// BFS depuis le var-node v, retourne le check-node non-voisin de v +// qui est le plus éloigné (ou le moins chargé en cas d'égalité) +fn peg_find_best_check( + v: usize, + var_to_chk: &[Vec], + chk_to_var: &[Vec], + n_chk: usize, +) -> usize { + use std::collections::VecDeque; + let current_neighbors: &[usize] = &var_to_chk[v]; + let mut dist_chk = vec![usize::MAX; n_chk]; + let mut dist_var = vec![usize::MAX; var_to_chk.len()]; + dist_var[v] = 0; + let mut queue: VecDeque<(bool, usize, usize)> = VecDeque::new(); + queue.push_back((true, v, 0)); + let mut max_dist = 0; + let mut reachable_non_neighbors: Vec<(usize, usize)> = Vec::new(); // (dist, chk) + + while let Some((is_var, node, dist)) = queue.pop_front() { + if is_var { + for &c in &var_to_chk[node] { + if dist_chk[c] == usize::MAX { + dist_chk[c] = dist + 1; + max_dist = max_dist.max(dist + 1); + queue.push_back((false, c, dist + 1)); + } + } + } else { + for &u in &chk_to_var[node] { + if dist_var[u] == usize::MAX { + dist_var[u] = dist + 1; + queue.push_back((true, u, dist + 1)); + } + } + } + } + + // Parmi les check-nodes non-voisins de v, choisir le plus éloigné + // (ou le moins chargé en cas d'égalité) + let mut best = (0usize, 0usize, usize::MAX); // (dist, chk_id, charge) + for c in 0..n_chk { + if current_neighbors.contains(&c) { + continue; + } + let d = dist_chk[c]; + let load = chk_to_var[c].len(); + if d > best.0 || (d == best.0 && load < best.2) { + best = (d, c, load); + } + } + // Fallback : check le moins chargé si aucun atteignable + if best.0 == 0 { + (0..n_chk).min_by_key(|&c| chk_to_var[c].len()).unwrap_or(0) + } else { + best.1 + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// +// #[test] +// fn test_gallager_regular_degrees() { +// let params = LdpcParams { +// n: 12, +// k: 6, +// topology: CodeTopology::Regular { wc: 2, wr: 4 }, +// generation: GenerationMethod::Gallager, +// seed: Some(42), +// }; +// let code = LdpcCode::new(params).unwrap(); +// for col in 0..code.n() { +// assert_eq!(code.h.col_weight(col), 2); +// } +// } +// +// #[test] +// fn test_mackay_neal_no_4cycle() { +// let params = LdpcParams { +// n: 20, +// k: 10, +// topology: CodeTopology::Regular { wc: 3, wr: 6 }, +// generation: GenerationMethod::MacKayNeal { max_attempts: 1000 }, +// seed: Some(0), +// }; +// let code = LdpcCode::new(params).unwrap(); +// assert!(!code.graph.has_4_cycle()); +// } +// +// #[test] +// fn test_peg_girth_at_least_6() { +// let params = LdpcParams { +// n: 30, +// k: 15, +// topology: CodeTopology::Regular { wc: 3, wr: 6 }, +// generation: GenerationMethod::Peg, +// seed: None, +// }; +// let code = LdpcCode::new(params).unwrap(); +// assert!(code.girth() >= 6); +// } +// } diff --git a/srctmp/decoder.rs b/srctmp/decoder.rs new file mode 100644 index 0000000..14370e1 --- /dev/null +++ b/srctmp/decoder.rs @@ -0,0 +1,446 @@ +use crate::code::LdpcCode; +use crate::graph::TannerGraph; +use crate::matrix::SparseMatrixGF2; +use crate::{BitVec, Gf2, Llr}; + +// Résultat + +#[derive(Debug, Clone)] +pub enum DecoderResult { + Converged(BitVec), + MaxIterationsReached(BitVec), + Failure, +} + +impl DecoderResult { + pub fn codeword(&self) -> Option<&BitVec> { + match self { + DecoderResult::Converged(c) | DecoderResult::MaxIterationsReached(c) => Some(c), + DecoderResult::Failure => None, + } + } + pub fn is_success(&self) -> bool { + matches!(self, DecoderResult::Converged(_)) + } +} + +// Configuration + +#[derive(Debug, Clone)] +pub struct DecoderConfig { + pub max_iterations: usize, + pub early_stopping: bool, +} + +impl Default for DecoderConfig { + fn default() -> Self { + Self { + max_iterations: 50, + early_stopping: true, + } + } +} + +// Trait Decoder + +pub trait Decoder: Send + Sync { + fn decode(&self, channel_llr: &[Llr]) -> DecoderResult; + + fn decode_hard(&self, received: &[Gf2]) -> DecoderResult { + let llr: Vec = received + .iter() + .map(|&b| if b == 0 { 1.0 } else { -1.0 }) + .collect(); + self.decode(&llr) + } +} + +// Primitives GF(2) et LLR + +#[inline] +pub fn hard_decision(llr: Llr) -> Gf2 { + if llr >= 0.0 { + 0 + } else { + 1 + } +} + +pub fn compute_syndrome(h: &SparseMatrixGF2, c: &[Gf2]) -> Vec { + h.multiply_vec(c) +} + +// phi(x) = -ln(tanh(|x|/2)) involution de Sum-Product +// phi(phi(x)) = x +#[inline] +fn phi(x: Llr) -> Llr { + let ax = x.abs().max(1e-10); + -((ax / 2.0).tanh().ln()) +} + +// Mises à jour des nœuds + +// Mise à jour Sum-Product du nœud de contrôle +// R_{c→v} = φ(∑_{v'≠v} φ(|L_{v'→c}|)) × sign(∏_{v'≠v} sign(L_{v'→c})) +fn check_node_update_sp(incoming: &[Llr], out: &mut [Llr]) { + let phi_sum: Llr = incoming.iter().map(|&l| phi(l.abs())).sum(); + let sign_prod: Llr = incoming.iter().map(|&l| l.signum()).product(); + for (j, (&l, r)) in incoming.iter().zip(out.iter_mut()).enumerate() { + let phi_excl = phi_sum - phi(l.abs()); + let sign_excl = sign_prod * l.signum(); + *r = sign_excl * phi(phi_excl); + } +} + +// Mise à jour Min-Sum avec facteur de normalisation α +// Approxime φ(∑ φ(|L|)) ≈ min(|L|). +// alpha in [0.75, 0.875] compense le biais de Min-Sum brut +fn check_node_update_ms(incoming: &[Llr], out: &mut [Llr], alpha: Llr) { + let sign_prod: Llr = incoming.iter().map(|&l| l.signum()).product(); + // Précalcul des deux plus petites valeurs absolues + let mut min1 = Llr::INFINITY; + let mut min2 = Llr::INFINITY; + let mut min1_idx = 0; + for (j, &l) in incoming.iter().enumerate() { + let al = l.abs(); + if al < min1 { + min2 = min1; + min1 = al; + min1_idx = j; + } else if al < min2 { + min2 = al; + } + } + for (j, (&l, r)) in incoming.iter().zip(out.iter_mut()).enumerate() { + let min_excl = if j == min1_idx { min2 } else { min1 }; + let sign_excl = sign_prod * l.signum(); + *r = alpha * sign_excl * min_excl; + } +} + +// Mise à jour du nœud de variable +// L_{v→c} = L_channel(v) + ∑_{c'≠c} R_{c'→v} +fn variable_node_update(ch_llr: Llr, incoming_c2v: &[Llr], out: &mut [Llr]) { + let total: Llr = ch_llr + incoming_c2v.iter().sum::(); + for ((&r, o)) in incoming_c2v.iter().zip(out.iter_mut()) { + *o = total - r; + } +} + +#[inline] +fn posterior_llr(ch_llr: Llr, c2v_msgs: &[Llr]) -> Llr { + ch_llr + c2v_msgs.iter().sum::() +} + +// Messages internes +// Indexés par (check_id, position_dans_liste_voisins) accès O(1) + +struct Messages { + v2c: Vec>, // v2c[c][j] : var_neighbor(c)[j] -> check c + c2v: Vec>, // c2v[c][j] : check c -> var_neighbor(c)[j] +} + +impl Messages { + fn new(graph: &TannerGraph) -> Self { + let v2c = (0..graph.n_chk) + .map(|c| vec![0.0; graph.chk_degree(c)]) + .collect(); + let c2v = (0..graph.n_chk) + .map(|c| vec![0.0; graph.chk_degree(c)]) + .collect(); + Self { v2c, c2v } + } +} + +// Table de correspondance, pour chaque (var, check), index dans la liste de voisins +// Précalculée une fois à la construction du décodeur +struct EdgeIndex { + // var_pos_in_chk[c][j] : position de var_neighbor(c)[j] dans var_to_chk[v] + var_pos_in_chk: Vec>, + // chk_pos_in_var[v][i] : position de var_neighbor(v)[i] dans chk_to_var[c] + chk_pos_in_var: Vec>, +} + +impl EdgeIndex { + fn build(graph: &TannerGraph) -> Self { + let var_pos_in_chk = (0..graph.n_chk) + .map(|c| { + graph + .chk_neighbors(c) + .iter() + .map(|&v| graph.var_neighbors(v).iter().position(|&x| x == c).unwrap()) + .collect() + }) + .collect(); + let chk_pos_in_var = (0..graph.n_var) + .map(|v| { + graph + .var_neighbors(v) + .iter() + .map(|&c| graph.chk_neighbors(c).iter().position(|&x| x == v).unwrap()) + .collect() + }) + .collect(); + Self { + var_pos_in_chk, + chk_pos_in_var, + } + } +} + +// Bit-Flipping + +pub struct BitFlippingDecoder { + graph: TannerGraph, + h: SparseMatrixGF2, + config: DecoderConfig, +} + +impl BitFlippingDecoder { + pub fn new(code: &LdpcCode, config: DecoderConfig) -> Self { + Self { + graph: code.graph.clone(), + h: code.h.clone(), + config, + } + } +} + +impl Decoder for BitFlippingDecoder { + fn decode(&self, channel_llr: &[Llr]) -> DecoderResult { + let mut bits: Vec = channel_llr.iter().map(|&l| hard_decision(l)).collect(); + + for _iter in 0..self.config.max_iterations { + let syndrome = compute_syndrome(&self.h, &bits); + if self.config.early_stopping && syndrome.iter().all(|&s| s == 0) { + return DecoderResult::Converged(bits); + } + let mut unsatisfied = vec![0usize; self.graph.n_var]; + for c in 0..self.graph.n_chk { + if syndrome[c] == 1 { + for &v in self.graph.chk_neighbors(c) { + unsatisfied[v] += 1; + } + } + } + let mut flipped = false; + for v in 0..self.graph.n_var { + if unsatisfied[v] > self.graph.var_degree(v) / 2 { + bits[v] ^= 1; + flipped = true; + } + } + if !flipped { + break; + } + } + + let synd = compute_syndrome(&self.h, &bits); + if synd.iter().all(|&s| s == 0) { + DecoderResult::Converged(bits) + } else { + DecoderResult::MaxIterationsReached(bits) + } + } +} + +// Noyau BP partagé par SP et MinSum + +fn bp_decode( + graph: &TannerGraph, + h: &SparseMatrixGF2, + config: &DecoderConfig, + channel_llr: &[Llr], + edge_idx: &EdgeIndex, + check_update: F, +) -> DecoderResult +where + F: Fn(&[Llr], &mut [Llr]), +{ + let mut msgs = Messages::new(graph); + + // Initialisation : v2c[c][j] = L_channel(var_neighbor(c)[j]) + for c in 0..graph.n_chk { + for (j, &v) in graph.chk_neighbors(c).iter().enumerate() { + msgs.v2c[c][j] = channel_llr[v]; + } + } + + for _iter in 0..config.max_iterations { + // Mise à jour des check-nodes + for c in 0..graph.n_chk { + let v2c = msgs.v2c[c].clone(); + check_update(&v2c, &mut msgs.c2v[c]); + } + + // Mise à jour des var-nodes + for v in 0..graph.n_var { + let neighbors = graph.var_neighbors(v); + // Rassembler les c2v entrants sur ce var-node + let incoming: Vec = neighbors + .iter() + .enumerate() + .map(|(i, &c)| { + let j = edge_idx.chk_pos_in_var[v][i]; + msgs.c2v[c][j] + }) + .collect(); + let mut new_v2c = vec![0.0; neighbors.len()]; + variable_node_update(channel_llr[v], &incoming, &mut new_v2c); + for (i, &c) in neighbors.iter().enumerate() { + let j = edge_idx.chk_pos_in_var[v][i]; + msgs.v2c[c][j] = new_v2c[i]; + } + } + + // Hard décision + arrêt + if config.early_stopping { + let bits = make_decision(graph, &msgs, channel_llr, edge_idx); + if compute_syndrome(h, &bits).iter().all(|&s| s == 0) { + return DecoderResult::Converged(bits); + } + } + } + + let bits = make_decision(graph, &msgs, channel_llr, edge_idx); + let synd = compute_syndrome(h, &bits); + if synd.iter().all(|&s| s == 0) { + DecoderResult::Converged(bits) + } else { + DecoderResult::MaxIterationsReached(bits) + } +} + +fn make_decision( + graph: &TannerGraph, + msgs: &Messages, + channel_llr: &[Llr], + edge_idx: &EdgeIndex, +) -> Vec { + (0..graph.n_var) + .map(|v| { + let incoming: Vec = graph + .var_neighbors(v) + .iter() + .enumerate() + .map(|(i, &c)| { + let j = edge_idx.chk_pos_in_var[v][i]; + msgs.c2v[c][j] + }) + .collect(); + hard_decision(posterior_llr(channel_llr[v], &incoming)) + }) + .collect() +} + +// Sum-Product + +pub struct SumProductDecoder { + graph: TannerGraph, + h: SparseMatrixGF2, + config: DecoderConfig, + edge_idx: EdgeIndex, +} + +impl SumProductDecoder { + pub fn new(code: &LdpcCode, config: DecoderConfig) -> Self { + let edge_idx = EdgeIndex::build(&code.graph); + Self { + graph: code.graph.clone(), + h: code.h.clone(), + config, + edge_idx, + } + } +} + +impl Decoder for SumProductDecoder { + fn decode(&self, channel_llr: &[Llr]) -> DecoderResult { + bp_decode( + &self.graph, + &self.h, + &self.config, + channel_llr, + &self.edge_idx, + |incoming, out| check_node_update_sp(incoming, out), + ) + } +} + +// Min-Sum + +pub struct MinSumDecoder { + graph: TannerGraph, + h: SparseMatrixGF2, + config: DecoderConfig, + scaling_factor: Llr, + edge_idx: EdgeIndex, +} + +impl MinSumDecoder { + pub fn new(code: &LdpcCode, config: DecoderConfig, scaling_factor: Llr) -> Self { + let edge_idx = EdgeIndex::build(&code.graph); + Self { + graph: code.graph.clone(), + h: code.h.clone(), + config, + scaling_factor, + edge_idx, + } + } +} + +impl Decoder for MinSumDecoder { + fn decode(&self, channel_llr: &[Llr]) -> DecoderResult { + let alpha = self.scaling_factor; + bp_decode( + &self.graph, + &self.h, + &self.config, + channel_llr, + &self.edge_idx, + move |incoming, out| check_node_update_ms(incoming, out, alpha), + ) + } +} + +// Factory + +#[derive(Debug, Clone)] +pub enum DecoderMethod { + BitFlipping, + SumProduct, + MinSum { scaling_factor: Llr }, +} + +pub fn build_decoder( + code: &LdpcCode, + method: DecoderMethod, + config: DecoderConfig, +) -> Box { + match method { + DecoderMethod::BitFlipping => Box::new(BitFlippingDecoder::new(code, config)), + DecoderMethod::SumProduct => Box::new(SumProductDecoder::new(code, config)), + DecoderMethod::MinSum { scaling_factor } => { + Box::new(MinSumDecoder::new(code, config, scaling_factor)) + } + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// +// #[test] +// fn test_phi_is_involution() { +// let x = 2.5_f64; +// assert!((phi(phi(x)) - x).abs() < 1e-9); +// } +// +// #[test] +// fn test_hard_decision() { +// assert_eq!(hard_decision(0.1), 0); +// assert_eq!(hard_decision(-0.1), 1); +// assert_eq!(hard_decision(0.0), 0); +// } +// } diff --git a/srctmp/encoder.rs b/srctmp/encoder.rs new file mode 100644 index 0000000..9ea240e --- /dev/null +++ b/srctmp/encoder.rs @@ -0,0 +1,202 @@ +use crate::code::{LdpcCode, SystematicForm}; +use crate::matrix::{DenseMatrixGF2, SparseMatrixGF2}; +use crate::{BitVec, Gf2, LdpcError, Result}; + +// Trait Encoder + +pub trait Encoder: Send + Sync { + fn encode(&self, message: &[Gf2]) -> Result; + fn message_len(&self) -> usize; + fn codeword_len(&self) -> usize; + + fn check_input(&self, msg: &[Gf2]) -> Result<()> { + if msg.len() != self.message_len() { + return Err(LdpcError::DimensionMismatch { + expected: self.message_len(), + got: msg.len(), + }); + } + Ok(()) + } +} + +// Méthode d'encodage + +#[derive(Debug, Clone)] +pub enum EncodingMethod { + // Via G = [I_k | P]. Complexité O(n^2). Point de départ simple + Systematic, + // Richardson-Urbanke. Complexité O(n + g^2), g = gap << n + RichardsonUrbanke, +} + +// Encodeur systématique +// c = m * G = [m | m*P] +// Les bits de parité p = m * P sont calculés par multiplication dense +// Le mot de code est ensuite réordonné selon la permutation inverse + +pub struct SystematicEncoder { + k: usize, + n: usize, + g: DenseMatrixGF2, + perm_inv: Vec, +} + +impl SystematicEncoder { + pub fn new(code: &mut LdpcCode) -> Result { + code.compute_systematic_form()?; + let sf = code.systematic_form.as_ref().unwrap(); + Ok(Self { + k: code.k(), + n: code.n(), + g: sf.g.clone(), + perm_inv: sf.col_perm_inv.clone(), + }) + } +} + +impl Encoder for SystematicEncoder { + fn encode(&self, message: &[Gf2]) -> Result { + self.check_input(message)?; + let c_perm = self.g.multiply_vec(message); + // Réordonner les bits selon la permutation inverse + let mut c = vec![0u8; self.n]; + for (i, &ci) in c_perm.iter().enumerate() { + c[self.perm_inv[i]] = ci; + } + Ok(c) + } + + fn message_len(&self) -> usize { + self.k + } + fn codeword_len(&self) -> usize { + self.n + } +} + +// Encodeur Richardson-Urbanke +// H est réarrangée par permutations en : +// +// ┌ ┐ +// │ A B T │ +// H = │ │ T = triangulaire inférieure, φ = -ET^{-1}B + D +// │ C D E │ +// └ ┘ +// +// Encodage en deux étapes : +// 1. p₂ = φ^{-1} * (-ET^{-1}*As - Cs) [dense g×g, g = gap] +// 2. p₁ = T^{-1} * (As + Bp₂) [substitution arrière] +// Complexité totale : O(n + g²) + +pub struct RichardsonUrbankeEncoder { + k: usize, + n: usize, + a: SparseMatrixGF2, + b: SparseMatrixGF2, + c: SparseMatrixGF2, + d: SparseMatrixGF2, + e: SparseMatrixGF2, + // T^{-1} : résolution par substitution arrière (T triangulaire inférieure) + t_inv: DenseMatrixGF2, + // φ^{-1} : petit (g×g), calculé une fois offline + phi_inv: DenseMatrixGF2, + col_perm_inv: Vec, +} + +impl RichardsonUrbankeEncoder { + pub fn new(code: &LdpcCode) -> Result { + // TODO: implémenter le pré-traitement de H par permutations de lignes/colonnes + // pour identifier les blocs A, B, C, D, E, T et calculer phi^{-1}. + todo!("Pré-traitement Richardson-Urbanke") + } + + // Substitution arrière sur T triangulaire inférieure : résout T*x = b en GF(2) + fn backward_substitution(t_inv: &DenseMatrixGF2, b: &[Gf2]) -> Vec { + let n = b.len(); + let mut x = vec![0u8; n]; + for i in 0..n { + let mut s = b[i]; + for j in 0..i { + s ^= t_inv.get(i, j) & x[j]; + } + x[i] = s; + } + x + } +} + +impl Encoder for RichardsonUrbankeEncoder { + fn encode(&self, message: &[Gf2]) -> Result { + self.check_input(message)?; + // Étape 1 : As + let a_s = self.a.multiply_vec(message); + // Étape 2 : p₂ = φ^{-1} * (E*T^{-1}*As + Cs) + let t_inv_as = Self::backward_substitution(&self.t_inv, &a_s); + let e_t_inv_as = self.e.multiply_vec(&t_inv_as); + let c_s = self.c.multiply_vec(message); + let rhs: Vec = e_t_inv_as + .iter() + .zip(c_s.iter()) + .map(|(&a, &b)| a ^ b) + .collect(); + let p2 = self.phi_inv.multiply_vec(&rhs); + // Étape 3 : p₁ = T^{-1} * (As + Bp₂) + let b_p2 = self.b.multiply_vec(&p2); + let sum: Vec = a_s.iter().zip(b_p2.iter()).map(|(&a, &b)| a ^ b).collect(); + let p1 = Self::backward_substitution(&self.t_inv, &sum); + // Assemblage et dépermutation + let mut c_perm: Vec = message + .iter() + .chain(p1.iter()) + .chain(p2.iter()) + .cloned() + .collect(); + let mut c = vec![0u8; self.n]; + for (i, &ci) in c_perm.iter().enumerate() { + c[self.col_perm_inv[i]] = ci; + } + Ok(c) + } + + fn message_len(&self) -> usize { + self.k + } + fn codeword_len(&self) -> usize { + self.n + } +} + +// Factory + +pub fn build_encoder(code: &mut LdpcCode, method: EncodingMethod) -> Result> { + match method { + EncodingMethod::Systematic => Ok(Box::new(SystematicEncoder::new(code)?)), + EncodingMethod::RichardsonUrbanke => Ok(Box::new(RichardsonUrbankeEncoder::new(code)?)), + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// use crate::code::{CodeTopology, GenerationMethod, LdpcCode, LdpcParams}; +// +// fn make_code() -> LdpcCode { +// LdpcCode::new(LdpcParams { +// n: 12, +// k: 6, +// topology: CodeTopology::Regular { wc: 2, wr: 4 }, +// generation: GenerationMethod::Gallager, +// seed: Some(42), +// }) +// .unwrap() +// } +// +// #[test] +// fn test_encoder_wrong_input_length_errors() { +// let mut code = make_code(); +// let enc = SystematicEncoder::new(&mut code).unwrap(); +// let result = enc.encode(&vec![0u8; 5]); // devrait être 6 +// assert!(result.is_err()); +// } +// } diff --git a/srctmp/graph.rs b/srctmp/graph.rs new file mode 100644 index 0000000..6738cfd --- /dev/null +++ b/srctmp/graph.rs @@ -0,0 +1,171 @@ +use crate::matrix::SparseMatrixGF2; +use std::collections::VecDeque; + +// Graphe de Tanner + +#[derive(Debug, Clone)] +pub struct TannerGraph { + pub n_var: usize, + pub n_chk: usize, + var_to_chk: Vec>, + chk_to_var: Vec>, +} + +impl TannerGraph { + pub fn from_matrix(h: &SparseMatrixGF2) -> Self { + let n_var = h.cols; + let n_chk = h.rows; + let chk_to_var: Vec> = (0..n_chk).map(|c| h.row_neighbors(c).to_vec()).collect(); + let mut var_to_chk = vec![vec![]; n_var]; + for c in 0..n_chk { + for &v in &chk_to_var[c] { + var_to_chk[v].push(c); + } + } + Self { + n_var, + n_chk, + var_to_chk, + chk_to_var, + } + } + + pub fn var_neighbors(&self, v: usize) -> &[usize] { + &self.var_to_chk[v] + } + pub fn chk_neighbors(&self, c: usize) -> &[usize] { + &self.chk_to_var[c] + } + pub fn var_degree(&self, v: usize) -> usize { + self.var_to_chk[v].len() + } + pub fn chk_degree(&self, c: usize) -> usize { + self.chk_to_var[c].len() + } + + // Calcule le girth par BFS depuis chaque noeud de variable + // O(n * (n + m)) + pub fn girth(&self) -> usize { + let mut min_girth = usize::MAX; + for start in 0..self.n_var { + if let Some(g) = self.bfs_girth_from_var(start) { + min_girth = min_girth.min(g); + if min_girth == 4 { + return 4; + } // minimum possible + } + } + min_girth + } + + // Détection rapide de cycles-4, deux var-nodes partagent >= check-nodes + pub fn has_4_cycle(&self) -> bool { + for v1 in 0..self.n_var { + for v2 in (v1 + 1)..self.n_var { + let common = self.var_to_chk[v1] + .iter() + .filter(|c| self.var_to_chk[v2].contains(c)) + .count(); + if common >= 2 { + return true; + } + } + } + false + } + + // Girth local depuis un noeud de variable v (pour PEG). + pub fn local_girth_from_var(&self, v: usize) -> usize { + self.bfs_girth_from_var(v).unwrap_or(usize::MAX) + } + + // BFS depuis le noeud de variable start, + // retourne la longueur du court cycle passant par ce noeud (None si aucun cycle) + fn bfs_girth_from_var(&self, start: usize) -> Option { + // dist_var[v] = distance depuis start jusqu'au var-node v + // dist_chk[c] = distance depuis start jusqu'au check-node c + let mut dist_var = vec![usize::MAX; self.n_var]; + let mut dist_chk = vec![usize::MAX; self.n_chk]; + dist_var[start] = 0; + + // File : (is_var: bool, index, distance) + let mut queue: VecDeque<(bool, usize, usize)> = VecDeque::new(); + queue.push_back((true, start, 0)); + let mut shortest = None; + + while let Some((is_var, node, dist)) = queue.pop_front() { + if is_var { + for &c in self.var_neighbors(node) { + if dist_chk[c] == usize::MAX { + dist_chk[c] = dist + 1; + queue.push_back((false, c, dist + 1)); + } else { + // Cycle trouvé + let cycle_len = dist + 1 + dist_chk[c]; + shortest = Some(shortest.map_or(cycle_len, |s: usize| s.min(cycle_len))); + } + } + } else { + for &v in self.chk_neighbors(node) { + if v == start && dist + 1 >= 2 { + let cycle_len = dist + 1; + shortest = Some(shortest.map_or(cycle_len, |s: usize| s.min(cycle_len))); + continue; + } + if dist_var[v] == usize::MAX { + dist_var[v] = dist + 1; + queue.push_back((true, v, dist + 1)); + } + } + } + } + shortest + } + + pub fn var_degree_distribution(&self) -> Vec { + let max_deg = self.var_to_chk.iter().map(|v| v.len()).max().unwrap_or(0); + let mut counts = vec![0usize; max_deg + 1]; + for v in 0..self.n_var { + counts[self.var_degree(v)] += 1; + } + counts + .iter() + .map(|&c| c as f64 / self.n_var as f64) + .collect() + } + + pub fn is_regular(&self) -> bool { + let d0 = self.var_degree(0); + let c0 = self.chk_degree(0); + self.var_to_chk.iter().all(|v| v.len() == d0) + && self.chk_to_var.iter().all(|c| c.len() == c0) + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// use crate::matrix::SparseMatrixGF2; +// +// fn simple_h() -> SparseMatrixGF2 { +// SparseMatrixGF2::from_dense(&vec![ +// vec![1u8, 1, 0, 1, 0], +// vec![0, 1, 1, 0, 1], +// vec![1, 0, 1, 0, 1], +// ]) +// } +// +// #[test] +// fn test_construction_from_matrix() { +// let h = simple_h(); +// let g = TannerGraph::from_matrix(&h); +// assert_eq!(g.n_var, 5); +// assert_eq!(g.n_chk, 3); +// } +// +// #[test] +// fn test_var_degree() { +// let g = TannerGraph::from_matrix(&simple_h()); +// assert_eq!(g.var_degree(0), 2); +// } +// } diff --git a/srctmp/lib.rs b/srctmp/lib.rs new file mode 100644 index 0000000..6dca682 --- /dev/null +++ b/srctmp/lib.rs @@ -0,0 +1,38 @@ +pub mod channel; +pub mod code; +pub mod decoder; +pub mod encoder; +pub mod graph; +pub mod matrix; + +pub type Gf2 = u8; +pub type Llr = f64; +pub type BitVec = Vec; + +#[derive(Debug, thiserror::Error)] +pub enum LdpcError { + #[error("Paramètres invalides : {0}")] + InvalidParameters(String), + + #[error("Matrice singulière : impossible d'inverser")] + SingularMatrix, + + #[error("Génération échouée après {attempts} tentatives")] + GenerationFailed { attempts: usize }, + + #[error("Dimension incorrecte : attendu {expected}, reçu {got}")] + DimensionMismatch { expected: usize, got: usize }, + + #[error("Le vecteur fourni n'est pas un mot de code valide")] + InvalidCodeword, + + #[error("Paramètre hors plage : {0}")] + OutOfRange(String), +} + +pub type Result = std::result::Result; + +pub use channel::Channel; +pub use code::{CodeTopology, GenerationMethod, LdpcCode, LdpcParams}; +pub use decoder::{Decoder, DecoderConfig, DecoderMethod, DecoderResult}; +pub use encoder::{Encoder, EncodingMethod}; diff --git a/srctmp/main.rs b/srctmp/main.rs new file mode 100644 index 0000000..490b7f9 --- /dev/null +++ b/srctmp/main.rs @@ -0,0 +1,60 @@ +use ldpc::{ + channel::{AwgnChannel, Channel}, + code::{CodeTopology, GenerationMethod, LdpcCode, LdpcParams}, + decoder::{build_decoder, DecoderConfig, DecoderMethod}, + encoder::{build_encoder, EncodingMethod}, +}; +use rand::{Rng, SeedableRng}; + +fn main() -> ldpc::Result<()> { + // ── 1. Générer le code ──────────────────────────────────────────────────── + let params = LdpcParams { + n: 100, + k: 50, + topology: CodeTopology::Regular { wc: 3, wr: 6 }, + generation: GenerationMethod::MacKayNeal { max_attempts: 500 }, + seed: Some(42), + }; + let mut code = LdpcCode::new(params)?; + println!( + "Code LDPC : n={}, k={}, taux={:.2}, girth={}", + code.n(), + code.k(), + code.rate(), + code.girth() + ); + + // ── 2. Encodeur ─────────────────────────────────────────────────────────── + let encoder = build_encoder(&mut code, EncodingMethod::Systematic)?; + + // ── 3. Canal ────────────────────────────────────────────────────────────── + let channel = AwgnChannel::new(2.0, code.rate())?; // 2 dB Eb/N0 + println!("Capacité AWGN ≈ {:.4} bits/utilisation", channel.capacity()); + + // ── 4. Décodeur ─────────────────────────────────────────────────────────── + let decoder = build_decoder(&code, DecoderMethod::SumProduct, DecoderConfig::default()); + + // ── 5. Simulation ───────────────────────────────────────────────────────── + let mut rng = rand::rngs::StdRng::seed_from_u64(123); + let n_trials = 100; + let mut errors = 0usize; + + for _ in 0..n_trials { + // Message aléatoire + let message: Vec = (0..code.k()).map(|_| rng.gen::() & 1).collect(); + // Encodage + let codeword = encoder.encode(&message)?; + // Transmission AWGN + let received_llr = channel.transmit(&codeword, &mut rng); + // Décodage + let result = decoder.decode(&received_llr); + // Vérification + if !result.is_success() { + errors += 1; + } + } + + let ber = errors as f64 / n_trials as f64; + println!("FER sur {} essais à 2dB : {:.2}%", n_trials, ber * 100.0); + Ok(()) +} diff --git a/srctmp/matrix.rs b/srctmp/matrix.rs new file mode 100644 index 0000000..9dd1f7a --- /dev/null +++ b/srctmp/matrix.rs @@ -0,0 +1,295 @@ +use crate::{Gf2, Result}; + +// Matrice creuse GF(2) — format CSR + CSC dual + +#[derive(Debug, Clone)] +pub struct SparseMatrixGF2 { + pub rows: usize, + pub cols: usize, + // CSR — accès ligne i : col_idx[ row_ptr[i] .. row_ptr[i+1] ] + row_ptr: Vec, + col_idx: Vec, + // CSC — accès col j : row_idx[ col_ptr[j] .. col_ptr[j+1] ] + col_ptr: Vec, + row_idx: Vec, +} + +impl SparseMatrixGF2 { + pub fn zeros(rows: usize, cols: usize) -> Self { + Self { + rows, + cols, + row_ptr: vec![0; rows + 1], + col_idx: vec![], + col_ptr: vec![0; cols + 1], + row_idx: vec![], + } + } + + // Depuis une liste de (row, col) indiquant les positions des 1s. + // Trie les entrées et construit CSR + CSC en un seul passage. + pub fn from_positions(rows: usize, cols: usize, mut ones: Vec<(usize, usize)>) -> Self { + // Construction CSR + ones.sort_unstable(); + let mut row_ptr = vec![0usize; rows + 1]; + let mut col_idx = Vec::with_capacity(ones.len()); + for &(r, c) in &ones { + row_ptr[r + 1] += 1; + col_idx.push(c); + } + for i in 0..rows { + row_ptr[i + 1] += row_ptr[i]; + } + // Construction CSC + let mut col_sorted = ones.clone(); + col_sorted.sort_unstable_by_key(|&(r, c)| (c, r)); + let mut col_ptr = vec![0usize; cols + 1]; + let mut row_idx = Vec::with_capacity(ones.len()); + for &(r, c) in &col_sorted { + col_ptr[c + 1] += 1; + row_idx.push(r); + } + for j in 0..cols { + col_ptr[j + 1] += col_ptr[j]; + } + Self { + rows, + cols, + row_ptr, + col_idx, + col_ptr, + row_idx, + } + } + + pub fn from_dense(dense: &[Vec]) -> Self { + let rows = dense.len(); + let cols = if rows > 0 { dense[0].len() } else { 0 }; + let ones: Vec<(usize, usize)> = dense + .iter() + .enumerate() + .flat_map(|(r, row)| { + row.iter() + .enumerate() + .filter(|(_, &v)| v == 1) + .map(move |(c, _)| (r, c)) + }) + .collect(); + Self::from_positions(rows, cols, ones) + } + + pub fn get(&self, row: usize, col: usize) -> Gf2 { + let slice = self.row_neighbors(row); + if slice.binary_search(&col).is_ok() { + 1 + } else { + 0 + } + } + + // Indices des colonnes où la ligne `row` vaut 1 (voisins check→var) + pub fn row_neighbors(&self, row: usize) -> &[usize] { + &self.col_idx[self.row_ptr[row]..self.row_ptr[row + 1]] + } + + // Indices des lignes où la colonne `col` vaut 1 (voisins var→check) + pub fn col_neighbors(&self, col: usize) -> &[usize] { + &self.row_idx[self.col_ptr[col]..self.col_ptr[col + 1]] + } + + pub fn row_weight(&self, row: usize) -> usize { + self.row_ptr[row + 1] - self.row_ptr[row] + } + + pub fn col_weight(&self, col: usize) -> usize { + self.col_ptr[col + 1] - self.col_ptr[col] + } + + pub fn nnz(&self) -> usize { + self.col_idx.len() + } + + pub fn density(&self) -> f64 { + self.nnz() as f64 / (self.rows * self.cols) as f64 + } + + // Produit H * x mod 2 (calcul du syndrome : s = H * c^T) + pub fn multiply_vec(&self, x: &[Gf2]) -> Vec { + (0..self.rows) + .map(|r| { + self.row_neighbors(r) + .iter() + .map(|&c| x[c]) + .fold(0u8, |acc, b| acc ^ b) + }) + .collect() + } + + pub fn transpose(&self) -> Self { + Self { + rows: self.cols, + cols: self.rows, + row_ptr: self.col_ptr.clone(), + col_idx: self.row_idx.clone(), + col_ptr: self.row_ptr.clone(), + row_idx: self.col_idx.clone(), + } + } + + // Vérifie si deux colonnes partagent >= 2 positions de 1 -> cycle-4 détecté + pub fn columns_share_two_ones(&self, c1: usize, c2: usize) -> bool { + let n1 = self.col_neighbors(c1); + let n2 = self.col_neighbors(c2); + let mut common = 0usize; + let (mut i, mut j) = (0, 0); + while i < n1.len() && j < n2.len() { + match n1[i].cmp(&n2[j]) { + std::cmp::Ordering::Less => i += 1, + std::cmp::Ordering::Greater => j += 1, + std::cmp::Ordering::Equal => { + common += 1; + if common >= 2 { + return true; + } + i += 1; + j += 1; + } + } + } + false + } + + pub fn to_dense(&self) -> Vec> { + let mut out = vec![vec![0u8; self.cols]; self.rows]; + for r in 0..self.rows { + for &c in self.row_neighbors(r) { + out[r][c] = 1; + } + } + out + } +} + +// Matrice dense GF(2) +// Utilisée pour la matrice génératrice G et les calculs de Gauss-Jordan +// G = [I | P], P est souvent dense + +#[derive(Debug, Clone)] +pub struct DenseMatrixGF2 { + pub rows: usize, + pub cols: usize, + data: Vec>, +} + +impl DenseMatrixGF2 { + pub fn zeros(rows: usize, cols: usize) -> Self { + Self { + rows, + cols, + data: vec![vec![0u8; cols]; rows], + } + } + + pub fn identity(n: usize) -> Self { + let mut m = Self::zeros(n, n); + for i in 0..n { + m.data[i][i] = 1; + } + m + } + + pub fn get(&self, row: usize, col: usize) -> Gf2 { + self.data[row][col] + } + pub fn set(&mut self, row: usize, col: usize, val: Gf2) { + self.data[row][col] = val; + } + + // Addition de deux lignes dans GF(2) + pub fn row_add(&mut self, dst: usize, src: usize) { + for j in 0..self.cols { + self.data[dst][j] ^= self.data[src][j]; + } + } + + pub fn row_swap(&mut self, r1: usize, r2: usize) { + self.data.swap(r1, r2); + } + + pub fn multiply_vec(&self, x: &[Gf2]) -> Vec { + self.data + .iter() + .map(|row| { + row.iter() + .zip(x.iter()) + .fold(0u8, |acc, (&a, &b)| acc ^ (a & b)) + }) + .collect() + } + + pub fn into_sparse(self) -> SparseMatrixGF2 { + SparseMatrixGF2::from_dense(&self.data) + } + + // Retourne la permutation de colonnes appliquée et le rang + // self est sous forme échelonnée réduite + pub fn gauss_jordan_gf2(&mut self) -> (Vec, usize) { + let mut perm: Vec = (0..self.cols).collect(); + let mut pivot_row = 0; + for col in 0..self.cols { + // Chercher un pivot dans la colonne courante + let pivot = (pivot_row..self.rows).find(|&r| self.data[r][col] == 1); + let Some(p) = pivot else { continue }; + self.row_swap(pivot_row, p); + // Éliminer dans toutes les autres lignes + for r in 0..self.rows { + if r != pivot_row && self.data[r][col] == 1 { + self.row_add(r, pivot_row); + } + } + pivot_row += 1; + if pivot_row == self.rows { + break; + } + } + (perm, pivot_row) + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// +// #[test] +// fn test_from_dense_roundtrip() { +// let dense = vec![vec![1u8, 0, 1, 0], vec![0, 1, 0, 1], vec![1, 1, 0, 0]]; +// let sparse = SparseMatrixGF2::from_dense(&dense); +// assert_eq!(sparse.to_dense(), dense); +// } +// +// #[test] +// fn test_multiply_vec_gf2() { +// let dense = vec![vec![1u8, 1, 0], vec![0, 1, 1]]; +// let h = SparseMatrixGF2::from_dense(&dense); +// let x = vec![1u8, 1, 1]; +// let s = h.multiply_vec(&x); +// // ligne 0 : 1^1^0 = 0, ligne 1 : 0^1^1 = 0 +// assert_eq!(s, vec![0u8, 0]); +// } +// +// #[test] +// fn test_transpose_double_is_identity() { +// let dense = vec![vec![1u8, 0, 1], vec![0, 1, 0]]; +// let h = SparseMatrixGF2::from_dense(&dense); +// assert_eq!(h.transpose().transpose().to_dense(), dense); +// } +// +// #[test] +// fn test_columns_share_two_ones() { +// // Colonnes 0 et 1 partagent les lignes 0 et 1 → cycle-4 +// let dense = vec![vec![1u8, 1, 0], vec![1, 1, 0], vec![0, 0, 1]]; +// let h = SparseMatrixGF2::from_dense(&dense); +// assert!(h.columns_share_two_ones(0, 1)); +// assert!(!h.columns_share_two_ones(0, 2)); +// } +// }