Added test_holo_from_correlations

test.py

@@ -344,3 +344,62 @@ def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
         res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
     res1 /= n
     assert_(_l2error(res0, res1) < 1e-4)
+
+@pmp('nxdirty', (32,16))
+@pmp('nydirty', (32,64))
+@pmp("nrow", (100,1000))
+@pmp("nchan", (1,))
+@pmp('epsilon', (1e-3,1e-6))
+def test_holo_from_correlations(nxdirty, nydirty, nchan, nrow, epsilon):
+    np.random.seed(420)
+    uvw = np.random.randn(nrow, 3)
+    uvw[:, 2] = 0
+    cell = 0.75/(2*np.abs(uvw).max())
+    conf = ng.GridderConfig(nxdirty=nxdirty,
+                            nydirty=nydirty,
+                            epsilon=epsilon,
+                            pixsize_x=cell,
+                            pixsize_y=cell)
+    f0 = 3e8
+    freq = f0 + np.arange(nchan)*(f0/nchan)
+    bl = ng.Baselines(coord=uvw, freq=freq)
+    flags = np.zeros((nrow, nchan), dtype=np.bool)
+    idx = ng.getIndices(bl, conf, flags)
+
+    # random vector to apply to
+    tmp = (np.random.randn(2*nxdirty, 2*nydirty) + 
+           1.0j*np.random.randn(2*nxdirty, 2*nydirty))
+    
+    # result of implicit operator
+    res1 = ng.apply_holo(bl, conf, idx, tmp)
+
+    # get explicit operator
+    W = conf.W()
+    holo_d = {}
+    for i in range(0, W):
+        for j in range(0, W):
+            # add entry to dict
+            holo_d[(i,j)] = np.fft.fftshift(ng.get_correlations(bl, conf, idx, du=i, dv=j)).flatten()
+            if i:
+                holo_d[(i,-j)] = np.fft.fftshift(ng.get_correlations(bl, conf, idx, du=i, dv=-j)).flatten()
+
+
+    # apply to random vector
+    Ntot = 2*nxdirty * 2*nydirty
+    x = np.fft.fftshift(tmp).flatten()
+    res2 = np.zeros(Ntot, dtype=np.complex128)
+    for key in holo_d:
+        i, j = key
+        shift = i*2*nydirty + j
+        diag = holo_d[key]
+        if shift:
+            # upper diagonal
+            res2[0:Ntot - shift] += diag[0:Ntot - shift] * x[shift::]
+            # lower diagonal
+            res2[shift::] += diag[0:Ntot - shift] * x[0:Ntot - shift]
+        else:
+            res2 += diag * x
+
+    res2 = np.fft.ifftshift(res2.reshape(2*nxdirty, 2*nydirty))
+
+    assert_array_almost_equal(res1, res2, decimal=8)
\ No newline at end of file