Massive performance boost possible in FittableImageModel

[krachyon]

Description

I'm currently messing around with the limitations of fitting larger groups of sources where the fit gets painfully slow. I noticed that 80% percent of the time spent was in the evaluate() method of the EPSF-model, digging into this more, the call to the interpolator seemed really slow when compared to hand-written linear interpolation.

It turns out it could be massively sped up by using a different way of calling RectBivariateSpline, ev() seems to be made for querying unstructured points, but with the grid=True argument in __call__() a regular grid of points can be interpolated much more effectively, see snippet below.

What to do about it?

I'm not sure how to properly make use of this though, using a sparse grid as input is a bit at odds with the interface provided by the modelling framework. The fitter expects arrays of same size for x,y and function value and one could possibly want to use the model as a plain function-object so I'm not sure about the ramifications for other usecases.
Also the logic to select the used pixels in BasicPsfPhotometry.do_photometry() would have to change and I'm not sure what the most elegant solution or possible unpleasant side effects would be, hence why this isn't a PR.

I'm starting to think that using the model abstraction is more of a roadblock for this usecase as it also requires the whole parameter-name matching machinery.
I'm in the process of create a variant of the PsfPhotometry class that uses the optimizer more directly to set bounds on parameters and to avoid the penalty of a large call-stack through the compound-model evaluation layers (see also #1217) and think this is actually a bit easier to implement and reason about. I'll link it here or create a PR when I have something presentable.

Demo

from scipy.interpolate import RectBivariateSpline
import numpy as np
import timeit

rng = np.random.default_rng(seed=10)


def benchmark(data_size, query_size):
    dummy_data = rng.uniform(0, 1, (data_size, data_size))

    y_input, x_input = np.ogrid[-300:300:data_size*1j, -300:300:data_size*1j]

    interpolator = RectBivariateSpline(x_input, y_input, dummy_data, kx=3, ky=3)

    y_grid, x_grid = np.mgrid[-200:200:query_size*1j, -200:200:query_size*1j]
    y_sparse, x_sparse = np.ogrid[-200:200:query_size*1j, -200:200:query_size*1j]

    def use_ev():
        return interpolator.ev(x_grid, y_grid)

    def use_call():
        return interpolator(x_sparse, y_sparse, grid=True)

    N = 5
    rep = 5
    ev = np.array(timeit.Timer(use_ev).repeat(repeat=rep, number=N))/N
    call = np.array(timeit.Timer(use_call).repeat(repeat=rep, number=N))/N
    print(f'{data_size=}, {query_size=}')
    print(f'using ev: {np.mean(ev)} sec/iter; std: {np.std(ev)}')
    print(f'using call: {np.mean(call)} sec/iter; std: {np.std(call)}')
    print(f'ratio: {np.mean(ev)/np.mean(call)}')
    print('\n')

benchmark(501, 700)
benchmark(201, 900)
benchmark(1024, 512)

output:

data_size=501, query_size=700
using ev: 0.3547264989600262 sec/iter; std: 0.018836194139951606
using call: 0.01233470320001288 sec/iter; std: 5.7584694313710254e-05
ratio: 28.75841381896046

data_size=201, query_size=900
using ev: 0.3189859828400222 sec/iter; std: 0.005763221323098215
using call: 0.020830674280023234 sec/iter; std: 0.0006746705740365767
ratio: 15.313281680273406

data_size=1024, query_size=512
using ev: 0.33410044439997366 sec/iter; std: 0.00031986051640275724
using call: 0.007152011239995773 sec/iter; std: 2.2031428508788802e-05
ratio: 46.71419453755936