Commit 7189cb31 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

tweak Python helpers; add upsampling to libsharp

parent 6a74493d
......@@ -129,12 +129,6 @@ std::unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (size_t nrings, size
std::unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t nphi, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat);
/*! Old name for sharp_make_fejer1_geom_info()
\ingroup geominfogroup */
static inline std::unique_ptr<sharp_geom_info> sharp_make_ecp_geom_info (size_t nrings, size_t nphi, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{ return sharp_make_fejer1_geom_info (nrings, nphi, phi0, stride_lon, stride_lat); }
/*! Creates a geometry information describing an ECP map with \a nrings
iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
pixel in each ring is \a phi0 (in radians). The index difference between
......
#ifndef MRUTIL_PYBIND_UTILS_H
#define MRUTIL_PYBIND_UTILS_H
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include "mr_util/mav.h"
namespace mr {
namespace detail_pybind {
namespace py = pybind11;
template<typename T> bool isPyarr(const py::object &obj)
{ return py::isinstance<py::array_t<T>>(obj); }
template<typename T> py::array_t<T> toPyarr(const py::object &obj)
{
auto tmp = obj.cast<py::array_t<T>>();
MR_assert(tmp.is(obj), "error during array conversion");
return tmp;
}
shape_t copy_shape(const py::array &arr)
{
shape_t res(size_t(arr.ndim()));
for (size_t i=0; i<res.size(); ++i)
res[i] = size_t(arr.shape(int(i)));
return res;
}
template<typename T> stride_t copy_strides(const py::array &arr)
{
stride_t res(size_t(arr.ndim()));
constexpr auto st = ptrdiff_t(sizeof(T));
for (size_t i=0; i<res.size(); ++i)
{
auto tmp = arr.strides(int(i));
MR_assert((tmp/st)*st==tmp, "bad stride");
res[i] = tmp/st;
}
return res;
}
template<size_t ndim> std::array<size_t, ndim> copy_fixshape(const py::array &arr)
{
MR_assert(size_t(arr.ndim())==ndim, "incorrect number of dimensions");
std::array<size_t, ndim> res;
for (size_t i=0; i<ndim; ++i)
res[i] = size_t(arr.shape(int(i)));
return res;
}
template<typename T, size_t ndim> std::array<ptrdiff_t, ndim> copy_fixstrides(const py::array &arr)
{
MR_assert(size_t(arr.ndim())==ndim, "incorrect number of dimensions");
std::array<ptrdiff_t, ndim> res;
constexpr auto st = ptrdiff_t(sizeof(T));
for (size_t i=0; i<ndim; ++i)
{
auto tmp = arr.strides(int(i));
MR_assert((tmp/st)*st==tmp, "bad stride");
res[i] = tmp/st;
}
return res;
}
template<typename T> py::array_t<T> make_Pyarr(const shape_t &dims)
{ return py::array_t<T>(dims); }
template<typename T> py::array_t<T> get_optional_Pyarr(py::object &arr_,
const shape_t &dims)
{
if (arr_.is_none()) return py::array_t<T>(dims);
MR_assert(isPyarr<T>(arr_), "incorrect data type");
auto tmp = toPyarr<T>(arr_);
MR_assert(dims.size()==size_t(tmp.ndim()), "dimension mismatch");
for (size_t i=0; i<dims.size(); ++i)
MR_assert(dims[i]==size_t(tmp.shape(int(i))), "dimension mismatch");
return tmp;
}
template<typename T> py::array_t<T> get_optional_const_Pyarr(
const py::object &arr_, const shape_t &dims)
{
if (arr_.is_none()) return py::array_t<T>(shape_t(dims.size(), 0));
MR_assert(isPyarr<T>(arr_), "incorrect data type");
auto tmp = toPyarr<T>(arr_);
MR_assert(dims.size()==size_t(tmp.ndim()), "dimension mismatch");
for (size_t i=0; i<dims.size(); ++i)
MR_assert(dims[i]==size_t(tmp.shape(int(i))), "dimension mismatch");
return tmp;
}
template<typename T> fmav<T> to_fmav(py::object &obj)
{
auto arr = toPyarr<T>(obj);
return fmav<T>(reinterpret_cast<T *>(arr.mutable_data()),
copy_shape(arr), copy_strides<T>(arr));
}
template<typename T> cfmav<T> to_cfmav(const py::object &obj)
{
auto arr = toPyarr<T>(obj);
return cfmav<T>(reinterpret_cast<const T *>(arr.data()),
copy_shape(arr), copy_strides<T>(arr));
}
template<typename T, size_t ndim> mav<T,ndim> to_mav(py::array &obj)
{
auto arr = toPyarr<T>(obj);
return mav<T,ndim>(reinterpret_cast<T *>(arr.mutable_data()),
copy_fixshape<ndim>(arr), copy_fixstrides<T,ndim>(arr));
}
template<typename T, size_t ndim> cmav<T,ndim> to_cmav(const py::array &obj)
{
auto arr = toPyarr<T>(obj);
return cmav<T,ndim>(reinterpret_cast<const T *>(arr.data()),
copy_fixshape<ndim>(arr), copy_fixstrides<T,ndim>(arr));
}
}
using detail_pybind::isPyarr;
using detail_pybind::make_Pyarr;
using detail_pybind::get_optional_Pyarr;
using detail_pybind::get_optional_const_Pyarr;
using detail_pybind::to_fmav;
using detail_pybind::to_cfmav;
using detail_pybind::to_mav;
using detail_pybind::to_cmav;
}
#endif
......@@ -12,4 +12,5 @@ include mr_util/simd.h
include mr_util/cmplx.h
include mr_util/unity_roots.h
include mr_util/error_handling.h
include mr_util/pybind_utils.h
include mr_util/threading.cc
......@@ -16,120 +16,54 @@
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* Copyright (C) 2019 Max-Planck-Society
/* Copyright (C) 2019-2020 Max-Planck-Society
Author: Martin Reinecke */
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include "mr_util/pybind_utils.h"
#include "gridder_cxx.h"
using namespace std;
using namespace gridder;
using namespace mr;
using gridder::detail::idx_t;
namespace py = pybind11;
namespace {
auto None = py::none();
template<typename T>
using pyarr = py::array_t<T, 0>;
template<typename T> bool isPytype(const py::array &arr)
{
auto t1=arr.dtype();
auto t2=pybind11::dtype::of<T>();
auto k1=t1.kind();
auto k2=t2.kind();
auto s1=t1.itemsize();
auto s2=t2.itemsize();
return (k1==k2)&&(s1==s2);
}
template<typename T> pyarr<T> getPyarr(const py::array &arr, const string &name)
{
auto t1=arr.dtype();
auto t2=pybind11::dtype::of<T>();
auto k1=t1.kind();
auto k2=t2.kind();
auto s1=t1.itemsize();
auto s2=t2.itemsize();
MR_assert((k1==k2)&&(s1==s2),
"type mismatch for array '", name, "': expected '", k2, s2,
"', but got '", k1, s1, "'.");
return arr.cast<pyarr<T>>();
}
template<typename T> pyarr<T> makeArray(const vector<size_t> &shape)
{ return pyarr<T>(shape); }
template<typename T> pyarr<T> providePotentialArray(const py::object &in,
const string &name, const vector<size_t> &shape)
{
if (in.is_none())
return makeArray<T>(vector<size_t>(shape.size(), 0));
return getPyarr<T>(in.cast<py::array>(), name);
}
template<size_t ndim, typename T> mav<T,ndim> make_mav(pyarr<T> &in)
{
MR_assert(ndim==in.ndim(), "dimension mismatch");
array<size_t,ndim> dims;
array<ptrdiff_t,ndim> str;
for (size_t i=0; i<ndim; ++i)
{
dims[i]=in.shape(i);
str[i]=in.strides(i)/sizeof(T);
MR_assert(str[i]*ptrdiff_t(sizeof(T))==in.strides(i), "weird strides");
}
return mav<T, ndim>(in.mutable_data(),dims,str);
}
template<size_t ndim, typename T> cmav<T,ndim> make_cmav(const pyarr<T> &in)
{
MR_assert(ndim==in.ndim(), "dimension mismatch");
array<size_t,ndim> dims;
array<ptrdiff_t,ndim> str;
for (size_t i=0; i<ndim; ++i)
{
dims[i]=in.shape(i);
str[i]=in.strides(i)/sizeof(T);
MR_assert(str[i]*ptrdiff_t(sizeof(T))==in.strides(i), "weird strides");
}
return cmav<T, ndim>(in.data(),dims,str);
}
template<typename T> py::array ms2dirty_general2(const py::array &uvw_,
const py::array &freq_, const py::array &ms_, const py::object &wgt_,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu,
size_t nv, double epsilon, bool do_wstacking, size_t nthreads,
size_t verbosity)
{
auto uvw = getPyarr<double>(uvw_, "uvw");
auto uvw2 = make_cmav<2>(uvw);
auto freq = getPyarr<double>(freq_, "freq");
auto freq2 = make_cmav<1>(freq);
auto ms = getPyarr<complex<T>>(ms_, "ms");
auto ms2 = make_cmav<2>(ms);
auto wgt = providePotentialArray<T>(wgt_, "wgt", {ms2.shape(0),ms2.shape(1)});
auto wgt2 = make_cmav<2>(wgt);
auto dirty = makeArray<T>({npix_x,npix_y});
auto dirty2 = make_mav<2>(dirty);
auto uvw = to_cmav<double,2>(uvw_);
auto freq = to_cmav<double,1>(freq_);
auto ms = to_cmav<complex<T>,2>(ms_);
auto wgt = get_optional_const_Pyarr<T>(wgt_, {ms.shape(0),ms.shape(1)});
auto wgt2 = to_cmav<T,2>(wgt);
auto dirty = make_Pyarr<T>({npix_x,npix_y});
auto dirty2 = to_mav<T,2>(dirty);
{
py::gil_scoped_release release;
ms2dirty_general(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y, nu, nv, epsilon,do_wstacking,
nthreads,dirty2,verbosity);
ms2dirty_general(uvw,freq,ms,wgt2,pixsize_x,pixsize_y,nu,nv,epsilon,
do_wstacking,nthreads,dirty2,verbosity);
}
return move(dirty);
}
py::array Pyms2dirty_general(const py::array &uvw,
const py::array &freq, const py::array &ms, const py::object &wgt,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu,
size_t nv, double epsilon, bool do_wstacking, size_t nthreads,
size_t verbosity)
{
if (isPytype<complex<float>>(ms))
if (isPyarr<complex<float>>(ms))
return ms2dirty_general2<float>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<complex<double>>(ms))
if (isPyarr<complex<double>>(ms))
return ms2dirty_general2<double>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
MR_fail("type matching failed: 'ms' has neither type 'c8' nor 'c16'");
......@@ -186,20 +120,17 @@ template<typename T> py::array dirty2ms_general2(const py::array &uvw_,
double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
auto uvw = getPyarr<double>(uvw_, "uvw");
auto uvw2 = make_cmav<2>(uvw);
auto freq = getPyarr<double>(freq_, "freq");
auto freq2 = make_cmav<1>(freq);
auto dirty = getPyarr<T>(dirty_, "dirty");
auto dirty2 = make_cmav<2>(dirty);
auto wgt = providePotentialArray<T>(wgt_, "wgt", {uvw2.shape(0),freq2.shape(0)});
auto wgt2 = make_cmav<2>(wgt);
auto ms = makeArray<complex<T>>({uvw2.shape(0),freq2.shape(0)});
auto ms2 = make_mav<2>(ms);
auto uvw = to_cmav<double,2>(uvw_);
auto freq = to_cmav<double,1>(freq_);
auto dirty = to_cmav<T,2>(dirty_);
auto wgt = get_optional_const_Pyarr<T>(wgt_, {uvw.shape(0),freq.shape(0)});
auto wgt2 = to_cmav<T,2>(wgt);
auto ms = make_Pyarr<complex<T>>({uvw.shape(0),freq.shape(0)});
auto ms2 = to_mav<complex<T>,2>(ms);
{
py::gil_scoped_release release;
dirty2ms_general(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,nu, nv,epsilon,do_wstacking,
nthreads,ms2,verbosity);
dirty2ms_general(uvw,freq,dirty,wgt2,pixsize_x,pixsize_y,nu,nv,epsilon,
do_wstacking,nthreads,ms2,verbosity);
}
return move(ms);
}
......@@ -208,10 +139,10 @@ py::array Pydirty2ms_general(const py::array &uvw,
double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<float>(dirty))
if (isPyarr<float>(dirty))
return dirty2ms_general2<float>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<double>(dirty))
if (isPyarr<double>(dirty))
return dirty2ms_general2<double>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
MR_fail("type matching failed: 'dirty' has neither type 'f4' nor 'f8'");
......
......@@ -46,6 +46,7 @@ def get_extension_modules():
'mr_util/mav.h',
'mr_util/cmplx.h',
'mr_util/unity_roots.h',
'mr_util/pybind_utils.h',
'gridder_cxx.h',
'setup.py'],
include_dirs=include_dirs,
......
......@@ -9,4 +9,5 @@ include mr_util/simd.h
include mr_util/cmplx.h
include mr_util/unity_roots.h
include mr_util/error_handling.h
include mr_util/pybind_utils.h
include mr_util/threading.cc
......@@ -6,7 +6,7 @@
/*
* Python interface.
*
* Copyright (C) 2019 Max-Planck-Society
* Copyright (C) 2019-2020 Max-Planck-Society
* Copyright (C) 2019 Peter Bell
* \author Martin Reinecke
* \author Peter Bell
......@@ -17,6 +17,7 @@
#include <pybind11/stl.h>
#include "mr_util/fft.h"
#include "mr_util/pybind_utils.h"
namespace {
......@@ -24,6 +25,9 @@ using mr::shape_t;
using mr::stride_t;
using mr::fmav;
using mr::cfmav;
using mr::to_fmav;
using mr::to_cfmav;
using mr::get_optional_Pyarr;
using std::size_t;
using std::ptrdiff_t;
......@@ -41,27 +45,6 @@ using f64 = double;
using flong = ldbl_t;
auto None = py::none();
shape_t copy_shape(const py::array &arr)
{
shape_t res(size_t(arr.ndim()));
for (size_t i=0; i<res.size(); ++i)
res[i] = size_t(arr.shape(int(i)));
return res;
}
template<typename T> stride_t copy_strides(const py::array &arr)
{
stride_t res(size_t(arr.ndim()));
constexpr auto st = ptrdiff_t(sizeof(T));
for (size_t i=0; i<res.size(); ++i)
{
auto tmp = arr.strides(int(i));
MR_assert((tmp/st)*st==tmp, "bad stride");
res[i] = tmp/st;
}
return res;
}
shape_t makeaxes(const py::array &in, const py::object &axes)
{
if (axes.is_none())
......@@ -111,34 +94,13 @@ template<typename T> T norm_fct(int inorm, const shape_t &shape,
return norm_fct<T>(inorm, N);
}
template<typename T> py::array_t<T> prepare_output(py::object &out_,
const shape_t &dims)
{
if (out_.is_none()) return py::array_t<T>(dims);
auto tmp = out_.cast<py::array_t<T>>();
if (!tmp.is(out_)) // a new object was created during casting
throw std::runtime_error("unexpected data type for output array");
return tmp;
}
template<typename T> fmav<T> to_fmav(py::array &arr)
{
return fmav<T>(reinterpret_cast<T *>(arr.mutable_data()),
copy_shape(arr), copy_strides<T>(arr));
}
template<typename T> cfmav<T> to_cfmav(const py::array &arr)
{
return cfmav<T>(reinterpret_cast<const T *>(arr.data()),
copy_shape(arr), copy_strides<T>(arr));
}
template<typename T> py::array c2c_internal(const py::array &in,
const py::object &axes_, bool forward, int inorm, py::object &out_,
size_t nthreads)
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<std::complex<T>>(in);
auto out = prepare_output<std::complex<T>>(out_, ain.shape());
auto out = get_optional_Pyarr<std::complex<T>>(out_, ain.shape());
auto aout = to_fmav<std::complex<T>>(out);
{
py::gil_scoped_release release;
......@@ -154,7 +116,7 @@ template<typename T> py::array c2c_sym_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<T>(in);
auto out = prepare_output<std::complex<T>>(out_, ain.shape());
auto out = get_optional_Pyarr<std::complex<T>>(out_, ain.shape());
auto aout = to_fmav<std::complex<T>>(out);
{
py::gil_scoped_release release;
......@@ -192,7 +154,7 @@ template<typename T> py::array r2c_internal(const py::array &in,
auto ain = to_cfmav<T>(in);
auto dims_out(ain.shape());
dims_out[axes.back()] = (dims_out[axes.back()]>>1)+1;
auto out = prepare_output<std::complex<T>>(out_, dims_out);
auto out = get_optional_Pyarr<std::complex<T>>(out_, dims_out);
auto aout = to_fmav<std::complex<T>>(out);
{
py::gil_scoped_release release;
......@@ -215,7 +177,7 @@ template<typename T> py::array r2r_fftpack_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<T>(in);
auto out = prepare_output<T>(out_, ain.shape());
auto out = get_optional_Pyarr<T>(out_, ain.shape());
auto aout = to_fmav<T>(out);
{
py::gil_scoped_release release;
......@@ -239,7 +201,7 @@ template<typename T> py::array dct_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<T>(in);
auto out = prepare_output<T>(out_, ain.shape());
auto out = get_optional_Pyarr<T>(out_, ain.shape());
auto aout = to_fmav<T>(out);
{
py::gil_scoped_release release;
......@@ -265,7 +227,7 @@ template<typename T> py::array dst_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<T>(in);
auto out = prepare_output<T>(out_, ain.shape());
auto out = get_optional_Pyarr<T>(out_, ain.shape());
auto aout = to_fmav<T>(out);
{
py::gil_scoped_release release;
......@@ -297,7 +259,7 @@ template<typename T> py::array c2r_internal(const py::array &in,
if ((lastsize/2) + 1 != ain.shape(axis))
throw std::invalid_argument("bad lastsize");
dims_out[axis] = lastsize;
auto out = prepare_output<T>(out_, dims_out);
auto out = get_optional_Pyarr<T>(out_, dims_out);
auto aout = to_fmav<T>(out);
{
py::gil_scoped_release release;
......@@ -319,7 +281,7 @@ template<typename T> py::array separable_hartley_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<T>(in);
auto out = prepare_output<T>(out_, ain.shape());
auto out = get_optional_Pyarr<T>(out_, ain.shape());
auto aout = to_fmav<T>(out);
{
py::gil_scoped_release release;
......@@ -341,7 +303,7 @@ template<typename T> py::array genuine_hartley_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);
auto ain = to_cfmav<T>(in);
auto out = prepare_output<T>(out_, ain.shape());
auto out = get_optional_Pyarr<T>(out_, ain.shape());
auto aout = to_fmav<T>(out);
{
py::gil_scoped_release release;
......
......@@ -36,10 +36,19 @@ else:
def get_extension_modules():
return [Extension('pypocketfft',
language='c++',
sources=['pypocketfft.cc','mr_util/threading.cc'],
depends=['mr_util/useful_macros.h', 'mr_util/fft.h', 'mr_util/fft1d.h', 'mr_util/mav.h', 'mr_util/threading.h',
'mr_util/aligned_array.h', 'mr_util/simd.h',
'mr_util/cmplx.h', 'mr_util/unity_roots.h', 'mr_util/error_handling.h',
sources=['pypocketfft.cc',
'mr_util/threading.cc'],
depends=['mr_util/useful_macros.h',
'mr_util/fft.h',
'mr_util/fft1d.h',
'mr_util/mav.h',
'mr_util/threading.h',
'mr_util/aligned_array.h',
'mr_util/simd.h',
'mr_util/cmplx.h',
'mr_util/unity_roots.h',
'mr_util/error_handling.h',
'mr_util/pybind_utils.h',
'setup.py'],
include_dirs=include_dirs,
define_macros=define_macros,
......
......@@ -13,6 +13,7 @@ include mr_util/constants.h
include mr_util/unity_roots.h
include mr_util/error_handling.h
include mr_util/useful_macros.h
include mr_util/pybind_utils.h
include libsharp2/sharp.h
include libsharp2/sharp_internal.h
include libsharp2/sharp_geomhelpers.h
......
......@@ -37,6 +37,9 @@
#include "libsharp2/sharp_almhelpers.h"
#include "mr_util/string_utils.h"
#include "mr_util/error_handling.h"
#include "mr_util/mav.h"
#include "mr_util/fft.h"
#include "mr_util/constants.h"
using namespace std;
using namespace mr;
......@@ -45,6 +48,7 @@ namespace py = pybind11;
namespace {
using a_d = py::array_t<double>;
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>;
......@@ -65,32 +69,53 @@ template<typename T> class py_sharpjob
", mmax=" + dataToString(mmax_) + ", npix=", dataToString(npix_) +".>";
}
void set_Gauss_geometry(int64_t nrings, int64_t nphi)
void set_gauss_geometry(int64_t nrings, int64_t nphi)
{
MR_assert((nrings>0)&&(nphi>0),"bad grid dimensions");
npix_=nrings*nphi;
ginfo = sharp_make_gauss_geom_info (nrings, nphi, 0., 1, nphi);
}
void set_Healpix_geometry(int64_t nside)
void set_healpix_geometry(int64_t nside)
{
MR_assert(nside>0,"bad Nside value");
npix_=12*nside*nside;
ginfo = sharp_make_healpix_geom_info (nside, 1);
}
void set_ECP_geometry(int64_t nrings, int64_t nphi)
void set_fejer1_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_ecp_geom_info (nrings, nphi, 0., 1, nphi);
ginfo = sharp_make_fejer1_geom_info (nrings, nphi, 0., 1, nphi);
}
void set_DH_geometry(int64_t nrings, int64_t nphi)
void set_fejer2_geometry(int64_t nrings, int64_t nphi)