diff --git a/Cargo.lock b/Cargo.lock
index 975e855..89b4412 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -690,7 +690,7 @@ dependencies = [
 
 [[package]]
 name = "pymoors"
-version = "0.1.1-rc1"
+version = "0.1.1"
 dependencies = [
  "criterion",
  "ndarray 0.16.1",
diff --git a/Cargo.toml b/Cargo.toml
index 9a374ee..1816c1b 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "pymoors"
-version = "0.1.1-rc1"
+version = "0.1.1"
 edition = "2021"
 
 [lib]
diff --git a/python/pymoors/_pymoors.pyi b/python/pymoors/_pymoors.pyi
index 6116931..16ecd5a 100644
--- a/python/pymoors/_pymoors.pyi
+++ b/python/pymoors/_pymoors.pyi
@@ -325,6 +325,44 @@ class Nsga3:
 
     def run(self) -> None: ...
 
+class _RNsg2Kwargs(_MooAlgorithmKwargs, total=False):
+    reference_points: TwoDArray
+    epsilon: float
+
+class RNsga2:
+    """
+    Implementation of RNsga2 (Reference-based NSGA-II).
+
+    RNsga2 is a variant of NSGA-II that incorporates reference points into the selection process
+    to enhance diversity in many-objective optimization problems. By integrating a reference
+    point–based ranking into the survival operator, RNsga2 aims to obtain a well-distributed
+    approximation of the Pareto front even in high-dimensional objective spaces.
+
+    References:
+        Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
+            A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
+            IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
+
+        Deb, K., & Jain, H. (2014).
+            An Evolutionary Many-Objective Optimization Algorithm Using Reference-Point-Based
+            Nondominated Sorting Approach, Part I: Solving Problems with Box Constraints.
+            IEEE Transactions on Evolutionary Computation, 18(4), 577–601.
+
+    """
+
+    def __init__(self, **kwargs: Unpack[_RNsg2Kwargs]) -> None: ...
+    @property
+    def population(self) -> Population:
+        """
+        Returns the current population of individuals.
+
+        Returns:
+            Population: The current population.
+        """
+        ...
+
+    def run(self) -> None: ...
+
 # Custom Errors
 
 class NoFeasibleIndividualsError(BaseException):
diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index 434fbb8..c7a5441 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -23,6 +23,7 @@ mod macros;
 pub mod nsga2;
 pub mod nsga3;
 pub mod py_errors;
+pub mod rnsga2;
 
 #[derive(Debug)]
 pub enum MultiObjectiveAlgorithmError {
@@ -167,7 +168,7 @@ impl MultiObjectiveAlgorithm {
         )
         .expect("Failed to concatenate current population genes with offspring genes");
         // Build fronts from the combined genes.
-        let fronts = self.evaluator.build_fronts(combined_genes);
+        let mut fronts = self.evaluator.build_fronts(combined_genes);
 
         // Check if there are no feasible individuals
         if fronts.is_empty() {
@@ -175,7 +176,7 @@ impl MultiObjectiveAlgorithm {
         }
 
         // Select the new population
-        self.population = self.survivor.operate(&fronts, self.pop_size);
+        self.population = self.survivor.operate(&mut fronts, self.pop_size);
         Ok(())
     }
 
diff --git a/src/algorithms/nsga2.rs b/src/algorithms/nsga2.rs
index 28f6d88..14ff3c5 100644
--- a/src/algorithms/nsga2.rs
+++ b/src/algorithms/nsga2.rs
@@ -1,8 +1,6 @@
-use crate::define_multiobj_pyclass;
-use crate::operators::selection::RankAndCrowdingSelection;
-use crate::operators::survival::RankCrowdingSurvival;
 use pyo3::prelude::*;
 
+use crate::define_multiobj_pyclass;
 use crate::helpers::functions::{
     create_population_constraints_closure, create_population_fitness_closure,
 };
@@ -10,6 +8,8 @@ use crate::helpers::parser::{
     unwrap_crossover_operator, unwrap_duplicates_cleaner, unwrap_mutation_operator,
     unwrap_sampling_operator,
 };
+use crate::operators::selection::RankAndCrowdingSelection;
+use crate::operators::survival::RankCrowdingSurvival;
 
 // Define the NSGA-II algorithm using the macro
 define_multiobj_pyclass!(Nsga2, PyNsga2, "Nsga2");
diff --git a/src/algorithms/rnsga2.rs b/src/algorithms/rnsga2.rs
new file mode 100644
index 0000000..0fe3086
--- /dev/null
+++ b/src/algorithms/rnsga2.rs
@@ -0,0 +1,121 @@
+use pyo3::prelude::*;
+
+use crate::define_multiobj_pyclass;
+use crate::helpers::functions::{
+    create_population_constraints_closure, create_population_fitness_closure,
+};
+use crate::helpers::parser::{
+    unwrap_crossover_operator, unwrap_duplicates_cleaner, unwrap_mutation_operator,
+    unwrap_sampling_operator,
+};
+use crate::operators::selection::{DiversityComparison, RankAndCrowdingSelection};
+use crate::operators::survival::RankReferencePointsSurvival;
+
+use numpy::{PyArray2, PyArrayMethods};
+
+// Define the NSGA-II algorithm using the macro
+define_multiobj_pyclass!(RNsga2, PyRNsga2, "RNsga2");
+
+// Implement PyO3 methods
+#[pymethods]
+impl PyRNsga2 {
+    #[new]
+    #[pyo3(signature = (
+        reference_points,
+        sampler,
+        crossover,
+        mutation,
+        fitness_fn,
+        n_vars,
+        pop_size,
+        n_offsprings,
+        num_iterations,
+        epsilon = 0.001,
+        mutation_rate=0.1,
+        crossover_rate=0.9,
+        keep_infeasible=false,
+        verbose=true,
+        duplicates_cleaner=None,
+        constraints_fn=None,
+        lower_bound=None,
+        upper_bound=None
+    ))]
+    pub fn py_new<'py>(
+        reference_points: &Bound<'py, PyArray2<f64>>,
+        sampler: PyObject,
+        crossover: PyObject,
+        mutation: PyObject,
+        fitness_fn: PyObject,
+        n_vars: usize,
+        pop_size: usize,
+        n_offsprings: usize,
+        num_iterations: usize,
+        epsilon: f64,
+        mutation_rate: f64,
+        crossover_rate: f64,
+        keep_infeasible: bool,
+        verbose: bool,
+        duplicates_cleaner: Option<PyObject>,
+        constraints_fn: Option<PyObject>,
+        // Optional lower bound for each gene.
+        lower_bound: Option<f64>,
+        // Optional upper bound for each gene.
+        upper_bound: Option<f64>,
+    ) -> PyResult<Self> {
+        // Unwrap the genetic operators
+        let sampler_box = unwrap_sampling_operator(sampler)?;
+        let crossover_box = unwrap_crossover_operator(crossover)?;
+        let mutation_box = unwrap_mutation_operator(mutation)?;
+        let duplicates_box = if let Some(py_obj) = duplicates_cleaner {
+            Some(unwrap_duplicates_cleaner(py_obj)?)
+        } else {
+            None
+        };
+
+        // Build the MANDATORY population-level fitness closure
+        let fitness_closure = create_population_fitness_closure(fitness_fn)?;
+
+        // Build OPTIONAL population-level constraints closure
+        let constraints_closure = if let Some(py_obj) = constraints_fn {
+            Some(create_population_constraints_closure(py_obj)?)
+        } else {
+            None
+        };
+
+        // Convert PyArray2 to Array2
+        let reference_points_array = reference_points.to_owned_array();
+
+        // Create an instance of the selection/survival struct
+        let selector_box = Box::new(RankAndCrowdingSelection::new_with_comparison(
+            DiversityComparison::Minimize,
+        ));
+        let survivor_box = Box::new(RankReferencePointsSurvival::new(
+            reference_points_array,
+            epsilon,
+        ));
+
+        // Create the Rust struct
+        let rs_obj = RNsga2::new(
+            sampler_box,
+            crossover_box,
+            mutation_box,
+            selector_box,
+            survivor_box,
+            duplicates_box,
+            fitness_closure,
+            n_vars,
+            pop_size,
+            n_offsprings,
+            num_iterations,
+            mutation_rate,
+            crossover_rate,
+            keep_infeasible,
+            verbose,
+            constraints_closure,
+            lower_bound,
+            upper_bound,
+        )?;
+
+        Ok(Self { inner: rs_obj })
+    }
+}
diff --git a/src/non_dominated_sorting/distance.rs b/src/diversity_metrics/crowding.rs
similarity index 100%
rename from src/non_dominated_sorting/distance.rs
rename to src/diversity_metrics/crowding.rs
diff --git a/src/diversity_metrics/mod.rs b/src/diversity_metrics/mod.rs
new file mode 100644
index 0000000..9bdb9ce
--- /dev/null
+++ b/src/diversity_metrics/mod.rs
@@ -0,0 +1,5 @@
+pub mod crowding;
+pub mod reference;
+
+pub use crowding::crowding_distance;
+pub use reference::{reference_points_rank_distance, weighted_distance_matrix};
diff --git a/src/diversity_metrics/reference.rs b/src/diversity_metrics/reference.rs
new file mode 100644
index 0000000..88db56d
--- /dev/null
+++ b/src/diversity_metrics/reference.rs
@@ -0,0 +1,243 @@
+use crate::genetic::PopulationFitness;
+use numpy::ndarray::{Array1, Array2, Axis};
+use std::cmp::Ordering; // Assume PopulationFitness is defined as Array2<f64> elsewhere.
+
+// ---------------------------------------------------------------------------
+// Auxiliary Functions for Distance and Ranking Computations
+// ---------------------------------------------------------------------------
+
+/// Computes the ideal point from a fitness matrix.
+/// Each element of the returned array is the minimum value along the corresponding column.
+pub fn get_nideal(population_fitness: &PopulationFitness) -> Array1<f64> {
+    population_fitness.fold_axis(Axis(0), f64::INFINITY, |a, &b| a.min(b))
+}
+
+/// Computes the nadir point from a fitness matrix.
+/// Each element of the returned array is the maximum value along the corresponding column.
+pub fn get_nadir(population_fitness: &PopulationFitness) -> Array1<f64> {
+    population_fitness.fold_axis(Axis(0), f64::NEG_INFINITY, |a, &b| a.max(b))
+}
+
+/// Computes the weighted, normalized Euclidean distance between two objective vectors `f1` and `f2`.
+/// Normalization is performed using the provided ideal (`nideal`) and nadir (`nadir`) points.
+/// If for any objective the range (nadir - nideal) is zero, the normalized difference is set to 0.0.
+pub fn weighted_normalized_euclidean_distance(
+    f1: &Array1<f64>,
+    f2: &Array1<f64>,
+    weights: &Array1<f64>,
+    nideal: &Array1<f64>,
+    nadir: &Array1<f64>,
+) -> f64 {
+    let diff = f1 - f2;
+    let mut sum_sq = 0.0;
+    for j in 0..diff.len() {
+        let range = nadir[j] - nideal[j];
+        let normalized_diff = if range == 0.0 { 0.0 } else { diff[j] / range };
+        sum_sq += weights[j] * normalized_diff * normalized_diff;
+    }
+    sum_sq.sqrt()
+}
+
+/// Computes the ranking for a single column of distances.
+/// Given a slice of distances, it returns a vector of the same length where each element
+/// represents the order (i.e., rank) of the corresponding solution in the sorted order (ascending).
+///
+/// For example, if the distances in a column are [0.3, 0.1, 0.5],
+/// the sorted order is [0.1, 0.3, 0.5] and the corresponding ranks are:
+/// - The solution with 0.1 gets rank 0,
+/// - The solution with 0.3 gets rank 1,
+/// - The solution with 0.5 gets rank 2.
+fn compute_rank_for_column(distances: &[f64]) -> Vec<usize> {
+    let n = distances.len();
+    let mut indices: Vec<usize> = (0..n).collect();
+    indices.sort_by(|&i, &j| {
+        distances[i]
+            .partial_cmp(&distances[j])
+            .unwrap_or(Ordering::Equal)
+    });
+    let mut ranks = vec![0; n];
+    for (order, &i) in indices.iter().enumerate() {
+        ranks[i] = order;
+    }
+    ranks
+}
+
+/// Given a matrix of distances (with shape (n, P), where n is the number of solutions
+/// in the front and P is the number of reference points), computes the ranking for each solution
+/// as:
+///
+///     r_i = min_{k=1..P} r_{ik}
+///
+/// where r_{ik} is the rank (order) of the distance d_{ik} among all solutions for reference point k.
+fn ranking_by_distances(distances: &Array2<f64>) -> Array1<usize> {
+    let (n, p) = distances.dim();
+    let mut rank_matrix: Vec<Vec<usize>> = Vec::with_capacity(p);
+    for k in 0..p {
+        let col = distances.index_axis(Axis(1), k);
+        let col_vec = col.to_vec();
+        let ranks = compute_rank_for_column(&col_vec);
+        rank_matrix.push(ranks);
+    }
+    let mut final_ranks = vec![usize::MAX; n];
+    for i in 0..n {
+        for k in 0..p {
+            final_ranks[i] = final_ranks[i].min(rank_matrix[k][i]);
+        }
+    }
+    Array1::from(final_ranks)
+}
+
+/// Computes the weighted distance matrix for a front.
+///
+/// # Parameters
+/// - `fitness`: A PopulationFitness (Array2<f64>) representing the objectives of each solution in the front.
+/// - `reference`: A matrix of reference points (Array2<f64>). In the RNSGA2 survival process, this function is used
+///   with the reference points to compute the distance between each solution and each reference point.
+///
+/// # Process
+/// 1. Determines the number of solutions (n) and the number of objectives (m).
+/// 2. Computes the ideal and nadir points using `get_nideal` and `get_nadir`.
+/// 3. Sets the weights as 1/m for each objective.
+/// 4. Computes the weighted, normalized Euclidean distance between each solution in the fitness matrix
+///    and each reference point.
+///
+/// # Returns
+/// Returns a distance matrix (Array2<f64>) with dimensions (n, r), where r is the number of rows of the reference matrix.
+pub fn weighted_distance_matrix(
+    fitness: &PopulationFitness,
+    reference: &Array2<f64>,
+) -> Array2<f64> {
+    let (n, m) = fitness.dim();
+    let r = reference.nrows();
+    let ideal = get_nideal(fitness);
+    let nadir = get_nadir(fitness);
+    // Set weights as 1/m for each objective.
+    let weights = Array1::from_elem(m, 1.0 / (m as f64));
+    let mut dist_matrix = Array2::<f64>::zeros((n, r));
+
+    for i in 0..n {
+        let f_i = fitness.row(i).to_owned();
+        for j in 0..r {
+            let r_j = reference.row(j).to_owned();
+            let d = weighted_normalized_euclidean_distance(&f_i, &r_j, &weights, &ideal, &nadir);
+            dist_matrix[[i, j]] = d;
+        }
+    }
+    dist_matrix
+}
+
+/// Computes the ranking based on reference points.
+/// It uses `weighted_distance_matrix` to obtain the distance matrix between each solution
+/// in the front and the reference points, then calculates the final ranking with `ranking_by_distances`.
+///
+/// # Parameters
+/// - `fitness`: A PopulationFitness matrix (Array2<f64>) containing the solutions of the front.
+/// - `reference`: A matrix of reference points (Array2<f64>).
+///
+/// # Returns
+/// Returns an Array1<f64> where each element represents the ranking for each solution.
+/// This ranking is used as the diversity measure in the survival process of the RNSGA2 algorithm.
+pub fn reference_points_rank_distance(
+    fitness: &PopulationFitness,
+    reference: &Array2<f64>,
+) -> Array1<f64> {
+    let distances = weighted_distance_matrix(fitness, reference);
+    let ranks = ranking_by_distances(&distances);
+    ranks.mapv(|r| r as f64)
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ndarray::array;
+
+    #[test]
+    fn test_get_nideal() {
+        // Example fitness matrix (3 solutions, 2 objectives)
+        let fitness = array![[1.0, 4.0], [2.0, 3.0], [0.5, 5.0]];
+        let ideal = get_nideal(&fitness);
+        // For each objective, the expected minimum value:
+        // First objective: min(1.0, 2.0, 0.5) = 0.5
+        // Second objective: min(4.0, 3.0, 5.0) = 3.0
+        assert_eq!(ideal, array![0.5, 3.0]);
+    }
+
+    #[test]
+    fn test_get_nadir() {
+        // Example fitness matrix (3 solutions, 2 objectives)
+        let fitness = array![[1.0, 4.0], [2.0, 3.0], [0.5, 5.0]];
+        let nadir = get_nadir(&fitness);
+        // For each objective, the expected maximum value:
+        // First objective: max(1.0, 2.0, 0.5) = 2.0
+        // Second objective: max(4.0, 3.0, 5.0) = 5.0
+        assert_eq!(nadir, array![2.0, 5.0]);
+    }
+
+    #[test]
+    fn test_weighted_normalized_euclidean_distance() {
+        let f1 = array![1.0, 2.0, 3.0];
+        let f2 = array![2.0, 3.0, 4.0];
+        // Assume 3 objectives; therefore, weights are 1/3 for each.
+        let weights = Array1::from_elem(3, 1.0 / 3.0);
+        // Define arbitrary ideal and nadir points for normalization:
+        let ideal = array![0.0, 0.0, 0.0];
+        let nadir = array![10.0, 10.0, 10.0];
+        let distance = weighted_normalized_euclidean_distance(&f1, &f2, &weights, &ideal, &nadir);
+        // The normalized differences will be (1/10, 1/10, 1/10)
+        // Sum of squares: 3*(0.01) = 0.03, sqrt(0.03) ≈ 0.1732.
+        assert!((distance - 0.1732).abs() < 1e-4);
+    }
+
+    #[test]
+    fn test_compute_rank_for_column() {
+        let distances = vec![0.3, 0.1, 0.5, 0.2];
+        let ranks = compute_rank_for_column(&distances);
+        // Sorted order: 0.1 (index 1), 0.2 (index 3), 0.3 (index 0), 0.5 (index 2)
+        // Expected ranks: [2, 0, 3, 1]
+        assert_eq!(ranks, vec![2, 0, 3, 1]);
+    }
+
+    #[test]
+    fn test_ranking_by_distances() {
+        // Example distance matrix (3 solutions, 2 reference points)
+        let distances = array![[0.3, 0.4], [0.1, 0.2], [0.5, 0.6]];
+        // For column 0: sorted order is [0.1 (index 1), 0.3 (index 0), 0.5 (index 2)] → ranks: [1, 0, 2]
+        // For column 1: sorted order is [0.2 (index 1), 0.4 (index 0), 0.6 (index 2)] → ranks: [1, 0, 2]
+        // Final ranking for each solution is the minimum across columns:
+        // Solution 0: min(1,1)=1, Solution 1: min(0,0)=0, Solution 2: min(2,2)=2.
+        let final_ranks = ranking_by_distances(&distances);
+        assert_eq!(final_ranks, array![1, 0, 2]);
+    }
+
+    #[test]
+    fn test_weighted_distance_matrix() {
+        // Example fitness matrix (3 solutions, 2 objectives)
+        let fitness = array![[1.0, 4.0], [2.0, 3.0], [0.5, 5.0]];
+        // Reference matrix: for this test we define it with 2 reference points.
+        let reference = array![[1.0, 4.0], [2.0, 3.0]];
+        // Calculate the distance matrix: dimensions (n, r) = (3, 2)
+        let dist_matrix = weighted_distance_matrix(&fitness, &reference);
+        // For example, the distance between solution 0 and reference point 0:
+        // Since fitness[0] == reference[0], the distance should be 0.
+        assert!((dist_matrix[[0, 0]] - 0.0).abs() < 1e-6);
+        // Verify that the matrix has the expected dimensions.
+        assert_eq!(dist_matrix.dim(), (3, 2));
+    }
+
+    #[test]
+    fn test_reference_points_rank_distance() {
+        // Example fitness matrix (3 solutions, 2 objectives)
+        let fitness = array![[1.0, 4.0], [2.0, 3.0], [0.5, 5.0]];
+        // Reference matrix: for this test we use the same matrix as fitness,
+        // so each solution will have at least one identical reference point (distance 0).
+        let reference = fitness.clone();
+        let rank_vector = reference_points_rank_distance(&fitness, &reference);
+        // Since each solution is compared with itself (distance 0) in at least one reference point,
+        // the resulting ranking should be 0 for each solution.
+        let expected = array![0.0, 0.0, 0.0];
+        assert_eq!(rank_vector, expected);
+    }
+}
diff --git a/src/evaluator.rs b/src/evaluator.rs
index 3a0ee8c..bd8f79e 100644
--- a/src/evaluator.rs
+++ b/src/evaluator.rs
@@ -1,6 +1,6 @@
 use crate::{
     genetic::{Fronts, Population, PopulationConstraints, PopulationFitness, PopulationGenes},
-    non_dominated_sorting::{crowding_distance, fast_non_dominated_sorting},
+    non_dominated_sorting::fast_non_dominated_sorting,
 };
 use numpy::ndarray::{Array1, Axis};
 
@@ -19,7 +19,7 @@ pub struct Evaluator {
 
 impl Evaluator {
     /// Creates a new `Evaluator` with a fitness function, an optional constraints function,
-    /// a flag to keep infeasible individuals, and optional lower/upper bounds for the genes.
+    /// a flag to keep infeasible individuals, and optional lower/upper bounds.
     pub fn new(
         fitness_fn: Box<dyn Fn(&PopulationGenes) -> PopulationFitness>,
         constraints_fn: Option<Box<dyn Fn(&PopulationGenes) -> PopulationConstraints>>,
@@ -63,7 +63,7 @@ impl Evaluator {
             let n = genes.nrows();
             let mut feasible_indices: Vec<usize> = (0..n).collect();
 
-            // Filter individuals that do not satisfy the constraints function (if provided)
+            // Filter individuals that do not satisfy the constraints function (if provided).
             if let Some(ref c) = constraints {
                 feasible_indices = feasible_indices
                     .into_iter()
@@ -88,45 +88,34 @@ impl Evaluator {
                     .collect();
             }
 
-            // Filter all relevant arrays (genes, fitness, and constraints if present)
+            // Filter all relevant arrays (genes, fitness, and constraints if present).
             genes = genes.select(Axis(0), &feasible_indices);
             fitness = fitness.select(Axis(0), &feasible_indices);
             constraints = constraints.map(|c_array| c_array.select(Axis(0), &feasible_indices));
         }
 
-        // Return an empty Vec if no feasible individuals remain after filtering
+        // Return an empty Vec if no feasible individuals remain after filtering.
         if genes.nrows() == 0 {
-            return Vec::new(); // Possible empty result due to constraints
+            return Vec::new();
         }
 
         let sorted_fronts = fast_non_dominated_sorting(&fitness);
         let mut results: Fronts = Vec::new();
 
-        // For each front (rank = front_index), extract the sub-population.
+        // For each front (with rank = front_index), extract the sub-population.
         for (front_index, indices) in sorted_fronts.iter().enumerate() {
-            // Slice out the genes and fitness for just these individuals.
             let front_genes = genes.select(Axis(0), &indices[..]);
             let front_fitness = fitness.select(Axis(0), &indices[..]);
-
-            // If constraints exist, slice them out too.
             let front_constraints = constraints
                 .as_ref()
                 .map(|c| c.select(Axis(0), &indices[..]));
 
-            // Calculate crowding distance for just these individuals.
-            let cd_front = crowding_distance(&front_fitness);
-
             // Create a rank Array1 (one rank value per individual in the front).
             let rank_arr = Array1::from_elem(indices.len(), front_index);
 
             // Build a `Population` representing this entire front.
-            let population_front = Population {
-                genes: front_genes,
-                fitness: front_fitness,
-                constraints: front_constraints,
-                rank: rank_arr,
-                crowding_distance: cd_front,
-            };
+            let population_front =
+                Population::new(front_genes, front_fitness, front_constraints, rank_arr);
 
             results.push(population_front);
         }
@@ -167,6 +156,7 @@ mod tests {
 
     #[test]
     fn test_evaluator_evaluate_fitness() {
+        // Using the Crowding metric in this test.
         let evaluator = Evaluator::new(Box::new(fitness_fn), None, true, None, None);
 
         let population_genes = array![[1.0, 2.0], [3.0, 4.0], [0.0, 0.0]];
@@ -191,9 +181,9 @@ mod tests {
         );
 
         let population_genes = array![
-            [1.0, 2.0], // Feasible (sum=3, sum-10=-7; genes ≥ 0)
-            [3.0, 4.0], // Feasible (sum=7, sum-10=-3; genes ≥ 0)
-            [5.0, 6.0]  // Infeasible (sum=11, sum-10=1>0)
+            [1.0, 2.0], // Feasible (sum=3, 3-10=-7; genes ≥ 0)
+            [3.0, 4.0], // Feasible (sum=7, 7-10=-3; genes ≥ 0)
+            [5.0, 6.0]  // Infeasible (sum=11, 11-10=1>0)
         ];
 
         if let Some(constraints_array) = evaluator.evaluate_constraints(&population_genes) {
@@ -240,18 +230,13 @@ mod tests {
             "Total individuals across all fronts should be 3."
         );
 
-        // Check that each front has matching shapes for rank and crowding distance.
+        // Check that in each front the length of the rank array matches.
         for (front_index, front_pop) in fronts.iter().enumerate() {
             let n = front_pop.genes.nrows();
             assert_eq!(
                 front_pop.rank.len(),
                 n,
-                "Rank length should match number of individuals in front"
-            );
-            assert_eq!(
-                front_pop.crowding_distance.len(),
-                n,
-                "Crowding distance length should match individuals in front"
+                "Rank length should match number of individuals in the front"
             );
             for &r in front_pop.rank.iter() {
                 assert_eq!(
@@ -263,7 +248,7 @@ mod tests {
                 assert_eq!(
                     c.nrows(),
                     n,
-                    "Constraints rows should match number of individuals"
+                    "Number of rows in constraints should match number of individuals"
                 );
             }
         }
@@ -272,7 +257,7 @@ mod tests {
     #[test]
     fn test_evaluator_build_fronts_with_infeasible_and_bounds() {
         // Use constraints function and bounds. Also, keep_infeasible is false.
-        // The bounds here require that every gene must be between 0 and 5.
+        // The bounds require each gene to be between 0 and 5.
         let evaluator = Evaluator::new(
             Box::new(fitness_fn),
             Some(Box::new(constraints_fn)),
@@ -284,8 +269,8 @@ mod tests {
         // Create a population with 3 individuals:
         //   - [1.0, 2.0]: Feasible (sum=3, all genes between 0 and 5)
         //   - [3.0, 4.0]: Feasible (sum=7, all genes between 0 and 5)
-        //   - [6.0, 1.0]: Infeasible (fails upper bound, since 6.0 > 5.0)
-        let population_genes = array![[1.0, 2.0], [3.0, 4.0], [6.0, 1.0],];
+        //   - [6.0, 1.0]: Infeasible (6.0 > 5.0)
+        let population_genes = array![[1.0, 2.0], [3.0, 4.0], [6.0, 1.0]];
 
         // Build fronts.
         let fronts = evaluator.build_fronts(population_genes);
diff --git a/src/genetic.rs b/src/genetic.rs
index cdaca48..01c4a33 100644
--- a/src/genetic.rs
+++ b/src/genetic.rs
@@ -5,75 +5,55 @@ use numpy::ndarray::{concatenate, Array1, Array2, ArrayViewMut1, Axis};
 pub type Genes = Array1<f64>;
 pub type GenesMut<'a> = ArrayViewMut1<'a, f64>;
 
-/// The `Parents` type represents the input for a binary genetic operator, such as a crossover operator.
-/// It is a tuple of two 2-dimensional arrays (`Array2<f64>`) where `Parents.0[i]` will be operated with
-/// `Parents.1[i]` for each `i` in the population size. Each array corresponds to one "parent" group
-/// participating in the operation.
-pub type Parents = (Array1<f64>, Array1<f64>);
-
-/// The `Children` type defines the output of a binary genetic operator, such as the crossover operator.
-/// It is a tuple of two 2-dimensional arrays (`Array2<f64>`) where each array represents the resulting
-/// offspring derived from the corresponding parent arrays in `Parents`.
-pub type Children = (Array1<f64>, Array1<f64>);
-
-/// The `PopulationGenes` type defines the current set of individuals in the population.
-/// It is represented as a 2-dimensional array (`Array2<f64>`), where each row corresponds to an individual.
-pub type PopulationGenes = Array2<f64>;
-
-/// Fitness associated to one Genes
-pub type IndividualFitness = Array1<f64>;
-/// PopulationGenes Fitness
-pub type PopulationFitness = Array2<f64>;
-
-pub type IndividualConstraints = Array1<f64>;
-
-pub type PopulationConstraints = Array2<f64>;
-
+/// Represents an individual with genes, fitness, constraints (if any),
+/// rank, and an optional diversity metric.
 pub struct Individual {
     pub genes: Genes,
-    pub fitness: IndividualFitness,
-    pub constraints: Option<IndividualConstraints>,
+    pub fitness: Array1<f64>,
+    pub constraints: Option<Array1<f64>>,
     pub rank: usize,
-    pub crowding_distance: f64,
+    pub diversity_metric: Option<f64>,
 }
 
 impl Individual {
     pub fn new(
         genes: Genes,
-        fitness: IndividualFitness,
-        constraints: Option<IndividualConstraints>,
+        fitness: Array1<f64>,
+        constraints: Option<Array1<f64>>,
         rank: usize,
-        crowding_distance: f64,
+        diversity_metric: Option<f64>,
     ) -> Self {
         Self {
             genes,
             fitness,
             constraints,
             rank,
-            crowding_distance,
+            diversity_metric,
         }
     }
 
     pub fn is_feasible(&self) -> bool {
         match &self.constraints {
             None => true,
-            Some(c) => {
-                let sum: f64 = c.iter().sum();
-                sum <= 0.0
-            }
+            Some(c) => c.iter().sum::<f64>() <= 0.0,
         }
     }
 }
 
-/// The `Population` struct containing genes, fitness, rank, and crowding distance.
-/// `rank` and `crowding_distance` are optional and may be set during the process.
+/// Type aliases to work with populations.
+pub type PopulationGenes = Array2<f64>;
+pub type PopulationFitness = Array2<f64>;
+pub type PopulationConstraints = Array2<f64>;
+
+/// The `Population` struct contains genes, fitness, constraints (if any),
+/// rank, and optionally a diversity metric vector.
 #[derive(Debug)]
 pub struct Population {
     pub genes: PopulationGenes,
     pub fitness: PopulationFitness,
     pub constraints: Option<PopulationConstraints>,
     pub rank: Array1<usize>,
-    pub crowding_distance: Array1<f64>,
+    pub diversity_metric: Option<Array1<f64>>,
 }
 
 impl Clone for Population {
@@ -83,79 +63,103 @@ impl Clone for Population {
             fitness: self.fitness.clone(),
             constraints: self.constraints.clone(),
             rank: self.rank.clone(),
-            crowding_distance: self.crowding_distance.clone(),
+            diversity_metric: self.diversity_metric.clone(),
         }
     }
 }
 
 impl Population {
-    /// Creates a new `Population` instance with the given genes and fitness.
-    /// `rank` and `crowding_distance` are initially `None`.
+    /// Creates a new `Population` instance with the given genes, fitness, constraints, and rank.
+    /// The `diversity_metric` field is set to `None` by default.
     pub fn new(
         genes: PopulationGenes,
         fitness: PopulationFitness,
         constraints: Option<PopulationConstraints>,
         rank: Array1<usize>,
-        crowding_distance: Array1<f64>,
     ) -> Self {
         Self {
             genes,
             fitness,
             constraints,
             rank,
-            crowding_distance,
+            diversity_metric: None, // Initialized to None by default.
         }
     }
 
     /// Retrieves an `Individual` from the population by index.
     pub fn get(&self, idx: usize) -> Individual {
         let constraints = self.constraints.as_ref().map(|c| c.row(idx).to_owned());
-
+        let diversity = self.diversity_metric.as_ref().map(|dm| dm[idx]);
         Individual::new(
             self.genes.row(idx).to_owned(),
             self.fitness.row(idx).to_owned(),
             constraints,
             self.rank[idx],
-            self.crowding_distance[idx],
+            diversity,
         )
     }
 
     /// Returns a new `Population` containing only the individuals at the specified indices.
-    /// Indices may be repeated, resulting in repeated individuals in the new population.
     pub fn selected(&self, indices: &[usize]) -> Population {
         let genes = self.genes.select(Axis(0), indices);
         let fitness = self.fitness.select(Axis(0), indices);
         let rank = self.rank.select(Axis(0), indices);
-        let crowding_distance = self.crowding_distance.select(Axis(0), indices);
-
+        let diversity_metric = self
+            .diversity_metric
+            .as_ref()
+            .map(|dm| dm.select(Axis(0), indices));
         let constraints = self
             .constraints
             .as_ref()
             .map(|c| c.select(Axis(0), indices));
 
-        Population::new(genes, fitness, constraints, rank, crowding_distance)
+        Population::new(genes, fitness, constraints, rank).with_diversity(diversity_metric)
     }
 
-    /// Returns the number of individuals in this population.
+    /// Returns the number of individuals in the population.
     pub fn len(&self) -> usize {
         self.genes.nrows()
     }
+
     /// Returns a new `Population` containing only the individuals with rank = 0.
     pub fn best(&self) -> Population {
         let indices: Vec<usize> = self
             .rank
             .iter()
             .enumerate()
-            .filter_map(|(i, &rank)| if rank == 0 { Some(i) } else { None })
+            .filter_map(|(i, &r)| if r == 0 { Some(i) } else { None })
             .collect();
         self.selected(&indices)
     }
+
+    /// Auxiliary method to chain the assignment of `diversity_metric` in the `selected` method.
+    fn with_diversity(mut self, diversity_metric: Option<Array1<f64>>) -> Self {
+        self.diversity_metric = diversity_metric;
+        self
+    }
+
+    /// Updates the population's `diversity_metric` field.
+    ///
+    /// This method validates that the provided `diversity` vector has the same number of elements
+    /// as individuals in the population. If not, it returns an error.
+    pub fn set_diversity(&mut self, diversity: Array1<f64>) -> Result<(), String> {
+        if diversity.len() != self.len() {
+            return Err(format!(
+                "The diversity vector has length {} but the population contains {} individuals.",
+                diversity.len(),
+                self.len()
+            ));
+        }
+        self.diversity_metric = Some(diversity);
+        Ok(())
+    }
 }
 
+/// Type alias for a vector of `Population` representing multiple fronts.
 pub type Fronts = Vec<Population>;
 
-/// An extension trait for the `Fronts` type to add a `.flatten()` method
-/// that combines multiple fronts into a single `Population`.
+/// An extension trait for `Fronts` that adds a `.to_population()` method
+/// which flattens multiple fronts into a single `Population`.
 pub trait FrontsExt {
     fn to_population(&self) -> Population;
 }
@@ -167,19 +171,25 @@ impl FrontsExt for Vec<Population> {
         }
 
         let has_constraints = self[0].constraints.is_some();
+        let has_diversity = self[0].diversity_metric.is_some();
 
         let mut genes_views = Vec::new();
         let mut fitness_views = Vec::new();
         let mut rank_views = Vec::new();
-        let mut cd_views = Vec::new();
+        let mut diversity_views = Vec::new();
         let mut constraints_views = Vec::new();
 
         for front in self.iter() {
             genes_views.push(front.genes.view());
             fitness_views.push(front.fitness.view());
             rank_views.push(front.rank.view());
-            cd_views.push(front.crowding_distance.view());
-
+            if has_diversity {
+                let dm = front
+                    .diversity_metric
+                    .as_ref()
+                    .expect("Inconsistent diversity_metric among fronts");
+                diversity_views.push(dm.view());
+            }
             if has_constraints {
                 let c = front
                     .constraints
@@ -193,11 +203,17 @@ impl FrontsExt for Vec<Population> {
             concatenate(Axis(0), &genes_views[..]).expect("Error concatenating genes");
         let merged_fitness =
             concatenate(Axis(0), &fitness_views[..]).expect("Error concatenating fitness");
-
-        // **Concatenate** (Axis(0)) for 1D arrays rank & cd:
         let merged_rank =
-            concatenate(Axis(0), &rank_views[..]).expect("Error concatenating rank arrays"); // 1D result
-        let merged_cd = concatenate(Axis(0), &cd_views[..]).expect("Error concatenating cd arrays"); // 1D result
+            concatenate(Axis(0), &rank_views[..]).expect("Error concatenating rank arrays");
+
+        let merged_diversity = if has_diversity {
+            Some(
+                concatenate(Axis(0), &diversity_views[..])
+                    .expect("Error concatenating diversity arrays"),
+            )
+        } else {
+            None
+        };
 
         let merged_constraints = if has_constraints {
             Some(
@@ -213,7 +229,7 @@ impl FrontsExt for Vec<Population> {
             fitness: merged_fitness,
             constraints: merged_constraints,
             rank: merged_rank,
-            crowding_distance: merged_cd,
+            diversity_metric: merged_diversity,
         }
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index ee48fc6..91da108 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -7,6 +7,7 @@ pub mod genetic;
 pub mod prelude;
 
 mod algorithms;
+pub mod diversity_metrics;
 mod helpers;
 pub mod non_dominated_sorting;
 mod operators;
@@ -15,6 +16,7 @@ use pyo3::prelude::*;
 
 pub use algorithms::nsga2::PyNsga2;
 pub use algorithms::nsga3::PyNsga3;
+pub use algorithms::rnsga2::PyRNsga2;
 pub use algorithms::py_errors::NoFeasibleIndividualsError;
 pub use helpers::duplicates::{PyCloseDuplicatesCleaner, PyExactDuplicatesCleaner};
 pub use operators::py_operators::{
@@ -30,6 +32,7 @@ fn _pymoors(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     // Add classes from algorithms
     m.add_class::<PyNsga2>()?;
     m.add_class::<PyNsga3>()?;
+    m.add_class::<PyRNsga2>()?;
 
     // Add classes from operators
     m.add_class::<PyBitFlipMutation>()?;
diff --git a/src/non_dominated_sorting/mod.rs b/src/non_dominated_sorting/mod.rs
index c630521..77a4b77 100644
--- a/src/non_dominated_sorting/mod.rs
+++ b/src/non_dominated_sorting/mod.rs
@@ -1,5 +1,3 @@
-mod distance;
 mod dominator;
 
-pub use distance::crowding_distance;
 pub use dominator::fast_non_dominated_sorting;
diff --git a/src/operators/mod.rs b/src/operators/mod.rs
index 0d23d16..97373e2 100644
--- a/src/operators/mod.rs
+++ b/src/operators/mod.rs
@@ -164,6 +164,7 @@ pub trait CrossoverOperator: GeneticOperator {
 }
 
 // Enum to represent the result of a tournament duel.
+#[derive(Debug, PartialEq, Eq)]
 pub enum DuelResult {
     LeftWins,
     RightWins,
@@ -260,5 +261,5 @@ pub trait SelectionOperator: GeneticOperator {
 
 pub trait SurvivalOperator: GeneticOperator {
     /// Selects the individuals that will survive to the next generation.
-    fn operate(&self, fronts: &Fronts, n_survive: usize) -> Population;
+    fn operate(&self, fronts: &mut Fronts, n_survive: usize) -> Population;
 }
diff --git a/src/operators/selection/mod.rs b/src/operators/selection/mod.rs
index df6070a..bf43a74 100644
--- a/src/operators/selection/mod.rs
+++ b/src/operators/selection/mod.rs
@@ -2,4 +2,4 @@ pub mod random_tournament;
 pub mod rank_and_crowding_tournament;
 
 pub use random_tournament::RandomSelection;
-pub use rank_and_crowding_tournament::RankAndCrowdingSelection;
+pub use rank_and_crowding_tournament::{DiversityComparison, RankAndCrowdingSelection};
diff --git a/src/operators/selection/rank_and_crowding_tournament.rs b/src/operators/selection/rank_and_crowding_tournament.rs
index c851995..6c1e24a 100644
--- a/src/operators/selection/rank_and_crowding_tournament.rs
+++ b/src/operators/selection/rank_and_crowding_tournament.rs
@@ -1,17 +1,34 @@
-use std::fmt::Debug;
-
 use crate::genetic::Individual;
 use crate::operators::{DuelResult, GeneticOperator, SelectionOperator};
 use rand::RngCore;
+use std::fmt::Debug;
 
-// TODO: Enable pressure. Currently is fixed in 2
+/// Controls how the diversity (crowding) metric is compared during tournament selection.
+#[derive(Clone, Debug)]
+pub enum DiversityComparison {
+    /// Larger diversity (crowding) metric is preferred.
+    Maximize,
+    /// Smaller diversity (crowding) metric is preferred.
+    Minimize,
+}
 
 #[derive(Clone, Debug)]
-pub struct RankAndCrowdingSelection {}
+pub struct RankAndCrowdingSelection {
+    diversity_comparison: DiversityComparison,
+}
 
 impl RankAndCrowdingSelection {
+    /// Creates a new RankAndCrowdingSelection with the default diversity comparison (Maximize).
     pub fn new() -> Self {
-        Self {}
+        Self {
+            diversity_comparison: DiversityComparison::Maximize,
+        }
+    }
+    /// Creates a new RankAndCrowdingSelection with the specified diversity comparison.
+    pub fn new_with_comparison(comparison: DiversityComparison) -> Self {
+        Self {
+            diversity_comparison: comparison,
+        }
     }
 }
 
@@ -22,42 +39,55 @@ impl GeneticOperator for RankAndCrowdingSelection {
 }
 
 impl SelectionOperator for RankAndCrowdingSelection {
-    /// Runs the tournament selection on the given population and returns a vector of winners.
-    /// This example assumes binary tournaments (pressure = 2)
+    /// Runs tournament selection on the given population and returns the duel result.
+    /// This example assumes binary tournaments (pressure = 2).
     fn tournament_duel(
         &self,
         p1: &Individual,
         p2: &Individual,
-        rng: &mut dyn RngCore,
+        _rng: &mut dyn RngCore,
     ) -> DuelResult {
-        // get feasibility
+        // Check feasibility.
         let p1_feasible = p1.is_feasible();
         let p2_feasible = p2.is_feasible();
-        // get rank
+        // Retrieve rank.
         let p1_rank = p1.rank;
         let p2_rank = p2.rank;
-        // get cd
-        let p1_cd = p1.crowding_distance;
-        let p2_cd = p2.crowding_distance;
+        // Retrieve diversity (crowding) metric.
+        let p1_cd = p1.diversity_metric;
+        let p2_cd = p2.diversity_metric;
 
         let winner = if p1_feasible && !p2_feasible {
             DuelResult::LeftWins
         } else if p2_feasible && !p1_feasible {
             DuelResult::RightWins
         } else {
-            // Both feasible or both infeasible
+            // Both are either feasible or infeasible.
             if p1_rank < p2_rank {
                 DuelResult::LeftWins
             } else if p2_rank < p1_rank {
                 DuelResult::RightWins
             } else {
-                // Same rank, compare crowding distance
-                if p1_cd > p2_cd {
-                    DuelResult::LeftWins
-                } else if p1_cd < p2_cd {
-                    DuelResult::RightWins
-                } else {
-                    DuelResult::Tie
+                // When ranks are equal, compare the diversity metric.
+                match self.diversity_comparison {
+                    DiversityComparison::Maximize => {
+                        if p1_cd > p2_cd {
+                            DuelResult::LeftWins
+                        } else if p1_cd < p2_cd {
+                            DuelResult::RightWins
+                        } else {
+                            DuelResult::Tie
+                        }
+                    }
+                    DiversityComparison::Minimize => {
+                        if p1_cd < p2_cd {
+                            DuelResult::LeftWins
+                        } else if p1_cd > p2_cd {
+                            DuelResult::RightWins
+                        } else {
+                            DuelResult::Tie
+                        }
+                    }
                 }
             }
         };
@@ -69,28 +99,62 @@ impl SelectionOperator for RankAndCrowdingSelection {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::genetic::Population;
-    use numpy::ndarray::{arr1, arr2, Array2};
+    use crate::genetic::{Individual, Population};
+    use crate::operators::{DuelResult, SelectionOperator};
+    use ndarray::{arr1, arr2, Array2};
     use rand::prelude::*;
 
+    #[test]
+    fn test_default_diversity_comparison_maximize() {
+        let selector = RankAndCrowdingSelection::new();
+        match selector.diversity_comparison {
+            DiversityComparison::Maximize => assert!(true),
+            DiversityComparison::Minimize => panic!("Default should be Maximize"),
+        }
+    }
+
+    #[test]
+    fn test_tournament_duel_maximize() {
+        // Create two individuals using the actual Individual type.
+        // Both individuals have the same rank (0) but different diversity metrics.
+        // In Maximize mode, the individual with the higher diversity (10.0) should win.
+        let p1 = Individual::new(arr1(&[1.0, 2.0]), arr1(&[0.5]), None, 0, Some(10.0));
+        let p2 = Individual::new(arr1(&[3.0, 4.0]), arr1(&[0.6]), None, 0, Some(5.0));
+        let selector = RankAndCrowdingSelection::new(); // Default: Maximize
+        let mut rng = thread_rng();
+        let result = selector.tournament_duel(&p1, &p2, &mut rng);
+        assert_eq!(result, DuelResult::LeftWins);
+    }
+
+    #[test]
+    fn test_tournament_duel_minimize() {
+        // Create two individuals using the actual Individual type.
+        // Both individuals have the same rank (0) but different diversity metrics.
+        // In Minimize mode, the individual with the lower diversity (5.0) should win.
+        let p1 = Individual::new(arr1(&[1.0, 2.0]), arr1(&[0.5]), None, 0, Some(10.0));
+        let p2 = Individual::new(arr1(&[3.0, 4.0]), arr1(&[0.6]), None, 0, Some(5.0));
+        let selector = RankAndCrowdingSelection::new_with_comparison(DiversityComparison::Minimize);
+        let mut rng = thread_rng();
+        let result = selector.tournament_duel(&p1, &p2, &mut rng);
+        assert_eq!(result, DuelResult::RightWins);
+    }
+
     #[test]
     fn test_tournament_selection_no_constraints_basic() {
         // For a population of 4:
         // Rank: [0, 1, 0, 1]
-        // CD: [10.0, 5.0, 9.0, 1.0]
+        // Diversity (CD): [10.0, 5.0, 9.0, 1.0]
         let genes = arr2(&[[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]]);
         let fitness = arr2(&[[0.5], [0.6], [0.7], [0.8]]);
         let constraints = None;
         let rank = arr1(&[0, 1, 0, 1]);
-        let cd = arr1(&[10.0, 5.0, 9.0, 1.0]);
 
-        let population = Population::new(genes, fitness, constraints, rank, cd);
+        let population = Population::new(genes, fitness, constraints, rank);
 
-        // n_crossovers = 2
-        // total_needed = 2 * 2 * 2 = 8 participants → 4 tournaments → 4 winners.
+        // n_crossovers = 2 → total_needed = 8 participants → 4 tournaments → 4 winners.
         // After splitting: pop_a = 2 winners, pop_b = 2 winners.
         let n_crossovers = 2;
-        let selector = RankAndCrowdingSelection {};
+        let selector = RankAndCrowdingSelection::new();
         let mut rng = thread_rng();
         let (pop_a, pop_b) = selector.operate(&population, n_crossovers, &mut rng);
 
@@ -107,15 +171,12 @@ mod tests {
         let fitness = arr2(&[[0.5], [0.6]]);
         let constraints = Some(arr2(&[[-1.0, 0.0], [1.0, 1.0]]));
         let rank = arr1(&[0, 0]);
-        let cd = arr1(&[5.0, 10.0]);
+        let population = Population::new(genes, fitness, constraints, rank);
 
-        let population = Population::new(genes, fitness, constraints, rank, cd);
-
-        // n_crossovers = 1
-        // total_needed = 1 * 2 * 2 = 4 participants → 2 tournaments → 2 winners total.
+        // n_crossovers = 1 → total_needed = 4 participants → 2 tournaments → 2 winners total.
         // After splitting: pop_a = 1 winner, pop_b = 1 winner.
         let n_crossovers = 1;
-        let selector = RankAndCrowdingSelection {};
+        let selector = RankAndCrowdingSelection::new();
         let mut rng = thread_rng();
         let (pop_a, pop_b) = selector.operate(&population, n_crossovers, &mut rng);
 
@@ -126,25 +187,24 @@ mod tests {
 
     #[test]
     fn test_tournament_selection_same_rank_and_cd() {
-        // If two individuals have the same rank and the same CD,
-        // the second one wins in the event of a tie.
+        // If two individuals have the same rank and the same crowding distance,
+        // the tournament duel should result in a tie.
         let genes = arr2(&[[1.0, 2.0], [3.0, 4.0]]);
         let fitness = arr2(&[[0.5], [0.6]]);
         let constraints = None;
         let rank = arr1(&[0, 0]);
-        let cd = arr1(&[10.0, 10.0]);
 
-        let population = Population::new(genes, fitness, constraints, rank, cd);
+        let population = Population::new(genes, fitness, constraints, rank);
 
-        // n_crossovers = 1
-        // total_needed = 4 participants → 2 tournaments → 2 winners
-        // After splitting: pop_a = 1 winner, pop_b = 1 winner
+        // n_crossovers = 1 → total_needed = 4 participants → 2 tournaments → 2 winners.
+        // After splitting: pop_a = 1 winner, pop_b = 1 winner.
         let n_crossovers = 1;
-        let selector = RankAndCrowdingSelection {};
+        let selector = RankAndCrowdingSelection::new();
         let mut rng = thread_rng();
         let (pop_a, pop_b) = selector.operate(&population, n_crossovers, &mut rng);
 
-        // In a tie, the second individual chosen in the tournament wins.
+        // In a tie, the overall selection process must eventually choose winners.
+        // For this test, we ensure that each subpopulation has 1 individual.
         assert_eq!(pop_a.len(), 1);
         assert_eq!(pop_b.len(), 1);
     }
@@ -154,7 +214,6 @@ mod tests {
         // Large population test to ensure stability.
         let pop_size = 100;
         let n_genes = 5;
-
         let genes = Array2::from_shape_fn((pop_size, n_genes), |(i, _)| i as f64);
         let fitness = Array2::from_shape_fn((pop_size, 1), |(i, _)| i as f64 / 100.0);
         let constraints = None;
@@ -163,17 +222,12 @@ mod tests {
         let rank_vec: Vec<usize> = (0..pop_size).map(|_| rng.gen_range(0..5)).collect();
         let rank = arr1(&rank_vec);
 
-        let cd_vec: Vec<f64> = (0..pop_size).map(|_| rng.gen_range(0.0..10.0)).collect();
-        let cd = arr1(&cd_vec);
+        let population = Population::new(genes, fitness, constraints, rank);
 
-        let population = Population::new(genes, fitness, constraints, rank, cd);
-
-        // n_crossovers = 50
-        // total_needed = 50 * 2 * 2 = 200 participants → 100 tournaments → 100 winners.
+        // n_crossovers = 50 → total_needed = 200 participants → 100 tournaments → 100 winners.
         // After splitting: pop_a = 50 winners, pop_b = 50 winners.
         let n_crossovers = 50;
-        let selector = RankAndCrowdingSelection {};
-        let mut rng = thread_rng();
+        let selector = RankAndCrowdingSelection::new();
         let (pop_a, pop_b) = selector.operate(&population, n_crossovers, &mut rng);
 
         assert_eq!(pop_a.len(), 50);
@@ -183,23 +237,21 @@ mod tests {
     #[test]
     fn test_tournament_selection_single_tournament() {
         // One crossover:
-        // total_needed = 1 * 2 * 2 = 4 participants → 2 tournaments → 2 winners
-        // After splitting: pop_a = 1, pop_b = 1
+        // total_needed = 4 participants → 2 tournaments → 2 winners.
+        // After splitting: pop_a = 1, pop_b = 1.
         let genes = arr2(&[[10.0, 20.0], [30.0, 40.0]]);
         let fitness = arr2(&[[1.0], [2.0]]);
         let constraints = None;
         let rank = arr1(&[0, 1]);
-        let cd = arr1(&[5.0, 1.0]);
 
-        let population = Population::new(genes, fitness, constraints, rank, cd);
+        let population = Population::new(genes, fitness, constraints, rank);
 
         let n_crossovers = 1;
-        let selector = RankAndCrowdingSelection {};
+        let selector = RankAndCrowdingSelection::new();
         let mut rng = thread_rng();
         let (pop_a, pop_b) = selector.operate(&population, n_crossovers, &mut rng);
 
-        // The first individual has a better rank, so it should be one of the winners.
-        // The second winner would be from the second tournament.
+        // The individual with the better rank should win one tournament.
         assert_eq!(pop_a.len(), 1);
         assert_eq!(pop_b.len(), 1);
     }
diff --git a/src/operators/survival/mod.rs b/src/operators/survival/mod.rs
index ebc3036..64fcd2f 100644
--- a/src/operators/survival/mod.rs
+++ b/src/operators/survival/mod.rs
@@ -1,5 +1,7 @@
 pub mod rank_and_crowding;
+pub mod rank_reference_points;
 pub mod reference_points;
 
 pub use rank_and_crowding::RankCrowdingSurvival;
+pub use rank_reference_points::RankReferencePointsSurvival;
 pub use reference_points::ReferencePointsSurvival;
diff --git a/src/operators/survival/rank_and_crowding.rs b/src/operators/survival/rank_and_crowding.rs
index 0ac174a..5784445 100644
--- a/src/operators/survival/rank_and_crowding.rs
+++ b/src/operators/survival/rank_and_crowding.rs
@@ -1,6 +1,8 @@
+use std::fmt::Debug;
+
+use crate::diversity_metrics::crowding_distance;
 use crate::genetic::{Fronts, FrontsExt, Population};
 use crate::operators::{GeneticOperator, SurvivalOperator};
-use std::fmt::Debug;
 
 #[derive(Clone, Debug)]
 pub struct RankCrowdingSurvival;
@@ -18,14 +20,19 @@ impl RankCrowdingSurvival {
 }
 
 impl SurvivalOperator for RankCrowdingSurvival {
-    fn operate(&self, fronts: &Fronts, n_survive: usize) -> Population {
+    fn operate(&self, fronts: &mut Fronts, n_survive: usize) -> Population {
         // We will collect sub-populations (fronts) that survive here
         let mut chosen_fronts: Vec<Population> = Vec::new();
         let mut n_survivors = 0;
 
-        for front in fronts.iter() {
+        for front in fronts.iter_mut() {
             let front_size = front.len();
-
+            // Compute crowding distance
+            let cd = crowding_distance(&front.fitness);
+            // Set the crowding distance to the front population
+            front
+                .set_diversity(cd)
+                .expect("Failed to set diversity metric");
             // If this entire front fits into the survivor count
             if n_survivors + front_size <= n_survive {
                 chosen_fronts.push(front.clone());
@@ -35,7 +42,7 @@ impl SurvivalOperator for RankCrowdingSurvival {
                 let remaining = n_survive - n_survivors;
                 if remaining > 0 {
                     // Sort by crowding distance (descending)
-                    let cd = front.crowding_distance.clone();
+                    let cd = front.diversity_metric.clone().unwrap();
                     let mut indices: Vec<usize> = (0..front_size).collect();
                     indices.sort_by(|&i, &j| {
                         cd[j]
@@ -70,14 +77,12 @@ mod tests {
         let fitness = arr2(&[[0.1], [0.2], [0.3]]);
         let constraints: Option<Array2<f64>> = None;
         let rank = arr1(&[0, 0, 0]);
-        let crowding_distance = arr1(&[10.0, 5.0, 7.0]);
 
         let population = Population::new(
             genes.clone(),
             fitness.clone(),
             constraints.clone(),
             rank.clone(),
-            crowding_distance.clone(),
         );
         let mut fronts: Fronts = vec![population];
 
@@ -98,14 +103,12 @@ mod tests {
         let fitness = arr2(&[[0.1], [0.2], [0.3]]);
         let constraints: Option<Array2<f64>> = None;
         let rank = arr1(&[0, 0, 0]);
-        let crowding_distance = arr1(&[10.0, 5.0, 7.0]);
 
         let population = Population::new(
             genes.clone(),
             fitness.clone(),
             constraints.clone(),
             rank.clone(),
-            crowding_distance.clone(),
         );
         let mut fronts: Fronts = vec![population];
 
@@ -131,20 +134,17 @@ mod tests {
         let front1_fitness = arr2(&[[0.1], [0.2]]);
         let front1_constraints: Option<Array2<f64>> = None;
         let front1_rank = arr1(&[0, 0]);
-        let front1_cd = arr1(&[8.0, 9.0]);
 
         let front2_genes = arr2(&[[4.0, 5.0], [6.0, 7.0], [8.0, 9.0]]);
         let front2_fitness = arr2(&[[0.3], [0.4], [0.5]]);
         let front2_constraints: Option<Array2<f64>> = None;
         let front2_rank = arr1(&[1, 1, 1]);
-        let front2_cd = arr1(&[3.0, 10.0, 1.0]);
 
         let population1 = Population::new(
             front1_genes,
             front1_fitness,
             front1_constraints,
             front1_rank,
-            front1_cd,
         );
 
         let population2 = Population::new(
@@ -152,7 +152,6 @@ mod tests {
             front2_fitness,
             front2_constraints,
             front2_rank,
-            front2_cd,
         );
 
         let mut fronts: Fronts = vec![population1, population2];
diff --git a/src/operators/survival/rank_reference_points.rs b/src/operators/survival/rank_reference_points.rs
new file mode 100644
index 0000000..8c1e9f2
--- /dev/null
+++ b/src/operators/survival/rank_reference_points.rs
@@ -0,0 +1,141 @@
+use std::cmp::Ordering;
+use std::fmt::Debug;
+
+use ndarray::{Array1, Array2};
+
+use crate::diversity_metrics::{reference_points_rank_distance, weighted_distance_matrix};
+use crate::genetic::{Fronts, FrontsExt, Population};
+use crate::operators::{GeneticOperator, SurvivalOperator};
+
+#[derive(Clone, Debug)]
+pub struct RankReferencePointsSurvival {
+    reference_points: Array2<f64>,
+    epsilon: f64,
+}
+
+impl GeneticOperator for RankReferencePointsSurvival {
+    fn name(&self) -> String {
+        "RankReferencePointsSurvival".to_string()
+    }
+}
+
+impl RankReferencePointsSurvival {
+    pub fn new(reference_points: Array2<f64>, epsilon: f64) -> Self {
+        Self {
+            reference_points,
+            epsilon,
+        }
+    }
+}
+
+/// This function computes the weighted, normalized Euclidean distance matrix between each solution
+/// in the front (fitness matrix) and a set of reference points.
+/// (It is already defined in your diversity_metrics module.)
+/// Here we assume it is available as:
+///    weighted_distance_matrix(fitness: &PopulationFitness, reference: &Array2<f64>) -> Array2<f64>
+/// and that reference_points_rank_distance calls it to compute a ranking.
+///
+/// The splitting front procedure below uses weighted_distance_matrix(fitness, fitness)
+/// to compute the internal distances among solutions.
+
+impl SurvivalOperator for RankReferencePointsSurvival {
+    fn operate(&self, fronts: &mut Fronts, n_survive: usize) -> Population {
+        // We will collect sub-populations (fronts) that survive here.
+        let mut chosen_fronts: Vec<Population> = Vec::new();
+        let mut n_survivors = 0;
+
+        for front in fronts.iter_mut() {
+            let front_size = front.len();
+            if n_survivors + front_size <= n_survive {
+                // Compute the reference ranking for the entire front using the given reference points.
+                let ranking_by_distance =
+                    reference_points_rank_distance(&front.fitness, &self.reference_points);
+                // Set the diversity metric (crowding distance) for the whole front.
+                front
+                    .set_diversity(ranking_by_distance)
+                    .expect("Failed to set diversity metric");
+                // If the entire front fits, add it as is.
+                chosen_fronts.push(front.clone());
+                n_survivors += front_size;
+            } else {
+                // Only part of this front will survive (splitting front).
+                let remaining = n_survive - n_survivors;
+                if remaining > 0 {
+                    let front_size = front.len();
+                    // Compute the internal distance matrix among solutions in the front.
+                    // Here, we use the fitness matrix as both arguments so that we get pairwise distances.
+                    let mut internal_dist =
+                        weighted_distance_matrix(&front.fitness, &front.fitness);
+                    // Set the diagonal to infinity so that a solution is not considered close to itself.
+                    for i in 0..front_size {
+                        internal_dist[[i, i]] = f64::INFINITY;
+                    }
+                    // Get the reference ranking for the front (using the reference points).
+                    let reference_ranking =
+                        reference_points_rank_distance(&front.fitness, &self.reference_points);
+                    let reference_ranking_vec = reference_ranking.to_vec();
+
+                    // Create set S: all indices of solutions in the front,
+                    // sorted in ascending order by the reference ranking.
+                    let mut s: Vec<usize> = (0..front_size).collect();
+                    s.sort_by(|&i, &j| {
+                        reference_ranking_vec[i]
+                            .partial_cmp(&reference_ranking_vec[j])
+                            .unwrap_or(Ordering::Equal)
+                    });
+                    // Initialize a crowding vector with the initial reference ranking values.
+                    let mut crowding = reference_ranking_vec.clone();
+                    // Define a penalty Δ as the ceiling of half the front size.
+                    let delta = (front_size as f64 / 2.0).ceil();
+
+                    // Iterative procedure:
+                    // While S is not empty, select the best (lowest ranked) solution and penalize close neighbors.
+                    while !s.is_empty() {
+                        let i_star = s[0];
+                        s.remove(0); // Remove i_star from S.
+                                     // Identify group G: all remaining solutions j in S such that the internal distance
+                                     // between i_star and j is less than epsilon.
+                        let group: Vec<usize> = s
+                            .iter()
+                            .cloned()
+                            .filter(|&j| internal_dist[[i_star, j]] < self.epsilon)
+                            .collect();
+                        // Penalize each solution j in group by adding Δ to its reference ranking.
+                        for &j in group.iter() {
+                            crowding[j] = reference_ranking_vec[j] + delta;
+                        }
+                        // Remove all indices in G from S.
+                        s.retain(|&x| !group.contains(&x));
+                    }
+                    // Now, sort the indices of the front by the updated (modified) crowding values (ascending).
+                    let mut sorted_indices: Vec<usize> = (0..front_size).collect();
+                    sorted_indices.sort_by(|&i, &j| {
+                        crowding[i]
+                            .partial_cmp(&crowding[j])
+                            .unwrap_or(Ordering::Equal)
+                    });
+                    // Select the first 'remaining' indices.
+                    let selected_indices: Vec<usize> =
+                        sorted_indices.into_iter().take(remaining).collect();
+                    // Create a new sub-population from the splitting front using only the selected individuals.
+                    let mut selected_front = front.selected(&selected_indices);
+                    // Set the updated crowding (diversity) for the selected sub-front.
+                    let selected_crowding = Array1::from(
+                        selected_indices
+                            .iter()
+                            .map(|&i| crowding[i])
+                            .collect::<Vec<f64>>(),
+                    );
+                    selected_front
+                        .set_diversity(selected_crowding)
+                        .expect("Failed to set updated diversity metric");
+                    chosen_fronts.push(selected_front);
+                }
+                // No more slots remain; break out.
+                break;
+            }
+        }
+        // Finally, combine the chosen fronts into a single population.
+        chosen_fronts.to_population()
+    }
+}
diff --git a/src/operators/survival/reference_points.rs b/src/operators/survival/reference_points.rs
index 0b5f4ac..4ce2834 100644
--- a/src/operators/survival/reference_points.rs
+++ b/src/operators/survival/reference_points.rs
@@ -25,7 +25,7 @@ impl ReferencePointsSurvival {
 }
 
 impl SurvivalOperator for ReferencePointsSurvival {
-    fn operate(&self, fronts: &Fronts, n_survive: usize) -> Population {
+    fn operate(&self, fronts: &mut Fronts, n_survive: usize) -> Population {
         // Initialize a vector to store selected fronts (populations)
         let mut chosen_fronts: Vec<Population> = Vec::new();
         let mut n_survivors = 0;
@@ -282,11 +282,10 @@ mod tests {
             fitness.clone(),
             None,
             Array1::from(vec![1, 1, 1, 1, 1, 1]),
-            Array1::from(vec![0.5, 0.6, 0.7, 0.8, 0.9, 1.0]),
         );
 
         // Define fronts (for simplicity, all individuals are in the first front)
-        let fronts = vec![population.clone()];
+        let mut fronts = vec![population.clone()];
 
         // Generate reference points.
         let n_ref_points = 3;
@@ -298,7 +297,7 @@ mod tests {
 
         // Perform survival operation to select 4 individuals.
         let n_survive = 4;
-        let new_population = nsga3_survival.operate(&fronts, n_survive);
+        let new_population = nsga3_survival.operate(&mut fronts, n_survive);
 
         // Assert that the new population has 4 individuals.
         assert_eq!(new_population.len(), n_survive);
@@ -319,7 +318,6 @@ mod tests {
             fitness.clone(),
             None,
             Array1::from(vec![1, 1, 1]),
-            Array1::from(vec![0.5, 0.6, 0.7]),
         );
 
         let normalized = normalize_front(&population);
@@ -354,7 +352,6 @@ mod tests {
             fitness.clone(),
             None,
             Array1::from(vec![1, 1, 1]),
-            Array1::from(vec![0.5, 0.6, 0.7]),
         );
 
         let reference_points = array![[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
@@ -389,7 +386,6 @@ mod tests {
             fitness.clone(),
             None,
             Array1::from(vec![1, 1, 1, 1]),
-            Array1::from(vec![0.5, 0.6, 0.7, 0.8]),
         );
 
         // Define assignments of individuals to reference points.
diff --git a/tests/algorithms/test_small_real_biobj.py b/tests/algorithms/test_small_real_biobj.py
index f6763bb..ea3d9c4 100644
--- a/tests/algorithms/test_small_real_biobj.py
+++ b/tests/algorithms/test_small_real_biobj.py
@@ -10,6 +10,9 @@
 )
 from pymoors.typing import TwoDArray
 
+# FIXME: Once RNsga2 is exposed to the users, change the import
+from pymoors._pymoors import RNsga2
+
 
 def f1(x: float, y: float) -> float:
     """Objective 1: x^2 + y^2."""
@@ -79,7 +82,7 @@ def test_small_real_biobjective_nsag2():
         num_iterations=50,
         mutation_rate=0.1,
         crossover_rate=0.9,
-        duplicates_cleaner=CloseDuplicatesCleaner(epsilon=1e-16),
+        duplicates_cleaner=CloseDuplicatesCleaner(epsilon=1e-5),
         keep_infeasible=False,
     )
     algorithm.run()
@@ -90,3 +93,48 @@ def test_small_real_biobjective_nsag2():
         assert i.genes[0] == pytest.approx(i.genes[1], abs=0.2)
 
     assert len(final_population) == 200
+
+
+def test_small_real_biobjective_rnsga2():
+    """
+    Test a 2D real-valued problem with RNsga2 using reference points and an epsilon value.
+    The objectives are:
+        f1 = x^2 + y^2
+        f2 = (x-1)^2 + (y-1)^2
+    with x, y in [0, 1].
+
+    The true Pareto front is approximately: { (x,y) in [0,1]^2 : x ≈ y }.
+    We provide reference points uniformly along the line x=y in [0,1] and set epsilon=0.01.
+
+    """
+    # Define reference points uniformly along the line x=y
+    r = np.linspace(0, 0.2, 2)
+    reference_points = np.column_stack((r, r))  # shape (11, 2)
+
+    algorithm = RNsga2(
+        sampler=RandomSamplingFloat(min=0.0, max=1.0),
+        crossover=SimulatedBinaryCrossover(distribution_index=2),
+        mutation=GaussianMutation(gene_mutation_rate=0.1, sigma=0.05),
+        fitness_fn=fitness_biobjective,
+        constraints_fn=constraints_biobjective,
+        n_vars=2,  # Two variables: x and y
+        pop_size=200,
+        n_offsprings=200,
+        num_iterations=50,
+        mutation_rate=0.1,
+        crossover_rate=0.9,
+        duplicates_cleaner=CloseDuplicatesCleaner(epsilon=1e-8),
+        keep_infeasible=False,
+        reference_points=reference_points,
+        epsilon=1,
+    )
+    algorithm.run()
+
+    final_population = algorithm.population
+    best = final_population.best
+    for individual in best:
+        # Since the Pareto front is approximately x = y, test that gene[0] is near gene[1]
+        assert individual.genes[0] == pytest.approx(individual.genes[1], abs=0.2)
+
+    # Also, the final population size should be exactly 200.
+    assert len(final_population) == 200
diff --git a/uv.lock b/uv.lock
index 15713e3..35c21bf 100644
--- a/uv.lock
+++ b/uv.lock
@@ -689,7 +689,7 @@ wheels = [
 
 [[package]]
 name = "pymoors"
-version = "0.1.1rc1"
+version = "0.1.1"
 source = { editable = "." }
 dependencies = [
     { name = "numpy" },