This is somewhat tricky as Lancszos interpolation used in the pipeline is buried quite deep and within C layer of the code. On the other hand, I have not been able to find  readily available  Lancszos interpolation codes in python that would reliably reproduce the results of the pipeline. I have found an implementation in torch which I could strip and modify, but it is not obvious it is worth the time and effort given the relatively small difference between currently used bilinear interpolation in my code and Lancszos code - and the fact that I will probably need to refactor my code to be LSST code based at some point in the future. I have to consider this a little bit further. 

It should be trivial to implement exactly the same Lanczos kernel in python (although I'm surprised that it's not in scipy).  I can do that if needs be.

Indeed, I have been able to implement it on my own (see small example of shift with maaaany wiggles). I searched quite extensively and was not able to identify Lanszscos (capable to going to 5) implementation that was readily usable (for example, no mention in scipy https://docs.scipy.org/doc/scipy/reference/search.html?q=lanczos&check_keywords=yes&area=default , only Lanczos fitting algorithm / one could construct it with scipy.special.sinc). I am just expressing some worry rewriting everything with my implementation which might be prone to errors and the amount of time that might take given other priorities.

After discussion with Robert on Slack, I have made slight modifications to the algorithm so that it can fully ran under LSST imports. There were some minor changes, which probably arise from slight differences in python version and other imports, but which do not change the final image. The only larger difference is that I was using scikit-image libarary for one resizing operation, which was trivially replaced.
I have not yet tested if I could run the fitting algorithm also under LSST code (the algorithm which searches for best parameters of the fit), which may have other dependencies not in LSST at the moment (most likely there might be some problems in the way how I do multiprocesing?). I would have to check that.

For the same test image in the focus, running on tiger-sumire, I am getting:

without LSST stack: `840 ms ± 3.69 ms per loop(mean ± std. dev. of 7 runs, 1 loop each)`
with LSST stack: `1.29 s ± 106 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)`

I ran %timeit many times and the result are roughly the same each time. I am confused my much lower variance in `without LSST' result. This result suggests that the slow down is around 50% for in-focus images. I have to investigate where the slowdown is coming from (i.e., is it from FFT calculation and using/not using of MKL libraries)

After some consideration, and talking with Josh Meyers, and looking how the pipeline code is set up, I think it is easiest and most convenient to use galsim.InterpolatedImage functionality. This makes it possible to use Lanczos5 as an interpolant and plays nicely with other functionality that I have + the pipeline as well. I have tested that the shift of 0.5 pixel with the pipeline function on an oversampled image in focus (with fuction lsst.afw.math import offsetImage, which is invoked in the pipeline) creates exactly the same change as when using galsim. This is shown in comparison_galsim_pipeline.png

That sounds fine, but won't you still have to integrate with the "real" Rubin pipeline soon?

Yes, indeed. I think we need to talk so that I can understand how this will be achieved, i.e., how will the psf object be constructed. But note, as you are aware, that galsim is LSST package. And I do need some functionality to do fitting of defocused and focused images which I dont think it would be easily run through ``real'' pipeline - but in this way it will be using the same algorithms as the ``real'' pipeline.

I have incorporated galsim Lanczos5 procedures at all relevant places in my module, internal version 0.47a. Closing this ticket, and evaluating in the next release of the solutions.