Update the code for n-d #4

scarlehoff · 2023-07-06T08:48:21Z

I've updated the code so that it is prepared to handle n-d interpolation.

My idea is to have a Vector of 1d-vectors for the input points and a nd-array for the values of the grid. The reason why I've changed array to vector is that, since we know (or enforce during construction if need be) at least for LHAPDF that the input is sorted 1d arrays, using vectors allow to use some of the nice vector methods that are considerably faster than my implementation of the search in ndarray. Indeed, this version is now faster than LHAPDF :)

This is also much easier to generalize to n-dimensions.

Other changes that I've made is to comment out the imports to the scatter part of the code, because clippy was making me crazy.

ndinterp/src/grid.rs

ndinterp_capi/src/lib.rs

ndinterp/src/grid.rs

scarlehoff · 2023-07-06T11:43:18Z

Now the code should work for any number of dimensions. It would be nice to be able to somehow join the D: dimension and the const N: usize parameters (since given one you know the other) but I haven't figured out how to do it...

scarlehoff · 2023-07-06T12:15:14Z

-        let p0 = yl * (2. * t3 - 3. * t2 + 1.);
-        let p1 = yu * (-2. * t3 + 3. * t2);
-        let m0 = dydxl * (t3 - 2. * t2 + t);
+        let p0 = yl * (2.0f64.mul_add(t3, -3. * t2) + 1.);
+        let p1 = yu * (-2.0f64).mul_add(t3, 3. * t2);
+        let m0 = dydxl * (2.0f64.mul_add(-t2, t3) + t);

If I do cargo clippy --fix it produces changes like these. Does it really make sense to apply them? (why?) imho it just makes things difficult to read?

alecandido · 2023-07-06T12:47:30Z

If I do cargo clippy --fix it produces changes like these. Does it really make sense to apply them? (why?) imho it just makes things difficult to read?

I guess it's disambiguating the kind of operation, or something like. Do we have a reference to the specific lint?

scarlehoff · 2023-07-06T12:51:14Z

It points to this https://rust-lang.github.io/rust-clippy/master/index.html#/suboptimal_flops
but even if they were suboptimal, I don't think it's impact will be noticeable (and furthermore, this is the kind of optimisation that I would be disappointed if the compiler is not able to do by itself) and I personally think that makes the code less readable.

But if it is the "rustacean" way to do it I can learn to live with it x)

cschwan · 2023-07-06T13:28:54Z

It points to this https://rust-lang.github.io/rust-clippy/master/index.html#/suboptimal_flops but even if they were suboptimal, I don't think it's impact will be noticeable (and furthermore, this is the kind of optimisation that I would be disappointed if the compiler is not able to do by itself) and I personally think that makes the code less readable.

But if it is the "rustacean" way to do it I can learn to live with it x)

I think it's possible the compiler will not do it or will only do it under very specific circumstances because it would require reordering floating-point operations, which usually isn't done because FLOPS aren't associative.

On Intel CPUs mul_add should use FMA operations, see also here: https://doc.rust-lang.org/std/primitive.f64.html#method.mul_add.

scarlehoff · 2023-07-06T13:31:27Z

But then clippy is taking a way too aggressive stance!

cschwan · 2023-07-06T13:44:45Z

Well, remind yourself where its name comes from 😄. But if it annoys you, just disable the lint.

scarlehoff · 2023-07-06T14:01:12Z

Some of the lints I will remove/ignore because I feel they are a bit extreme. But this particular one I want to know what you guys usually do.

alecandido · 2023-07-07T17:26:59Z

Some of the lints I will remove/ignore because I feel they are a bit extreme. But this particular one I want to know what you guys usually do.

To be fair, it seems like I've never done enough math in a single line in Rust to invoke this lint.
On the other hand, I understand the optimization, and it's not stupid (maybe we could save the f64 suffix), but it's also hardware-specific, and most likely not critical.

I would try to benchmark once just for the sake of it (since benchmarks are already in place), and then you could discard in favor of simplicity.

scarlehoff · 2023-07-07T17:48:23Z

No, I'm sure it's clever. But a bit extreme.

I the benchmark the differences are within the variance.

alecandido · 2023-07-07T18:16:45Z

Ok, then you could forget that lint. I bet it won't affect much even in the 2D case.

scarlehoff · 2023-07-07T18:17:58Z

I need to configure the silencing of specific lints because my OCD won't allow me to leave a warning alive :_

alecandido · 2023-07-07T18:24:45Z

Isn't this one?
https://github.com/rust-lang/rust-clippy#allowingdenying-lints

e.g. #![deny(clippy::single_match, clippy::box_vec)]

I guess the proper one to be: #![deny(clippy::suboptimal_flops)]

alecandido · 2023-07-07T18:30:29Z

Moreover, it seems you could also do it per-function:

allow/warn/deny can be limited to a single function or module using #[allow(...)], etc.

https://github.com/rust-lang/rust-clippy#allowingdenying-lints

I guess it should work as an attribute:

#[deny(clippy::suboptimal_flops)]
fn interpolate(&self, query: f64) -> Result<f64, InterpolationError> {
    ...
}

But I haven't tried myself, so I'm not sure.

scarlehoff · 2023-07-07T18:35:13Z

I'll play around with it next week.

For the rest, if you are happy with the "generalization" in this branch, next week* I'll fix/silence clippy warnings and implement the 2D version.

*I intend on doing this on the train on Thursday more specifically. I might even pass by the department since I'll have a 3 or 4 hours train-change in Milan :)

ndinterp/src/grid.rs

alecandido · 2023-07-08T11:31:47Z

ndinterp/src/grid.rs

+    pub fn closest_below<const N: usize>(
+        &self,
+        input_query: &[f64],
+    ) -> Result<[usize; N], InterpolationError> {


Here, I believe you'd like to use something like D::NDIM.

Unfortunately, D::NDIM: Option<usize>, so you'd need something like .unwrap(), but executed at compile time.

https://docs.rs/ndarray/0.15.6/ndarray/trait.Dimension.html#associatedconstant.NDIM

A simpler alternative might be to split the bulk of the function, and then implement thin wrappers around (one per Dimension).
Even the implementors of Dimension itself are just 7 (+ 1).

Otherwise, you should be able to use a const function to unwrap it, defining an associated constant, and try to use the associated constant in the signature.

https://users.rust-lang.org/t/compile-time-const-unwrapping/51619

EDIT: Unfortunately, const panic is not stabilized, because of an issue with Rust 2021 rust-lang/rust#85194, maybe you could check if const pattern matching is allowed

A const function to unwrap IxN at compile time you mean?

Indeed, in order to turn it into a usize (if you can not panic, you can always turn IxDyn into a 0, it's useless anyhow, and it should not happen, so I wouldn't care that much).

I think this cannot be done. Or at least I cannot find a way of passing the D to the const function:

generic parameters may not be used in const operations type parameter may not be used in const expressions

and in any case, if I use a function that could depend on one of the parameters:

const parameters may only be used as standalone arguments

Not sure whether that means there's no way of doing what I want to do (without creating a path per dimension, at that point I rather leave the N)

I'd really like to solve this issue. But if it's not simple, we could postpone.

I'd like to as well, but I haven't been able to...

ndinterp/src/grid/cubic.rs

ndinterp_capi/src/lib.rs

scarlehoff · 2023-07-11T12:46:02Z

I've decided to remove pedantic because even the clippy docs talk about false positives I don't want to sprinkle allows around the entire code. For some reason that I cannot fully understand the warnings for ndinterp_capi/Cargo.toml are not fixed even all the metadata is there...

ndinterp_capi/Cargo.toml

alecandido · 2023-07-11T15:26:57Z

@scarlehoff whenever you're ready, feel free to merge

scarlehoff · 2023-07-16T09:45:07Z

@alecandido have a look when you have time to the last commit

The motivation for the GridSlice is that at various stages of other-D interpolation one needs to perform interpolation in D-1 version of the grid. Since for now we are doing 1 and 2 I've decided GridSlice is going to be limited to 1D (but it shouldn't be to painful to extend).
As a very nice side-effect, the explicit access to everything as a view allows for some optimization by the rust compiler and now is even faster than before.

After this I'll start with the 2D interpolation, I promise :P

alecandido · 2023-07-16T11:58:16Z

ndinterp/src/grid.rs

-impl Grid<Ix1> {
+/// A grid slice is always 1-Dimensional
+#[derive(Debug)]
+pub struct GridSlice<'a> {


Eventually, you started playing with lifetimes...

I don't understand why the compiler is not able to deduce all of them though.

There are rules to infer references lifetimes for methods, but even those are limited (to simplest but also most common cases, so they are very useful). You definitely want to infer as much as possible, but consistency is hard to impose when you have many rules (I recently read about attempts to introduce linear types, and they are a mess...).

I'm not sure, but I believe the reason you need to make them explicit should be that there is an alternative way of doing, e.g. having more lifetimes and some of them explicitly specified (like generics), or something like that. Maybe?

alecandido · 2023-07-16T12:00:37Z

As a very nice side-effect, the explicit access to everything as a view allows for some optimization by the rust compiler and now is even faster than before.

Maybe the optimization boils down to fewer allocations. Could it be?

scarlehoff · 2023-07-16T12:34:14Z

If you are happy with the GridSlice (you can see the use-case in #5) I'd merge this as it is.

Maybe the optimization boils down to fewer allocations. Could it be?

If it is due to that I cannot see where the extra allocations happen.

alecandido · 2023-07-16T12:51:35Z

If you are happy with the GridSlice (you can see the use-case in #5) I'd merge this as it is.

I'm going to rereview quickly, but I guess we can merge.

If it is due to that I cannot see where the extra allocations happen.

I assumed you were using GridSlice in place of Grid, so you were avoiding some local copies...

scarlehoff · 2023-07-16T15:56:41Z

I assumed you were using GridSlice in place of Grid

For the 1D case this is (almost) literally the case, but the only difference is that now the compiler has a bit more information: "From here onwards nothing will be touched". For the 2D case GridSlice is a 1D slice of Grid.

alecandido · 2023-07-16T21:24:00Z

ndinterp/src/grid.rs

-    pub values: Array1<f64>,
+pub struct Grid<D: Dimension> {
+    /// Arrays with the input vectors (x_i)
+    pub xgrid: Vec<Vec<f64>>,


Here an even nicer option (for later on) would be something like Vec<[f64; N]>, with N: usize determined from D.

But it is again the same problem we already discussed before...

alecandido · 2023-07-16T21:25:36Z

ndinterp/src/grid.rs

+#[derive(Debug)]
+pub struct GridSlice<'a> {
+    /// A reference to one of the input vectors of the grid
+    pub x: &'a Vec<f64>,


Are you sure you need a reference to a vector, and not a slice?

(with slice I mean &'a [f64], as per standard library)

Depends on what you mean with “need”, that might also work (and be more correct?), but x is a reference to one of the quantities in the xgrid which are vectors hence the choice.

If you need to use vector-specific operations, then yes, you need the reference to a vector.
But if you need to reference just a collection of points, the slice is more general, and it relies less on the underlying structure (while still be plenty of possible methods, https://doc.rust-lang.org/std/primitive.slice.html).

However, if it has to be one entire element (the full x coordinates) maybe it makes sense to keep as it is. I'd suggest instead to give it a name:

type XPoint = Vec<f64>;

(if you're able to find a better name, you're welcome). It's just an alias, for greater clarity, but here a newtype would be too much (I believe).

Ah, no. I do need a vector because of the search function.

But it is not a point, it is actually the x grid for \vec{y} = f(\vec{x}). I could call them xgrid, I decided to change it to x since it is not exposed to the user and only used internally (as opposed to the xgrid which needs to be filled by the user).

But I can use xgrid as well here if you prefer (but I'll do the change directly in #5)

"better name" was referring to the name of the type alias I was proposing. x for the attribute is just fine.

The reason behind the alias is that an XPoint is isomorphic to a Vec<f64>, but not all Vec<f64> are XPoints. It doesn't help with automatic check, because an alias is interchangeable with the aliased type, but it helps the developer reading it (to have help also from the compiler, we could use a newtype, but as I said it seems to me that it's not worth, at the moment)

alecandido · 2023-07-16T22:45:51Z

ndinterp/src/grid.rs

-    pub fn derivative_at(&self, index: usize) -> f64 {
-        let dx = self.input[index] - self.input[index - 1];
-        let dy = self.values[index] - self.values[index - 1];
+    pub fn derivative_at(&'a self, index: usize) -> f64 {


The annoying part of lifetimes is that you then need to specify them everywhere. Luckily not in calls, but in all types yes...

alecandido · 2023-07-16T22:49:11Z

ndinterp/src/grid/cubic.rs

+impl<'a> GridSlice<'a> {
+    /// Implements utilities for a GridSlice that can be used by cubic interpolation Nd
+    /// Takes as input the value being queried and its index within the given slice
+    fn cubic_interpolate_1d(&'a self, query: f64, idx: usize) -> Result<f64, InterpolationError> {


Does it make sense to call it 1d? A slice is always 1D... (or are you already planning to extend beyond?)

Who knows… might be necessary in the future.

scarlehoff · 2023-07-20T13:41:55Z

I'll merge this (so that further discussion can happen in #5 or other future changes)

scarlehoff added 2 commits July 6, 2023 09:54

update code to prepare for n-d

aac5264

fix double free

bd1a08d

scarlehoff requested review from cschwan and alecandido July 6, 2023 08:48

add comments and tests

bf93fb3

scarlehoff commented Jul 6, 2023

View reviewed changes

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

cschwan reviewed Jul 6, 2023

View reviewed changes

ndinterp_capi/src/lib.rs Outdated Show resolved Hide resolved

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

make closest below be dimension independent

4ac6fda

alecandido reviewed Jul 8, 2023

View reviewed changes

scarlehoff added 3 commits July 11, 2023 13:51

put clippy warnings to defaults

0513987

fix warnings

c2633b9

apply review comments

84dbc07

alecandido reviewed Jul 11, 2023

View reviewed changes

ndinterp_capi/Cargo.toml Outdated Show resolved Hide resolved

scarlehoff added 2 commits July 12, 2023 09:42

inherit from main cargo.toml

7088c68

create a gridslice struct

2c0dae3

alecandido reviewed Jul 16, 2023

View reviewed changes

scarlehoff merged commit 381fa8f into main Jul 20, 2023

scarlehoff deleted the prepare_for_ndinterpolation branch July 20, 2023 13:42

Update the code for n-d #4

Update the code for n-d #4

Conversation

scarlehoff commented Jul 6, 2023

scarlehoff commented Jul 6, 2023

scarlehoff commented Jul 6, 2023 • edited Loading

alecandido commented Jul 6, 2023

scarlehoff commented Jul 6, 2023 • edited Loading

cschwan commented Jul 6, 2023 • edited Loading

scarlehoff commented Jul 6, 2023 • edited Loading

cschwan commented Jul 6, 2023

scarlehoff commented Jul 6, 2023

alecandido commented Jul 7, 2023

scarlehoff commented Jul 7, 2023

alecandido commented Jul 7, 2023

scarlehoff commented Jul 7, 2023

alecandido commented Jul 7, 2023 • edited Loading

alecandido commented Jul 7, 2023

scarlehoff commented Jul 7, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alecandido Jul 8, 2023 • edited Loading

Choose a reason for hiding this comment

scarlehoff Jul 11, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

scarlehoff commented Jul 11, 2023

alecandido commented Jul 11, 2023

scarlehoff commented Jul 16, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alecandido commented Jul 16, 2023 • edited Loading

scarlehoff commented Jul 16, 2023 • edited Loading

alecandido commented Jul 16, 2023

scarlehoff commented Jul 16, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

scarlehoff commented Jul 20, 2023

scarlehoff commented Jul 6, 2023 •

edited

Loading

scarlehoff commented Jul 6, 2023 •

edited

Loading

cschwan commented Jul 6, 2023 •

edited

Loading

scarlehoff commented Jul 6, 2023 •

edited

Loading

alecandido commented Jul 7, 2023 •

edited

Loading

alecandido Jul 8, 2023 •

edited

Loading

scarlehoff Jul 11, 2023 •

edited

Loading

alecandido commented Jul 16, 2023 •

edited

Loading

scarlehoff commented Jul 16, 2023 •

edited

Loading