dfm / python-nufft Goto Github PK

View Code? Open in Web Editor NEW

58.0 5.0 12.0 111 KB

Python bindings to a subset of the NUFFT algorithm

License: Other

Python 19.08% Fortran 80.92%

python-nufft's Introduction

Python-NUFFT

Python bindings to a subset of the NUFFT algorithm. 1D, 2D, and 3D cases are implemented.

Usage

The documentation can be found on ReadTheDocs.

To install, run python setup.py install. Then, to evaluate a type-3 FT in 1D, use nufft.nufft1d3. Assuming that you have a time series in t and y and you want to evaluate it at (angular) frequencies f:

import nufft
ft = nufft.nufft1d3(t, y, f)

You can specify your required precision using eps=1e-15. The default is 1e-15.

Authors and License

Python bindings by Dan Foreman-Mackey, Thomas Arildsen, and Marc T. Henry de Frahan but the code that actually does the work is from the Greengard lab at NYU (see the website). The Fortran code is BSD licensed and the Python bindings are MIT licensed.

python-nufft's People

Contributors

Stargazers

Watchers

Forkers

nanaln vidalalcala thomasa afcarl deepdumbo mikaem freesiemens dimitar-petrov cedmrnl i-a-morozov

python-nufft's Issues

Publish wheels to pypi

Fortran compilation error under gcc 10 on OS X

I tried to install nufft using gfortran which comes with gcc version 10, but due to the following errors it failed:

   16 |       call dcftb1 (n,c,wsave,wsave(iw1),wsave(iw2))
      |                                                                        
Error: Type mismatch in argument 'ifac' at (1); passed REAL(8) to INTEGER(4)
src/nufft/dfftpack.f:27:72:

   27 |       call dcftf1 (n,c,wsave,wsave(iw1),wsave(iw2))
      |                                                                        
Error: Type mismatch in argument 'ifac' at (1); passed REAL(8) to INTEGER(4)
src/nufft/dfftpack.f:38:72:

   38 |       call dcfti1 (n,wsave(iw1),wsave(iw2))
      |                                                                        
Error: Type mismatch in argument 'ifac' at (1); passed REAL(8) to INTEGER(4)

However if I downgrade gcc to version 6 it compiles.

Is this project still active?

I have some questions I would like to ask about the code.

Code worked once, now core dumped

I had a code which worked, now it does not, and I think it has to do something with the fortran part of the package. It is described here:
https://stackoverflow.com/questions/52392709/code-worked-once-now-core-dumped-from-nufft-algorithm-what-went-wrong
Is this only me, or is it a general issue?
Cheers:)

Installation Issue: No module named 'nufft._nufft'

I have Ubuntu 18.04 with gfortran, and I'm using conda python 3.6.5. I ran the following to install:

git clone https://github.com/dfm/python-nufft.git
cd python-nufft
python setup.py install

However, I got the following error when trying to import:

Python 3.6.5 |Anaconda, Inc.| (default, Apr 29 2018, 16:14:56) 
[GCC 7.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import nufft
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File ".../python-nufft/nufft/__init__.py", line 21, in <module>
    from .nufft import nufft1d1freqs, nufft1d1, nufft1d2, nufft1d3, nufft2d1, nufft2d2, nufft2d3, nufft3d1, nufft3d2, nufft3d3
  File ".../python-nufft/nufft/nufft.py", line 18, in <module>
    from ._nufft import (
ModuleNotFoundError: No module named 'nufft._nufft'

Solution:
I replaced .../python-nufft/nufft/ with .../python-nufft/build/lib.linux-x86_64-3.6/nufft/ and now everything works.

issues importing in python3

I don't see python3 support or lack of support mentioned anywhere, although that may be the root of this problem. In any case, I can import and use nufft with no problems in python2.7, but running in python3 gives the following error. Any suggestions how to fix this are appreciated.


Python 3.5.2 (default, Nov 23 2017, 16:37:01) 
[GCC 5.4.0 20160609] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import nufft
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.7/dist-packages/nufft/__init__.py", line 21, in <module>
    from .nufft import nufft1d1freqs, nufft1d1, nufft1d2, nufft1d3, nufft2d1, nufft2d2, nufft2d3, nufft3d1, nufft3d2, nufft3d3
  File "/usr/local/lib/python2.7/dist-packages/nufft/nufft.py", line 18, in <module>
    from ._nufft import (
ImportError: /usr/local/lib/python2.7/dist-packages/nufft/_nufft.so: undefined symbol: PyCObject_Type

benchmarks

Can you set up and document a few benchmarks to show that it is far faster than Lomb-Scargle? Those benchmarks could be described in the readme or documentation.

computing the FT for non-uniform u,v coordinates

I would like to compute the real and imaginary parts of V(u,v) from I(x,y) (a 2D image with a known pixel scale) for non-uniformly spaced u,v coordinates, where,

V(u,v) = sum_i^I sum_j^J I(x_i,y_j) exp(-2 * pi * i * (x_i * u + y_j * v))

where the u,v data are in units of wavelengths (or inverse radians) and x,y are in units of radians.

Is there a way to perform this calculation with nufft or does it only works for when x, y are non-uniformly space?

Below is an example script of the discrete FT (which becomes very slow for a large number of u,v points) that I want to compare against nufft.

import numpy as np
import matplotlib.pyplot as plt
from astropy import units

def z(x, y, sigma_x, sigma_y):
    return 1.0 / (2.0 * np.pi * sigma_x * sigma_y) * np.exp(-(x**2.0 / (2.0 * sigma_x**2) + y**2.0 / (2.0 * sigma_y**2.0)))


size = 5.0 # The angular size of the image in arcsec.
n_pixels = int(2**6.0)

# compute the pixel scale
pixel_scale = size / n_pixels

# generate the x,y regular grid
x_rad = np.linspace(
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)
y_rad = np.linspace(
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)
x_rad, y_rad = np.meshgrid(
    x_rad,
    y_rad
)

# make an image (example 1)
sigma_x = 0.5 # in units of arcsec
sigma_y = 0.5 # in units of arcsec
image = z(
    x=x_rad,
    y=y_rad,
    sigma_x=sigma_x * units.arcsec.to(units.rad),
    sigma_y=sigma_y * units.arcsec.to(units.rad)
)

# Generate the non-uniform u,v coordinates (in units of inverse radians or wavelengths)
N = 100
u = np.random.uniform(
    -1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    +1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    N)
v = np.random.uniform(
    -1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    +1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    N)

# Generate visibilities using dFT
visibilities = np.zeros(
    shape=(u.shape[-1]),
    dtype="complex"
)
image_flatten = np.ndarray.flatten(image)
x_rad_flatten = np.ndarray.flatten(x_rad)
y_rad_flatten = np.ndarray.flatten(y_rad)
for j in range(u.shape[-1]):
    for i in range(image_flatten.shape[-1]):
        visibilities[j] += image_flatten[i] * np.exp(-2j * np.pi * (x_rad_flatten[i] * u[j] + y_rad_flatten[i] * v[j]))
visibilities_real__dFT = visibilities.real
visibilities_imag__dFT = visibilities.imag

plt.figure()
plt.scatter(
    visibilities_real__dFT,
    visibilities_imag__dFT,
    marker="o",
    alpha=0.85,
    label="dFT"
)
plt.xlabel("Real", fontsize=15)
plt.ylabel("Imag", fontsize=15)
plt.legend()
plt.show()