refactor: simplify convolution function signatures

- change the signature of the `convolve` method, simplifying the underlying logic. Input/output grid are not arguments anymore, and only the grid size is passed.
- adapting the testing to this new signature

Author: mmolari

packages/treetime-convolution/src/algos/algos.rs

@@ -64,11 +64,5 @@ impl ConvolutionAlgorithm {
 pub trait Algo: Send + Sync {
   fn name(&self) -> &'static str;
 
-  fn convolve(
-    &self,
-    input_grid: &Array1<f64>,
-    f_values: &Array1<f64>,
-    g_values: &Array1<f64>,
-    output_grid: &Array1<f64>,
-  ) -> Result<Array1<f64>, Report>;
+  fn convolve(&self, dx: f64, f_values: &Array1<f64>, g_values: &Array1<f64>) -> Result<Array1<f64>, Report>;
 }

packages/treetime-convolution/src/algos/ndarray_conv/ndarray_conv.rs

@@ -2,37 +2,13 @@ use crate::algos::algos::Algo;
 use eyre::Report;
 use ndarray::Array1;
 use ndarray_conv::{ConvExt, ConvMode, PaddingMode};
-use ndarray_interp::interp1d::Interp1DBuilder;
-use treetime_utils::ndarray::is_uniform_grid;
 
 /// Convolution using ndarray_conv library for uniform grids
-pub fn convolve_ndarray_conv(
-  input_grid: &Array1<f64>,
-  f_values: &Array1<f64>,
-  g_values: &Array1<f64>,
-  output_grid: &Array1<f64>,
-) -> Result<Array1<f64>, Report> {
-  debug_assert!(is_uniform_grid(input_grid), "input_grid must be uniform");
-  debug_assert!(is_uniform_grid(output_grid), "output_grid must be uniform");
-
-  let grid_spacing = input_grid[1] - input_grid[0];
-
+/// takes as input the grid spacing, and arrays with function values
+pub fn convolve_ndarray_conv(dx: f64, f_values: &Array1<f64>, g_values: &Array1<f64>) -> Result<Array1<f64>, Report> {
   let discrete_conv = f_values.conv(g_values, ConvMode::Full, PaddingMode::Zeros)?;
-  let continuous_conv = &discrete_conv * grid_spacing;
-
-  let conv_result_len = discrete_conv.len();
-  let conv_result_min = input_grid[0] + input_grid[0];
-  let conv_result_max = input_grid[input_grid.len() - 1] + input_grid[input_grid.len() - 1];
-  let conv_result_grid = Array1::linspace(conv_result_min, conv_result_max, conv_result_len);
-
-  let temp_interp = Interp1DBuilder::new(continuous_conv).x(conv_result_grid).build()?;
-
-  let mut result = Array1::zeros(output_grid.len());
-  for (i, &x) in output_grid.iter().enumerate() {
-    result[i] = temp_interp.interp_scalar(x).unwrap_or(0.0);
-  }
-
-  Ok(result)
+  let continuous_conv = &discrete_conv * dx;
+  Ok(continuous_conv)
 }
 
 pub struct NdarrayAlgo;
@@ -42,13 +18,7 @@ impl Algo for NdarrayAlgo {
     "ndarray-conv"
   }
 
-  fn convolve(
-    &self,
-    input_grid: &Array1<f64>,
-    f_values: &Array1<f64>,
-    g_values: &Array1<f64>,
-    output_grid: &Array1<f64>,
-  ) -> Result<Array1<f64>, Report> {
-    convolve_ndarray_conv(input_grid, f_values, g_values, output_grid)
+  fn convolve(&self, dx: f64, f_values: &Array1<f64>, g_values: &Array1<f64>) -> Result<Array1<f64>, Report> {
+    convolve_ndarray_conv(dx, f_values, g_values)
   }
 }

packages/treetime-convolution/src/algos/ndarray_conv_fft/ndarray_conv_fft.rs

@@ -2,36 +2,17 @@ use crate::algos::algos::Algo;
 use eyre::Report;
 use ndarray::Array1;
 use ndarray_conv::{ConvFFTExt, ConvMode, PaddingMode};
-use ndarray_interp::interp1d::Interp1DBuilder;
-use treetime_utils::ndarray::is_uniform_grid;
 
+/// Convolution using ndarray_conv library with FFT for uniform grids.
+/// Takes as input the grid spacing, and arrays with function values
 pub fn convolve_ndarray_conv_fft(
-  input_grid: &Array1<f64>,
+  dx: f64,
   f_values: &Array1<f64>,
   g_values: &Array1<f64>,
-  output_grid: &Array1<f64>,
 ) -> Result<Array1<f64>, Report> {
-  debug_assert!(is_uniform_grid(input_grid), "input_grid must be uniform");
-  debug_assert!(is_uniform_grid(output_grid), "output_grid must be uniform");
-
-  let grid_spacing = input_grid[1] - input_grid[0];
-
   let discrete_conv = f_values.conv_fft(g_values, ConvMode::Full, PaddingMode::Zeros)?;
-  let continuous_conv = &discrete_conv * grid_spacing;
-
-  let conv_result_len = discrete_conv.len();
-  let conv_result_min = input_grid[0] + input_grid[0];
-  let conv_result_max = input_grid[input_grid.len() - 1] + input_grid[input_grid.len() - 1];
-  let conv_result_grid = Array1::linspace(conv_result_min, conv_result_max, conv_result_len);
-
-  let temp_interp = Interp1DBuilder::new(continuous_conv).x(conv_result_grid).build()?;
-
-  let mut result = Array1::zeros(output_grid.len());
-  for (i, &x) in output_grid.iter().enumerate() {
-    result[i] = temp_interp.interp_scalar(x).unwrap_or(0.0);
-  }
-
-  Ok(result)
+  let continuous_conv = &discrete_conv * dx;
+  Ok(continuous_conv)
 }
 
 /// Ndarray FFT-based convolution algorithm
@@ -42,13 +23,7 @@ impl Algo for NdarrayConvFftAlgo {
     "ndarray-conv-fft"
   }
 
-  fn convolve(
-    &self,
-    input_grid: &Array1<f64>,
-    f_values: &Array1<f64>,
-    g_values: &Array1<f64>,
-    output_grid: &Array1<f64>,
-  ) -> Result<Array1<f64>, Report> {
-    convolve_ndarray_conv_fft(input_grid, f_values, g_values, output_grid)
+  fn convolve(&self, dx: f64, f_values: &Array1<f64>, g_values: &Array1<f64>) -> Result<Array1<f64>, Report> {
+    convolve_ndarray_conv_fft(dx, f_values, g_values)
   }
 }

packages/treetime-convolution/src/algos/riemann/riemann.rs

@@ -1,28 +1,15 @@
 use crate::algos::algos::Algo;
 use eyre::Report;
-use itertools::izip;
 use ndarray::Array1;
-use ndarray_interp::interp1d::Interp1DBuilder;
 
-/// Convolution using Riemann sum integration
-pub fn convolve_riemann(
-  input_grid: &Array1<f64>,
-  f_values: &Array1<f64>,
-  g_values: &Array1<f64>,
-  output_grid: &Array1<f64>,
-) -> Result<Array1<f64>, Report> {
-  let g_interp = Interp1DBuilder::new(g_values.view()).x(input_grid.view()).build()?;
+/// Riemann convolution algorithm on uniform grids
+pub fn convolve_riemann(dx: f64, f_values: &Array1<f64>, g_values: &Array1<f64>) -> Result<Array1<f64>, Report> {
+  let mut result = Array1::zeros(f_values.len() + g_values.len() - 1);
 
-  let ds = input_grid[1] - input_grid[0];
-
-  let mut result = Array1::zeros(output_grid.len());
-  for (output_idx, &x_eval) in output_grid.iter().enumerate() {
-    let mut sum = 0.0;
-    for (&x_input, &f_at_x_input) in izip!(input_grid, f_values) {
-      let g_at_shifted = g_interp.interp_scalar(x_eval - x_input).unwrap_or(0.0);
-      sum += f_at_x_input * g_at_shifted;
+  for (i, &f_val) in f_values.iter().enumerate() {
+    for (j, &g_val) in g_values.iter().enumerate() {
+      result[i + j] += f_val * g_val * dx;
     }
-    result[output_idx] = sum * ds;
   }
 
   Ok(result)
@@ -35,13 +22,7 @@ impl Algo for RiemannAlgo {
     "riemann"
   }
 
-  fn convolve(
-    &self,
-    input_grid: &Array1<f64>,
-    f_values: &Array1<f64>,
-    g_values: &Array1<f64>,
-    output_grid: &Array1<f64>,
-  ) -> Result<Array1<f64>, Report> {
-    convolve_riemann(input_grid, f_values, g_values, output_grid)
+  fn convolve(&self, dx: f64, f_values: &Array1<f64>, g_values: &Array1<f64>) -> Result<Array1<f64>, Report> {
+    convolve_riemann(dx, f_values, g_values)
   }
 }

packages/treetime-convolution/src/testing/framework/test_case.rs

@@ -14,8 +14,4 @@ pub trait TestCase: Clone + Send + Sync + Serialize {
   fn input_grid_domain(&self) -> (f64, f64);
 
   fn input_grid_n_points(&self) -> usize;
-
-  fn output_grid_domain(&self) -> (f64, f64);
-
-  fn output_grid_n_points(&self) -> usize;
 }

packages/treetime-convolution/src/testing/run.rs

@@ -262,19 +262,20 @@ fn run_test<S: TestSuite>(
   let input_grid_n_points = test_case.input_grid_n_points();
   let input_grid = Array1::linspace(input_grid_min, input_grid_max, input_grid_n_points);
 
-  let (output_grid_min, output_grid_max) = test_case.output_grid_domain();
-  let output_grid_n_points = test_case.output_grid_n_points();
-  let output_grid = Array1::linspace(output_grid_min, output_grid_max, output_grid_n_points);
+  let (evaluation_grid_min, evaluation_grid_max) = (input_grid_min * 2.0, input_grid_max * 2.0);
+  let evaluation_grid_n_points = 2 * input_grid_n_points - 1;
+  let evaluation_grid = Array1::linspace(evaluation_grid_min, evaluation_grid_max, evaluation_grid_n_points);
 
   let f_values = suite.create_f(test_case, &input_grid)?;
   let g_values = suite.create_g(test_case, &input_grid)?;
 
-  let actual_values = algo.convolve(&input_grid, &f_values, &g_values, &output_grid)?;
-  let expected_values = suite.analytical_convolution(test_case, &output_grid)?;
+  let dx = input_grid[1] - input_grid[0];
+  let actual_values = algo.convolve(dx, &f_values, &g_values)?;
+  let expected_values = suite.analytical_convolution(test_case, &evaluation_grid)?;
 
   let execution_time = start_time.elapsed().as_secs_f64() * 1000.0;
 
-  let metrics = ConvolutionMetrics::new(&output_grid, &actual_values, &expected_values, execution_time)?;
+  let metrics = ConvolutionMetrics::new(&evaluation_grid, &actual_values, &expected_values, execution_time)?;
 
   Ok(TestResult {
     algorithm: algo.name().to_owned(),
@@ -284,7 +285,7 @@ fn run_test<S: TestSuite>(
     f_y_values: f_values,
     g_x_values: input_grid,
     g_y_values: g_values,
-    evaluation_grid: output_grid,
+    evaluation_grid,
     actual_values,
     expected_values,
     metrics,

packages/treetime-convolution/src/testing/test_suites/exponential.rs

@@ -51,8 +51,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 2.0,
         input_grid_domain: (-1.0, 10.0),
         input_grid_n_points: 1101,
-        output_grid_domain: (-1.0, 10.0),
-        output_grid_n_points: 1101,
       },
       ExponentialTestCase {
         name: "moderate_coarse_grid".to_owned(),
@@ -65,8 +63,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 0.8,
         input_grid_domain: (0.0, 20.0),
         input_grid_n_points: 201,
-        output_grid_domain: (0.0, 30.0),
-        output_grid_n_points: 301,
       },
       ExponentialTestCase {
         name: "tight_truncation".to_owned(),
@@ -80,8 +76,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 2.0,
         input_grid_domain: (0.0, 5.0),
         input_grid_n_points: 501,
-        output_grid_domain: (0.0, 6.0),
-        output_grid_n_points: 601,
       },
       ExponentialTestCase {
         name: "baseline_distinct_rates".to_owned(),
@@ -94,8 +88,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 2.0,
         input_grid_domain: (0.0, 10.0),
         input_grid_n_points: 1001,
-        output_grid_domain: (0.0, 17.0),
-        output_grid_n_points: 1701,
       },
       ExponentialTestCase {
         name: "equal_rates_limit".to_owned(),
@@ -107,8 +99,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 1.5,
         input_grid_domain: (0.0, 12.0),
         input_grid_n_points: 1201,
-        output_grid_domain: (0.0, 24.0),
-        output_grid_n_points: 2401,
       },
       ExponentialTestCase {
         name: "near_equal_rates".to_owned(),
@@ -121,8 +111,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 1.0 + 1e-6,
         input_grid_domain: (0.0, 30.0),
         input_grid_n_points: 15001,
-        output_grid_domain: (0.0, 60.0),
-        output_grid_n_points: 30001,
       },
       ExponentialTestCase {
         name: "fast_slow_decay".to_owned(),
@@ -135,8 +123,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 0.1,
         input_grid_domain: (0.0, 200.0),
         input_grid_n_points: 200001,
-        output_grid_domain: (0.0, 202.0),
-        output_grid_n_points: 202001,
       },
       ExponentialTestCase {
         name: "slow_fast_decay".to_owned(),
@@ -149,8 +135,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 10.0,
         input_grid_domain: (0.0, 200.0),
         input_grid_n_points: 200001,
-        output_grid_domain: (0.0, 202.0),
-        output_grid_n_points: 202001,
       },
       ExponentialTestCase {
         name: "fine_grid_reference".to_owned(),
@@ -163,8 +147,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 0.8,
         input_grid_domain: (0.0, 20.0),
         input_grid_n_points: 10001,
-        output_grid_domain: (0.0, 40.0),
-        output_grid_n_points: 20001,
       },
       ExponentialTestCase {
         name: "long_horizon_tail_precision".to_owned(),
@@ -176,8 +158,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 0.04,
         input_grid_domain: (0.0, 750.0),
         input_grid_n_points: 75001,
-        output_grid_domain: (0.0, 1350.0),
-        output_grid_n_points: 135001,
       },
       ExponentialTestCase {
         name: "large_range_underflow_guard".to_owned(),
@@ -190,8 +170,6 @@ impl TestSuite for ExponentialTestSuite {
         b: 0.7,
         input_grid_domain: (0.0, 400.0),
         input_grid_n_points: 20001,
-        output_grid_domain: (0.0, 800.0),
-        output_grid_n_points: 40001,
       },
       ExponentialTestCase {
         name: "extreme_near_equality".to_owned(),
@@ -204,8 +182,8 @@ impl TestSuite for ExponentialTestSuite {
         b: 1.0 + 1e-12,
         input_grid_domain: (0.0, 80.0),
         input_grid_n_points: 16001,
-        output_grid_domain: (0.0, 160.0),
-        output_grid_n_points: 32001,
+        
+        
       },
       ExponentialTestCase {
         name: "very_fast_decays".to_owned(),
@@ -217,8 +195,8 @@ impl TestSuite for ExponentialTestSuite {
         b: 30.0,
         input_grid_domain: (0.0, 0.5),
         input_grid_n_points: 1001,
-        output_grid_domain: (0.0, 1.0),
-        output_grid_n_points: 2001,
+        
+        
       },
     ]
   }
@@ -236,8 +214,8 @@ pub struct ExponentialTestCase {
   pub b: f64,
   pub input_grid_domain: (f64, f64),
   pub input_grid_n_points: usize,
-  pub output_grid_domain: (f64, f64),
-  pub output_grid_n_points: usize,
+  
+  
 }
 
 impl TestCase for ExponentialTestCase {
@@ -269,11 +247,4 @@ impl TestCase for ExponentialTestCase {
     self.input_grid_n_points
   }
 
-  fn output_grid_domain(&self) -> (f64, f64) {
-    self.output_grid_domain
-  }
-
-  fn output_grid_n_points(&self) -> usize {
-    self.output_grid_n_points
-  }
 }