diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..afd0b22 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,499 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "adler2" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" + +[[package]] +name = "bincode" +version = "1.3.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b1f45e9417d87227c7a56d22e471c6206462cba514c7590c09aff4cf6d1ddcad" +dependencies = [ + "serde", +] + +[[package]] +name = "bit_field" +version = "0.10.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e4b40c7323adcfc0a41c4b88143ed58346ff65a288fc144329c5c45e05d70c6" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "bytemuck" +version = "1.25.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8efb64bd706a16a1bdde310ae86b351e4d21550d98d056f22f8a7f7a2183fec" + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cfg-if" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" + +[[package]] +name = "color_quant" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d7b894f5411737b7867f4827955924d7c254fc9f4d91a6aad6b097804b1018b" + +[[package]] +name = "crc32fast" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9481c1c90cbf2ac953f07c8d4a58aa3945c425b7185c9154d67a65e4230da511" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "crossbeam-deque" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +dependencies = [ + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28" + +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + +[[package]] +name = "either" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" + +[[package]] +name = "exr" +version = "1.74.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4300e043a56aa2cb633c01af81ca8f699a321879a7854d3896a0ba89056363be" +dependencies = [ + "bit_field", + "half", + "lebe", + "miniz_oxide", + "rayon-core", + "smallvec", + "zune-inflate", +] + +[[package]] +name = "fdeflate" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e6853b52649d4ac5c0bd02320cddc5ba956bdb407c4b75a2c6b75bf51500f8c" +dependencies = [ + "simd-adler32", +] + +[[package]] +name = "flate2" +version = "1.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "843fba2746e448b37e26a819579957415c8cef339bf08564fe8b7ddbd959573c" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + +[[package]] +name = "getrandom" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + +[[package]] +name = "gif" +version = "0.13.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ae047235e33e2829703574b54fdec96bfbad892062d97fed2f76022287de61b" +dependencies = [ + "color_quant", + "weezl", +] + +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "cfg-if", + "crunchy", + "zerocopy", +] + +[[package]] +name = "image" +version = "0.24.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5690139d2f55868e080017335e4b94cb7414274c74f1669c84fb5feba2c9f69d" +dependencies = [ + "bytemuck", + "byteorder", + "color_quant", + "exr", + "gif", + "jpeg-decoder", + "num-traits", + "png", + "qoi", + "tiff", +] + +[[package]] +name = "jpeg-decoder" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "00810f1d8b74be64b13dbf3db89ac67740615d6c891f0e7b6179326533011a07" +dependencies = [ + "rayon", +] + +[[package]] +name = "ldpc" +version = "0.1.0" +dependencies = [ + "approx", + "bincode", + "image", + "rand", + "rand_distr", + "serde", + "thiserror", +] + +[[package]] +name = "lebe" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7a79a3332a6609480d7d0c9eab957bca6b455b91bb84e66d19f5ff66294b85b8" + +[[package]] +name = "libc" +version = "0.2.186" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68ab91017fe16c622486840e4c83c9a37afeff978bd239b5293d61ece587de66" + +[[package]] +name = "libm" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" + +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", + "simd-adler32", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", + "libm", +] + +[[package]] +name = "png" +version = "0.17.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "82151a2fc869e011c153adc57cf2789ccb8d9906ce52c0b39a6b5697749d7526" +dependencies = [ + "bitflags", + "crc32fast", + "fdeflate", + "flate2", + "miniz_oxide", +] + +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + +[[package]] +name = "proc-macro2" +version = "1.0.106" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8fd00f0bb2e90d81d1044c2b32617f68fcb9fa3bb7640c23e9c748e53fb30934" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "qoi" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f6d64c71eb498fe9eae14ce4ec935c555749aef511cca85b5568910d6e48001" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "quote" +version = "1.0.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "41f2619966050689382d2b44f664f4bc593e129785a36d6ee376ddf37259b924" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rand" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ca0ecfa931c29007047d1bc58e623ab12e5590e8c7cc53200d5202b69266d8a" +dependencies = [ + "libc", + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rand_distr" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31" +dependencies = [ + "num-traits", + "rand", +] + +[[package]] +name = "rayon" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + +[[package]] +name = "serde" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a8e94ea7f378bd32cbbd37198a4a91436180c5bb472411e48b5ec2e2124ae9e" +dependencies = [ + "serde_core", + "serde_derive", +] + +[[package]] +name = "serde_core" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "41d385c7d4ca58e59fc732af25c3983b67ac852c1a25000afe1175de458b67ad" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "simd-adler32" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "703d5c7ef118737c72f1af64ad2f6f8c5e1921f818cdcb97b8fe6fc69bf66214" + +[[package]] +name = "smallvec" +version = "1.15.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "67b1b7a3b5fe4f1376887184045fcf45c69e92af734b7aaddc05fb777b6fbd03" + +[[package]] +name = "syn" +version = "2.0.117" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e665b8803e7b1d2a727f4023456bbbbe74da67099c585258af0ad9c5013b9b99" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "tiff" +version = "0.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba1310fcea54c6a9a4fd1aad794ecc02c31682f6bfbecdf460bf19533eed1e3e" +dependencies = [ + "flate2", + "jpeg-decoder", + "weezl", +] + +[[package]] +name = "unicode-ident" +version = "1.0.24" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75" + +[[package]] +name = "wasi" +version = "0.11.1+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b" + +[[package]] +name = "weezl" +version = "0.1.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a28ac98ddc8b9274cb41bb4d9d4d5c425b6020c50c46f25559911905610b4a88" + +[[package]] +name = "zerocopy" +version = "0.8.48" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eed437bf9d6692032087e337407a86f04cd8d6a16a37199ed57949d415bd68e9" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.48" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "70e3cd084b1788766f53af483dd21f93881ff30d7320490ec3ef7526d203bad4" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "zune-inflate" +version = "0.2.54" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73ab332fe2f6680068f3582b16a24f90ad7096d5d39b974d1c0aff0125116f02" +dependencies = [ + "simd-adler32", +] diff --git a/Cargo.toml b/Cargo.toml index 1cd198f..924849d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,15 @@ [package] -name = "ldpc" +name = "ldpc" version = "0.1.0" -edition = "2024" +edition = "2021" [dependencies] +image = "0.24" +rand = "0.8" +rand_distr = "0.4.3" +thiserror = "1.0" +serde = { version = "1.0", features = ["derive"] } +bincode = "1.3" + +[dev-dependencies] +approx = "0.5" diff --git a/cache_ldpc_MN_n1944_k972.bin b/cache_ldpc_MN_n1944_k972.bin new file mode 100644 index 0000000..dfa3dc3 Binary files /dev/null and b/cache_ldpc_MN_n1944_k972.bin differ diff --git a/cache_ldpc_MN_n4096_k2048.bin b/cache_ldpc_MN_n4096_k2048.bin new file mode 100644 index 0000000..9ce431e Binary files /dev/null and b/cache_ldpc_MN_n4096_k2048.bin differ diff --git a/decoded_out.png b/decoded_out.png new file mode 100644 index 0000000..56cd6b9 Binary files /dev/null and b/decoded_out.png differ diff --git a/noisy_out.png b/noisy_out.png new file mode 100644 index 0000000..fa7cb55 Binary files /dev/null and b/noisy_out.png differ diff --git a/src/benchmark.rs b/src/benchmark.rs new file mode 100644 index 0000000..d800e56 --- /dev/null +++ b/src/benchmark.rs @@ -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 = (0..code.k()).map(|_| rng.gen::() & 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 { + 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 { + 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); + } + } + } +} diff --git a/src/channel.rs b/src/channel.rs new file mode 100644 index 0000000..0b69b74 --- /dev/null +++ b/src/channel.rs @@ -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; + 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 { + 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 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 + } +} diff --git a/src/code.rs b/src/code.rs new file mode 100644 index 0000000..bc21084 --- /dev/null +++ b/src/code.rs @@ -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, +} + +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, + // Permutation inverse => reformer le mot de code + pub col_perm_inv: Vec, +} + +// 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, +} + +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)? + } + }; + 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) + 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) -> 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 { + 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 = (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 = 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 { + 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 = (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 = 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)) +} diff --git a/src/decoder.rs b/src/decoder.rs new file mode 100644 index 0000000..82153e8 --- /dev/null +++ b/src/decoder.rs @@ -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 = 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) +} + +#[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::(); + 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) (O(1)) +struct Messages { + v2c: Vec>, + c2v: Vec>, +} + +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>, + 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) + } + } +} + +// BP +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); + + // 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 = 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)) + } + } +} diff --git a/src/encoder.rs b/src/encoder.rs new file mode 100644 index 0000000..6193801 --- /dev/null +++ b/src/encoder.rs @@ -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; + fn message_len(&self) -> usize; + fn codeword_len(&self) -> usize; + fn extract_message(&self, codeword: &[Gf2]) -> Vec; + + 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, + col_perm: 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_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 { + 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 { + 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> { + match method { + EncodingMethod::Systematic => Ok(Box::new(SystematicEncoder::new(code)?)), + } +} diff --git a/src/graph.rs b/src/graph.rs new file mode 100644 index 0000000..ca49da6 --- /dev/null +++ b/src/graph.rs @@ -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>, + 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 + 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 { + 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 { + 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) + } +} diff --git a/src/image_sim.rs b/src/image_sim.rs new file mode 100644 index 0000000..0deb354 --- /dev/null +++ b/src/image_sim.rs @@ -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 { + 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 { + 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 = 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::, _>::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::, _>::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(()) +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..0968aab --- /dev/null +++ b/src/lib.rs @@ -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; + +#[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/src/main.rs b/src/main.rs index e7a11a9..908a617 100644 --- a/src/main.rs +++ b/src/main.rs @@ -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(()) } diff --git a/src/matrix.rs b/src/matrix.rs new file mode 100644 index 0000000..9df3ad9 --- /dev/null +++ b/src/matrix.rs @@ -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, + 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 + 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 { + (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 +// 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>, +} + +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 { + 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) { + let m = self.rows; + let n = self.cols; + let mut col_perm: Vec = (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) + } +} 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)); +// } +// } diff --git a/test.png b/test.png new file mode 100644 index 0000000..a9e8789 Binary files /dev/null and b/test.png differ