multiplicatively convolutional fast integral transforms
Switch branches/tags
Nothing to show
Clone or download
Latest commit e58e21b Sep 17, 2018
Permalink
Failed to load latest commit information.
doc Edit Sep 17, 2018
mcfit Bump version Sep 14, 2018
.gitignore Reset github commit Nov 15, 2017
LICENSE Initial commit Mar 10, 2017
README.rst Edit Sep 17, 2018
runtests.py patch from Yu Mar 18, 2017
setup.py Add __version__ Jun 26, 2018

README.rst

Multiplicatively Convolutional Fast Integral Transforms

mcfit computes integral transforms of the form

G(y) = \int_0^\infty F(x) K(xy) \frac{dx}x

where F(x) is the input function, G(y) is the output function, and K(xy) is the integral kernel. One is free to scale all three functions by a power law

g(y) = \int_0^\infty f(x) k(xy) \frac{dx}x

where f(x)=x^{-q}F(x), g(y)=y^q G(y), and k(t)=t^q K(t). And q is a tilt parameter serving to shift power of x between the input function and the kernel.

mcfit implements the FFTLog algorithm. The idea is to take advantage of the convolution theorem in \ln x and \ln y. It approximates the input function with truncated Fourier series over one period of the periodic approximant, and use the exact Fourier transform of the kernel. One can calculate the latter analytically as a Mellin transform. This algorithm is optimal when the input function is smooth in \ln x, and is ideal for oscillatory kernels with input spanning a wide range in \ln x.

Examples

One can perform the following pair of Hankel transforms

e^{-y} &= \int_0^\infty (1+x^2)^{-\frac32} J_0(xy) x dx \\
(1+y^2)^{-\frac32} &= \int_0^\infty e^{-y} J_0(xy) x dx

easily as follows

def F_fun(x): return 1 / (1 + x*x)**1.5
def G_fun(y): return numpy.exp(-y)

from mcfit import Hankel

x = numpy.logspace(-3, 3, num=60, endpoint=False)
F = F_fun(x)
H = Hankel(x)
y, G = H(F)
numpy.allclose(G, G_fun(y), rtol=1e-8, atol=1e-8)

y = numpy.logspace(-4, 2, num=60, endpoint=False)
G = G_fun(y)
H_inv = Hankel(y)
x, F = H_inv(G)
numpy.allclose(F, F_fun(x), rtol=1e-10, atol=1e-10)

Cosmologists often need to transform a power spectrum to its correlation function

from mcfit import P2xi
k, P = numpy.loadtxt('P.txt', unpack=True)
r, xi = P2xi(k)(P)

and the other way around

from mcfit import xi2P
r, xi = numpy.loadtxt('xi.txt', unpack=True)
k, P = xi2P(r)(xi)

Similarly for the quadrupoles

k, P2 = numpy.loadtxt('P2.txt', unpack=True)
r, xi2 = P2xi(k, l=2)(P2)

Also useful to the cosmologists is the tool below that computes the variance of the overdensity field as a function of radius, from which \sigma_8 can be interpolated.

R, var = TophatVar(k)(P)
from scipy.interpolate import CubicSpline
varR = CubicSpline(R, var)
sigma8 = numpy.sqrt(varR(8))