diff --git a/README.rst b/README.rst index bf9ea27e..e9036f26 100644 --- a/README.rst +++ b/README.rst @@ -55,9 +55,8 @@ Manual Installation hipercam is written in Python3; it does not support Python2.x It relies on multiple third-party packages, as described in the next -section. At minimum, you should ensure that Cython and trm.pgplot are -installed. Once you have, get the hipercam pipeline software itself -using:: +section. At minimum, you should ensure that trm.pgplot is installed. +Once you have, get the hipercam pipeline software itself using:: git clone https://github.com/HiPERCAM/hipercam.git @@ -110,31 +109,28 @@ For development, you can install the package in editable mode:: This will install the package in development mode, so changes to the source code are immediately reflected without needing to reinstall. -Note: The package includes Cython extensions that need to be compiled, -so you'll need a C compiler and the build dependencies (Cython, numpy) +Note: The package includes compiled extensions that need to be built, +so you'll need a C++ compiler and the build dependencies (numpy, pybind11) available during the build process. Third-Party Modules =================== -Apart from Cython and trm.pgplot, I hope that most of the extras will -get automatically installed if necessary by pip. So, if you have -Cython and PGPLOT ready, you might as well try ``pip install -. --user`` here and now. +Apart from trm.pgplot, I hope that most of the extras will get +automatically installed if necessary by pip. So, if you have PGPLOT ready, +you might as well try ``pip install . --user`` here and now. If something seems amiss, here are details of the third-party packages which you can either install via pip or by looking for the packages in -your O/S package manager. e.g. under fedora, Cython appears as -``python3-Cython``. +your O/S package manager. astropy : astronomical Python package with lots of useful stuff. - Cython : - C-extensions for Python. Widely used package used to interface - to C-libraries and to enable faster code when critical. It is - needed at the setup stage so it might have to be installed first - rather than relying on pip finding it, although I could be wrong. + pybind11 : + Lightweight C++ binding library used to expose compiled helpers for + the fitting and support routines. It is needed at the setup stage so + it may need to be installed before building the package. fitsio : Provides fairly direct access to FITS through the cfitsio library. diff --git a/hipercam/fitting.cpp b/hipercam/fitting.cpp new file mode 100644 index 00000000..0de6313c --- /dev/null +++ b/hipercam/fitting.cpp @@ -0,0 +1,1183 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef _OPENMP +#include +#endif + +namespace py = pybind11; + +// Portable wrapper for compiler-specific restrict support +#if defined(__clang__) || defined(__GNUC__) +#define RESTRICT __restrict__ +#elif defined(_MSC_VER) +#define RESTRICT __restrict +#else +#define RESTRICT +#endif + +namespace { + +// Helper function to generate sub-pixel offsets for binning. +// This creates a vector of offsets for each bin and sub-bin, centered around +// zero, to be used in the sub-pixellation loops in the profile evaluation and +// derivatives. +std::vector make_subpixel_offsets(int bin, int ndiv) { + std::vector offsets; + if (ndiv <= 0) { + return offsets; + } + offsets.reserve(static_cast(bin * ndiv)); + double inv_ndiv = 1.0 / static_cast(ndiv); + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + for (int ibin = 0; ibin < bin; ++ibin) { + double base_off = ibin - (bin - 1) / 2.0 - soff; + for (int isub = 0; isub < ndiv; ++isub) { + offsets.push_back(base_off + isub * inv_ndiv); + } + } + + return offsets; +} + +// Helper function to calculate the Moffat alpha parameter. +double calc_moffat_alpha(double fwhm, double beta) { + double tbeta = std::max(0.01, beta); + return 4.0 * (std::pow(2.0, 1.0 / tbeta) - 1.0) / (fwhm * fwhm); +} + +// Evaluate one Moffat model value at a single pixel coordinate. +// This is used by the selected-pixel residual path to avoid building +// full 2D model arrays when only masked pixels are needed. +double moffat_value_at(double x_val, double y_val, double height, double xcen, + double ycen, double alpha, double tbeta, + const std::vector &x_offsets, + const std::vector &y_offsets) { + if (!x_offsets.empty() && !y_offsets.empty()) { + double prof = 0.0; + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + prof += std::pow(1.0 + alpha * rsq, -tbeta); + } + } + return height * inv_nadd * prof; + } + + double dx = x_val - xcen; + double dy = y_val - ycen; + double rsq = dx * dx + dy * dy; + return height * std::pow(1.0 + alpha * rsq, -tbeta); +} + +// Evaluate one Gaussian model value at a single pixel coordinate. +// This is used by the selected-pixel residual path to avoid building +// full 2D model arrays when only masked pixels are needed. +double gaussian_value_at(double x_val, double y_val, double height, double xcen, + double ycen, double alpha, + const std::vector &x_offsets, + const std::vector &y_offsets) { + if (!x_offsets.empty() && !y_offsets.empty()) { + double prof = 0.0; + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + prof += std::exp(-alpha * rsq); + } + } + return height * inv_nadd * prof; + } + + double dx = x_val - xcen; + double dy = y_val - ycen; + double rsq = dx * dx + dy * dy; + return height * std::exp(-alpha * rsq); +} + +// Calculate Moffat derivatives at one pixel coordinate. +// The outputs follow the same normalization as dmoffat so the fit path +// stays numerically equivalent. +void moffat_derivs_at(double x_val, double y_val, double height, double xcen, + double ycen, double alpha, double tbeta, + double dfwhm_coeff, double dbeta_coeff, + const std::vector &x_offsets, + const std::vector &y_offsets, bool comp_dfwhm, + bool comp_dbeta, double &dheight, double &dxcen, + double &dycen, double &dfwhm, double &dbeta) { + double two_alpha_tbeta = 2.0 * alpha * tbeta; + dheight = 0.0; + dxcen = 0.0; + dycen = 0.0; + dfwhm = 0.0; + dbeta = 0.0; + + if (!x_offsets.empty() && !y_offsets.empty()) { + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + + // Instead of checking comp_dfwhm and comp_dbeta inside the innermost loop, + // we split into four separate loops here to save time on the checks. + if (comp_dfwhm && comp_dbeta) { + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + double denom = 1.0 + alpha * rsq; + // Compute denom^-tbeta once, then derive save1 from it to avoid + // paying for a second pow call in the loop. + double dh = std::pow(denom, -tbeta); + double save1 = height * dh / denom; + double save2 = save1 * rsq; + dheight += dh; + dxcen += two_alpha_tbeta * dx * save1; + dycen += two_alpha_tbeta * dy * save1; + dfwhm += dfwhm_coeff * save2; + double log_denom = std::log(denom); + dbeta += (-log_denom * height * dh + dbeta_coeff * save2); + } + } + } else if (comp_dfwhm) { + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + double denom = 1.0 + alpha * rsq; + double dh = std::pow(denom, -tbeta); + double save1 = height * dh / denom; + double save2 = save1 * rsq; + dheight += dh; + dxcen += two_alpha_tbeta * dx * save1; + dycen += two_alpha_tbeta * dy * save1; + dfwhm += dfwhm_coeff * save2; + } + } + } else if (comp_dbeta) { + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + double denom = 1.0 + alpha * rsq; + double dh = std::pow(denom, -tbeta); + double save1 = height * dh / denom; + double save2 = save1 * rsq; + dheight += dh; + dxcen += two_alpha_tbeta * dx * save1; + dycen += two_alpha_tbeta * dy * save1; + double log_denom = std::log(denom); + dbeta += (-log_denom * height * dh + dbeta_coeff * save2); + } + } + } else { + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + double denom = 1.0 + alpha * rsq; + double dh = std::pow(denom, -tbeta); + double save1 = height * dh / denom; + dheight += dh; + dxcen += two_alpha_tbeta * dx * save1; + dycen += two_alpha_tbeta * dy * save1; + } + } + } + + dheight *= inv_nadd; + dxcen *= inv_nadd; + dycen *= inv_nadd; + if (comp_dfwhm) { + dfwhm *= inv_nadd; + } + if (comp_dbeta) { + dbeta *= inv_nadd; + } + return; + } + + double dx = x_val - xcen; + double dy = y_val - ycen; + double rsq = dx * dx + dy * dy; + double denom = 1.0 + alpha * rsq; + dheight = std::pow(denom, -tbeta); + double save1 = height * dheight / denom; + double save2 = save1 * rsq; + dxcen = two_alpha_tbeta * dx * save1; + dycen = two_alpha_tbeta * dy * save1; + if (comp_dfwhm) { + dfwhm = dfwhm_coeff * save2; + } + if (comp_dbeta) { + double log_denom = std::log(denom); + dbeta = (-log_denom * height * dheight + dbeta_coeff * save2); + } +} + +// Calculate Gaussian derivatives at one pixel coordinate. +// The outputs follow the same normalization as dgaussian so the fit path +// stays numerically equivalent. +void gaussian_derivs_at(double x_val, double y_val, double height, double xcen, + double ycen, double alpha, double two_alpha_height, + double dfwhm_coeff, + const std::vector &x_offsets, + const std::vector &y_offsets, bool comp_dfwhm, + double &dheight, double &dxcen, double &dycen, + double &dfwhm) { + dheight = 0.0; + dxcen = 0.0; + dycen = 0.0; + dfwhm = 0.0; + + if (!x_offsets.empty() && !y_offsets.empty()) { + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + + // Same comp_dfwhm check optimization as in moffat_derivs_at, + // except here we don't have comp_dbeta to worry about. + if (comp_dfwhm) { + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + double dh = std::exp(-alpha * rsq); + dheight += dh; + dxcen += two_alpha_height * dh * dx; + dycen += two_alpha_height * dh * dy; + dfwhm += dfwhm_coeff * dh * rsq; + } + } + } else { + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + double dh = std::exp(-alpha * rsq); + dheight += dh; + dxcen += two_alpha_height * dh * dx; + dycen += two_alpha_height * dh * dy; + } + } + } + + dheight *= inv_nadd; + dxcen *= inv_nadd; + dycen *= inv_nadd; + if (comp_dfwhm) { + dfwhm *= inv_nadd; + } + return; + } + + double dx = x_val - xcen; + double dy = y_val - ycen; + double rsq = dx * dx + dy * dy; + + dheight = std::exp(-alpha * rsq); + dxcen = two_alpha_height * dheight * dx; + dycen = two_alpha_height * dheight * dy; + if (comp_dfwhm) { + dfwhm = dfwhm_coeff * dheight * rsq; + } +} + +// Shared shape validation for the selected-pixel residual/Jacobian entry +// points, which all require a set of same-shaped 2D input arrays. +// py::buffer_info is move-only, so callers pass pointers rather than values. +void check_matching_2d_shapes( + std::initializer_list infos) { + const py::buffer_info &first = **infos.begin(); + bool ok = first.ndim == 2; + for (const py::buffer_info *info : infos) { + ok = ok && info->ndim == 2 && info->shape[0] == first.shape[0] && + info->shape[1] == first.shape[1]; + } + if (!ok) { + throw std::runtime_error( + "Input arrays have invalid dimensions or mismatched shapes"); + } +} + +// Companion check for the 1D ok_indices array passed alongside the 2D +// inputs validated by check_matching_2d_shapes. +void check_ok_indices_1d(const py::buffer_info &ok_info) { + if (ok_info.ndim != 1) { + throw std::runtime_error( + "Input arrays have invalid dimensions or mismatched shapes"); + } +} + +// Bounds-check every selected pixel index against the flattened pixel count. +void check_ok_indices_range(const std::int64_t *ok_ptr, std::size_t n_ok, + std::size_t n_pixels) { + for (std::size_t i = 0; i < n_ok; ++i) { + std::int64_t idx = ok_ptr[i]; + if (idx < 0 || static_cast(idx) >= n_pixels) { + throw std::runtime_error("ok_indices contains out-of-range values"); + } + } +} + +} // namespace + +py::array_t moffat_resid(py::array_t x, py::array_t y, + py::array_t data, + py::array_t sigma, + py::array_t ok_indices, + double sky, double height, double xcen, + double ycen, double fwhm, double beta, + int xbin, int ybin, int ndiv) { + + // Compute residuals only at valid indices provided by Python. This + // bypasses full-frame residual assembly and boolean masking in Python. + + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + py::buffer_info data_info = data.request(); + py::buffer_info sigma_info = sigma.request(); + py::buffer_info ok_info = ok_indices.request(); + + check_matching_2d_shapes({&x_info, &y_info, &data_info, &sigma_info}); + check_ok_indices_1d(ok_info); + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + const double *data_ptr = static_cast(data_info.ptr); + const double *sigma_ptr = static_cast(sigma_info.ptr); + const std::int64_t *ok_ptr = static_cast(ok_info.ptr); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + std::size_t n_ok = ok_info.shape[0]; + + double tbeta = std::max(0.01, beta); + double alpha = calc_moffat_alpha(fwhm, beta); + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + + py::array_t result(static_cast(n_ok)); + py::buffer_info result_info = result.request(); + double *result_ptr = static_cast(result_info.ptr); + + check_ok_indices_range(ok_ptr, n_ok, n_pixels); + + // Release the GIL now the Python operations are complete. + // pybind11 buffer operations (request, array creation) require the GIL, + // so we can only release it after extracting all pointers and dimensions. + py::gil_scoped_release release; + + for (std::size_t i = 0; i < n_ok; ++i) { + std::int64_t idx = ok_ptr[i]; + + double model = + sky + moffat_value_at(x_ptr[idx], y_ptr[idx], height, xcen, ycen, alpha, + tbeta, x_offsets, y_offsets); + result_ptr[i] = (data_ptr[idx] - model) / sigma_ptr[idx]; + } + + return result; +} + +py::array_t dmoffat_jac(py::array_t x, py::array_t y, + py::array_t sigma, + py::array_t ok_indices, + double sky, double height, double xcen, + double ycen, double fwhm, double beta, int xbin, + int ybin, int ndiv, bool comp_dfwhm, + bool comp_dbeta, const std::vector &inds) { + + // sky affects residuals but not derivatives; suppress unused-param warning + (void)sky; + + // Build Jacobian rows directly for selected pixels and selected parameter + // columns (`inds`), matching the same derivative ordering used by Python. + + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + py::buffer_info sigma_info = sigma.request(); + py::buffer_info ok_info = ok_indices.request(); + + check_matching_2d_shapes({&x_info, &y_info, &sigma_info}); + check_ok_indices_1d(ok_info); + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + const double *sigma_ptr = static_cast(sigma_info.ptr); + const std::int64_t *ok_ptr = static_cast(ok_info.ptr); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + std::size_t n_ok = ok_info.shape[0]; + std::size_t n_par = inds.size(); + + double tbeta = std::max(0.01, beta); + double alpha = calc_moffat_alpha(fwhm, beta); + double two_alpha_tbeta = 2.0 * alpha * tbeta; + double dfwhm_coeff = two_alpha_tbeta / fwhm; + double dbeta_coeff = + 4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / tbeta / (fwhm * fwhm); + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + + py::array_t result( + {static_cast(n_ok), static_cast(n_par)}); + py::buffer_info result_info = result.request(); + double *result_ptr = static_cast(result_info.ptr); + + check_ok_indices_range(ok_ptr, n_ok, n_pixels); + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + for (std::size_t i = 0; i < n_ok; ++i) { + std::int64_t idx = ok_ptr[i]; + + double dheight, dxcen, dycen, dfwhm, dbeta; + moffat_derivs_at(x_ptr[idx], y_ptr[idx], height, xcen, ycen, alpha, tbeta, + dfwhm_coeff, dbeta_coeff, x_offsets, y_offsets, comp_dfwhm, + comp_dbeta, dheight, dxcen, dycen, dfwhm, dbeta); + + double d0 = 1.0; + double d1 = dheight; + double d2 = dxcen; + double d3 = dycen; + double d4; + double d5; + + // Keep compatibility with the legacy dmoffat API that duplicates + // placeholders when derivatives are not requested. + if (comp_dfwhm && comp_dbeta) { + d4 = dfwhm; + d5 = dbeta; + } else if (comp_dfwhm) { + d4 = dfwhm; + d5 = dfwhm; + } else if (comp_dbeta) { + d4 = dbeta; + d5 = dbeta; + } else { + d4 = dycen; + d5 = dycen; + } + + double derivs[6] = {d0, d1, d2, d3, d4, d5}; + double inv_sigma = -1.0 / sigma_ptr[idx]; + + for (std::size_t j = 0; j < n_par; ++j) { + int ind = inds[j]; + if (ind < 0 || ind > 5) { + throw std::runtime_error("inds contains out-of-range derivative index"); + } + result_ptr[i * n_par + j] = derivs[ind] * inv_sigma; + } + } + + return result; +} + +py::array_t gaussian_resid(py::array_t x, py::array_t y, + py::array_t data, + py::array_t sigma, + py::array_t ok_indices, + double sky, double height, double xcen, + double ycen, double fwhm, int xbin, int ybin, + int ndiv) { + + // Gaussian equivalent of moffat_resid: selected-pixel residuals only. + + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + py::buffer_info data_info = data.request(); + py::buffer_info sigma_info = sigma.request(); + py::buffer_info ok_info = ok_indices.request(); + + check_matching_2d_shapes({&x_info, &y_info, &data_info, &sigma_info}); + check_ok_indices_1d(ok_info); + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + const double *data_ptr = static_cast(data_info.ptr); + const double *sigma_ptr = static_cast(sigma_info.ptr); + const std::int64_t *ok_ptr = static_cast(ok_info.ptr); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + std::size_t n_ok = ok_info.shape[0]; + + double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + + py::array_t result(static_cast(n_ok)); + py::buffer_info result_info = result.request(); + double *result_ptr = static_cast(result_info.ptr); + + check_ok_indices_range(ok_ptr, n_ok, n_pixels); + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + for (std::size_t i = 0; i < n_ok; ++i) { + std::int64_t idx = ok_ptr[i]; + + double model = sky + gaussian_value_at(x_ptr[idx], y_ptr[idx], height, xcen, + ycen, alpha, x_offsets, y_offsets); + result_ptr[i] = (data_ptr[idx] - model) / sigma_ptr[idx]; + } + + return result; +} + +py::array_t dgaussian_jac(py::array_t x, py::array_t y, + py::array_t sigma, + py::array_t ok_indices, + double sky, double height, double xcen, + double ycen, double fwhm, int xbin, int ybin, + int ndiv, bool comp_dfwhm, + const std::vector &inds) { + + // sky affects residuals but not derivatives; suppress unused-param warning + (void)sky; + + // Selected-pixel Jacobian assembly for Gaussian fits. + + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + py::buffer_info sigma_info = sigma.request(); + py::buffer_info ok_info = ok_indices.request(); + + check_matching_2d_shapes({&x_info, &y_info, &sigma_info}); + check_ok_indices_1d(ok_info); + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + const double *sigma_ptr = static_cast(sigma_info.ptr); + const std::int64_t *ok_ptr = static_cast(ok_info.ptr); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + std::size_t n_ok = ok_info.shape[0]; + std::size_t n_par = inds.size(); + + double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); + double two_alpha_height = 2.0 * alpha * height; + double dfwhm_coeff = two_alpha_height / fwhm; + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + + py::array_t result( + {static_cast(n_ok), static_cast(n_par)}); + py::buffer_info result_info = result.request(); + double *result_ptr = static_cast(result_info.ptr); + + check_ok_indices_range(ok_ptr, n_ok, n_pixels); + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + for (std::size_t i = 0; i < n_ok; ++i) { + std::int64_t idx = ok_ptr[i]; + + double dheight, dxcen, dycen, dfwhm; + gaussian_derivs_at(x_ptr[idx], y_ptr[idx], height, xcen, ycen, alpha, + two_alpha_height, dfwhm_coeff, x_offsets, y_offsets, + comp_dfwhm, dheight, dxcen, dycen, dfwhm); + + double d0 = 1.0; + double d1 = dheight; + double d2 = dxcen; + double d3 = dycen; + double d4 = comp_dfwhm ? dfwhm : dycen; + + double derivs[5] = {d0, d1, d2, d3, d4}; + double inv_sigma = -1.0 / sigma_ptr[idx]; + + for (std::size_t j = 0; j < n_par; ++j) { + int ind = inds[j]; + if (ind < 0 || ind > 4) { + throw std::runtime_error("inds contains out-of-range derivative index"); + } + result_ptr[i * n_par + j] = derivs[ind] * inv_sigma; + } + } + + return result; +} + +// C++ implementation of the Moffat profile function +py::array_t moffat(py::array_t x, py::array_t y, + double sky, double height, double xcen, double ycen, + double fwhm, double beta, int xbin, int ybin, + int ndiv) { + + // Get input array dimensions and data pointers + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || + x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1]) { + throw std::runtime_error("Input arrays must be 2D and have the same shape"); + } + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + + // Create output array with same shape as input + py::array_t result = py::array_t(x_info.shape); + py::buffer_info result_info = result.request(); + double *RESTRICT result_ptr = static_cast(result_info.ptr); + + // Calculate Moffat profile parameters + double tbeta = std::max(0.01, beta); + double alpha = calc_moffat_alpha(fwhm, beta); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + if (ndiv > 0) { + // With sub-pixellation + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + +// OpenMP parallelization of for loop with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + + // Loop over all pixels + for (std::size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + double x_val = x_ptr[pixel_idx]; + double y_val = y_ptr[pixel_idx]; + double prof = 0.0; + + // Reuse precomputed sub-pixel offsets across all pixels. + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + prof += std::pow(1.0 + alpha * rsq, -tbeta); + } + } + + result_ptr[pixel_idx] = sky + height * inv_nadd * prof; + } + } else { + // Fast calculation at pixel centers + for (std::size_t i = 0; i < n_pixels; ++i) { + double dx = x_ptr[i] - xcen; + double dy = y_ptr[i] - ycen; + double rsq = dx * dx + dy * dy; + result_ptr[i] = sky + height * std::pow(1.0 + alpha * rsq, -tbeta); + } + } + + return result; +} + +// C++ implementation of the Moffat profile derivatives +std::vector> +dmoffat(py::array_t x, py::array_t y, double sky, double height, + double xcen, double ycen, double fwhm, double beta, int xbin, int ybin, + int ndiv, bool comp_dfwhm, bool comp_dbeta) { + + // sky affects residuals but not derivatives; suppress unused-param warning + (void)sky; + + // Get input array dimensions and data pointers + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || + x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1]) { + throw std::runtime_error("Input arrays must be 2D and have the same shape"); + } + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + + // Create output arrays with same shape as input + std::vector> result; + result.reserve(6); + + // Always need dsky, dheight, dxcen, dycen + py::array_t dsky = py::array_t(x_info.shape); + py::array_t dheight = py::array_t(x_info.shape); + py::array_t dxcen = py::array_t(x_info.shape); + py::array_t dycen = py::array_t(x_info.shape); + + // Optionally need dfwhm and dbeta + py::array_t dfwhm; + py::array_t dbeta; + + if (comp_dfwhm) { + dfwhm = py::array_t(x_info.shape); + } + + if (comp_dbeta) { + dbeta = py::array_t(x_info.shape); + } + + // Get data pointers + py::buffer_info dsky_info = dsky.request(); + py::buffer_info dheight_info = dheight.request(); + py::buffer_info dxcen_info = dxcen.request(); + py::buffer_info dycen_info = dycen.request(); + + double *RESTRICT dsky_ptr = static_cast(dsky_info.ptr); + double *RESTRICT dheight_ptr = static_cast(dheight_info.ptr); + double *RESTRICT dxcen_ptr = static_cast(dxcen_info.ptr); + double *RESTRICT dycen_ptr = static_cast(dycen_info.ptr); + + py::buffer_info dfwhm_info; + py::buffer_info dbeta_info; + double *RESTRICT dfwhm_ptr = nullptr; + double *RESTRICT dbeta_ptr = nullptr; + + if (comp_dfwhm) { + dfwhm_info = dfwhm.request(); + dfwhm_ptr = static_cast(dfwhm_info.ptr); + } + + if (comp_dbeta) { + dbeta_info = dbeta.request(); + dbeta_ptr = static_cast(dbeta_info.ptr); + } + + // Calculate Moffat profile parameters + double tbeta = std::max(0.01, beta); + double alpha = calc_moffat_alpha(fwhm, beta); + double two_alpha_tbeta = 2.0 * alpha * tbeta; + double dfwhm_coeff = two_alpha_tbeta / fwhm; + double dbeta_coeff = + 4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / tbeta / (fwhm * fwhm); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + // Initialize dsky to ones (derivative of sky is always 1) + std::fill_n(dsky_ptr, n_pixels, 1.0); + + if (ndiv > 0) { + // With sub-pixellation + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + +// OpenMP parallelization of for loop with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + + // Loop over all pixels + for (std::size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + double x_val = x_ptr[pixel_idx]; + double y_val = y_ptr[pixel_idx]; + // Use _sum suffix to avoid shadowing the outer py::array_t declarations. + double dheight_sum = 0.0; + double dxcen_sum = 0.0; + double dycen_sum = 0.0; + double dfwhm_sum = 0.0; + double dbeta_sum = 0.0; + + // Reuse precomputed sub-pixel offsets across all pixels. + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + + double denom = 1.0 + alpha * rsq; + // Derivatives + // Keep the ordering here aligned with moffat_derivs_at: compute + // denom^-tbeta once, then derive save1 from it to avoid a second pow. + double dh = std::pow(denom, -tbeta); + double save1 = height * dh / denom; + double save2 = save1 * rsq; + dheight_sum += dh; + dxcen_sum += two_alpha_tbeta * dx * save1; + dycen_sum += two_alpha_tbeta * dy * save1; + + if (comp_dfwhm) { + dfwhm_sum += dfwhm_coeff * save2; + } + + if (comp_dbeta) { + double log_denom = std::log(denom); + dbeta_sum += (-log_denom * height * dh + dbeta_coeff * save2); + } + } + } + + dheight_ptr[pixel_idx] = dheight_sum * inv_nadd; + dxcen_ptr[pixel_idx] = dxcen_sum * inv_nadd; + dycen_ptr[pixel_idx] = dycen_sum * inv_nadd; + + if (comp_dfwhm) { + dfwhm_ptr[pixel_idx] = dfwhm_sum * inv_nadd; + } + + if (comp_dbeta) { + dbeta_ptr[pixel_idx] = dbeta_sum * inv_nadd; + } + } + } else { + // Fast calculation at pixel centers + for (std::size_t i = 0; i < n_pixels; ++i) { + double dx = x_ptr[i] - xcen; + double dy = y_ptr[i] - ycen; + double rsq = dx * dx + dy * dy; + + double denom = 1.0 + alpha * rsq; + dheight_ptr[i] = std::pow(denom, -tbeta); + double save1 = height * dheight_ptr[i] / denom; + double save2 = save1 * rsq; + + // Derivatives + dxcen_ptr[i] = two_alpha_tbeta * dx * save1; + dycen_ptr[i] = two_alpha_tbeta * dy * save1; + + if (comp_dfwhm) { + dfwhm_ptr[i] = dfwhm_coeff * save2; + } + + if (comp_dbeta) { + double log_denom = std::log(denom); + dbeta_ptr[i] = + (-log_denom * height * dheight_ptr[i] + dbeta_coeff * save2); + } + } + } + + // Add arrays to result vector + result.push_back(dsky); + result.push_back(dheight); + result.push_back(dxcen); + result.push_back(dycen); + + if (comp_dfwhm && comp_dbeta) { + result.push_back(dfwhm); + result.push_back(dbeta); + } else if (comp_dfwhm) { + result.push_back(dfwhm); + result.push_back(dfwhm); // duplicate to match Python API + } else if (comp_dbeta) { + result.push_back(dbeta); + result.push_back(dbeta); // duplicate to match Python API + } else { + result.push_back(dycen); // duplicate to match Python API + result.push_back(dycen); // duplicate to match Python API + } + + return result; +} + +// C++ implementation of the Gaussian profile function +py::array_t gaussian(py::array_t x, py::array_t y, + double sky, double height, double xcen, + double ycen, double fwhm, int xbin, int ybin, + int ndiv) { + + // Get input array dimensions and data pointers + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || + x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1]) { + throw std::runtime_error("Input arrays must be 2D and have the same shape"); + } + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + + // Create output array with same shape as input + py::array_t result = py::array_t(x_info.shape); + py::buffer_info result_info = result.request(); + double *RESTRICT result_ptr = static_cast(result_info.ptr); + + // Calculate Gaussian profile parameter + double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + if (ndiv > 0) { + // With sub-pixellation + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + +// OpenMP parallelization of for loop with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + + // Loop over all pixels + for (std::size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + double x_val = x_ptr[pixel_idx]; + double y_val = y_ptr[pixel_idx]; + double prof = 0.0; + + // Reuse precomputed sub-pixel offsets across all pixels. + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + prof += std::exp(-alpha * rsq); + } + } + + result_ptr[pixel_idx] = sky + height * inv_nadd * prof; + } + } else { + // Fast calculation at pixel centers + for (std::size_t i = 0; i < n_pixels; ++i) { + double dx = x_ptr[i] - xcen; + double dy = y_ptr[i] - ycen; + double rsq = dx * dx + dy * dy; + result_ptr[i] = sky + height * std::exp(-alpha * rsq); + } + } + + return result; +} + +// C++ implementation of the Gaussian profile derivatives +std::vector> +dgaussian(py::array_t x, py::array_t y, double sky, + double height, double xcen, double ycen, double fwhm, int xbin, + int ybin, int ndiv, bool comp_dfwhm) { + + // sky affects residuals but not derivatives; suppress unused-param warning + (void)sky; + + // Get input array dimensions and data pointers + py::buffer_info x_info = x.request(); + py::buffer_info y_info = y.request(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || + x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1]) { + throw std::runtime_error("Input arrays must be 2D and have the same shape"); + } + + const double *x_ptr = static_cast(x_info.ptr); + const double *y_ptr = static_cast(y_info.ptr); + + // Create output arrays with same shape as input + std::vector> result; + result.reserve(5); + + // Always need dsky, dheight, dxcen, dycen + py::array_t dsky = py::array_t(x_info.shape); + py::array_t dheight = py::array_t(x_info.shape); + py::array_t dxcen = py::array_t(x_info.shape); + py::array_t dycen = py::array_t(x_info.shape); + + // Optionally need dfwhm + py::array_t dfwhm; + + if (comp_dfwhm) { + dfwhm = py::array_t(x_info.shape); + } + + // Get data pointers + py::buffer_info dsky_info = dsky.request(); + py::buffer_info dheight_info = dheight.request(); + py::buffer_info dxcen_info = dxcen.request(); + py::buffer_info dycen_info = dycen.request(); + + double *RESTRICT dsky_ptr = static_cast(dsky_info.ptr); + double *RESTRICT dheight_ptr = static_cast(dheight_info.ptr); + double *RESTRICT dxcen_ptr = static_cast(dxcen_info.ptr); + double *RESTRICT dycen_ptr = static_cast(dycen_info.ptr); + + py::buffer_info dfwhm_info; + double *RESTRICT dfwhm_ptr = nullptr; + + if (comp_dfwhm) { + dfwhm_info = dfwhm.request(); + dfwhm_ptr = static_cast(dfwhm_info.ptr); + } + + // Calculate Gaussian profile parameter + double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); + double two_alpha_height = 2.0 * alpha * height; + double dfwhm_coeff = two_alpha_height / fwhm; + + std::size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + // Initialize dsky to ones (derivative of sky is always 1) + std::fill_n(dsky_ptr, n_pixels, 1.0); + + if (ndiv > 0) { + // With sub-pixellation + const std::vector x_offsets = make_subpixel_offsets(xbin, ndiv); + const std::vector y_offsets = make_subpixel_offsets(ybin, ndiv); + double inv_nadd = + 1.0 / static_cast(x_offsets.size() * y_offsets.size()); + +// OpenMP parallelization of for loop with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + + // Loop over all pixels + for (std::size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + double x_val = x_ptr[pixel_idx]; + double y_val = y_ptr[pixel_idx]; + // Use _sum suffix to avoid shadowing the outer py::array_t declarations. + double dheight_sum = 0.0; + double dxcen_sum = 0.0; + double dycen_sum = 0.0; + double dfwhm_sum = 0.0; + + // Reuse precomputed sub-pixel offsets across all pixels. + for (double yoff : y_offsets) { + for (double xoff : x_offsets) { + double dx = x_val + xoff - xcen; + double dy = y_val + yoff - ycen; + double rsq = dx * dx + dy * dy; + + // Gaussian value + double dh = std::exp(-alpha * rsq); + dheight_sum += dh; + dxcen_sum += two_alpha_height * dh * dx; + dycen_sum += two_alpha_height * dh * dy; + + if (comp_dfwhm) { + dfwhm_sum += dfwhm_coeff * dh * rsq; + } + } + } + + dheight_ptr[pixel_idx] = dheight_sum * inv_nadd; + dxcen_ptr[pixel_idx] = dxcen_sum * inv_nadd; + dycen_ptr[pixel_idx] = dycen_sum * inv_nadd; + + if (comp_dfwhm) { + dfwhm_ptr[pixel_idx] = dfwhm_sum * inv_nadd; + } + } + } else { + // Fast calculation at pixel centers + for (std::size_t i = 0; i < n_pixels; ++i) { + double dx = x_ptr[i] - xcen; + double dy = y_ptr[i] - ycen; + double rsq = dx * dx + dy * dy; + + // Gaussian value + double dh = std::exp(-alpha * rsq); + dheight_ptr[i] = dh; + dxcen_ptr[i] = two_alpha_height * dh * dx; + dycen_ptr[i] = two_alpha_height * dh * dy; + + if (comp_dfwhm) { + dfwhm_ptr[i] = dfwhm_coeff * dh * rsq; + } + } + } + + // Add arrays to result vector + result.push_back(dsky); + result.push_back(dheight); + result.push_back(dxcen); + result.push_back(dycen); + + if (comp_dfwhm) { + result.push_back(dfwhm); + } else { + result.push_back(dycen); // duplicate to match Python API + } + + return result; +} + +PYBIND11_MODULE(_fitting_cpp, m) { + m.doc() = "C++ implementation of profile fitting functions"; + + m.def("moffat", &moffat, "C++ implementation of Moffat profile", py::arg("x"), + py::arg("y"), py::arg("sky"), py::arg("height"), py::arg("xcen"), + py::arg("ycen"), py::arg("fwhm"), py::arg("beta"), py::arg("xbin"), + py::arg("ybin"), py::arg("ndiv")); + + m.def("dmoffat", &dmoffat, "C++ implementation of Moffat profile derivatives", + py::arg("x"), py::arg("y"), py::arg("sky"), py::arg("height"), + py::arg("xcen"), py::arg("ycen"), py::arg("fwhm"), py::arg("beta"), + py::arg("xbin"), py::arg("ybin"), py::arg("ndiv"), + py::arg("comp_dfwhm"), py::arg("comp_dbeta")); + + m.def("moffat_resid", &moffat_resid, + "C++ implementation of Moffat residuals at selected pixels", + py::arg("x"), py::arg("y"), py::arg("data"), py::arg("sigma"), + py::arg("ok_indices"), py::arg("sky"), py::arg("height"), + py::arg("xcen"), py::arg("ycen"), py::arg("fwhm"), py::arg("beta"), + py::arg("xbin"), py::arg("ybin"), py::arg("ndiv")); + + m.def("dmoffat_jac", &dmoffat_jac, + "C++ implementation of Moffat residual Jacobian at selected pixels", + py::arg("x"), py::arg("y"), py::arg("sigma"), py::arg("ok_indices"), + py::arg("sky"), py::arg("height"), py::arg("xcen"), py::arg("ycen"), + py::arg("fwhm"), py::arg("beta"), py::arg("xbin"), py::arg("ybin"), + py::arg("ndiv"), py::arg("comp_dfwhm"), py::arg("comp_dbeta"), + py::arg("inds")); + + m.def("gaussian", &gaussian, "C++ implementation of Gaussian profile", + py::arg("x"), py::arg("y"), py::arg("sky"), py::arg("height"), + py::arg("xcen"), py::arg("ycen"), py::arg("fwhm"), py::arg("xbin"), + py::arg("ybin"), py::arg("ndiv")); + + m.def("dgaussian", &dgaussian, + "C++ implementation of Gaussian profile derivatives", py::arg("x"), + py::arg("y"), py::arg("sky"), py::arg("height"), py::arg("xcen"), + py::arg("ycen"), py::arg("fwhm"), py::arg("xbin"), py::arg("ybin"), + py::arg("ndiv"), py::arg("comp_dfwhm")); + + m.def("gaussian_resid", &gaussian_resid, + "C++ implementation of Gaussian residuals at selected pixels", + py::arg("x"), py::arg("y"), py::arg("data"), py::arg("sigma"), + py::arg("ok_indices"), py::arg("sky"), py::arg("height"), + py::arg("xcen"), py::arg("ycen"), py::arg("fwhm"), py::arg("xbin"), + py::arg("ybin"), py::arg("ndiv")); + + m.def("dgaussian_jac", &dgaussian_jac, + "C++ implementation of Gaussian residual Jacobian at selected pixels", + py::arg("x"), py::arg("y"), py::arg("sigma"), py::arg("ok_indices"), + py::arg("sky"), py::arg("height"), py::arg("xcen"), py::arg("ycen"), + py::arg("fwhm"), py::arg("xbin"), py::arg("ybin"), py::arg("ndiv"), + py::arg("comp_dfwhm"), py::arg("inds")); +} diff --git a/hipercam/fitting.py b/hipercam/fitting.py index 285c4466..1f830745 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -4,12 +4,27 @@ Moffat profiles plus constants. """ +import importlib + from numba import jit import numpy as np from scipy.optimize import least_squares from .core import * from .window import * -from . import support + +try: + _fitting_cpp = importlib.import_module("._fitting_cpp", __package__) + _moffat_cpp = _fitting_cpp.moffat + _moffat_resid = _fitting_cpp.moffat_resid + _dmoffat_cpp = _fitting_cpp.dmoffat + _dmoffat_jac = _fitting_cpp.dmoffat_jac + _gaussian_cpp = _fitting_cpp.gaussian + _gaussian_resid = _fitting_cpp.gaussian_resid + _dgaussian_cpp = _fitting_cpp.dgaussian + _dgaussian_jac = _fitting_cpp.dgaussian_jac + FITTING_CCP_AVAILABLE = True +except ImportError: + FITTING_CCP_AVAILABLE = False __all__ = ("combFit", "fitMoffat", "fitGaussian", "moffat", "gaussian") @@ -30,7 +45,8 @@ def combFit( beta_fix, thresh, ndiv=0, - max_nfev=None + max_nfev=None, + ls_tol=1e-8, ): """Fits a stellar profile in a :class:Window using either a 2D Gaussian or Moffat profile. This is a convenience wrapper of fitMoffat and @@ -96,6 +112,15 @@ def combFit( will slow things. To simply evaluate the profile once at the centre of each pixel in `wind`, set ndiv = 0. + max_nfev : int or None + maximum number of function evaluations during fits. + Passed directly to scipy.optimize.least_squares. + + ls_tol : float or None + tolerance for least squares termination. + Used to set ftol, xtol and gtol in scipy.optimize.least_squares. + + Returns:: (pars, epars, extras) where:: @@ -130,7 +155,7 @@ def combFit( (fit, X, Y, chisq, nok, nrej, npar, nfev) ) = fitGaussian( wind, sigma, sky, height, x, y, fwhm, fwhm_min, fwhm_fix, - thresh, ndiv, max_nfev + thresh, ndiv, max_nfev, ls_tol ) elif method == "m": @@ -141,7 +166,7 @@ def combFit( (fit, X, Y, chisq, nok, nrej, npar, nfev) ) = fitMoffat( wind, sigma, sky, height, x, y, fwhm, fwhm_min, fwhm_fix, beta, - beta_max, beta_fix, thresh, ndiv, max_nfev + beta_max, beta_fix, thresh, ndiv, max_nfev, ls_tol ) else: @@ -199,6 +224,7 @@ def fitMoffat( thresh, ndiv, max_nfev=None, + ls_tol=1e-8, ): """Fits the profile of one target in a Window with a symmetric 2D Moffat profile plus a constant "c + h/(1+alpha**2)**beta" where r is the distance @@ -289,8 +315,12 @@ def fitMoffat( `wind`, set ndiv = 0. max_nfev : int or None - maximum number of function evaluations during fits. Passed - direct to least_squares. + maximum number of function evaluations during fits. + Passed directly to scipy.optimize.least_squares. + + ls_tol : float or None + tolerance for least squares termination. + Used to set ftol, xtol and gtol in scipy.optimize.least_squares. Returns:: tuple @@ -357,7 +387,14 @@ def fitMoffat( # carry out fit res = least_squares( - mfit.fun, param, jac=mfit.jac, method="lm", max_nfev=max_nfev + mfit.fun, + param, + jac=mfit.jac, + method="lm", + max_nfev=max_nfev, + ftol=ls_tol, + xtol=ls_tol, + gtol=ls_tol, ) if not res.success: raise HipercamError(res.message) @@ -405,7 +442,7 @@ def fitMoffat( # first fit carried out with higher threshold for safety sigma[ok & (np.abs(resid) > 2*sfac*thresh)] *= -1 else: - sigma[ok & (np.abs(resid) > 2*sfac*thresh)] *= -1 + sigma[ok & (np.abs(resid) > sfac*thresh)] *= -1 # check whether any have been rejected ok = mfit.mask & (sigma > 0) @@ -451,7 +488,6 @@ def fitMoffat( ) -@jit(nopython=True, cache=True) def moffat(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv): """ Returns a numpy array corresponding to the ordinate grids in xy @@ -461,6 +497,10 @@ def moffat(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv): beta. As beta becomes large, this tends to a Gaussian shape but has more extended wings at low beta. + This function will use the fastest available implementation: + 1. C++ implementation (if available) + 2. Numba-accelerated Python implementation (fallback) + Parameters: x : 2D numpy array @@ -506,10 +546,18 @@ def moffat(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv): obviously will slow things. To simply evaluate the profile once at the centre of each pixel, set ndiv = 0. - Returns:: 2D numpy array containg the Moffat profile plus constant evaluated + Returns:: 2D numpy array containing the Moffat profile plus constant evaluated on the ordinate grids in xy. """ + if FITTING_CCP_AVAILABLE: + return _moffat_cpp(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv) + + return _moffat_numba(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv) + + +@jit(nopython=True, cache=True) +def _moffat_numba(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv): tbeta = max(0.01, beta) alpha = 4 * (2 ** (1.0 / tbeta) - 1) / fwhm ** 2 @@ -545,7 +593,6 @@ def moffat(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv): return height * (1 + alpha * rsq) ** (-tbeta) + sky -@jit(nopython=True, cache=True) def dmoffat( x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv, comp_dfwhm, comp_dbeta ): @@ -557,6 +604,10 @@ def dmoffat( a Gaussian shape but has more extended wings at low beta. The partial derivatives are in the order of the parameters. + This function will use the fastest available implementation: + 1. C++ implementation (if available) + 2. Numba-accelerated Python implementation (fallback) + Parameters: x : 2D numpy array @@ -618,6 +669,44 @@ def dmoffat( numba just-in-time compiler function better. """ + if FITTING_CCP_AVAILABLE: + return _dmoffat_cpp( + x, + y, + sky, + height, + xcen, + ycen, + fwhm, + beta, + xbin, + ybin, + ndiv, + comp_dfwhm, + comp_dbeta, + ) + + return _dmoffat_numba( + x, + y, + sky, + height, + xcen, + ycen, + fwhm, + beta, + xbin, + ybin, + ndiv, + comp_dfwhm, + comp_dbeta, + ) + + +@jit(nopython=True, cache=True) +def _dmoffat_numba( + x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv, comp_dfwhm, comp_dbeta +): tbeta = max(0.01, beta) alpha = 4 * (2 ** (1 / tbeta) - 1) / fwhm ** 2 @@ -757,7 +846,7 @@ def _mask(wind, x, y): class Mfit: """Object providing 'fun' and 'jac' methods for least_squares for - Mofffat models. Eight operating modes each of which allows the + Moffat models. Eight operating modes each of which allows the following to be free [else not]. Sky assumed = 0 when not free. mode == 'sfb' : sky, FWHM, beta @@ -788,7 +877,6 @@ def __init__(self, wind, sigma, ndiv, mode, fwhm=None, beta=None): ndiv : int pixel sub-division factor. See comments in fitMoffat """ - self.sigma = sigma x = wind.x(np.arange(wind.nx)) y = wind.y(np.arange(wind.ny)) self.x, self.y = np.meshgrid(x, y) @@ -796,21 +884,41 @@ def __init__(self, wind, sigma, ndiv, mode, fwhm=None, beta=None): self.xbin = wind.xbin self.ybin = wind.ybin self.ndiv = ndiv + self.sigma = sigma self.mask = _mask(wind, self.x, self.y) + self.ok = self.mask & (self.sigma > 0) + self.ok_indices = np.flatnonzero(self.ok.ravel()).astype(np.int64) + self.set_mode(mode, fwhm, beta) def set_mode(self, mode, fwhm, beta): """Set the operation mode with some light checks""" if mode not in ("sfb", "sb", "sf", "s", "fb", "b", "f", ""): raise HipercamError("invalid mode = {:s}".format(mode)) - - if (mode.find("f") == -1 and fwhm is None) or ( - mode.find("b") == -1 and beta is None - ): - raise HipercamError("invalid mode / fwhm / beta combination") self.mode = mode + self.fwhm = fwhm self.beta = beta + self.comp_fwhm = mode.find("f") > -1 + self.comp_beta = mode.find("b") > -1 + if (self.fwhm is None and not self.comp_fwhm) or ( + self.beta is None and not self.comp_beta + ): + raise HipercamError("invalid mode / fwhm / beta combination") + + # Precompute derivative indices (moved out of jac() for performance) + if mode == "sfb": + self.inds = (0, 1, 2, 3, 4, 5) + elif mode == "sb" or mode == "sf": + self.inds = (0, 1, 2, 3, 4) + elif mode == "s": + self.inds = (0, 1, 2, 3) + elif mode == "fb": + self.inds = (1, 2, 3, 4, 5) + elif mode == "b" or mode == "f": + self.inds = (1, 2, 3, 4) + elif mode == "": + self.inds = (1, 2, 3) def get_par(self, param): """Gets parameters (sky, height, xcen, ycen, fwhm, beta) according to @@ -858,13 +966,13 @@ def set_par(self, sky, height, xcen, ycen, fwhm, beta): Argument:: - param : 1D array - unpacks to (sky, height, xcen, ycen, fwhm, beta) where: - 'sky' is the background per pixel; 'height' is the - central height of the Moffat function; 'xcen' and 'ycen' - are the ordinates of its centre in unbinned CCD pixels - with (1,1) at the left corner of the physical imagine - area; 'fwhm' is the FWHM in unbinned + param : 1D array + unpacks to (sky, height, xcen, ycen, fwhm, beta) where: + 'sky' is the background per pixel; 'height' is the + central height of the Moffat function; 'xcen' and 'ycen' + are the ordinates of its centre in unbinned CCD pixels + with (1,1) at the left corner of the physical imagine + area; 'fwhm' is the FWHM in unbinned """ if self.mode == "sfb": @@ -878,7 +986,7 @@ def set_par(self, sky, height, xcen, ycen, fwhm, beta): elif self.mode == "fb": return (height, xcen, ycen, fwhm, beta) elif self.mode == "b": - reurn(height, xcen, ycen, beta) + return (height, xcen, ycen, beta) elif self.mode == "f": return (height, xcen, ycen, fwhm) elif self.mode == "": @@ -924,42 +1032,24 @@ def get_epar(self, param): return (skye, heighte, xcene, ycene, fwhme, betae) - def fun(self, param): - """ - Returns 1D array of normalised residuals. See the model - method for a description of the argument 'param' + def model(self, param): """ - mod = self.model(param) - diff = (self.data - mod) / self.sigma - ok = self.mask & (self.sigma > 0) - return diff[ok].ravel() + Returns 2D array with model given a parameter vector. - def jac(self, param): - """ - Returns list of 1D arrays of the partial derivatives of - the normalised residuals with respect to the variable - parameters. + Argument:: + + param : 1D array + parameter vector, with values that depend upon the mode, + but could include some or all of (sky, height, xcen, ycen, + fwhm, beta) where 'sky' is the background per pixel; + 'height' is the central height of the Moffat function; + 'xcen' and 'ycen' are the ordinates of its centre in + unbinned CCD pixels with (1,1) at the left corner of the + physical imaging area; 'fwhm' is the FWHM in unbinned + pixels; 'beta' is the Moffat exponent. """ sky, height, xcen, ycen, fwhm, beta = self.get_par(param) - - comp_fwhm = self.mode.find("f") > -1 - comp_beta = self.mode.find("b") > -1 - - # work out which derivatives to bother with - if self.mode == "sfb": - inds = (0, 1, 2, 3, 4, 5) - elif self.mode == "sb" or self.mode == "sf": - inds = (0, 1, 2, 3, 4) - elif self.mode == "s": - inds = (0, 1, 2, 3) - elif self.mode == "fb": - inds = (1, 2, 3, 4, 5) - elif self.mode == "b" or self.mode == "f": - inds = (1, 2, 3, 4) - elif self.mode == "": - inds = (1, 2, 3) - - derivs = dmoffat( + return moffat( self.x, self.y, sky, @@ -971,33 +1061,17 @@ def jac(self, param): self.xbin, self.ybin, self.ndiv, - comp_fwhm, - comp_beta, ) - ok = self.mask & (self.sigma > 0) - return np.column_stack( - [(-derivs[ind][ok] / self.sigma[ok]).ravel() for ind in inds] - ) - - def model(self, param): + def dmodel(self, param): """ - Returns 2D array with model given a parameter vector. + Returns list of 2D arrays of the partial derivatives of + the model with respect to the variable parameters. - Argument:: - - param : 1D array - parameter vector, with values that depend upon the mode, - but could include some or all of (sky, height, xcen, ycen, - fwhm, beta) where 'sky' is the background per pixel; - 'height' is the central height of the Moffat function; - 'xcen' and 'ycen' are the ordinates of its centre in - unbinned CCD pixels with (1,1) at the left corner of the - physical imaging area; 'fwhm' is the FWHM in unbinned - pixels; 'beta' is the Moffat exponent. + See the model method for a description of the argument 'param'. """ sky, height, xcen, ycen, fwhm, beta = self.get_par(param) - return moffat( + return dmoffat( self.x, self.y, sky, @@ -1009,8 +1083,74 @@ def model(self, param): self.xbin, self.ybin, self.ndiv, + self.comp_fwhm, + self.comp_beta, ) + def fun(self, param): + """ + Returns 1D array of normalised residuals. + See the model method for a description of the argument 'param'. + + Used by scipy.optimize.least_squares. + """ + if FITTING_CCP_AVAILABLE: + sky, height, xcen, ycen, fwhm, beta = self.get_par(param) + return _moffat_resid( + self.x, + self.y, + self.data, + self.sigma, + self.ok_indices, + sky, + height, + xcen, + ycen, + fwhm, + beta, + self.xbin, + self.ybin, + self.ndiv, + ) + + mod = self.model(param) + diff = (self.data - mod) / self.sigma + return diff[self.ok].ravel() + + def jac(self, param): + """ + Returns list of 1D arrays of the partial derivatives of + the normalised residuals with respect to the variable + parameters. + + Used by scipy.optimize.least_squares. + """ + if FITTING_CCP_AVAILABLE: + sky, height, xcen, ycen, fwhm, beta = self.get_par(param) + return _dmoffat_jac( + self.x, + self.y, + self.sigma, + self.ok_indices, + sky, + height, + xcen, + ycen, + fwhm, + beta, + self.xbin, + self.ybin, + self.ndiv, + self.comp_fwhm, + self.comp_beta, + self.inds, + ) + + derivs = self.dmodel(param) + arr = np.stack([derivs[ind][self.ok] for ind in self.inds], axis=-1) + return -arr / self.sigma[self.ok, None] + + ########################################## # @@ -1031,7 +1171,8 @@ def fitGaussian( fwhm_fix, thresh, ndiv, - max_nfev=0, + max_nfev=None, + ls_tol=1e-8, ): """Fits the profile of one target in an Window with a 2D symmetric Gaussian profile "c + h*exp(-alpha*r**2)" where r is the distance from the centre @@ -1107,9 +1248,14 @@ def fitGaussian( simply evaluate the profile once at the centre of each pixel in `wind`, set ndiv = 0. - max_nfev : int - maximum number of function evaluations during fits. Passed directly - to leastsq. + max_nfev : int or None + maximum number of function evaluations during fits. + Passed directly to scipy.optimize.least_squares. + + ls_tol : float or None + tolerance for least squares termination. + Used to set ftol, xtol and gtol in scipy.optimize.least_squares. + Returns:: tuple of tuples @@ -1163,11 +1309,18 @@ def fitGaussian( # carry out fit res = least_squares( - gfit.fun, param, jac=gfit.jac, method="lm", max_nfev=max_nfev + gfit.fun, + param, + jac=gfit.jac, + method="lm", + max_nfev=max_nfev, + ftol=ls_tol, + xtol=ls_tol, + gtol=ls_tol, ) if not res.success: raise HipercamError(res.message) - nfev += nfev + nfev += res.nfev # get Jacobian J = np.matrix(res.jac) @@ -1253,7 +1406,7 @@ def fitGaussian( extras ) -@jit(nopython=True, cache=True) + def gaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): """Returns a numpy array corresponding to the ordinate grids in xy set to a symmetric 2D Gaussian plus a constant. The profile is essentially defined @@ -1261,6 +1414,9 @@ def gaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): and alpha is set to give the desired FWHM, but account is taken of the finite size of the pixels by summing over multiple points in each one. + This function will use the optimized C++ implementation if available, + falling back to a Numba-accelerated Python implementation if not. + Arguments:: x : 2D numpy array @@ -1307,6 +1463,14 @@ def gaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): on the ordinate grids in xy. """ + if FITTING_CCP_AVAILABLE: + return _gaussian_cpp(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv) + + return _gaussian_numba(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv) + + +@jit(nopython=True, cache=True) +def _gaussian_numba(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): alpha = 4.0 * np.log(2.0) / fwhm ** 2 @@ -1340,7 +1504,6 @@ def gaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): return sky + height * np.exp(-alpha * rsq) -@jit(nopython=True, cache=True) def dgaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm): """Returns a list of four or five numpy arrays corresponding to the ordinate grids in xy set to the partial derivatives of a symmetric 2D Gaussian plus @@ -1348,6 +1511,9 @@ def dgaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) distance from the centre and alpha is set to give the desired FWHM. The partial derivatives are in the order of the parameters. + This function will use the optimized C++ implementation if available, + falling back to a Numba-accelerated Python implementation if not. + Arguments:: x : 2D numpy array @@ -1400,6 +1566,14 @@ def dgaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) appear in the function call. """ + if FITTING_CCP_AVAILABLE: + return _dgaussian_cpp(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) + + return _dgaussian_numba(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) + + +@jit(nopython=True, cache=True) +def _dgaussian_numba(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm): alpha = 4.0 * np.log(2.0) / fwhm ** 2 dsky = np.ones_like(x) @@ -1468,6 +1642,7 @@ def dgaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) else: return (dsky, dheight, dxcen, dycen, dycen) + class Gfit: """Object providing 'fun' and 'jac' methods for least_squares for Gaussian models. Four operating modes each of which allows the @@ -1503,7 +1678,6 @@ def __init__(self, wind, sigma, ndiv, mode, fwhm=None): FWHM in unbinned pixels """ - self.sigma = sigma x = wind.x(np.arange(wind.nx)) y = wind.y(np.arange(wind.ny)) self.x, self.y = np.meshgrid(x, y) @@ -1511,18 +1685,32 @@ def __init__(self, wind, sigma, ndiv, mode, fwhm=None): self.xbin = wind.xbin self.ybin = wind.ybin self.ndiv = ndiv + self.sigma = sigma self.mask = _mask(wind, self.x, self.y) + self.ok = self.mask & (self.sigma > 0) + self.ok_indices = np.flatnonzero(self.ok.ravel()).astype(np.int64) self.set_mode(mode, fwhm) def set_mode(self, mode, fwhm): """Set the operation mode with some light checks""" if mode not in ("sf", "s", "f", ""): raise HipercamError("invalid mode = {}".format(mode)) - - if mode.find("f") == -1 and fwhm is None: - raise HipercamError("invalid mode / fwhm combination") self.mode = mode + self.fwhm = fwhm + self.comp_fwhm = mode.find("f") > -1 + if fwhm is None and not self.comp_fwhm: + raise HipercamError("invalid mode / fwhm combination") + + # Precompute derivative indices (moved out of jac() for performance) + if mode == "sf": + self.inds = (0, 1, 2, 3, 4) + elif mode == "s": + self.inds = (0, 1, 2, 3) + elif mode == "f": + self.inds = (1, 2, 3, 4) + elif mode == "": + self.inds = (1, 2, 3) def get_par(self, param): """Gets parameters (sky, height, xcen, ycen, fwhm) according to @@ -1588,39 +1776,24 @@ def get_epar(self, param): return (skye, heighte, xcene, ycene, fwhme) - def fun(self, param): - """ - Returns 1D array of normalised residuals. See the model - method for a description of the argument 'param' + def model(self, param): """ - mod = self.model(param) - diff = (self.data - mod) / self.sigma - ok = self.mask & (self.sigma > 0) - return diff[ok].ravel() + Returns 2D array with model given a parameter vector. - def jac(self, param): - """ - Returns list of 1D arrays of the partial derivatives of - the normalised residuals with respect to the variable - parameters. + Argument:: + + param : 1D array + parameter vector, with values that depend upon the mode, + but could include some or all of (sky, height, xcen, ycen, + fwhm) where 'sky' is the background per pixel; + 'height' is the central height of the gaussian function; + 'xcen' and 'ycen' are the ordinates of its centre in + unbinned CCD pixels with (1,1) at the left corner of the + physical imaging area; 'fwhm' is the FWHM in unbinned + pixels. """ sky, height, xcen, ycen, fwhm = self.get_par(param) - - comp_fwhm = self.mode.find("f") > -1 - - # work out which derivatives to bother with - if self.mode == "sf": - inds = (0, 1, 2, 3, 4) - elif self.mode == "s": - inds = (0, 1, 2, 3) - elif self.mode == "f": - inds = (1, 2, 3, 4) - elif self.mode == "": - inds = (1, 2, 3) - else: - raise HipercamError("invalid mode") - - derivs = dgaussian( + return gaussian( self.x, self.y, sky, @@ -1631,32 +1804,17 @@ def jac(self, param): self.xbin, self.ybin, self.ndiv, - comp_fwhm, - ) - - ok = self.mask & (self.sigma > 0) - return np.column_stack( - [(-derivs[ind][ok] / self.sigma[ok]).ravel() for ind in inds] ) - def model(self, param): + def dmodel(self, param): """ - Returns 2D array with model given a parameter vector. + Returns list of 2D arrays of the partial derivatives of + the model with respect to the variable parameters. - Argument:: - - param : 1D array - parameter vector, with values that depend upon the mode, - but could include some or all of (sky, height, xcen, ycen, - fwhm) where 'sky' is the background per pixel; - 'height' is the central height of the gaussian function; - 'xcen' and 'ycen' are the ordinates of its centre in - unbinned CCD pixels with (1,1) at the left corner of the - physical imaging area; 'fwhm' is the FWHM in unbinned - pixels. + See the model method for a description of the argument 'param'. """ sky, height, xcen, ycen, fwhm = self.get_par(param) - return gaussian( + return dgaussian( self.x, self.y, sky, @@ -1667,5 +1825,65 @@ def model(self, param): self.xbin, self.ybin, self.ndiv, + self.comp_fwhm, ) + def fun(self, param): + """ + Returns 1D array of normalised residuals. + See the model method for a description of the argument 'param'. + + Used by scipy.optimize.least_squares. + """ + if FITTING_CCP_AVAILABLE: + sky, height, xcen, ycen, fwhm = self.get_par(param) + return _gaussian_resid( + self.x, + self.y, + self.data, + self.sigma, + self.ok_indices, + sky, + height, + xcen, + ycen, + fwhm, + self.xbin, + self.ybin, + self.ndiv, + ) + + mod = self.model(param) + diff = (self.data - mod) / self.sigma + return diff[self.ok].ravel() + + def jac(self, param): + """ + Returns list of 1D arrays of the partial derivatives of + the normalised residuals with respect to the variable + parameters. + + Used by scipy.optimize.least_squares. + """ + if FITTING_CCP_AVAILABLE: + sky, height, xcen, ycen, fwhm = self.get_par(param) + return _dgaussian_jac( + self.x, + self.y, + self.sigma, + self.ok_indices, + sky, + height, + xcen, + ycen, + fwhm, + self.xbin, + self.ybin, + self.ndiv, + self.comp_fwhm, + self.inds, + ) + + derivs = self.dmodel(param) + arr = np.stack([derivs[ind][self.ok] for ind in self.inds], axis=-1) + return -arr / self.sigma[self.ok, None] diff --git a/hipercam/support.cpp b/hipercam/support.cpp new file mode 100644 index 00000000..40e35f7c --- /dev/null +++ b/hipercam/support.cpp @@ -0,0 +1,136 @@ +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +namespace { + +// Type alias for the 3D float array +using Array3F = py::array_t; + +// Compute the average and standard deviation of a 3D float array along the +// first axis, ignoring outliers beyond a specified sigma threshold. Returns a +// tuple of (avg, std, num), where avg and std are 2D arrays of shape (ny, nx) +// and num is a 2D array of shape (ny, nx) containing the number of valid pixels +// used in the computation for each (y, x) position. +py::tuple avgstd_impl(const py::array_t &cube, float sigma) { + if (cube.ndim() != 3) { + throw std::runtime_error("cube must be 3-dimensional"); + } + if (sigma <= 1.0f) { + throw std::runtime_error("sigma must be greater than 1"); + } + + auto buf = cube.request(); + const std::size_t nf = static_cast(buf.shape[0]); + const std::size_t ny = static_cast(buf.shape[1]); + const std::size_t nx = static_cast(buf.shape[2]); + const float *data = static_cast(buf.ptr); + + py::array_t avg({ny, nx}); + py::array_t stddev({ny, nx}); + py::array_t num({ny, nx}); + + auto avg_mut = avg.mutable_unchecked<2>(); + auto std_mut = stddev.mutable_unchecked<2>(); + auto num_mut = num.mutable_unchecked<2>(); + + std::vector vals(nf); + std::vector ok(nf, 1); + + for (std::size_t iy = 0; iy < ny; ++iy) { + for (std::size_t ix = 0; ix < nx; ++ix) { + for (std::size_t iz = 0; iz < nf; ++iz) { + vals[iz] = data[iz * ny * nx + iy * nx + ix]; + ok[iz] = 1; + } + + std::size_t ncur = nf; + while (true) { + if (ncur == 0) { + avg_mut(iy, ix) = 0.0f; + std_mut(iy, ix) = 0.0f; + num_mut(iy, ix) = 0; + break; + } + + double sum = 0.0; + std::size_t nused = 0; + for (std::size_t iz = 0; iz < nf; ++iz) { + if (!ok[iz]) { + continue; + } + sum += vals[iz]; + ++nused; + } + + double tavg = nused > 0 ? sum / static_cast(nused) : 0.0; + double sumsq = 0.0; + for (std::size_t iz = 0; iz < nf; ++iz) { + if (!ok[iz]) { + continue; + } + const double diff = vals[iz] - tavg; + sumsq += diff * diff; + } + + double tstd = + nused > 1 ? std::sqrt(sumsq / static_cast(nused - 1)) : 0.0; + const double thresh = sigma * tstd; + + std::vector new_ok(ok.begin(), ok.end()); + std::size_t nnew = 0; + for (std::size_t iz = 0; iz < nf; ++iz) { + const bool keep = ok[iz] && std::fabs(vals[iz] - tavg) <= thresh; + new_ok[iz] = keep ? 1 : 0; + if (keep) { + ++nnew; + } + } + + if (nnew == ncur) { + double final_sum = 0.0; + for (std::size_t iz = 0; iz < nf; ++iz) { + if (!new_ok[iz]) { + continue; + } + final_sum += vals[iz]; + } + const double final_avg = + nnew > 0 ? final_sum / static_cast(nnew) : 0.0; + double final_sumsq = 0.0; + for (std::size_t iz = 0; iz < nf; ++iz) { + if (!new_ok[iz]) { + continue; + } + const double diff = vals[iz] - final_avg; + final_sumsq += diff * diff; + } + const double final_std = + nnew > 1 ? std::sqrt(final_sumsq / static_cast(nnew - 1)) + : 0.0; + avg_mut(iy, ix) = static_cast(final_avg); + std_mut(iy, ix) = static_cast(final_std); + num_mut(iy, ix) = static_cast(nnew); + break; + } + + ok.swap(new_ok); + ncur = nnew; + } + } + } + + return py::make_tuple(avg, stddev, num); +} + +} // namespace + +PYBIND11_MODULE(_support_cpp, m) { + m.def("avgstd", &avgstd_impl, py::arg("cube"), py::arg("sigma")); +} diff --git a/hipercam/support.py b/hipercam/support.py new file mode 100644 index 00000000..7f935b35 --- /dev/null +++ b/hipercam/support.py @@ -0,0 +1,90 @@ +"""Support routines for image combination and profile evaluation. + +The implementation uses a pybind11-backed extension when available and falls +back to a pure-NumPy implementation otherwise. +""" + +from __future__ import annotations +import importlib + +import numpy as np + +try: + _support_cpp = importlib.import_module("._support_cpp", __package__) + _avgstd_cpp = _support_cpp.avgstd + SUPPORT_CPP_AVAILABLE = True +except ImportError: + SUPPORT_CPP_AVAILABLE = False + +__all__ = ["avgstd"] + + +def avgstd(cube, sigma): + """Compute clipped mean, standard deviation and accepted-frame counts. + + Parameters + ---------- + cube : numpy.ndarray + Three-dimensional array of shape ``(nf, ny, nx)`` with ``float32`` data. + sigma : float + Rejection threshold for the iterative clipped-mean routine. + + Returns + ------- + tuple + ``(avg, std, num)`` arrays with shapes ``(ny, nx)``. + """ + + cube_arr = np.asarray(cube, dtype=np.float32) + if cube_arr.ndim != 3: + raise ValueError("cube must be a 3D array") + if sigma <= 1.0: + raise ValueError("sigma must be greater than 1") + + if SUPPORT_CPP_AVAILABLE: + return _avgstd_cpp(cube_arr, float(sigma)) + + return _avgstd_numpy(cube_arr, float(sigma)) + + +def _avgstd_numpy(cube_arr, sigma): + """Compute clipped mean, standard deviation and accepted-frame counts.""" + nf, ny, nx = cube_arr.shape + avg = np.empty((ny, nx), dtype=np.float32) + std = np.empty((ny, nx), dtype=np.float32) + num = np.empty((ny, nx), dtype=np.int32) + + for iy in range(ny): + for ix in range(nx): + vals = cube_arr[:, iy, ix].astype(np.float64) + ok = np.ones(nf, dtype=bool) + ncur = nf + + while True: + if ncur <= 0: + avg[iy, ix] = 0.0 + std[iy, ix] = 0.0 + num[iy, ix] = 0 + break + + tavg = vals[ok].mean() + if ncur > 1: + tstd = vals[ok].std(ddof=1) + else: + tstd = 0.0 + + thresh = sigma * tstd + new_ok = ok.copy() + new_ok[np.abs(vals - tavg) > thresh] = False + + if np.all(new_ok == ok): + kept = vals[ok] + avg[iy, ix] = float(kept.mean()) + std[iy, ix] = float(kept.std(ddof=1)) if kept.size > 1 else 0.0 + num[iy, ix] = int(kept.size) + break + + ok = new_ok + ncur = int(np.count_nonzero(ok)) + + return avg, std, num diff --git a/hipercam/support.pyx b/hipercam/support.pyx deleted file mode 100644 index 4cd06df2..00000000 --- a/hipercam/support.pyx +++ /dev/null @@ -1,223 +0,0 @@ -""" -Cython routines -""" - -import numpy as np -cimport numpy as np -cimport cython - -from libc.math cimport sqrt - -FTYPE = np.float32 -ctypedef np.float32_t FTYPE_t - -DTYPE = np.float64 -ctypedef np.float64_t DTYPE_t - -# Next two decorators reduce time from 108ms to 75ms -# in the test I carried out. -@cython.cdivision(True) -@cython.boundscheck(False) -def avgstd(np.ndarray[FTYPE_t, ndim=3] cube, float sigma): - """ - Given a 3D cube of dimensions (nf,ny,nx) which can be thought of as 'nf' - frames each of dimension (ny,nx), avgstd computes the mean, standard - deviation and number of frames used to calculate them for each pixel after - rejection at 'sigma' standard deviations across the first dimension. It - returns three 2D frames each of dimension (ny,nx). The rejection is - carried out in a cycle until no more rejection occurs. - - This routine is heavily cythonised and runs much faster than a pure-python - implementation, but that comes with a restriction on the data type of - 'cube' as listed below. - - Arguments:: - - cube : (3D numpy array, dtype=np.float32) - the set of frames to process. - - sigma : (float) - the rejection threshold. 3 to 4 goodish value. Must be > 1. - - Returns (avg,std,num) the average, standard deviation and number - of contributing frames. - """ - - # Dimensions and indices - cdef unsigned int nf = cube.shape[0] - cdef unsigned int ny = cube.shape[1] - cdef unsigned int nx = cube.shape[2] - cdef unsigned int ix, iy, iz - - # Sanity checks - assert (cube.dtype == FTYPE) and (sigma > 1.) - - # Output arrays: average, standard deviation and numbers of frames - cdef np.ndarray[FTYPE_t, ndim=2] avg = np.empty((ny,nx), dtype=FTYPE) - cdef np.ndarray[FTYPE_t, ndim=2] std = np.empty((ny,nx), dtype=FTYPE) - cdef np.ndarray[int, ndim=2] nfm = np.empty((ny,nx), dtype=np.int32) - - # Temporary arrays for storing values and rejection flags - cdef np.ndarray[FTYPE_t, ndim=1] vals = np.empty((nf,), dtype=FTYPE) - cdef np.ndarray[int, ndim=1] ok = np.empty((nf,), dtype=np.int32) - - # Other temporaries - cdef FTYPE_t tavg, tstd, thresh - cdef unsigned int num, numt - cdef double sum - - # The main loop starts. Anything inside this is done nx*ny - # times and better be fast... - for iy in xrange(ny): - for ix in xrange(nx): - - # extract values of given pixel for all frames, set flags to 1 - for iz in xrange(nf): - vals[iz] = cube[iz,iy,ix] - ok[iz] = 1 - - # calculate initial values of tavg, tstd and num - num = nf - sum = 0. - for iz in xrange(nf): - sum += vals[iz] - tavg = sum / num - sum = 0. - for iz in xrange(nf): - sum += (vals[iz]-tavg)**2 - if num > 1: - tstd = sqrt(sum/(num-1)) - else: - tstd = 0. - break - - # Now the rejection loop. Keep rejecting if - # number of rejections is positive - while True: - # pre-compute rejection threshold - thresh = sigma*tstd - sum = 0. - numt = 0 - for iz in xrange(nf): - if abs(vals[iz]-tavg) > thresh: - ok[iz] = 0 - - if ok[iz]: - sum += vals[iz] - numt += 1 - - tavg = sum / numt - sum = 0. - for iz in xrange(nf): - if ok[iz]: - sum += (vals[iz]-tavg)**2 - if numt > 1: - tstd = sqrt(sum/(numt-1)) - else: - tstd = 0. - break - - if numt == num: - break - num = numt - - # store the final numbers - avg[iy,ix] = tavg - std[iy,ix] = tstd - nfm[iy,ix] = num - - # return the frames - return (avg,std,nfm) - -@cython.cdivision(True) -@cython.boundscheck(False) -def gaussian( - np.ndarray[DTYPE_t, ndim=2] x, np.ndarray[DTYPE_t, ndim=2] y, - double sky, double height, double xcen, double ycen, double fwhm, - int xbin, int ybin, int ndiv): - """Returns a numpy array corresponding to the ordinate grids in xy set to a - symmetric 2D Gaussian plus a constant. The profile is essentially defined - by sky + height*exp(-alpha*r**2) where r is the distance from the centre - and alpha is set to give the desired FWHM, but account is taken of the - finite size of the pixels by summing over multiple points in each one. - - Arguments:: - - x : 2D numpy array - the X ordinates over which you want to compute the - profile. They should be measured in term of unbinned pixels. - - y : 2D numpy array - the Y ordinates over which you want to compute the - profile. They should be measured in term of unbinned pixels. - - sky : float - sky background - - height : float - height of central peak - - xcen : float - X-ordinate of centre - - ycen : float - Y-ordinate of centre - - fwhm : float - FWHM of the profile (unbinned pixels) - - xbin : int - X-size of pixels in terms of unbinned pixels, i.e. the binning factor in X - - ybin : int - Y-size of pixels in terms of unbinned pixels, i.e. the binning factor in Y - - ndiv : int - Parameter controlling treatment of sub-pixellation. If > 0, every - pixel will be sub-divided first into unbinned pixels, and then each - unbinned pixels will be split into a square array of ndiv by ndiv - points. The profile will be evaluated and averaged over each of these - points. Thus if the pixels are binned xbin by ybin, there will be - xbin*ybin*ndiv**2 evaluations per pixel. This is to cope with the - case where the seeing is becoming small compared to the pixels, but - obviously will slow things. To simply evaluate the profile once at - the centre of each pixel, set ndiv = 0. - - Returns:: 2D numpy array containing the Gaussian plus constant evaluated - on the ordinate grids in xy. - - """ - - # Dimensions and indices - cdef unsigned int ny = x.shape[0] - cdef unsigned int nx = x.shape[1] - cdef double alpha = 4.*np.log(2.)/fwhm**2 - cdef np.ndarray[DTYPE_t, ndim=2] rsq = np.empty((ny,nx), dtype=DTYPE) - cdef np.ndarray[DTYPE_t, ndim=2] prof = np.zeros((ny,nx), dtype=DTYPE) - cdef unsigned int ix, iy, isx, isy - cdef double xoff, yoff, xsoff, ysoff, soff - - if ndiv > 0: - # Complicated case with sub-pixellation allowed for - soff = (ndiv-1)/(2*ndiv) - for iy in range(ybin): - # loop over unbinned pixels in Y - yoff = iy-(ybin-1)/2 - soff - for ix in range(xbin): - # loop over unbinned pixels in X - xoff = ix-(xbin-1)/2 - soff - for isy in range(ndiv): - # loop over sub-pixels in y - ysoff = yoff + isy/ndiv - for isx in range(ndiv): - # loop over sub-pixels in x - xsoff = xoff + isx/ndiv - rsq = (x+xsoff-xcen)**2+(y+ysoff-ycen)**2 - prof += np.exp(-alpha*rsq) - - return sky+(height/xbin/ybin/ndiv**2)*prof - - else: - # Fast as possible, compute profile at pixel centres only - rsq = (x-xcen)**2+(y-ycen)**2 - return sky+height*np.exp(-alpha*rsq) diff --git a/hipercam/tests/fitting_test.py b/hipercam/tests/fitting_test.py new file mode 100644 index 00000000..0a957980 --- /dev/null +++ b/hipercam/tests/fitting_test.py @@ -0,0 +1,383 @@ +#!/usr/bin/env python3 +"""Test and benchmark the Numba and C++ fitting functions.""" + +import argparse +import time + +import numpy as np + +from hipercam import fitting + + +def _assert_allclose(name: str, a: np.ndarray, b: np.ndarray, rtol: float, atol: float) -> None: + """Assert that two arrays are close, raising an error with details if not.""" + if not np.allclose(a, b, rtol=rtol, atol=atol, equal_nan=True): + diff = np.max(np.abs(a - b)) + rel = np.max(np.abs(a - b) / np.maximum(np.abs(b), 1e-300)) + raise AssertionError( + f"{name} mismatch: max_abs={diff:.3e}, max_rel={rel:.3e}, " + f"rtol={rtol:.1e}, atol={atol:.1e}" + ) + + +def _make_case(rng: np.random.Generator) -> dict[str, object]: + """Make a single random test case for fitting functions.""" + # Make a simple 2D grid of x and y coordinates + ny = int(rng.integers(8, 24)) + nx = int(rng.integers(8, 24)) + x1 = np.linspace(100.25, 100.25 + 0.6 * (nx - 1), nx) + y1 = np.linspace(200.75, 200.75 + 0.7 * (ny - 1), ny) + x, y = np.meshgrid(x1, y1) + + return { + "x": x, + "y": y, + "sky": float(rng.uniform(-100.0, 500.0)), + "height": float(rng.uniform(1.0, 5e4)), + "xcen": float(rng.uniform(x.min() - 1.0, x.max() + 1.0)), + "ycen": float(rng.uniform(y.min() - 1.0, y.max() + 1.0)), + "fwhm": float(rng.uniform(0.5, 8.0)), + "beta": float(rng.uniform(0.05, 8.0)), + "xbin": int(rng.integers(1, 5)), + "ybin": int(rng.integers(1, 5)), + "ndiv": int(rng.integers(0, 4)), + } + + +def _make_cases(ncases: int, seed: int) -> list[dict[str, object]]: + """Make a list of random test cases for fitting functions.""" + rng = np.random.default_rng(seed) + return [_make_case(rng) for _ in range(ncases)] + + +def run_comparison(case_id: int, case: dict[str, object], rtol: float, atol: float) -> None: + """Run a single test case, comparing C++ and Numba implementations.""" + if not fitting.FITTING_CCP_AVAILABLE: + raise RuntimeError("hipercam.fitting_cpp is not available. Build/install extension first.") + + # Moffat profile + results_numba = fitting._moffat_numba( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['beta'], + case['xbin'], + case['ybin'], + case['ndiv'] + ) + results_cpp = fitting._moffat_cpp( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['beta'], + case['xbin'], + case['ybin'], + case['ndiv'] + ) + _assert_allclose(f"case {case_id} moffat", results_numba, results_cpp, rtol, atol) + + # Moffat derivatives (all flag combinations) + for comp_dfwhm in (False, True): + for comp_dbeta in (False, True): + results_numba = fitting._dmoffat_numba( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['beta'], + case['xbin'], + case['ybin'], + case['ndiv'], + comp_dfwhm, + comp_dbeta, + ) + results_cpp = fitting._dmoffat_cpp( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['beta'], + case['xbin'], + case['ybin'], + case['ndiv'], + comp_dfwhm, + comp_dbeta, + ) + _assert_allclose( + f"case {case_id} dmoffat[dfwhm={comp_dfwhm},dbeta={comp_dbeta}]", + results_numba, + results_cpp, + rtol, + atol, + ) + + # Gaussian profile + results_numba = fitting._gaussian_numba( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['xbin'], + case['ybin'], + case['ndiv'] + ) + results_cpp = fitting._gaussian_cpp( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['xbin'], + case['ybin'], + case['ndiv'] + ) + _assert_allclose(f"case {case_id} gaussian", results_numba, results_cpp, rtol, atol) + + # Gaussian derivatives + for comp_dfwhm in (False, True): + results_numba = fitting._dgaussian_numba( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['xbin'], + case['ybin'], + case['ndiv'], + comp_dfwhm + ) + results_cpp = fitting._dgaussian_cpp( + case['x'], + case['y'], + case['sky'], + case['height'], + case['xcen'], + case['ycen'], + case['fwhm'], + case['xbin'], + case['ybin'], + case['ndiv'], + comp_dfwhm + ) + _assert_allclose( + f"case {case_id} dgaussian[dfwhm={comp_dfwhm}]", + results_numba, + results_cpp, + rtol, + atol, + ) + + +def _time_function( + callable_fn, cases: list[dict[str, object]], nruns: int, *args, **kwargs +) -> np.ndarray: + """Time a function over multiple runs with given cases and return the elapsed times.""" + times = [] + for _ in range(nruns): + start = time.perf_counter() + for case in cases: + callable_fn(case, *args, **kwargs) + times.append(time.perf_counter() - start) + return [time / len(cases) for time in times] # Return time per case + + +def _compare_functions(function_numba, function_cpp, cases, nruns, *args, **kwargs): + """Compare the timing of Numba and C++ functions over multiple runs.""" + times_numba = _time_function(function_numba, cases, nruns, *args, **kwargs) + times_cpp = _time_function(function_cpp, cases, nruns, *args, **kwargs) + print( + f"{np.mean(times_numba) * 1000:5.3f}±{np.std(times_numba) * 1000:5.3f} ms " + f"{np.mean(times_cpp) * 1000:5.3f}±{np.std(times_cpp) * 1000:5.3f} ms " + f"{np.mean(times_numba) / np.mean(times_cpp):5.2f}x" + ) + + +def run_benchmarks(cases: list[dict[str, object]], nruns: int) -> None: + """Run benchmarks comparing C++ and Numba implementations of fitting kernels.""" + print(f"\nBenchmarking {len(cases)} cases over {nruns} runs:") + print("-" * 85) + print(f"{'function':32s} {'Numba':18s} {'C++':18s} {'Speedup':8s}") + print("-" * 85) + + # Run the comparisons + # Moffat kernel + def moffat_numba(case): + return fitting._moffat_numba( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["beta"], + case["xbin"], + case["ybin"], + case["ndiv"], + ) + def moffat_cpp(case): + return fitting._moffat_cpp( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["beta"], + case["xbin"], + case["ybin"], + case["ndiv"], + ) + print(f"{'moffat':32s}", end=" ") + _compare_functions(moffat_numba, moffat_cpp, cases, nruns) + + # Gaussian kernel + def gaussian_numba(case): + return fitting._gaussian_numba( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["xbin"], + case["ybin"], + case["ndiv"], + ) + def gaussian_cpp(case): + return fitting._gaussian_cpp( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["xbin"], + case["ybin"], + case["ndiv"], + ) + print(f"{'gaussian':32s}", end=" ") + _compare_functions(gaussian_numba, gaussian_cpp, cases, nruns) + + # Moffat derivatives + def dmoffat_numba(case, comp_dfwhm, comp_dbeta): + return fitting._dmoffat_numba( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["beta"], + case["xbin"], + case["ybin"], + case["ndiv"], + comp_dfwhm, + comp_dbeta, + ) + def dmoffat_cpp(case, comp_dfwhm, comp_dbeta): + return fitting._dmoffat_cpp( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["beta"], + case["xbin"], + case["ybin"], + case["ndiv"], + comp_dfwhm, + comp_dbeta, + ) + for comp_dfwhm in (False, True): + for comp_dbeta in (False, True): + print(f"{f'dmoffat[dfwhm={comp_dfwhm},dbeta={comp_dbeta}]':32s}", end=" ") + _compare_functions(dmoffat_numba, dmoffat_cpp, cases, nruns, comp_dfwhm, comp_dbeta) + + # Gaussian derivatives + def dgaussian_numba(case, comp_dfwhm): + return fitting._dgaussian_numba( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["xbin"], + case["ybin"], + case["ndiv"], + comp_dfwhm, + ) + def dgaussian_cpp(case, comp_dfwhm): + return fitting._dgaussian_cpp( + case["x"], + case["y"], + case["sky"], + case["height"], + case["xcen"], + case["ycen"], + case["fwhm"], + case["xbin"], + case["ybin"], + case["ndiv"], + comp_dfwhm, + ) + for comp_dfwhm in (False, True): + print(f"{f'dgaussian[dfwhm={comp_dfwhm}]':32s}", end=" ") + _compare_functions(dgaussian_numba, dgaussian_cpp, cases, nruns, comp_dfwhm) + + +def main() -> int: + parser = argparse.ArgumentParser( + description="Compare fitting timing between the Numpy and C++ implementations" + ) + parser.add_argument("--repeats", type=int, default=100) + parser.add_argument("--cases", type=int, default=20) + parser.add_argument("--rtol", type=float, default=1e-9) + parser.add_argument("--atol", type=float, default=1e-11) + parser.add_argument("--seed", type=int, default=12345) + args = parser.parse_args() + + # Make the random test cases + cases = _make_cases(args.cases, args.seed) + + # Run the equivalence checks for each case + for case_id, case in enumerate(cases, start=1): + run_comparison(case_id, case, args.rtol, args.atol) + print(f"All checks passed: {args.cases} cases, rtol={args.rtol:.1e}, atol={args.atol:.1e}") + + # Run benchmarks (using a new set of cases with a different seed) + bench_cases = _make_cases(args.cases, args.seed + 1) + run_benchmarks(bench_cases, args.repeats) + + +if __name__ == "__main__": + main() diff --git a/hipercam/tests/support_test.py b/hipercam/tests/support_test.py new file mode 100644 index 00000000..fb3037cc --- /dev/null +++ b/hipercam/tests/support_test.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +"""Test and benchmark the Numpy and C++ avgstd implementations.""" + +import argparse +import time + +import numpy as np + +from hipercam import support + + +def main() -> None: + parser = argparse.ArgumentParser( + description="Compare avgstd timing between the Numpy and C++ implementations" + ) + parser.add_argument("--repeats", type=int, default=10) + parser.add_argument("--frames", type=int, default=32) + parser.add_argument("--height", type=int, default=64) + parser.add_argument("--width", type=int, default=64) + parser.add_argument("--sigma", type=float, default=3.0) + parser.add_argument("--seed", type=int, default=12345, help="RNG seed") + args = parser.parse_args() + + # Build a random data cube for testing + rng = np.random.default_rng(seed=args.seed) + cube = rng.normal(loc=0.0, scale=1.0, size=(args.frames, args.height, args.width)) + cube = cube.astype(np.float32) + print(f"Cube shape: {cube.shape}") + + # Call the functions (this also ensures the C++ extension is loaded and compiled for timing) + numpy_avg, numpy_std, numpy_num = support._avgstd_numpy(cube, args.sigma) + cpp_avg, cpp_std, cpp_num = support._avgstd_cpp(cube, args.sigma) + + # Raise an error if the results differ + if not ( + np.allclose(cpp_avg, numpy_avg) + and np.allclose(cpp_std, numpy_std) + and np.array_equal(cpp_num, numpy_num) + ): + raise RuntimeError("C++ and Numpy avgstd implementations produced different results") + else: + print("C++ and Numpy avgstd implementations produced the same results") + + # Now time the functions and print the results + numpy_times = [] + for _ in range(args.repeats): + start = time.perf_counter() + support._avgstd_numpy(cube, args.sigma) + numpy_times.append(time.perf_counter() - start) + + cpp_times = [] + for _ in range(args.repeats): + start = time.perf_counter() + support._avgstd_cpp(cube, args.sigma) + cpp_times.append(time.perf_counter() - start) + + numpy_time = float(np.mean(numpy_times)) + cpp_time = float(np.mean(cpp_times)) + print(f"Numpy avgstd mean time over {args.repeats} runs: {numpy_time * 1000:.2f} ms") + print(f"C++ avgstd mean time over {args.repeats} runs: {cpp_time * 1000:.2f} ms") + print(f"Speedup: {numpy_time / cpp_time:.2f}x") + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index c230cae5..947ec0a7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=61.0", "setuptools-scm", "wheel", "Cython", "numpy"] +requires = ["setuptools>=61.0", "setuptools-scm", "wheel", "numpy", "pybind11>=2.6.0"] build-backend = "setuptools.build_meta" [project] @@ -19,16 +19,16 @@ classifiers = [ "Programming Language :: Python :: 3", "Programming Language :: Python :: 3.6", ] -requires-python = ">=3.6" +requires-python = ">=3.9" dependencies = [ "astropy", - "Cython", "fitsio", "keyring", "matplotlib", "numba", "numpy", "pandas", + "pybind11>=2.6.0", "requests", "sep>=1.4", "trm.cline", diff --git a/setup.py b/setup.py index 38f2ec2f..76baf5bc 100755 --- a/setup.py +++ b/setup.py @@ -1,25 +1,54 @@ """ -Minimal setup.py for Cython extension support. +Minimal setup.py for pybind11 extension support. All other metadata is in pyproject.toml. """ import os +import sys + import numpy as np -from Cython.Build import cythonize from setuptools import setup -from setuptools.extension import Extension +import pybind11 +from pybind11.setup_helpers import Pybind11Extension, build_ext + +# pybind11 extension for profile fitting +use_openmp = os.environ.get("HIPERCAM_USE_OPENMP", "1") != "0" +openmp_args = [] +openmp_link_args = [] +if use_openmp and sys.platform.startswith("linux"): + openmp_args = ["-fopenmp"] + openmp_link_args = ["-fopenmp"] -# cython support routine -extension = [ - Extension( - "hipercam.support", - [os.path.join("hipercam", "support.pyx")], - libraries=["m"], - include_dirs=[np.get_include()], - extra_compile_args=["-fno-strict-aliasing"], +pybind11_extensions = [ + Pybind11Extension( + "hipercam._fitting_cpp", + ["hipercam/fitting.cpp"], + include_dirs=[ + np.get_include(), + pybind11.get_include(), + pybind11.get_include(user=True), + ], + language="c++", + extra_compile_args=["-std=c++11", "-O3", "-ffast-math", "-march=native"] + + openmp_args, + extra_link_args=openmp_link_args, + ), + Pybind11Extension( + "hipercam._support_cpp", + ["hipercam/support.cpp"], + include_dirs=[ + np.get_include(), + pybind11.get_include(), + pybind11.get_include(user=True), + ], + language="c++", + extra_compile_args=["-std=c++11", "-O3", "-ffast-math", "-march=native"] + + openmp_args, + extra_link_args=openmp_link_args, ), ] setup( - ext_modules=cythonize(extension), + ext_modules=pybind11_extensions, + cmdclass={"build_ext": build_ext}, )