From 9aaf87833d273bd5c4acd7463a51ffd30881deaa Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Mon, 27 May 2019 11:11:43 +0200 Subject: [PATCH 01/21] first try --- demos/bench_gridder.py | 7 +-- nifty5/library/gridder.py | 93 ++++++++++++++++----------------- test/test_operators/test_nft.py | 15 ++---- 3 files changed, 51 insertions(+), 64 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 5dd68c1c..0a0a7f67 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -27,12 +27,9 @@ for ii in range(10, 23): img = ift.from_global_data(uvspace, img) t0 = time() - GM = ift.GridderMaker(uvspace, eps=1e-7) - idx = GM.getReordering(uv) - uv = uv[idx] - vis = vis[idx] + GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv) vis = ift.from_global_data(visspace, vis) - op = GM.getFull(uv).adjoint + op = GM.getFull().adjoint t1 = time() op(img).to_global_data() t2 = time() diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index a1f1dcda..0d416483 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -26,52 +26,45 @@ from ..sugar import from_global_data, makeDomain class GridderMaker(object): - def __init__(self, domain, eps=2e-13): - from nifty_gridder import get_w - domain = makeDomain(domain) - if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or - not len(domain.shape) == 2): - raise ValueError("need domain with exactly one 2D RGSpace") - nu, nv = domain.shape - if nu % 2 != 0 or nv % 2 != 0: - raise ValueError("dimensions must be even") - nu2, nv2 = 2*nu, 2*nv - w = get_w(eps) - nsafe = (w+1)//2 - nu2 = max([nu2, 2*nsafe]) - nv2 = max([nv2, 2*nsafe]) - - oversampled_domain = RGSpace( - [nu2, nv2], distances=[1, 1], harmonic=False) - - self._eps = eps - self._rest = _RestOperator(domain, oversampled_domain, eps) - - def getReordering(self, uv): - from nifty_gridder import peanoindex - nu2, nv2 = self._rest._domain.shape - return peanoindex(uv, nu2, nv2) - - def getGridder(self, uv): - return RadioGridder(self._rest.domain, self._eps, uv) + def __init__(self, dirty_domain, uv, eps=2e-13): + import nifty_gridder + dirty_domain = makeDomain(dirty_domain) + if (len(dirty_domain) != 1 or + not isinstance(dirty_domain[0], RGSpace) or + not len(dirty_domain.shape) == 2): + raise ValueError("need dirty_domain with exactly one 2D RGSpace") + bl = nifty_gridder.Baselines(uv, np.array([1.])); + nxdirty, nydirty = dirty_domain.shape + gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) + nu = gconf.Nu() + nv = gconf.Nv() + idx = bl.getIndices() + idx = gconf.reorderIndices(idx, bl) + + grid_domain = RGSpace([nu, nv], distances=[1, 1], harmonic=False) + + self._rest = _RestOperator(dirty_domain, grid_domain, gconf) + self._gridder = RadioGridder(grid_domain, bl, gconf, idx) + + def getGridder(self): + return self._gridder def getRest(self): return self._rest - def getFull(self, uv): - return self.getRest() @ self.getGridder(uv) + def getFull(self): + return self.getRest() @ self._gridder class _RestOperator(LinearOperator): - def __init__(self, domain, oversampled_domain, eps): - from nifty_gridder import correction_factors - self._domain = makeDomain(oversampled_domain) - self._target = domain - nu, nv = domain.shape - nu2, nv2 = oversampled_domain.shape - - fu = correction_factors(nu2, nu//2+1, eps) - fv = correction_factors(nv2, nv//2+1, eps) + def __init__(self, dirty_domain, grid_domain, gconf): + import nifty_gridder + self._domain = makeDomain(grid_domain) + self._target = makeDomain(dirty_domain) + self._gconf = gconf + fu = gconf.U_corrections() + fv = gconf.V_corrections() + nu, nv = dirty_domain.shape # compute deconvolution operator rng = np.arange(nu) k = np.minimum(rng, nu-rng) @@ -82,6 +75,7 @@ class _RestOperator(LinearOperator): self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): + import nifty_gridder self._check_input(x, mode) nu, nv = self._target.shape res = x.to_global_data() @@ -103,20 +97,21 @@ class _RestOperator(LinearOperator): class RadioGridder(LinearOperator): - def __init__(self, target, eps, uv): - self._domain = DomainTuple.make( - UnstructuredDomain((uv.shape[0],))) - self._target = DomainTuple.make(target) + def __init__(self, grid_domain, bl, gconf, idx): + self._domain = DomainTuple.make(UnstructuredDomain((idx.shape[0],))) + self._target = DomainTuple.make(grid_domain) + self._bl = bl + self._gconf = gconf + self._idx = idx self._capability = self.TIMES | self.ADJOINT_TIMES - self._eps = float(eps) - self._uv = uv # FIXME: should we write-protect this? def apply(self, x, mode): - from nifty_gridder import to_grid, from_grid + import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: - nu2, nv2 = self._target.shape - res = to_grid(self._uv, x.to_global_data(), nu2, nv2, self._eps) + res = nifty_gridder.ms2grid( + self._bl, self._gconf, self._idx, x.to_global_data().reshape((-1,1))) else: - res = from_grid(self._uv, x.to_global_data(), self._eps) + res = nifty_gridder.grid2ms( + self._bl, self._gconf, self._idx, x.to_global_data()) return from_global_data(self._tgt(mode), res) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 0e137e9c..c4946a16 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -39,13 +39,10 @@ def test_gridding(nu, nv, N, eps): vis = np.random.randn(N) + 1j*np.random.randn(N) # Nifty - GM = ift.GridderMaker(ift.RGSpace((nu, nv)), eps=eps) - # re-order for performance - idx = GM.getReordering(uv) - uv, vis = uv[idx], vis[idx] + GM = ift.GridderMaker(ift.RGSpace((nu, nv)), eps=eps, uv=uv) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) - Op = GM.getFull(uv) + Op = GM.getFull() pynu = Op(vis2).to_global_data() # DFT x, y = np.meshgrid( @@ -63,14 +60,12 @@ def test_gridding(nu, nv, N, eps): def test_build(nu, nv, N, eps): dom = ift.RGSpace([nu, nv]) uv = np.random.rand(N, 2) - 0.5 - GM = ift.GridderMaker(dom, eps=eps) + GM = ift.GridderMaker(dom, eps=eps, uv=uv) # re-order for performance - idx = GM.getReordering(uv) - uv = uv[idx] - R0 = GM.getGridder(uv) + R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 - RF = GM.getFull(uv) + RF = GM.getFull() # Consistency checks flt = np.float64 -- GitLab From 3b88ba48b78d1a751b94b715efe9482e233da1d5 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Mon, 27 May 2019 22:22:37 +0200 Subject: [PATCH 02/21] seems to work. needs lots of polishing --- demos/bench_gridder.py | 11 +++++---- nifty5/library/gridder.py | 43 +++++++++------------------------ test/test_operators/test_nft.py | 15 ++++++------ 3 files changed, 26 insertions(+), 43 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 0a0a7f67..99459b2c 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -9,13 +9,13 @@ np.random.seed(40) N0s, a0s, b0s, c0s = [], [], [], [] -for ii in range(10, 23): - nu = 1024 - nv = 1024 +for ii in range(10, 26): + nu = 2048 + nv = 2048 N = int(2**ii) print('N = {}'.format(N)) - uv = np.random.rand(N, 2) - 0.5 + uvw = np.random.rand(N, 3) - 0.5 vis = np.random.randn(N) + 1j*np.random.randn(N) uvspace = ift.RGSpace((nu, nv)) @@ -27,7 +27,7 @@ for ii in range(10, 23): img = ift.from_global_data(uvspace, img) t0 = time() - GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv) + GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, channel_fact=np.array([1.])) vis = ift.from_global_data(visspace, vis) op = GM.getFull().adjoint t1 = time() @@ -35,6 +35,7 @@ for ii in range(10, 23): t2 = time() op.adjoint(vis).to_global_data() t3 = time() + print(t2-t1, t3-t2) N0s.append(N) a0s.append(t1 - t0) b0s.append(t2 - t1) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index 0d416483..c7129498 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -26,20 +26,21 @@ from ..sugar import from_global_data, makeDomain class GridderMaker(object): - def __init__(self, dirty_domain, uv, eps=2e-13): + def __init__(self, dirty_domain, uvw, channel_fact, eps=2e-13): import nifty_gridder dirty_domain = makeDomain(dirty_domain) if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace) or not len(dirty_domain.shape) == 2): raise ValueError("need dirty_domain with exactly one 2D RGSpace") - bl = nifty_gridder.Baselines(uv, np.array([1.])); + if channel_fact.ndim != 1: + raise ValueError("channel_fact must be a 1D array") + bl = nifty_gridder.Baselines(uvw, channel_fact); nxdirty, nydirty = dirty_domain.shape gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) nu = gconf.Nu() nv = gconf.Nv() - idx = bl.getIndices() - idx = gconf.reorderIndices(idx, bl) + idx = nifty_gridder.getIndices(bl, gconf) grid_domain = RGSpace([nu, nv], distances=[1, 1], harmonic=False) @@ -62,43 +63,23 @@ class _RestOperator(LinearOperator): self._domain = makeDomain(grid_domain) self._target = makeDomain(dirty_domain) self._gconf = gconf - fu = gconf.U_corrections() - fv = gconf.V_corrections() - nu, nv = dirty_domain.shape - # compute deconvolution operator - rng = np.arange(nu) - k = np.minimum(rng, nu-rng) - self._deconv_u = np.roll(fu[k], -nu//2).reshape((-1, 1)) - rng = np.arange(nv) - k = np.minimum(rng, nv-rng) - self._deconv_v = np.roll(fv[k], -nv//2).reshape((1, -1)) self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): import nifty_gridder self._check_input(x, mode) - nu, nv = self._target.shape res = x.to_global_data() if mode == self.TIMES: - res = hartley(res) - res = np.roll(res, (nu//2, nv//2), axis=(0, 1)) - res = res[:nu, :nv] - res *= self._deconv_u - res *= self._deconv_v + res = self._gconf.grid2dirty(res) else: - res = res*self._deconv_u - res *= self._deconv_v - nu2, nv2 = self._domain.shape - res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant', - constant_values=0) - res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1)) - res = hartley(res) + res = self._gconf.dirty2grid(res) return from_global_data(self._tgt(mode), res) class RadioGridder(LinearOperator): def __init__(self, grid_domain, bl, gconf, idx): - self._domain = DomainTuple.make(UnstructuredDomain((idx.shape[0],))) + self._domain = DomainTuple.make(UnstructuredDomain( + (idx.shape[0],))) self._target = DomainTuple.make(grid_domain) self._bl = bl self._gconf = gconf @@ -109,9 +90,9 @@ class RadioGridder(LinearOperator): import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: - res = nifty_gridder.ms2grid( - self._bl, self._gconf, self._idx, x.to_global_data().reshape((-1,1))) + res = nifty_gridder.vis2grid( + self._bl, self._gconf, self._idx, x.to_global_data()) else: - res = nifty_gridder.grid2ms( + res = nifty_gridder.grid2vis( self._bl, self._gconf, self._idx, x.to_global_data()) return from_global_data(self._tgt(mode), res) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index c4946a16..0e157ee3 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -35,11 +35,13 @@ def _l2error(a, b): @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) def test_gridding(nu, nv, N, eps): - uv = np.random.rand(N, 2) - 0.5 - vis = np.random.randn(N) + 1j*np.random.randn(N) + uvw = np.random.rand(N, 3) - 0.5 + ms = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) + # FIXME temporary! + vis = np.ones((N,))+1j*np.ones((N,)) # Nifty - GM = ift.GridderMaker(ift.RGSpace((nu, nv)), eps=eps, uv=uv) + GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, channel_fact=np.array([1.]), eps=eps) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) Op = GM.getFull() @@ -49,7 +51,7 @@ def test_gridding(nu, nv, N, eps): *[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij') dft = pynu*0. for i in range(N): - dft += (vis[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))).real + dft += (vis[i]*np.exp(2j*np.pi*(x*uvw[i, 0] + y*uvw[i, 1]))).real assert_(_l2error(dft, pynu) < eps) @@ -59,9 +61,8 @@ def test_gridding(nu, nv, N, eps): @pmp('N', [1, 10, 100]) def test_build(nu, nv, N, eps): dom = ift.RGSpace([nu, nv]) - uv = np.random.rand(N, 2) - 0.5 - GM = ift.GridderMaker(dom, eps=eps, uv=uv) - # re-order for performance + uvw = np.random.rand(N, 3) - 0.5 + GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=np.array([1.]), eps=eps) R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 -- GitLab From a9b8d316767a867344f711613b2fcbbacc7992f0 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Mon, 27 May 2019 23:14:37 +0200 Subject: [PATCH 03/21] fixes --- nifty5/library/gridder.py | 5 ++++- test/test_operators/test_nft.py | 4 +--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index c7129498..c25097f6 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -90,9 +90,12 @@ class RadioGridder(LinearOperator): import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: + x = x.to_global_data().reshape((-1,1)) + x = self._bl.ms2vis(x, self._idx) res = nifty_gridder.vis2grid( - self._bl, self._gconf, self._idx, x.to_global_data()) + self._bl, self._gconf, self._idx, x) else: res = nifty_gridder.grid2vis( self._bl, self._gconf, self._idx, x.to_global_data()) + res = self._bl.vis2ms(res, self._idx).reshape((-1,)) return from_global_data(self._tgt(mode), res) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 0e157ee3..9ec56a9f 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -36,9 +36,7 @@ def _l2error(a, b): @pmp('N', [1, 10, 100]) def test_gridding(nu, nv, N, eps): uvw = np.random.rand(N, 3) - 0.5 - ms = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) - # FIXME temporary! - vis = np.ones((N,))+1j*np.ones((N,)) + vis = (np.random.randn(N) + 1j*np.random.randn(N)) # Nifty GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, channel_fact=np.array([1.]), eps=eps) -- GitLab From aea10faa66d5a704d0ffab386c8ab07fc9499bf0 Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Tue, 28 May 2019 10:05:16 +0200 Subject: [PATCH 04/21] Cleanup --- nifty5/library/gridder.py | 24 ++++++++---------------- 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index c25097f6..e5bf958a 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -15,12 +15,9 @@ # # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. -import numpy as np - from ..domain_tuple import DomainTuple from ..domains.rg_space import RGSpace from ..domains.unstructured_domain import UnstructuredDomain -from ..fft import hartley from ..operators.linear_operator import LinearOperator from ..sugar import from_global_data, makeDomain @@ -29,13 +26,12 @@ class GridderMaker(object): def __init__(self, dirty_domain, uvw, channel_fact, eps=2e-13): import nifty_gridder dirty_domain = makeDomain(dirty_domain) - if (len(dirty_domain) != 1 or - not isinstance(dirty_domain[0], RGSpace) or - not len(dirty_domain.shape) == 2): + if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace) + or not len(dirty_domain.shape) == 2): raise ValueError("need dirty_domain with exactly one 2D RGSpace") if channel_fact.ndim != 1: raise ValueError("channel_fact must be a 1D array") - bl = nifty_gridder.Baselines(uvw, channel_fact); + bl = nifty_gridder.Baselines(uvw, channel_fact) nxdirty, nydirty = dirty_domain.shape gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) nu = gconf.Nu() @@ -59,14 +55,12 @@ class GridderMaker(object): class _RestOperator(LinearOperator): def __init__(self, dirty_domain, grid_domain, gconf): - import nifty_gridder self._domain = makeDomain(grid_domain) self._target = makeDomain(dirty_domain) self._gconf = gconf self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): - import nifty_gridder self._check_input(x, mode) res = x.to_global_data() if mode == self.TIMES: @@ -78,8 +72,7 @@ class _RestOperator(LinearOperator): class RadioGridder(LinearOperator): def __init__(self, grid_domain, bl, gconf, idx): - self._domain = DomainTuple.make(UnstructuredDomain( - (idx.shape[0],))) + self._domain = DomainTuple.make(UnstructuredDomain((idx.shape[0],))) self._target = DomainTuple.make(grid_domain) self._bl = bl self._gconf = gconf @@ -90,12 +83,11 @@ class RadioGridder(LinearOperator): import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: - x = x.to_global_data().reshape((-1,1)) + x = x.to_global_data().reshape((-1, 1)) x = self._bl.ms2vis(x, self._idx) - res = nifty_gridder.vis2grid( - self._bl, self._gconf, self._idx, x) + res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x) else: - res = nifty_gridder.grid2vis( - self._bl, self._gconf, self._idx, x.to_global_data()) + res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx, + x.to_global_data()) res = self._bl.vis2ms(res, self._idx).reshape((-1,)) return from_global_data(self._tgt(mode), res) -- GitLab From d05c02bbae512214004173d5dbd121074d3ae589 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 28 May 2019 12:11:06 +0200 Subject: [PATCH 05/21] adjust --- nifty5/library/gridder.py | 8 ++++---- test/test_operators/test_nft.py | 3 ++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index c25097f6..93d7c65a 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -35,7 +35,9 @@ class GridderMaker(object): raise ValueError("need dirty_domain with exactly one 2D RGSpace") if channel_fact.ndim != 1: raise ValueError("channel_fact must be a 1D array") - bl = nifty_gridder.Baselines(uvw, channel_fact); + bl = nifty_gridder.Baselines( + uvw, channel_fact, + np.zeros((uvw.shape[0], channel_fact.shape[0]), dtype=np.bool)) nxdirty, nydirty = dirty_domain.shape gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) nu = gconf.Nu() @@ -59,14 +61,12 @@ class GridderMaker(object): class _RestOperator(LinearOperator): def __init__(self, dirty_domain, grid_domain, gconf): - import nifty_gridder self._domain = makeDomain(grid_domain) self._target = makeDomain(dirty_domain) self._gconf = gconf self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): - import nifty_gridder self._check_input(x, mode) res = x.to_global_data() if mode == self.TIMES: @@ -90,7 +90,7 @@ class RadioGridder(LinearOperator): import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: - x = x.to_global_data().reshape((-1,1)) + x = x.to_global_data().reshape((-1, 1)) x = self._bl.ms2vis(x, self._idx) res = nifty_gridder.vis2grid( self._bl, self._gconf, self._idx, x) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 9ec56a9f..4503e79b 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -39,7 +39,8 @@ def test_gridding(nu, nv, N, eps): vis = (np.random.randn(N) + 1j*np.random.randn(N)) # Nifty - GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, channel_fact=np.array([1.]), eps=eps) + GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, + channel_fact=np.array([1.]), eps=eps) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) Op = GM.getFull() -- GitLab From f22c1e8427b165415018e65dca8c8c0a53c0f3a3 Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Tue, 28 May 2019 14:03:24 +0200 Subject: [PATCH 06/21] Add flags --- nifty5/library/gridder.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index e5bf958a..873abe71 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -23,7 +23,7 @@ from ..sugar import from_global_data, makeDomain class GridderMaker(object): - def __init__(self, dirty_domain, uvw, channel_fact, eps=2e-13): + def __init__(self, dirty_domain, uvw, channel_fact, flags, eps=2e-13): import nifty_gridder dirty_domain = makeDomain(dirty_domain) if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace) @@ -31,17 +31,18 @@ class GridderMaker(object): raise ValueError("need dirty_domain with exactly one 2D RGSpace") if channel_fact.ndim != 1: raise ValueError("channel_fact must be a 1D array") - bl = nifty_gridder.Baselines(uvw, channel_fact) + bl = nifty_gridder.Baselines(uvw, channel_fact, flags) nxdirty, nydirty = dirty_domain.shape gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) nu = gconf.Nu() nv = gconf.Nv() - idx = nifty_gridder.getIndices(bl, gconf) + self._idx = nifty_gridder.getIndices(bl, gconf) + self._bl = bl grid_domain = RGSpace([nu, nv], distances=[1, 1], harmonic=False) self._rest = _RestOperator(dirty_domain, grid_domain, gconf) - self._gridder = RadioGridder(grid_domain, bl, gconf, idx) + self._gridder = RadioGridder(grid_domain, bl, gconf, self._idx) def getGridder(self): return self._gridder @@ -52,6 +53,9 @@ class GridderMaker(object): def getFull(self): return self.getRest() @ self._gridder + def ms2vis(self, x): + return self._bl.ms2vis(x, self._idx) + class _RestOperator(LinearOperator): def __init__(self, dirty_domain, grid_domain, gconf): -- GitLab From 2c9cf99feee8a499d1a07fbd426de706a3827a79 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 28 May 2019 14:22:21 +0200 Subject: [PATCH 07/21] adjust to new gridder interface --- demos/bench_gridder.py | 4 +++- test/test_operators/test_nft.py | 6 ++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 99459b2c..56c70e00 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -27,7 +27,9 @@ for ii in range(10, 26): img = ift.from_global_data(uvspace, img) t0 = time() - GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, channel_fact=np.array([1.])) + GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, + channel_fact=np.array([1.]), + flags=np.zeros((N,1), dtype=np.bool)) vis = ift.from_global_data(visspace, vis) op = GM.getFull().adjoint t1 = time() diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 4503e79b..d525dde7 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -40,7 +40,8 @@ def test_gridding(nu, nv, N, eps): # Nifty GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, - channel_fact=np.array([1.]), eps=eps) + channel_fact=np.array([1.]), eps=eps, + flags=np.zeros((N,1), dtype=np.bool)) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) Op = GM.getFull() @@ -61,7 +62,8 @@ def test_gridding(nu, nv, N, eps): def test_build(nu, nv, N, eps): dom = ift.RGSpace([nu, nv]) uvw = np.random.rand(N, 3) - 0.5 - GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=np.array([1.]), eps=eps) + GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=np.array([1.]), eps=eps, + flags=np.zeros((N,1), dtype=np.bool)) R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 -- GitLab From f90e6228c3d90b93df5fcc8b6497271c3f4ea480 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 28 May 2019 14:23:22 +0200 Subject: [PATCH 08/21] cosmetics --- demos/bench_gridder.py | 2 +- test/test_operators/test_nft.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 56c70e00..8289f3ff 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -29,7 +29,7 @@ for ii in range(10, 26): t0 = time() GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, channel_fact=np.array([1.]), - flags=np.zeros((N,1), dtype=np.bool)) + flags=np.zeros((N, 1), dtype=np.bool)) vis = ift.from_global_data(visspace, vis) op = GM.getFull().adjoint t1 = time() diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index d525dde7..9ae1a906 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -41,7 +41,7 @@ def test_gridding(nu, nv, N, eps): # Nifty GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, channel_fact=np.array([1.]), eps=eps, - flags=np.zeros((N,1), dtype=np.bool)) + flags=np.zeros((N, 1), dtype=np.bool)) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) Op = GM.getFull() @@ -63,7 +63,7 @@ def test_build(nu, nv, N, eps): dom = ift.RGSpace([nu, nv]) uvw = np.random.rand(N, 3) - 0.5 GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=np.array([1.]), eps=eps, - flags=np.zeros((N,1), dtype=np.bool)) + flags=np.zeros((N, 1), dtype=np.bool)) R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 -- GitLab From 2bc1bff11a5b56061aefdf509273b18cee3835bd Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 28 May 2019 14:24:47 +0200 Subject: [PATCH 09/21] use correct nifty_gridder branch --- Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index afe594fb..8cbd63bb 100644 --- a/Dockerfile +++ b/Dockerfile @@ -13,7 +13,7 @@ RUN apt-get update && apt-get install -y \ python3-mpi4py python3-matplotlib \ # more optional NIFTy dependencies && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ - && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \ + && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git@powergrid \ && pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft.git \ && pip3 install jupyter \ && rm -rf /var/lib/apt/lists/* -- GitLab From 8e9c5244419a2b8bd5f0ba2221a9c24f2dd5e62c Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Tue, 28 May 2019 16:53:54 +0200 Subject: [PATCH 10/21] Add test of channel_facts --- test/test_operators/test_nft.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 9ae1a906..23a2f456 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -34,13 +34,14 @@ def _l2error(a, b): @pmp('nu', [12, 128]) @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) -def test_gridding(nu, nv, N, eps): +@pmp('channel_fact', [1, 1.2]) +def test_gridding(nu, nv, N, eps, channel_fact): uvw = np.random.rand(N, 3) - 0.5 vis = (np.random.randn(N) + 1j*np.random.randn(N)) # Nifty GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, - channel_fact=np.array([1.]), eps=eps, + channel_fact=np.array([channel_fact]), eps=eps, flags=np.zeros((N, 1), dtype=np.bool)) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) @@ -51,7 +52,7 @@ def test_gridding(nu, nv, N, eps): *[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij') dft = pynu*0. for i in range(N): - dft += (vis[i]*np.exp(2j*np.pi*(x*uvw[i, 0] + y*uvw[i, 1]))).real + dft += (vis[i]*np.exp(2j*np.pi*(x*uvw[i, 0] + y*uvw[i, 1])*channel_fact)).real assert_(_l2error(dft, pynu) < eps) -- GitLab From 859c79befeda1faba3c5f7ccbd4d4c63ebb9f035 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 28 May 2019 18:04:54 +0200 Subject: [PATCH 11/21] fix gridder --- demos/bench_gridder.py | 4 ++-- nifty5/library/gridder.py | 8 ++++---- test/test_operators/test_nft.py | 11 +++++++---- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 8289f3ff..50398860 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -16,11 +16,11 @@ for ii in range(10, 26): print('N = {}'.format(N)) uvw = np.random.rand(N, 3) - 0.5 - vis = np.random.randn(N) + 1j*np.random.randn(N) + vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) uvspace = ift.RGSpace((nu, nv)) - visspace = ift.UnstructuredDomain(N) + visspace = ift.UnstructuredDomain((N,1)) img = np.random.randn(nu*nv) img = img.reshape((nu, nv)) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index 377150dc..86dafc84 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -77,7 +77,8 @@ class _RestOperator(LinearOperator): class RadioGridder(LinearOperator): def __init__(self, grid_domain, bl, gconf, idx): - self._domain = DomainTuple.make(UnstructuredDomain((idx.shape[0],))) + self._domain = DomainTuple.make( + UnstructuredDomain((bl.Nrows(),bl.Nchannels()))) self._target = DomainTuple.make(grid_domain) self._bl = bl self._gconf = gconf @@ -88,11 +89,10 @@ class RadioGridder(LinearOperator): import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: - x = x.to_global_data().reshape((-1, 1)) - x = self._bl.ms2vis(x, self._idx) + x = self._bl.ms2vis(x.to_global_data(), self._idx) res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x) else: res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx, x.to_global_data()) - res = self._bl.vis2ms(res, self._idx).reshape((-1,)) + res = self._bl.vis2ms(res, self._idx) return from_global_data(self._tgt(mode), res) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 23a2f456..0cca0f33 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -37,7 +37,7 @@ def _l2error(a, b): @pmp('channel_fact', [1, 1.2]) def test_gridding(nu, nv, N, eps, channel_fact): uvw = np.random.rand(N, 3) - 0.5 - vis = (np.random.randn(N) + 1j*np.random.randn(N)) + vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) # Nifty GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, @@ -60,11 +60,14 @@ def test_gridding(nu, nv, N, eps, channel_fact): @pmp('nu', [12, 128]) @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) -def test_build(nu, nv, N, eps): +@pmp('cfact', [np.array([1.]), np.array([0.3, 0.5, 2.3])]) +def test_build(nu, nv, N, eps, cfact): dom = ift.RGSpace([nu, nv]) uvw = np.random.rand(N, 3) - 0.5 - GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=np.array([1.]), eps=eps, - flags=np.zeros((N, 1), dtype=np.bool)) + flags=np.zeros((N, cfact.shape[0]), dtype=np.bool) + flags[0,0]=True + GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=cfact, eps=eps, + flags=flags) R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 -- GitLab From 369354fc5e2929f3cfd8e52661c3dc7d0ad00213 Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Wed, 29 May 2019 13:48:01 +0200 Subject: [PATCH 12/21] Add subtraction mode to adder --- nifty5/operators/adder.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nifty5/operators/adder.py b/nifty5/operators/adder.py index ccee05ba..cefaea03 100644 --- a/nifty5/operators/adder.py +++ b/nifty5/operators/adder.py @@ -28,12 +28,15 @@ class Adder(Operator): field : Field or MultiField The field by which the input is shifted. """ - def __init__(self, field): + def __init__(self, field, neg=False): if not isinstance(field, (Field, MultiField)): raise TypeError self._field = field self._domain = self._target = field.domain + self._neg = bool(neg) def apply(self, x): self._check_input(x) + if self._neg: + return x - self._field return x + self._field -- GitLab From 8a8e675141494e2a943739c2d2aceee81c6caa0d Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Wed, 29 May 2019 14:07:06 +0200 Subject: [PATCH 13/21] Gaussian energy takes inverse covariance as input This is a performance tweak. Since inverse diagonal operators are automatically converted to a pure diagonal operator with the inverse on its diagonal, this commit saves memory. --- demos/getting_started_3.py | 3 ++- nifty5/operators/energy_operators.py | 19 +++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py index a4f6458f..00875a34 100644 --- a/demos/getting_started_3.py +++ b/demos/getting_started_3.py @@ -109,7 +109,8 @@ if __name__ == '__main__': minimizer = ift.NewtonCG(ic_newton) # Set up likelihood and information Hamiltonian - likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response) + likelihood = ift.GaussianEnergy(mean=data, + inverse_covariance=N.inverse)(signal_response) H = ift.StandardHamiltonian(likelihood, ic_sampling) initial_mean = ift.MultiField.full(H.domain, 0.) diff --git a/nifty5/operators/energy_operators.py b/nifty5/operators/energy_operators.py index 17a97d80..aad60f4a 100644 --- a/nifty5/operators/energy_operators.py +++ b/nifty5/operators/energy_operators.py @@ -110,8 +110,8 @@ class GaussianEnergy(EnergyOperator): ---------- mean : Field Mean of the Gaussian. Default is 0. - covariance : LinearOperator - Covariance of the Gaussian. Default is the identity operator. + inverse_covariance : LinearOperator + Inverse covariance of the Gaussian. Default is the identity operator. domain : Domain, DomainTuple, tuple of Domain or MultiDomain Operator domain. By default it is inferred from `mean` or `covariance` if specified @@ -121,28 +121,27 @@ class GaussianEnergy(EnergyOperator): At least one of the arguments has to be provided. """ - def __init__(self, mean=None, covariance=None, domain=None): + def __init__(self, mean=None, inverse_covariance=None, domain=None): if mean is not None and not isinstance(mean, (Field, MultiField)): raise TypeError - if covariance is not None and not isinstance(covariance, - LinearOperator): + if inverse_covariance is not None and not isinstance(inverse_covariance, LinearOperator): raise TypeError self._domain = None if mean is not None: self._checkEquivalence(mean.domain) - if covariance is not None: - self._checkEquivalence(covariance.domain) + if inverse_covariance is not None: + self._checkEquivalence(inverse_covariance.domain) if domain is not None: self._checkEquivalence(domain) if self._domain is None: raise ValueError("no domain given") self._mean = mean - if covariance is None: + if inverse_covariance is None: self._op = SquaredNormOperator(self._domain).scale(0.5) else: - self._op = QuadraticFormOperator(covariance.inverse) - self._icov = None if covariance is None else covariance.inverse + self._op = QuadraticFormOperator(inverse_covariance) + self._icov = None if inverse_covariance is None else inverse_covariance def _checkEquivalence(self, newdom): newdom = makeDomain(newdom) -- GitLab From ef3f25d271c4083c132e54a7466ca7cbe21c3189 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Wed, 29 May 2019 15:19:27 +0200 Subject: [PATCH 14/21] adjust to new interface --- nifty5/library/gridder.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index 86dafc84..adb9e995 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -32,12 +32,12 @@ class GridderMaker(object): raise ValueError("need dirty_domain with exactly one 2D RGSpace") if channel_fact.ndim != 1: raise ValueError("channel_fact must be a 1D array") - bl = nifty_gridder.Baselines(uvw, channel_fact, flags) + bl = nifty_gridder.Baselines(uvw, channel_fact) nxdirty, nydirty = dirty_domain.shape gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) nu = gconf.Nu() nv = gconf.Nv() - self._idx = nifty_gridder.getIndices(bl, gconf) + self._idx = nifty_gridder.getIndices(bl, gconf, flags) self._bl = bl grid_domain = RGSpace([nu, nv], distances=[1, 1], harmonic=False) -- GitLab From 3be21948bf6826ecb5ef7240228dd4eb8653fee0 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 4 Jun 2019 16:56:09 +0200 Subject: [PATCH 15/21] test with channels --- demos/bench_gridder.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 50398860..58f53b87 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -10,17 +10,19 @@ np.random.seed(40) N0s, a0s, b0s, c0s = [], [], [], [] for ii in range(10, 26): - nu = 2048 - nv = 2048 + nu = 1024 + nv = 1024 N = int(2**ii) print('N = {}'.format(N)) - - uvw = np.random.rand(N, 3) - 0.5 - vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) + nchan=16 + nrow=N//nchan + fct = 1.+0.00001*np.arange(nchan) + uvw = np.random.rand(nrow, 3) - 0.5 + vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((nrow,nchan)) uvspace = ift.RGSpace((nu, nv)) - visspace = ift.UnstructuredDomain((N,1)) + visspace = ift.UnstructuredDomain((N//nchan,nchan)) img = np.random.randn(nu*nv) img = img.reshape((nu, nv)) @@ -28,8 +30,8 @@ for ii in range(10, 26): t0 = time() GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, - channel_fact=np.array([1.]), - flags=np.zeros((N, 1), dtype=np.bool)) + channel_fact=fct, + flags=np.zeros((N//nchan, nchan), dtype=np.bool)) vis = ift.from_global_data(visspace, vis) op = GM.getFull().adjoint t1 = time() -- GitLab From 28c6c7b60f46b565dc47718c7545db098da7d297 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 18 Jun 2019 20:14:10 +0200 Subject: [PATCH 16/21] switch to better_params branch --- Dockerfile | 2 +- demos/bench_gridder.py | 6 ++++-- nifty5/library/gridder.py | 10 +++++----- test/test_operators/test_nft.py | 28 ++++++++++++++++++---------- 4 files changed, 28 insertions(+), 18 deletions(-) diff --git a/Dockerfile b/Dockerfile index 8cbd63bb..f1b81e61 100644 --- a/Dockerfile +++ b/Dockerfile @@ -13,7 +13,7 @@ RUN apt-get update && apt-get install -y \ python3-mpi4py python3-matplotlib \ # more optional NIFTy dependencies && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ - && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git@powergrid \ + && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git@better_params \ && pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft.git \ && pip3 install jupyter \ && rm -rf /var/lib/apt/lists/* diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 58f53b87..5e107800 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -10,13 +10,15 @@ np.random.seed(40) N0s, a0s, b0s, c0s = [], [], [], [] for ii in range(10, 26): + fovx = 0.0001 + fovy = 0.0002 nu = 1024 nv = 1024 N = int(2**ii) print('N = {}'.format(N)) nchan=16 nrow=N//nchan - fct = 1.+0.00001*np.arange(nchan) + freq = 1e9+1e6*np.arange(nchan) uvw = np.random.rand(nrow, 3) - 0.5 vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((nrow,nchan)) @@ -30,7 +32,7 @@ for ii in range(10, 26): t0 = time() GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, - channel_fact=fct, + freq=freq, fovx=fovx, fovy=fovy, flags=np.zeros((N//nchan, nchan), dtype=np.bool)) vis = ift.from_global_data(visspace, vis) op = GM.getFull().adjoint diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index adb9e995..4e872571 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -24,17 +24,17 @@ import numpy as np class GridderMaker(object): - def __init__(self, dirty_domain, uvw, channel_fact, flags, eps=2e-13): + def __init__(self, dirty_domain, uvw, freq, fovx, fovy, flags, eps=2e-13): import nifty_gridder dirty_domain = makeDomain(dirty_domain) if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace) or not len(dirty_domain.shape) == 2): raise ValueError("need dirty_domain with exactly one 2D RGSpace") - if channel_fact.ndim != 1: - raise ValueError("channel_fact must be a 1D array") - bl = nifty_gridder.Baselines(uvw, channel_fact) + if freq.ndim != 1: + raise ValueError("freq must be a 1D array") + bl = nifty_gridder.Baselines(uvw, freq) nxdirty, nydirty = dirty_domain.shape - gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) + gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, fovx, fovy) nu = gconf.Nu() nv = gconf.Nv() self._idx = nifty_gridder.getIndices(bl, gconf, flags) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index 0cca0f33..d8b7d166 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -29,19 +29,24 @@ pmp = pytest.mark.parametrize def _l2error(a, b): return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2)) +speedOfLight = 299792458. @pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13]) @pmp('nu', [12, 128]) @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) -@pmp('channel_fact', [1, 1.2]) -def test_gridding(nu, nv, N, eps, channel_fact): - uvw = np.random.rand(N, 3) - 0.5 +@pmp('freq', [1e9]) +def test_gridding(nu, nv, N, eps, freq): + fovx = 0.0001 + fovy = 0.0002 + uvw = (np.random.rand(N, 3) - 0.5) + uvw[:,0] /= fovx*freq/speedOfLight + uvw[:,1] /= fovy*freq/speedOfLight vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) # Nifty GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, - channel_fact=np.array([channel_fact]), eps=eps, + freq=np.array([freq]), eps=eps, fovx=fovx, fovy=fovy, flags=np.zeros((N, 1), dtype=np.bool)) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) @@ -50,9 +55,11 @@ def test_gridding(nu, nv, N, eps, channel_fact): # DFT x, y = np.meshgrid( *[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij') + x *= fovx*freq/speedOfLight + y *= fovy*freq/speedOfLight dft = pynu*0. for i in range(N): - dft += (vis[i]*np.exp(2j*np.pi*(x*uvw[i, 0] + y*uvw[i, 1])*channel_fact)).real + dft += (vis[i]*np.exp(2j*np.pi*(x*uvw[i, 0] + y*uvw[i, 1]))).real assert_(_l2error(dft, pynu) < eps) @@ -60,14 +67,15 @@ def test_gridding(nu, nv, N, eps, channel_fact): @pmp('nu', [12, 128]) @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) -@pmp('cfact', [np.array([1.]), np.array([0.3, 0.5, 2.3])]) -def test_build(nu, nv, N, eps, cfact): +@pmp('freq', [np.array([1e9]), np.array([1e9, 2e9, 2.5e9])]) +def test_build(nu, nv, N, eps, freq): dom = ift.RGSpace([nu, nv]) + fov = np.pi/180/60 uvw = np.random.rand(N, 3) - 0.5 - flags=np.zeros((N, cfact.shape[0]), dtype=np.bool) + flags=np.zeros((N, freq.shape[0]), dtype=np.bool) flags[0,0]=True - GM = ift.GridderMaker(dom, uvw=uvw, channel_fact=cfact, eps=eps, - flags=flags) + GM = ift.GridderMaker(dom, uvw=uvw, freq=freq, eps=eps, + flags=flags, fovx=fov, fovy=fov) R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 -- GitLab From 437c39947de21c782f34164ea28a4499d53bca81 Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Fri, 21 Jun 2019 11:00:39 +0200 Subject: [PATCH 17/21] Add more information to __repr__ in LogRGSpaces --- nifty5/domains/log_rg_space.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nifty5/domains/log_rg_space.py b/nifty5/domains/log_rg_space.py index 65cfb46a..4cd66c03 100644 --- a/nifty5/domains/log_rg_space.py +++ b/nifty5/domains/log_rg_space.py @@ -80,8 +80,8 @@ class LogRGSpace(StructuredDomain): return np.array(self._t_0) def __repr__(self): - return ("LogRGSpace(shape={}, harmonic={})".format( - self.shape, self.harmonic)) + return ("LogRGSpace(shape={}, bindistances={}, t_0={}, harmonic={})".format( + self.shape, self.bindistances, self.t_0, self.harmonic)) def get_default_codomain(self): """Returns a :class:`LogRGSpace` object representing the (position or -- GitLab From 2d9e651174951f964feb16c7df95c117e2dba66e Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Fri, 21 Jun 2019 11:42:33 +0200 Subject: [PATCH 18/21] revert to old gridder interface (almost) --- demos/bench_gridder.py | 16 +++++----------- nifty5/library/gridder.py | 30 ++++++++++++++++++----------- test/test_operators/test_nft.py | 34 +++++++++++---------------------- 3 files changed, 35 insertions(+), 45 deletions(-) diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 5e107800..48ae4f3b 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -10,30 +10,24 @@ np.random.seed(40) N0s, a0s, b0s, c0s = [], [], [], [] for ii in range(10, 26): - fovx = 0.0001 - fovy = 0.0002 nu = 1024 nv = 1024 N = int(2**ii) print('N = {}'.format(N)) - nchan=16 - nrow=N//nchan - freq = 1e9+1e6*np.arange(nchan) - uvw = np.random.rand(nrow, 3) - 0.5 - vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((nrow,nchan)) + + uv = np.random.rand(N, 2) - 0.5 + vis = np.random.randn(N) + 1j*np.random.randn(N) uvspace = ift.RGSpace((nu, nv)) - visspace = ift.UnstructuredDomain((N//nchan,nchan)) + visspace = ift.UnstructuredDomain(N) img = np.random.randn(nu*nv) img = img.reshape((nu, nv)) img = ift.from_global_data(uvspace, img) t0 = time() - GM = ift.GridderMaker(uvspace, eps=1e-7, uvw=uvw, - freq=freq, fovx=fovx, fovy=fovy, - flags=np.zeros((N//nchan, nchan), dtype=np.bool)) + GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv) vis = ift.from_global_data(visspace, vis) op = GM.getFull().adjoint t1 = time() diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index 4e872571..aff1b7c4 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -24,20 +24,28 @@ import numpy as np class GridderMaker(object): - def __init__(self, dirty_domain, uvw, freq, fovx, fovy, flags, eps=2e-13): + def __init__(self, dirty_domain, uv, eps=2e-13): import nifty_gridder dirty_domain = makeDomain(dirty_domain) if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace) or not len(dirty_domain.shape) == 2): raise ValueError("need dirty_domain with exactly one 2D RGSpace") - if freq.ndim != 1: - raise ValueError("freq must be a 1D array") - bl = nifty_gridder.Baselines(uvw, freq) + if uv.ndim != 2: + raise ValueError("uv must be a 2D array") + if uv.shape[1] != 2: + raise ValueError("second dimension of uv must have length 2") + # wasteful hack to adjust to shape required by nifty_gridder + uvw = np.empty((uv.shape[0],3), dtype=np.float64) + uvw[:,0:2] = uv + uvw[:,2] = 0. + speedOfLight = 299792458. + bl = nifty_gridder.Baselines(uvw, np.array([speedOfLight])) nxdirty, nydirty = dirty_domain.shape - gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, fovx, fovy) - nu = gconf.Nu() - nv = gconf.Nv() - self._idx = nifty_gridder.getIndices(bl, gconf, flags) + nxd, nyd = dirty_domain.shape + gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.) + nu, nv = gconf.Nu(), gconf.Nv() + self._idx = nifty_gridder.getIndices( + bl, gconf, np.zeros((uv.shape[0],1),dtype=np.bool)) self._bl = bl grid_domain = RGSpace([nu, nv], distances=[1, 1], harmonic=False) @@ -78,7 +86,7 @@ class _RestOperator(LinearOperator): class RadioGridder(LinearOperator): def __init__(self, grid_domain, bl, gconf, idx): self._domain = DomainTuple.make( - UnstructuredDomain((bl.Nrows(),bl.Nchannels()))) + UnstructuredDomain((bl.Nrows()))) self._target = DomainTuple.make(grid_domain) self._bl = bl self._gconf = gconf @@ -89,10 +97,10 @@ class RadioGridder(LinearOperator): import nifty_gridder self._check_input(x, mode) if mode == self.TIMES: - x = self._bl.ms2vis(x.to_global_data(), self._idx) + x = self._bl.ms2vis(x.to_global_data().reshape((-1, 1)), self._idx) res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x) else: res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx, x.to_global_data()) - res = self._bl.vis2ms(res, self._idx) + res = self._bl.vis2ms(res, self._idx).reshape((-1,)) return from_global_data(self._tgt(mode), res) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index d8b7d166..114675bb 100644 --- a/test/test_operators/test_nft.py +++ b/test/test_operators/test_nft.py @@ -29,37 +29,30 @@ pmp = pytest.mark.parametrize def _l2error(a, b): return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2)) -speedOfLight = 299792458. @pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13]) @pmp('nu', [12, 128]) @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) -@pmp('freq', [1e9]) -def test_gridding(nu, nv, N, eps, freq): - fovx = 0.0001 - fovy = 0.0002 - uvw = (np.random.rand(N, 3) - 0.5) - uvw[:,0] /= fovx*freq/speedOfLight - uvw[:,1] /= fovy*freq/speedOfLight - vis = (np.random.randn(N) + 1j*np.random.randn(N)).reshape((-1,1)) +def test_gridding(nu, nv, N, eps): + uv = np.random.rand(N, 2) - 0.5 + vis = np.random.randn(N) + 1j*np.random.randn(N) # Nifty - GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uvw=uvw, - freq=np.array([freq]), eps=eps, fovx=fovx, fovy=fovy, - flags=np.zeros((N, 1), dtype=np.bool)) + GM = ift.GridderMaker(ift.RGSpace((nu, nv)), uv=uv, eps=eps) vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis) Op = GM.getFull() pynu = Op(vis2).to_global_data() + import matplotlib.pyplot as plt + plt.imshow(pynu) + plt.show() # DFT x, y = np.meshgrid( *[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij') - x *= fovx*freq/speedOfLight - y *= fovy*freq/speedOfLight dft = pynu*0. for i in range(N): - dft += (vis[i]*np.exp(2j*np.pi*(x*uvw[i, 0] + y*uvw[i, 1]))).real + dft += (vis[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))).real assert_(_l2error(dft, pynu) < eps) @@ -67,15 +60,10 @@ def test_gridding(nu, nv, N, eps, freq): @pmp('nu', [12, 128]) @pmp('nv', [4, 12, 128]) @pmp('N', [1, 10, 100]) -@pmp('freq', [np.array([1e9]), np.array([1e9, 2e9, 2.5e9])]) -def test_build(nu, nv, N, eps, freq): +def test_build(nu, nv, N, eps): dom = ift.RGSpace([nu, nv]) - fov = np.pi/180/60 - uvw = np.random.rand(N, 3) - 0.5 - flags=np.zeros((N, freq.shape[0]), dtype=np.bool) - flags[0,0]=True - GM = ift.GridderMaker(dom, uvw=uvw, freq=freq, eps=eps, - flags=flags, fovx=fov, fovy=fov) + uv = np.random.rand(N, 2) - 0.5 + GM = ift.GridderMaker(dom, uv=uv, eps=eps) R0 = GM.getGridder() R1 = GM.getRest() R = R1@R0 -- GitLab From 8b1423c66788474872e83ac9be1de6d06cbbf2fb Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Fri, 21 Jun 2019 13:31:49 +0200 Subject: [PATCH 19/21] Make convention for uv in gridder compatible with NIFTy spaces --- nifty5/library/gridder.py | 4 ++++ test/test_operators/test_nft.py | 11 ++++++----- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index aff1b7c4..e01ef6b3 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -34,10 +34,14 @@ class GridderMaker(object): raise ValueError("uv must be a 2D array") if uv.shape[1] != 2: raise ValueError("second dimension of uv must have length 2") + dstx, dsty = dirty_domain[0].distances # wasteful hack to adjust to shape required by nifty_gridder uvw = np.empty((uv.shape[0],3), dtype=np.float64) uvw[:,0:2] = uv uvw[:,2] = 0. + # Scale uv such that 0 Date: Fri, 21 Jun 2019 14:08:41 +0200 Subject: [PATCH 20/21] typo --- nifty5/library/gridder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index e01ef6b3..17f505f6 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -39,7 +39,7 @@ class GridderMaker(object): uvw = np.empty((uv.shape[0],3), dtype=np.float64) uvw[:,0:2] = uv uvw[:,2] = 0. - # Scale uv such that 0 Date: Mon, 24 Jun 2019 13:15:57 +0200 Subject: [PATCH 21/21] more tests, cleanups --- nifty5/library/gridder.py | 16 +++++++-------- test/test_operators/test_nft.py | 35 ++++++++++++++++++++++++++++++--- 2 files changed, 40 insertions(+), 11 deletions(-) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index 17f505f6..c9d33ed9 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -36,23 +36,23 @@ class GridderMaker(object): raise ValueError("second dimension of uv must have length 2") dstx, dsty = dirty_domain[0].distances # wasteful hack to adjust to shape required by nifty_gridder - uvw = np.empty((uv.shape[0],3), dtype=np.float64) - uvw[:,0:2] = uv - uvw[:,2] = 0. + uvw = np.empty((uv.shape[0], 3), dtype=np.float64) + uvw[:, 0:2] = uv + uvw[:, 2] = 0. # Scale uv such that 0