From ac2d93cba576c1fbfa5acf8453bb89463a2db67e Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 10 Nov 2025 15:53:14 +0000 Subject: [PATCH 1/2] simplify base flex structure --- README.md | 4 +- setup.py | 13 -- src/flex/__init__.py | 4 +- src/flex/flexcython.py | 195 -------------------------- src/flex/flexnumba.py | 256 ----------------------------------- src/flex/laguerre_cython.pyx | 55 -------- 6 files changed, 3 insertions(+), 524 deletions(-) delete mode 100644 src/flex/flexcython.py delete mode 100644 src/flex/flexnumba.py delete mode 100644 src/flex/laguerre_cython.pyx diff --git a/README.md b/README.md index 9593fd2..170b2aa 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ **Fourier-Laguerre expansions for images of galaxies.** -[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://github.com/michael-petersen/flex/blob/main/LICENSE) +[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://github.com/ObservationalExpansions/flex/blob/main/LICENSE) @@ -11,7 +11,7 @@ Installation of `flex` currently proceeds from local builds after cloning this repository: ``` -git clone https://github.com/michael-petersen/flex.git +git clone https://github.com/ObservationalExpansions/flex.git ``` ``` diff --git a/setup.py b/setup.py index b5bffb6..d4bf07b 100644 --- a/setup.py +++ b/setup.py @@ -1,22 +1,9 @@ from setuptools import setup, Extension -from Cython.Build import cythonize import numpy as np -extensions = [ - Extension( - name="flex.laguerre_cython", - sources=["src/flex/laguerre_cython.pyx"], - include_dirs=[np.get_include()], - ) -] setup( name="flex", packages=["flex"], package_dir={"": "src"}, - ext_modules=cythonize( - extensions, - language_level="3", - annotate=True, # generates .html for optimisation insight - ), ) diff --git a/src/flex/__init__.py b/src/flex/__init__.py index 3d6eec5..a554454 100644 --- a/src/flex/__init__.py +++ b/src/flex/__init__.py @@ -7,9 +7,7 @@ for further details and usage instructions. """ from .flexbase import FLEX -from .flexnumba import FLEXY -from .flexcython import FLEXC from importlib.metadata import version __version__ = version("flex") -__all__ = ["FLEX", "FLEXY", "FLEXC"] +__all__ = ["FLEX"] diff --git a/src/flex/flexcython.py b/src/flex/flexcython.py deleted file mode 100644 index c139a2c..0000000 --- a/src/flex/flexcython.py +++ /dev/null @@ -1,195 +0,0 @@ - - -import numpy as np - -# for the Laguerre polynomials -from .laguerre_cython import laguerre_eval - -class FLEXC: - """ - FLEX class for calculating Laguerre basis amplitudes. - - This class provides methods for calculating Laguerre basis amplitudes based on Weinberg & Petersen (2021). - - Parameters: - rscl (float): Scale parameter for the Laguerre basis. - mass (array-like): Mass values for particles. - phi (array-like): Angular phi values. - velocity (array-like): Velocity values. - R (array-like): Radial values. - mmax (int): Maximum order parameter for m. - nmax (int): Maximum order parameter for n. - - Methods: - gamma_n(nrange, rscl): Calculate the Laguerre alpha=1 normalisation. - G_n(R, nrange, rscl): Calculate the Laguerre basis. - n_m(): Calculate the angular normalisation. - laguerre_amplitudes(): Calculate Laguerre amplitudes for the given parameters. - laguerre_reconstruction(rr, pp): Calculate Laguerre reconstruction. - - Attributes: - rscl (float): Scale parameter for the Laguerre basis. - mass (array-like): Mass values for particles. - phi (array-like): Angular phi values. - velocity (array-like): Velocity values. - R (array-like): Radial values. - mmax (int): Maximum order parameter for m. - nmax (int): Maximum order parameter for n. - coscoefs (array-like): Cosine coefficients. - sincoefs (array-like): Sine coefficients. - reconstruction (array-like): Laguerre reconstruction result. - """ - - def __init__(self, rscl, mmax, nmax, R, phi, mass=1., velocity=1.,newaxis=False): - """ - Initialize the LaguerreAmplitudes instance with parameters. - - Args: - rscl (float): Scale parameter for the Laguerre basis. - mmax (int): Maximum Fourier harmonic order. - nmax (int): Maximum Laguerre order. - R (array-like): Radial values. - velocity (array-like): Velocity values. - mass (integer or array-like): Mass values for particles. - phi (integer or array-like): Angular phi values. - - """ - - # check for input validity - if not isinstance(rscl, (int, float)): - raise ValueError("rscl must be a scalar value.") - if not isinstance(mmax, int) or mmax < 0: - raise ValueError("mmax must be a non-negative integer.") - if not isinstance(nmax, int) or nmax < 0: - raise ValueError("nmax must be a non-negative integer.") - if not isinstance(R, (np.ndarray, list)): - raise ValueError("R must be an array-like structure.") - if not isinstance(phi, (np.ndarray, list)): - raise ValueError("phi must be an array-like structure.") - if not isinstance(mass, (int, float, np.ndarray, list)): - raise ValueError("mass must be a scalar or array-like structure.") - if not isinstance(velocity, (int, float, np.ndarray, list)): - raise ValueError("velocity must be a scalar or array-like structure.") - - self.rscl = rscl - self.mmax = mmax - self.nmax = nmax - self.R = R - self.phi = phi - self.mass = mass - self.velocity = velocity - - # run the amplitude calculation - if newaxis: - self.laguerre_amplitudes_newaxis() - else: - # default behaviour - self.laguerre_amplitudes() - - def _gamma_n(self,nrange, rscl): - """ - Calculate the Laguerre alpha=1 normalisation. - - Args: - nrange (array-like): Range of order parameters. - rscl (float): Scale parameter for the Laguerre basis. - - Returns: - array-like: Laguerre alpha=1 normalisation values. - """ - return (rscl / 2.) * np.sqrt(nrange + 1.) - - def _G_n(self,R, nrange, rscl): - """ - Calculate the Laguerre basis. - - Args: - R (array-like): Radial values. - nrange (array-like): Range of order parameters. - rscl (float): Scale parameter for the Laguerre basis. - - Returns: - array-like: Laguerre basis values. - """ - laguerrevalues = np.array([laguerre_eval(n, 2 * R / rscl)/self._gamma_n(n, rscl) for n in nrange]) - return np.exp(-R / rscl) * laguerrevalues - - def _n_m(self): - """ - Calculate the angular normalisation. - - Returns: - array-like: Angular normalisation values. - """ - deltam0 = np.zeros(self.mmax+1) - - deltam0[0] = 1.0 - - return np.power((deltam0 + 1) * np.pi / 2.,-1.) - - def laguerre_amplitudes_newaxis(self): - """ - Calculate Laguerre amplitudes for the given parameters. - - Returns: - tuple: Tuple containing the cosine and sine amplitudes. - """ - - G_j = self._G_n(self.R, np.arange(0,self.nmax,1), self.rscl) - - nmvals = self._n_m() - cosm = np.array([nmvals[m]*np.cos(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - sinm = np.array([nmvals[m]*np.sin(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - - # broadcast to sum values - self.coscoefs = np.nansum(cosm[:, np.newaxis, :] * G_j[np.newaxis, :, :] * self.mass * self.velocity,axis=2) - self.sincoefs = np.nansum(sinm[:, np.newaxis, :] * G_j[np.newaxis, :, :] * self.mass * self.velocity,axis=2) - - - def laguerre_amplitudes(self): - """ - Calculate Laguerre amplitudes for the given parameters. - - Returns: - none. Attributes are added containing the cosine and sine amplitudes. - """ - - G_j = self._G_n(self.R, np.arange(0,self.nmax,1), self.rscl) - - nmvals = self._n_m() - cosm = np.array([nmvals[m]*np.cos(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - sinm = np.array([nmvals[m]*np.sin(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - - #if scalar: - if np.isscalar(self.mass) and np.isscalar(self.velocity): - scale = self.mass * self.velocity # scalar - self.coscoefs = scale * np.einsum('mn,jn->mj', cosm, G_j) - self.sincoefs = scale * np.einsum('mn,jn->mj', sinm, G_j) - else: - # vector case - self.coscoefs = np.einsum('mn,jn,n->mj', cosm, G_j, self.mass * self.velocity) - self.sincoefs = np.einsum('mn,jn,n->mj', sinm, G_j, self.mass * self.velocity) - - def laguerre_reconstruction(self, rr, pp): - """ - Reconstruct a function using Laguerre amplitudes. - - Args: - rr (array-like): Radial values. - pp (array-like): Angular phi values. - - This method reconstructs a function using the Laguerre amplitudes calculated with the `laguerre_amplitudes` method. - - Returns: - array-like: The reconstructed function values. - """ - nmvals = self._n_m() - G_j = self._G_n(rr, np.arange(0, self.nmax, 1), self.rscl) - - fftotal = 0. - for m in range(0, self.mmax+1): - for n in range(0, self.nmax): - fftotal += self.coscoefs[m, n] * np.cos(m * pp) * G_j[n] - fftotal += self.sincoefs[m, n] * np.sin(m * pp) * G_j[n] - - self.reconstruction = fftotal * 0.5 diff --git a/src/flex/flexnumba.py b/src/flex/flexnumba.py deleted file mode 100644 index 9131b9a..0000000 --- a/src/flex/flexnumba.py +++ /dev/null @@ -1,256 +0,0 @@ -import numpy as np - -import numpy as np -from numba import njit - -@njit -def laguerre_eval(n, alpha, x_vals): - result = np.empty(len(x_vals)) - for i in range(len(x_vals)): - x = x_vals[i] - if n == 0: - result[i] = 1.0 - elif n == 1: - result[i] = 1.0 + alpha - x - else: - Lnm2 = 1.0 - Lnm1 = 1.0 + alpha - x - for k in range(2, n + 1): - L = ((2 * k - 1 + alpha - x) * Lnm1 - (k - 1 + alpha) * Lnm2) / k - Lnm2 = Lnm1 - Lnm1 = L - result[i] = Lnm1 - return result - - -@njit -def laguerre_eval_one(n, x_vals): - """ - Compute the Laguerre polynomial L_n^1(x) for a single n and multiple - x values. - - Parameters: - n (int): Order of the Laguerre polynomial. - x_vals (1D array): Input x values. - Returns: - 1D array of Laguerre polynomial values for each x in x_vals. - """ - result = np.empty(len(x_vals)) - for i in range(len(x_vals)): - x = x_vals[i] - if n == 0: - result[i] = 1.0 - elif n == 1: - result[i] = 2.0 - x - else: - Lnm2 = 1.0 - Lnm1 = 2.0 - x - for k in range(2, n + 1): - L = ((2 * k - x) * Lnm1 - k * Lnm2) / k - Lnm2 = Lnm1 - Lnm1 = L - result[i] = Lnm1 - return result - -@njit -def laguerre_all_orders_one(nmax, x_vals): - """ - Compute L_n^1(x) for all n = 0..nmax and all x in x_vals. - - Parameters: - nmax (int): Maximum order n. - x_vals (1D array): Input x values. - - Returns: - 2D array of shape (nmax+1, len(x_vals)) with L_n^1(x). - """ - nx = len(x_vals) - result = np.empty((nmax + 1, nx)) - - for i in range(nx): - x = x_vals[i] - result[0, i] = 1.0 - if nmax >= 1: - result[1, i] = 2.0 - x - for n in range(2, nmax + 1): - result[n, i] = (2 * n - x) * result[n - 1, i]/ n - result[n - 2, i] - - return result - -class FLEXY: - """ - FLEX class for calculating Laguerre basis amplitudes. - - This class provides methods for calculating Laguerre basis amplitudes based on Weinberg & Petersen (2021). - - Parameters: - rscl (float): Scale parameter for the Laguerre basis. - mass (array-like): Mass values for particles. - phi (array-like): Angular phi values. - velocity (array-like): Velocity values. - R (array-like): Radial values. - mmax (int): Maximum order parameter for m. - nmax (int): Maximum order parameter for n. - - Methods: - gamma_n(nrange, rscl): Calculate the Laguerre alpha=1 normalisation. - G_n(R, nrange, rscl): Calculate the Laguerre basis. - n_m(): Calculate the angular normalisation. - laguerre_amplitudes(): Calculate Laguerre amplitudes for the given parameters. - laguerre_reconstruction(rr, pp): Calculate Laguerre reconstruction. - - Attributes: - rscl (float): Scale parameter for the Laguerre basis. - mass (array-like): Mass values for particles. - phi (array-like): Angular phi values. - velocity (array-like): Velocity values. - R (array-like): Radial values. - mmax (int): Maximum order parameter for m. - nmax (int): Maximum order parameter for n. - coscoefs (array-like): Cosine coefficients. - sincoefs (array-like): Sine coefficients. - reconstruction (array-like): Laguerre reconstruction result. - """ - - def __init__(self, rscl, mmax, nmax, R, phi, mass=1., velocity=1.): - """ - Initialize the LaguerreAmplitudes instance with parameters. - - Args: - rscl (float): Scale parameter for the Laguerre basis. - mmax (int): Maximum Fourier harmonic order. - nmax (int): Maximum Laguerre order. - R (array-like): Radial values. - velocity (array-like): Velocity values. - mass (integer or array-like): Mass values for particles. - phi (integer or array-like): Angular phi values. - - """ - self.rscl = rscl - self.mmax = mmax - self.nmax = nmax - self.R = R - self.phi = phi - self.mass = mass - self.velocity = velocity - - - # run the amplitude calculation - self.laguerre_amplitudes() - - def _gamma_n(self,nrange, rscl): - """ - Calculate the Laguerre alpha=1 normalisation. - - Args: - nrange (array-like): Range of order parameters. - rscl (float): Scale parameter for the Laguerre basis. - - Returns: - array-like: Laguerre alpha=1 normalisation values. - """ - return (rscl / 2.) * np.sqrt(nrange + 1.) - - def _G_n(self, R, nrange, rscl): - """ - Calculate the Laguerre basis. - - Args: - R (array-like): Radial values. - nrange (array-like): Range of order parameters. - rscl (float): Scale parameter for the Laguerre basis. - - Returns: - array-like: Laguerre basis values. - """ - R = np.asarray(R) - x = 2 * R / rscl - gamma_values = np.array([self._gamma_n(n, rscl) for n in nrange]) - - # Preallocate output array for speed - laguerrevalues = np.empty((len(nrange), len(R))) - - for i, n in enumerate(nrange): - laguerrevalues[i] = laguerre_eval_one(n, x) / gamma_values[i] - - return np.exp(-R / rscl) * laguerrevalues - - def _n_m(self): - """ - Calculate the angular normalisation. - - Returns: - array-like: Angular normalisation values. - """ - deltam0 = np.zeros(self.mmax+1) - deltam0[0] = 1.0 - #return np.power((deltam0 + 1) * np.pi / 2., -0.5) - return np.power((deltam0 + 1) * np.pi / 2.,-1.) - #return np.power(deltam0+1,-1.0)#np.power((deltam0 + 1), -0.5) - - def laguerre_amplitudes_newaxis(self): - """ - Calculate Laguerre amplitudes for the given parameters. - - Returns: - tuple: Tuple containing the cosine and sine amplitudes. - """ - - G_j = self._G_n(self.R, np.arange(0,self.nmax,1), self.rscl) - - nmvals = self._n_m() - cosm = np.array([nmvals[m]*np.cos(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - sinm = np.array([nmvals[m]*np.sin(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - - # broadcast to sum values - self.coscoefs = np.nansum(cosm[:, np.newaxis, :] * G_j[np.newaxis, :, :] * self.mass * self.velocity,axis=2) - self.sincoefs = np.nansum(sinm[:, np.newaxis, :] * G_j[np.newaxis, :, :] * self.mass * self.velocity,axis=2) - - - def laguerre_amplitudes(self): - """ - Calculate Laguerre amplitudes for the given parameters. - - Returns: - none. Attributes are added containing the cosine and sine amplitudes. - """ - - G_j = self._G_n(self.R, np.arange(0,self.nmax,1), self.rscl) - - nmvals = self._n_m() - cosm = np.array([nmvals[m]*np.cos(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - sinm = np.array([nmvals[m]*np.sin(m*self.phi) for m in np.arange(0,self.mmax+1,1)]) - - #if scalar: - if np.isscalar(self.mass) and np.isscalar(self.velocity): - scale = self.mass * self.velocity # scalar - self.coscoefs = scale * np.einsum('mn,jn->mj', cosm, G_j) - self.sincoefs = scale * np.einsum('mn,jn->mj', sinm, G_j) - else: - # vector case - self.coscoefs = np.einsum('mn,jn,n->mj', cosm, G_j, self.mass * self.velocity) - self.sincoefs = np.einsum('mn,jn,n->mj', sinm, G_j, self.mass * self.velocity) - - def laguerre_reconstruction(self, rr, pp): - """ - Reconstruct a function using Laguerre amplitudes. - - Args: - rr (array-like): Radial values. - pp (array-like): Angular phi values. - - This method reconstructs a function using the Laguerre amplitudes calculated with the `laguerre_amplitudes` method. - - Returns: - array-like: The reconstructed function values. - """ - nmvals = self._n_m() - G_j = self._G_n(rr, np.arange(0, self.nmax, 1), self.rscl) - - fftotal = 0. - for m in range(0, self.mmax+1): - for n in range(0, self.nmax): - fftotal += self.coscoefs[m, n] * np.cos(m * pp) * G_j[n] - fftotal += self.sincoefs[m, n] * np.sin(m * pp) * G_j[n] - - self.reconstruction = fftotal * 0.5 #/ np.pi diff --git a/src/flex/laguerre_cython.pyx b/src/flex/laguerre_cython.pyx deleted file mode 100644 index d175bb6..0000000 --- a/src/flex/laguerre_cython.pyx +++ /dev/null @@ -1,55 +0,0 @@ -import numpy as np -cimport numpy as np -cimport cython - -@cython.boundscheck(False) -@cython.wraparound(False) -def laguerre_eval(int n, np.ndarray[double, ndim=1] x_vals): - cdef Py_ssize_t i - cdef int k - cdef double x, L, Lnm1, Lnm2 - cdef int size = x_vals.shape[0] - cdef np.ndarray[double, ndim=1] result = np.empty(size) - - for i in range(size): - x = x_vals[i] - if n == 0: - result[i] = 1.0 - elif n == 1: - result[i] = 2.0 - x - else: - Lnm2 = 1.0 - Lnm1 = 2.0 - x - for k in range(2, n + 1): - L = ((2 * k - x) * Lnm1 - k * Lnm2) / k - Lnm2 = Lnm1 - Lnm1 = L - result[i] = Lnm1 - return result - - -@cython.boundscheck(False) -@cython.wraparound(False) -def laguerre_eval_generic(int n, double alpha, np.ndarray[double, ndim=1] x_vals): - cdef Py_ssize_t i - cdef int k - cdef double x, L, Lnm1, Lnm2 - cdef int size = x_vals.shape[0] - cdef np.ndarray[double, ndim=1] result = np.empty(size) - - for i in range(size): - x = x_vals[i] - if n == 0: - result[i] = 1.0 - elif n == 1: - result[i] = 1.0 + alpha - x - else: - Lnm2 = 1.0 - Lnm1 = 1.0 + alpha - x - for k in range(2, n + 1): - L = ((2 * k - 1 + alpha - x) * Lnm1 - (k - 1 + alpha) * Lnm2) / k - Lnm2 = Lnm1 - Lnm1 = L - result[i] = Lnm1 - return result - From 01fb2d4f7d088229b51f427c528690269b4cbad4 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Mon, 10 Nov 2025 16:00:02 +0000 Subject: [PATCH 2/2] Update setup.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index d4bf07b..b8d1303 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,4 @@ -from setuptools import setup, Extension +from setuptools import setup import numpy as np