diff --git a/Cargo.toml b/Cargo.toml
index 1479091..4debf7d 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -11,6 +11,7 @@ readme = "README.md"
 
 [features]
 progress_bar = ["dep:indicatif"]
+quicklog = ["dep:thiserror"]
 
 [dependencies]
 serde_json = "1.0.115"
@@ -21,7 +22,7 @@ rayon = "1.10.0"
 indicatif = { version = "0.17.8", optional = true }
 ndarray-npy ="0.8.1"
 itertools = "0.12.1"
-thiserror = "1.0.58"
+thiserror = { version = "1.0.58", optional = true }
 
 [dev-dependencies]
 criterion = "0.5.1"
@@ -29,6 +30,11 @@ ndarray-rand = "0.14.0"
 anyhow = "1.0.81"
 muscat = { path = ".", features = ["progress_bar"] }
 
+[[example]]
+name = "snr"
+path = "examples/snr.rs"
+required-features = ["quicklog"]
+
 [[bench]]
 name = "cpa"
 harness = false
diff --git a/benches/cpa.rs b/benches/cpa.rs
index 5dd2845..1aff2cf 100644
--- a/benches/cpa.rs
+++ b/benches/cpa.rs
@@ -29,7 +29,7 @@ pub fn leakage_model_normal(value: ArrayView1<usize>, guess: usize) -> usize {
     hw(sbox((value[1] ^ guess) as u8) as usize)
 }
 
-fn cpa_normal_sequential(leakages: &Array2<f64>, plaintexts: &Array2<u8>) -> cpa_normal::Cpa {
+fn cpa_normal_sequential(leakages: &Array2<f64>, plaintexts: &Array2<u8>) -> Cpa {
     let chunk_size = 500;
 
     let mut cpa =
diff --git a/benches/dpa.rs b/benches/dpa.rs
index 9aad207..6fadbd1 100644
--- a/benches/dpa.rs
+++ b/benches/dpa.rs
@@ -6,8 +6,8 @@ use ndarray_rand::rand::{rngs::StdRng, SeedableRng};
 use ndarray_rand::rand_distr::Uniform;
 use ndarray_rand::RandomExt;
 
-fn selection_function(metadata: Array1<u8>, guess: usize) -> usize {
-    sbox(metadata[1] ^ guess as u8).into()
+fn selection_function(metadata: Array1<u8>, guess: usize) -> bool {
+    usize::from(sbox(metadata[1] ^ guess as u8)) & 1 == 1
 }
 
 fn dpa_sequential(leakages: &Array2<f32>, plaintexts: &Array2<u8>) -> Dpa {
diff --git a/examples/cpa.rs b/examples/cpa.rs
index 17eaa08..b2357a1 100644
--- a/examples/cpa.rs
+++ b/examples/cpa.rs
@@ -4,6 +4,7 @@ use muscat::cpa_normal::CpaProcessor;
 use muscat::leakage::{hw, sbox};
 use muscat::util::{progress_bar, read_array2_from_npy_file, save_array2};
 use ndarray::*;
+use ndarray_npy::write_npy;
 use rayon::iter::{ParallelBridge, ParallelIterator};
 use std::time;
 
@@ -50,9 +51,9 @@ fn cpa() -> Result<()> {
         );
 
     let cpa = cpa_parallel.finalize();
-    println!("Guessed key = {}", cpa.pass_guess());
+    println!("Guessed key = {}", cpa.best_guess());
 
-    save_array2("results/corr.npy", cpa.pass_corr_array().view())?;
+    save_array2("results/corr.npy", cpa.corr())?;
 
     Ok(())
 }
@@ -70,30 +71,38 @@ fn success() -> Result<()> {
 
     let mut cpa = CpaProcessor::new(size, batch, guess_range, leakage_model);
 
-    cpa.success_traces(rank_traces);
+    let mut rank = Array1::zeros(guess_range);
+    let mut processed_traces = 0;
     for i in (0..nfiles).progress() {
         let dir_l = format!("{folder}/l/{i}.npy");
         let dir_p = format!("{folder}/p/{i}.npy");
         let leakages = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
         let plaintext = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
-        let len_leakages = leakages.shape()[0];
-        for row in (0..len_leakages).step_by(batch) {
+        for row in (0..leakages.shape()[0]).step_by(batch) {
             let range_samples = start_sample..end_sample;
-            let range_rows: std::ops::Range<usize> = row..row + batch;
-            let range_metadat = 0..plaintext.shape()[1];
+            let range_rows = row..row + batch;
+            let range_metadata = 0..plaintext.shape()[1];
             let sample_traces = leakages
                 .slice(s![range_rows.clone(), range_samples])
                 .map(|l| *l as f32);
-            let sample_metadata = plaintext.slice(s![range_rows, range_metadat]);
-            cpa.update_success(sample_traces.view(), sample_metadata);
+            let sample_metadata = plaintext.slice(s![range_rows, range_metadata]);
+
+            cpa.update(sample_traces.view(), sample_metadata);
+            processed_traces += sample_traces.len();
+            if processed_traces % rank_traces == 0 {
+                // rank can be saved to get its evolution
+                rank = cpa.finalize().rank();
+            }
         }
     }
 
     let cpa = cpa.finalize();
-    println!("Guessed key = {}", cpa.pass_guess());
+    println!("Guessed key = {}", cpa.best_guess());
+
+    println!("{:?}", rank);
 
     // save corr key curves in npy
-    save_array2("results/success.npy", cpa.pass_rank().view())?;
+    write_npy("results/success.npy", &cpa.rank().map(|&x| x as u64))?;
 
     Ok(())
 }
diff --git a/examples/cpa_partioned.rs b/examples/cpa_partioned.rs
index e1ad877..4a2122c 100644
--- a/examples/cpa_partioned.rs
+++ b/examples/cpa_partioned.rs
@@ -46,14 +46,14 @@ fn cpa() -> Result<()> {
         })
         .reduce(
             || CpaProcessor::new(size, guess_range, target_byte, leakage_model),
-            |a: CpaProcessor, b| a + b,
+            |a, b| a + b,
         );
 
     let cpa_result = cpa.finalize();
-    println!("Guessed key = {}", cpa_result.pass_guess());
+    println!("Guessed key = {}", cpa_result.best_guess());
 
     // save corr key curves in npy
-    save_array("../results/corr.npy", &cpa_result.pass_corr_array())?;
+    save_array("../results/corr.npy", &cpa_result.corr())?;
 
     Ok(())
 }
diff --git a/examples/dpa.rs b/examples/dpa.rs
index 370c829..f450344 100644
--- a/examples/dpa.rs
+++ b/examples/dpa.rs
@@ -3,16 +3,15 @@ use indicatif::ProgressIterator;
 use muscat::dpa::DpaProcessor;
 use muscat::leakage::sbox;
 use muscat::util::read_array2_from_npy_file;
-use ndarray::{s, Array1, Array2};
+use ndarray::{s, Array1};
 use rayon::iter::{ParallelBridge, ParallelIterator};
 
 // traces format
 type FormatTraces = f64;
 type FormatMetadata = u8;
 
-// leakage model
-pub fn leakage_model(value: Array1<FormatMetadata>, guess: usize) -> usize {
-    (sbox((value[1] as usize ^ guess) as u8)) as usize
+pub fn selection_function(value: Array1<FormatMetadata>, guess: usize) -> bool {
+    (sbox((value[1] as usize ^ guess) as u8)) as usize & 1 == 1
 }
 
 fn dpa() -> Result<()> {
@@ -23,21 +22,21 @@ fn dpa() -> Result<()> {
     let folder = String::from("../../data/cw");
     let dir_l = format!("{folder}/leakages.npy");
     let dir_p = format!("{folder}/plaintexts.npy");
-    let leakages: Array2<FormatTraces> = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
-    let plaintext: Array2<FormatMetadata> = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
+    let leakages = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
+    let plaintext = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
     let len_traces = 20000; //leakages.shape()[0];
-    let mut dpa = DpaProcessor::new(size, guess_range, leakage_model);
+    let mut dpa_proc = DpaProcessor::new(size, guess_range, selection_function);
     for i in (0..len_traces).progress() {
         let tmp_trace = leakages
             .row(i)
             .slice(s![start_sample..end_sample])
-            .map(|t| *t as f32);
+            .mapv(|t| t as f32);
         let tmp_metadata = plaintext.row(i);
-        dpa.update(tmp_trace.view(), tmp_metadata.to_owned());
+        dpa_proc.update(tmp_trace.view(), tmp_metadata.to_owned());
     }
-    let dpa = dpa.finalize();
-    println!("Guessed key = {:02x}", dpa.pass_guess());
-    // let corr = dpa.pass_corr_array();
+    let dpa = dpa_proc.finalize();
+    println!("Guessed key = {:02x}", dpa.best_guess());
+    // let corr = dpa.corr();
 
     Ok(())
 }
@@ -51,25 +50,30 @@ fn dpa_success() -> Result<()> {
     let folder = String::from("../../data/cw");
     let dir_l = format!("{folder}/leakages.npy");
     let dir_p = format!("{folder}/plaintexts.npy");
-    let leakages: Array2<FormatTraces> = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
-    let plaintext: Array2<FormatMetadata> = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
+    let leakages = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
+    let plaintext = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
     let len_traces = leakages.shape()[0];
-    let mut dpa = DpaProcessor::new(size, guess_range, leakage_model);
+    let mut dpa_proc = DpaProcessor::new(size, guess_range, selection_function);
     let rank_traces: usize = 100;
-    dpa.assign_rank_traces(rank_traces);
 
+    let mut rank = Array1::zeros(guess_range);
     for i in (0..len_traces).progress() {
         let tmp_trace = leakages
             .row(i)
             .slice(s![start_sample..end_sample])
-            .map(|t| *t as f32);
+            .mapv(|t| t as f32);
         let tmp_metadata = plaintext.row(i).to_owned();
-        dpa.update_success(tmp_trace.view(), tmp_metadata);
+        dpa_proc.update(tmp_trace.view(), tmp_metadata);
+
+        if i % rank_traces == 0 {
+            // rank can be saved to get its evolution
+            rank = dpa_proc.finalize().rank();
+        }
     }
 
-    let dpa = dpa.finalize();
-    println!("Guessed key = {:02x}", dpa.pass_guess());
-    // let succss = dpa.pass_rank().to_owned();
+    let dpa = dpa_proc.finalize();
+    println!("Guessed key = {:02x}", dpa.best_guess());
+    println!("{:?}", rank);
 
     Ok(())
 }
@@ -83,21 +87,22 @@ fn dpa_parallel() -> Result<()> {
     let folder = String::from("../../data/cw");
     let dir_l = format!("{folder}/leakages.npy");
     let dir_p = format!("{folder}/plaintexts.npy");
-    let leakages: Array2<FormatTraces> = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
-    let plaintext: Array2<FormatMetadata> = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
+    let leakages = read_array2_from_npy_file::<FormatTraces>(&dir_l)?;
+    let plaintext = read_array2_from_npy_file::<FormatMetadata>(&dir_p)?;
     let len_traces = 20000; //leakages.shape()[0];
     let batch = 2500;
-    let mut dpa = (0..len_traces)
+    let dpa = (0..len_traces)
         .step_by(batch)
         .par_bridge()
-        .map(|range_rows: usize| {
+        .map(|range_rows| {
             let tmp_leakages = leakages
                 .slice(s![range_rows..range_rows + batch, start_sample..end_sample])
-                .map(|l| *l as f32);
+                .mapv(|l| l as f32);
             let tmp_metadata = plaintext
                 .slice(s![range_rows..range_rows + batch, ..])
                 .to_owned();
-            let mut dpa_inner = DpaProcessor::new(size, guess_range, leakage_model);
+
+            let mut dpa_inner = DpaProcessor::new(size, guess_range, selection_function);
             for i in 0..batch {
                 let trace = tmp_leakages.row(i);
                 let metadata = tmp_metadata.row(i).to_owned();
@@ -106,13 +111,13 @@ fn dpa_parallel() -> Result<()> {
             dpa_inner
         })
         .reduce(
-            || DpaProcessor::new(size, guess_range, leakage_model),
+            || DpaProcessor::new(size, guess_range, selection_function),
             |x, y| x + y,
-        );
+        )
+        .finalize();
 
-    let dpa = dpa.finalize();
-    println!("{:2x}", dpa.pass_guess());
-    // let corr = dpa.pass_corr_array();
+    println!("{:2x}", dpa.best_guess());
+    // let corr = dpa.corr();
 
     Ok(())
 }
diff --git a/examples/rank.rs b/examples/rank.rs
index a8b3cfd..a8e9acc 100644
--- a/examples/rank.rs
+++ b/examples/rank.rs
@@ -39,7 +39,7 @@ fn rank() -> Result<()> {
                 .par_bridge()
                 .fold(
                     || CpaProcessor::new(size, guess_range, target_byte, leakage_model),
-                    |mut r: CpaProcessor, n| {
+                    |mut r, n| {
                         r.update(
                             l_sample.row(n).map(|l| *l as usize).view(),
                             p_sample.row(n).map(|p| *p as usize).view(),
@@ -58,7 +58,7 @@ fn rank() -> Result<()> {
     let rank = rank.finalize();
 
     // save rank key curves in npy
-    save_array("../results/rank.npy", &rank.pass_rank())?;
+    save_array("../results/rank.npy", &rank.rank().map(|&x| x as u64))?;
 
     Ok(())
 }
diff --git a/src/cpa.rs b/src/cpa.rs
index 1e67d78..5718d03 100644
--- a/src/cpa.rs
+++ b/src/cpa.rs
@@ -1,26 +1,27 @@
-use crate::util::max_per_row;
-use ndarray::{concatenate, s, Array1, Array2, ArrayView1, ArrayView2, Axis};
+use crate::util::{argmax_by, argsort_by, max_per_row};
+use ndarray::{s, Array1, Array2, ArrayView1, ArrayView2, Axis};
 use rayon::{
     iter::ParallelBridge,
     prelude::{IntoParallelIterator, ParallelIterator},
 };
 use std::{iter::zip, ops::Add};
 
-/// Computes the [`Cpa`] of the given traces using [`CpaProcessor`].
+/// Compute the [`Cpa`] of the given traces using [`CpaProcessor`].
 ///
 /// # Panics
 /// - Panic if `leakages.shape()[0] != plaintexts.shape()[0]`
 /// - Panic if `chunk_size` is 0.
-pub fn cpa<T>(
+pub fn cpa<T, F>(
     leakages: ArrayView2<T>,
     plaintexts: ArrayView2<T>,
     guess_range: usize,
     target_byte: usize,
-    leakage_func: fn(usize, usize) -> usize,
+    leakage_func: F,
     chunk_size: usize,
 ) -> Cpa
 where
     T: Into<usize> + Copy + Sync,
+    F: Fn(usize, usize) -> usize + Send + Sync + Copy,
 {
     assert_eq!(leakages.shape()[0], plaintexts.shape()[0]);
     assert!(chunk_size > 0);
@@ -46,43 +47,51 @@ where
     .finalize()
 }
 
+/// Result of the CPA[^1] on some traces.
+///
+/// [^1]: <https://www.iacr.org/archive/ches2004/31560016/31560016.pdf>
+#[derive(Debug)]
 pub struct Cpa {
-    /// Guess range upper excluded bound
-    guess_range: usize,
     /// Pearson correlation coefficients
-    corr: Array2<f32>,
-    /// Max pearson correlation coefficients
-    max_corr: Array1<f32>,
-    rank_slice: Array2<f32>,
+    pub(crate) corr: Array2<f32>,
 }
 
 impl Cpa {
-    pub fn pass_rank(&self) -> ArrayView2<f32> {
-        self.rank_slice.view()
+    /// Rank guesses.
+    pub fn rank(&self) -> Array1<usize> {
+        let rank = argsort_by(&self.max_corr().to_vec()[..], f32::total_cmp);
+
+        Array1::from_vec(rank)
     }
 
-    pub fn pass_corr_array(&self) -> ArrayView2<f32> {
+    /// Return the Pearson correlation coefficients.
+    pub fn corr(&self) -> ArrayView2<f32> {
         self.corr.view()
     }
 
-    pub fn pass_guess(&self) -> usize {
-        let mut init_value = 0.0;
-        let mut guess = 0;
-
-        for i in 0..self.guess_range {
-            if self.max_corr[i] > init_value {
-                init_value = self.max_corr[i];
-                guess = i;
-            }
-        }
+    /// Return the guess with the highest Pearson correlation coefficient.
+    pub fn best_guess(&self) -> usize {
+        argmax_by(self.max_corr().view(), f32::total_cmp)
+    }
 
-        guess
+    /// Return the maximum Pearson correlation coefficient for each guess.
+    pub fn max_corr(&self) -> Array1<f32> {
+        max_per_row(self.corr.view())
     }
 }
 
-pub struct CpaProcessor {
+/// A processor that computes the [`Cpa`] of the given traces.
+///
+/// It implements algorithm 4 from [^1].
+///
+/// [^1]: <https://eprint.iacr.org/2013/794.pdf>
+pub struct CpaProcessor<F>
+where
+    F: Fn(usize, usize) -> usize,
+{
     /// Number of samples per trace
-    len_samples: usize,
+    num_samples: usize,
+    /// Target byte index in a block
     target_byte: usize,
     /// Guess range upper excluded bound
     guess_range: usize,
@@ -94,49 +103,53 @@ pub struct CpaProcessor {
     guess_sum_leakages: Array1<usize>,
     /// Sum of square of traces per key guess
     guess_sum_squares_leakages: Array1<usize>,
-    a_l: Array2<usize>,
+    /// Sum of traces per plaintext used
+    /// See 4.3 in <https://eprint.iacr.org/2013/794.pdf>
+    plaintext_sum_leakages: Array2<usize>,
     /// Leakage model
-    leakage_func: fn(usize, usize) -> usize,
+    leakage_func: F,
     /// Number of traces processed
-    len_leakages: usize,
+    num_traces: usize,
 }
 
-/* This class implements the CPA shown in this paper: https://eprint.iacr.org/2013/794.pdf */
-impl CpaProcessor {
+impl<F> CpaProcessor<F>
+where
+    F: Fn(usize, usize) -> usize + Sync,
+{
     pub fn new(
-        size: usize,
+        num_samples: usize,
         guess_range: usize,
         target_byte: usize,
-        leakage_func: fn(usize, usize) -> usize,
+        leakage_func: F,
     ) -> Self {
         Self {
-            len_samples: size,
+            num_samples,
             target_byte,
             guess_range,
-            sum_leakages: Array1::zeros(size),
-            sum_squares_leakages: Array1::zeros(size),
+            sum_leakages: Array1::zeros(num_samples),
+            sum_squares_leakages: Array1::zeros(num_samples),
             guess_sum_leakages: Array1::zeros(guess_range),
             guess_sum_squares_leakages: Array1::zeros(guess_range),
-            a_l: Array2::zeros((guess_range, size)),
+            plaintext_sum_leakages: Array2::zeros((guess_range, num_samples)),
             leakage_func,
-            len_leakages: 0,
+            num_traces: 0,
         }
     }
 
     /// # Panics
-    /// Panic in debug if `trace.shape()[0] != self.len_samples`.
+    /// Panic in debug if `trace.shape()[0] != self.num_samples`.
     pub fn update<T>(&mut self, trace: ArrayView1<T>, plaintext: ArrayView1<T>)
     where
         T: Into<usize> + Copy,
     {
-        debug_assert_eq!(trace.shape()[0], self.len_samples);
+        debug_assert_eq!(trace.shape()[0], self.num_samples);
 
-        /* This function updates the main arrays of the CPA, as shown in Alg. 4
-        in the paper.*/
-
-        for i in 0..self.len_samples {
+        let partition = plaintext[self.target_byte].into();
+        for i in 0..self.num_samples {
             self.sum_leakages[i] += trace[i].into();
             self.sum_squares_leakages[i] += trace[i].into() * trace[i].into();
+
+            self.plaintext_sum_leakages[[partition, i]] += trace[i].into();
         }
 
         for guess in 0..self.guess_range {
@@ -145,91 +158,85 @@ impl CpaProcessor {
             self.guess_sum_squares_leakages[guess] += value * value;
         }
 
-        let partition = plaintext[self.target_byte].into();
-        for i in 0..self.len_samples {
-            self.a_l[[partition, i]] += trace[i].into();
-        }
-
-        self.len_leakages += 1;
+        self.num_traces += 1;
     }
 
+    /// Finalize the calculation after feeding the overall traces.
     pub fn finalize(&self) -> Cpa {
-        /* This function finalizes the calculation after feeding the
-        overall traces */
-        let mut p = Array2::zeros((self.guess_range, self.guess_range));
+        let mut modeled_leakages = Array1::zeros(self.guess_range);
+
+        let mut corr = Array2::zeros((self.guess_range, self.num_samples));
         for guess in 0..self.guess_range {
-            for x in 0..self.guess_range {
-                p[[x, guess]] = (self.leakage_func)(x, guess);
+            for u in 0..self.guess_range {
+                modeled_leakages[u] = (self.leakage_func)(u, guess);
             }
-        }
 
-        let mut corr = Array2::zeros((self.guess_range, self.len_samples));
-        for guess in 0..self.guess_range {
-            let mean_key = self.guess_sum_leakages[guess] as f32 / self.len_leakages as f32;
+            let mean_key = self.guess_sum_leakages[guess] as f32 / self.num_traces as f32;
             let mean_squares_key =
-                self.guess_sum_squares_leakages[guess] as f32 / self.len_leakages as f32;
+                self.guess_sum_squares_leakages[guess] as f32 / self.num_traces as f32;
             let var_key = mean_squares_key - (mean_key * mean_key);
 
-            /* Parallel operation using multi-threading */
-            let tmp: Vec<f32> = (0..self.len_samples)
+            let guess_corr: Vec<_> = (0..self.num_samples)
                 .into_par_iter()
-                .map(|x| {
-                    let mean_leakages = self.sum_leakages[x] as f32 / self.len_leakages as f32;
-                    let summult = self.sum_mult(self.a_l.slice(s![.., x]), p.slice(s![.., guess]));
-                    let upper1 = summult as f32 / self.len_leakages as f32;
-                    let upper = upper1 - (mean_key * mean_leakages);
+                .map(|u| {
+                    let mean_leakages = self.sum_leakages[u] as f32 / self.num_traces as f32;
+
+                    let cov = self.sum_mult(
+                        self.plaintext_sum_leakages.slice(s![.., u]),
+                        modeled_leakages.view(),
+                    );
+                    let cov = cov as f32 / self.num_traces as f32 - (mean_key * mean_leakages);
 
                     let mean_squares_leakages =
-                        self.sum_squares_leakages[x] as f32 / self.len_leakages as f32;
+                        self.sum_squares_leakages[u] as f32 / self.num_traces as f32;
                     let var_leakages = mean_squares_leakages - (mean_leakages * mean_leakages);
-                    let lower = f32::sqrt(var_key * var_leakages);
-
-                    f32::abs(upper / lower)
+                    f32::abs(cov / f32::sqrt(var_key * var_leakages))
                 })
                 .collect();
 
             #[allow(clippy::needless_range_loop)]
-            for u in 0..self.len_samples {
-                corr[[guess, u]] = tmp[u];
+            for u in 0..self.num_samples {
+                corr[[guess, u]] = guess_corr[u];
             }
         }
 
-        let max_corr = max_per_row(corr.view());
-
-        let mut rank_slice = Array2::zeros((self.guess_range, 1));
-        rank_slice = concatenate![
-            Axis(1),
-            rank_slice,
-            max_corr
-                .clone()
-                .into_shape((max_corr.shape()[0], 1))
-                .unwrap()
-        ];
+        Cpa { corr }
+    }
 
-        Cpa {
-            guess_range: self.guess_range,
-            corr,
-            max_corr,
-            rank_slice,
-        }
+    fn sum_mult(&self, a: ArrayView1<usize>, b: ArrayView1<usize>) -> usize {
+        a.dot(&b)
     }
 
-    fn sum_mult(&self, a: ArrayView1<usize>, b: ArrayView1<usize>) -> i32 {
-        a.dot(&b) as i32
+    /// Determine if two [`CpaProcessor`] are compatible for addition.
+    ///
+    /// If they were created with the same parameters, they are compatible.
+    ///
+    /// Note: [`CpaProcessor::leakage_func`] cannot be checked for equality, but they must have the
+    /// same leakage functions in order to be compatible.
+    fn is_compatible_with(&self, other: &Self) -> bool {
+        self.num_samples == other.num_samples
+            && self.target_byte == other.target_byte
+            && self.guess_range == other.guess_range
     }
 }
 
-impl Add for CpaProcessor {
+impl<F> Add for CpaProcessor<F>
+where
+    F: Fn(usize, usize) -> usize + Sync,
+{
     type Output = Self;
 
+    /// Merge computations of two [`CpaProcessor`]. Processors need to be compatible to be merged
+    /// together, otherwise it can panic or yield incoherent result (see
+    /// [`CpaProcessor::is_compatible_with`]).
+    ///
+    /// # Panics
+    /// Panics in debug if the processors are not compatible.
     fn add(self, rhs: Self) -> Self::Output {
-        debug_assert_eq!(self.target_byte, rhs.target_byte);
-        debug_assert_eq!(self.guess_range, rhs.guess_range);
-        debug_assert_eq!(self.len_samples, rhs.len_samples);
-        debug_assert_eq!(self.leakage_func, rhs.leakage_func);
+        debug_assert!(self.is_compatible_with(&rhs));
 
         Self {
-            len_samples: self.len_samples,
+            num_samples: self.num_samples,
             target_byte: self.target_byte,
             guess_range: self.guess_range,
             sum_leakages: self.sum_leakages + rhs.sum_leakages,
@@ -237,9 +244,9 @@ impl Add for CpaProcessor {
             guess_sum_leakages: self.guess_sum_leakages + rhs.guess_sum_leakages,
             guess_sum_squares_leakages: self.guess_sum_squares_leakages
                 + rhs.guess_sum_squares_leakages,
-            a_l: self.a_l + rhs.a_l,
+            plaintext_sum_leakages: self.plaintext_sum_leakages + rhs.plaintext_sum_leakages,
             leakage_func: self.leakage_func,
-            len_leakages: self.len_leakages + rhs.len_leakages,
+            num_traces: self.num_traces + rhs.num_traces,
         }
     }
 }
diff --git a/src/cpa_normal.rs b/src/cpa_normal.rs
index ee5aea2..3bee39a 100644
--- a/src/cpa_normal.rs
+++ b/src/cpa_normal.rs
@@ -1,24 +1,25 @@
-use ndarray::{concatenate, Array1, Array2, ArrayView1, ArrayView2, Axis};
+use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
 use rayon::iter::{ParallelBridge, ParallelIterator};
 use std::{iter::zip, ops::Add};
 
-use crate::util::max_per_row;
+use crate::cpa::Cpa;
 
-/// Computes the [`Cpa`] of the given traces using [`CpaProcessor`].
+/// Compute the [`Cpa`] of the given traces using [`CpaProcessor`].
 ///
 /// # Panics
 /// - Panic if `leakages.shape()[0] != plaintexts.shape()[0]`
 /// - Panic if `chunk_size` is 0.
-pub fn cpa<T, U>(
+pub fn cpa<T, U, F>(
     leakages: ArrayView2<T>,
     plaintexts: ArrayView2<U>,
     guess_range: usize,
-    leakage_func: fn(ArrayView1<usize>, usize) -> usize,
+    leakage_func: F,
     chunk_size: usize,
 ) -> Cpa
 where
     T: Into<f32> + Copy + Sync,
     U: Into<usize> + Copy + Sync,
+    F: Fn(ArrayView1<usize>, usize) -> usize + Send + Sync + Copy,
 {
     assert_eq!(leakages.shape()[0], plaintexts.shape()[0]);
     assert!(chunk_size > 0);
@@ -41,41 +42,15 @@ where
     .finalize()
 }
 
-pub struct Cpa {
-    /// Guess range upper excluded bound
-    guess_range: usize,
-    corr: Array2<f32>,
-    max_corr: Array1<f32>,
-    rank_slice: Array2<f32>,
-}
-
-impl Cpa {
-    pub fn pass_rank(&self) -> ArrayView2<f32> {
-        self.rank_slice.view()
-    }
-
-    pub fn pass_corr_array(&self) -> ArrayView2<f32> {
-        self.corr.view()
-    }
-
-    pub fn pass_guess(&self) -> usize {
-        let mut init_value = 0.0;
-        let mut guess = 0;
-
-        for i in 0..self.guess_range {
-            if self.max_corr[i] > init_value {
-                init_value = self.max_corr[i];
-                guess = i;
-            }
-        }
-
-        guess
-    }
-}
-
-pub struct CpaProcessor {
+/// A processor that computes the [`Cpa`] of the given traces.
+///
+/// [^1]: <https://www.iacr.org/archive/ches2004/31560016/31560016.pdf>
+pub struct CpaProcessor<F>
+where
+    F: Fn(ArrayView1<usize>, usize) -> usize,
+{
     /// Number of samples per trace
-    len_samples: usize,
+    num_samples: usize,
     /// Guess range upper excluded bound
     guess_range: usize,
     /// Sum of traces
@@ -88,61 +63,54 @@ pub struct CpaProcessor {
     guess_sum2_leakages: Array1<f32>,
     values: Array2<f32>,
     cov: Array2<f32>,
-    rank_traces: usize, // Number of traces to calculate succes rate
     /// Batch size
     batch_size: usize,
     /// Leakage model
-    leakage_func: fn(ArrayView1<usize>, usize) -> usize,
+    leakage_func: F,
     /// Number of traces processed
-    len_leakages: usize,
+    num_traces: usize,
 }
 
-/* This class implements the CPA algorithm shown in:
-https://www.iacr.org/archive/ches2004/31560016/31560016.pdf */
-
-impl CpaProcessor {
-    pub fn new(
-        size: usize,
-        batch_size: usize,
-        guess_range: usize,
-        leakage_func: fn(ArrayView1<usize>, usize) -> usize,
-    ) -> Self {
+impl<F> CpaProcessor<F>
+where
+    F: Fn(ArrayView1<usize>, usize) -> usize,
+{
+    pub fn new(num_samples: usize, batch_size: usize, guess_range: usize, leakage_func: F) -> Self {
         Self {
-            len_samples: size,
+            num_samples,
             guess_range,
-            sum_leakages: Array1::zeros(size),
-            sum2_leakages: Array1::zeros(size),
+            sum_leakages: Array1::zeros(num_samples),
+            sum2_leakages: Array1::zeros(num_samples),
             guess_sum_leakages: Array1::zeros(guess_range),
             guess_sum2_leakages: Array1::zeros(guess_range),
             values: Array2::zeros((batch_size, guess_range)),
-            cov: Array2::zeros((guess_range, size)),
-            rank_traces: 0,
+            cov: Array2::zeros((guess_range, num_samples)),
             batch_size,
             leakage_func,
-            len_leakages: 0,
+            num_traces: 0,
         }
     }
 
     /// # Panics
     /// - Panic in debug if `trace_batch.shape()[0] != plaintext_batch.shape()[0]`.
-    /// - Panic in debug if `trace_batch.shape()[1] != self.len_samples`.
+    /// - Panic in debug if `trace_batch.shape()[1] != self.num_samples`.
     pub fn update<T, U>(&mut self, trace_batch: ArrayView2<T>, plaintext_batch: ArrayView2<U>)
     where
         T: Into<f32> + Copy,
         U: Into<usize> + Copy,
     {
         debug_assert_eq!(trace_batch.shape()[0], plaintext_batch.shape()[0]);
-        debug_assert_eq!(trace_batch.shape()[1], self.len_samples);
+        debug_assert_eq!(trace_batch.shape()[1], self.num_samples);
 
         /* This function updates the internal arrays of the CPA
         It accepts trace_batch and plaintext_batch to update them*/
-        let trace_batch = trace_batch.map(|&t| t.into());
-        let plaintext_batch = plaintext_batch.map(|&m| m.into());
+        let trace_batch = trace_batch.mapv(|t| t.into());
+        let plaintext_batch = plaintext_batch.mapv(|m| m.into());
 
         self.update_values(trace_batch.view(), plaintext_batch.view(), self.guess_range);
         self.update_key_leakages(trace_batch.view(), self.guess_range);
 
-        self.len_leakages += self.batch_size;
+        self.num_traces += self.batch_size;
     }
 
     fn update_values(
@@ -163,7 +131,7 @@ impl CpaProcessor {
     }
 
     fn update_key_leakages(&mut self, trace: ArrayView2<f32>, guess_range: usize) {
-        for i in 0..self.len_samples {
+        for i in 0..self.num_samples {
             self.sum_leakages[i] += trace.column(i).sum(); // trace[i] as usize;
             self.sum2_leakages[i] += trace.column(i).dot(&trace.column(i)); // (trace[i] * trace[i]) as usize;
         }
@@ -176,59 +144,17 @@ impl CpaProcessor {
         }
     }
 
-    pub fn update_success<T, U>(
-        &mut self,
-        trace_batch: ArrayView2<T>,
-        plaintext_batch: ArrayView2<U>,
-    ) -> Option<Cpa>
-    where
-        T: Into<f32> + Copy,
-        U: Into<usize> + Copy,
-    {
-        /* This function updates the main arrays of the CPA for the success rate*/
-        self.update(trace_batch, plaintext_batch);
-
-        // WARN: if self.rank_traces == 0 this function will panic (division by zero)
-        // WARN: if self.rank_traces is not divisible by self.batch_size this branch will never be
-        // taken
-        if self.len_leakages % self.rank_traces == 0 {
-            let mut cpa = self.finalize();
-
-            if self.len_leakages == self.rank_traces {
-                cpa.rank_slice = cpa
-                    .max_corr
-                    .clone()
-                    .into_shape((cpa.max_corr.shape()[0], 1))
-                    .unwrap();
-            } else {
-                cpa.rank_slice = concatenate![
-                    Axis(1),
-                    cpa.rank_slice,
-                    cpa.max_corr
-                        .clone()
-                        .into_shape((cpa.max_corr.shape()[0], 1))
-                        .unwrap()
-                ];
-            }
-
-            Some(cpa)
-        } else {
-            None
-        }
-    }
-
+    /// Finalize the calculation after feeding the overall traces.
     pub fn finalize(&self) -> Cpa {
-        /* This function finalizes the calculation after
-        feeding all stored acc arrays */
-        let cov_n = self.cov.clone() / self.len_leakages as f32;
-        let avg_keys = self.guess_sum_leakages.clone() / self.len_leakages as f32;
-        let std_key = self.guess_sum2_leakages.clone() / self.len_leakages as f32;
-        let avg_leakages = self.sum_leakages.clone() / self.len_leakages as f32;
-        let std_leakages = self.sum2_leakages.clone() / self.len_leakages as f32;
+        let cov_n = self.cov.clone() / self.num_traces as f32;
+        let avg_keys = self.guess_sum_leakages.clone() / self.num_traces as f32;
+        let std_key = self.guess_sum2_leakages.clone() / self.num_traces as f32;
+        let avg_leakages = self.sum_leakages.clone() / self.num_traces as f32;
+        let std_leakages = self.sum2_leakages.clone() / self.num_traces as f32;
 
-        let mut corr = Array2::zeros((self.guess_range, self.len_samples));
+        let mut corr = Array2::zeros((self.guess_range, self.num_samples));
         for i in 0..self.guess_range {
-            for x in 0..self.len_samples {
+            for x in 0..self.num_samples {
                 let numerator = cov_n[[i, x]] - (avg_keys[i] * avg_leakages[x]);
 
                 let denominator_1 = std_key[i] - (avg_keys[i] * avg_keys[i]);
@@ -240,31 +166,39 @@ impl CpaProcessor {
             }
         }
 
-        let max_corr = max_per_row(corr.view());
-
-        Cpa {
-            guess_range: self.guess_range,
-            corr,
-            max_corr,
-            rank_slice: Array2::zeros((self.guess_range, 1)),
-        }
+        Cpa { corr }
     }
 
-    pub fn success_traces(&mut self, traces_no: usize) {
-        self.rank_traces = traces_no;
+    /// Determine if two [`CpaProcessor`] are compatible for addition.
+    ///
+    /// If they were created with the same parameters, they are compatible.
+    ///
+    /// Note: [`CpaProcessor::leakage_func`] cannot be checked for equality, but they must have the
+    /// same leakage functions in order to be compatible.
+    fn is_compatible_with(&self, other: &Self) -> bool {
+        self.num_samples == other.num_samples
+            && self.batch_size == other.batch_size
+            && self.guess_range == other.guess_range
     }
 }
 
-impl Add for CpaProcessor {
+impl<F> Add for CpaProcessor<F>
+where
+    F: Fn(ArrayView1<usize>, usize) -> usize,
+{
     type Output = Self;
 
+    /// Merge computations of two [`CpaProcessor`]. Processors need to be compatible to be merged
+    /// together, otherwise it can panic or yield incoherent result (see
+    /// [`CpaProcessor::is_compatible_with`]).
+    ///
+    /// # Panics
+    /// Panics in debug if the processors are not compatible.
     fn add(self, rhs: Self) -> Self::Output {
-        debug_assert_eq!(self.len_samples, rhs.len_samples);
-        debug_assert_eq!(self.batch_size, rhs.batch_size);
-        debug_assert_eq!(self.guess_range, rhs.guess_range);
+        debug_assert!(self.is_compatible_with(&rhs));
 
         Self {
-            len_samples: self.len_samples,
+            num_samples: self.num_samples,
             guess_range: self.guess_range,
             sum_leakages: self.sum_leakages + rhs.sum_leakages,
             sum2_leakages: self.sum2_leakages + rhs.sum2_leakages,
@@ -272,10 +206,9 @@ impl Add for CpaProcessor {
             guess_sum2_leakages: self.guess_sum2_leakages + rhs.guess_sum2_leakages,
             values: self.values + rhs.values,
             cov: self.cov + rhs.cov,
-            rank_traces: self.rank_traces,
             batch_size: self.batch_size,
             leakage_func: self.leakage_func,
-            len_leakages: self.len_leakages + rhs.len_leakages,
+            num_traces: self.num_traces + rhs.num_traces,
         }
     }
 }
diff --git a/src/dpa.rs b/src/dpa.rs
index 41083ab..cc06ee3 100644
--- a/src/dpa.rs
+++ b/src/dpa.rs
@@ -1,21 +1,24 @@
-use ndarray::{concatenate, Array1, Array2, ArrayView1, ArrayView2, Axis};
+use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
 use rayon::iter::{ParallelBridge, ParallelIterator};
-use std::{iter::zip, ops::Add};
+use std::{iter::zip, marker::PhantomData, ops::Add};
 
-use crate::util::max_per_row;
+use crate::util::{argmax_by, argsort_by, max_per_row};
 
+/// Compute the [`Dpa`] of the given traces using [`DpaProcessor`].
+///
 /// # Panics
-/// Panics if `chunk_size` is not strictly positive.
-pub fn dpa<M, T>(
+/// Panic if `chunk_size` is not strictly positive.
+pub fn dpa<M, T, F>(
     leakages: ArrayView2<T>,
     metadata: ArrayView1<M>,
     guess_range: usize,
-    leakage_func: fn(M, usize) -> usize,
+    selection_function: F,
     chunk_size: usize,
 ) -> Dpa
 where
     T: Into<f32> + Copy + Sync,
-    M: Clone + Sync,
+    M: Clone + Send + Sync,
+    F: Fn(M, usize) -> bool + Send + Sync + Copy,
 {
     assert!(chunk_size > 0);
 
@@ -25,7 +28,7 @@ where
     )
     .par_bridge()
     .fold(
-        || DpaProcessor::new(leakages.shape()[1], guess_range, leakage_func),
+        || DpaProcessor::new(leakages.shape()[1], guess_range, selection_function),
         |mut dpa, (leakages_chunk, metadata_chunk)| {
             for i in 0..leakages_chunk.shape()[0] {
                 dpa.update(leakages_chunk.row(i), metadata_chunk[i].clone());
@@ -39,190 +42,162 @@ where
     .finalize()
 }
 
+/// Result of the DPA[^1] on some traces.
+///
+/// [^1]: <https://paulkocher.com/doc/DifferentialPowerAnalysis.pdf>
+#[derive(Debug)]
 pub struct Dpa {
-    /// Guess range upper excluded bound
-    guess_range: usize,
-    corr: Array2<f32>,
-    max_corr: Array1<f32>,
+    differential_curves: Array2<f32>,
 }
 
 impl Dpa {
-    pub fn pass_corr_array(&self) -> ArrayView2<f32> {
-        self.corr.view()
+    /// Return the rank of guesses
+    pub fn rank(&self) -> Array1<usize> {
+        let rank = argsort_by(&self.max_differential_curves().to_vec()[..], f32::total_cmp);
+
+        Array1::from_vec(rank)
     }
 
-    pub fn pass_guess(&self) -> usize {
-        let mut init_value = 0.0;
-        let mut guess = 0;
+    /// Return the differential curves
+    pub fn differential_curves(&self) -> ArrayView2<f32> {
+        self.differential_curves.view()
+    }
 
-        for i in 0..self.guess_range {
-            if self.max_corr[i] > init_value {
-                init_value = self.max_corr[i];
-                guess = i;
-            }
-        }
+    /// Return the guess with the highest differential peak.
+    pub fn best_guess(&self) -> usize {
+        argmax_by(self.max_differential_curves().view(), f32::total_cmp)
+    }
 
-        guess
+    /// Return the maximum differential peak for each guess.
+    pub fn max_differential_curves(&self) -> Array1<f32> {
+        max_per_row(self.differential_curves.view())
     }
 }
 
-pub struct DpaProcessor<M> {
+/// A processor that computes the [`Dpa`] of the given traces.
+///
+/// [^1]: <https://paulkocher.com/doc/DifferentialPowerAnalysis.pdf>
+/// [^2]: <https://web.mit.edu/6.857/OldStuff/Fall03/ref/kocher-DPATechInfo.pdf>
+pub struct DpaProcessor<M, F>
+where
+    F: Fn(M, usize) -> bool,
+{
     /// Number of samples per trace
-    len_samples: usize,
+    num_samples: usize,
     /// Guess range upper excluded bound
     guess_range: usize,
-    /// Sum of traces for which the selection function equals 0
+    /// Sum of traces for which the selection function equals false
     sum_0: Array2<f32>,
-    /// Sum of traces for which the selection function equals 1
+    /// Sum of traces for which the selection function equals true
     sum_1: Array2<f32>,
-    /// Number of traces processed for which the selection function equals 0
+    /// Number of traces processed for which the selection function equals false
     count_0: Array1<usize>,
-    /// Number of traces processed for which the selection function equals 1
+    /// Number of traces processed for which the selection function equals true
     count_1: Array1<usize>,
-    rank_slice: Array2<f32>,
-    rank_traces: usize, // Number of traces to calculate succes rate
-    /// Selection function
-    leakage_func: fn(M, usize) -> usize,
+    selection_function: F,
     /// Number of traces processed
-    len_leakages: usize,
+    num_traces: usize,
+    _metadata: PhantomData<M>,
 }
 
-/* This class implements the DPA algorithm shown in:
-https://paulkocher.com/doc/DifferentialPowerAnalysis.pdf
-https://web.mit.edu/6.857/OldStuff/Fall03/ref/kocher-DPATechInfo.pdf */
-
-impl<M: Clone> DpaProcessor<M> {
-    pub fn new(size: usize, guess_range: usize, f: fn(M, usize) -> usize) -> Self {
+impl<M, F> DpaProcessor<M, F>
+where
+    M: Clone,
+    F: Fn(M, usize) -> bool,
+{
+    pub fn new(num_samples: usize, guess_range: usize, selection_function: F) -> Self {
         Self {
-            len_samples: size,
+            num_samples,
             guess_range,
-            sum_0: Array2::zeros((guess_range, size)),
-            sum_1: Array2::zeros((guess_range, size)),
+            sum_0: Array2::zeros((guess_range, num_samples)),
+            sum_1: Array2::zeros((guess_range, num_samples)),
             count_0: Array1::zeros(guess_range),
             count_1: Array1::zeros(guess_range),
-            rank_slice: Array2::zeros((guess_range, 1)),
-            rank_traces: 0,
-            leakage_func: f,
-            len_leakages: 0,
+            selection_function,
+            num_traces: 0,
+            _metadata: PhantomData,
         }
     }
 
     /// # Panics
-    /// Panic in debug if `trace.shape()[0] != self.len_samples`.
+    /// Panic in debug if `trace.shape()[0] != self.num_samples`.
     pub fn update<T>(&mut self, trace: ArrayView1<T>, metadata: M)
     where
         T: Into<f32> + Copy,
     {
-        debug_assert_eq!(trace.shape()[0], self.len_samples);
+        debug_assert_eq!(trace.shape()[0], self.num_samples);
 
-        /* This function updates the internal arrays of the DPA
-        It accepts trace_batch and plaintext_batch to update them*/
         for guess in 0..self.guess_range {
-            let index = (self.leakage_func)(metadata.clone(), guess);
-            if index & 1 == 1 {
-                // classification is performed based on the lsb
-                // let tmp_row: Array1<f32> = self.sum_1.row(guess).to_owned() + tmp_trace.clone();
-                // self.sum_1.row_mut(guess).assign(&tmp_row);
-                for i in 0..self.len_samples {
+            if (self.selection_function)(metadata.clone(), guess) {
+                for i in 0..self.num_samples {
                     self.sum_1[[guess, i]] += trace[i].into();
                 }
                 self.count_1[guess] += 1;
             } else {
-                // let tmp_row: Array1<f32> = self.sum_0.row(guess).to_owned() + tmp_trace.clone();
-                // self.sum_0.row_mut(guess).assign(&tmp_row);
-
-                for i in 0..self.len_samples {
+                for i in 0..self.num_samples {
                     self.sum_0[[guess, i]] += trace[i].into();
                 }
                 self.count_0[guess] += 1;
             }
         }
-        self.len_leakages += 1;
-    }
 
-    pub fn update_success<T>(&mut self, trace_batch: ArrayView1<T>, plaintext_batch: M)
-    where
-        T: Into<f32> + Copy,
-    {
-        /* This function updates the main arrays of the DPA for the success rate*/
-        self.update(trace_batch, plaintext_batch);
-
-        if self.len_leakages % self.rank_traces == 0 {
-            let dpa = self.finalize();
-
-            if self.len_leakages == self.rank_traces {
-                self.rank_slice = dpa
-                    .max_corr
-                    .clone()
-                    .into_shape((dpa.max_corr.shape()[0], 1))
-                    .unwrap();
-            } else {
-                self.rank_slice = concatenate![
-                    Axis(1),
-                    self.rank_slice,
-                    dpa.max_corr
-                        .clone()
-                        .into_shape((dpa.max_corr.shape()[0], 1))
-                        .unwrap()
-                ];
-            }
-        }
+        self.num_traces += 1;
     }
 
-    pub fn assign_rank_traces(&mut self, value: usize) {
-        self.rank_traces = value;
-    }
-
-    pub fn finalize(&mut self) -> Dpa {
-        /* This function finalizes the calculation after feeding all stored acc arrays */
-        let mut tmp_avg_0 = Array2::zeros((self.guess_range, self.len_samples));
-        let mut tmp_avg_1 = Array2::zeros((self.guess_range, self.len_samples));
+    /// Finalizes the calculation after feeding the overall traces.
+    pub fn finalize(&self) -> Dpa {
+        let mut differential_curves = Array2::zeros((self.guess_range, self.num_samples));
+        for guess in 0..self.guess_range {
+            for i in 0..self.num_samples {
+                let mean_0 = self.sum_0[[guess, i]] / self.count_0[guess] as f32;
+                let mean_1 = self.sum_1[[guess, i]] / self.count_1[guess] as f32;
 
-        for row in 0..self.guess_range {
-            let tmp_row_0 = self.sum_0.row(row).to_owned() / self.count_0[row] as f32;
-            let tmp_row_1 = self.sum_1.row(row).to_owned() / self.count_1[row] as f32;
-            tmp_avg_0.row_mut(row).assign(&tmp_row_0);
-            tmp_avg_1.row_mut(row).assign(&tmp_row_1);
+                differential_curves[[guess, i]] = f32::abs(mean_0 - mean_1);
+            }
         }
 
-        let corr = (tmp_avg_0 - tmp_avg_1).map(|e| f32::abs(*e));
-        let max_corr = max_per_row(corr.view());
-
         Dpa {
-            guess_range: self.guess_range,
-            corr,
-            max_corr,
+            differential_curves,
         }
     }
 
-    pub fn success_traces(&mut self, traces_no: usize) {
-        self.rank_traces = traces_no;
-    }
-
-    pub fn pass_rank(&self) -> ArrayView2<f32> {
-        self.rank_slice.view()
+    /// Determine if two [`DpaProcessor`] are compatible for addition.
+    ///
+    /// If they were created with the same parameters, they are compatible.
+    ///
+    /// Note: [`DpaProcessor::selection_function`] cannot be checked for equality, but they must
+    /// have the same selection functions in order to be compatible.
+    fn is_compatible_with(&self, other: &Self) -> bool {
+        self.num_samples == other.num_samples && self.guess_range == other.guess_range
     }
 }
 
-impl<M> Add for DpaProcessor<M> {
+impl<M, F> Add for DpaProcessor<M, F>
+where
+    F: Fn(M, usize) -> bool,
+    M: Clone,
+{
     type Output = Self;
 
+    /// Merge computations of two [`DpaProcessor`]. Processors need to be compatible to be merged
+    /// together, otherwise it can panic or yield incoherent result (see
+    /// [`DpaProcessor::is_compatible_with`]).
+    ///
+    /// # Panics
+    /// Panics in debug if the processors are not compatible.
     fn add(self, rhs: Self) -> Self::Output {
-        debug_assert_eq!(self.len_samples, rhs.len_samples);
-        debug_assert_eq!(self.guess_range, rhs.guess_range);
-        debug_assert_eq!(self.leakage_func, rhs.leakage_func);
+        debug_assert!(self.is_compatible_with(&rhs));
 
         Self {
-            len_samples: self.len_samples,
+            num_samples: self.num_samples,
             guess_range: self.guess_range,
             sum_0: self.sum_0 + rhs.sum_0,
             sum_1: self.sum_1 + rhs.sum_1,
             count_0: self.count_0 + rhs.count_0,
             count_1: self.count_1 + rhs.count_1,
-            rank_slice: self.rank_slice,
-            rank_traces: self.rank_traces,
-            leakage_func: self.leakage_func,
-            len_leakages: self.len_leakages + rhs.len_leakages,
+            selection_function: self.selection_function,
+            num_traces: self.num_traces + rhs.num_traces,
+            _metadata: PhantomData,
         }
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index b33811c..505c3aa 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -4,6 +4,8 @@ pub mod dpa;
 pub mod leakage;
 pub mod preprocessors;
 pub mod processors;
-pub mod quicklog;
 pub mod trace;
 pub mod util;
+
+#[cfg(feature = "quicklog")]
+pub mod quicklog;
diff --git a/src/preprocessors.rs b/src/preprocessors.rs
index 0ba3e42..e6fd7cc 100644
--- a/src/preprocessors.rs
+++ b/src/preprocessors.rs
@@ -6,6 +6,7 @@ use crate::processors::MeanVar;
 
 /// Computes the centered product of "order" leakage samples
 /// Used particularly when performing high-order SCA
+#[derive(Debug)]
 pub struct CenteredProduct {
     /// Sum of traces
     acc: Array1<i64>,
@@ -49,14 +50,14 @@ impl CenteredProduct {
     /// Compute the mean
     pub fn finalize(&mut self) {
         if self.count != 0 {
-            self.mean = self.acc.map(|&x| x as f64 / self.count as f64)
+            self.mean = self.acc.mapv(|x| x as f64 / self.count as f64)
         }
         self.processed = true
     }
 
     /// Apply the processing to an input trace
     /// The centered product substract the mean of the traces and then perform products between every input time samples
-    pub fn apply<T: Into<f64> + Copy>(&mut self, trace: ArrayView1<T>) -> Array1<f64> {
+    pub fn apply<T: Into<f64> + Copy>(&self, trace: ArrayView1<T>) -> Array1<f64> {
         // First we substract the mean trace
         let centered_trace: Array1<f64> = trace.mapv(|x| x.into()) - &self.mean;
         let length_out_trace: usize = self.intervals.iter().map(|x| x.len()).product();
@@ -79,6 +80,7 @@ impl CenteredProduct {
 }
 
 /// Elevates parts of a trace to a certain power
+#[derive(Debug)]
 pub struct Power {
     intervals: Vec<Range<i32>>,
     power: i32,
@@ -110,6 +112,7 @@ impl Power {
 }
 
 /// Standardization of the traces by removing the mean and scaling to unit variance
+#[derive(Debug)]
 pub struct StandardScaler {
     /// meanVar processor
     meanvar: MeanVar,
@@ -140,7 +143,7 @@ impl StandardScaler {
     }
 
     /// Apply the processing to an input trace
-    pub fn apply<T: Into<f64> + Copy>(&mut self, trace: ArrayView1<T>) -> Array1<f64> {
+    pub fn apply<T: Into<f64> + Copy>(&self, trace: ArrayView1<T>) -> Array1<f64> {
         (trace.mapv(|x| x.into()) - &self.mean) / &self.std
     }
 }
@@ -200,7 +203,7 @@ mod tests {
 
         for (i, t) in traces.iter().enumerate() {
             assert_eq!(
-                processor2.apply(t.view()).map(|x| round_to_2_digits(*x)),
+                processor2.apply(t.view()).mapv(round_to_2_digits),
                 expected_results[i]
             );
         }
@@ -228,15 +231,15 @@ mod tests {
         ];
 
         assert_eq!(
-            processor1.process(t.view()).map(|x| round_to_2_digits(*x)),
+            processor1.process(t.view()).mapv(round_to_2_digits),
             expected_results[0]
         );
         assert_eq!(
-            processor2.process(t.view()).map(|x| round_to_2_digits(*x)),
+            processor2.process(t.view()).mapv(round_to_2_digits),
             expected_results[1]
         );
         assert_eq!(
-            processor3.process(t.view()).map(|x| round_to_2_digits(*x)),
+            processor3.process(t.view()).mapv(round_to_2_digits),
             expected_results[2]
         );
     }
@@ -277,7 +280,7 @@ mod tests {
 
         for (i, t) in traces.iter().enumerate() {
             assert_eq!(
-                processor.apply(t.view()).map(|x| round_to_2_digits(*x)),
+                processor.apply(t.view()).mapv(round_to_2_digits),
                 expected_results[i]
             );
         }
diff --git a/src/processors.rs b/src/processors.rs
index 6b75829..96d1415 100644
--- a/src/processors.rs
+++ b/src/processors.rs
@@ -4,7 +4,7 @@ use rayon::iter::{ParallelBridge, ParallelIterator};
 use std::{iter::zip, ops::Add};
 
 /// Processes traces to calculate mean and variance.
-#[derive(Clone)]
+#[derive(Debug, Clone)]
 pub struct MeanVar {
     /// Sum of traces
     sum: Array1<i64>,
@@ -33,7 +33,7 @@ impl MeanVar {
     /// # Panics
     /// Panics in debug if the length of the trace is different form the size of [`MeanVar`].
     pub fn process<T: Into<i64> + Copy>(&mut self, trace: ArrayView1<T>) {
-        debug_assert!(trace.len() == self.sum.len());
+        debug_assert!(trace.len() == self.size());
 
         for i in 0..self.sum.len() {
             let x = trace[i].into();
@@ -49,7 +49,7 @@ impl MeanVar {
     pub fn mean(&self) -> Array1<f64> {
         let count = self.count as f64;
 
-        self.sum.map(|&x| x as f64 / count)
+        self.sum.mapv(|x| x as f64 / count)
     }
 
     /// Calculates and returns traces variance.
@@ -61,16 +61,36 @@ impl MeanVar {
             .collect()
     }
 
+    /// Returns the trace size handled.
+    pub fn size(&self) -> usize {
+        self.sum.len()
+    }
+
     /// Returns the number of traces processed.
     pub fn count(&self) -> usize {
         self.count
     }
+
+    /// Determine if two [`MeanVar`] are compatible for addition.
+    ///
+    /// If they were created with the same parameters, they are compatible.
+    fn is_compatible_with(&self, other: &Self) -> bool {
+        self.size() == other.size()
+    }
 }
 
 impl Add for MeanVar {
     type Output = Self;
 
+    /// Merge computations of two [`MeanVar`]. Processors need to be compatible to be merged
+    /// together, otherwise it can panic or yield incoherent result (see
+    /// [`MeanVar::is_compatible_with`]).
+    ///
+    /// # Panics
+    /// Panics in debug if the processors are not compatible.
     fn add(self, rhs: Self) -> Self::Output {
+        debug_assert!(self.is_compatible_with(&rhs));
+
         Self {
             sum: self.sum + rhs.sum,
             sum_squares: self.sum_squares + rhs.sum_squares,
@@ -117,7 +137,7 @@ where
 }
 
 /// Processes traces to calculate the Signal-to-Noise Ratio.
-#[derive(Clone)]
+#[derive(Debug, Clone)]
 pub struct Snr {
     mean_var: MeanVar,
     /// Sum of traces per class
@@ -133,11 +153,11 @@ impl Snr {
     ///
     /// * `size` - Size of the input traces
     /// * `classes` - Number of classes
-    pub fn new(size: usize, classes: usize) -> Self {
+    pub fn new(size: usize, num_classes: usize) -> Self {
         Self {
             mean_var: MeanVar::new(size),
-            classes_sum: Array2::zeros((classes, size)),
-            classes_count: Array1::zeros(classes),
+            classes_sum: Array2::zeros((num_classes, size)),
+            classes_count: Array1::zeros(num_classes),
         }
     }
 
@@ -146,11 +166,11 @@ impl Snr {
     /// # Panics
     /// Panics in debug if the length of the trace is different from the size of [`Snr`].
     pub fn process<T: Into<i64> + Copy>(&mut self, trace: ArrayView1<T>, class: usize) {
-        debug_assert!(trace.len() == self.classes_sum.shape()[1]);
+        debug_assert!(trace.len() == self.size());
 
         self.mean_var.process(trace);
 
-        for i in 0..self.classes_sum.shape()[1] {
+        for i in 0..self.size() {
             self.classes_sum[[class, i]] += trace[i].into();
         }
 
@@ -160,11 +180,11 @@ impl Snr {
     /// Returns Signal-to-Noise Ratio of the traces.
     /// SNR = V[E[L|X]] / E[V[L|X]]
     pub fn snr(&self) -> Array1<f64> {
-        let size = self.classes_sum.shape()[1];
-        let classes = self.classes_sum.shape()[0];
+        let size = self.size();
+        let num_classes = self.num_classes();
 
         let mut acc: Array1<f64> = Array1::zeros(size);
-        for class in 0..classes {
+        for class in 0..num_classes {
             if self.classes_count[class] == 0 {
                 continue;
             }
@@ -181,12 +201,37 @@ impl Snr {
         let velx = (acc / self.mean_var.count as f64) - mean.mapv(|x| x.powi(2));
         1f64 / (var / velx - 1f64)
     }
+
+    /// Returns the trace size handled
+    pub fn size(&self) -> usize {
+        self.classes_sum.shape()[1]
+    }
+
+    /// Returns the number of classes handled.
+    pub fn num_classes(&self) -> usize {
+        self.classes_count.len()
+    }
+
+    /// Determine if two [`Snr`] are compatible for addition.
+    ///
+    /// If they were created with the same parameters, they are compatible.
+    fn is_compatible_with(&self, other: &Self) -> bool {
+        self.size() == other.size() && self.num_classes() == other.num_classes()
+    }
 }
 
 impl Add for Snr {
     type Output = Self;
 
+    /// Merge computations of two [`Snr`]. Processors need to be compatible to be merged
+    /// together, otherwise it can panic or yield incoherent result (see
+    /// [`Snr::is_compatible_with`]).
+    ///
+    /// # Panics
+    /// Panics in debug if the processors are not compatible.
     fn add(self, rhs: Self) -> Self::Output {
+        debug_assert!(self.is_compatible_with(&rhs));
+
         Self {
             mean_var: self.mean_var + rhs.mean_var,
             classes_sum: self.classes_sum + rhs.classes_sum,
@@ -196,6 +241,7 @@ impl Add for Snr {
 }
 
 /// Processes traces to calculate Welch's T-Test.
+#[derive(Debug)]
 pub struct TTest {
     mean_var_1: MeanVar,
     mean_var_2: MeanVar,
@@ -222,9 +268,9 @@ impl TTest {
     /// * `class` - Indicates to which of the two partitions the given trace belongs.
     ///
     /// # Panics
-    /// Panics in debug if `trace.len() != self.mean_var_1.sum.len()`.
+    /// Panics in debug if `trace.len() != self.size()`.
     pub fn process<T: Into<i64> + Copy>(&mut self, trace: ArrayView1<T>, class: bool) {
-        debug_assert!(trace.len() == self.mean_var_1.sum.len());
+        debug_assert!(trace.len() == self.size());
 
         if class {
             self.mean_var_2.process(trace);
@@ -244,12 +290,32 @@ impl TTest {
             .mapv(f64::sqrt);
         q / d
     }
+
+    /// Returns the trace size handled.
+    pub fn size(&self) -> usize {
+        self.mean_var_1.size()
+    }
+
+    /// Determine if two [`TTest`] are compatible for addition.
+    ///
+    /// If they were created with the same parameters, they are compatible.
+    fn is_compatible_with(&self, other: &Self) -> bool {
+        self.size() == other.size()
+    }
 }
 
 impl Add for TTest {
     type Output = Self;
 
+    /// Merge computations of two [`TTest`]. Processors need to be compatible to be merged
+    /// together, otherwise it can panic or yield incoherent result (see
+    /// [`TTest::is_compatible_with`]).
+    ///
+    /// # Panics
+    /// Panics in debug if the processors are not compatible.
     fn add(self, rhs: Self) -> Self::Output {
+        debug_assert!(self.is_compatible_with(&rhs));
+
         Self {
             mean_var_1: self.mean_var_1 + rhs.mean_var_1,
             mean_var_2: self.mean_var_2 + rhs.mean_var_2,
diff --git a/src/util.rs b/src/util.rs
index 1736745..04c9670 100644
--- a/src/util.rs
+++ b/src/util.rs
@@ -1,9 +1,11 @@
 //! Convenient utility functions.
 
-use std::{io::Read, path::Path};
+use std::{cmp::Ordering, io::Read, path::Path};
 
-use ndarray::{Array, Array1, Array2, ArrayView2, Axis};
-use ndarray_npy::{read_npy, write_npy, ReadNpyError, ReadableElement, WriteNpyError};
+use ndarray::{Array, Array1, Array2, ArrayView1, ArrayView2, Axis};
+use ndarray_npy::{
+    read_npy, write_npy, ReadNpyError, ReadableElement, WritableElement, WriteNpyError,
+};
 use npyz::{Deserialize, NpyFile};
 
 #[cfg(feature = "progress_bar")]
@@ -58,7 +60,10 @@ pub fn read_array2_from_npy_file<T: ReadableElement>(
     read_npy(path)
 }
 
-pub fn save_array2(path: impl AsRef<Path>, array: ArrayView2<f32>) -> Result<(), WriteNpyError> {
+pub fn save_array2<T>(path: impl AsRef<Path>, array: ArrayView2<T>) -> Result<(), WriteNpyError>
+where
+    T: WritableElement,
+{
     write_npy(path, &array)
 }
 
@@ -78,3 +83,31 @@ pub fn max_per_row(arr: ArrayView2<f32>) -> Array1<f32> {
         })
         .collect()
 }
+
+/// Return the indices that would sort the given array with a comparison function.
+pub fn argsort_by<T, F>(data: &[T], compare: F) -> Vec<usize>
+where
+    F: Fn(&T, &T) -> Ordering,
+{
+    let mut indices: Vec<usize> = (0..data.len()).collect();
+
+    indices.sort_by(|&a, &b| compare(&data[a], &data[b]));
+
+    indices
+}
+
+/// Return the index of the maximum value in the given array.
+pub fn argmax_by<T, F>(array: ArrayView1<T>, compare: F) -> usize
+where
+    F: Fn(&T, &T) -> Ordering,
+{
+    let mut idx_max = 0;
+
+    for i in 0..array.len() {
+        if compare(&array[i], &array[idx_max]).is_gt() {
+            idx_max = i;
+        }
+    }
+
+    idx_max
+}