From 3200749fa6c714db982e1f98dec5a0ff97b6cff4 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Fri, 13 Jun 2025 11:50:41 +0100 Subject: [PATCH 01/11] Bug fixing --- hipercam/fitting.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/hipercam/fitting.py b/hipercam/fitting.py index 285c4466..09a50b62 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -405,7 +405,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) @@ -878,7 +878,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 == "": @@ -1167,7 +1167,7 @@ def fitGaussian( ) if not res.success: raise HipercamError(res.message) - nfev += nfev + nfev += res.nfev # get Jacobian J = np.matrix(res.jac) @@ -1668,4 +1668,3 @@ def model(self, param): self.ybin, self.ndiv, ) - From c5913dd01cd0897b6ee3255d9af63e53fbe5d2f2 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Wed, 29 Oct 2025 15:50:07 +0000 Subject: [PATCH 02/11] Optimisations in fit classes --- hipercam/fitting.py | 256 +++++++++++++++++++++++--------------------- 1 file changed, 136 insertions(+), 120 deletions(-) diff --git a/hipercam/fitting.py b/hipercam/fitting.py index 09a50b62..3a1eacbb 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -757,7 +757,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 +788,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 +795,39 @@ 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.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 +875,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": @@ -924,42 +941,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 +970,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 +992,34 @@ 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. + """ + 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. + """ + 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] + + ########################################## # @@ -1468,6 +1477,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 +1513,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 +1520,31 @@ 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.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 +1610,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 +1638,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. - - Argument:: + Returns list of 2D arrays of the partial derivatives of + the model with respect to the variable parameters. - 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,4 +1659,28 @@ 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. + """ + 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. + """ + 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] From 5c6bcb0d84194772962d925bd4dfa21de38c190b Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Thu, 14 May 2026 16:18:05 +0100 Subject: [PATCH 03/11] Add C++ versions of fitting functions --- hipercam/fitting.py | 99 ++++++- hipercam/fitting_cpp.cpp | 564 +++++++++++++++++++++++++++++++++++++++ pyproject.toml | 5 +- setup.py | 27 +- 4 files changed, 684 insertions(+), 11 deletions(-) create mode 100644 hipercam/fitting_cpp.cpp diff --git a/hipercam/fitting.py b/hipercam/fitting.py index 3a1eacbb..5e6ef9c6 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -9,7 +9,11 @@ from scipy.optimize import least_squares from .core import * from .window import * -from . import support + +try: + from . import fitting_cpp +except ImportError: + fitting_cpp = None __all__ = ("combFit", "fitMoffat", "fitGaussian", "moffat", "gaussian") @@ -451,7 +455,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 +464,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 +513,22 @@ 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_cpp is not None: + return fitting_cpp.moffat( + x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv + ) + else: + 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 +564,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 +575,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 +640,44 @@ def dmoffat( numba just-in-time compiler function better. """ + if fitting_cpp is not None: + return fitting_cpp.dmoffat( + x, + y, + sky, + height, + xcen, + ycen, + fwhm, + beta, + xbin, + ybin, + ndiv, + comp_dfwhm, + comp_dbeta, + ) + else: + 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 @@ -1262,7 +1322,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 @@ -1270,6 +1330,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 @@ -1316,6 +1379,16 @@ def gaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): on the ordinate grids in xy. """ + if fitting_cpp is not None: + return fitting_cpp.gaussian( + x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv + ) + else: + 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 @@ -1349,7 +1422,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 @@ -1357,6 +1429,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 @@ -1409,6 +1484,18 @@ def dgaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) appear in the function call. """ + if fitting_cpp is not None: + return fitting_cpp.dgaussian( + x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm + ) + else: + 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) diff --git a/hipercam/fitting_cpp.cpp b/hipercam/fitting_cpp.cpp new file mode 100644 index 00000000..ed38becf --- /dev/null +++ b/hipercam/fitting_cpp.cpp @@ -0,0 +1,564 @@ +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +// Helper function to calculate the Moffat alpha parameter +inline 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); +} + +/** + * C++ implementation of the Moffat profile function + */ +py::array_t moffat_cpp(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 *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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + double norm = height / xbin / ybin / (ndiv * ndiv); + + if (ndiv > 0) { + // With sub-pixellation + std::fill_n(result_ptr, n_pixels, 0.0); + + // Mean offset within sub-pixels + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + // Loop over all pixels + for (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; + + // Loop over sub-pixels + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy / static_cast(ndiv); + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx / static_cast(ndiv); + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + prof += (height / xbin / ybin / (ndiv * ndiv)) * + std::pow(1.0 + alpha * rsq, -tbeta); + } + } + } + } + + result_ptr[pixel_idx] = sky + prof; + } + } else { + // Fast calculation at pixel centers + for (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] = height * std::pow(1.0 + alpha * rsq, -tbeta) + sky; + } + } + + return result; +} + +/** + * C++ implementation of the Moffat profile derivatives + */ +std::vector> +dmoffat_cpp(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) { + + // 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; + + // 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 *dsky_ptr = static_cast(dsky_info.ptr); + double *dheight_ptr = static_cast(dheight_info.ptr); + double *dxcen_ptr = static_cast(dxcen_info.ptr); + double *dycen_ptr = static_cast(dycen_info.ptr); + + py::buffer_info dfwhm_info; + py::buffer_info dbeta_info; + double *dfwhm_ptr = nullptr; + double *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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + // 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 + std::fill_n(dheight_ptr, n_pixels, 0.0); + std::fill_n(dxcen_ptr, n_pixels, 0.0); + std::fill_n(dycen_ptr, n_pixels, 0.0); + + if (comp_dfwhm) { + std::fill_n(dfwhm_ptr, n_pixels, 0.0); + } + + if (comp_dbeta) { + std::fill_n(dbeta_ptr, n_pixels, 0.0); + } + + // Mean offset within sub-pixels + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + // Loop over all pixels + for (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]; + + // Loop over sub-pixels + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy / static_cast(ndiv); + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx / static_cast(ndiv); + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + + double denom = 1.0 + alpha * rsq; + double save1 = height * std::pow(denom, -tbeta - 1.0); + double save2 = save1 * rsq; + + // Derivatives + double dh = std::pow(denom, -tbeta); + dheight_ptr[pixel_idx] += dh; + dxcen_ptr[pixel_idx] += (2.0 * alpha * tbeta) * dx * save1; + dycen_ptr[pixel_idx] += (2.0 * alpha * tbeta) * dy * save1; + + if (comp_dfwhm) { + dfwhm_ptr[pixel_idx] += (2.0 * alpha * tbeta / fwhm) * save2; + } + + if (comp_dbeta) { + double log_denom = std::log(denom); + dbeta_ptr[pixel_idx] += + (-log_denom * height * dh + + (4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / tbeta / + (fwhm * fwhm)) * + save2); + } + } + } + } + } + } + + // Normalize by number of evaluations + double nadd = xbin * ybin * ndiv * ndiv; + for (size_t i = 0; i < n_pixels; ++i) { + dheight_ptr[i] /= nadd; + dxcen_ptr[i] /= nadd; + dycen_ptr[i] /= nadd; + + if (comp_dfwhm) { + dfwhm_ptr[i] /= nadd; + } + + if (comp_dbeta) { + dbeta_ptr[i] /= nadd; + } + } + } else { + // Fast calculation at pixel centers + for (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; + double save1 = height * std::pow(denom, -tbeta - 1.0); + double save2 = save1 * rsq; + + // Derivatives + dheight_ptr[i] = std::pow(denom, -tbeta); + dxcen_ptr[i] = (2.0 * alpha * tbeta) * dx * save1; + dycen_ptr[i] = (2.0 * alpha * tbeta) * dy * save1; + + if (comp_dfwhm) { + dfwhm_ptr[i] = (2.0 * alpha * tbeta / fwhm) * save2; + } + + if (comp_dbeta) { + double log_denom = std::log(denom); + dbeta_ptr[i] = (-log_denom * height * dheight_ptr[i] + + (4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / + tbeta / (fwhm * fwhm)) * + 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_cpp(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 *result_ptr = static_cast(result_info.ptr); + + // Calculate Gaussian profile parameter + double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + if (ndiv > 0) { + // With sub-pixellation + std::fill_n(result_ptr, n_pixels, 0.0); + + // Mean offset within sub-pixels + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + // Loop over all pixels + for (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; + + // Loop over sub-pixels + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy / static_cast(ndiv); + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx / static_cast(ndiv); + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + prof += std::exp(-alpha * rsq); + } + } + } + } + + result_ptr[pixel_idx] = + sky + (height / xbin / ybin / (ndiv * ndiv)) * prof; + } + } else { + // Fast calculation at pixel centers + for (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_cpp(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) { + + // 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; + + // 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 *dsky_ptr = static_cast(dsky_info.ptr); + double *dheight_ptr = static_cast(dheight_info.ptr); + double *dxcen_ptr = static_cast(dxcen_info.ptr); + double *dycen_ptr = static_cast(dycen_info.ptr); + + py::buffer_info dfwhm_info; + double *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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + + // 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 + std::fill_n(dheight_ptr, n_pixels, 0.0); + std::fill_n(dxcen_ptr, n_pixels, 0.0); + std::fill_n(dycen_ptr, n_pixels, 0.0); + + if (comp_dfwhm) { + std::fill_n(dfwhm_ptr, n_pixels, 0.0); + } + + // Mean offset within sub-pixels + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + // Loop over all pixels + for (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]; + + // Loop over sub-pixels + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy / static_cast(ndiv); + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx / static_cast(ndiv); + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + + // Gaussian value + double dh = std::exp(-alpha * rsq); + dheight_ptr[pixel_idx] += dh; + dxcen_ptr[pixel_idx] += (2.0 * alpha * height) * dh * dx; + dycen_ptr[pixel_idx] += (2.0 * alpha * height) * dh * dy; + + if (comp_dfwhm) { + dfwhm_ptr[pixel_idx] += + (2.0 * alpha * height / fwhm) * dh * rsq; + } + } + } + } + } + } + + // Normalize by number of evaluations + double nadd = xbin * ybin * ndiv * ndiv; + for (size_t i = 0; i < n_pixels; ++i) { + dheight_ptr[i] /= nadd; + dxcen_ptr[i] /= nadd; + dycen_ptr[i] /= nadd; + + if (comp_dfwhm) { + dfwhm_ptr[i] /= nadd; + } + } + } else { + // Fast calculation at pixel centers + for (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] = (2.0 * alpha * height) * dh * dx; + dycen_ptr[i] = (2.0 * alpha * height) * dh * dy; + + if (comp_dfwhm) { + dfwhm_ptr[i] = (2.0 * alpha * height / fwhm) * 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_cpp, "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_cpp, + "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("gaussian", &gaussian_cpp, "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_cpp, + "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")); +} diff --git a/pyproject.toml b/pyproject.toml index c230cae5..a1d6f20e 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", "Cython", "numpy", "pybind11>=2.6.0"] build-backend = "setuptools.build_meta" [project] @@ -19,7 +19,7 @@ classifiers = [ "Programming Language :: Python :: 3", "Programming Language :: Python :: 3.6", ] -requires-python = ">=3.6" +requires-python = ">=3.9" dependencies = [ "astropy", "Cython", @@ -29,6 +29,7 @@ dependencies = [ "numba", "numpy", "pandas", + "pybind11>=2.6.0", "requests", "sep>=1.4", "trm.cline", diff --git a/setup.py b/setup.py index 38f2ec2f..327d3196 100755 --- a/setup.py +++ b/setup.py @@ -1,16 +1,21 @@ """ -Minimal setup.py for Cython extension support. +Minimal setup.py for Cython and pybind11 extension support. All other metadata is in pyproject.toml. """ import os + +# need for Cython and pybind11 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 + # cython support routine -extension = [ +cython_extensions = [ Extension( "hipercam.support", [os.path.join("hipercam", "support.pyx")], @@ -20,6 +25,22 @@ ), ] +# pybind11 extension for profile fitting +pybind11_extensions = [ + Pybind11Extension( + "hipercam.fitting_cpp", + ["hipercam/fitting_cpp.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"], + ), +] + setup( - ext_modules=cythonize(extension), + ext_modules=cythonize(cython_extensions) + pybind11_extensions, + cmdclass={"build_ext": build_ext}, ) From 723335c1d061afbe7b52a9afea898ac7c6b28982 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Thu, 14 May 2026 16:53:53 +0100 Subject: [PATCH 04/11] Some optimisations --- hipercam/fitting_cpp.cpp | 134 ++++++++++++++++++--------------------- 1 file changed, 61 insertions(+), 73 deletions(-) diff --git a/hipercam/fitting_cpp.cpp b/hipercam/fitting_cpp.cpp index ed38becf..09d9a299 100644 --- a/hipercam/fitting_cpp.cpp +++ b/hipercam/fitting_cpp.cpp @@ -45,11 +45,11 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - double norm = height / xbin / ybin / (ndiv * ndiv); - if (ndiv > 0) { // With sub-pixellation std::fill_n(result_ptr, n_pixels, 0.0); + double norm = height / xbin / ybin / (ndiv * ndiv); + double inv_ndiv = 1.0 / static_cast(ndiv); // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); @@ -66,14 +66,13 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, for (int ix = 0; ix < xbin; ++ix) { double xoff = ix - (xbin - 1) / 2.0 - soff; for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy / static_cast(ndiv); + double ysoff = yoff + isy * inv_ndiv; for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx / static_cast(ndiv); + double xsoff = xoff + isx * inv_ndiv; double dx = x_val + xsoff - xcen; double dy = y_val + ysoff - ycen; double rsq = dx * dx + dy * dy; - prof += (height / xbin / ybin / (ndiv * ndiv)) * - std::pow(1.0 + alpha * rsq, -tbeta); + prof += norm * std::pow(1.0 + alpha * rsq, -tbeta); } } } @@ -117,6 +116,7 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, // 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); @@ -165,6 +165,10 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, // 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); size_t n_pixels = x_info.shape[0] * x_info.shape[1]; @@ -173,17 +177,8 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, if (ndiv > 0) { // With sub-pixellation - std::fill_n(dheight_ptr, n_pixels, 0.0); - std::fill_n(dxcen_ptr, n_pixels, 0.0); - std::fill_n(dycen_ptr, n_pixels, 0.0); - - if (comp_dfwhm) { - std::fill_n(dfwhm_ptr, n_pixels, 0.0); - } - - if (comp_dbeta) { - std::fill_n(dbeta_ptr, n_pixels, 0.0); - } + double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); + double inv_ndiv = 1.0 / static_cast(ndiv); // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); @@ -192,6 +187,11 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, for (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 dheight = 0.0; + double dxcen = 0.0; + double dycen = 0.0; + double dfwhm = 0.0; + double dbeta = 0.0; // Loop over sub-pixels for (int iy = 0; iy < ybin; ++iy) { @@ -199,9 +199,9 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, for (int ix = 0; ix < xbin; ++ix) { double xoff = ix - (xbin - 1) / 2.0 - soff; for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy / static_cast(ndiv); + double ysoff = yoff + isy * inv_ndiv; for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx / static_cast(ndiv); + double xsoff = xoff + isx * inv_ndiv; double dx = x_val + xsoff - xcen; double dy = y_val + ysoff - ycen; double rsq = dx * dx + dy * dy; @@ -212,41 +212,33 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, // Derivatives double dh = std::pow(denom, -tbeta); - dheight_ptr[pixel_idx] += dh; - dxcen_ptr[pixel_idx] += (2.0 * alpha * tbeta) * dx * save1; - dycen_ptr[pixel_idx] += (2.0 * alpha * tbeta) * dy * save1; + dheight += dh; + dxcen += two_alpha_tbeta * dx * save1; + dycen += two_alpha_tbeta * dy * save1; if (comp_dfwhm) { - dfwhm_ptr[pixel_idx] += (2.0 * alpha * tbeta / fwhm) * save2; + dfwhm += dfwhm_coeff * save2; } if (comp_dbeta) { double log_denom = std::log(denom); - dbeta_ptr[pixel_idx] += - (-log_denom * height * dh + - (4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / tbeta / - (fwhm * fwhm)) * - save2); + dbeta += (-log_denom * height * dh + dbeta_coeff * save2); } } } } } - } - // Normalize by number of evaluations - double nadd = xbin * ybin * ndiv * ndiv; - for (size_t i = 0; i < n_pixels; ++i) { - dheight_ptr[i] /= nadd; - dxcen_ptr[i] /= nadd; - dycen_ptr[i] /= nadd; + dheight_ptr[pixel_idx] = dheight * inv_nadd; + dxcen_ptr[pixel_idx] = dxcen * inv_nadd; + dycen_ptr[pixel_idx] = dycen * inv_nadd; if (comp_dfwhm) { - dfwhm_ptr[i] /= nadd; + dfwhm_ptr[pixel_idx] = dfwhm * inv_nadd; } if (comp_dbeta) { - dbeta_ptr[i] /= nadd; + dbeta_ptr[pixel_idx] = dbeta * inv_nadd; } } } else { @@ -262,19 +254,17 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, // Derivatives dheight_ptr[i] = std::pow(denom, -tbeta); - dxcen_ptr[i] = (2.0 * alpha * tbeta) * dx * save1; - dycen_ptr[i] = (2.0 * alpha * tbeta) * dy * save1; + dxcen_ptr[i] = two_alpha_tbeta * dx * save1; + dycen_ptr[i] = two_alpha_tbeta * dy * save1; if (comp_dfwhm) { - dfwhm_ptr[i] = (2.0 * alpha * tbeta / fwhm) * save2; + dfwhm_ptr[i] = dfwhm_coeff * save2; } if (comp_dbeta) { double log_denom = std::log(denom); - dbeta_ptr[i] = (-log_denom * height * dheight_ptr[i] + - (4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / - tbeta / (fwhm * fwhm)) * - save2); + dbeta_ptr[i] = + (-log_denom * height * dheight_ptr[i] + dbeta_coeff * save2); } } } @@ -336,6 +326,8 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, if (ndiv > 0) { // With sub-pixellation std::fill_n(result_ptr, n_pixels, 0.0); + double norm = height / xbin / ybin / (ndiv * ndiv); + double inv_ndiv = 1.0 / static_cast(ndiv); // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); @@ -352,9 +344,9 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, for (int ix = 0; ix < xbin; ++ix) { double xoff = ix - (xbin - 1) / 2.0 - soff; for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy / static_cast(ndiv); + double ysoff = yoff + isy * inv_ndiv; for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx / static_cast(ndiv); + double xsoff = xoff + isx * inv_ndiv; double dx = x_val + xsoff - xcen; double dy = y_val + ysoff - ycen; double rsq = dx * dx + dy * dy; @@ -364,8 +356,7 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, } } - result_ptr[pixel_idx] = - sky + (height / xbin / ybin / (ndiv * ndiv)) * prof; + result_ptr[pixel_idx] = sky + norm * prof; } } else { // Fast calculation at pixel centers @@ -403,6 +394,7 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, // 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); @@ -438,6 +430,8 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, // 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; size_t n_pixels = x_info.shape[0] * x_info.shape[1]; @@ -446,13 +440,8 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, if (ndiv > 0) { // With sub-pixellation - std::fill_n(dheight_ptr, n_pixels, 0.0); - std::fill_n(dxcen_ptr, n_pixels, 0.0); - std::fill_n(dycen_ptr, n_pixels, 0.0); - - if (comp_dfwhm) { - std::fill_n(dfwhm_ptr, n_pixels, 0.0); - } + double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); + double inv_ndiv = 1.0 / static_cast(ndiv); // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); @@ -461,6 +450,10 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, for (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 dheight = 0.0; + double dxcen = 0.0; + double dycen = 0.0; + double dfwhm = 0.0; // Loop over sub-pixels for (int iy = 0; iy < ybin; ++iy) { @@ -468,38 +461,33 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, for (int ix = 0; ix < xbin; ++ix) { double xoff = ix - (xbin - 1) / 2.0 - soff; for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy / static_cast(ndiv); + double ysoff = yoff + isy * inv_ndiv; for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx / static_cast(ndiv); + double xsoff = xoff + isx * inv_ndiv; double dx = x_val + xsoff - xcen; double dy = y_val + ysoff - ycen; double rsq = dx * dx + dy * dy; // Gaussian value double dh = std::exp(-alpha * rsq); - dheight_ptr[pixel_idx] += dh; - dxcen_ptr[pixel_idx] += (2.0 * alpha * height) * dh * dx; - dycen_ptr[pixel_idx] += (2.0 * alpha * height) * dh * dy; + dheight += dh; + dxcen += two_alpha_height * dh * dx; + dycen += two_alpha_height * dh * dy; if (comp_dfwhm) { - dfwhm_ptr[pixel_idx] += - (2.0 * alpha * height / fwhm) * dh * rsq; + dfwhm += dfwhm_coeff * dh * rsq; } } } } } - } - // Normalize by number of evaluations - double nadd = xbin * ybin * ndiv * ndiv; - for (size_t i = 0; i < n_pixels; ++i) { - dheight_ptr[i] /= nadd; - dxcen_ptr[i] /= nadd; - dycen_ptr[i] /= nadd; + dheight_ptr[pixel_idx] = dheight * inv_nadd; + dxcen_ptr[pixel_idx] = dxcen * inv_nadd; + dycen_ptr[pixel_idx] = dycen * inv_nadd; if (comp_dfwhm) { - dfwhm_ptr[i] /= nadd; + dfwhm_ptr[pixel_idx] = dfwhm * inv_nadd; } } } else { @@ -512,11 +500,11 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, // Gaussian value double dh = std::exp(-alpha * rsq); dheight_ptr[i] = dh; - dxcen_ptr[i] = (2.0 * alpha * height) * dh * dx; - dycen_ptr[i] = (2.0 * alpha * height) * dh * dy; + dxcen_ptr[i] = two_alpha_height * dh * dx; + dycen_ptr[i] = two_alpha_height * dh * dy; if (comp_dfwhm) { - dfwhm_ptr[i] = (2.0 * alpha * height / fwhm) * dh * rsq; + dfwhm_ptr[i] = dfwhm_coeff * dh * rsq; } } } From 215ce1d77ab95cf00e24dceae9b5edbdb3e7135a Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Thu, 14 May 2026 17:30:30 +0100 Subject: [PATCH 05/11] Optimise residual/jacobian in C++ --- hipercam/fitting.py | 80 ++++++ hipercam/fitting_cpp.cpp | 532 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 612 insertions(+) diff --git a/hipercam/fitting.py b/hipercam/fitting.py index 5e6ef9c6..87a815dd 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -858,6 +858,8 @@ def __init__(self, wind, sigma, ndiv, mode, fwhm=None, beta=None): 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): @@ -1063,6 +1065,25 @@ def fun(self, param): Used by scipy.optimize.least_squares. """ + if fitting_cpp is not None and hasattr(fitting_cpp, "moffat_resid"): + sky, height, xcen, ycen, fwhm, beta = self.get_par(param) + return fitting_cpp.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() @@ -1075,6 +1096,27 @@ def jac(self, param): Used by scipy.optimize.least_squares. """ + if fitting_cpp is not None and hasattr(fitting_cpp, "dmoffat_jac"): + sky, height, xcen, ycen, fwhm, beta = self.get_par(param) + return fitting_cpp.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] @@ -1610,6 +1652,7 @@ def __init__(self, wind, sigma, ndiv, mode, fwhm=None): 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): @@ -1756,6 +1799,24 @@ def fun(self, param): Used by scipy.optimize.least_squares. """ + if fitting_cpp is not None and hasattr(fitting_cpp, "gaussian_resid"): + sky, height, xcen, ycen, fwhm = self.get_par(param) + return fitting_cpp.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() @@ -1768,6 +1829,25 @@ def jac(self, param): Used by scipy.optimize.least_squares. """ + if fitting_cpp is not None and hasattr(fitting_cpp, "dgaussian_jac"): + sky, height, xcen, ycen, fwhm = self.get_par(param) + return fitting_cpp.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/fitting_cpp.cpp b/hipercam/fitting_cpp.cpp index 09d9a299..b93cd824 100644 --- a/hipercam/fitting_cpp.cpp +++ b/hipercam/fitting_cpp.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -13,6 +14,508 @@ inline double calc_moffat_alpha(double fwhm, double 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. +inline double moffat_value_at(double x_val, double y_val, double height, + double xcen, double ycen, double alpha, + double tbeta, int xbin, int ybin, int ndiv) { + if (ndiv > 0) { + double prof = 0.0; + double inv_ndiv = 1.0 / static_cast(ndiv); + double norm = height / xbin / ybin / (ndiv * ndiv); + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy * inv_ndiv; + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx * inv_ndiv; + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + prof += norm * std::pow(1.0 + alpha * rsq, -tbeta); + } + } + } + } + return 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); +} + +inline double gaussian_value_at(double x_val, double y_val, double height, + double xcen, double ycen, double alpha, + int xbin, int ybin, int ndiv) { + if (ndiv > 0) { + double prof = 0.0; + double inv_ndiv = 1.0 / static_cast(ndiv); + double norm = height / xbin / ybin / (ndiv * ndiv); + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy * inv_ndiv; + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx * inv_ndiv; + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + prof += std::exp(-alpha * rsq); + } + } + } + } + return norm * prof; + } + + double dx = x_val - xcen; + double dy = y_val - ycen; + double rsq = dx * dx + dy * dy; + return height * std::exp(-alpha * rsq); +} + +// Derivatives at one pixel coordinate. The outputs follow the same +// normalization as dmoffat_cpp and dgaussian_cpp so the fit path stays +// numerically equivalent. +inline 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, int xbin, int ybin, int ndiv, + 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 (ndiv > 0) { + double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); + double inv_ndiv = 1.0 / static_cast(ndiv); + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy * inv_ndiv; + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx * inv_ndiv; + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - ycen; + double rsq = dx * dx + dy * dy; + + double denom = 1.0 + alpha * rsq; + double save1 = height * std::pow(denom, -tbeta - 1.0); + double save2 = save1 * rsq; + + double dh = std::pow(denom, -tbeta); + dheight += dh; + 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 * dh + dbeta_coeff * save2); + } + } + } + } + } + + 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; + double save1 = height * std::pow(denom, -tbeta - 1.0); + double save2 = save1 * rsq; + + dheight = std::pow(denom, -tbeta); + 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); + } +} + +inline 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, + int xbin, int ybin, int ndiv, bool comp_dfwhm, + double &dheight, double &dxcen, double &dycen, + double &dfwhm) { + dheight = 0.0; + dxcen = 0.0; + dycen = 0.0; + dfwhm = 0.0; + + if (ndiv > 0) { + double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); + double inv_ndiv = 1.0 / static_cast(ndiv); + double soff = (ndiv - 1.0) / (2.0 * ndiv); + + for (int iy = 0; iy < ybin; ++iy) { + double yoff = iy - (ybin - 1) / 2.0 - soff; + for (int ix = 0; ix < xbin; ++ix) { + double xoff = ix - (xbin - 1) / 2.0 - soff; + for (int isy = 0; isy < ndiv; ++isy) { + double ysoff = yoff + isy * inv_ndiv; + for (int isx = 0; isx < ndiv; ++isx) { + double xsoff = xoff + isx * inv_ndiv; + double dx = x_val + xsoff - xcen; + double dy = y_val + ysoff - 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; + if (comp_dfwhm) { + dfwhm += dfwhm_coeff * dh * rsq; + } + } + } + } + } + + 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; + } +} + +py::array_t +moffat_resid_cpp(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(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || data_info.ndim != 2 || + sigma_info.ndim != 2 || ok_info.ndim != 1 || + x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1] || + x_info.shape[0] != data_info.shape[0] || + x_info.shape[1] != data_info.shape[1] || + x_info.shape[0] != sigma_info.shape[0] || + x_info.shape[1] != sigma_info.shape[1]) { + throw std::runtime_error( + "Input arrays have invalid dimensions or mismatched shapes"); + } + + 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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + size_t n_ok = ok_info.shape[0]; + + double tbeta = std::max(0.01, beta); + double alpha = calc_moffat_alpha(fwhm, beta); + + py::array_t result(static_cast(n_ok)); + py::buffer_info result_info = result.request(); + double *result_ptr = static_cast(result_info.ptr); + + for (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"); + } + + double model = sky + moffat_value_at(x_ptr[idx], y_ptr[idx], height, xcen, + ycen, alpha, tbeta, xbin, ybin, ndiv); + result_ptr[i] = (data_ptr[idx] - model) / sigma_ptr[idx]; + } + + return result; +} + +py::array_t dmoffat_jac_cpp( + 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) { + + (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(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || sigma_info.ndim != 2 || + ok_info.ndim != 1 || x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1] || + x_info.shape[0] != sigma_info.shape[0] || + x_info.shape[1] != sigma_info.shape[1]) { + throw std::runtime_error( + "Input arrays have invalid dimensions or mismatched shapes"); + } + + 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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + size_t n_ok = ok_info.shape[0]; + 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); + + 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); + + for (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"); + } + + double dheight, dxcen, dycen, dfwhm, dbeta; + moffat_derivs_at(x_ptr[idx], y_ptr[idx], height, xcen, ycen, alpha, tbeta, + dfwhm_coeff, dbeta_coeff, xbin, ybin, ndiv, 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 (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_cpp(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_cpp: 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(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || data_info.ndim != 2 || + sigma_info.ndim != 2 || ok_info.ndim != 1 || + x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1] || + x_info.shape[0] != data_info.shape[0] || + x_info.shape[1] != data_info.shape[1] || + x_info.shape[0] != sigma_info.shape[0] || + x_info.shape[1] != sigma_info.shape[1]) { + throw std::runtime_error( + "Input arrays have invalid dimensions or mismatched shapes"); + } + + 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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + size_t n_ok = ok_info.shape[0]; + + double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); + + py::array_t result(static_cast(n_ok)); + py::buffer_info result_info = result.request(); + double *result_ptr = static_cast(result_info.ptr); + + for (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"); + } + + double model = sky + gaussian_value_at(x_ptr[idx], y_ptr[idx], height, xcen, + ycen, alpha, xbin, ybin, ndiv); + result_ptr[i] = (data_ptr[idx] - model) / sigma_ptr[idx]; + } + + return result; +} + +py::array_t dgaussian_jac_cpp( + 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) { + + (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(); + + if (x_info.ndim != 2 || y_info.ndim != 2 || sigma_info.ndim != 2 || + ok_info.ndim != 1 || x_info.shape[0] != y_info.shape[0] || + x_info.shape[1] != y_info.shape[1] || + x_info.shape[0] != sigma_info.shape[0] || + x_info.shape[1] != sigma_info.shape[1]) { + throw std::runtime_error( + "Input arrays have invalid dimensions or mismatched shapes"); + } + + 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); + + size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + size_t n_ok = ok_info.shape[0]; + 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; + + 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); + + for (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"); + } + + double dheight, dxcen, dycen, dfwhm; + gaussian_derivs_at(x_ptr[idx], y_ptr[idx], height, xcen, ycen, alpha, + two_alpha_height, dfwhm_coeff, xbin, ybin, ndiv, + 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 (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 */ @@ -539,6 +1042,21 @@ PYBIND11_MODULE(fitting_cpp, m) { py::arg("ybin"), py::arg("ndiv"), py::arg("comp_dfwhm"), py::arg("comp_dbeta")); + m.def("moffat_resid", &moffat_resid_cpp, + "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_cpp, + "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_cpp, "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"), @@ -549,4 +1067,18 @@ PYBIND11_MODULE(fitting_cpp, m) { 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_cpp, + "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_cpp, + "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")); } From 6e99f02397e41bcf8f28349e67758bafa50d7972 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Fri, 15 May 2026 11:03:36 +0100 Subject: [PATCH 06/11] Add GIL release and OpenMP parallelisation --- hipercam/fitting_cpp.cpp | 59 ++++++++++++++++++++++++++++++++++++++++ setup.py | 12 +++++++- 2 files changed, 70 insertions(+), 1 deletion(-) diff --git a/hipercam/fitting_cpp.cpp b/hipercam/fitting_cpp.cpp index b93cd824..9c4b9d7d 100644 --- a/hipercam/fitting_cpp.cpp +++ b/hipercam/fitting_cpp.cpp @@ -5,6 +5,9 @@ #include #include #include +#ifdef _OPENMP +#include +#endif namespace py = pybind11; @@ -548,6 +551,15 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + // Now the Python operations are complete, release GIL for parallel + // computation from here on. + // pybind11 buffer operations (request, array creation) require the GIL, + // so we can only release it after extracting all pointers and dimensions. + // This allows the OpenMP parallelization in the loops below to run + // without Python thread contention, while still allowing Python to run + // other tasks in parallel if needed. + py::gil_scoped_release release; + if (ndiv > 0) { // With sub-pixellation std::fill_n(result_ptr, n_pixels, 0.0); @@ -557,6 +569,11 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); +// OpenMP parallelization of the outer loop over pixels with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + // Loop over all pixels for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { double x_val = x_ptr[pixel_idx]; @@ -675,6 +692,15 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + // Now the Python operations are complete, release GIL for parallel + // computation from here on. + // pybind11 buffer operations (request, array creation) require the GIL, + // so we can only release it after extracting all pointers and dimensions. + // This allows the OpenMP parallelization in the loops below to run + // without Python thread contention, while still allowing Python to run + // other tasks in parallel if needed. + py::gil_scoped_release release; + // Initialize dsky to ones (derivative of sky is always 1) std::fill_n(dsky_ptr, n_pixels, 1.0); @@ -686,6 +712,11 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); +// OpenMP parallelization of the outer loop over pixels with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + // Loop over all pixels for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { double x_val = x_ptr[pixel_idx]; @@ -826,6 +857,15 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + // Now the Python operations are complete, release GIL for parallel + // computation from here on. + // pybind11 buffer operations (request, array creation) require the GIL, + // so we can only release it after extracting all pointers and dimensions. + // This allows the OpenMP parallelization in the loops below to run + // without Python thread contention, while still allowing Python to run + // other tasks in parallel if needed. + py::gil_scoped_release release; + if (ndiv > 0) { // With sub-pixellation std::fill_n(result_ptr, n_pixels, 0.0); @@ -835,6 +875,11 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); +// OpenMP parallelization of the outer loop over pixels with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + // Loop over all pixels for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { double x_val = x_ptr[pixel_idx]; @@ -938,6 +983,15 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + // Now the Python operations are complete, release GIL for parallel + // computation from here on. + // pybind11 buffer operations (request, array creation) require the GIL, + // so we can only release it after extracting all pointers and dimensions. + // This allows the OpenMP parallelization in the loops below to run + // without Python thread contention, while still allowing Python to run + // other tasks in parallel if needed. + py::gil_scoped_release release; + // Initialize dsky to ones (derivative of sky is always 1) std::fill_n(dsky_ptr, n_pixels, 1.0); @@ -949,6 +1003,11 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, // Mean offset within sub-pixels double soff = (ndiv - 1.0) / (2.0 * ndiv); +// OpenMP parallelization of the outer loop over pixels with SIMD vectorization. +#ifdef _OPENMP +#pragma omp parallel for simd +#endif + // Loop over all pixels for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { double x_val = x_ptr[pixel_idx]; diff --git a/setup.py b/setup.py index 327d3196..108d61b4 100755 --- a/setup.py +++ b/setup.py @@ -4,6 +4,7 @@ """ import os +import sys # need for Cython and pybind11 import numpy as np @@ -26,6 +27,13 @@ ] # 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"] + pybind11_extensions = [ Pybind11Extension( "hipercam.fitting_cpp", @@ -36,7 +44,9 @@ pybind11.get_include(user=True), ], language="c++", - extra_compile_args=["-std=c++11", "-O3", "-ffast-math"], + extra_compile_args=["-std=c++11", "-O3", "-ffast-math", "-march=native"] + + openmp_args, + extra_link_args=openmp_link_args, ), ] From 6a7a374ea0077a9f5f9417ad1af4dbc32df35dbb Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Fri, 15 May 2026 16:55:03 +0100 Subject: [PATCH 07/11] More optimisations --- hipercam/fitting_cpp.cpp | 669 +++++++++++++++++++++------------------ 1 file changed, 356 insertions(+), 313 deletions(-) diff --git a/hipercam/fitting_cpp.cpp b/hipercam/fitting_cpp.cpp index 9c4b9d7d..dba14e39 100644 --- a/hipercam/fitting_cpp.cpp +++ b/hipercam/fitting_cpp.cpp @@ -11,7 +11,39 @@ namespace py = pybind11; -// Helper function to calculate the Moffat alpha parameter +// 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 + +// 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. +inline 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. inline 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); @@ -22,30 +54,23 @@ inline double calc_moffat_alpha(double fwhm, double beta) { // full 2D model arrays when only masked pixels are needed. inline double moffat_value_at(double x_val, double y_val, double height, double xcen, double ycen, double alpha, - double tbeta, int xbin, int ybin, int ndiv) { - if (ndiv > 0) { + 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_ndiv = 1.0 / static_cast(ndiv); - double norm = height / xbin / ybin / (ndiv * ndiv); - double soff = (ndiv - 1.0) / (2.0 * ndiv); - - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - prof += norm * std::pow(1.0 + alpha * rsq, -tbeta); - } - } + 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 prof; + return height * inv_nadd * prof; } double dx = x_val - xcen; @@ -54,32 +79,27 @@ inline double moffat_value_at(double x_val, double y_val, double height, 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. inline double gaussian_value_at(double x_val, double y_val, double height, double xcen, double ycen, double alpha, - int xbin, int ybin, int ndiv) { - if (ndiv > 0) { + const std::vector &x_offsets, + const std::vector &y_offsets) { + if (!x_offsets.empty() && !y_offsets.empty()) { double prof = 0.0; - double inv_ndiv = 1.0 / static_cast(ndiv); - double norm = height / xbin / ybin / (ndiv * ndiv); - double soff = (ndiv - 1.0) / (2.0 * ndiv); - - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - prof += std::exp(-alpha * rsq); - } - } + 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 norm * prof; + return height * inv_nadd * prof; } double dx = x_val - xcen; @@ -88,16 +108,16 @@ inline double gaussian_value_at(double x_val, double y_val, double height, return height * std::exp(-alpha * rsq); } -// Derivatives at one pixel coordinate. The outputs follow the same -// normalization as dmoffat_cpp and dgaussian_cpp so the fit path stays -// numerically equivalent. -inline 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, int xbin, int ybin, int ndiv, - bool comp_dfwhm, bool comp_dbeta, double &dheight, - double &dxcen, double &dycen, double &dfwhm, - double &dbeta) { +// Calculate Moffat derivatives at one pixel coordinate. +// The outputs follow the same normalization as dmoffat_cpp so the fit path +// stays numerically equivalent. +inline 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; @@ -105,41 +125,77 @@ inline void moffat_derivs_at(double x_val, double y_val, double height, dfwhm = 0.0; dbeta = 0.0; - if (ndiv > 0) { - double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); - double inv_ndiv = 1.0 / static_cast(ndiv); - double soff = (ndiv - 1.0) / (2.0 * ndiv); - - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - - double denom = 1.0 + alpha * rsq; - double save1 = height * std::pow(denom, -tbeta - 1.0); - double save2 = save1 * rsq; - - double dh = std::pow(denom, -tbeta); - dheight += dh; - 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 * dh + dbeta_coeff * save2); - } - } + 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; } } } @@ -160,10 +216,9 @@ inline void moffat_derivs_at(double x_val, double y_val, double height, double dy = y_val - ycen; double rsq = dx * dx + dy * dy; double denom = 1.0 + alpha * rsq; - double save1 = height * std::pow(denom, -tbeta - 1.0); - double save2 = save1 * 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) { @@ -175,42 +230,50 @@ inline void moffat_derivs_at(double x_val, double y_val, double height, } } +// Calculate Gaussian derivatives at one pixel coordinate. +// The outputs follow the same normalization as dgaussian_cpp so the fit path +// stays numerically equivalent. inline 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, - int xbin, int ybin, int ndiv, bool comp_dfwhm, - double &dheight, double &dxcen, double &dycen, - double &dfwhm) { + 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 (ndiv > 0) { - double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); - double inv_ndiv = 1.0 / static_cast(ndiv); - double soff = (ndiv - 1.0) / (2.0 * ndiv); - - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - 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; - if (comp_dfwhm) { - dfwhm += dfwhm_coeff * dh * rsq; - } - } + 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; } } } @@ -275,6 +338,8 @@ moffat_resid_cpp(py::array_t x, py::array_t y, 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(); @@ -285,9 +350,19 @@ moffat_resid_cpp(py::array_t x, py::array_t y, if (idx < 0 || static_cast(idx) >= n_pixels) { throw std::runtime_error("ok_indices contains out-of-range values"); } + } - double model = sky + moffat_value_at(x_ptr[idx], y_ptr[idx], height, xcen, - ycen, alpha, tbeta, xbin, ybin, ndiv); + // 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 (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]; } @@ -300,6 +375,7 @@ py::array_t dmoffat_jac_cpp( 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 @@ -334,6 +410,8 @@ py::array_t dmoffat_jac_cpp( 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)}); @@ -345,10 +423,17 @@ py::array_t dmoffat_jac_cpp( if (idx < 0 || static_cast(idx) >= n_pixels) { throw std::runtime_error("ok_indices contains out-of-range values"); } + } + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + for (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, xbin, ybin, ndiv, comp_dfwhm, + dfwhm_coeff, dbeta_coeff, x_offsets, y_offsets, comp_dfwhm, comp_dbeta, dheight, dxcen, dycen, dfwhm, dbeta); double d0 = 1.0; @@ -426,6 +511,8 @@ gaussian_resid_cpp(py::array_t x, py::array_t y, 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(); @@ -436,9 +523,16 @@ gaussian_resid_cpp(py::array_t x, py::array_t y, if (idx < 0 || static_cast(idx) >= n_pixels) { throw std::runtime_error("ok_indices contains out-of-range values"); } + } + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + for (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, xbin, ybin, ndiv); + ycen, alpha, x_offsets, y_offsets); result_ptr[i] = (data_ptr[idx] - model) / sigma_ptr[idx]; } @@ -451,6 +545,7 @@ py::array_t dgaussian_jac_cpp( 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. @@ -481,6 +576,8 @@ py::array_t dgaussian_jac_cpp( 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)}); @@ -492,10 +589,17 @@ py::array_t dgaussian_jac_cpp( if (idx < 0 || static_cast(idx) >= n_pixels) { throw std::runtime_error("ok_indices contains out-of-range values"); } + } + + // Release the GIL now the Python operations are complete. + py::gil_scoped_release release; + + for (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, xbin, ybin, ndiv, + two_alpha_height, dfwhm_coeff, x_offsets, y_offsets, comp_dfwhm, dheight, dxcen, dycen, dfwhm); double d0 = 1.0; @@ -519,9 +623,7 @@ py::array_t dgaussian_jac_cpp( return result; } -/** - * C++ implementation of the Moffat profile function - */ +// C++ implementation of the Moffat profile function py::array_t moffat_cpp(py::array_t x, py::array_t y, double sky, double height, double xcen, double ycen, double fwhm, double beta, int xbin, @@ -543,7 +645,7 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, // 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 *result_ptr = static_cast(result_info.ptr); + double *RESTRICT result_ptr = static_cast(result_info.ptr); // Calculate Moffat profile parameters double tbeta = std::max(0.01, beta); @@ -551,25 +653,17 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - // Now the Python operations are complete, release GIL for parallel - // computation from here on. - // pybind11 buffer operations (request, array creation) require the GIL, - // so we can only release it after extracting all pointers and dimensions. - // This allows the OpenMP parallelization in the loops below to run - // without Python thread contention, while still allowing Python to run - // other tasks in parallel if needed. + // Release the GIL now the Python operations are complete. py::gil_scoped_release release; if (ndiv > 0) { // With sub-pixellation - std::fill_n(result_ptr, n_pixels, 0.0); - double norm = height / xbin / ybin / (ndiv * ndiv); - double inv_ndiv = 1.0 / static_cast(ndiv); - - // Mean offset within sub-pixels - double soff = (ndiv - 1.0) / (2.0 * ndiv); + 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 the outer loop over pixels with SIMD vectorization. +// OpenMP parallelization of for loop with SIMD vectorization. #ifdef _OPENMP #pragma omp parallel for simd #endif @@ -580,25 +674,17 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, double y_val = y_ptr[pixel_idx]; double prof = 0.0; - // Loop over sub-pixels - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - prof += norm * std::pow(1.0 + alpha * rsq, -tbeta); - } - } + // 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 + prof; + result_ptr[pixel_idx] = sky + height * inv_nadd * prof; } } else { // Fast calculation at pixel centers @@ -606,21 +692,22 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, double dx = x_ptr[i] - xcen; double dy = y_ptr[i] - ycen; double rsq = dx * dx + dy * dy; - result_ptr[i] = height * std::pow(1.0 + alpha * rsq, -tbeta) + sky; + result_ptr[i] = sky + height * std::pow(1.0 + alpha * rsq, -tbeta); } } return result; } -/** - * C++ implementation of the Moffat profile derivatives - */ +// C++ implementation of the Moffat profile derivatives std::vector> dmoffat_cpp(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(); @@ -662,15 +749,15 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, py::buffer_info dxcen_info = dxcen.request(); py::buffer_info dycen_info = dycen.request(); - double *dsky_ptr = static_cast(dsky_info.ptr); - double *dheight_ptr = static_cast(dheight_info.ptr); - double *dxcen_ptr = static_cast(dxcen_info.ptr); - double *dycen_ptr = static_cast(dycen_info.ptr); + 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 *dfwhm_ptr = nullptr; - double *dbeta_ptr = nullptr; + double *RESTRICT dfwhm_ptr = nullptr; + double *RESTRICT dbeta_ptr = nullptr; if (comp_dfwhm) { dfwhm_info = dfwhm.request(); @@ -692,13 +779,7 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - // Now the Python operations are complete, release GIL for parallel - // computation from here on. - // pybind11 buffer operations (request, array creation) require the GIL, - // so we can only release it after extracting all pointers and dimensions. - // This allows the OpenMP parallelization in the loops below to run - // without Python thread contention, while still allowing Python to run - // other tasks in parallel if needed. + // Release the GIL now the Python operations are complete. py::gil_scoped_release release; // Initialize dsky to ones (derivative of sky is always 1) @@ -706,13 +787,12 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, if (ndiv > 0) { // With sub-pixellation - double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); - double inv_ndiv = 1.0 / static_cast(ndiv); - - // Mean offset within sub-pixels - double soff = (ndiv - 1.0) / (2.0 * ndiv); + 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 the outer loop over pixels with SIMD vectorization. +// OpenMP parallelization of for loop with SIMD vectorization. #ifdef _OPENMP #pragma omp parallel for simd #endif @@ -721,58 +801,52 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, for (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 dheight = 0.0; - double dxcen = 0.0; - double dycen = 0.0; - double dfwhm = 0.0; - double dbeta = 0.0; - - // Loop over sub-pixels - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - - double denom = 1.0 + alpha * rsq; - double save1 = height * std::pow(denom, -tbeta - 1.0); - double save2 = save1 * rsq; - - // Derivatives - double dh = std::pow(denom, -tbeta); - dheight += dh; - 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 * dh + dbeta_coeff * save2); - } - } + // 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 * inv_nadd; - dxcen_ptr[pixel_idx] = dxcen * inv_nadd; - dycen_ptr[pixel_idx] = dycen * inv_nadd; + 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 * inv_nadd; + dfwhm_ptr[pixel_idx] = dfwhm_sum * inv_nadd; } if (comp_dbeta) { - dbeta_ptr[pixel_idx] = dbeta * inv_nadd; + dbeta_ptr[pixel_idx] = dbeta_sum * inv_nadd; } } } else { @@ -783,11 +857,11 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, double rsq = dx * dx + dy * dy; double denom = 1.0 + alpha * rsq; - double save1 = height * std::pow(denom, -tbeta - 1.0); + dheight_ptr[i] = std::pow(denom, -tbeta); + double save1 = height * dheight_ptr[i] / denom; double save2 = save1 * rsq; // Derivatives - dheight_ptr[i] = std::pow(denom, -tbeta); dxcen_ptr[i] = two_alpha_tbeta * dx * save1; dycen_ptr[i] = two_alpha_tbeta * dy * save1; @@ -826,9 +900,7 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, return result; } -/** - * C++ implementation of the Gaussian profile function - */ +// C++ implementation of the Gaussian profile function py::array_t gaussian_cpp(py::array_t x, py::array_t y, double sky, double height, double xcen, double ycen, double fwhm, int xbin, int ybin, @@ -850,32 +922,24 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, // 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 *result_ptr = static_cast(result_info.ptr); + double *RESTRICT result_ptr = static_cast(result_info.ptr); // Calculate Gaussian profile parameter double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - // Now the Python operations are complete, release GIL for parallel - // computation from here on. - // pybind11 buffer operations (request, array creation) require the GIL, - // so we can only release it after extracting all pointers and dimensions. - // This allows the OpenMP parallelization in the loops below to run - // without Python thread contention, while still allowing Python to run - // other tasks in parallel if needed. + // Release the GIL now the Python operations are complete. py::gil_scoped_release release; if (ndiv > 0) { // With sub-pixellation - std::fill_n(result_ptr, n_pixels, 0.0); - double norm = height / xbin / ybin / (ndiv * ndiv); - double inv_ndiv = 1.0 / static_cast(ndiv); - - // Mean offset within sub-pixels - double soff = (ndiv - 1.0) / (2.0 * ndiv); + 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 the outer loop over pixels with SIMD vectorization. +// OpenMP parallelization of for loop with SIMD vectorization. #ifdef _OPENMP #pragma omp parallel for simd #endif @@ -886,25 +950,17 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, double y_val = y_ptr[pixel_idx]; double prof = 0.0; - // Loop over sub-pixels - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - prof += std::exp(-alpha * rsq); - } - } + // 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 + norm * prof; + result_ptr[pixel_idx] = sky + height * inv_nadd * prof; } } else { // Fast calculation at pixel centers @@ -919,14 +975,15 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, return result; } -/** - * C++ implementation of the Gaussian profile derivatives - */ +// C++ implementation of the Gaussian profile derivatives std::vector> dgaussian_cpp(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(); @@ -963,13 +1020,13 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, py::buffer_info dxcen_info = dxcen.request(); py::buffer_info dycen_info = dycen.request(); - double *dsky_ptr = static_cast(dsky_info.ptr); - double *dheight_ptr = static_cast(dheight_info.ptr); - double *dxcen_ptr = static_cast(dxcen_info.ptr); - double *dycen_ptr = static_cast(dycen_info.ptr); + 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 *dfwhm_ptr = nullptr; + double *RESTRICT dfwhm_ptr = nullptr; if (comp_dfwhm) { dfwhm_info = dfwhm.request(); @@ -983,13 +1040,7 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - // Now the Python operations are complete, release GIL for parallel - // computation from here on. - // pybind11 buffer operations (request, array creation) require the GIL, - // so we can only release it after extracting all pointers and dimensions. - // This allows the OpenMP parallelization in the loops below to run - // without Python thread contention, while still allowing Python to run - // other tasks in parallel if needed. + // Release the GIL now the Python operations are complete. py::gil_scoped_release release; // Initialize dsky to ones (derivative of sky is always 1) @@ -997,13 +1048,12 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, if (ndiv > 0) { // With sub-pixellation - double inv_nadd = 1.0 / static_cast(xbin * ybin * ndiv * ndiv); - double inv_ndiv = 1.0 / static_cast(ndiv); - - // Mean offset within sub-pixels - double soff = (ndiv - 1.0) / (2.0 * ndiv); + 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 the outer loop over pixels with SIMD vectorization. +// OpenMP parallelization of for loop with SIMD vectorization. #ifdef _OPENMP #pragma omp parallel for simd #endif @@ -1012,44 +1062,37 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, for (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 dheight = 0.0; - double dxcen = 0.0; - double dycen = 0.0; - double dfwhm = 0.0; - - // Loop over sub-pixels - for (int iy = 0; iy < ybin; ++iy) { - double yoff = iy - (ybin - 1) / 2.0 - soff; - for (int ix = 0; ix < xbin; ++ix) { - double xoff = ix - (xbin - 1) / 2.0 - soff; - for (int isy = 0; isy < ndiv; ++isy) { - double ysoff = yoff + isy * inv_ndiv; - for (int isx = 0; isx < ndiv; ++isx) { - double xsoff = xoff + isx * inv_ndiv; - double dx = x_val + xsoff - xcen; - double dy = y_val + ysoff - ycen; - double rsq = dx * dx + dy * dy; - - // Gaussian value - double dh = std::exp(-alpha * rsq); - dheight += dh; - dxcen += two_alpha_height * dh * dx; - dycen += two_alpha_height * dh * dy; - - if (comp_dfwhm) { - dfwhm += dfwhm_coeff * dh * rsq; - } - } + // 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 * inv_nadd; - dxcen_ptr[pixel_idx] = dxcen * inv_nadd; - dycen_ptr[pixel_idx] = dycen * inv_nadd; + 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 * inv_nadd; + dfwhm_ptr[pixel_idx] = dfwhm_sum * inv_nadd; } } } else { From 9a9cd5dbef80506fc3f98f4be50311ed578b3a98 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Tue, 11 Aug 2026 16:31:57 +0100 Subject: [PATCH 08/11] Move fully from Cython to Pybind11 --- README.rst | 28 ++- hipercam/{fitting_cpp.cpp => fitting.cpp} | 40 ++-- hipercam/fitting.py | 4 +- hipercam/support.cpp | 136 +++++++++++++ hipercam/support.py | 85 +++++++++ hipercam/support.pyx | 223 ---------------------- pyproject.toml | 3 +- setup.py | 36 ++-- 8 files changed, 274 insertions(+), 281 deletions(-) rename hipercam/{fitting_cpp.cpp => fitting.cpp} (97%) create mode 100644 hipercam/support.cpp create mode 100644 hipercam/support.py delete mode 100644 hipercam/support.pyx 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.cpp b/hipercam/fitting.cpp similarity index 97% rename from hipercam/fitting_cpp.cpp rename to hipercam/fitting.cpp index dba14e39..029ff353 100644 --- a/hipercam/fitting_cpp.cpp +++ b/hipercam/fitting.cpp @@ -109,7 +109,7 @@ inline double gaussian_value_at(double x_val, double y_val, double height, } // Calculate Moffat derivatives at one pixel coordinate. -// The outputs follow the same normalization as dmoffat_cpp so the fit path +// The outputs follow the same normalization as dmoffat so the fit path // stays numerically equivalent. inline void moffat_derivs_at(double x_val, double y_val, double height, double xcen, @@ -231,7 +231,7 @@ moffat_derivs_at(double x_val, double y_val, double height, double xcen, } // Calculate Gaussian derivatives at one pixel coordinate. -// The outputs follow the same normalization as dgaussian_cpp so the fit path +// The outputs follow the same normalization as dgaussian so the fit path // stays numerically equivalent. inline void gaussian_derivs_at(double x_val, double y_val, double height, double xcen, double ycen, double alpha, @@ -300,7 +300,7 @@ inline void gaussian_derivs_at(double x_val, double y_val, double height, } py::array_t -moffat_resid_cpp(py::array_t x, py::array_t y, +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, @@ -369,7 +369,7 @@ moffat_resid_cpp(py::array_t x, py::array_t y, return result; } -py::array_t dmoffat_jac_cpp( +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, @@ -475,13 +475,13 @@ py::array_t dmoffat_jac_cpp( } py::array_t -gaussian_resid_cpp(py::array_t x, py::array_t y, +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_cpp: selected-pixel residuals only. + // Gaussian equivalent of moffat_resid: selected-pixel residuals only. py::buffer_info x_info = x.request(); py::buffer_info y_info = y.request(); @@ -539,7 +539,7 @@ gaussian_resid_cpp(py::array_t x, py::array_t y, return result; } -py::array_t dgaussian_jac_cpp( +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, @@ -624,7 +624,7 @@ py::array_t dgaussian_jac_cpp( } // C++ implementation of the Moffat profile function -py::array_t moffat_cpp(py::array_t x, py::array_t y, +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) { @@ -701,7 +701,7 @@ py::array_t moffat_cpp(py::array_t x, py::array_t y, // C++ implementation of the Moffat profile derivatives std::vector> -dmoffat_cpp(py::array_t x, py::array_t y, double sky, +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) { @@ -901,7 +901,7 @@ dmoffat_cpp(py::array_t x, py::array_t y, double sky, } // C++ implementation of the Gaussian profile function -py::array_t gaussian_cpp(py::array_t x, py::array_t y, +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) { @@ -977,7 +977,7 @@ py::array_t gaussian_cpp(py::array_t x, py::array_t y, // C++ implementation of the Gaussian profile derivatives std::vector> -dgaussian_cpp(py::array_t x, py::array_t y, double sky, +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) { @@ -1129,29 +1129,29 @@ dgaussian_cpp(py::array_t x, py::array_t y, double sky, return result; } -PYBIND11_MODULE(fitting_cpp, m) { +PYBIND11_MODULE(_fitting_cpp, m) { m.doc() = "C++ implementation of profile fitting functions"; - m.def("moffat", &moffat_cpp, "C++ implementation of Moffat profile", + 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_cpp, + 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_cpp, + 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_cpp, + 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"), @@ -1159,25 +1159,25 @@ PYBIND11_MODULE(fitting_cpp, m) { py::arg("ndiv"), py::arg("comp_dfwhm"), py::arg("comp_dbeta"), py::arg("inds")); - m.def("gaussian", &gaussian_cpp, "C++ implementation of Gaussian profile", + 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_cpp, + 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_cpp, + 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_cpp, + 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"), diff --git a/hipercam/fitting.py b/hipercam/fitting.py index 87a815dd..aac71495 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -4,6 +4,8 @@ Moffat profiles plus constants. """ +import importlib + from numba import jit import numpy as np from scipy.optimize import least_squares @@ -11,7 +13,7 @@ from .window import * try: - from . import fitting_cpp + fitting_cpp = importlib.import_module("._fitting_cpp", __package__) except ImportError: fitting_cpp = None diff --git a/hipercam/support.cpp b/hipercam/support.cpp new file mode 100644 index 00000000..6af12489 --- /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 std({ny, nx}); + py::array_t num({ny, nx}); + + auto avg_mut = avg.mutable_unchecked<2>(); + auto std_mut = std.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, std, 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..9276f7a2 --- /dev/null +++ b/hipercam/support.py @@ -0,0 +1,85 @@ +"""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__) +except ImportError: + support_cpp = None + +__all__ = ["avgstd", "gaussian"] + + +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 is not None: + avg, std, num = support_cpp.avgstd(cube_arr, float(sigma)) + return avg, std, num + + # Fallback to pure-NumPy implementation if the extension is not available + 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/pyproject.toml b/pyproject.toml index a1d6f20e..947ec0a7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=61.0", "setuptools-scm", "wheel", "Cython", "numpy", "pybind11>=2.6.0"] +requires = ["setuptools>=61.0", "setuptools-scm", "wheel", "numpy", "pybind11>=2.6.0"] build-backend = "setuptools.build_meta" [project] @@ -22,7 +22,6 @@ classifiers = [ requires-python = ">=3.9" dependencies = [ "astropy", - "Cython", "fitsio", "keyring", "matplotlib", diff --git a/setup.py b/setup.py index 108d61b4..76baf5bc 100755 --- a/setup.py +++ b/setup.py @@ -1,31 +1,16 @@ """ -Minimal setup.py for Cython and pybind11 extension support. +Minimal setup.py for pybind11 extension support. All other metadata is in pyproject.toml. """ import os import sys -# need for Cython and pybind11 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 - -# cython support routine -cython_extensions = [ - Extension( - "hipercam.support", - [os.path.join("hipercam", "support.pyx")], - libraries=["m"], - include_dirs=[np.get_include()], - extra_compile_args=["-fno-strict-aliasing"], - ), -] - # pybind11 extension for profile fitting use_openmp = os.environ.get("HIPERCAM_USE_OPENMP", "1") != "0" openmp_args = [] @@ -36,8 +21,21 @@ pybind11_extensions = [ Pybind11Extension( - "hipercam.fitting_cpp", - ["hipercam/fitting_cpp.cpp"], + "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(), @@ -51,6 +49,6 @@ ] setup( - ext_modules=cythonize(cython_extensions) + pybind11_extensions, + ext_modules=pybind11_extensions, cmdclass={"build_ext": build_ext}, ) From 154f2ebb346d618f4cbcf2951e5118ecf19a2be1 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Wed, 12 Aug 2026 14:20:26 +0100 Subject: [PATCH 09/11] Add comparison test scripts for C++ --- hipercam/fitting.py | 99 +++++---- hipercam/support.py | 19 +- hipercam/tests/fitting_test.py | 383 +++++++++++++++++++++++++++++++++ hipercam/tests/support_test.py | 65 ++++++ 4 files changed, 509 insertions(+), 57 deletions(-) create mode 100644 hipercam/tests/fitting_test.py create mode 100644 hipercam/tests/support_test.py diff --git a/hipercam/fitting.py b/hipercam/fitting.py index aac71495..e05ca6e6 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -13,9 +13,18 @@ from .window import * try: - fitting_cpp = importlib.import_module("._fitting_cpp", __package__) + _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_cpp = None + FITTING_CCP_AVAILABLE = False __all__ = ("combFit", "fitMoffat", "fitGaussian", "moffat", "gaussian") @@ -519,14 +528,10 @@ def moffat(x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv): on the ordinate grids in xy. """ - if fitting_cpp is not None: - return fitting_cpp.moffat( - x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv - ) - else: - return _moffat_numba( - x, y, sky, height, xcen, ycen, fwhm, beta, xbin, ybin, ndiv - ) + 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) @@ -642,24 +647,8 @@ def dmoffat( numba just-in-time compiler function better. """ - if fitting_cpp is not None: - return fitting_cpp.dmoffat( - x, - y, - sky, - height, - xcen, - ycen, - fwhm, - beta, - xbin, - ybin, - ndiv, - comp_dfwhm, - comp_dbeta, - ) - else: - return _dmoffat_numba( + if FITTING_CCP_AVAILABLE: + return _dmoffat_cpp( x, y, sky, @@ -675,6 +664,22 @@ def dmoffat( 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( @@ -1067,9 +1072,9 @@ def fun(self, param): Used by scipy.optimize.least_squares. """ - if fitting_cpp is not None and hasattr(fitting_cpp, "moffat_resid"): + if FITTING_CCP_AVAILABLE: sky, height, xcen, ycen, fwhm, beta = self.get_par(param) - return fitting_cpp.moffat_resid( + return _moffat_resid( self.x, self.y, self.data, @@ -1098,9 +1103,9 @@ def jac(self, param): Used by scipy.optimize.least_squares. """ - if fitting_cpp is not None and hasattr(fitting_cpp, "dmoffat_jac"): + if FITTING_CCP_AVAILABLE: sky, height, xcen, ycen, fwhm, beta = self.get_par(param) - return fitting_cpp.dmoffat_jac( + return _dmoffat_jac( self.x, self.y, self.sigma, @@ -1423,12 +1428,10 @@ def gaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv): on the ordinate grids in xy. """ - if fitting_cpp is not None: - return fitting_cpp.gaussian( - x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv - ) - else: - return _gaussian_numba(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv) + 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) @@ -1528,14 +1531,10 @@ def dgaussian(x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm) appear in the function call. """ - if fitting_cpp is not None: - return fitting_cpp.dgaussian( - x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm - ) - else: - return _dgaussian_numba( - x, y, sky, height, xcen, ycen, fwhm, xbin, ybin, ndiv, comp_dfwhm - ) + 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) @@ -1801,9 +1800,9 @@ def fun(self, param): Used by scipy.optimize.least_squares. """ - if fitting_cpp is not None and hasattr(fitting_cpp, "gaussian_resid"): + if FITTING_CCP_AVAILABLE: sky, height, xcen, ycen, fwhm = self.get_par(param) - return fitting_cpp.gaussian_resid( + return _gaussian_resid( self.x, self.y, self.data, @@ -1831,9 +1830,9 @@ def jac(self, param): Used by scipy.optimize.least_squares. """ - if fitting_cpp is not None and hasattr(fitting_cpp, "dgaussian_jac"): + if FITTING_CCP_AVAILABLE: sky, height, xcen, ycen, fwhm = self.get_par(param) - return fitting_cpp.dgaussian_jac( + return _dgaussian_jac( self.x, self.y, self.sigma, diff --git a/hipercam/support.py b/hipercam/support.py index 9276f7a2..7f935b35 100644 --- a/hipercam/support.py +++ b/hipercam/support.py @@ -10,11 +10,13 @@ import numpy as np try: - support_cpp = importlib.import_module("._support_cpp", __package__) + _support_cpp = importlib.import_module("._support_cpp", __package__) + _avgstd_cpp = _support_cpp.avgstd + SUPPORT_CPP_AVAILABLE = True except ImportError: - support_cpp = None + SUPPORT_CPP_AVAILABLE = False -__all__ = ["avgstd", "gaussian"] +__all__ = ["avgstd"] def avgstd(cube, sigma): @@ -39,11 +41,14 @@ def avgstd(cube, sigma): if sigma <= 1.0: raise ValueError("sigma must be greater than 1") - if support_cpp is not None: - avg, std, num = support_cpp.avgstd(cube_arr, float(sigma)) - return avg, std, num + if SUPPORT_CPP_AVAILABLE: + return _avgstd_cpp(cube_arr, float(sigma)) - # Fallback to pure-NumPy implementation if the extension is not available + 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) 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() From 245c52be1cc5aaf00c795cb16c7adcab51210487 Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Tue, 25 Aug 2026 14:02:46 +0100 Subject: [PATCH 10/11] Add tol options --- hipercam/fitting.py | 57 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 11 deletions(-) diff --git a/hipercam/fitting.py b/hipercam/fitting.py index e05ca6e6..1f830745 100644 --- a/hipercam/fitting.py +++ b/hipercam/fitting.py @@ -45,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 @@ -111,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:: @@ -145,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": @@ -156,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: @@ -214,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 @@ -304,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 @@ -372,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) @@ -1149,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 @@ -1225,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 @@ -1281,7 +1309,14 @@ 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) From b7bf49beddd4e226f7822938a9b5819b2b9e6c4c Mon Sep 17 00:00:00 2001 From: Martin Dyer Date: Tue, 25 Aug 2026 14:45:01 +0100 Subject: [PATCH 11/11] Final C++ cleaning --- hipercam/fitting.cpp | 313 +++++++++++++++++++++---------------------- hipercam/support.cpp | 6 +- 2 files changed, 158 insertions(+), 161 deletions(-) diff --git a/hipercam/fitting.cpp b/hipercam/fitting.cpp index 029ff353..0de6313c 100644 --- a/hipercam/fitting.cpp +++ b/hipercam/fitting.cpp @@ -1,9 +1,11 @@ #include #include +#include #include #include #include #include +#include #include #ifdef _OPENMP #include @@ -20,16 +22,18 @@ namespace py = pybind11; #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. -inline std::vector make_subpixel_offsets(int bin, int ndiv) { +std::vector make_subpixel_offsets(int bin, int ndiv) { std::vector offsets; if (ndiv <= 0) { return offsets; } - offsets.reserve(static_cast(bin * ndiv)); + offsets.reserve(static_cast(bin * ndiv)); double inv_ndiv = 1.0 / static_cast(ndiv); double soff = (ndiv - 1.0) / (2.0 * ndiv); @@ -44,7 +48,7 @@ inline std::vector make_subpixel_offsets(int bin, int ndiv) { } // Helper function to calculate the Moffat alpha parameter. -inline double calc_moffat_alpha(double fwhm, double beta) { +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); } @@ -52,11 +56,10 @@ inline double calc_moffat_alpha(double fwhm, double beta) { // 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. -inline 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) { +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 = @@ -82,10 +85,10 @@ inline double moffat_value_at(double x_val, double y_val, double height, // 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. -inline 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) { +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 = @@ -111,13 +114,13 @@ inline double gaussian_value_at(double x_val, double y_val, double height, // Calculate Moffat derivatives at one pixel coordinate. // The outputs follow the same normalization as dmoffat so the fit path // stays numerically equivalent. -inline 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) { +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; @@ -233,13 +236,13 @@ moffat_derivs_at(double x_val, double y_val, double height, double xcen, // Calculate Gaussian derivatives at one pixel coordinate. // The outputs follow the same normalization as dgaussian so the fit path // stays numerically equivalent. -inline 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) { +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; @@ -299,12 +302,52 @@ inline void gaussian_derivs_at(double x_val, double y_val, double height, } } -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) { +// 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. @@ -315,17 +358,8 @@ moffat_resid(py::array_t x, py::array_t y, py::buffer_info sigma_info = sigma.request(); py::buffer_info ok_info = ok_indices.request(); - if (x_info.ndim != 2 || y_info.ndim != 2 || data_info.ndim != 2 || - sigma_info.ndim != 2 || ok_info.ndim != 1 || - x_info.shape[0] != y_info.shape[0] || - x_info.shape[1] != y_info.shape[1] || - x_info.shape[0] != data_info.shape[0] || - x_info.shape[1] != data_info.shape[1] || - x_info.shape[0] != sigma_info.shape[0] || - x_info.shape[1] != sigma_info.shape[1]) { - throw std::runtime_error( - "Input arrays have invalid dimensions or mismatched shapes"); - } + 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); @@ -333,8 +367,8 @@ moffat_resid(py::array_t x, py::array_t y, const double *sigma_ptr = static_cast(sigma_info.ptr); const std::int64_t *ok_ptr = static_cast(ok_info.ptr); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - size_t n_ok = ok_info.shape[0]; + 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); @@ -345,19 +379,14 @@ moffat_resid(py::array_t x, py::array_t y, py::buffer_info result_info = result.request(); double *result_ptr = static_cast(result_info.ptr); - for (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"); - } - } + 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 (size_t i = 0; i < n_ok; ++i) { + for (std::size_t i = 0; i < n_ok; ++i) { std::int64_t idx = ok_ptr[i]; double model = @@ -369,11 +398,13 @@ moffat_resid(py::array_t x, py::array_t y, 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) { +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; @@ -386,23 +417,17 @@ py::array_t dmoffat_jac( py::buffer_info sigma_info = sigma.request(); py::buffer_info ok_info = ok_indices.request(); - if (x_info.ndim != 2 || y_info.ndim != 2 || sigma_info.ndim != 2 || - ok_info.ndim != 1 || x_info.shape[0] != y_info.shape[0] || - x_info.shape[1] != y_info.shape[1] || - x_info.shape[0] != sigma_info.shape[0] || - x_info.shape[1] != sigma_info.shape[1]) { - throw std::runtime_error( - "Input arrays have invalid dimensions or mismatched shapes"); - } + 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); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - size_t n_ok = ok_info.shape[0]; - size_t n_par = inds.size(); + 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); @@ -418,17 +443,12 @@ py::array_t dmoffat_jac( py::buffer_info result_info = result.request(); double *result_ptr = static_cast(result_info.ptr); - for (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"); - } - } + 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 (size_t i = 0; i < n_ok; ++i) { + for (std::size_t i = 0; i < n_ok; ++i) { std::int64_t idx = ok_ptr[i]; double dheight, dxcen, dycen, dfwhm, dbeta; @@ -462,7 +482,7 @@ py::array_t dmoffat_jac( double derivs[6] = {d0, d1, d2, d3, d4, d5}; double inv_sigma = -1.0 / sigma_ptr[idx]; - for (size_t j = 0; j < n_par; ++j) { + 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"); @@ -474,12 +494,13 @@ py::array_t dmoffat_jac( 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) { +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. @@ -489,17 +510,8 @@ gaussian_resid(py::array_t x, py::array_t y, py::buffer_info sigma_info = sigma.request(); py::buffer_info ok_info = ok_indices.request(); - if (x_info.ndim != 2 || y_info.ndim != 2 || data_info.ndim != 2 || - sigma_info.ndim != 2 || ok_info.ndim != 1 || - x_info.shape[0] != y_info.shape[0] || - x_info.shape[1] != y_info.shape[1] || - x_info.shape[0] != data_info.shape[0] || - x_info.shape[1] != data_info.shape[1] || - x_info.shape[0] != sigma_info.shape[0] || - x_info.shape[1] != sigma_info.shape[1]) { - throw std::runtime_error( - "Input arrays have invalid dimensions or mismatched shapes"); - } + 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); @@ -507,8 +519,8 @@ gaussian_resid(py::array_t x, py::array_t y, const double *sigma_ptr = static_cast(sigma_info.ptr); const std::int64_t *ok_ptr = static_cast(ok_info.ptr); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - size_t n_ok = ok_info.shape[0]; + 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); @@ -518,17 +530,12 @@ gaussian_resid(py::array_t x, py::array_t y, py::buffer_info result_info = result.request(); double *result_ptr = static_cast(result_info.ptr); - for (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"); - } - } + 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 (size_t i = 0; i < n_ok; ++i) { + 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, @@ -539,11 +546,13 @@ gaussian_resid(py::array_t x, py::array_t y, 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) { +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; @@ -555,23 +564,17 @@ py::array_t dgaussian_jac( py::buffer_info sigma_info = sigma.request(); py::buffer_info ok_info = ok_indices.request(); - if (x_info.ndim != 2 || y_info.ndim != 2 || sigma_info.ndim != 2 || - ok_info.ndim != 1 || x_info.shape[0] != y_info.shape[0] || - x_info.shape[1] != y_info.shape[1] || - x_info.shape[0] != sigma_info.shape[0] || - x_info.shape[1] != sigma_info.shape[1]) { - throw std::runtime_error( - "Input arrays have invalid dimensions or mismatched shapes"); - } + 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); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; - size_t n_ok = ok_info.shape[0]; - size_t n_par = inds.size(); + 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; @@ -584,17 +587,12 @@ py::array_t dgaussian_jac( py::buffer_info result_info = result.request(); double *result_ptr = static_cast(result_info.ptr); - for (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"); - } - } + 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 (size_t i = 0; i < n_ok; ++i) { + for (std::size_t i = 0; i < n_ok; ++i) { std::int64_t idx = ok_ptr[i]; double dheight, dxcen, dycen, dfwhm; @@ -611,7 +609,7 @@ py::array_t dgaussian_jac( double derivs[5] = {d0, d1, d2, d3, d4}; double inv_sigma = -1.0 / sigma_ptr[idx]; - for (size_t j = 0; j < n_par; ++j) { + 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"); @@ -625,9 +623,9 @@ py::array_t dgaussian_jac( // 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) { + 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(); @@ -651,7 +649,7 @@ py::array_t moffat(py::array_t x, py::array_t y, double tbeta = std::max(0.01, beta); double alpha = calc_moffat_alpha(fwhm, beta); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + 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; @@ -669,7 +667,7 @@ py::array_t moffat(py::array_t x, py::array_t y, #endif // Loop over all pixels - for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + 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; @@ -688,7 +686,7 @@ py::array_t moffat(py::array_t x, py::array_t y, } } else { // Fast calculation at pixel centers - for (size_t i = 0; i < n_pixels; ++i) { + 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; @@ -701,9 +699,9 @@ py::array_t moffat(py::array_t x, py::array_t y, // 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) { +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; @@ -777,7 +775,7 @@ dmoffat(py::array_t x, py::array_t y, double sky, double dbeta_coeff = 4.0 * std::log(2.0) * std::pow(2.0, 1.0 / tbeta) / tbeta / (fwhm * fwhm); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + 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; @@ -798,7 +796,7 @@ dmoffat(py::array_t x, py::array_t y, double sky, #endif // Loop over all pixels - for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + 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. @@ -851,7 +849,7 @@ dmoffat(py::array_t x, py::array_t y, double sky, } } else { // Fast calculation at pixel centers - for (size_t i = 0; i < n_pixels; ++i) { + 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; @@ -902,9 +900,9 @@ dmoffat(py::array_t x, py::array_t y, double sky, // 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) { + 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(); @@ -927,7 +925,7 @@ py::array_t gaussian(py::array_t x, py::array_t y, // Calculate Gaussian profile parameter double alpha = 4.0 * std::log(2.0) / (fwhm * fwhm); - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + 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; @@ -945,7 +943,7 @@ py::array_t gaussian(py::array_t x, py::array_t y, #endif // Loop over all pixels - for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + 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; @@ -964,7 +962,7 @@ py::array_t gaussian(py::array_t x, py::array_t y, } } else { // Fast calculation at pixel centers - for (size_t i = 0; i < n_pixels; ++i) { + 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; @@ -978,8 +976,8 @@ py::array_t gaussian(py::array_t x, py::array_t y, // 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) { + 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; @@ -1038,7 +1036,7 @@ dgaussian(py::array_t x, py::array_t y, double sky, double two_alpha_height = 2.0 * alpha * height; double dfwhm_coeff = two_alpha_height / fwhm; - size_t n_pixels = x_info.shape[0] * x_info.shape[1]; + 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; @@ -1059,7 +1057,7 @@ dgaussian(py::array_t x, py::array_t y, double sky, #endif // Loop over all pixels - for (size_t pixel_idx = 0; pixel_idx < n_pixels; ++pixel_idx) { + 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. @@ -1097,7 +1095,7 @@ dgaussian(py::array_t x, py::array_t y, double sky, } } else { // Fast calculation at pixel centers - for (size_t i = 0; i < n_pixels; ++i) { + 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; @@ -1132,17 +1130,16 @@ dgaussian(py::array_t x, py::array_t y, double sky, 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"), + 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"), py::arg("comp_dfwhm"), - py::arg("comp_dbeta")); + 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", diff --git a/hipercam/support.cpp b/hipercam/support.cpp index 6af12489..40e35f7c 100644 --- a/hipercam/support.cpp +++ b/hipercam/support.cpp @@ -33,11 +33,11 @@ py::tuple avgstd_impl(const py::array_t &cube, float sigma) { const float *data = static_cast(buf.ptr); py::array_t avg({ny, nx}); - py::array_t std({ny, nx}); + py::array_t stddev({ny, nx}); py::array_t num({ny, nx}); auto avg_mut = avg.mutable_unchecked<2>(); - auto std_mut = std.mutable_unchecked<2>(); + auto std_mut = stddev.mutable_unchecked<2>(); auto num_mut = num.mutable_unchecked<2>(); std::vector vals(nf); @@ -126,7 +126,7 @@ py::tuple avgstd_impl(const py::array_t &cube, float sigma) { } } - return py::make_tuple(avg, std, num); + return py::make_tuple(avg, stddev, num); } } // namespace