Commit fc601ac0 authored by Martin Reinecke's avatar Martin Reinecke

experimental speedup

parent 9e84bd2f
......@@ -766,27 +766,65 @@ class GridderConfig
{
checkShape(dirty.shape(), {nx_dirty, ny_dirty});
checkShape(dirty2.shape(), {nx_dirty, ny_dirty});
double x0 = -0.5*nx_dirty*psx,
y0 = -0.5*ny_dirty*psy;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (size_t i=0; i<=nx_dirty/2; ++i)
#if 0
if ((nx_dirty==ny_dirty) && (psx==psy)) // extra symmetry
{
double p0 = -0.5*nx_dirty*psx;
#pragma omp parallel for num_threads(nthreads) schedule(dynamic,10)
for (size_t i=0; i<=nx_dirty/2; ++i)
{
double fx = p0+i*psx;
fx *= fx;
for (size_t j=0; j<=i; ++j)
{
double fy = p0+j*psy;
auto ws = complex<T>(wscreen(fx, fy*fy, w, adjoint, divide_by_n));
dirty2(i,j) = dirty(i,j)*ws; // lower left
if (j!=i) dirty2(j,i) = dirty(j,i)*ws; // lower left
size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
if ((i>0)&&(i<i2))
{
dirty2(i2,j) = dirty(i2,j)*ws; // lower right
if (i!=j) dirty2(j,i2) = dirty(j,i2)*ws; // lower right
if ((j>0)&&(j<j2))
{
dirty2(i2,j2) = dirty(i2,j2)*ws; // upper right
if (i!=j) dirty2(j2,i2) = dirty(j2,i2)*ws; // upper right
}
}
if ((j>0)&&(j<j2))
{
dirty2(i,j2) = dirty(i,j2)*ws; // upper left
if (i!=j) dirty2(j2,i) = dirty(j2,i)*ws; // upper left
}
}
}
}
else
#endif
{
double fx = x0+i*psx;
fx *= fx;
for (size_t j=0; j<=ny_dirty/2; ++j)
double x0 = -0.5*nx_dirty*psx,
y0 = -0.5*ny_dirty*psy;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (size_t i=0; i<=nx_dirty/2; ++i)
{
double fy = y0+j*psy;
auto ws = complex<T>(wscreen(fx, fy*fy, w, adjoint, divide_by_n));
dirty2(i,j) = dirty(i,j)*ws; // lower left
size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
if ((i>0)&&(i<i2))
double fx = x0+i*psx;
fx *= fx;
for (size_t j=0; j<=ny_dirty/2; ++j)
{
dirty2(i2,j) = dirty(i2,j)*ws; // lower right
double fy = y0+j*psy;
auto ws = complex<T>(wscreen(fx, fy*fy, w, adjoint, divide_by_n));
dirty2(i,j) = dirty(i,j)*ws; // lower left
size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
if ((i>0)&&(i<i2))
{
dirty2(i2,j) = dirty(i2,j)*ws; // lower right
if ((j>0)&&(j<j2))
dirty2(i2,j2) = dirty(i2,j2)*ws; // upper right
}
if ((j>0)&&(j<j2))
dirty2(i2,j2) = dirty(i2,j2)*ws; // upper right
dirty2(i,j2) = dirty(i,j2)*ws; // upper left
}
if ((j>0)&&(j<j2))
dirty2(i,j2) = dirty(i,j2)*ws; // upper left
}
}
}
......
......@@ -460,15 +460,18 @@ def test_wplane_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
@pmp("nrow", (10, 100))
@pmp("nchan", (1, 10))
@pmp("fov", (1, 10))
@pmp("yscale", (0.5, 1.))
@pmp("epsilon", (1e-2, 1e-7, 2e-12))
def test_wgridder_against_wdft(nxdirty, nydirty, nchan, nrow, fov, epsilon):
def test_wgridder_against_wdft(nxdirty, nydirty, yscale, nchan, nrow, fov, epsilon):
np.random.seed(40)
pixsize = fov*np.pi/180/nxdirty
pixsize_x = pixsize
pixsize_y = yscale*pixsize
conf = ng.GridderConfig(nxdirty=nxdirty,
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=pixsize,
pixsize_y=0.5*pixsize)
pixsize_x=pixsize_x,
pixsize_y=pixsize_y)
speedoflight, f0 = 299792458., 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
......@@ -478,8 +481,9 @@ def test_wgridder_against_wdft(nxdirty, nydirty, nchan, nrow, fov, epsilon):
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
vis = bl.ms2vis(ms, idx)
res0 = ng.vis2dirty(bl, conf, idx, vis, do_wstacking=True)
res1 = explicit_gridder(uvw, freq, ms, None, nxdirty, nydirty, pixsize,
0.5*pixsize, True)
res1 = explicit_gridder(uvw, freq, ms, None, nxdirty, nydirty, pixsize_x,
pixsize_y, True)
print(_l2error(res0, res1)/epsilon)
assert_(_l2error(res0, res1) < epsilon)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment