/*
 *  This file is part of libsharp2.
 *
 *  libsharp2 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.
 *
 *  libsharp2 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 libsharp2; 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
 */

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

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

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

#include "libsharp2/sharp_cxx.h"
#include "libsharp2/sharp_legendre_roots.h"

using namespace std;

namespace py = pybind11;

namespace {

template<typename T> using pyarr = py::array_t<T>;
template<typename T> using pyarr_c
  = py::array_t<T, py::array::c_style | py::array::forcecast>;

using a_d = py::array_t<double>;
using a_c = py::array_t<complex<double>>;
using a_i = py::array_t<int64_t>;
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>;

void myassert(bool cond, const char *msg)
  { if (!cond) throw runtime_error(msg); }

template<typename T> class py_sharpjob
  {
  private:
    sharp_cxxjob<T> job;
    int64_t lmax_, mmax_, npix_;

  public:
    py_sharpjob () : lmax_(0), mmax_(0), npix_(0) {}

    void set_Gauss_geometry(int64_t nrings, int64_t nphi)
      {
      myassert((nrings>0)&&(nphi>0),"bad grid dimensions");
      npix_=nrings*nphi;
      job.set_Gauss_geometry(nrings, nphi);
      }
    void set_Healpix_geometry(int64_t nside)
      {
      myassert(nside>0,"bad Nside value");
      npix_=12*nside*nside;
      job.set_Healpix_geometry(nside);
      }
    void set_triangular_alm_info (int64_t lmax, int64_t mmax)
      {
      myassert(mmax>=0,"negative mmax");
      myassert(mmax<=lmax,"mmax must not be larger than lmax");
      lmax_=lmax; mmax_=mmax;
      job.set_triangular_alm_info (lmax,mmax);
      }

    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
      {
      myassert(npix_>0,"no map geometry specified");
      myassert (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>();
      job.alm2map(&ar[0],&mr[0],false);
      return map;
      }
    a_c_c alm2map_adjoint (const a_d_c &map) const
      {
      myassert(npix_>0,"no map geometry specified");
      myassert (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>();
      job.alm2map_adjoint(&mr[0],&ar[0],false);
      return alm;
      }
    a_c_c map2alm (const a_d_c &map) const
      {
      myassert(npix_>0,"no map geometry specified");
      myassert (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>();
      job.map2alm(&mr[0],&ar[0],false);
      return alm;
      }
    a_d_c alm2map_spin (const a_c_c &alm, int64_t spin) const
      {
      myassert(npix_>0,"no map geometry specified");
      auto ar=alm.unchecked<2>();
      myassert((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>();
      job.alm2map_spin(&ar(0,0),&ar(1,0),&mr(0,0),&mr(1,0),spin,false);
      return map;
      }
    a_c_c map2alm_spin (const a_d_c &map, int64_t spin) const
      {
      myassert(npix_>0,"no map geometry specified");
      auto mr=map.unchecked<2>();
      myassert ((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>();
      job.map2alm_spin(&mr(0,0),&mr(1,0),&ar(0,0),&ar(1,0),spin,false);
      return alm;
      }
  };

a_d_c GL_weights(int64_t nlat, int64_t nlon)
  {
  constexpr double twopi=6.283185307179586476925286766559005768394;
  a_d_c res(nlat);
  auto rr=res.mutable_unchecked<1>();
  vector<double> dummy_roots(nlat);
  sharp_legendre_roots(nlat, dummy_roots.data(), &rr[0]);
  for (size_t i=0; i<size_t(rr.shape(0)); ++i)
    rr[i]*=twopi/nlon;
  return res;
  }

a_d_c GL_thetas(int64_t nlat)
  {
  a_d_c res(nlat);
  auto rr=res.mutable_unchecked<1>();
  vector<double> dummy_weights(nlat);
  sharp_legendre_roots(nlat, &rr[0], dummy_weights.data());
  for (size_t i=0; i<size_t(rr.shape(0)); ++i)
    rr[i]=acos(-rr[i]);
  return res;
  }


const char *pysharp_DS = R"""(
Python interface for libsharp2

All angles are interpreted as radians.
The theta coordinate is measured as co-latitude, ranging from 0 (North Pole)
to pi (South Pole).

Error conditions are reported by raising exceptions.
)""";


} // unnamed namespace

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

  m.doc() = pysharp_DS;

  py::class_<py_sharpjob<double>> (m, "sharpjob_d")
    .def(py::init<>())
    .def("set_Gauss_geometry", &py_sharpjob<double>::set_Gauss_geometry,
      "nrings"_a,"nphi"_a)
    .def("set_Healpix_geometry", &py_sharpjob<double>::set_Healpix_geometry,
      "nside"_a)
    .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)
    ;

  m.def("GL_weights",&GL_weights, "nlat"_a, "nlon"_a);
  m.def("GL_thetas",&GL_thetas, "nlat"_a);
  }