Cubic interpolation 1d

[scarlehoff]

First implementation of the cubic_1d interpolation algorithm that LHAPDF uses for alpha_s. At the moment I put everything in src/grid.rs because I am an inhuman beast, but also because I was not sure about the organization.

Maybe having the Grid struct in src/grid.rs and then inside src/grid a file for every interpolation method that can be implemented for the struct Grid? That way the user (i.e., partons) can decide whether to call cubic_1d or magic_interpolation_1d depending on the necessities of the interpolation.

In any case... there are 4 pieces of code:

• pub struct Grid ...
• impl Grid { pub fn derivative_at ...
• pub trait Interpolation { fn cubic_1d...
• impl Interpolation for Grid ... { fn cubic_1d..

What is the rustonic (rustacean?) way of organising this? I'm guessing the first two would go together and then in a different file (maybe inside src/grid/interpolation.rs) the trait Interpolation is defined? But where would the specific functions implemented for the trait? Would interpolation.rs act as a header file? (I'm poisoned by C/C++, please liberate me).

I've also added an alphas example which just runs over a few values of alphas(q) and compares them to the LHAPDF result.

Still to do:

From the point of view of the functionalities

• [x] What to do at the borders (right now it will just crash when taking the point at i+1
• [x] (decide) what to do in extrapolation (raise an exception? Should I be returning a Result? say no? (extrapolate? Should we "extrapol" as well as "interpol"?) Maybe we implement some different options and then partons can decide whether to do extrapolate_error or extrapolate_constant or whatever.

From the point of view of the code

• [ ] Clean it up. I'm sure there are better ways to find the index...
• [ ] Vectorize it? The 1d case is simple enough that it might be worth it -for me- learning how to do it here before even starting with 2d... This greatly depends on me understanding well enough the iter->enumerate->map

[alecandido]

Maybe having the Grid struct in src/grid.rs and then inside src/grid a file for every interpolation method that can be implemented for the struct Grid?

Consider that src/grid.rs is the equivalent of src/grid/__init__.rs. That said it might be fine (in principle in Rust, you don't even have the problem of circular imports, since it's not happening at runtime).

[alecandido]

What is the rustonic (rustacean?)

Someone even said "Rustic, I'm a native speaker", but I've seen someone else replying "Then you suck at your own language" (since rustic is an actual word, with a different meaning).

My answer to this dilemma is just idiomatic in Rust. And if it's clear from the context, idiomatic is already sufficient. Though boring...

P.S.: a Rustacean is a developer writing in Rust, so essentially you; the word comes after Ferris the Crab https://rustacean.net/

[alecandido]

What is the rustonic (rustacean?) way of organising this? I'm guessing the first two would go together and then in a different file (maybe inside src/grid/interpolation.rs) the trait Interpolation is defined? But where would the specific functions implemented for the trait? Would interpolation.rs act as a header file? (I'm poisoned by C/C++, please liberate me).

I believe the correct way of thinking about this in Rust is just visibility, not the definition order (like in C/C++, but also in Python):

This works:

struct Come {}

impl Ciao for Come {
    fn ciao(self) -> String {
        "come va?".to_owned()
    }
}

trait Ciao {
    fn ciao(self) -> String;
}

fn main() {
    println!("{}", (Come {}).ciao());
}

For as long as the object you want to reference is in a visible scope, it doesn't matter if it's before or after, the compiler will work it out for you :)

So, feel free to implement the functions in the files inside src/grid/.

However, consider that every trait can be implemented only once per struct, and the Interpolation trait would work better if it's unique (even though generic).

To me, the best option is to implement one struct for interpolation method you want to define. If Grid is common to all of them, include it as a member. If it is the only member, adopt the newtype pattern https://doc.rust-lang.org/rust-by-example/generics/new_types.html

(if instead you might have something like "kinds" of grids, we could consider collecting the related interpolators in enums, but this we can worry later on, when we will have a bunch of these interpolation functions)

[alecandido]

From the point of view of the functionalities

Return Result errors every time you are not able to do what you're expected to do. You can define specific type of errors in a simple way with the thiserror crate, to allow specialized handling down below.

However, this is an interpolation library, extrapolation is beyond the scope, it will be up to the library user :+1: (we can change perspective at some point, but I'd try to keep a minimal scope for the time being) (and in particular I would just raise an error in partons as the default, and for some time there won't be much more than the default)

[scarlehoff]

This is ready for (first) review, as it is now functionally complete (but with decisions to be taken/discussed).

At the moment everything happening in interpolation follows what LHAPDF does. As per the above discussion, I think that this is the best option until we can reproduce all LHAPDF-interpolation behaviour. Once that's done we can improve on it/add new features.

In extrapolation I decided I want to try that thiserror library that I found while looking through pineappl and have created an Error for the extrapolation above and below the max/min. The idea is that when that happens this error is returned and it is the job of the caller to decide what to do there (this is called ndinterp and not ndextrap after all). We might want to also add "extrapolation utilities" at some point (and we will need them either here or in partons if we want to be lhapdf-feature-complete) but for the time being this is enough.

When reviewing the code please complain to your heart's content. There's a bit of cargo-cult-coding in the error handling btw, might be correct but I don't think I fully comprehend what I'm doing.

[alecandido]

The idea is that when that happens this error is returned and it is the job of the caller to decide what to do there (this is called ndinterp and not ndextrap after all).

This is sensible. In principle, the caller can even avoid getting the error, but if he doesn't, then an error will be.

We might want to also add "extrapolation utilities" at some point (and we will need them either here or in partons if we want to be lhapdf-feature-complete) but for the time being this is enough.

This is what LHAPDF is doing, since (in principle) you could select the interpolator and extrapolator (e.g. the "error extrapolator" should be available, or even simple saturation).

However, to improve over LHAPDF choices, if ever, there will also be an extrapolation crate. But (most likely) I'm not the one who's going to write. From a physics perspective, extrapolation is black magic at best (otherwise just random numbers). More a hope than an ansatz.

There's a bit of cargo-cult-coding in the error handling btw, might be correct but I don't think I fully comprehend what I'm doing.

In this context, cargo is ambiguous :P

[scarlehoff]

So, what should the next steps be?

I think it would be a good idea to try and use this already in partons as it is developed in order to find possible pitfalls as early as possible.

Personally I'd like to have a go at vectorization but, depending on the stage at which partons is, it might make more sense to go from here to 2d interpolation in order to be able to use PDF grids asap.

[alecandido]

So, what should the next steps be?

The bicubic, definitely (and I'll finish reviewing this today, but I suspect is just missing a green button click).

I think it would be a good idea to try and use this already in partons as it is developed in order to find possible pitfalls as early as possible.

Yes and no: Partons is already supporting for grids, because they are a more fundamental feature than $\alpha_s$, but missing the support for this one. Since tests are working and sensible, I would accept this, but not sure if implementing $\alpha_s$ is really top-priority for Partons (there are other clean-ups ongoing, and limited manpower available).

Personally I'd like to have a go at vectorization but, depending on the stage at which partons is, it might make more sense to go from here to 2d interpolation in order to be able to use PDF grids asap.

Indeed, 2d then vectorization would be appreciated :P

[scarlehoff]

Not that alpha_s should be top priority but rather having things being used will help discover possible issues (and which are probably relevant also for the pdf)

[alecandido]

Not that alpha_s should be top priority but rather having things being used will help discover possible issues (and which are probably relevant also for the pdf)

I know the feeling, and I agree to the spirit. It's on my to-do list.

P.S.: you could also try to author it yourself; I have a specific idea on how to do it, but if you have a first draft, I can also work on that (just to say that if you have time, don't feel blocked by myself, you should have full write access to the repo :+1:)

[scarlehoff]

just to say that if you have time, don't feel blocked by myself, you should have full write access to the repo

I think I do, but since this is not time sensitive in any way I prefer if you do it. Once I feel comfortable enough with rust that I think I can find my own mistakes I will also start touching partons.

I'm treating this as a toy I can break and rebuild (I know I know, PRs and version control, but that takes time from everyone involved)

[alecandido]

@scarlehoff for me you can merge this whenever you wish :)

[scarlehoff]

I want to add a benchmark with the LHAPDF interpolation in the example before merging this and starting with the PDF interpol. Then it'd be done.

[alecandido]

A benchmark in the sense of performances?

If that's the case, remember that Cargo has a dedicated command for benchmarks, cargo bench, and a default directory benches/ (at the same level of src/). Moreover, there are dedicated tools for writing benchmarks, like Criterion. I would suggest taking a look at this page: https://nnethercote.github.io/perf-book/benchmarking.html

In case you want an example of how to use Criterion, we already have a pretty basic one in quad: https://github.com/AleCandido/quad/blob/main/quad/benches/my_benchmark.rs

[scarlehoff]

I'll merge this since it makes more sense to separate the benchmarks into another pull request.