I have a formula that is completely defined in numpy and scipy and I want to define it using theano ops for use in pymc3 and I want to be able to use HMC or NUTS.
The inputs are v \sin i which is a float, f which are frequencies from an associated Fourier transform, and F_\lambda which are flux densities (light energy per density per wavelength)
Here is the formula:
This is essentially a broadening kernel used in astronomy. We want to implement this convolution by multiplying the kernel in the Fourier domain, so the python implementation looks like this
import numpy as np from scipy.special import j1 def rotational_broaden(wavelengths, fluxes, vsini): # This is the frequency spacing with maximum information dv = calculate_dv(wavelengths) freq = np.fft.rfftfreq(fluxes.shape[-1], d=dv) flux_ff = np.fft.rfft(fluxes) # ignore the 0th frequency ub = 2. * np.pi * vsini * freq[1:] # Calculate the stellar broadening kernel (Gray 2008) sb = j1(ub) / ub - 3 * np.cos(ub) / (2 * ub ** 2) + 3. * np.sin(ub) / (2 * ub ** 3) flux_ff[:, 1:] *= sb flux_final = np.fft.irfft(flux_ff, n=fluxes.shape[-1]) return flux_final
I am sort of able to incorporate this using theano, for instance there exists
tt.j1. I’m not sure how to work with the Fourier transforms, though, the theano implementations don’t match the numpy versions. Also is the issue with the frequency spacing and only defining the kernel from
freq[1:] to avoid f=0 in the denominator of the kernel.