pyinterpol_ng.cc 7.23 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/*
 *  Copyright (C) 2020 Max-Planck-Society
 *  Author: Martin Reinecke
 */

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include "interpol_ng.h"

using namespace std;
using namespace mr;

namespace py = pybind11;

namespace {

template<typename T> class PyInterpolator: public Interpolator<T>
  {
  protected:
    using Interpolator<T>::lmax;
    using Interpolator<T>::kmax;
22
    using Interpolator<T>::ncomp;
23
24
25
    using Interpolator<T>::interpol;
    using Interpolator<T>::deinterpol;
    using Interpolator<T>::getSlm;
Martin Reinecke's avatar
Martin Reinecke committed
26

27
vector<Alm<complex<T>>> makevec(const py::array &inp, int64_t lmax, int64_t kmax)
28
  {
29
  auto inp2 = to_mav<complex<T>,2>(inp);
30
31
32
33
34
  vector<Alm<complex<T>>> res;
  for (size_t i=0; i<inp2.shape(1); ++i)
    res.push_back(Alm<complex<T>>(inp2.template subarray<1>({0,i},{inp2.shape(0),0}),lmax, kmax));
  return res;
  }
35
36
37
38
39
40
41
42
43
void makevec_v(py::array &inp, int64_t lmax, int64_t kmax, vector<Alm<complex<T>>> &res)
  {
  auto inp2 = to_mav<complex<T>,2>(inp, true);
  for (size_t i=0; i<inp2.shape(1); ++i)
    {
    auto xtmp = inp2.template subarray<1>({0,i},{inp2.shape(0),0});
    res.emplace_back(xtmp, lmax, kmax);
    }
  }
Martin Reinecke's avatar
Martin Reinecke committed
44
  public:
45
    PyInterpolator(const py::array &slm, const py::array &blm,
46
      bool separate, int64_t lmax, int64_t kmax, double epsilon, int nthreads=0)
47
48
      : Interpolator<T>(makevec(slm, lmax, lmax),
                        makevec(blm, lmax, kmax),
49
50
51
                        separate, epsilon, nthreads) {}
    PyInterpolator(int64_t lmax, int64_t kmax, int64_t ncomp_, double epsilon, int nthreads=0)
      : Interpolator<T>(lmax, kmax, ncomp_, epsilon, nthreads) {}
52
    py::array pyinterpol(const py::array &ptg) const
Martin Reinecke's avatar
Martin Reinecke committed
53
54
      {
      auto ptg2 = to_mav<T,2>(ptg);
55
      auto res = make_Pyarr<double>({ptg2.shape(0),ncomp});
56
      auto res2 = to_mav<double,2>(res,true);
57
      interpol(ptg2, res2);
Martin Reinecke's avatar
Martin Reinecke committed
58
59
60
      return res;
      }

61
    void pydeinterpol(const py::array &ptg, const py::array &data)
Martin Reinecke's avatar
Martin Reinecke committed
62
63
      {
      auto ptg2 = to_mav<T,2>(ptg);
64
      auto data2 = to_mav<T,2>(data);
65
      deinterpol(ptg2, data2);
Martin Reinecke's avatar
Martin Reinecke committed
66
      }
67
    py::array pygetSlm(const py::array &blm_)
Martin Reinecke's avatar
Martin Reinecke committed
68
      {
69
70
      auto blm = makevec(blm_, lmax, kmax);
      auto res = make_Pyarr<complex<T>>({Alm_Base::Num_Alms(lmax, lmax),blm.size()});
71
72
      vector<Alm<complex<double>>> slm;
      makevec_v(res, lmax, lmax, slm);
73
      getSlm(blm, slm);
Martin Reinecke's avatar
Martin Reinecke committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
      return res;
      }
  };

#if 1
template<typename T> py::array pyrotate_alm(const py::array &alm_, int64_t lmax,
  double psi, double theta, double phi)
  {
  auto a1 = to_mav<complex<T>,1>(alm_);
  auto alm = make_Pyarr<complex<T>>({a1.shape(0)});
  auto a2 = to_mav<complex<T>,1>(alm,true);
  for (size_t i=0; i<a1.shape(0); ++i) a2.v(i)=a1(i);
  auto tmp = Alm<complex<T>>(a2,lmax,lmax);
  rotate_alm(tmp, psi, theta, phi);
  return alm;
  }
#endif

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
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
)""";

Martin Reinecke's avatar
Martin Reinecke committed
217
218
219
220
221
222
} // unnamed namespace

PYBIND11_MODULE(pyinterpol_ng, m)
  {
  using namespace pybind11::literals;

223
224
225
  m.doc() = pyinterpol_ng_DS;

  py::class_<PyInterpolator<double>> (m, "PyInterpolator", pyinterpolator_DS)
226
227
    .def(py::init<const py::array &, const py::array &, bool, int64_t, int64_t, double, int>(),
      initnormal_DS, "sky"_a, "beam"_a, "separate"_a, "lmax"_a, "kmax"_a, "epsilon"_a,
228
      "nthreads"_a)
229
230
    .def(py::init<int64_t, int64_t, int64_t, double, int>(), initadjoint_DS,
      "lmax"_a, "kmax"_a, "ncomp"_a, "epsilon"_a, "nthreads"_a)
231
232
233
    .def ("interpol", &PyInterpolator<double>::pyinterpol, interpol_DS, "ptg"_a)
    .def ("deinterpol", &PyInterpolator<double>::pydeinterpol, deinterpol_DS, "ptg"_a, "data"_a)
    .def ("getSlm", &PyInterpolator<double>::pygetSlm, getSlm_DS, "blmT"_a);
Martin Reinecke's avatar
Martin Reinecke committed
234
235
236
237
238
#if 1
  m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
    "phi"_a);
#endif
  }