search

moble / scri

Python/numba code for manipulating time-dependent functions of spin-weighted spherical harmonics on future null infinity
MIT License
18 stars 20 forks source link

Enhancement: Allow reuse/precomputation of interpolating splines in ModesTimeSeries (and abd) #73

Open duetosymmetry opened 2 years ago

duetosymmetry commented 2 years ago

Currently, calling interpolate, derivative, or antiderivative on an instance of ModesTimeSeries will construct a scipy.interpolate.CubicSpline interpolant, and then throw it out when the function returns. Building the interpolant is, AFAIK, most of the time in performing the interpolation/derivative/antiderivative. This means that calculations can't benefit from a speedup by reusing the interpolant.

I see two approaches to reusing the spline interpolant.

One, make it the user's problem; add an optional keyword argument spline to those three routines, and the user is responsible for passing in the correct interpolant (if they pass in garbage, they get back garbage). Along with this change also add a function to compute and return the interpolant.

Two, try to handle saving/reusing the interpolant inside the object. This would mean adding a member variable to store the interpolant, computing it the first time it's requested, and reusing it whenever needed.

The advantage of the second approach is that it will speedup any code that currently could gain a speedup by reuse.

The disadvantage of the second approach is that the spline may become invalid if somebody modifies the data of the object, and I don't know if it's possible to (cheaply) determine if the spline has become invalid.

Comments? @moble @keefemitman

keefemitman commented 2 years ago

I think the 2nd option would probably run into issues more often than we expect. For example, people often will just change the time member variable after reading in data as an MTS. But this would require the spline to be recreated, which is an extra computation that really shouldn't be present in such a simple task.

Because this is really only relevant in speeding up frame transformations (I can't think of other instances where we'd want to use this feature?) I think it's probably best if we just make an optional keyword argument spline for these routines.

But I'm not necessarily set on this and will defer to Mike.

moble commented 2 years ago

I actually implemented the second version at some point — maybe in sxs, but the principle is the same. Basically, I made a function that copied the original object (and its data), attached the precomputed spline, and changed the various interpolation, dot, antiderivative methods to use that spline. (Some times, I'll admit that python's flexibility is nice.)

As for modifying the data, that's not a problem because numpy lets you set the writeable flags to False. So this new object's t, data, and frame couldn't be modified; any code that tried would raise an error. That restriction also isn't a problem, because the functions for which this is useful (e.g., BMS transformation or doing the algebra in all these charge calculations) don't modify waveforms in place.

The downside of this approach is, of course, the extra memory it requires. Incidentally, the slowest part of building a spline is making this weird tridiagonal matrix from the time data, and then taking its LU decomposition. The interpolated data are basically put into a linear solver with that LU decomposition, which is relatively fast. So in my Julia code, I just store that LU decomposition, and for each mode I redo the linear solve, interpolate, and move onto the next mode. That uses a lot less memory, and is just as fast as storing the whole spline for all the modes.