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
be306136
Commit
be306136
authored
Aug 26, 2019
by
Martin Reinecke
Browse files
remove real-only methods
parent
c32c14c7
Changes
3
Hide whitespace changes
Inline
Side-by-side
gridder_cxx.h
View file @
be306136
...
...
@@ -31,35 +31,32 @@
#include
"pocketfft_hdronly.h"
#include
<array>
template
<
typename
T
,
size_t
ndim
,
bool
c_contiguous
>
class
mav
namespace
gridder
{
using
namespace
std
;
template
<
typename
T
,
size_t
ndim
>
class
mav
{
static_assert
((
ndim
>
0
)
&&
(
ndim
<
3
),
"only supports 1D and 2D arrays"
);
private:
T
*
d
;
std
::
array
<
size_t
,
ndim
>
shp
;
std
::
array
<
ptrdiff_t
,
ndim
>
str
;
array
<
size_t
,
ndim
>
shp
;
array
<
ptrdiff_t
,
ndim
>
str
;
public:
mav
(
T
*
d_
,
const
std
::
array
<
size_t
,
ndim
>
&
shp_
,
const
std
::
array
<
ptrdiff_t
,
ndim
>
&
str_
)
:
d
(
d_
),
shp
(
shp_
),
str
(
str_
)
{
static_assert
(
!
c_contiguous
,
"this type does not accept strides"
);
}
mav
(
T
*
d_
,
const
std
::
array
<
size_t
,
ndim
>
&
shp_
)
mav
(
T
*
d_
,
const
array
<
size_t
,
ndim
>
&
shp_
,
const
array
<
ptrdiff_t
,
ndim
>
&
str_
)
:
d
(
d_
),
shp
(
shp_
),
str
(
str_
)
{}
mav
(
T
*
d_
,
const
array
<
size_t
,
ndim
>
&
shp_
)
:
d
(
d_
),
shp
(
shp_
)
{
static_assert
(
c_contiguous
,
"this type requires strides"
);
str
[
ndim
-
1
]
=
1
;
for
(
size_t
d
=
2
;
d
<=
ndim
;
++
d
)
str
[
ndim
-
d
]
=
str
[
ndim
-
d
+
1
]
*
shp
[
ndim
-
d
+
1
];
}
operator
mav
<
const
T
,
ndim
,
true
>
()
const
{
static_assert
(
c_contiguous
);
return
mav
<
const
T
,
ndim
,
true
>
(
d
,
shp
);
}
operator
mav
<
const
T
,
ndim
,
false
>
()
const
{
return
mav
<
const
T
,
ndim
,
false
>
(
d
,
shp
,
str
);
}
operator
mav
<
const
T
,
ndim
>
()
const
{
return
mav
<
const
T
,
ndim
>
(
d
,
shp
,
str
);
}
T
&
operator
[](
size_t
i
)
{
return
operator
()(
i
);
}
const
T
&
operator
[](
size_t
i
)
const
...
...
@@ -67,25 +64,25 @@ template<typename T, size_t ndim, bool c_contiguous> class mav
T
&
operator
()(
size_t
i
)
{
static_assert
(
ndim
==
1
,
"ndim must be 1"
);
return
c_contiguous
?
d
[
i
]
:
d
[
str
[
0
]
*
i
];
return
d
[
str
[
0
]
*
i
];
}
const
T
&
operator
()(
size_t
i
)
const
{
static_assert
(
ndim
==
1
,
"ndim must be 1"
);
return
c_contiguous
?
d
[
i
]
:
d
[
str
[
0
]
*
i
];
return
d
[
str
[
0
]
*
i
];
}
T
&
operator
()(
size_t
i
,
size_t
j
)
{
static_assert
(
ndim
==
2
,
"ndim must be 2"
);
return
c_contiguous
?
d
[
str
[
0
]
*
i
+
j
]
:
d
[
str
[
0
]
*
i
+
str
[
1
]
*
j
];
return
d
[
str
[
0
]
*
i
+
str
[
1
]
*
j
];
}
const
T
&
operator
()(
size_t
i
,
size_t
j
)
const
{
static_assert
(
ndim
==
2
,
"ndim must be 2"
);
return
c_contiguous
?
d
[
str
[
0
]
*
i
+
j
]
:
d
[
str
[
0
]
*
i
+
str
[
1
]
*
j
];
return
d
[
str
[
0
]
*
i
+
str
[
1
]
*
j
];
}
size_t
shape
(
size_t
i
)
const
{
return
shp
[
i
];
}
const
std
::
array
<
size_t
,
ndim
>
&
shape
()
const
{
return
shp
;
}
const
array
<
size_t
,
ndim
>
&
shape
()
const
{
return
shp
;
}
size_t
size
()
const
{
size_t
res
=
1
;
...
...
@@ -94,20 +91,11 @@ template<typename T, size_t ndim, bool c_contiguous> class mav
}
size_t
stride
(
size_t
i
)
const
{
return
str
[
i
];
}
T
*
data
()
{
// static_assert(c_contiguous, "type is not C contiguous");
return
d
;
}
{
return
d
;
}
const
T
*
data
()
const
{
// static_assert(c_contiguous, "type is not C contiguous");
return
d
;
}
{
return
d
;
}
};
template
<
typename
T
,
size_t
ndim
>
using
cmav
=
mav
<
T
,
ndim
,
true
>
;
template
<
typename
T
,
size_t
ndim
>
using
smav
=
mav
<
T
,
ndim
,
false
>
;
//
// basic utilities
//
...
...
@@ -115,11 +103,11 @@ template<typename T, size_t ndim> using smav = mav<T,ndim,false>;
void
myassert
(
bool
cond
,
const
char
*
msg
)
{
if
(
cond
)
return
;
throw
std
::
runtime_error
(
msg
);
throw
runtime_error
(
msg
);
}
template
<
size_t
ndim
>
void
checkShape
(
const
std
::
array
<
size_t
,
ndim
>
&
shp1
,
const
std
::
array
<
size_t
,
ndim
>
&
shp2
)
(
const
array
<
size_t
,
ndim
>
&
shp1
,
const
array
<
size_t
,
ndim
>
&
shp2
)
{
for
(
size_t
i
=
0
;
i
<
ndim
;
++
i
)
myassert
(
shp1
[
i
]
==
shp2
[
i
],
"shape mismatch"
);
...
...
@@ -138,27 +126,12 @@ template<typename T> inline T fmodulo (T v1, T v2)
static
size_t
nthreads
=
1
;
constexpr
auto
set_nthreads_DS
=
R"""(
Specifies the number of threads to be used by the module
Parameters
==========
nthreads: int
the number of threads to be used. Must be >=1.
)"""
;
void
set_nthreads
(
size_t
nthreads_
)
{
myassert
(
nthreads_
>=
1
,
"nthreads must be >= 1"
);
nthreads
=
nthreads_
;
}
constexpr
auto
get_nthreads_DS
=
R"""(
Returns the number of threads used by the module
Returns
=======
int : the number of threads used by the module
)"""
;
size_t
get_nthreads
()
{
return
nthreads
;
}
...
...
@@ -169,7 +142,7 @@ size_t get_nthreads()
static
inline
double
one_minus_x2
(
double
x
)
{
return
(
fabs
(
x
)
>
0.1
)
?
(
1.
+
x
)
*
(
1.
-
x
)
:
1.
-
x
*
x
;
}
void
legendre_prep
(
int
n
,
std
::
vector
<
double
>
&
x
,
std
::
vector
<
double
>
&
w
)
void
legendre_prep
(
int
n
,
vector
<
double
>
&
x
,
vector
<
double
>
&
w
)
{
constexpr
double
pi
=
3.141592653589793238462643383279502884197
;
constexpr
double
eps
=
3e-14
;
...
...
@@ -229,7 +202,7 @@ void legendre_prep(int n, std::vector<double> &x, std::vector<double> &w)
size_t
get_supp
(
double
epsilon
)
{
static
const
std
::
vector
<
double
>
maxmaperr
{
1e8
,
0.32
,
0.021
,
6.2e-4
,
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
};
...
...
@@ -237,70 +210,7 @@ size_t get_supp(double epsilon)
for
(
size_t
i
=
1
;
i
<
maxmaperr
.
size
();
++
i
)
if
(
epssq
>
maxmaperr
[
i
])
return
i
;
throw
std
::
runtime_error
(
"requested epsilon too small - minimum is 2e-13"
);
}
template
<
typename
T
>
void
complex2hartley
(
const
smav
<
const
std
::
complex
<
T
>
,
2
>
&
grid
,
smav
<
T
,
2
>
&
grid2
)
{
myassert
(
grid
.
shape
()
==
grid2
.
shape
(),
"shape mismatch"
);
size_t
nu
=
grid
.
shape
(
0
),
nv
=
grid
.
shape
(
1
);
#pragma omp parallel for num_threads(nthreads)
for
(
size_t
u
=
0
;
u
<
nu
;
++
u
)
{
size_t
xu
=
(
u
==
0
)
?
0
:
nu
-
u
;
for
(
size_t
v
=
0
;
v
<
nv
;
++
v
)
{
size_t
xv
=
(
v
==
0
)
?
0
:
nv
-
v
;
grid2
(
u
,
v
)
+=
T
(
0.5
)
*
(
grid
(
u
,
v
).
real
()
+
grid
(
u
,
v
).
imag
()
+
grid
(
xu
,
xv
).
real
()
-
grid
(
xu
,
xv
).
imag
());
}
}
}
template
<
typename
T
>
void
hartley2complex
(
const
smav
<
const
T
,
2
>
&
grid
,
smav
<
std
::
complex
<
T
>
,
2
>
&
grid2
)
{
myassert
(
grid
.
shape
()
==
grid2
.
shape
(),
"shape mismatch"
);
size_t
nu
=
grid
.
shape
(
0
),
nv
=
grid
.
shape
(
1
);
#pragma omp parallel for num_threads(nthreads)
for
(
size_t
u
=
0
;
u
<
nu
;
++
u
)
{
size_t
xu
=
(
u
==
0
)
?
0
:
nu
-
u
;
for
(
size_t
v
=
0
;
v
<
nv
;
++
v
)
{
size_t
xv
=
(
v
==
0
)
?
0
:
nv
-
v
;
T
v1
=
T
(
0.5
)
*
grid
(
u
,
v
);
T
v2
=
T
(
0.5
)
*
grid
(
xu
,
xv
);
grid2
(
u
,
v
)
=
std
::
complex
<
T
>
(
v1
+
v2
,
v1
-
v2
);
}
}
}
template
<
typename
T
>
void
hartley2_2D
(
const
smav
<
const
T
,
2
>
&
in
,
smav
<
T
,
2
>
&
out
)
{
myassert
(
in
.
shape
()
==
out
.
shape
(),
"shape mismatch"
);
size_t
nu
=
in
.
shape
(
0
),
nv
=
in
.
shape
(
1
),
sz
=
sizeof
(
T
);
pocketfft
::
stride_t
str
{
ptrdiff_t
(
sz
*
nv
),
ptrdiff_t
(
sz
)};
auto
d_i
=
in
.
data
();
auto
ptmp
=
out
.
data
();
pocketfft
::
r2r_separable_hartley
({
nu
,
nv
},
str
,
str
,
{
0
,
1
},
d_i
,
ptmp
,
T
(
1
),
nthreads
);
#pragma omp parallel for num_threads(nthreads)
for
(
size_t
i
=
1
;
i
<
(
nu
+
1
)
/
2
;
++
i
)
for
(
size_t
j
=
1
;
j
<
(
nv
+
1
)
/
2
;
++
j
)
{
T
a
=
ptmp
[
i
*
nv
+
j
];
T
b
=
ptmp
[(
nu
-
i
)
*
nv
+
j
];
T
c
=
ptmp
[
i
*
nv
+
nv
-
j
];
T
d
=
ptmp
[(
nu
-
i
)
*
nv
+
nv
-
j
];
ptmp
[
i
*
nv
+
j
]
=
T
(
0.5
)
*
(
a
+
b
+
c
-
d
);
ptmp
[(
nu
-
i
)
*
nv
+
j
]
=
T
(
0.5
)
*
(
a
+
b
+
d
-
c
);
ptmp
[
i
*
nv
+
nv
-
j
]
=
T
(
0.5
)
*
(
a
+
c
+
d
-
b
);
ptmp
[(
nu
-
i
)
*
nv
+
nv
-
j
]
=
T
(
0.5
)
*
(
b
+
c
+
d
-
a
);
}
throw
runtime_error
(
"requested epsilon too small - minimum is 2e-13"
);
}
template
<
typename
T
>
class
EC_Kernel
...
...
@@ -322,7 +232,7 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T>
protected:
static
constexpr
T
pi
=
T
(
3.141592653589793238462643383279502884197
L
);
int
p
;
std
::
vector
<
T
>
x
,
wgt
,
psi
;
vector
<
T
>
x
,
wgt
,
psi
;
using
EC_Kernel
<
T
>::
supp
;
public:
...
...
@@ -347,10 +257,10 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T>
};
/* Compute correction factors for the ES gridding kernel
This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
std
::
vector
<
double
>
correction_factors
(
size_t
n
,
size_t
nval
,
size_t
supp
)
vector
<
double
>
correction_factors
(
size_t
n
,
size_t
nval
,
size_t
supp
)
{
EC_Kernel_with_correction
<
double
>
kernel
(
supp
);
std
::
vector
<
double
>
res
(
nval
);
vector
<
double
>
res
(
nval
);
double
xn
=
1.
/
n
;
#pragma omp parallel for schedule(static) num_threads(nthreads)
for
(
size_t
k
=
0
;
k
<
nval
;
++
k
)
...
...
@@ -370,12 +280,12 @@ template<typename T> struct UVW
template
<
typename
T
>
class
Baselines
{
protected:
std
::
vector
<
UVW
<
T
>>
coord
;
std
::
vector
<
T
>
f_over_c
;
vector
<
UVW
<
T
>>
coord
;
vector
<
T
>
f_over_c
;
size_t
nrows
,
nchan
;
public:
Baselines
(
const
s
mav
<
const
T
,
2
>
&
coord_
,
const
s
mav
<
const
T
,
1
>
&
freq
)
Baselines
(
const
mav
<
const
T
,
2
>
&
coord_
,
const
mav
<
const
T
,
1
>
&
freq
)
{
constexpr
double
speedOfLight
=
299792458.
;
myassert
(
coord_
.
shape
(
1
)
==
3
,
"dimension mismatch"
);
...
...
@@ -401,7 +311,7 @@ template<typename T> class Baselines
size_t
Nrows
()
const
{
return
nrows
;
}
size_t
Nchannels
()
const
{
return
nchan
;
}
void
effectiveUVW
(
const
s
mav
<
const
uint32_t
,
1
>
&
idx
,
s
mav
<
T
,
2
>
&
res
)
const
void
effectiveUVW
(
const
mav
<
const
uint32_t
,
1
>
&
idx
,
mav
<
T
,
2
>
&
res
)
const
{
size_t
nvis
=
idx
.
shape
(
0
);
myassert
(
res
.
shape
(
0
)
==
nvis
,
"shape mismatch"
);
...
...
@@ -415,8 +325,8 @@ template<typename T> class Baselines
}
}
template
<
typename
T2
>
void
ms2vis
(
const
s
mav
<
const
T2
,
2
>
&
ms
,
const
s
mav
<
const
uint32_t
,
1
>
&
idx
,
s
mav
<
T2
,
1
>
&
vis
)
const
template
<
typename
T2
>
void
ms2vis
(
const
mav
<
const
T2
,
2
>
&
ms
,
const
mav
<
const
uint32_t
,
1
>
&
idx
,
mav
<
T2
,
1
>
&
vis
)
const
{
myassert
(
ms
.
shape
(
0
)
==
nrows
,
"shape mismatch"
);
myassert
(
ms
.
shape
(
1
)
==
nchan
,
"shape mismatch"
);
...
...
@@ -432,8 +342,8 @@ template<typename T> class Baselines
}
}
template
<
typename
T2
>
void
vis2ms
(
const
s
mav
<
const
T2
,
1
>
&
vis
,
const
s
mav
<
const
uint32_t
,
1
>
&
idx
,
s
mav
<
T2
,
2
>
&
ms
)
const
template
<
typename
T2
>
void
vis2ms
(
const
mav
<
const
T2
,
1
>
&
vis
,
const
mav
<
const
uint32_t
,
1
>
&
idx
,
mav
<
T2
,
2
>
&
ms
)
const
{
size_t
nvis
=
vis
.
shape
(
0
);
myassert
(
idx
.
shape
(
0
)
==
nvis
,
"shape mismatch"
);
...
...
@@ -457,11 +367,10 @@ template<typename T> class GridderConfig
double
eps
,
psx
,
psy
;
size_t
supp
,
nsafe
,
nu
,
nv
;
T
beta
;
std
::
vector
<
T
>
cfu
,
cfv
;
vector
<
T
>
cfu
,
cfv
;
std
::
complex
<
T
>
wscreen
(
double
x
,
double
y
,
double
w
,
bool
adjoint
)
const
complex
<
T
>
wscreen
(
double
x
,
double
y
,
double
w
,
bool
adjoint
)
const
{
using
namespace
std
;
constexpr
double
pi
=
3.141592653589793238462643383279502884197
;
#if 1
double
eps
=
sqrt
(
x
+
y
);
...
...
@@ -482,7 +391,7 @@ template<typename T> class GridderConfig
:
nx_dirty
(
nxdirty
),
ny_dirty
(
nydirty
),
eps
(
epsilon
),
psx
(
pixsize_x
),
psy
(
pixsize_y
),
supp
(
get_supp
(
epsilon
)),
nsafe
((
supp
+
1
)
/
2
),
nu
(
std
::
max
(
2
*
nsafe
,
2
*
nx_dirty
)),
nv
(
std
::
max
(
2
*
nsafe
,
2
*
ny_dirty
)),
nu
(
max
(
2
*
nsafe
,
2
*
nx_dirty
)),
nv
(
max
(
2
*
nsafe
,
2
*
ny_dirty
)),
beta
(
2.3
*
supp
),
cfu
(
nx_dirty
),
cfv
(
ny_dirty
)
{
...
...
@@ -514,78 +423,30 @@ template<typename T> class GridderConfig
size_t
Nsafe
()
const
{
return
nsafe
;
}
T
Beta
()
const
{
return
beta
;
}
void
grid2dirty
(
const
smav
<
const
T
,
2
>
&
grid
,
smav
<
T
,
2
>
&
grid2
)
const
{
checkShape
(
grid
.
shape
(),
{
nu
,
nv
});
std
::
vector
<
T
>
tmpdat
(
nu
*
nv
);
auto
tmp
=
smav
<
T
,
2
>
(
tmpdat
.
data
(),{
nu
,
nv
},{
nv
,
1
});
hartley2_2D
<
T
>
(
grid
,
tmp
);
checkShape
(
grid2
.
shape
(),
{
nx_dirty
,
ny_dirty
});
for
(
size_t
i
=
0
;
i
<
nx_dirty
;
++
i
)
for
(
size_t
j
=
0
;
j
<
ny_dirty
;
++
j
)
{
size_t
i2
=
nu
-
nx_dirty
/
2
+
i
;
if
(
i2
>=
nu
)
i2
-=
nu
;
size_t
j2
=
nv
-
ny_dirty
/
2
+
j
;
if
(
j2
>=
nv
)
j2
-=
nv
;
grid2
(
i
,
j
)
=
tmp
(
i2
,
j2
)
*
cfu
[
i
]
*
cfv
[
j
];
}
}
#if 0
pyarr_c<T> apply_taper(const pyarr_c<T> &img, bool divide) const
{
checkArray(img, "img", {nx_dirty, ny_dirty});
auto pin = img.data();
auto res = makeArray<T>({nx_dirty, ny_dirty});
auto pout = res.mutable_data();
void
apply_taper
(
const
mav
<
const
T
,
2
>
&
img
,
mav
<
T
,
2
>
&
img2
,
bool
divide
)
const
{
py::gil_scoped_release release;
checkShape
(
img
.
shape
(),
{
nx_dirty
,
ny_dirty
});
checkShape
(
img2
.
shape
(),
{
nx_dirty
,
ny_dirty
});
if
(
divide
)
for
(
size_t
i
=
0
;
i
<
nx_dirty
;
++
i
)
for
(
size_t
j
=
0
;
j
<
ny_dirty
;
++
j
)
pout[ny_dirty*i + j] = pin[ny_dirty*i + j]
/(cfu[i]*cfv[j]);
img2
(
i
,
j
)
=
img
(
i
,
j
)
/
(
cfu
[
i
]
*
cfv
[
j
]);
else
for
(
size_t
i
=
0
;
i
<
nx_dirty
;
++
i
)
for
(
size_t
j
=
0
;
j
<
ny_dirty
;
++
j
)
pout[ny_dirty*i + j] = pin[ny_dirty*i + j]*cfu[i]*cfv[j];
}
return res;
}
pyarr_c<complex<T>> grid2dirty_c(const pyarr_c<complex<T>> &grid) const
{
checkArray(grid, "grid", {nu, nv});
auto tmp = makeArray<complex<T>>({nu, nv});
auto ptmp = tmp.mutable_data();
pocketfft::c2c({nu,nv},{grid.strides(0),grid.strides(1)},
{tmp.strides(0), tmp.strides(1)}, {0,1}, pocketfft::BACKWARD,
grid.data(), tmp.mutable_data(), T(1), nthreads);
auto res = makeArray<complex<T>>({nx_dirty, ny_dirty});
auto pout = res.mutable_data();
{
py::gil_scoped_release release;
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
{
size_t i2 = nu-nx_dirty/2+i;
if (i2>=nu) i2-=nu;
size_t j2 = nv-ny_dirty/2+j;
if (j2>=nv) j2-=nv;
pout[ny_dirty*i + j] = ptmp[nv*i2+j2]*cfu[i]*cfv[j];
}
}
return res;
img2
(
i
,
j
)
=
img
(
i
,
j
)
*
cfu
[
i
]
*
cfv
[
j
];
}
pyarr_c<T> dirty2grid(const pyarr_c<T> &dirty) const
{
checkArray(dirty, "dirty", {nx_dirty, ny_dirty});
auto pdirty = dirty.data();
auto tmp = makeArray<T>({nu, nv});
auto ptmp = tmp.mutable_data();
void
grid2dirty_c
(
const
mav
<
const
complex
<
T
>
,
2
>
&
grid
,
mav
<
complex
<
T
>
,
2
>
&
dirty
)
const
{
py::gil_scoped_release release;
for (size_t i=0; i<nu*nv; ++i)
ptmp[i] = 0.;
checkShape
(
grid
.
shape
(),
{
nu
,
nv
});
checkShape
(
dirty
.
shape
(),
{
nx_dirty
,
ny_dirty
});
vector
<
complex
<
T
>>
tmpdat
(
nu
*
nv
);
auto
tmp
=
mav
<
complex
<
T
>
,
2
>
(
tmpdat
.
data
(),{
nu
,
nv
},{
nv
,
1
});
constexpr
auto
sc
=
sizeof
(
complex
<
T
>
);
pocketfft
::
c2c
({
nu
,
nv
},{
grid
.
stride
(
0
)
*
sc
,
grid
.
stride
(
1
)
*
sc
},
{
tmp
.
stride
(
0
)
*
sc
,
tmp
.
stride
(
1
)
*
sc
},
{
0
,
1
},
pocketfft
::
BACKWARD
,
grid
.
data
(),
tmp
.
data
(),
T
(
1
),
nthreads
);
for
(
size_t
i
=
0
;
i
<
nx_dirty
;
++
i
)
for
(
size_t
j
=
0
;
j
<
ny_dirty
;
++
j
)
{
...
...
@@ -593,12 +454,10 @@ template<typename T> class GridderConfig
if
(
i2
>=
nu
)
i2
-=
nu
;
size_t
j2
=
nv
-
ny_dirty
/
2
+
j
;
if
(
j2
>=
nv
)
j2
-=
nv
;
ptmp[nv*i2+j2] = pdirty[ny_dirty*i + j]
*cfu[i]*cfv[j];
dirty
(
i
,
j
)
=
tmp
(
i2
,
j2
)
*
cfu
[
i
]
*
cfv
[
j
];
}
}
hartley2_2D<T>(tmp, tmp);
return tmp;
}
#if 0
pyarr_c<complex<T>> dirty2grid_c(const pyarr_c<complex<T>> &dirty) const
{
checkArray(dirty, "dirty", {nx_dirty, ny_dirty});
...
...
@@ -677,7 +536,7 @@ template<typename T> class GridderConfig
constexpr
int
logsquare
=
4
;
template
<
typename
T
,
typename
T2
=
std
::
complex
<
T
>
>
class
Helper
template
<
typename
T
,
typename
T2
=
complex
<
T
>
>
class
Helper
{
private:
const
GridderConfig
<
T
>
&
gconf
;
...
...
@@ -689,7 +548,7 @@ template<typename T, typename T2=std::complex<T>> class Helper
int
iu0
,
iv0
;
// start index of the current visibility
int
bu0
,
bv0
;
// start index of the current buffer
std
::
vector
<
T2
>
rbuf
,
wbuf
;
vector
<
T2
>
rbuf
,
wbuf
;
void
dump
()
const
{
...
...
@@ -731,7 +590,7 @@ template<typename T, typename T2=std::complex<T>> class Helper
public:
const
T2
*
p0r
;
T2
*
p0w
;
std
::
vector
<
T
>
kernel
;
vector
<
T
>
kernel
;
Helper
(
const
GridderConfig
<
T
>
&
gconf_
,
const
T2
*
grid_r_
,
T2
*
grid_w_
)
:
gconf
(
gconf_
),
nu
(
gconf
.
Nu
()),
nv
(
gconf
.
Nv
()),
nsafe
(
gconf
.
Nsafe
()),
...
...
@@ -775,4 +634,6 @@ template<typename T, typename T2=std::complex<T>> class Helper
}
};
}
// namespace gridder
#endif
nifty_gridder.cc
View file @
be306136
...
...
@@ -25,6 +25,7 @@
#include
"gridder_cxx.h"
using
namespace
std
;
using
namespace
gridder
;
namespace
py
=
pybind11
;
...
...
@@ -39,15 +40,27 @@ auto None = py::none();
//
// Start of real gridder functionality
//
constexpr
auto
set_nthreads_DS
=
R"""(
Specifies the number of threads to be used by the module
Parameters
==========
nthreads: int
the number of threads to be used. Must be >=1.
)"""
;
constexpr
auto
get_nthreads_DS
=
R"""(
Returns the number of threads used by the module
Returns
=======
int : the number of threads used by the module
)"""
;
template
<
typename
T
>
using
pyarr
=
py
::
array_t
<
T
>
;
// The "_c" suffix here stands for "C memory order, contiguous"
template
<
typename
T
>
using
pyarr_c
=
py
::
array_t
<
T
,
py
::
array
::
c_style
|
py
::
array
::
forcecast
>
;
template
<
typename
T
>
pyarr
_c
<
T
>
makeArray
(
const
vector
<
size_t
>
&
shape
)
{
return
pyarr
_c
<
T
>
(
shape
);
}
template
<
typename
T
>
pyarr
<
T
>
makeArray
(
const
vector
<
size_t
>
&
shape
)
{
return
pyarr
<
T
>
(
shape
);
}
void
checkArray
(
const
py
::
array
&
arr
,
const
char
*
aname
,
const
vector
<
size_t
>
&
shape
)
...
...
@@ -94,7 +107,7 @@ template<typename T> pyarr<T> providePotentialArray(const py::object &in,
return
tmp_
;
}
template
<
typename
T
>
pyarr
_c
<
T
>
provideCArray
(
py
::
object
&
in
,
template
<
typename
T
>
pyarr
<
T
>
provideCArray
(
py
::
object
&
in
,
const
vector
<
size_t
>
&
shape
)
{
if
(
in
.
is_none
())
...
...
@@ -106,19 +119,12 @@ template<typename T> pyarr_c<T> provideCArray(py::object &in,
tmp
[
i
]
=
T
(
0
);
return
tmp_
;
}
auto
tmp_
=
in
.
cast
<
pyarr
_c
<
T
>>
();
auto
tmp_
=
in
.
cast
<
pyarr
<
T
>>
();
checkArray
(
tmp_
,
"temporary"
,
shape
);
return
tmp_
;
}
template
<
size_t
ndim
,
typename
T
>
cmav
<
T
,
ndim
>
make_cmav
(
pyarr_c
<
T
>
&
in
)
{
myassert
(
ndim
==
in
.
ndim
(),
"dimension mismatch"
);
array
<
size_t
,
ndim
>
dims
;
for
(
size_t
i
=
0
;
i
<
ndim
;
++
i
)
dims
[
i
]
=
in
.
shape
(
i
);
return
cmav
<
T
,
ndim
>
(
in
.
mutable_data
(),
dims
);
}
template
<
size_t
ndim
,
typename
T
>
smav
<
T
,
ndim
>
make_smav
(
pyarr
<
T
>
&
in
)
template
<
size_t
ndim
,
typename
T
>
mav
<
T
,
ndim
>
make_mav
(
pyarr
<
T
>
&
in
)
{
myassert
(
ndim
==
in
.
ndim
(),
"dimension mismatch"
);
array
<
size_t
,
ndim
>
dims
;
...
...
@@ -129,9 +135,9 @@ template<size_t ndim, typename T> smav<T,ndim> make_smav(pyarr<T> &in)
str
[
i
]
=
in
.
strides
(
i
)
/
sizeof
(
T
);
myassert
(
str
[
i
]
*
ptrdiff_t
(
sizeof
(
T
))
==
in
.
strides
(
i
),
"weird strides"
);
}
return
s
mav
<
T
,
ndim
>
(
in
.
mutable_data
(),
dims
,
str
);
return
mav
<
T
,
ndim
>
(
in
.
mutable_data
(),
dims
,
str
);
}
template
<
size_t
ndim
,
typename
T
>
s
mav
<
T
,
ndim
>
make_
smav
(
pyarr
_c
<
T
>
&
in
)
template
<
size_t
ndim
,
typename
T
>
mav
<
const
T
,
ndim
>
make_
const_mav
(
const
pyarr
<
T
>
&
in
)
{
myassert
(
ndim
==
in
.
ndim
(),
"dimension mismatch"
);
array
<
size_t
,
ndim
>
dims
;
...
...
@@ -142,77 +148,7 @@ template<size_t ndim, typename T> smav<T,ndim> make_smav(pyarr_c<T> &in)
str
[
i
]
=
in
.
strides
(
i
)
/
sizeof
(
T
);
myassert
(
str
[
i
]
*
ptrdiff_t
(
sizeof
(
T
))
==
in
.
strides
(
i
),
"weird strides"
);
}
return
smav
<
T
,
ndim
>
(
in
.
mutable_data
(),
dims
,
str
);
}
template
<
size_t
ndim
,
typename
T
>
cmav
<
const
T
,
ndim
>
make_const_cmav
(
const
pyarr_c
<
T
>
&
in
)
{
myassert
(
ndim
==
in
.
ndim
(),
"dimension mismatch"
);
array
<
size_t
,
ndim
>
dims
;
for
(
size_t
i
=
0
;
i
<
ndim
;
++
i
)
dims
[
i
]
=
in
.
shape
(
i
);
return
cmav
<
const
T
,
ndim
>
(
in
.
data
(),
dims
);
}
template
<
size_t
ndim
,
typename
T
>
smav
<
const
T
,
ndim
>
make_const_smav
(
const
pyarr
<
T
>
&
in
)
{
myassert
(
ndim
==
in
.
ndim
(),
"dimension mismatch"
);
array
<
size_t
,
ndim
>
dims
;
array
<
ptrdiff_t
,
ndim
>
str
;
for
(
size_t
i
=
0
;
i
<
ndim
;
++
i
)
{
dims
[
i
]
=
in
.
shape
(
i
);
str
[
i
]
=
in
.
strides
(
i
)
/
sizeof
(
T
);
myassert
(
str
[
i
]
*
ptrdiff_t
(
sizeof
(
T
))
==
in
.
strides
(
i
),
"weird strides"
);
}
return
smav
<
const
T
,
ndim
>
(
in
.
data
(),
dims
,
str
);
}
template
<
size_t
ndim
,
typename
T
>
smav
<
const
T
,
ndim
>
make_const_smav
(
const
pyarr_c
<
T
>
&
in
)
{
myassert
(
ndim
==
in
.
ndim
(),
"dimension mismatch"
);
array
<
size_t
,
ndim
>
dims
;
array
<
ptrdiff_t
,
ndim
>
str
;
for
(
size_t
i
=
0
;
i
<
ndim
;
++
i
)
{
dims
[
i
]
=
in
.
shape
(
i
);
str
[
i
]
=
in
.
strides
(
i
)
/
sizeof
(
T
);
myassert
(
str
[
i
]
*
ptrdiff_t
(
sizeof
(
T
))
==
in
.
strides
(
i
),
"weird strides"
);
}
return
smav
<
const
T
,
ndim
>
(
in
.
data
(),
dims
,
str
);
}
template
<
typename
T
>
pyarr_c
<
T
>
complex2hartley
(
const
pyarr_c
<
complex
<
T
>>
&
grid_
,
py
::
object
&
grid_in
)
{
auto
grid
=
make_const_smav
<
2
>
(
grid_
);
size_t
nu
=
grid
.
shape
(
0
),
nv
=
grid
.
shape
(
1
);
auto
res
=
provideCArray
<
T
>
(
grid_in
,
{
nu
,
nv
});
auto
grid2
=
make_smav
<
2
>
(
res
);
{
py
::
gil_scoped_release
release
;
complex2hartley
(
grid
,
grid2
);
}
return
res
;
}
template
<
typename
T
>
pyarr_c
<
complex
<
T
>>
hartley2complex
(
const
pyarr_c
<
T
>
&
grid_
)
{
auto
grid
=
make_const_smav
<
2
>
(
grid_
);
size_t
nu
=
grid
.
shape
(
0
),
nv
=
grid
.
shape
(
1
);
auto
res
=
makeArray
<
complex
<
T
>>
({
nu
,
nv
});
auto
grid2
=
make_smav
<
2
>
(
res
);
{