Remove empty entries from luminosity function

pineappl/src/empty_subgrid.rs

@@ -1,6 +1,7 @@
 //! TODO
 
 use super::grid::Ntuple;
+use super::sparse_array3::SparseArray3;
 use super::subgrid::{Subgrid, SubgridEnum};
 use either::Either;
 use serde::{Deserialize, Serialize};
@@ -66,6 +67,10 @@ impl Subgrid for EmptySubgridV1 {
     fn iter(&self) -> Box<dyn Iterator<Item = ((usize, usize, usize), &f64)> + '_> {
         Box::new(iter::empty())
     }
+
+    fn export(&self) -> Cow<SparseArray3<f64>> {
+        unreachable!();
+    }
 }
 
 #[cfg(test)]

pineappl/src/grid.rs

@@ -12,7 +12,7 @@ use float_cmp::approx_eq;
 use git_version::git_version;
 use itertools::Itertools;
 use lz_fear::{framed::DecompressionError::WrongMagic, LZ4FrameReader};
-use ndarray::{Array3, Dimension};
+use ndarray::{s, Array3, Dimension};
 use rustc_hash::FxHashMap;
 use serde::{Deserialize, Serialize};
 use std::borrow::Cow;
@@ -940,7 +940,8 @@ impl Grid {
             .key_values()
             .map_or(true, |map| map["initial_state_1"] == map["initial_state_2"])
         {
-            self.symmetrize();
+            self.symmetrize_lumi();
+            self.optimize_lumi();
         }
 
         for subgrid in self.subgrids.iter_mut() {
@@ -967,7 +968,39 @@ impl Grid {
         }
     }
 
-    fn symmetrize(&mut self) {
+    fn optimize_lumi(&mut self) {
+        let mut keep_lumi_indices = vec![];
+        let mut new_lumi_entries = vec![];
+
+        for (lumi, entry) in self.lumi.iter().enumerate() {
+            let slice = self.subgrids.slice(s![.., .., lumi]);
+
+            if !slice.iter().all(|subgrid| subgrid.is_empty()) {
+                keep_lumi_indices.push(lumi);
+                new_lumi_entries.push(entry.clone());
+            }
+        }
+
+        let new_subgrids = Array3::from_shape_fn(
+            (
+                self.orders.len(),
+                self.bin_info().bins(),
+                keep_lumi_indices.len(),
+            ),
+            |(order, bin, new_lumi)| {
+                mem::replace(
+                    &mut self.subgrids[[order, bin, keep_lumi_indices[new_lumi]]],
+                    EmptySubgridV1::default().into(),
+                )
+            },
+        );
+
+        self.lumi = new_lumi_entries;
+        self.subgrids = new_subgrids;
+    }
+
+    // TODO: simplify the method, because `optimize_lumi` already removes empty entries
+    fn symmetrize_lumi(&mut self) {
         let mut indices: Vec<usize> = (0..self.lumi.len()).rev().collect();
         let mut pairs: Vec<(usize, usize)> = Vec::new();
         let mut not_symmetrized: Vec<usize> = Vec::new();

pineappl/src/import_only_subgrid.rs

@@ -187,6 +187,10 @@ impl Subgrid for ImportOnlySubgridV1 {
     fn iter(&self) -> Box<dyn Iterator<Item = ((usize, usize, usize), &f64)> + '_> {
         Box::new(self.array.indexed_iter())
     }
+
+    fn export(&self) -> Cow<SparseArray3<f64>> {
+        Cow::Borrowed(&self.array)
+    }
 }
 
 impl From<&LagrangeSubgridV2> for ImportOnlySubgridV1 {

pineappl/src/lagrange_subgrid.rs

@@ -353,6 +353,10 @@ impl Subgrid for LagrangeSubgridV1 {
         .into()
     }
 
+    fn export(&self) -> Cow<SparseArray3<f64>> {
+        todo!();
+    }
+
     fn iter(&self) -> Box<dyn Iterator<Item = ((usize, usize, usize), &f64)> + '_> {
         self.grid.as_ref().map_or_else(
             || Box::new(iter::empty()) as Box<dyn Iterator<Item = ((usize, usize, usize), &f64)>>,
@@ -732,6 +736,10 @@ impl Subgrid for LagrangeSubgridV2 {
             },
         )
     }
+
+    fn export(&self) -> Cow<SparseArray3<f64>> {
+        todo!();
+    }
 }
 
 /// Subgrid which uses Lagrange-interpolation, but also stores its contents in a space-efficient
@@ -974,6 +982,10 @@ impl Subgrid for LagrangeSparseSubgridV1 {
     fn iter(&self) -> Box<dyn Iterator<Item = ((usize, usize, usize), &f64)> + '_> {
         Box::new(self.array.indexed_iter())
     }
+
+    fn export(&self) -> Cow<SparseArray3<f64>> {
+        Cow::Borrowed(&self.array)
+    }
 }
 
 impl From<&LagrangeSubgridV1> for LagrangeSparseSubgridV1 {

pineappl/src/ntuple_subgrid.rs

@@ -1,6 +1,7 @@
 //! Provides an implementation of the `Grid` trait with n-tuples.
 
 use super::grid::Ntuple;
+use super::sparse_array3::SparseArray3;
 use super::subgrid::{Subgrid, SubgridEnum};
 use either::Either;
 use serde::{Deserialize, Serialize};
@@ -90,6 +91,10 @@ impl Subgrid for NtupleSubgridV1 {
     fn iter(&self) -> Box<dyn Iterator<Item = ((usize, usize, usize), &f64)>> {
         unimplemented!();
     }
+
+    fn export(&self) -> Cow<SparseArray3<f64>> {
+        unimplemented!();
+    }
 }
 
 #[cfg(test)]

pineappl/src/sparse_array3.rs

@@ -8,7 +8,7 @@ use std::slice::{Iter, IterMut};
 
 /// Struct for a sparse three-dimensional array, which is optimized for the sparsity of
 /// interpolation grids.
-#[derive(Deserialize, Serialize)]
+#[derive(Clone, Deserialize, Serialize)]
 pub struct SparseArray3<T> {
     entries: Vec<T>,
     indices: Vec<(usize, usize)>,

pineappl/src/subgrid.rs

@@ -5,6 +5,7 @@ use super::grid::Ntuple;
 use super::import_only_subgrid::ImportOnlySubgridV1;
 use super::lagrange_subgrid::{LagrangeSparseSubgridV1, LagrangeSubgridV1, LagrangeSubgridV2};
 use super::ntuple_subgrid::NtupleSubgridV1;
+use super::sparse_array3::SparseArray3;
 use either::Either;
 use enum_dispatch::enum_dispatch;
 use serde::{Deserialize, Serialize};
@@ -78,6 +79,9 @@ pub trait Subgrid {
     /// Fill the q2-slice with index `q2_slice` into `grid`.
     fn fill_q2_slice(&self, q2_slice: usize, grid: &mut [f64]);
 
+    /// Export the subgrid as a [`SparseArray3`].
+    fn export(&self) -> Cow<SparseArray3<f64>>;
+
     /// Assumes that the initial states for this grid are the same and uses this to optimize the
     /// grid by getting rid of almost half of the entries.
     fn symmetrize(&mut self);