Enhanced interpolation improvement

reichherzerp commented 3 years ago

Added comments and implemented consisting naming for functions and variables. Improved the readibility of some parts of the code.

antoniusf commented 3 years ago

Hey, I just pushed two commits: – fixing the bug where ipolType wasn't set – two of the constructors were duplicated somehow (maybe from the merge?) and the code didn't build for me. I removed the duplicates.

antoniusf commented 3 years ago

So I just checked the new SIMD version against the old version as well, and they are mostly equal (to about 1e-7, which looks like it's float precision, but isn't, since it's an absolute error. The maximum relative error is actually about 1e-3). I am a little bit confused by this, since I'm not sure why these errors were introduced (compiler optimizations maybeeee?), but otherwise the fields are pretty much equal and that's a good thing I think.

antoniusf commented 3 years ago

Okay so I don't think this is a problem, I just re-arranged the SIMD computation in exactly the same order I would expect the original formula to be evaluated, and the error went to 0. I think this means that we can safely assume that the above error is solely due to errors from inaccuracies in the floating point format. (As you know, floating point operations like addition and multiplication are not associative, strictly speaking. While the results are generally close to each other (and to the actually correct result), they are not guaranteed to be the exact same floating point number. Each step introduces a tiny little rounding error, which can accumulate and be amplified, thus leading to the differences I observed.)

antoniusf commented 3 years ago

I'm not sure if this is the right place for this, given that we mainly want to review Patrick's changes, but I don't have another place to put this right now so it's going to go in here.

Right now, the tricubicInterpolate functions use periodicGet to retrieve grid points while fulfilling boundary conditions – except that periodicGet only works for periodic boundary conditions (obviously). I think we just never got around to adding the same thing for reflective boundaries, but we should definitely do that before it gets merged into CRPropa master.
There are also no generic implementations for tricubicInterpolate. (There's one for double and one for Vector3f, but what about float or Vector3d?) I think we could probably just take the double version and modify it to work generically.

Again, both of these points are not things that have to be done in this PR (especially since they are not related to anything Patrick did but rather issues that have been present for a long time), but I wanted somewhere to write them down, and they came up while I reviewed the code.

reichherzerp commented 3 years ago

Okay so I don't think this is a problem, I just re-arranged the SIMD computation in exactly the same order I would expect the original formula to be evaluated, and the error went to 0. I think this means that we can safely assume that the above error is solely due to errors from inaccuracies in the floating point format. (As you know, floating point operations like addition and multiplication are not associative, strictly speaking. While the results are generally close to each other (and to the actually correct result), they are not guaranteed to be the exact same floating point number. Each step introduces a tiny little rounding error, which can accumulate and be amplified, thus leading to the differences I observed.)

Thanks for checking, this makes total sense! So you think that we can just let it as it is, or do you prefer your ordering? I don't mind.

antoniusf commented 3 years ago

I don't think there's a good reason to use my ordering, it's way less readable than what's in there right now. It might be worth considering whether it would be valuable to change the scalar implementation instead to have both implemenations (with SIMD and without) return the exact same values, or at least have them match within 1e-7. But for what it's worth, the TD13 optimization doesn't provide this either and in fact the fields TD13 produces are only comparable for small enough x.

reichherzerp commented 3 years ago

Thanks Antonius for your great review and sorry for taking so long to respond!

antoniusf commented 3 years ago

Thanks for taking the time to fix things!

LeanderSchlegel / CRPropa3

Enhanced interpolation improvement #1