Backend for ducc0.nufft?

[mreineck]

Hi,

do you think it might be worthwhile to provide a wrapper for ducc0.nufft (https://gitlab.mpcdf.mpg.de/mtr/ducc)?
This is currently packaged for Python only, but the backend is C++ and consists only of two global functions for carrying out type 1 and 2 non-uniform FFTs, so I think making this accessible from Julia should not be too difficult.

This implementation

• does not support explicit planning, making the interface much simpler
• is as fast or faster as a pre-planned FINUFFT transform nevertheless, and
• for the same requested accuracy, provides better actual accuracy than FINUFFT in almost all cases

Unfortunately I don't have any exprience with Julia myself (yet), so I was not able to create direct benchmarks against NFFT.jl.
If you are interested, I have created a small benchmark script, which tries to imtate the calculations you are showing in the "performance" section of your docs.

import ducc0
import numpy as np
from time import time


class Bench:
    def __init__(self, shape, npoints):
        self._shape = shape
        ndim = len(shape)
        self._coord = (2*np.pi*np.random.uniform(size=(npoints,ndim)) - np.pi).astype(np.float32)
        self._points = (np.random.uniform(size=npoints)-0.5
                 + 1j * np.random.uniform(size=npoints)-0.5j).astype(np.complex64)
        self._values = (np.random.uniform(size=shape)-0.5
                 + 1j * np.random.uniform(size=shape)-0.5j).astype(np.complex64)

        self._res_fiducial_1 = ducc0.nufft.nu2u(
            points=self._points.astype(np.complex128),
            coord=self._coord.astype(np.float64),
            forward=True,
            epsilon=3e-13,
            nthreads=0,
            verbosity=1,
            out=np.empty(shape, dtype=np.complex128))

        self._res_fiducial_2 = ducc0.nufft.u2nu(
            grid=self._values.astype(np.complex128),
            coord=self._coord.astype(np.float64),
            forward=True,
            epsilon=3e-13,
            nthreads=0,
            verbosity=1)

    def run(self, epsilon, singleprec, nthreads):
        rdtype = np.float32 if singleprec else np.float64
        dtype = np.complex64 if singleprec else np.complex128

        res={}
        res["shape"] = self._shape
        res["npoints"] = self._coord.shape[0]
        res["epsilon"] = epsilon
        res["nthreads"] = nthreads
        res["singleprec"] = singleprec

        shape = self._shape
        ndim = len(shape)
        npoints = self._coord.shape[0]
        coord = self._coord.astype(rdtype)
        points = self._points.astype(dtype)
        values = self._values.astype(dtype)

        out = np.ones(shape, dtype=dtype)
        t0 = time()
        res_ducc = ducc0.nufft.nu2u(points=points, coord=coord, forward=True, epsilon=epsilon, nthreads=nthreads, verbosity=1, out=out)
        res["ducc_full_1"] = time()-t0
        res["err_ducc_1"] = ducc0.misc.l2error(res_ducc, self._res_fiducial_1)

        out=np.ones(shape=(npoints,), dtype=dtype)
        t0 = time()
        res_ducc = ducc0.nufft.u2nu(grid=values, coord=coord, forward=True, epsilon=epsilon, nthreads=nthreads, verbosity=1)
        res["ducc_full_2"] = time()-t0
        res["err_ducc_2"] = ducc0.misc.l2error(res_ducc, self._res_fiducial_2)
        return res


def plot(res, fname):
    import matplotlib.pyplot as plt
    tducc1 = [r["ducc_full_1"] for r in res]
    tducc2 = [r["ducc_full_2"] for r in res]
    educc1 = np.array([r["err_ducc_1"] for r in res])
    educc2 = np.array([r["err_ducc_2"] for r in res])
    eps = np.array([r["epsilon"] for r in res])
    plt.xscale("log")
    plt.yscale("log")
    plt.plot(educc1,tducc1,label="ducc, type 1")
    plt.plot(educc2,tducc2,label="ducc, type 2")
    plt.title("shape={}, npoints={}, nthreads={}".format(res[0]["shape"], res[0]["npoints"], res[0]["nthreads"]))
    plt.xlabel("real error")
    plt.ylabel("t[s]")
    plt.legend()
    plt.savefig(fname)
    plt.close()


def runbench(shape, npoints, nthreads, fname, singleprec=False):
    res=[]
    mybench = Bench(shape, npoints)
    epslist = [4e-6, 1e-5, 1e-4, 1e-3, 1e-2] if singleprec else [4e-13, 1e-12, 1e-11, 1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]
    for eps in epslist:
        res.append(mybench.run(eps, singleprec, nthreads))
    plot(res, fname)


# NFFT.jl benchmarks
runbench((   512*512,),  512*512, 1, "bench_1d.png", False)
runbench((   512,512,),  512*512, 1, "bench_2d.png", False)
runbench((   64,64,64,), 64*64*64, 1, "bench_3d.png", False)

[tknopp]

Hi Martin,

yes in general having more more wrapped libraries is a good thing. If this is just for comparison with NFFT.jl then we can host that within NFFT.jl itself, see https://github.com/JuliaMath/NFFT.jl/tree/master/Wrappers
That would probably the best starting point anyway. But these wrappers are not meant to be "official" since I don't want this repository to support/maintain all NFFT wrappers.

If you want to provide real Julia bindings to ducc0 I would recommend that you create a dedicated Julia package and implement the AbstractNFFT interface. In that way, other users can exchange libraries.

Having said that, I can support you creating the Wrapper. The most important thing (in which I unfortunately don't have experience) is to create a "Binary Julia" package like this one: https://github.com/JuliaBinaryWrappers/finufft_jll.jl. You can see the build script of finufft here: https://github.com/JuliaPackaging/Yggdrasil/blob/master/F/finufft/build_tarballs.jl Yggdrasil is the repository that holds the build scripts and the binaries itself are then created using CI services.

Regarding C++, the best would be to provide a C Wrapper, since we can call that directly without any problems. When you are at that point, I can help you.

Independent of this, I would be interested in what you do different than other NFFT implementations:

• Is it the classical algorithm or is it a completely different approach?
• Do you use the same blocking technique like finufft?
• Do you use FFTW and if not is your FFT faster?
  In case that ducc0.nufft is faster than NFFT.jl it would be great if you could explain some tricks that you applied. We just submitted a paper that describes what NFFT.jl does: https://arxiv.org/abs/2208.00049

Side comment:
"does not support explicit planning, making the interface much simpler"
I never get why this makes an interface simpler. Grouping the memory in a struct plus using caching techniques is an established programming pattern.

Best

Tobi

[mreineck]

Hi Tobi,

many thanks for the quick reply!

I will have a look at the necessary steps on the C side to come up with a wrapper. It should not be too hard, probably similar to the FFTW way of describing multi-D arrays. The only slightly tedious thing is that all template specializations might have to be written down explicitly.

Is it the classical algorithm or is it a completely different approach?

It's the classical approach, using a slightly tweaked version of the finufft ES kernel. Kernels are evaluated via precomputed polynomials, so in principle any kernel could be used.

The special feature of the library is that, given the grid size, the number of nonuniform points and the desired accuracy, it will automatically determine a combination of kernel support, oversampling factor and concrete kernel function that results in near-optimal performance.

(More details at https://arxiv.org/abs/2010.10122; this describes a specialized application for radio interferometry, but most of the technical details also apply to the NFFT implementation.)

Do you use the same blocking technique like finufft?

Not quite. I subdivide the grid into linear/quadratic/cubic patches that fit into L1 cache (except for large kernel supports in 3D, where this is no longer possible). While processing a patch, I make a local copy of the grid data (even for the uniform-to-nonuniform case), which seems to help with parallelization.
For the NU points I compute a permutation order, such that all points belonging to one patch are processed together, but I don't actually rearrange the points, just process them in this order. (This causes somewhat nasty random memory accesses, but it seems better than actual reshuffling the data.'

Do you use FFTW and if not is your FFT faster?

I'm not using FFTW because I"d like to avoid external dependencies as much as possible.
My own FFT is faster than FFTW with FFTW_ESTIMATE, and often comparable to FFTW with FFTW_MEASURE, without FFTW's big planning overhead. For codes that are not totally dominated by FFTs this feels like a good compromise.

In case that ducc0.nufft is faster than NFFT.jl it would be great if you could explain some tricks that you applied. We just submitted a paper that describes what NFFT.jl does: https://arxiv.org/abs/2208.00049

I have seen the paper, but unfortunately didn't have enough time yet for a closer look; I certainly plan to read it carefully soon.
Concerning our own tricks, the most in-depth description is currently at https://arxiv.org/abs/2010.10122; however this does not cover technical details like template metaprogramming and SIMD usage. I'm more than happy to discuss details at any time though!

Side comment:
"does not support explicit planning, making the interface much simpler"
I never get why this makes an interface simpler. Grouping the memory in a struct plus using caching techniques is an established programming pattern.

It certainly is. But it means keeping global state in a library, which can cause all sorts of headaches. For example, you have to guarantee the thread safety of the planning functions. And everyone has slightly different ideas about the guarantees made by the library, as you have discovered, for example, with FINUFFT, if I interpret the comments in FINUFFT.jl correctly :-)
It can be done for ducc0.nufft, but I sincerely doubt it will be worth the hassle. Same with ducc0.fft...

Cheers,
Martin

[mreineck]

Sorry, I forgot to mention one thing:

If you want to run benchmarks, please be sure to use

pip3 install --no-binary ducc0 --user ducc0

or similar when installing the package to enforce local compilation. If you install the vanilla binary package, you will not benefit from CPU features like AVX etc.

[tknopp]

yes I did compile from source.

[tknopp]

Some initial benchmarks on the same computer indeed indicate that ducc0.nufft is quite fast. In 1D its slower (even slower than FINUFFT). In 2D it depends on the accuracy but I would say that overall NFFT.jl and ducc0.nufft have the same speed. In 3D, ducc0.nufft seems to be faster in the type2 transform and at a similar speed as the adjoint transform. This is quite impressive.

[mreineck]

Thanks for running the benchmarks! In 1D I expect the slowness is mostly due to the FFT part (1D FFTs are the weakest component of ducc0.fft). This could be tested by checking if the performance difference decreases with increasing number of nonuniform points.

[mreineck]

Also, in 1D, the preparation step is proportionally more expensive than in 2D and 3D. Did you test against NFFT.jl including or excluding precomputation?

[tknopp]

It is faster with and without preparation. But for a fair comparison it would be better if ducc0.fft would allow to cache the precomputation ;-)

[mreineck]

I won't be able to do that quickly ... but if you are interested in the potential speedup, you can set verbosity=1 in the calls and look at the timing values for "building index" and "parameter calculation". These are the only steps that can be done ahead of time.

[tknopp]

Oh wow, I just realize that you change the oversampling factor adaptively. This is very important information. Do you use some heuristic there? I would be very much interested in such a thing.

For the benchmarks shown on our homepage it would be great to have version with and without adaptive parameter selection. This allows to understand where the performance gain comes from.

[mreineck]

To determine the oversampling factor I'm using a simple cost model: the higher the oversampling factor, the more expensive the FFT, but the lower the kernel support and therefore the gridding/degridding cost. For all the possible oversampling factors I estimate the computation time using heuristic formulas and pick the one that is expected to run fastest.

You can already now limit the allowed oversampling factors using the sigma_min and sigma_max parameters. If you want to ensure a sigma of 2, it's probably best to use something like sigma_min=1.98, sigma_max=2.02, since I have factors tabulated in steps of 0.05.

It is possible that the code won't let you use a certain oversampling factor even though the kernel for it theoretically exists; this is because the deconvolution step might blow up the errors so much that the requested tolerance cannot be guaranteed any more. This is a somewhat obscure effect, but I have noticed it causing problems in real-world calculations.

[mreineck]

For the parameter selection code, have a look at https://gitlab.mpcdf.mpg.de/mtr/ducc/-/blob/ducc0/src/ducc0/nufft/nufft.h#L147.

[tknopp]

Thanks, this is definitely a TODO item for NFFT.jl

[tknopp]

What is quite interesting that NFFT is faster in the "spreading" routine (in 3D) but slower in "interpolation".
I have:
interpolation / spreading
0.6872 s / 0.4061 s NFFT.jl
0.1964 s / 0.5570 s ducc0

And this at the same oversampling factor kernel size. So there is definitely something we can learn when analyzing both libraries.

And your FFT faster than FFTW, which is quite a statement. It would be great if there would be a Julia package wrapping the FFT part of ducc0. Then we can make it possible to use different FFT implementations within NFFT.jl.

[tknopp]

Hmm, your first sentence seems to contradict your table ... which one is correct?

yes, sorry. I edit that.

[mreineck]

Apple M1 Max. Not sure what that has

Something else entirely (SVE), which I'm using even more blindly ;-)

It's a pity that the spreading part is the slow one; I probably have some idea how to fix that, but then AVX and friends will be slowed down, and I don"t lke the idea of having different code paths for different CPUs.

[tknopp]

Please don't take my benchmarks too serious. My Julia code is actually running through Rosetta2, i.e. the emulation layer. Not sure if ducc is compiled to native code here.

But this has convinced me that we should really work on a ducc0.nufft wrapper such that we can compare all libraries against each other.

[mreineck]

Do you have any experience with https://github.com/JuliaInterop/CxxWrap.jl? It seems very close to pybind11 which I'm using for the interface to Python. Using a wrapper like this might save quite a lot of time, and we wouldn't have to go through pure C interfaces with all their raw pointers.

[tknopp]

Experience not. But it is an established package and certainly a good option.

[mreineck]

OK, I fear that a nice Julia interface using CxxWrap.jl will take quite some time ... I have no idea how to communicate multi-D arrays properly.

For the moment, how about the following minimalistic C interface functions?

/*
ndim:    number of dimensions (1/2/3)
npoints: number of non-uniform points
shape:   points to a dense Julia array of shape(ndim,) containing the grid
         dimensions in Julia order
grid:    points to a dense Julia array of shape (2,shape)
         the leading dimension is for real and imaginary parts
coord:   points to a dense Julia array of shape(ndim,npoints)
forward: 0 ==> FFT exponent is  1
         1 ==> FFT exponent is -1
out:     points to a mutable dense Julia array of shape (2,npoints)
         the leading dimension is for real and imaginary parts
*/
void nufft_u2nu_julia_double (size_t ndim,
                              size_t npoints,
                              const size_t *shape,
                              const double *grid,
                              const double *coord,
                              int forward,
                              double epsilon,
                              size_t nthreads,
                              double *out,
                              size_t verbosity,
                              double sigma_min,
                              double sigma_max);
/*
ndim:    number of dimensions (1/2/3)
npoints: number of non-uniform points
shape:   points to a dense Julia array of shape(ndim,) containing the grid
         dimensions in Julia order
points:  points to a dense Julia array of shape (2,npoints)
         the leading dimension is for real and imaginary parts
coord:   points to a dense Julia array of shape(ndim,npoints)
forward: 0 ==> FFT exponent is  1
         1 ==> FFT exponent is -1
out:     points to a dense mutable Julia array of shape (2,shape)
         the leading dimension is for real and imaginary parts
*/

void nufft_nu2u_julia_double (size_t ndim,
                              size_t npoints,
                              const size_t *shape,
                              const double *points,
                              const double *coord,
                              int forward,
                              double epsilon,
                              size_t nthreads,
                              double *out,
                              size_t verbosity,
                              double sigma_min,
                              double sigma_max);

Sorry that I cannot use C99 complex data types, since they do not play well with the C++ standard; this is one of the very few situations where legal C99 is not legal C++ (as far as I understand). From the Julia standpoint it shouldn't make a difference though, its should be fine just to pass the pointers to the complex data.

I should be able to prepare everything necessary on the compiled side, including a simple Makefile which generates a shared library exporting these two functions.

(Equivalents for single precision data should be easy once we have these two functions working.)

[tknopp]

Yes, I have written the wrapper and can do this for ducc. But for finufft and nfft3 there are already packaged binaries that I could use them right away.

[mreineck]

Just to understand this better: is it sufficient to write an equivalent of the build_tarballs.jl that you linked, so that you can use it locally, or do we really need to have this completely experimental JLL registered at Yggdrasil, where it is precompiled for various platforms and distributed in binary form? This is pre-alpha interface code that will not be of use to anyone at the beginning, so making a public release feels a bit like environmental pollution ;-)

[mreineck]

Here is my current draft. I have no means to test it, so it's definitely still broken. The relevant tarball should be downloadable, however.

# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using BinaryBuilderBase

include(joinpath(@__DIR__, "..", "..", "platforms", "microarchitectures.jl"))

name = "ducc0"
version = v"0.26.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gitlab.mpcdf.mpg.de/mtr/ducc/-/archive/juliatest.tar.gz", "4106a45714351b152048483344d0d4741714175eea7027906be2b464406d5872")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ducc-*/julia
${CXX} ${CFLAGS} -O3 -ffast-math -I ../src/ ducc_julia.cc -Wfatal-errors -pthread -std=c++17 -fPIC -c
${CXX} ${CFLAGS} -O3 -march=native -o libducc_julia.${dlext} ducc_julia.o -Wfatal-errors -pthread -std=c++17 -shared -fPIC
install -Dvm 0755 "libducc_julia.${dlext}" "${libdir}/libducc_julia.${dlext}"
"""

# Expand for microarchitectures on x86_64 (library doesn't have CPU dispatching)
# Tests on Linux/x86_64 yielded a slow binary with avx512 for some reason, so disable that
platforms = expand_cxxstring_abis(expand_microarchitectures(supported_platforms(), ["x86_64", "avx", "avx2"]); skip=!Sys.iswindows)

augment_platform_block = """
    $(MicroArchitectures.augment)
    function augment_platform!(platform::Platform)
        # We augment only x86_64
        @static if Sys.ARCH === :x86_64
            augment_microarchitecture!(platform)
        else
            platform
        end
    end
    """

# The products that we will ensure are always built
products = [
    LibraryProduct("libducc_julia", :libducc_julia)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"8", julia_compat="1.6", augment_platform_block)

[tknopp]

I have no means to test that either. At least locally the script fails. But I think one needs to run it in Yggdrasil.

Just to understand this better: ... This is pre-alpha interface code that will not be of use to anyone at the beginning, so making a public release feels a bit like environmental pollution

That depends on your goals. I thought that you want to develop Julia bindings for ducc0 don't you? If that is not the case it obviously does not make sense to build JLL packages.

As I said previously, it would be better to build a dedicated ducc0.jl package and make it implement AbstractFFTs and AbstractNFFTs. Then ducc0.jl would also be a drop-in replacement for FFTW.jl.

And a small comment on the Julia C wrapper: It would IMHO make sense to make this just a C wrapper, which you could also test internally within the ducc0 unit tests. I don't see anything Julia specific there and it would make ducc0 accessible by many other programming languages.

[mreineck]

That depends on your goals. I thought that you want to develop Julia bindings for ducc0 don't you? If that is not the case it obviously does not make sense to build JLL packages.

OK, that's probably the crucial point. I'd be happy to support anyone who wants to develop the Julia bindings, by providing something on the C/C++ side that is easily callable from the Julia side. I don't speak a single word of Julia at the moment, and even though I'd very much like to change that, it will be a long way until I can take responsibility as a maintainer for a Julia package. This really should be done by someone with more expertise than myself.
I'll see if I can find someone with the necessary knowledge at my institute.

And a small comment on the Julia C wrapper: It would IMHO make sense to make this just a C wrapper, which you could also test internally within the ducc0 unit tests. I don't see anything Julia specific there and it would make ducc0 accessible by many other programming languages.

That is a very good point. It should be re-usable for wrappers in Fortran (yes, there has been a request for that ;-) and perhaps Rust etc.

[mreineck]

I have now released a new version which allows planned NFFTs. If only looking at the execution times of the pre-planned transforms, this provides significant speedups for 1D and 2D. In 3D the planning cost is dominated by the actual execution cost, so the effect is small there.

I'm attaching a few benchmark plots similar to the NFFT.jl ones (sorry, I still cannot compare directly with the Julia results, so Finufft will have to do for the moment).

[mreineck]

BTW, I asked around locally, but unfortunately I could not find anyone with the skills necessary to create the Julia side of the wrapping code. So I fear this must be postponed. I will keep the C wrapper up to date however and hope that someone can be found in the future!

[mreineck]

I tried to hack something together myself; if you are interested, please have a look at #107 !