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
Martin Reinecke
ducc
Commits
adaa0f6c
Commit
adaa0f6c
authored
Mar 24, 2020
by
Martin Reinecke
Browse files
better function names; add docstrings
parent
23c22b21
Pipeline
#71331
passed with stages
in 12 minutes and 35 seconds
Changes
2
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
pyinterpol_ng/interpol_ng.h
View file @
adaa0f6c
...
...
@@ -229,7 +229,7 @@ template<typename T> class Interpolator
cube
.
apply
([](
T
&
v
){
v
=
0.
;});
}
void
interpol
x
(
const
mav
<
T
,
2
>
&
ptg
,
mav
<
T
,
1
>
&
res
)
const
void
interpol
(
const
mav
<
T
,
2
>
&
ptg
,
mav
<
T
,
1
>
&
res
)
const
{
MR_assert
(
!
adjoint
,
"cannot be called in adjoint mode"
);
MR_assert
(
ptg
.
shape
(
0
)
==
res
.
shape
(
0
),
"dimension mismatch"
);
...
...
@@ -275,7 +275,7 @@ template<typename T> class Interpolator
});
}
void
deinterpol
x
(
const
mav
<
T
,
2
>
&
ptg
,
const
mav
<
T
,
1
>
&
data
)
void
deinterpol
(
const
mav
<
T
,
2
>
&
ptg
,
const
mav
<
T
,
1
>
&
data
)
{
MR_assert
(
adjoint
,
"can only be called in adjoint mode"
);
MR_assert
(
ptg
.
shape
(
0
)
==
data
.
shape
(
0
),
"dimension mismatch"
);
...
...
@@ -351,7 +351,7 @@ template<typename T> class Interpolator
}
});
}
void
getSlm
x
(
const
Alm
<
complex
<
T
>>
&
blmT
,
Alm
<
complex
<
T
>>
&
slmT
)
void
getSlm
(
const
Alm
<
complex
<
T
>>
&
blmT
,
Alm
<
complex
<
T
>>
&
slmT
)
{
MR_assert
(
adjoint
,
"can only be called in adjoint mode"
);
Alm
<
complex
<
T
>>
a1
(
lmax
,
lmax
),
a2
(
lmax
,
lmax
);
...
...
pyinterpol_ng/pyinterpol_ng.cc
View file @
adaa0f6c
...
...
@@ -19,9 +19,9 @@ template<typename T> class PyInterpolator: public Interpolator<T>
protected:
using
Interpolator
<
T
>::
lmax
;
using
Interpolator
<
T
>::
kmax
;
using
Interpolator
<
T
>::
interpol
x
;
using
Interpolator
<
T
>::
deinterpol
x
;
using
Interpolator
<
T
>::
getSlm
x
;
using
Interpolator
<
T
>::
interpol
;
using
Interpolator
<
T
>::
deinterpol
;
using
Interpolator
<
T
>::
getSlm
;
public:
PyInterpolator
(
const
py
::
array
&
slmT
,
const
py
::
array
&
blmT
,
...
...
@@ -31,28 +31,28 @@ template<typename T> class PyInterpolator: public Interpolator<T>
epsilon
,
nthreads
)
{}
PyInterpolator
(
int64_t
lmax
,
int64_t
kmax
,
double
epsilon
,
int
nthreads
=
0
)
:
Interpolator
<
T
>
(
lmax
,
kmax
,
epsilon
,
nthreads
)
{}
py
::
array
interpol
(
const
py
::
array
&
ptg
)
const
py
::
array
py
interpol
(
const
py
::
array
&
ptg
)
const
{
auto
ptg2
=
to_mav
<
T
,
2
>
(
ptg
);
auto
res
=
make_Pyarr
<
double
>
({
ptg2
.
shape
(
0
)});
auto
res2
=
to_mav
<
double
,
1
>
(
res
,
true
);
interpol
x
(
ptg2
,
res2
);
interpol
(
ptg2
,
res2
);
return
res
;
}
void
deinterpol
(
const
py
::
array
&
ptg
,
const
py
::
array
&
data
)
void
py
deinterpol
(
const
py
::
array
&
ptg
,
const
py
::
array
&
data
)
{
auto
ptg2
=
to_mav
<
T
,
2
>
(
ptg
);
auto
data2
=
to_mav
<
T
,
1
>
(
data
);
deinterpol
x
(
ptg2
,
data2
);
deinterpol
(
ptg2
,
data2
);
}
py
::
array
getSlm
(
const
py
::
array
&
blmT_
)
py
::
array
py
getSlm
(
const
py
::
array
&
blmT_
)
{
auto
res
=
make_Pyarr
<
complex
<
T
>>
({
Alm_Base
::
Num_Alms
(
lmax
,
lmax
)});
Alm
<
complex
<
T
>>
blmT
(
to_mav
<
complex
<
T
>
,
1
>
(
blmT_
,
false
),
lmax
,
kmax
);
auto
slmT_
=
to_mav
<
complex
<
T
>
,
1
>
(
res
,
true
);
Alm
<
complex
<
T
>>
slmT
(
slmT_
,
lmax
,
lmax
);
getSlm
x
(
blmT
,
slmT
);
getSlm
(
blmT
,
slmT
);
return
res
;
}
};
...
...
@@ -71,20 +71,148 @@ template<typename T> py::array pyrotate_alm(const py::array &alm_, int64_t lmax,
}
#endif
constexpr
const
char
*
pyinterpol_ng_DS
=
R"""(
Python interface for total convolution/interpolation library
All arrays containing spherical harmonic coefficients are assumed to have the
following format:
- values for m=0, l going from 0 to lmax
(these values must have an imaginary part of zero)
- values for m=1, l going from 1 to lmax
(these values can be fully complex)
- values for m=2, l going from 2 to lmax
- ...
- values for m=mmax, l going from mmax to lmax
Error conditions are reported by raising exceptions.
)"""
;
constexpr
const
char
*
pyinterpolator_DS
=
R"""(
Class encapsulating the convolution/interpolation functionality
The class can be configured for interpolation or for adjoint interpolation, by
means of two different constructors.
)"""
;
constexpr
const
char
*
initnormal_DS
=
R"""(
Constructor for interpolation mode
Parameters
----------
sky : numpy.ndarray(numpy.complex)
spherical harmonic coefficients of the sky
beam : numpy.ndarray(numpy.complex)
spherical harmonic coefficients of the sky
lmax : int
maximum l in the coefficient arays
mmax : int
maximum m in the sky coefficients
kmax : int
maximum azimuthal moment in the beam coefficients
epsilon : float
desired accuracy for the interpolation; a typical value is 1e-5
nthreads : the number of threads to use for computation
)"""
;
constexpr
const
char
*
initadjoint_DS
=
R"""(
Constructor for adjoint interpolation mode
Parameters
----------
lmax : int
maximum l in the coefficient arays
mmax : int
maximum m in the sky coefficients
kmax : int
maximum azimuthal moment in the beam coefficients
epsilon : float
desired accuracy for the interpolation; a typical value is 1e-5
nthreads : the number of threads to use for computation
)"""
;
constexpr
const
char
*
interpol_DS
=
R"""(
Computes the interpolated values for a given set of angle triplets
Parameters
----------
ptg : numpy.ndarray(numpy.float64) of shape(N,3)
theta, phi and psi angles (in radian) for N pointings
theta must be in the range [0; pi]
phi must be in the range [0; 2pi]
psi should be in the range [-2pi; 2pi]
Returns
-------
numpy.array(numpy.float64) of shape (N,)
the interpolated values
Notes
-----
- Can only be called in "normal" (i.e. not adjoint) mode
- repeated calls to this method are fine, but for good performance the
number of pointings passed per call should be as large as possible.
)"""
;
constexpr
const
char
*
deinterpol_DS
=
R"""(
Takes a set of angle triplets and interpolated values and spreads them onto the
data cube.
Parameters
----------
ptg : numpy.ndarray(numpy.float64) of shape(N,3)
theta, phi and psi angles (in radian) for N pointings
theta must be in the range [0; pi]
phi must be in the range [0; 2pi]
psi should be in the range [-2pi; 2pi]
data : numpy.ndarray(numpy.float64) of shape (N,)
the interpolated values
Notes
-----
- Can only be called in adjoint mode
- repeated calls to this method are fine, but for good performance the
number of pointings passed per call should be as large as possible.
)"""
;
constexpr
const
char
*
getSlm_DS
=
R"""(
Returns a set of sky spherical hamonic coefficients resulting from adjoint
interplation
Parameters
----------
blmT : numpy.array(numpy.complex)
spherical harmonic coefficients of the beam with lmax and kmax defined
in the constructor call
Returns
-------
numpy.array(numpy.complex):
spherical harmonic coefficients of the sky with lmax and mmax defined
in the constructor call
Notes
-----
- Can only be called in adjoint mode
- must be the last call to the object
)"""
;
}
// unnamed namespace
PYBIND11_MODULE
(
pyinterpol_ng
,
m
)
{
using
namespace
pybind11
::
literals
;
py
::
class_
<
PyInterpolator
<
double
>>
(
m
,
"PyInterpolator"
)
m
.
doc
()
=
pyinterpol_ng_DS
;
py
::
class_
<
PyInterpolator
<
double
>>
(
m
,
"PyInterpolator"
,
pyinterpolator_DS
)
.
def
(
py
::
init
<
const
py
::
array
&
,
const
py
::
array
&
,
int64_t
,
int64_t
,
double
,
int
>
(),
"sky"
_a
,
"beam"
_a
,
"lmax"
_a
,
"kmax"
_a
,
"epsilon"
_a
,
"nthreads"
_a
)
.
def
(
py
::
init
<
int64_t
,
int64_t
,
double
,
int
>
(),
initnormal_DS
,
"sky"
_a
,
"beam"
_a
,
"lmax"
_a
,
"kmax"
_a
,
"epsilon"
_a
,
"nthreads"
_a
)
.
def
(
py
::
init
<
int64_t
,
int64_t
,
double
,
int
>
(),
initadjoint_DS
,
"lmax"
_a
,
"kmax"
_a
,
"epsilon"
_a
,
"nthreads"
_a
)
.
def
(
"interpol"
,
&
PyInterpolator
<
double
>::
interpol
,
"ptg"
_a
)
.
def
(
"deinterpol"
,
&
PyInterpolator
<
double
>::
deinterpol
,
"ptg"
_a
,
"data"
_a
)
.
def
(
"getSlm"
,
&
PyInterpolator
<
double
>::
getSlm
,
"blmT"
_a
);
.
def
(
"interpol"
,
&
PyInterpolator
<
double
>::
py
interpol
,
interpol_DS
,
"ptg"
_a
)
.
def
(
"deinterpol"
,
&
PyInterpolator
<
double
>::
py
deinterpol
,
deinterpol_DS
,
"ptg"
_a
,
"data"
_a
)
.
def
(
"getSlm"
,
&
PyInterpolator
<
double
>::
py
getSlm
,
getSlm_DS
,
"blmT"
_a
);
#if 1
m
.
def
(
"rotate_alm"
,
&
pyrotate_alm
<
double
>
,
"alm"
_a
,
"lmax"
_a
,
"psi"
_a
,
"theta"
_a
,
"phi"
_a
);
...
...
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