From 9234e395e9b76dff1406706eb9f7c06f68960bc7 Mon Sep 17 00:00:00 2001
From: Vince Buffalo <vsbuffalo@gmail.com>
Date: Thu, 8 Aug 2024 20:28:37 -0700
Subject: [PATCH] Fixed small bug in interpolation.

The previous implementation was not using a first position of zero.
---
 src/numeric.rs |  1 +
 src/recmap.rs  | 18 +++++++++++++-----
 2 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/src/numeric.rs b/src/numeric.rs
index d13f748..e0e5e87 100644
--- a/src/numeric.rs
+++ b/src/numeric.rs
@@ -149,6 +149,7 @@ where
     Tx: PartialOrd + ToPrimitive + Copy + Debug,
     Ty: ToPrimitive + NumCast + Copy + Debug,
 {
+    // dbg!((&x, &y, &x0));
     assert!(x.len() == y.len());
     let index = search_sorted(x, x0);
     match index {
diff --git a/src/recmap.rs b/src/recmap.rs
index 753946c..da08af5 100644
--- a/src/recmap.rs
+++ b/src/recmap.rs
@@ -111,7 +111,7 @@ pub struct RateMap {
     pub ends: Vec<Position>,
     /// The n rates in Morgans, between each genomic position.
     pub rates: Vec<RateFloat>,
-    /// The n+1 cumulative map lengths at each genomic position.
+    /// The n cumulative map lengths at each genomic position.
     pub map_pos: Vec<RateFloat>,
 }
 
@@ -436,9 +436,17 @@ impl RecMap {
             .get(name)
             .ok_or(RecMapError::NoChrom(name.to_string()))?;
         let ends = &rate_map.ends;
-        let interp_result = interp1d(&ends[0..ends.len() - 1], 
-                                     &rate_map.map_pos, 
-                                     position);
+        // dbg!((&ends, &rate_map.map_pos));
+        assert_eq!(ends.len(), rate_map.map_pos.len() + 1);
+
+        // Create a new vector with 0 at the beginning, followed by map_pos
+        let extended_map_pos = std::iter::once(&0.0)
+            .chain(rate_map.map_pos.iter())
+            .cloned()
+            .collect::<Vec<f64>>();
+
+        // NOTE: This has changed recently <!> Still under review.
+        let interp_result = interp1d(&ends, &extended_map_pos, position);
         let interpolated_map_pos =
             interp_result.ok_or(RecMapError::LookupOutOfBounds(name.to_string(), position))?;
         Ok(interpolated_map_pos)
@@ -669,7 +677,7 @@ mod tests {
             &chr1_map.map_pos,
             &to_morgans(vec![0.0, 5.5, 13., 34.65]),
             1e-3,
-            );
+        );
     }
 
     #[test]