Generalize static node value detection

pineappl/src/empty_subgrid.rs

@@ -1,7 +1,7 @@
 //! TODO
 
 use super::interpolation::Interp;
-use super::subgrid::{Mu2, Stats, Subgrid, SubgridEnum, SubgridIndexedIter};
+use super::subgrid::{Stats, Subgrid, SubgridEnum, SubgridIndexedIter};
 use serde::{Deserialize, Serialize};
 use std::iter;
 
@@ -47,9 +47,7 @@ impl Subgrid for EmptySubgridV1 {
         }
     }
 
-    fn static_scale(&self) -> Option<Mu2> {
-        None
-    }
+    fn optimize_static_nodes(&mut self) {}
 }
 
 #[cfg(test)]
@@ -66,6 +64,7 @@ mod tests {
         subgrid.merge(&EmptySubgridV1.into(), None);
         subgrid.scale(2.0);
         subgrid.symmetrize(1, 2);
+        subgrid.optimize_static_nodes();
         assert_eq!(
             subgrid.stats(),
             Stats {
@@ -76,7 +75,6 @@ mod tests {
                 bytes_per_value: 0,
             }
         );
-        assert_eq!(subgrid.static_scale(), None);
     }
 
     #[test]

pineappl/src/grid.rs

@@ -850,23 +850,18 @@ impl Grid {
         }
     }
 
-    fn optimize_subgrid_type(&mut self, static_scale_detection: bool) {
+    fn optimize_subgrid_type(&mut self, optimize_static_nodes: bool) {
         for subgrid in &mut self.subgrids {
             match subgrid {
                 // replace empty subgrids of any type with `EmptySubgridV1`
                 _ if subgrid.is_empty() => {
                     *subgrid = EmptySubgridV1.into();
                 }
                 _ => {
-                    // TODO: this requires a `pub(crate)` in `InterpSubgridV1`; we should
-                    // replace this with a method
-                    if !static_scale_detection {
-                        if let SubgridEnum::InterpSubgridV1(subgrid) = subgrid {
-                            // disable static-scale detection
-                            subgrid.static_q2 = -1.0;
-                        }
+                    if optimize_static_nodes {
+                        subgrid.optimize_static_nodes();
                     }
-
+                    // TODO: check if we should remove this
                     *subgrid = ImportSubgridV1::from(&*subgrid).into();
                 }
             }

pineappl/src/import_subgrid.rs

@@ -3,8 +3,8 @@
 use super::evolution::EVOLVE_INFO_TOL_ULPS;
 use super::interpolation::Interp;
 use super::packed_array::PackedArray;
-use super::subgrid::{Mu2, Stats, Subgrid, SubgridEnum, SubgridIndexedIter};
-use float_cmp::{approx_eq, assert_approx_eq};
+use super::subgrid::{Stats, Subgrid, SubgridEnum, SubgridIndexedIter};
+use float_cmp::approx_eq;
 use itertools::izip;
 use serde::{Deserialize, Serialize};
 use std::mem;
@@ -130,17 +130,7 @@ impl Subgrid for ImportSubgridV1 {
         }
     }
 
-    fn static_scale(&self) -> Option<Mu2> {
-        if let &[static_scale] = self.node_values()[0].as_slice() {
-            Some(Mu2 {
-                ren: static_scale,
-                fac: static_scale,
-                frg: -1.0,
-            })
-        } else {
-            None
-        }
-    }
+    fn optimize_static_nodes(&mut self) {}
 }
 
 impl From<&SubgridEnum> for ImportSubgridV1 {
@@ -161,32 +151,18 @@ impl From<&SubgridEnum> for ImportSubgridV1 {
             },
         );
 
-        let mut new_node_values: Vec<_> = subgrid
+        let new_node_values: Vec<_> = subgrid
             .node_values()
             .iter()
             .zip(&ranges)
             .map(|(values, range)| values[range.clone()].to_vec())
             .collect();
-        let static_scale = if let Some(Mu2 { ren, fac, frg }) = subgrid.static_scale() {
-            assert_approx_eq!(f64, ren, fac, ulps = 4);
-            assert_approx_eq!(f64, frg, -1.0, ulps = 4);
-            new_node_values[0] = vec![fac];
-            true
-        } else {
-            false
-        };
 
         let mut array = PackedArray::new(new_node_values.iter().map(Vec::len).collect());
 
         for (mut indices, value) in subgrid.indexed_iter() {
-            for (idx, (index, range)) in indices.iter_mut().zip(&ranges).enumerate() {
-                // TODO: generalize static scale detection
-                *index = if static_scale && idx == 0 {
-                    // if there's a static scale we want every value to be added to same grid point
-                    0
-                } else {
-                    *index - range.start
-                };
+            for (index, range) in indices.iter_mut().zip(&ranges) {
+                *index = *index - range.start;
             }
 
             array[indices.as_slice()] += value;

pineappl/src/interp_subgrid.rs

@@ -2,7 +2,7 @@
 
 use super::interpolation::{self, Interp};
 use super::packed_array::PackedArray;
-use super::subgrid::{Mu2, Stats, Subgrid, SubgridEnum, SubgridIndexedIter};
+use super::subgrid::{Stats, Subgrid, SubgridEnum, SubgridIndexedIter};
 use float_cmp::approx_eq;
 use serde::{Deserialize, Serialize};
 use std::mem;
@@ -12,7 +12,7 @@ use std::mem;
 pub struct InterpSubgridV1 {
     array: PackedArray<f64>,
     interps: Vec<Interp>,
-    pub(crate) static_q2: f64,
+    static_nodes: Vec<Option<f64>>,
 }
 
 impl InterpSubgridV1 {
@@ -22,7 +22,7 @@ impl InterpSubgridV1 {
         Self {
             array: PackedArray::new(interps.iter().map(Interp::nodes).collect()),
             interps: interps.to_vec(),
-            static_q2: 0.0,
+            static_nodes: vec![Some(-1.0); interps.len()],
         }
     }
 }
@@ -32,14 +32,14 @@ impl Subgrid for InterpSubgridV1 {
         debug_assert_eq!(interps.len(), ntuple.len());
 
         if interpolation::interpolate(interps, ntuple, weight, &mut self.array) {
-            // TODO: make this more general
-            let q2 = ntuple[0];
-            if self.static_q2 == 0.0 {
-                self.static_q2 = q2;
-            } else if !approx_eq!(f64, self.static_q2, -1.0, ulps = 4)
-                && !approx_eq!(f64, self.static_q2, q2, ulps = 4)
-            {
-                self.static_q2 = -1.0;
+            for (value, previous_node) in ntuple.iter().zip(&mut self.static_nodes) {
+                if let Some(previous_value) = previous_node {
+                    if *previous_value < 0.0 {
+                        *previous_value = *value;
+                    } else if !approx_eq!(f64, *previous_value, *value, ulps = 4) {
+                        *previous_node = None;
+                    }
+                }
             }
         }
     }
@@ -110,12 +110,40 @@ impl Subgrid for InterpSubgridV1 {
         }
     }
 
-    fn static_scale(&self) -> Option<Mu2> {
-        (self.static_q2 > 0.0).then_some(Mu2 {
-            ren: self.static_q2,
-            fac: self.static_q2,
-            frg: -1.0,
-        })
+    fn optimize_static_nodes(&mut self) {
+        let mut new_array = PackedArray::new(
+            self.array
+                .shape()
+                .iter()
+                .zip(&self.static_nodes)
+                .map(|(&dim, static_node)| if static_node.is_some() { 1 } else { dim })
+                .collect(),
+        );
+
+        for (mut index, value) in self.array.indexed_iter() {
+            for (idx, static_node) in index.iter_mut().zip(&self.static_nodes) {
+                if static_node.is_some() {
+                    *idx = 0;
+                }
+            }
+            new_array[index.as_slice()] += value;
+        }
+
+        self.array = new_array;
+
+        for (static_node, interp) in self.static_nodes.iter_mut().zip(&mut self.interps) {
+            if let &mut Some(value) = static_node {
+                *interp = Interp::new(
+                    value,
+                    value,
+                    1,
+                    0,
+                    interp.reweight_meth(),
+                    interp.map(),
+                    interp.interp_meth(),
+                );
+            }
+        }
     }
 }
 
@@ -175,14 +203,6 @@ mod tests {
 
         assert!(!subgrid.is_empty());
         assert_eq!(subgrid.indexed_iter().count(), 4 * 4 * 4);
-        assert_eq!(
-            subgrid.static_scale(),
-            Some(Mu2 {
-                ren: 1000.0,
-                fac: 1000.0,
-                frg: -1.0
-            })
-        );
         assert_eq!(
             subgrid.stats(),
             Stats {
@@ -198,7 +218,6 @@ mod tests {
 
         assert!(!subgrid.is_empty());
         assert_eq!(subgrid.indexed_iter().count(), 2 * 4 * 4 * 4);
-        assert_eq!(subgrid.static_scale(), None);
         assert_eq!(
             subgrid.stats(),
             Stats {

pineappl/src/interpolation.rs

@@ -207,9 +207,13 @@ impl Interp {
     /// TODO
     #[must_use]
     pub fn node_values(&self) -> Vec<f64> {
-        (0..self.nodes)
-            .map(|node| self.map_y_to_x(self.gety(node)))
-            .collect()
+        if self.nodes == 1 {
+            vec![self.map_y_to_x(self.min)]
+        } else {
+            (0..self.nodes)
+                .map(|node| self.map_y_to_x(self.gety(node)))
+                .collect()
+        }
     }
 
     fn map_y_to_x(&self, y: f64) -> f64 {
@@ -249,6 +253,18 @@ impl Interp {
     pub const fn map(&self) -> Map {
         self.map
     }
+
+    /// TODO
+    #[must_use]
+    pub const fn interp_meth(&self) -> InterpMeth {
+        self.interp_meth
+    }
+
+    /// TODO
+    #[must_use]
+    pub const fn reweight_meth(&self) -> ReweightMeth {
+        self.reweight
+    }
 }
 
 /// TODO
@@ -682,7 +698,6 @@ mod tests {
 
         assert_eq!(node_values.len(), 1);
 
-        // TODO: the return value is not the one expected (90^2), because `deltay` is zero
-        assert!(node_values[0].is_nan());
+        assert_approx_eq!(f64, node_values[0], 90.0 * 90.0, ulps = 16);
     }
 }

pineappl/src/subgrid.rs

@@ -85,8 +85,8 @@ pub trait Subgrid {
     /// Return statistics for this subgrid.
     fn stats(&self) -> Stats;
 
-    /// Return the static (single) scale, if this subgrid has one.
-    fn static_scale(&self) -> Option<Mu2>;
+    /// TODO
+    fn optimize_static_nodes(&mut self);
 }
 
 /// Type to iterate over the non-zero contents of a subgrid. The tuple contains the indices of the