This commit is contained in:
2026-05-07 18:32:36 +02:00
parent 2fcd368e34
commit b7faade0c8
25 changed files with 4209 additions and 4 deletions

268
src/benchmark.rs Normal file
View File

@ -0,0 +1,268 @@
use crate::channel::{AwgnChannel, Channel};
use crate::code::{CodeTopology, GenerationMethod, LdpcCode, LdpcParams};
use crate::decoder::{build_decoder, DecoderConfig, DecoderMethod};
use crate::encoder::{build_encoder, EncodingMethod};
use crate::Result as LdpcResult;
use rand::{Rng, SeedableRng};
use std::fs::File;
use std::io::{Read, Write};
use std::path::Path;
use std::time::Instant;
pub fn run_simulation(mut code: LdpcCode) -> LdpcResult<()> {
println!("[*] Étape 1 : Construction Mathématique et Graphe de Tanner");
println!(
" - Dimensions : n={}, k={}, m={} (Taux R={:.3})",
code.n(),
code.k(),
code.m(),
code.rate()
);
let methode_nom = match code.params.generation {
GenerationMethod::MacKayNeal { .. } => "MacKay-Neal",
GenerationMethod::Gallager => "Gallager",
};
println!(" - Topologie : Régulier via {}", methode_nom);
println!(" - Densité H : {:.2}%", code.h.density() * 100.0);
println!(" - Girth : {}", code.girth());
println!(
" - Cycles-4 : {}",
if code.graph.has_4_cycle() {
"Présents (Problématique)"
} else {
"Aucun (Optimal)"
}
);
println!("\n[*] Étape 2 : Extraction de G^T et Instanciation de l'Encodeur");
let start_enc = Instant::now();
let encoder = build_encoder(&mut code, EncodingMethod::Systematic)?;
println!(
" - Forme systématique calculée en {:.2?}",
start_enc.elapsed()
);
println!("\n[*] Étape 3 : Instanciation des Décodeurs sur Graphe");
let config = DecoderConfig {
max_iterations: 50,
early_stopping: true,
};
let dec_sp = build_decoder(&code, DecoderMethod::SumProduct, config.clone());
let dec_ms = build_decoder(
&code,
DecoderMethod::MinSum {
scaling_factor: 0.8,
},
config.clone(),
);
let dec_bf = build_decoder(&code, DecoderMethod::BitFlipping, config.clone());
println!(" - Moteurs prêts : Sum-Product, Min-Sum (α=0.8), Bit-Flipping");
let snr_range = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0];
let n_trials = 100;
let mut rng = rand::rngs::StdRng::seed_from_u64(42);
println!(
"\n[*] Étape 4 : Simulation sur Canal AWGN ({} trames par SNR)",
n_trials
);
println!("{:-<115}", "");
println!(
"{:>8} | {:>9} || {:>10} | {:>10} || {:>10} | {:>10} || {:>10} | {:>10} |",
"SNR (dB)",
"Capacité",
"FER (SP)",
"BER (SP)",
"FER (MS)",
"BER (MS)",
"FER (BF)",
"BER (BF)"
);
println!("{:-<115}", "");
for &snr in &snr_range {
let channel = AwgnChannel::new(snr, code.rate())?;
let mut err_sp_frames = 0;
let mut err_sp_bits = 0;
let mut err_ms_frames = 0;
let mut err_ms_bits = 0;
let mut err_bf_frames = 0;
let mut err_bf_bits = 0;
for _ in 0..n_trials {
let message: Vec<u8> = (0..code.k()).map(|_| rng.gen::<u8>() & 1).collect();
let codeword = encoder.encode(&message)?;
let received_llr = channel.transmit(&codeword, &mut rng);
// SP
let res_sp = dec_sp.decode(&received_llr);
if let Some(decoded) = res_sp.codeword() {
let errs = count_bit_errors(&codeword, decoded);
if errs > 0 {
err_sp_frames += 1;
err_sp_bits += errs;
}
} else {
err_sp_frames += 1;
err_sp_bits += code.n();
}
// MS
let res_ms = dec_ms.decode(&received_llr);
if let Some(decoded) = res_ms.codeword() {
let errs = count_bit_errors(&codeword, decoded);
if errs > 0 {
err_ms_frames += 1;
err_ms_bits += errs;
}
} else {
err_ms_frames += 1;
err_ms_bits += code.n();
}
// BF
let res_bf = dec_bf.decode(&received_llr);
if let Some(decoded) = res_bf.codeword() {
let errs = count_bit_errors(&codeword, decoded);
if errs > 0 {
err_bf_frames += 1;
err_bf_bits += errs;
}
} else {
err_bf_frames += 1;
err_bf_bits += code.n();
}
}
let total_bits = (n_trials * code.n()) as f64;
let total_frames = n_trials as f64;
println!(
"{:>8.2} | {:>9.4} || {:>9.2}% | {:>9.2}% || {:>9.2}% | {:>9.2}% || {:>9.2}% | {:>9.2}% |",
snr, channel.capacity(),
(err_sp_frames as f64 / total_frames) * 100.0, (err_sp_bits as f64 / total_bits) * 100.0,
(err_ms_frames as f64 / total_frames) * 100.0, (err_ms_bits as f64 / total_bits) * 100.0,
(err_bf_frames as f64 / total_frames) * 100.0, (err_bf_bits as f64 / total_bits) * 100.0
);
}
println!("{:-<115}", "");
Ok(())
}
#[inline]
fn count_bit_errors(transmitted: &[u8], decoded: &[u8]) -> usize {
transmitted
.iter()
.zip(decoded.iter())
.filter(|(a, b)| a != b)
.count()
}
pub fn generate_valid_code(
n: usize,
k: usize,
wc: usize,
wr: usize,
generation_method: GenerationMethod,
) -> LdpcResult<LdpcCode> {
let mut attempt = 0;
let start_gen = Instant::now();
loop {
attempt += 1;
let params = LdpcParams {
n,
k,
topology: CodeTopology::Regular { wc, wr },
generation: generation_method.clone(),
seed: Some(rand::random()),
};
if let Ok(mut code) = LdpcCode::new(params) {
if code.compute_systematic_form().is_ok() {
if attempt > 1 {
println!(
" -> Matrice inversible obtenue après {} tentatives.",
attempt
);
}
println!(" - Génération : Terminée en {:.2?}", start_gen.elapsed());
return Ok(code);
}
}
}
}
pub fn get_or_generate_cached_code(
n: usize,
k: usize,
wc: usize,
wr: usize,
generation_method: GenerationMethod,
) -> LdpcResult<LdpcCode> {
let method_str = match generation_method {
GenerationMethod::MacKayNeal { .. } => "MN",
GenerationMethod::Gallager => "GAL",
};
let cache_filename = format!("cache_ldpc_{}_n{}_k{}.bin", method_str, n, k);
let path = Path::new(&cache_filename);
// Chargement
if path.exists() {
let start_load = Instant::now();
let mut file = File::open(&path).expect("Impossible d'ouvrir le fichier de cache");
let mut buffer = Vec::new();
file.read_to_end(&mut buffer).unwrap();
let mut code: LdpcCode = bincode::deserialize(&buffer).expect(
"Erreur de désérialisation du cache LDPC. Supprimez le fichier .bin et réessayez.",
);
// Reconstruction du Graphe de Tanner
code.graph = crate::graph::TannerGraph::from_matrix(&code.h);
println!(
" -> Matrice chargée depuis le cache en {:.2?}",
start_load.elapsed()
);
return Ok(code);
}
// Génération
println!(" -> Aucun cache trouvé. Génération en cours ...");
let mut attempt = 0;
let start_gen = Instant::now();
loop {
attempt += 1;
let params = LdpcParams {
n,
k,
topology: CodeTopology::Regular { wc, wr },
generation: generation_method.clone(),
seed: Some(rand::random()),
};
if let Ok(mut code) = LdpcCode::new(params) {
if code.compute_systematic_form().is_ok() {
println!(
" - Génération et Pivot de Gauss terminés en {:.2?}",
start_gen.elapsed()
);
// sauvegarde cache
let encoded = bincode::serialize(&code).expect("Échec de la sérialisation");
let mut file =
File::create(&path).expect("Impossible de créer le fichier de cache");
file.write_all(&encoded)
.expect("Impossible d'écrire sur le disque");
println!(" -> Matrice sauvegardée ({})", cache_filename);
return Ok(code);
}
}
}
}

74
src/channel.rs Normal file
View File

@ -0,0 +1,74 @@
use crate::{Gf2, LdpcError, Llr, Result};
// Trait Channel
pub trait Channel: Send + Sync {
fn transmit(&self, codeword: &[Gf2], rng: &mut impl rand::Rng) -> Vec<Llr>;
fn capacity(&self) -> f64;
}
// Canal AWGN
// BPSK mod 0 -> +1.0, 1 -> -1.0
// Signal recu y = x + n, n ~ N(0, sig^2)
// LLR optimal L(y) = 2y/sig^2
// sig^2 = 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<Self> {
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<Llr> {
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 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^2))]
sum += (1.0 + (-llr).exp()).log2();
}
1.0 - sum / n_samples as f64
}
}

329
src/code.rs Normal file
View File

@ -0,0 +1,329 @@
use crate::graph::TannerGraph;
use crate::matrix::{DenseMatrixGF2, SparseMatrixGF2};
use crate::{Gf2, LdpcError, Result};
use serde::{Deserialize, Serialize};
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct LdpcParams {
pub n: usize,
pub k: usize,
pub topology: CodeTopology,
pub generation: GenerationMethod,
pub seed: Option<u64>,
}
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, Serialize, Deserialize)]
pub enum CodeTopology {
// Chaque varnode a degré wc, chaque checknode a degré wr
// Condition nécessaire n * wc == m * wr
Regular { wc: usize, wr: usize },
// Potentielle implémentation Irregular (plus tard !)
}
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 ({}) != m*wr ({}) pour LDPC régulier",
params.n * wc,
params.m() * wr
)));
}
}
}
Ok(())
}
}
// Méthode de génération
#[derive(Debug, Clone, Serialize, Deserialize)]
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 cycle4 créé
MacKayNeal { max_attempts: usize },
}
// Forme systématique
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct SystematicForm {
// G = [I_k | P], dense
pub g: DenseMatrixGF2,
// Permutation de colonnes appliquée à H => [A | I_m]
pub col_perm: Vec<usize>,
// Permutation inverse => reformer le mot de code
pub col_perm_inv: Vec<usize>,
}
// Structure principale
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct LdpcCode {
pub params: LdpcParams,
pub h: SparseMatrixGF2,
#[serde(skip, default = "default_graph")]
pub graph: TannerGraph,
pub systematic_form: Option<SystematicForm>,
}
impl LdpcCode {
pub fn new(params: LdpcParams) -> Result<Self> {
params.validate()?;
let mut rng = build_rng(params.seed);
let h = match &params.generation {
GenerationMethod::Gallager => generate_gallager(&params, &mut rng)?,
GenerationMethod::MacKayNeal { max_attempts } => {
generate_mackay_neal(&params, *max_attempts, &mut rng)?
}
};
let graph = TannerGraph::from_matrix(&h);
Ok(Self {
params,
h,
graph,
systematic_form: None,
})
}
pub fn from_matrix(h: SparseMatrixGF2, k: usize) -> Result<Self> {
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)
pub fn compute_systematic_form(&mut self) -> Result<()> {
if self.systematic_form.is_some() {
return Ok(());
}
let mut dense = DenseMatrixGF2::zeros(self.m(), self.n());
let h_dense = self.h.to_dense();
for r in 0..self.m() {
for c in 0..self.n() {
dense.set(r, c, h_dense[r][c]);
}
}
let (col_perm, rank) = dense.systematize(self.k());
// Possibilité d'avoir rg < m donc juste rejeter et recommencer
if rank < self.m() {
return Err(LdpcError::SingularMatrix);
}
// Création de G^T (taille n,k) pour encodage c = G^T * m
let mut g_t = DenseMatrixGF2::zeros(self.n(), self.k());
// I_k (haut)
for i in 0..self.k() {
g_t.set(i, i, 1);
}
// (bas) Matrice A (extraite des k premières colonnes de la matrice dense pivotée)
for r in 0..self.m() {
for c in 0..self.k() {
g_t.set(self.k() + r, c, dense.get(r, c));
}
}
// Création de la permutation inverse pour réordonner le mot de code
let mut col_perm_inv = vec![0usize; self.n()];
for (i, &p) in col_perm.iter().enumerate() {
col_perm_inv[p] = i;
}
self.systematic_form = Some(SystematicForm {
g: g_t,
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)
}
}
fn build_rng(seed: Option<u64>) -> 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 = matrice régulière (ligne contient wr 1)
// H2..Hwc = permutations aléatoires de colonnes de H1
fn generate_gallager(params: &LdpcParams, rng: &mut impl rand::Rng) -> Result<SparseMatrixGF2> {
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();
for r in 0..rows_per_block {
for j in 0..wr {
ones.push((r, r * wr + j));
}
}
use rand::seq::SliceRandom;
for block in 1..wc {
let mut perm: Vec<usize> = (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]));
}
}
}
for b in 1..wc {
let r_source = b * rows_per_block;
let row_0_cols: Vec<usize> = ones
.iter()
.filter(|&&(r, _)| r == 0)
.map(|&(_, c)| c)
.collect();
if let Some(idx) = ones
.iter()
.position(|&(r, c)| r == r_source && !row_0_cols.contains(&c))
{
ones[idx].0 = 0;
}
}
Ok(SparseMatrixGF2::from_positions(m, n, ones))
}
// MacKay-Neal
// Ajoute les colonnes une à une avec poids wc
// Rejette colonne créant cycle4 (2 colonnes n'ont qu'un 1 en commun)
fn generate_mackay_neal(
params: &LdpcParams,
max_attempts: usize,
rng: &mut impl rand::Rng,
) -> Result<SparseMatrixGF2> {
let CodeTopology::Regular { wc, wr } = params.topology else {
return Err(LdpcError::InvalidParameters(
"MacKayNeal nécessite régulier".into(),
));
};
let n = params.n;
let m = params.m();
let mut ones: Vec<(usize, usize)> = Vec::new();
// Suivi du poids de chaque ligne
let mut row_weights = vec![0usize; m];
use rand::seq::SliceRandom;
for col in 0..n {
let mut placed = false;
for _attempt in 0..max_attempts {
let mut available_rows: Vec<usize> = (0..m).filter(|&r| row_weights[r] < wr).collect();
if available_rows.len() < wc {
break; // Plus de lignes dispo
}
available_rows.shuffle(rng);
let candidate = available_rows[..wc].to_vec();
// Vérifier cycle4
let mut ok = true;
let mut c2 = 0;
while c2 < col && ok {
let existing: Vec<usize> = 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));
row_weights[r] += 1;
}
placed = true;
break;
}
}
if !placed {
return Err(LdpcError::GenerationFailed {
attempts: max_attempts,
});
}
}
Ok(SparseMatrixGF2::from_positions(m, n, ones))
}
fn default_graph() -> TannerGraph {
TannerGraph::from_matrix(&SparseMatrixGF2::zeros(1, 1))
}

410
src/decoder.rs Normal file
View File

@ -0,0 +1,410 @@
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<Llr> = 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<Gf2> {
h.multiply_vec(c)
}
#[inline]
fn phi(x: Llr) -> Llr {
let ax = x.abs().max(1e-10);
-((ax / 2.0).tanh().ln())
}
// Mises à jour des noeuds
// Mise à jour Sum-Product du noeud de contrôle
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 a
// 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();
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 noeud de variable
fn variable_node_update(ch_llr: Llr, incoming_c2v: &[Llr], out: &mut [Llr]) {
let total: Llr = ch_llr + incoming_c2v.iter().sum::<Llr>();
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::<Llr>()
}
// Messages internes
// Indexés par (check_id, position_dans_liste_voisins) (O(1))
struct Messages {
v2c: Vec<Vec<Llr>>,
c2v: Vec<Vec<Llr>>,
}
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é une fois à construction du décodeur
struct EdgeIndex {
var_pos_in_chk: Vec<Vec<usize>>,
chk_pos_in_var: Vec<Vec<usize>>,
}
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<Gf2> = 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)
}
}
}
// BP
fn bp_decode<F>(
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);
// Init
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 {
// Maj des checknodes
for c in 0..graph.n_chk {
let v2c = msgs.v2c[c].clone();
check_update(&v2c, &mut msgs.c2v[c]);
}
// Maj des varnodes
for v in 0..graph.n_var {
let neighbors = graph.var_neighbors(v);
// Rassembler les c2v entrants sur ce varnode
let incoming: Vec<Llr> = 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<Gf2> {
(0..graph.n_var)
.map(|v| {
let incoming: Vec<Llr> = 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<dyn Decoder> {
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))
}
}
}

89
src/encoder.rs Normal file
View File

@ -0,0 +1,89 @@
use crate::code::LdpcCode;
use crate::matrix::DenseMatrixGF2;
use crate::{BitVec, Gf2, LdpcError, Result};
pub trait Encoder: Send + Sync {
fn encode(&self, message: &[Gf2]) -> Result<BitVec>;
fn message_len(&self) -> usize;
fn codeword_len(&self) -> usize;
fn extract_message(&self, codeword: &[Gf2]) -> Vec<Gf2>;
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(())
}
}
#[derive(Debug, Clone)]
pub enum EncodingMethod {
Systematic,
}
pub struct SystematicEncoder {
k: usize,
n: usize,
g_t: DenseMatrixGF2,
perm_inv: Vec<usize>,
col_perm: Vec<usize>,
}
impl SystematicEncoder {
pub fn new(code: &mut LdpcCode) -> Result<Self> {
code.compute_systematic_form()?;
let sf = code.systematic_form.as_ref().unwrap();
Ok(Self {
k: code.k(),
n: code.n(),
g_t: sf.g.clone(),
perm_inv: sf.col_perm_inv.clone(),
col_perm: sf.col_perm.clone(),
})
}
}
impl Encoder for SystematicEncoder {
fn encode(&self, message: &[Gf2]) -> Result<BitVec> {
self.check_input(message)?;
let c_perm = self.g_t.multiply_vec(message);
// Retablir l'ordre initial des bits selon la permutation de H
let mut c = vec![0u8; self.n];
// for (i, &ci) in c_perm.iter().enumerate() {
// c[self.perm_inv[i]] = ci;
// }
for i in 0..self.n {
c[i] = c_perm[self.perm_inv[i]];
}
Ok(c)
}
fn extract_message(&self, codeword: &[Gf2]) -> Vec<Gf2> {
let mut msg = vec![0u8; self.k];
for j in 0..self.k {
msg[j] = codeword[self.col_perm[j]];
}
msg
}
fn message_len(&self) -> usize {
self.k
}
fn codeword_len(&self) -> usize {
self.n
}
}
pub fn build_encoder(code: &mut LdpcCode, method: EncodingMethod) -> Result<Box<dyn Encoder>> {
match method {
EncodingMethod::Systematic => Ok(Box::new(SystematicEncoder::new(code)?)),
}
}

142
src/graph.rs Normal file
View File

@ -0,0 +1,142 @@
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<Vec<usize>>,
chk_to_var: Vec<Vec<usize>>,
}
impl TannerGraph {
pub fn from_matrix(h: &SparseMatrixGF2) -> Self {
let n_var = h.cols;
let n_chk = h.rows;
let chk_to_var: Vec<Vec<usize>> = (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
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
}
}
min_girth
}
// Détection cycles-4, 2 varnodes 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
}
pub fn local_girth_from_var(&self, v: usize) -> usize {
self.bfs_girth_from_var(v).unwrap_or(usize::MAX)
}
// retourne la longueur du court cycle passant par ce noeud (None si pas cycle)
fn bfs_girth_from_var(&self, start: usize) -> Option<usize> {
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, index, distance, parent_index)
let mut queue: VecDeque<(bool, usize, usize, usize)> = VecDeque::new();
queue.push_back((true, start, 0, usize::MAX));
let mut shortest = None;
while let Some((is_var, node, dist, parent)) = queue.pop_front() {
if is_var {
for &c in self.var_neighbors(node) {
if c == parent {
continue;
} // aller retour immédiat impossible
if dist_chk[c] == usize::MAX {
dist_chk[c] = dist + 1;
queue.push_back((false, c, dist + 1, node));
} else {
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 == parent {
continue;
} // aller retour immédiat impossible
if v == start {
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, node));
}
}
}
}
shortest
}
pub fn var_degree_distribution(&self) -> Vec<f64> {
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)
}
}

131
src/image_sim.rs Normal file
View File

@ -0,0 +1,131 @@
use crate::{
channel::{AwgnChannel, Channel},
decoder::Decoder,
encoder::Encoder,
Gf2, Result,
};
use image::{ImageBuffer, Rgb};
use rand::SeedableRng;
// Convertit un tableau d'octets en un flux de bits
pub fn bytes_to_bits(bytes: &[u8]) -> Vec<Gf2> {
let mut bits = Vec::with_capacity(bytes.len() * 8);
for &byte in bytes {
for i in (0..8).rev() {
bits.push((byte >> i) & 1);
}
}
bits
}
// Convertit un flux de bits en tableau d'octets
pub fn bits_to_bytes(bits: &[Gf2]) -> Vec<u8> {
let mut bytes = Vec::with_capacity(bits.len() / 8);
for chunk in bits.chunks(8) {
let mut byte = 0u8;
for (i, &bit) in chunk.iter().enumerate() {
byte |= bit << (7 - i);
}
bytes.push(byte);
}
bytes
}
// Transmet une image à travers le canal avec codage LDPC
pub fn transmit_image(
input_path: &str,
noisy_out_path: &str,
decoded_out_path: &str,
encoder: &dyn Encoder,
decoder: &dyn Decoder,
channel: &AwgnChannel,
) -> Result<()> {
println!("[*] Chargement de l'image : {}", input_path);
let img = image::open(input_path)
.expect("Erreur de chargement de l'image")
.to_rgb8();
let (width, height) = img.dimensions();
let raw_bytes = img.into_raw();
let mut bits = bytes_to_bits(&raw_bytes);
let original_bit_len = bits.len();
// Padding
let k = encoder.message_len();
let remainder = bits.len() % k;
if remainder != 0 {
bits.resize(bits.len() + (k - remainder), 0);
}
let num_blocks = bits.len() / k;
println!(" - Taille: {}x{} pixels", width, height);
println!(" - Blocs à transmettre (k={}): {}", k, num_blocks);
let mut rng = rand::rngs::StdRng::seed_from_u64(42);
let mut noisy_bits = Vec::with_capacity(num_blocks * k);
let mut decoded_bits = Vec::with_capacity(num_blocks * k);
let mut frame_errors = 0;
println!("[*] Transmission et Décodage en cours...");
for (i, block) in bits.chunks(k).enumerate() {
if i % 100 == 0 && i > 0 {
println!(" - Progession: {} / {} blocs...", i, num_blocks);
}
let codeword = encoder.encode(block)?;
let rx_llr = channel.transmit(&codeword, &mut rng);
// Sans correction LDPC
let hard_codeword: Vec<Gf2> = rx_llr
.iter()
.map(|&l| if l < 0.0 { 1 } else { 0 })
.collect();
let noisy_block = encoder.extract_message(&hard_codeword);
noisy_bits.extend_from_slice(&noisy_block);
// Décodage LDPC
let res = decoder.decode(&rx_llr);
if let Some(decoded_codeword) = res.codeword() {
let decoded_msg = encoder.extract_message(decoded_codeword);
decoded_bits.extend_from_slice(&decoded_msg);
if decoded_msg != block {
frame_errors += 1;
}
} else {
decoded_bits.extend_from_slice(&noisy_block);
frame_errors += 1;
}
}
println!(
"[*] Transmission terminée. FER : {:.2}%",
(frame_errors as f64 / num_blocks as f64) * 100.0
);
// Suppression du padding
noisy_bits.truncate(original_bit_len);
decoded_bits.truncate(original_bit_len);
// Reconstitution des images
let noisy_bytes = bits_to_bytes(&noisy_bits);
let decoded_bytes = bits_to_bytes(&decoded_bits);
let noisy_img = ImageBuffer::<Rgb<u8>, _>::from_raw(width, height, noisy_bytes)
.expect("Erreur de reconstruction de l'image bruitée");
noisy_img.save(noisy_out_path).unwrap();
let decoded_img = ImageBuffer::<Rgb<u8>, _>::from_raw(width, height, decoded_bytes)
.expect("Erreur de reconstruction de l'image décodée");
decoded_img.save(decoded_out_path).unwrap();
println!(
"[*] Images sauvegardées : {} et {}",
noisy_out_path, decoded_out_path
);
Ok(())
}

40
src/lib.rs Normal file
View File

@ -0,0 +1,40 @@
pub mod benchmark;
pub mod channel;
pub mod code;
pub mod decoder;
pub mod encoder;
pub mod graph;
pub mod image_sim;
pub mod matrix;
pub type Gf2 = u8;
pub type Llr = f64;
pub type BitVec = Vec<Gf2>;
#[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<T> = std::result::Result<T, LdpcError>;
pub use channel::Channel;
pub use code::{CodeTopology, GenerationMethod, LdpcCode, LdpcParams};
pub use decoder::{Decoder, DecoderConfig, DecoderMethod, DecoderResult};
pub use encoder::{Encoder, EncodingMethod};

View File

@ -1,3 +1,76 @@
fn main() {
println!("Hello, world!");
// use ldpc::benchmark::{generate_valid_code, run_simulation};
// use ldpc::code::GenerationMethod;
//
// fn main() -> ldpc::Result<()> {
// let n = 1944;
// let k = 972;
// let wc = 3;
// let wr = 6;
//
// println!("Benchmark: MacKayNeal vs Gallager");
// println!();
//
// println!("Test 1: Génération MacKayNeal\n");
// let code_mn = generate_valid_code(
// n,
// k,
// wc,
// wr,
// GenerationMethod::MacKayNeal { max_attempts: 1000 },
// )?;
// run_simulation(code_mn)?;
//
// println!("\nTest 2 : Génération Gallager\n");
// let code_gal = generate_valid_code(n, k, wc, wr, GenerationMethod::Gallager)?;
// run_simulation(code_gal)?;
//
// Ok(())
// }
use ldpc::benchmark::get_or_generate_cached_code;
use ldpc::channel::AwgnChannel;
use ldpc::code::GenerationMethod;
use ldpc::decoder::{build_decoder, DecoderConfig, DecoderMethod};
use ldpc::encoder::{build_encoder, EncodingMethod};
use ldpc::image_sim::transmit_image;
fn main() -> ldpc::Result<()> {
let n = 1296;
let k = 864;
let wc = 3;
let wr = 6;
println!("Transmission d'image via code LDPC");
let code_mn = get_or_generate_cached_code(
n,
k,
wc,
wr,
GenerationMethod::MacKayNeal { max_attempts: 5000 },
)?;
let mut code = code_mn;
let encoder = build_encoder(&mut code, EncodingMethod::Systematic)?;
let config = DecoderConfig {
max_iterations: 50,
early_stopping: true,
};
// Sum-Product
let decoder = build_decoder(&code, DecoderMethod::SumProduct, config);
let channel = AwgnChannel::new(2.0, code.rate())?;
transmit_image(
"test.png",
"noisy_out.png",
"decoded_out.png",
&*encoder,
&*decoder,
&channel,
)?;
Ok(())
}

284
src/matrix.rs Normal file
View File

@ -0,0 +1,284 @@
use crate::Gf2;
use serde::{Deserialize, Serialize};
// Matrice creuse format CSR + CSC
#[derive(Debug, Clone, Serialize, Deserialize)]
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<usize>,
col_idx: Vec<usize>,
// CSC accès col j : row_idx[col_ptr[j]..col_ptr[j+1]]
col_ptr: Vec<usize>,
row_idx: Vec<usize>,
}
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<Gf2>]) -> 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
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
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 (syndrome : s = H * c^T)
pub fn multiply_vec(&self, x: &[Gf2]) -> Vec<Gf2> {
(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<Vec<Gf2>> {
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
// Utilisée pour G et Gauss-Jordan
// G = [I | P], P dense
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct DenseMatrixGF2 {
pub rows: usize,
pub cols: usize,
data: Vec<Vec<Gf2>>,
}
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;
}
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 col_swap(&mut self, c1: usize, c2: usize) {
for r in 0..self.rows {
let tmp = self.data[r][c1];
self.data[r][c1] = self.data[r][c2];
self.data[r][c2] = tmp;
}
}
pub fn multiply_vec(&self, x: &[Gf2]) -> Vec<Gf2> {
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)
}
pub fn systematize(&mut self, k: usize) -> (Vec<usize>, usize) {
let m = self.rows;
let n = self.cols;
let mut col_perm: Vec<usize> = (0..n).collect();
let mut rank = 0;
for i in 0..m {
// Placer le pivot ligne i à la colonne target_c (former I_m à droite)
let target_c = k + i;
let mut pivot_found = false;
// pivot, chercher en premier dans les colonnes cibles restantes
// et après dans les colonnes de données (0..k)
for c_search in (target_c..n).chain(0..k) {
for r_search in i..m {
if self.data[r_search][c_search] == 1 {
// pivot à la position (i, target_c)
self.row_swap(i, r_search);
if c_search != target_c {
self.col_swap(target_c, c_search);
col_perm.swap(target_c, c_search);
}
pivot_found = true;
break;
}
}
if pivot_found {
break;
}
}
if pivot_found {
rank += 1;
// Elimination dans toutes les autres lignes => forcer la colonne à 0 (sauf le pivot à 1)
for r in 0..m {
if r != i && self.data[r][target_c] == 1 {
self.row_add(r, i);
}
}
}
}
(col_perm, rank)
}
}