totalconvolve.cc 9.57 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
/*
 *  This code is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This code is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this code; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

Martin Reinecke's avatar
Martin Reinecke committed
17
18
19
20
21
22
23
/*
 *  Copyright (C) 2020 Max-Planck-Society
 *  Author: Martin Reinecke
 */

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
24
#include "python/totalconvolve.h"
Martin Reinecke's avatar
Martin Reinecke committed
25

Martin Reinecke's avatar
Martin Reinecke committed
26
namespace ducc0 {
27

Martin Reinecke's avatar
Martin Reinecke committed
28
namespace detail_pymodule_totalconvolve {
29

Martin Reinecke's avatar
Martin Reinecke committed
30
31
32
33
34
35
36
37
38
using namespace std;

namespace py = pybind11;

template<typename T> class PyInterpolator: public Interpolator<T>
  {
  protected:
    using Interpolator<T>::lmax;
    using Interpolator<T>::kmax;
39
    using Interpolator<T>::ncomp;
40
41
42
    using Interpolator<T>::interpol;
    using Interpolator<T>::deinterpol;
    using Interpolator<T>::getSlm;
Martin Reinecke's avatar
Martin Reinecke committed
43

44
vector<Alm<complex<T>>> makevec(const py::array &inp, int64_t lmax, int64_t kmax)
45
  {
46
  auto inp2 = to_mav<complex<T>,2>(inp);
47
48
49
50
51
  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;
  }
52
53
54
55
56
57
58
59
60
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
61
  public:
62
    PyInterpolator(const py::array &slm, const py::array &blm,
63
      bool separate, int64_t lmax, int64_t kmax, T epsilon, T ofactor, int nthreads)
64
65
      : Interpolator<T>(makevec(slm, lmax, lmax),
                        makevec(blm, lmax, kmax),
66
67
68
69
70
71
                        separate, epsilon, ofactor, nthreads) {}
    PyInterpolator(int64_t lmax, int64_t kmax, int64_t ncomp_, T epsilon, T ofactor, int nthreads)
      : Interpolator<T>(lmax, kmax, ncomp_, epsilon, ofactor, nthreads) {}

    using Interpolator<T>::support;

72
    py::array pyinterpol(const py::array &ptg) const
Martin Reinecke's avatar
Martin Reinecke committed
73
74
      {
      auto ptg2 = to_mav<T,2>(ptg);
75
76
      auto res = make_Pyarr<T>({ptg2.shape(0),ncomp});
      auto res2 = to_mav<T,2>(res,true);
77
      interpol(ptg2, res2);
78
      return move(res);
Martin Reinecke's avatar
Martin Reinecke committed
79
80
      }

81
    void pydeinterpol(const py::array &ptg, const py::array &data)
Martin Reinecke's avatar
Martin Reinecke committed
82
83
      {
      auto ptg2 = to_mav<T,2>(ptg);
84
      auto data2 = to_mav<T,2>(data);
85
      deinterpol(ptg2, data2);
Martin Reinecke's avatar
Martin Reinecke committed
86
      }
87
    py::array pygetSlm(const py::array &blm_)
Martin Reinecke's avatar
Martin Reinecke committed
88
      {
89
90
      auto blm = makevec(blm_, lmax, kmax);
      auto res = make_Pyarr<complex<T>>({Alm_Base::Num_Alms(lmax, lmax),blm.size()});
91
      vector<Alm<complex<T>>> slm;
92
      makevec_v(res, lmax, lmax, slm);
93
      getSlm(blm, slm);
94
      return move(res);
Martin Reinecke's avatar
Martin Reinecke committed
95
96
97
      }
  };

Martin Reinecke's avatar
Martin Reinecke committed
98
constexpr const char *totalconvolve_DS = R"""(
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
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
----------
126
127
128
129
130
131
sky : numpy.ndarray((nalm_sky, ncomp), dtype=numpy.complex)
    spherical harmonic coefficients of the sky. ncomp can be 1 or 3.
beam : numpy.ndarray((nalm_beam, ncomp), dtype=numpy.complex)
    spherical harmonic coefficients of the beam. ncomp can be 1 or 3
separate : bool
    whether contributions of individual components should be added together.
132
133
134
135
136
137
lmax : int
    maximum l in the coefficient arays
kmax : int
    maximum azimuthal moment in the beam coefficients
epsilon : float
    desired accuracy for the interpolation; a typical value is 1e-5
138
139
140
141
142
ofactor : float
    oversampling factor to be used for the interpolation grids.
    Should be in the range [1.2; 2], a typical value is 1.5
    Increasing this factor makes (adjoint) convolution slower and
    increases memory consumption, but speeds up interpolation/deinterpolation.
143
144
145
146
147
148
149
150
151
152
153
154
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
kmax : int
    maximum azimuthal moment in the beam coefficients
155
156
157
ncomp : int
    the number of components which are going to input to `deinterpol`.
    Can be 1 or 3.
158
159
epsilon : float
    desired accuracy for the interpolation; a typical value is 1e-5
160
161
162
163
164
ofactor : float
    oversampling factor to be used for the interpolation grids.
    Should be in the range [1.2; 2], a typical value is 1.5
    Increasing this factor makes (adjoint) convolution slower and
    increases memory consumption, but speeds up interpolation/deinterpolation.
165
166
167
168
169
170
171
172
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
----------
173
ptg : numpy.ndarray((N, 3), dtype=numpy.float64)
174
175
176
177
178
179
180
    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
-------
181
numpy.array((N, n2), dtype=numpy.float64)
182
    the interpolated values
183
184
    n2 is either 1 (if separate=True was used in the constructor) or the
    second dimension of the input slm and blm arrays (otherwise)
185
186
187
188
189
190
191
192
193
194
195
196
197
198

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
----------
199
ptg : numpy.ndarray((N,3), dtype=numpy.float64)
200
201
202
203
    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]
204
data : numpy.ndarray((N, n2), dtype=numpy.float64)
205
    the interpolated values
206
    n2 must match the `ncomp` value specified in the constructor.
207
208
209
210
211
212
213
214
215
216

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
217
interpolation
218
219
220

Parameters
----------
221
beam : numpy.array(nalm_beam, nbeam), dtype=numpy.complex)
222
223
    spherical harmonic coefficients of the beam with lmax and kmax defined
    in the constructor call
224
    nbeam must match the ncomp specified in the constructor, unless ncomp was 1.
225
226
227

Returns
-------
228
229
numpy.array(nalm_sky, nbeam), dtype=numpy.complex):
    spherical harmonic coefficients of the sky with lmax defined
230
231
232
233
234
235
236
237
    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
238
void add_totalconvolve(py::module &msup)
Martin Reinecke's avatar
Martin Reinecke committed
239
240
  {
  using namespace pybind11::literals;
Martin Reinecke's avatar
Martin Reinecke committed
241
  auto m = msup.def_submodule("totalconvolve");
Martin Reinecke's avatar
Martin Reinecke committed
242

Martin Reinecke's avatar
Martin Reinecke committed
243
  m.doc() = totalconvolve_DS;
244

Martin Reinecke's avatar
Martin Reinecke committed
245
  using inter_d = PyInterpolator<double>;
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
246
  py::class_<inter_d> (m, "Interpolator", py::module_local(), pyinterpolator_DS)
247
248
    .def(py::init<const py::array &, const py::array &, bool, int64_t, int64_t, double, double, int>(),
      initnormal_DS, "sky"_a, "beam"_a, "separate"_a, "lmax"_a, "kmax"_a, "epsilon"_a, "ofactor"_a=1.5,
249
      "nthreads"_a=0)
250
251
    .def(py::init<int64_t, int64_t, int64_t, double, double, int>(), initadjoint_DS,
      "lmax"_a, "kmax"_a, "ncomp"_a, "epsilon"_a, "ofactor"_a=1.5, "nthreads"_a=0)
Martin Reinecke's avatar
Martin Reinecke committed
252
253
254
255
    .def ("interpol", &inter_d::pyinterpol, interpol_DS, "ptg"_a)
    .def ("deinterpol", &inter_d::pydeinterpol, deinterpol_DS, "ptg"_a, "data"_a)
    .def ("getSlm", &inter_d::pygetSlm, getSlm_DS, "beam"_a)
    .def ("support", &inter_d::support);
256
  using inter_f = PyInterpolator<float>;
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
257
  py::class_<inter_f> (m, "Interpolator_f", py::module_local(), pyinterpolator_DS)
258
259
260
261
262
263
264
265
266
    .def(py::init<const py::array &, const py::array &, bool, int64_t, int64_t, float, float, int>(),
      initnormal_DS, "sky"_a, "beam"_a, "separate"_a, "lmax"_a, "kmax"_a, "epsilon"_a, "ofactor"_a=1.5f,
      "nthreads"_a=0)
    .def(py::init<int64_t, int64_t, int64_t, float, float, int>(), initadjoint_DS,
      "lmax"_a, "kmax"_a, "ncomp"_a, "epsilon"_a, "ofactor"_a=1.5f, "nthreads"_a=0)
    .def ("interpol", &inter_f::pyinterpol, interpol_DS, "ptg"_a)
    .def ("deinterpol", &inter_f::pydeinterpol, deinterpol_DS, "ptg"_a, "data"_a)
    .def ("getSlm", &inter_f::pygetSlm, getSlm_DS, "beam"_a)
    .def ("support", &inter_f::support);
267
  m.def("epsilon_guess", &epsilon_guess, "support"_a, "ofactor"_a);
Martin Reinecke's avatar
Martin Reinecke committed
268
  }
269
270
271

}

Martin Reinecke's avatar
Martin Reinecke committed
272
using detail_pymodule_totalconvolve::add_totalconvolve;
273
274

}