Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
ift
nifty_gridder
Commits
8e2ad904
Commit
8e2ad904
authored
Oct 01, 2019
by
Martin Reinecke
Browse files
remove redundant code
parent
7cac7cb7
Changes
1
Hide whitespace changes
Inline
Side-by-side
test.py
View file @
8e2ad904
...
...
@@ -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
)
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment