Merge branch 'better_beta' into 'master'

Better beta

See merge request !24

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<double> 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<opt_beta.size(), "bad support size");
+      return opt_beta[supp];
+      }
+
     static size_t get_supp(double epsilon)
       {
-      static const vector<double> 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<double> 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; i<maxmaperr.size(); ++i)
         if (epssq>maxmaperr[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<typename T> void apply_wcorr(const GridderConfig &gconf,
     }
   }
 
-template<typename Serv> void wstack_common(
-  const GridderConfig &gconf, const Serv &srv, double &wmin,
-  double &dw, size_t &nplanes, vector<size_t> &nvis_plane,
-  vector<int> &minplane, size_t verbosity)
+template<typename Serv> 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<nvis; ++ipart)
@@ -1232,6 +1238,18 @@ template<typename Serv> void wstack_common(
     wmin = min(wmin,wval);
     wmax = max(wmax,wval);
     }
+  }
+
+template<typename Serv> void wstack_common(
+  const GridderConfig &gconf, const Serv &srv, double &wmin,
+  double &dw, size_t &nplanes, vector<size_t> &nvis_plane,
+  vector<int> &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(),

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)