Evolution results in "interpolations do not match"

[felixhekhorn]

Evolving e.g. HERA_NC_300GEV_EP-SIGMARED results in interpolations do not match raised from here:

pineappl/pineappl/src/grid.rs

Lines 464 to 465 in a1614db

	if self.interpolations() != other.interpolations() {
	return Err(Error::General("interpolations do not match".to_owned()));

using the debug script from below I can generate the following output:

$ pineappl evolve ../test/HERA_NC_300GEV_EP_SIGMARED.pineappl.lz4  ../test/HERA_NC_300GEV_EP_SIGMARED.tar ../test/fk.pineappl.lz4 NNPDF40_nnlo_as_01180
scales0: [Some(2.7224999999999997), None], scales1: [Some(300.00000000000006), None]
scales0: [Some(2.7224999999999997), None], scales1: [Some(45.00000000000001), None]
enter merge
merge
a 30000.00000000002, NaN, 2.571202181814284, 30000.00000000002, 50, 3, ApplGridH0, NoReweight, Lagrange
b 30000.00000000002, NaN, 2.571202181814284, 30000.00000000002, 50, 3, ApplGridH0, NoReweight, Lagrange
a==b -> false
a==b|min -> true
a==b|max -> true
a==b|rmin -> false
a==b|rmax -> true
a==b|nodes -> true
a==b|order -> true
a==b|map -> true
a==b|reweight_meth -> true
a==b|interp_meth -> true
a 0.00000019999999999999989, 0, 20.424947470398376, 1, 50, 4, ApplGridF2, ApplGridX, Lagrange
b 0.00000019999999999999989, 0, 20.424947470398376, 1, 50, 4, ApplGridF2, ApplGridX, Lagrange
a==b -> true
a==b|min -> true
a==b|max -> true
a==b|rmin -> true
a==b|rmax -> true
a==b|nodes -> true
a==b|order -> true
a==b|map -> true
a==b|reweight_meth -> true
a==b|interp_meth -> true
Error: interpolations do not match

(rmin is "real" min as saved in the struct and not a parsed version thereof)

so we are failing, because we set an interpolation bound to NaN instead of repeating the entry. I haven't dug deeper yet, but wanted to clear my mind before leaving the office and open already the issue, just in case someone already knows what to do (actually, let me tag @Radonirinaunimi and @cschwan to raise a bit the priority - sorry guys 😇 )

Things I don't understand:

• why is this a problem now (and not before)?
• I understood from @ecole41 + @andrpie this is only a problem for some datasets - why?
• why are there 3 Q2 printed? 45 + 300 + 30000 (all appear in the dataset)

If necessary I can provide the test data @ecole41 provided to me in the first place.

For the moment, I believe this is a PineAPPL bug and not an EKO bug.

---

diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs
index af45d6f0..3db82a11 100644
--- a/pineappl/src/grid.rs
+++ b/pineappl/src/grid.rs
@@ -460,6 +460,21 @@ impl Grid {
         if self.kinematics() != other.kinematics() {
             return Err(Error::General("kinematics do not match".to_owned()));
         }
+        println!("merge");
+        for (a,b) in self.interpolations().into_iter().zip(other.interpolations()) {
+            println!("a {}, {}, {}, {}, {}, {}, {:?}, {:?}, {:?}", a.min(),a.rmin(), a.rmax(), a.max(), a.nodes(), a.order(), a.map(), a.reweight_meth(), a.interp_meth());
+            println!("b {}, {}, {}, {}, {}, {}, {:?}, {:?}, {:?}", b.min(),b.rmin(), b.rmax(), b.max(), b.nodes(), b.order(), b.map(), b.reweight_meth(), a.interp_meth());
+            println!("a==b -> {}", a == b);
+            println!("a==b|min -> {}", a.min() == b.min());
+            println!("a==b|max -> {}", a.max() == b.max());
+            println!("a==b|rmin -> {}", a.rmin() == b.rmin());
+            println!("a==b|rmax -> {}", a.rmax() == b.rmax());
+            println!("a==b|nodes -> {}", a.nodes() == b.nodes());
+            println!("a==b|order -> {}", a.order() == b.order());
+            println!("a==b|map -> {}", a.map() == b.map());
+            println!("a==b|reweight_meth -> {}", a.reweight_meth() == b.reweight_meth());
+            println!("a==b|interp_meth -> {}", a.interp_meth() == b.interp_meth());
+        }
         // TODO: relax check if subgrid types don't use interpolation
         if self.interpolations() != other.interpolations() {
             return Err(Error::General("interpolations do not match".to_owned()));
@@ -1298,7 +1313,9 @@ impl Grid {
                 reference: self.reference.clone(),
             };
 
+            println!("scales0: {:?}, scales1: {:?}", scales0, scales1);
             if let Some(result) = &mut result {
+                println!("enter merge");
                 result.merge(evolved_slice)?;
             } else {
                 result = Some(evolved_slice);
diff --git a/pineappl/src/interpolation.rs b/pineappl/src/interpolation.rs
index 02211279..5cb564bd 100644
--- a/pineappl/src/interpolation.rs
+++ b/pineappl/src/interpolation.rs
@@ -62,7 +62,7 @@ fn lagrange_weights(i: usize, n: usize, u: f64) -> f64 {
 
 /// TODO
 #[repr(C)]
-#[derive(Clone, Copy, Deserialize, Eq, PartialEq, Serialize)]
+#[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize,)]
 pub enum ReweightMeth {
     /// TODO
     ApplGridX,
@@ -82,7 +82,7 @@ pub enum Map {
 
 /// TODO
 #[repr(C)]
-#[derive(Clone, Copy, Deserialize, Eq, PartialEq, Serialize)]
+#[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize)]
 pub enum InterpMeth {
     /// TODO
     Lagrange,
@@ -251,6 +251,9 @@ impl Interp {
         self.map_y_to_x(self.min).max(self.map_y_to_x(self.max))
     }
 
+    pub fn rmax(&self) -> f64 { self.max }
+    pub fn rmin(&self) -> f64 { self.min }
+
     /// TODO
     #[must_use]
     pub const fn map(&self) -> Map {

[Radonirinaunimi]

@felixhekhorn Could you please upload the EKO somewhere?

[cschwan]

How did you create an EKO that resulted in min = NaN? It should be a number.

[felixhekhorn]

@felixhekhorn Could you please upload the EKO somewhere?

https://nextcloud.jyu.fi/index.php/s/2GQ5iM9Wzoj57bG - this was provided to me by @ecole41

How did you create an EKO that resulted in min = NaN? It should be a number.

• this should refer to a scale, right? that should never be NaN I agree (but whatever the user requested)
• recall that the eko content (i.e. the 4D operator itself) might be NaN (and for this specific dataset it will happen I think), due to being in a non-perturbative region, which may result in NaN. This is not a problem since these bins are then cut out at fit level

[cschwan]

Is this below $\Lambda_\text{QCD} $? In that case the interpolation doesn't work.

[Radonirinaunimi]

Is this below Λ QCD ? In that case the interpolation doesn't work.

Yes, the problem is the following line:

pineappl/pineappl/src/interpolation.rs

Lines 40 to 42 in a1614db

	pub fn ftau0(q2: f64) -> f64 {
	(q2 / 0.0625).ln().ln()
	}

which if q2 is below 0.0625 (as that HERA dataset do) results in those NaN. Therefore, the bit-for-bit comparison of the fields of Interp with PartialEq fails.

[Radonirinaunimi]

I am not sure how this was not an issue before but we should definitely not rely on PartialEq for Interp. Instead, we should do something like:

use float_cmp::approx_eq;

impl PartialEq for Interp {
    fn eq(&self, other: &Self) -> bool {
        self.nodes == other.nodes
            && self.order == other.order
            && self.reweight == other.reweight
            && self.map == other.map
            && self.interp_meth == other.interp_meth
            && approx_eq!(f64, self.min, other.min, ulps = 8)
            && approx_eq!(f64, self.max, other.max, ulps = 8)
    }
}

impl Eq for Interp {}

I will open a PR for this.

[felixhekhorn]

I am not sure how this was not an issue before

this is indeed strange, these bins with very low Q2 (< $\Lambda_{QCD}$) have always been there. We were always happy to just get some NaNs (which is a reasonable value within pQCD) and cut them out later ...

Anyway, if this solves the issue, i.e. we can live with a NaN in the struct, I'm fine

[Radonirinaunimi]

I would guess it might be due to some CPU optimizations that compares the NaNs the same, but we can't know for sure unless the people who managed to correctly compute FK tables recently check this explicitly in the computer where they were produced.

[felixhekhorn]

just for pedantic reasons (which just comes to my mind): the relevant scale for PineAPPL is not $\Lambda_{QCD}$, but rather the hardcoded scale $Q_{min}^2 = 0.0625$ GeV². Instead, the relevant scale for EKO is $\Lambda_{QCD}$, which however depends on external things, such as PTO, $\alpha_s(M_Z)$, or $n_f$. This specific dataset contains a bin at $Q^2 = 0.045$ GeV², which is below either scale (but is of course measurable just fine).