Implement Cubic 2d interpolation as seen in LHAPDF #5

scarlehoff · 2023-07-16T12:32:18Z

This would require some caveats since in practice LHAPDF is defining a hierarchy on the axes and I'm not 100% sure it is symmetric. I would like to check that.

TODO:

~~Implement a Grid->GridSlice method that exists only for D > 2~~ This would generate unnecessary copies
C-api + benchmark just like with alpha_s
Check whether the hierarchy of the interpolation is symmetric if it is or if it isn't document the functions accordingly
Change Cubic1d to Cubic<Ix?>

Leaving this here so I don't forget, although technically is not related to this PR and I guess I should move it to an issue:
The extrapolation question arises again. In 1D is clear, if you are trying to extrapolate you can error out but now you might extrapolate only in 1D so you might still want the interpolation in the other. This would requirepartons to implement interpolation-extrapolation routines which would basically repeat some of the code here, so we might want to discuss at some point how to deal with that.
It might still be partons calling the proper extrapolation function from ndinterp, but I believe the functions themselves should be implemented here. So partons take the decisions and ndinter performs the operations.

scarlehoff · 2023-07-18T07:59:01Z

It's slower than LHAPDF (about 20% slower).

There are two possible improvements: vectorization and caching the derivatives of the PDF with respect to x for every q at construction time but I though it was preferable to have something working (20% slower is not much when talking about PDF interpolation) and deal with a possible improvement later since both those things should be transparent to the user.

ndinterp/src/grid/cubic.rs

alecandido · 2023-07-18T15:28:45Z

It's slower than LHAPDF (about 20% slower).

Do you have any idea about why? In principle, C++ and Rust are fully equivalent (they might produce the same LLVM IR), and you took the implementation from there, so I wonder where is the bottleneck.

However, 20% is perfectly acceptable, so we could even keep as it is...

scarlehoff · 2023-07-18T15:37:10Z

I haven’t looked in detail. This time I took it from there but not bug for bug, I just had a general look… LHAPDF might be doing some aggressive caching of the coefficients of the interpolation. If this is the difference this would be cheap to do when the grid is first created.

alecandido · 2023-07-19T09:27:42Z

ndinterp/src/grid.rs

+/// Together with the trait [`ToDimension`] this struct allows to convert a `usize` into a
+/// `Dimension` from the `ndarray` crate.
+pub struct DimensionHelper<const D: usize>;
+
+/// See documentation of [`DimensionHelper`].
+pub trait ToDimension {
+    /// A `Dimension` type from the `ndarray` crate.
+    type Dim: Dimension;
+}
+
+impl ToDimension for DimensionHelper<1> {
+    type Dim = Ix1;
+}
+
+impl ToDimension for DimensionHelper<2> {
+    type Dim = Ix2;
+}


Ok, this is not as I expected, but at least the repeated part is quite minimal!

Thanks @scarlehoff for the workaround!

I wonder if you can give a blanket implementation for usize, and just overwrite for 1 and 2 (i.e. saying that by default is Ix0). At least, we would reduce the boilerplate at call site, since all usize would implement the ToDimension trait.

I take it back: it's not possible (yet)

https://rust-lang.github.io/rfcs/1210-impl-specialization.html

scarlehoff · 2023-07-19T11:47:48Z

Thanks @cschwan I'll try to understand your changes.

Btw, for some reason the newest version of rust makes the code even slower. Now it's 40% worse than LHAPDF (from 20% before).
I just tried updated with rustup to see whether that would make a difference... I was hoping the difference to go in the opposite direction though... (so I need to read the changelog to see what has changed there).

cschwan · 2023-07-19T12:57:55Z

Hm, that's strange. It should make any difference, because I shouldn't have changed any (acutal) code.

scarlehoff · 2023-07-19T13:11:29Z

Nono, independently of your changes, I updated rust and the whole thing became slower.

alecandido · 2023-07-19T13:32:50Z

Nono, independently of your changes, I updated rust and the whole thing became slower.

This is also weird, a regression on performances of such a small amount of code is unexpected...
Are you sure it is not a statistical fluctuation?

scarlehoff · 2023-07-20T11:20:13Z

Ok, now this is ready for review.

RE performance, no, I tried a previous version of rust and the performance was recovered. I've seen in the changelog there was a few changes regarding the M1 mac son it might be just a problem for my laptop.

Please let me know if the gridNd_to_slice1d methods are ok or are confusing. Thinking on future improvements they are nice because they will allow me to have a cache of the derivatives (by generating all necessary slices together with the loading of the PDF) but I wont implement this in the short term since feature-parity with LHAPDF is much more necessary.

scarlehoff · 2023-07-20T14:15:44Z

Everything seems to work well. The results seem symmetric upon changing x1 and x2 so nothing weird with that. We can merge when you want.

alecandido

I left a few comments, but it's essentially fine. Merge whenever you wish.

ndinterp/src/grid.rs

alecandido · 2023-07-20T14:58:47Z

ndinterp/src/grid/cubic.rs

+
+/// Cubic interpolation
+#[derive(Debug)]
+pub struct Cubic<const D: usize>


To make it (possibly) even simpler, you could make it a 1-element tuple.

This is exactly the newtype pattern, even though Rust is missing a dedicated mechanism for that.

However, it's substantially the same, it would be truely simpler only if you could avoid the generics boilerplate somehow...

ndinterp/src/grid.rs

implement 2d cubic interpolation

0ba3d85

scarlehoff changed the base branch from main to prepare_for_ndinterpolation July 16, 2023 12:32

scarlehoff mentioned this pull request Jul 16, 2023

Update the code for n-d #4

Merged

scarlehoff added 3 commits July 17, 2023 10:36

change cubic1d to cubic<D> ; thanks @cschwan!

2686b4a

implement check against LHAPDF

bad19f1

clean up 2d interpolation

aecbf7d

alecandido reviewed Jul 18, 2023

View reviewed changes

ndinterp/src/grid/cubic.rs Outdated Show resolved Hide resolved

Convert the use of Dimension to usize

c50c07a

alecandido reviewed Jul 19, 2023

View reviewed changes

exploit the D to remove the N

23cb61e

scarlehoff added 2 commits July 19, 2023 16:50

simplify 2d interpolation

b8f3296

create a toslice method

2085cba

cschwan approved these changes Jul 20, 2023

View reviewed changes

Base automatically changed from prepare_for_ndinterpolation to main July 20, 2023 13:42

alecandido approved these changes Jul 20, 2023

View reviewed changes

restrict visibilities

c765f4c

alecandido reviewed Jul 21, 2023

View reviewed changes

ndinterp/src/grid.rs Show resolved Hide resolved

remove innecessary traitg

73d1715

alecandido merged commit f0f6311 into main Sep 9, 2023

alecandido deleted the cubic2d_lha branch September 9, 2023 08:30

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement Cubic 2d interpolation as seen in LHAPDF #5

Implement Cubic 2d interpolation as seen in LHAPDF #5

scarlehoff commented Jul 16, 2023 •

edited

Loading

scarlehoff commented Jul 18, 2023

alecandido commented Jul 18, 2023

scarlehoff commented Jul 18, 2023

alecandido Jul 19, 2023

alecandido Jul 19, 2023

scarlehoff commented Jul 19, 2023

cschwan commented Jul 19, 2023

scarlehoff commented Jul 19, 2023

alecandido commented Jul 19, 2023 •

edited

Loading

scarlehoff commented Jul 20, 2023

scarlehoff commented Jul 20, 2023 •

edited

Loading

alecandido left a comment

alecandido Jul 20, 2023 •

edited

Loading

Implement Cubic 2d interpolation as seen in LHAPDF #5

Implement Cubic 2d interpolation as seen in LHAPDF #5

Conversation

scarlehoff commented Jul 16, 2023 • edited Loading

scarlehoff commented Jul 18, 2023

alecandido commented Jul 18, 2023

scarlehoff commented Jul 18, 2023

alecandido Jul 19, 2023

Choose a reason for hiding this comment

alecandido Jul 19, 2023

Choose a reason for hiding this comment

scarlehoff commented Jul 19, 2023

cschwan commented Jul 19, 2023

scarlehoff commented Jul 19, 2023

alecandido commented Jul 19, 2023 • edited Loading

scarlehoff commented Jul 20, 2023

scarlehoff commented Jul 20, 2023 • edited Loading

alecandido left a comment

Choose a reason for hiding this comment

alecandido Jul 20, 2023 • edited Loading

Choose a reason for hiding this comment

scarlehoff commented Jul 16, 2023 •

edited

Loading

alecandido commented Jul 19, 2023 •

edited

Loading

scarlehoff commented Jul 20, 2023 •

edited

Loading

alecandido Jul 20, 2023 •

edited

Loading