diff --git a/gridder_cxx.h b/gridder_cxx.h index db1d99765d70ac6c5595fcc0f3a2527fe09254b5..fe6f687b90aa32c8a29a23f58af5019b30e29431 100644 --- a/gridder_cxx.h +++ b/gridder_cxx.h @@ -375,7 +375,7 @@ class ES_Kernel public: ES_Kernel(size_t supp_, size_t nthreads) - : beta(2.3*supp_), p(int(1.5*supp_+2)), supp(supp_) + : beta(get_beta(supp_)*supp_), p(int(1.5*supp_+2)), supp(supp_) { legendre_prep(2*p,x,wgt,nthreads); psi=x; @@ -392,16 +392,24 @@ class ES_Kernel tmp += wgt[i]*psi[i]*cos(pi*supp*v*x[i]); return 1./(supp*tmp); } + static double get_beta(size_t supp) + { + static const vector opt_beta {-1, 0.14, 1.70, 2.08, 2.205, 2.26, + 2.29, 2.307, 2.316, 2.3265, 2.3324, 2.282, 2.294, 2.304, 2.3138, 2.317}; + myassert(supp maxmaperr { 1e8, 0.32, 0.021, 6.2e-4, - 1.08e-5, 1.25e-7, 8.25e-10, 5.70e-12, 1.22e-13, 2.48e-15, 4.82e-17, - 6.74e-19, 5.41e-21, 4.41e-23, 7.88e-25, 3.9e-26 }; + static const vector maxmaperr { 1e8, 0.19, 2.98e-3, 5.98e-5, + 1.11e-6, 2.01e-8, 3.55e-10, 5.31e-12, 8.81e-14, 1.34e-15, 2.17e-17, + 2.12e-19, 2.88e-21, 3.92e-23, 8.21e-25, 7.13e-27 }; double epssq = epsilon*epsilon; for (size_t i=1; imaxmaperr[i]) return i; - myfail("requested epsilon too small - minimum is 2e-13"); + myfail("requested epsilon too small - minimum is 1e-13"); } }; @@ -552,7 +560,7 @@ class GridderConfig psx(pixsize_x), psy(pixsize_y), supp(ES_Kernel::get_supp(epsilon)), nsafe((supp+1)/2), nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)), - beta(2.3*supp), + beta(ES_Kernel::get_beta(supp)*supp), cfu(nx_dirty), cfv(ny_dirty), nthreads(nthreads_), ushift(supp*(-0.5)+1+nu), vshift(supp*(-0.5)+1+nv), maxiu0((nu+nsafe)-supp), maxiv0((nv+nsafe)-supp) @@ -1214,16 +1222,14 @@ template void apply_wcorr(const GridderConfig &gconf, } } -template void wstack_common( - const GridderConfig &gconf, const Serv &srv, double &wmin, - double &dw, size_t &nplanes, vector &nvis_plane, - vector &minplane, size_t verbosity) +template void wminmax(const GridderConfig &gconf, + const Serv &srv, double &wmin, double &wmax) { size_t nvis = srv.Nvis(); size_t nthreads = gconf.Nthreads(); - wmin=1e38; - double wmax=-1e38; + wmin= 1e38; + wmax=-1e38; // FIXME maybe this can be done more intelligently #pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax) for (size_t ipart=0; ipart void wstack_common( wmin = min(wmin,wval); wmax = max(wmax,wval); } + } + +template void wstack_common( + const GridderConfig &gconf, const Serv &srv, double &wmin, + double &dw, size_t &nplanes, vector &nvis_plane, + vector &minplane, size_t verbosity) + { + size_t nvis = srv.Nvis(); + size_t nthreads = gconf.Nthreads(); + double wmax; + + wminmax(gconf, srv, wmin, wmax); if (verbosity>0) cout << "Using " << nthreads << " threads" << endl; if (verbosity>0) cout << "W range: " << wmin << " to " << wmax << endl; double x0 = -0.5*gconf.Nxdirty()*gconf.Pixsize_x(), diff --git a/test.py b/test.py index 6b154ce9475f178723582b86d97b27bbd06827ba..2e8b107da242c8908af28d4f5519defb7ac857bc 100644 --- a/test.py +++ b/test.py @@ -15,7 +15,7 @@ def _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow): epsilon=epsilon, pixsize_x=0.368*pixsize, pixsize_y=pixsize) - speedoflight, f0 = 3e8, 1e9 + speedoflight, f0 = 299792458., 1e9 freq = f0 + np.arange(nchan)*(f0/nchan) uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight) baselines = ng.Baselines(coord=uvw, freq=freq) @@ -24,26 +24,6 @@ def _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow): return baselines, conf, idx -def _dft(uvw, vis, conf, apply_w): - shp = (conf.Nxdirty(), conf.Nydirty()) - x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in shp], - indexing='ij') - x *= conf.Pixsize_x() - y *= conf.Pixsize_y() - dirty = np.zeros(shp, dtype=np.complex128) - if apply_w: - n = np.sqrt(1-x**2-y**2) - nm1 = n-1 - for ii, vv in enumerate(vis): - phase = x*uvw[ii, 0] + y*uvw[ii, 1] - if apply_w: - phase -= uvw[ii, 2]*nm1 - dirty += (vv*np.exp(2j*np.pi*phase)).real - if apply_w: - return dirty/n - return dirty - - def _l2error(a, b): return np.sqrt(np.sum(np.abs(a-b)**2)/np.maximum(np.sum(np.abs(a)**2), np.sum(np.abs(b)**2))) @@ -90,7 +70,7 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, nrow, nchan, epsilon, singleprec np.random.seed(42) pixsizex = np.pi/180/60/nxdirty*0.2398 pixsizey = np.pi/180/60/nxdirty - speedoflight, f0 = 3e8, 1e9 + speedoflight, f0 = 299792458., 1e9 freq = f0 + np.arange(nchan)*(f0/nchan) uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizey*f0/speedoflight) ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) @@ -125,7 +105,7 @@ def test_ms2dirty_against_wdft(nxdirty, nydirty, nrow, nchan, epsilon, singlepre np.random.seed(40) pixsizex = fov*np.pi/180/nxdirty pixsizey = fov*np.pi/180/nydirty*1.1 - speedoflight, f0 = 3e8, 1e9 + speedoflight, f0 = 299792458., 1e9 freq = f0 + np.arange(nchan)*(f0/nchan) uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizex*f0/speedoflight) ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) @@ -185,10 +165,9 @@ def test_hoisted_grid_allocation(nxdirty, nydirty, nrow, nchan, epsilon): gconf = ng.GridderConfig(1024, 1024, 2e-13, 2.0, 2.0) freq = np.linspace(.856e9, 2*.856e9, nchan) f0 = freq[freq.shape[0]//2] - speedoflight = 3e8 pixsize = np.pi/180/60/nxdirty # assume 1 arcmin FOV - speedoflight = 3e8 + speedoflight = 299792458. uvw = (np.random.rand(nrow, 3)-0.5) / (pixsize*f0/speedoflight) baselines = ng.Baselines(coord=uvw, freq=freq) gconf = ng.GridderConfig(nxdirty=nxdirty, nydirty=nydirty, @@ -327,22 +306,6 @@ def test_no_correlation(nrow, nchan, epsilon, weight): ng.get_correlations(bl, conf, idx, uu, vv, wgt) -@pmp('epsilon', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13]) -@pmp('nxdirty', [2, 12, 128]) -@pmp('nydirty', [2, 4, 12, 128]) -@pmp("nrow", (10, 100)) -@pmp("nchan", (1, 10)) -def test_against_dft(nxdirty, nydirty, epsilon, nchan, nrow): - np.random.seed(42) - bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow) - ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) - vis = bl.ms2vis(ms, idx) - uvw = bl.effectiveuvw(idx) - res0 = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, idx, vis)).real - res1 = _dft(uvw, vis, conf, False) - assert_(_l2error(res0, res1) < epsilon) - - @pmp('nxdirty', [16, 64]) @pmp('nydirty', [64]) @pmp("nrow", (10, 100)) @@ -357,7 +320,7 @@ def test_wstack_against_wdft(nxdirty, nydirty, nchan, nrow, fov): epsilon=epsilon, pixsize_x=pixsize, pixsize_y=0.5*pixsize) - speedoflight, f0 = 3e8, 1e9 + speedoflight, f0 = 299792458., 1e9 freq = f0 + np.arange(nchan)*(f0/nchan) uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight) bl = ng.Baselines(coord=uvw, freq=freq) @@ -365,9 +328,9 @@ def test_wstack_against_wdft(nxdirty, nydirty, nchan, nrow, fov): idx = ng.getIndices(bl, conf, flags) ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) vis = bl.ms2vis(ms, idx) - uvw = bl.effectiveuvw(idx) res0 = np.zeros((nxdirty, nydirty), dtype=np.float64) - mi, ma = np.min(np.abs(uvw[:, 2])), np.max(np.abs(uvw[:, 2])) + mi = np.min(np.abs(uvw[:, 2]))*freq[0]/speedoflight + ma = np.max(np.abs(uvw[:, 2]))*freq[-1]/speedoflight nplanes = 10000 ws = mi + np.arange(nplanes)*(ma-mi)/(nplanes-1) for ii in range(len(ws)-1): @@ -378,7 +341,8 @@ def test_wstack_against_wdft(nxdirty, nydirty, nchan, nrow, fov): continue dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, iind, bl.ms2vis(ms, iind))) res0 += conf.apply_wscreen(dd, 0.5*(ws[ii+1]+ws[ii]), adjoint=True).real - res1 = _dft(uvw, vis, conf, True) + res1 = explicit_gridder(uvw, freq, ms, None, nxdirty, nydirty, pixsize, + 0.5*pixsize, True) assert_(_l2error(res0, res1) < 1e-4) # Very inaccurate @@ -396,7 +360,7 @@ def test_wplane_against_wdft(nxdirty, nydirty, nchan, nrow, fov): epsilon=epsilon, pixsize_x=pixsize, pixsize_y=0.2*pixsize) - speedoflight, f0 = 3e8, 1e9 + speedoflight, f0 = 299792458., 1e9 freq = f0 + np.arange(nchan)*(f0/nchan) uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight) uvw[:, 2] = uvw[0, 2] # use exactly one random w @@ -405,12 +369,12 @@ def test_wplane_against_wdft(nxdirty, nydirty, nchan, nrow, fov): idx = ng.getIndices(bl, conf, flags) ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) vis = bl.ms2vis(ms, idx) - uvw = bl.effectiveuvw(idx) - w = np.min(uvw[:, 2]) + w = np.min(uvw[:, 2])*freq[0]/speedoflight iind = ng.getIndices(bl, conf, flags) dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, iind, bl.ms2vis(ms, iind))) res0 = conf.apply_wscreen(dd, w, adjoint=True).real - res1 = _dft(uvw, vis, conf, True) + res1 = explicit_gridder(uvw, freq, ms, None, nxdirty, nydirty, pixsize, + 0.2*pixsize, True) assert_(_l2error(res0, res1) < epsilon) @@ -428,7 +392,7 @@ def test_wgridder_against_wdft(nxdirty, nydirty, nchan, nrow, fov, epsilon): epsilon=epsilon, pixsize_x=pixsize, pixsize_y=0.5*pixsize) - speedoflight, f0 = 3e8, 1e9 + speedoflight, f0 = 299792458., 1e9 freq = f0 + np.arange(nchan)*(f0/nchan) uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight) bl = ng.Baselines(coord=uvw, freq=freq) @@ -436,9 +400,9 @@ def test_wgridder_against_wdft(nxdirty, nydirty, nchan, nrow, fov, epsilon): idx = ng.getIndices(bl, conf, flags) ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) vis = bl.ms2vis(ms, idx) - uvw = bl.effectiveuvw(idx) res0 = ng.vis2dirty(bl, conf, idx, vis, do_wstacking=True) - res1 = _dft(uvw, vis, conf, True) + res1 = explicit_gridder(uvw, freq, ms, None, nxdirty, nydirty, pixsize, + 0.5*pixsize, True) assert_(_l2error(res0, res1) < epsilon) @@ -473,4 +437,4 @@ def test_holo_from_correlations(nxdirty, nydirty, nchan, nrow, epsilon): res0 += np.roll(d1[ii]*grid, du, axis=0) res0 += d2*grid res1 = ng.apply_holo(bl, conf, idx, grid) - assert_allclose(res0, res1, rtol=1e-13) + assert_allclose(res0, res1, rtol=1e-13, atol=1e-13)