sht.cc 7.91 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
/*
 *  This file is part of pyHealpix.
 *
 *  pyHealpix 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.
 *
 *  pyHealpix 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 pyHealpix; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *  For more information about HEALPix, see http://healpix.sourceforge.net
 */

/*
 *  pyHealpix is being developed at the Max-Planck-Institut fuer Astrophysik
 */

/*
 *  Copyright (C) 2017-2020 Max-Planck-Society
 *  Author: Martin Reinecke
 */

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <vector>
#include <complex>

35
36
37
38
39
40
41
#include "mr_util/sharp/sharp.h"
#include "mr_util/sharp/sharp_geomhelpers.h"
#include "mr_util/sharp/sharp_almhelpers.h"
#include "mr_util/infra/string_utils.h"
#include "mr_util/infra/error_handling.h"
#include "mr_util/math/constants.h"
#include "mr_util/bindings/pybind_utils.h"
42

43
44
namespace mr {

Martin Reinecke's avatar
Martin Reinecke committed
45
namespace detail_pymodule_sht {
46

47
48
49
50
using namespace std;

namespace py = pybind11;

51
using a_d = py::array_t<double>;
52
53
54
55
56
57
58
59
60
61
using a_d_c = py::array_t<double, py::array::c_style | py::array::forcecast>;
using a_c_c = py::array_t<complex<double>,
  py::array::c_style | py::array::forcecast>;

template<typename T> class py_sharpjob
  {
  private:
    unique_ptr<sharp_geom_info> ginfo;
    unique_ptr<sharp_alm_info> ainfo;
    int64_t lmax_, mmax_, npix_;
62
    int nthreads;
63
64

  public:
65
    py_sharpjob () : lmax_(0), mmax_(0), npix_(0), nthreads(1) {}
66
67
68
69
70
71
72

    string repr() const
      {
      return "<sharpjob_d: lmax=" + dataToString(lmax_) +
        ", mmax=" + dataToString(mmax_) + ", npix=", dataToString(npix_) +".>";
      }

73
74
    void set_nthreads(int64_t nthreads_)
      { nthreads = int(nthreads_); }
75
    void set_gauss_geometry(int64_t nrings, int64_t nphi)
76
77
78
79
80
      {
      MR_assert((nrings>0)&&(nphi>0),"bad grid dimensions");
      npix_=nrings*nphi;
      ginfo = sharp_make_gauss_geom_info (nrings, nphi, 0., 1, nphi);
      }
81
    void set_healpix_geometry(int64_t nside)
82
83
84
85
86
      {
      MR_assert(nside>0,"bad Nside value");
      npix_=12*nside*nside;
      ginfo = sharp_make_healpix_geom_info (nside, 1);
      }
87
    void set_fejer1_geometry(int64_t nrings, int64_t nphi)
88
89
90
91
      {
      MR_assert(nrings>0,"bad nrings value");
      MR_assert(nphi>0,"bad nphi value");
      npix_=nrings*nphi;
92
      ginfo = sharp_make_fejer1_geom_info (nrings, nphi, 0., 1, nphi);
93
      }
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    void set_fejer2_geometry(int64_t nrings, int64_t nphi)
      {
      MR_assert(nrings>0,"bad nrings value");
      MR_assert(nphi>0,"bad nphi value");
      npix_=nrings*nphi;
      ginfo = sharp_make_fejer2_geom_info (nrings, nphi, 0., 1, nphi);
      }
    void set_cc_geometry(int64_t nrings, int64_t nphi)
      {
      MR_assert(nrings>0,"bad nrings value");
      MR_assert(nphi>0,"bad nphi value");
      npix_=nrings*nphi;
      ginfo = sharp_make_cc_geom_info (nrings, nphi, 0., 1, nphi);
      }
    void set_dh_geometry(int64_t nrings, int64_t nphi)
109
110
111
112
113
114
      {
      MR_assert(nrings>1,"bad nrings value");
      MR_assert(nphi>0,"bad nphi value");
      npix_=nrings*nphi;
      ginfo = sharp_make_dh_geom_info (nrings, nphi, 0., 1, nphi);
      }
115
116
117
118
119
120
121
    void set_mw_geometry(int64_t nrings, int64_t nphi)
      {
      MR_assert(nrings>0,"bad nrings value");
      MR_assert(nphi>0,"bad nphi value");
      npix_=nrings*nphi;
      ginfo = sharp_make_mw_geom_info (nrings, nphi, 0., 1, nphi);
      }
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
    void set_triangular_alm_info (int64_t lmax, int64_t mmax)
      {
      MR_assert(mmax>=0,"negative mmax");
      MR_assert(mmax<=lmax,"mmax must not be larger than lmax");
      lmax_=lmax; mmax_=mmax;
      ainfo = sharp_make_triangular_alm_info(lmax,mmax,1);
      }

    int64_t n_alm() const
      { return ((mmax_+1)*(mmax_+2))/2 + (mmax_+1)*(lmax_-mmax_); }

    a_d_c alm2map (const a_c_c &alm) const
      {
      MR_assert(npix_>0,"no map geometry specified");
      MR_assert (alm.size()==n_alm(),
        "incorrect size of a_lm array");
      a_d_c map(npix_);
      auto mr=map.mutable_unchecked<1>();
      auto ar=alm.unchecked<1>();
141
      sharp_alm2map(&ar[0], &mr[0], *ginfo, *ainfo, 0, nthreads);
142
143
144
145
146
147
148
149
150
      return map;
      }
    a_c_c alm2map_adjoint (const a_d_c &map) const
      {
      MR_assert(npix_>0,"no map geometry specified");
      MR_assert (map.size()==npix_,"incorrect size of map array");
      a_c_c alm(n_alm());
      auto mr=map.unchecked<1>();
      auto ar=alm.mutable_unchecked<1>();
151
      sharp_map2alm(&ar[0], &mr[0], *ginfo, *ainfo, 0, nthreads);
152
153
154
155
156
157
158
159
160
      return alm;
      }
    a_c_c map2alm (const a_d_c &map) const
      {
      MR_assert(npix_>0,"no map geometry specified");
      MR_assert (map.size()==npix_,"incorrect size of map array");
      a_c_c alm(n_alm());
      auto mr=map.unchecked<1>();
      auto ar=alm.mutable_unchecked<1>();
161
      sharp_map2alm(&ar[0], &mr[0], *ginfo, *ainfo, SHARP_USE_WEIGHTS, nthreads);
162
163
164
165
166
167
168
169
170
171
      return alm;
      }
    a_d_c alm2map_spin (const a_c_c &alm, int64_t spin) const
      {
      MR_assert(npix_>0,"no map geometry specified");
      auto ar=alm.unchecked<2>();
      MR_assert((ar.shape(0)==2)&&(ar.shape(1)==n_alm()),
        "incorrect size of a_lm array");
      a_d_c map(vector<size_t>{2,size_t(npix_)});
      auto mr=map.mutable_unchecked<2>();
172
      sharp_alm2map_spin(spin, &ar(0,0), &ar(1,0), &mr(0,0), &mr(1,0), *ginfo, *ainfo, 0, nthreads);
173
174
175
176
177
178
179
180
181
182
      return map;
      }
    a_c_c map2alm_spin (const a_d_c &map, int64_t spin) const
      {
      MR_assert(npix_>0,"no map geometry specified");
      auto mr=map.unchecked<2>();
      MR_assert ((mr.shape(0)==2)&&(mr.shape(1)==npix_),
        "incorrect size of map array");
      a_c_c alm(vector<size_t>{2,size_t(n_alm())});
      auto ar=alm.mutable_unchecked<2>();
183
      sharp_map2alm_spin(spin, &ar(0,0), &ar(1,0), &mr(0,0), &mr(1,0), *ginfo, *ainfo, SHARP_USE_WEIGHTS, nthreads);
184
185
186
187
      return alm;
      }
  };

Martin Reinecke's avatar
Martin Reinecke committed
188
const char *sht_DS = R"""(
189
190
191
Python interface for some of the libsharp functionality

Error conditions are reported by raising exceptions.
Martin Reinecke's avatar
Martin Reinecke committed
192
)""";
193

Martin Reinecke's avatar
Martin Reinecke committed
194
void add_sht(py::module &msup)
195
196
  {
  using namespace pybind11::literals;
Martin Reinecke's avatar
Martin Reinecke committed
197
  auto m = msup.def_submodule("sht");
Martin Reinecke's avatar
Martin Reinecke committed
198
  m.doc() = sht_DS;
199

Martin Reinecke's avatar
Martin Reinecke committed
200
  py::class_<py_sharpjob<double>> (m, "sharpjob_d", py::module_local())
201
    .def(py::init<>())
202
    .def("set_nthreads", &py_sharpjob<double>::set_nthreads, "nthreads"_a)
203
    .def("set_gauss_geometry", &py_sharpjob<double>::set_gauss_geometry,
204
      "nrings"_a,"nphi"_a)
205
    .def("set_healpix_geometry", &py_sharpjob<double>::set_healpix_geometry,
206
      "nside"_a)
207
208
209
210
211
212
213
    .def("set_fejer1_geometry", &py_sharpjob<double>::set_fejer1_geometry,
      "nrings"_a, "nphi"_a)
    .def("set_fejer2_geometry", &py_sharpjob<double>::set_fejer2_geometry,
      "nrings"_a, "nphi"_a)
    .def("set_cc_geometry", &py_sharpjob<double>::set_cc_geometry,
      "nrings"_a, "nphi"_a)
    .def("set_dh_geometry", &py_sharpjob<double>::set_dh_geometry,
214
      "nrings"_a, "nphi"_a)
215
    .def("set_mw_geometry", &py_sharpjob<double>::set_mw_geometry,
216
      "nrings"_a, "nphi"_a)
217
218
219
220
221
222
223
224
    .def("set_triangular_alm_info",
      &py_sharpjob<double>::set_triangular_alm_info, "lmax"_a, "mmax"_a)
    .def("n_alm", &py_sharpjob<double>::n_alm)
    .def("alm2map", &py_sharpjob<double>::alm2map,"alm"_a)
    .def("alm2map_adjoint", &py_sharpjob<double>::alm2map_adjoint,"map"_a)
    .def("map2alm", &py_sharpjob<double>::map2alm,"map"_a)
    .def("alm2map_spin", &py_sharpjob<double>::alm2map_spin,"alm"_a,"spin"_a)
    .def("map2alm_spin", &py_sharpjob<double>::map2alm_spin,"map"_a,"spin"_a)
225
    .def("__repr__", &py_sharpjob<double>::repr);
226
  }
227
228
229

}

Martin Reinecke's avatar
Martin Reinecke committed
230
using detail_pymodule_sht::add_sht;
231
232
233

}