From e55cc6ef3990785fce3150179663e277ff41a23c Mon Sep 17 00:00:00 2001 From: parsa roshanak Date: Fri, 4 Jul 2025 01:08:15 +0330 Subject: [PATCH 1/2] 0.0.3 c1 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Redesigned `GLpoint`: Now uses a fast, single-pass recurrence to compute the Grünwald-Letnikov fractional derivative at the endpoint. This new implementation is both efficient and direct. - Legacy and alternative endpoint methods (`GLpoint_direct` and `GLpoint_via_GL`) have been moved to `special.py` for reference and comparison. --- changelog.md | 11 ++- differintP/core.py | 156 ++++++++++++++++++++++++++---------------- differintP/special.py | 89 ++++++++++++++++++++++++ pyproject.toml | 2 +- 4 files changed, 196 insertions(+), 62 deletions(-) create mode 100644 differintP/special.py diff --git a/changelog.md b/changelog.md index 4b39e4a..81f603e 100644 --- a/changelog.md +++ b/changelog.md @@ -109,4 +109,13 @@ | 1+e4 | 1.0374 ms | 9.7152 ms | | 1+e5 | 12.5892 ms | 108.2792 ms | | 1+e6 | 126.6893 ms | 995.7722 ms | -| 1+e7 | 1257.0316 ms| * | \ No newline at end of file +| 1+e7 | 1257.0316 ms| * | + + +### 0.0.3 + +- **Redesigned `GLpoint`:** Now uses a fast, single-pass recurrence to compute the Grünwald-Letnikov fractional derivative at the endpoint. This new implementation is both efficient and direct. +- Legacy and alternative endpoint methods (`GLpoint_direct` and `GLpoint_via_GL`) have been moved to `special.py` for reference and comparison. + + + diff --git a/differintP/core.py b/differintP/core.py index a470ba7..e24fb6e 100644 --- a/differintP/core.py +++ b/differintP/core.py @@ -30,7 +30,7 @@ def checkValues( domain_start: int | float, domain_end: int | float, num_points: int, - support_complex_alpha: bool = False + support_complex_alpha: bool = False, ) -> bool | None: """Type checking for valid inputs.""" @@ -158,6 +158,11 @@ def MittagLeffler( return np.sum(exp_vals / gamma_vals, axis=1) +######################################################################################### +########################################## GL ########################################### +######################################################################################### + + def GLcoeffs(alpha: float, n: int) -> np.ndarray: """Vectorized GL coefficient computation""" """ Computes the GL coefficient array of size n. @@ -180,55 +185,6 @@ def GLcoeffs(alpha: float, n: int) -> np.ndarray: return np.cumprod(factors) -def GLpoint( - alpha: float, - f_name: Callable | np.ndarray | list, - domain_start: float = 0.0, - domain_end: float = 1.0, - num_points: int = 100, -) -> float: - """Computes the GL fractional derivative of a function at a point. - - Parameters - ========== - alpha : float - The order of the differintegral to be computed. - f_name : function handle, lambda function, list, or 1d-array of - function values - This is the function that is to be differintegrated. - domain_start : float - The left-endpoint of the function domain. Default value is 0. - domain_end : float - The right-endpoint of the function domain; the point at which the - differintegral is being evaluated. Default value is 1. - num_points : integer - The number of points in the domain. Default value is 100. - - Examples: - >>> DF_poly = GLpoint(-0.5, lambda x: 3*x**2 - 9*x + 2) - >>> DF_sqrt = GLpoint(0.5, lambda x: np.sqrt(x), 0., 1., 100) - """ - # Flip the domain limits if they are in the wrong order. - if domain_start > domain_end: - domain_start, domain_end = domain_end, domain_start - - # Check inputs. - checkValues(alpha, domain_start, domain_end, num_points) - f_values, _ = functionCheck(f_name, domain_start, domain_end, num_points) - - # Calculate the GL differintegral, avoiding the explicit calculation of - # the gamma function. - GL_previous = f_values[1] - for index in range(2, num_points): - GL_current = ( - GL_previous * (num_points - alpha - index - 1) / (num_points - index) - + f_values[index] - ) - GL_previous = GL_current - - return GL_current * (num_points / (domain_end - domain_start)) ** alpha # type: ignore - - def GL( alpha: float, f_name: Callable | np.ndarray | list, @@ -285,6 +241,65 @@ def GL( return result +def GLpoint( + alpha: float, + f_name: Callable[[np.ndarray], np.ndarray] | list[float] | np.ndarray, + domain_start: float = 0.0, + domain_end: float = 1.0, + num_points: int = 100, +) -> float: + """ + Efficiently computes the Grünwald-Letnikov fractional derivative at the endpoint + using vectorized coefficient computation and dot product. + + Parameters + ---------- + alpha : float + The order of the fractional derivative. + f_name : Callable[[np.ndarray], np.ndarray] or Sequence[float] or np.ndarray + The function to differentiate (callable or array-like). + domain_start : float, optional + The starting value of the domain (default is 0.0). + domain_end : float, optional + The ending value of the domain (default is 1.0). + num_points : int, optional + Number of discretization points (default is 100). + + Returns + ------- + float + The Grünwald-Letnikov fractional derivative at the endpoint. + """ + if domain_start > domain_end: + domain_start, domain_end = domain_end, domain_start + + x = np.linspace(domain_start, domain_end, num_points) + if callable(f_name): + f_values = f_name(x) + else: + f_values = np.asarray(f_name) + if len(f_values) != num_points: + raise ValueError("Function array length doesn't match num_points") + + step = (domain_end - domain_start) / (num_points - 1) + step_power = step ** (-alpha) + k = num_points - 1 + + j = np.arange(0, k + 1) + denom = Gamma(j + 1) * Gamma(alpha - j + 1) + with np.errstate(divide="ignore", invalid="ignore"): + coeffs = Gamma(alpha + 1) / denom + coeffs *= (-1) ** j + coeffs[~np.isfinite(coeffs)] = 0 # set nan/inf to zero + + result = np.dot(coeffs, f_values[::-1]) + return step_power * result + + +######################################################################################### +######################################## GL Gpu ######################################### +######################################################################################### + def _gpu_GLcoeffs( alpha: float, n: int, @@ -342,6 +357,23 @@ def GL_gpu( return cupy_manager.cp.asnumpy(result) # Convert back to CPU # type: ignore +######################################################################################### +##################################### GLI - Crone ####################################### +######################################################################################### + +class GLIinterpolat: + """Class for computing interpolation of function values for the + improved GL algorithm. + + Using a class here helps avoid type flexibility for these constants. + """ + + def __init__(self, alpha): + # Determine coefficients for quadratic interpolation. + self.nxt = alpha * (2 + alpha) / 8 + self.crr = (4 - alpha * alpha) / 4 + self.prv = alpha * (alpha - 2) / 8 + def GLI( alpha: float, @@ -472,6 +504,10 @@ def CRONEfilter(siz, alpha): raise InputError(f_name, "f_name must have dimension <= 2") +######################################################################################### +########################################### RL ########################################## +######################################################################################### + def RLmatrix(alpha, N): """Vectorized RL coefficient matrix generation""" # Precompute all required powers @@ -671,18 +707,14 @@ def RL( return result -class GLIinterpolat: - """Class for computing interpolation of function values for the - improved GL algorithm. - Using a class here helps avoid type flexibility for these constants. - """ - def __init__(self, alpha): - # Determine coefficients for quadratic interpolation. - self.nxt = alpha * (2 + alpha) / 8 - self.crr = (4 - alpha * alpha) / 4 - self.prv = alpha * (alpha - 2) / 8 + + +######################################################################################### +######################################### Caputo ######################################## +######################################################################################### + def CaputoL1point( @@ -944,6 +976,10 @@ def CaputoFromRLpoint( return C +######################################################################################### +######################################## PCsolver ####################################### +######################################################################################### + def PCcoeffs(alpha, j, n): if 1 < alpha: if j == 0: diff --git a/differintP/special.py b/differintP/special.py new file mode 100644 index 0000000..549306c --- /dev/null +++ b/differintP/special.py @@ -0,0 +1,89 @@ +from typing import Callable + +import numpy as np + +from differintP import functionCheck, checkValues, GL + + +def GLpoint_via_GL( + alpha: float, + f_name: Callable[[np.ndarray], np.ndarray] | list[float] | np.ndarray, + domain_start: float = 0.0, + domain_end: float = 1.0, + num_points: int = 100, +) -> float: + """ + Efficiently computes the Grünwald-Letnikov fractional derivative at the endpoint + by evaluating the full array via the optimized GL function and returning the last value. + + Parameters + ---------- + alpha : float + The order of the fractional derivative. + f_name : Callable[[np.ndarray], np.ndarray] or Sequence[float] or np.ndarray + The function to differentiate (callable or array-like). + domain_start : float, optional + The starting value of the domain (default is 0.0). + domain_end : float, optional + The ending value of the domain (default is 1.0). + num_points : int, optional + Number of discretization points (default is 100). + + Returns + ------- + float + The Grünwald-Letnikov fractional derivative at the endpoint. + """ + values = GL(alpha, f_name, domain_start, domain_end, num_points) + return float(values[-1]) + + +def GLpoint_direct( + alpha: float, + f_name: Callable | np.ndarray | list, + domain_start: float = 0.0, + domain_end: float = 1.0, + num_points: int = 100, +) -> float: + """Efficiently computes the Grünwald-Letnikov fractional derivative + of a function at the right endpoint of the domain using a direct single-pass + recurrence relation. + + Parameters + ========== + alpha : float + The order of the differintegral to be computed. + f_name : function handle, lambda function, list, or 1d-array of + function values + This is the function that is to be differintegrated. + domain_start : float + The left-endpoint of the function domain. Default value is 0. + domain_end : float + The right-endpoint of the function domain; the point at which the + differintegral is being evaluated. Default value is 1. + num_points : integer + The number of points in the domain. Default value is 100. + + Examples: + >>> DF_poly = GLpoint(-0.5, lambda x: 3*x**2 - 9*x + 2) + >>> DF_sqrt = GLpoint(0.5, lambda x: np.sqrt(x), 0., 1., 100) + """ + # Flip the domain limits if they are in the wrong order. + if domain_start > domain_end: + domain_start, domain_end = domain_end, domain_start + + # Check inputs. + checkValues(alpha, domain_start, domain_end, num_points) + f_values, _ = functionCheck(f_name, domain_start, domain_end, num_points) + + # Calculate the GL differintegral, avoiding the explicit calculation of + # the gamma function. + GL_previous = f_values[1] + for index in range(2, num_points): + GL_current = ( + GL_previous * (num_points - alpha - index - 1) / (num_points - index) + + f_values[index] + ) + GL_previous = GL_current + + return GL_current * (num_points / (domain_end - domain_start)) ** alpha # type: ignore diff --git a/pyproject.toml b/pyproject.toml index 0626552..1d35e97 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "differintP" -version = "0.0.2" +version = "0.0.3" description = "Modern, pure Python fractional calculus library (fork of differint)" readme = "README.md" requires-python = ">=3.7" From a37b106f7706edc671b083275893871d35d5ca2e Mon Sep 17 00:00:00 2001 From: parsa roshanak Date: Fri, 4 Jul 2025 04:08:31 +0330 Subject: [PATCH 2/2] 0.0.3 c2 --- .gitignore | 2 + changelog.md | 40 ++++++++++++-- differintP/core.py | 135 +++++++++++++++++++++++++++------------------ pyproject.toml | 1 + tests/test_core.py | 12 ++-- 5 files changed, 126 insertions(+), 64 deletions(-) diff --git a/.gitignore b/.gitignore index 7bbc71c..34c869e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +bentchmark.ipynb + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/changelog.md b/changelog.md index 81f603e..62a8a74 100644 --- a/changelog.md +++ b/changelog.md @@ -73,13 +73,15 @@ # differintP Changelog +## 0.0.x + ### 0.0.2 (7/3/2025) - Optimized core functions: `GL`, `RL`, and `RLpoint` - Added GPU-accelerated `GL` (`GL_gpu`) via [CuPy](https://cupy.dev/) (optional dependency) - Modernized package setup and build system -#### GL: +[Benchmark](https://github.com/iparsw/differintC/blob/main/BENCHMARK.md) Result for GL: | count | differintP | differint | | ----- | ----------- | ----------- | @@ -91,7 +93,7 @@ | 1+e7 | * | * | -#### RL +[Benchmark](https://github.com/iparsw/differintC/blob/main/BENCHMARK.md) Result for RL | count | differintP | differint | | ----- | ----------- | ----------- | @@ -100,7 +102,7 @@ | 1+e4 | 1372.8949 ms| 3.8757 ms | -#### RLpoint +[Benchmark](https://github.com/iparsw/differintC/blob/main/BENCHMARK.md) Result for RLpoint | count | differintP | differint | | ----- | ----------- | ----------- | @@ -112,10 +114,40 @@ | 1+e7 | 1257.0316 ms| * | -### 0.0.3 +### 0.0.3 (7/4/2025) - **Redesigned `GLpoint`:** Now uses a fast, single-pass recurrence to compute the Grünwald-Letnikov fractional derivative at the endpoint. This new implementation is both efficient and direct. - Legacy and alternative endpoint methods (`GLpoint_direct` and `GLpoint_via_GL`) have been moved to `special.py` for reference and comparison. +[Benchmark](https://github.com/iparsw/differintC/blob/main/BENCHMARK.md) Result for GLpoint: + +| count | differintP | differint | Speedup Factor | +| ----- | ----------- | ------------ | -------------- | +| 1+e2 | 0.0148 ms | 0.0416 ms | x2.81 | +| 1+e3 | 0.0276 ms | 0.4213 ms | x15.26 | +| 1+e4 | 0.0397 ms | 3.8376 ms | x96.66 | +| 1+e5 | 0.3285 ms | 36.303 ms | x110.5 | +| 1+e6 | 5.267 ms | 385.9205 ms | x73.27 | +| 1+e7 | 51.2645 ms | 3674.3541 ms | x71.67 | + + +- **Optimized `GLI` implementation:** + + - Replaced class-based coefficient calculation with direct variable computation for improved clarity and performance. + - Vectorized the 3-point Lagrange interpolation step for all interior points. + - Switched from per-point convolution in a loop to a single global convolution for the entire array, greatly increasing speed. + - Introduced a hybrid approach: uses direct convolution for arrays smaller than 800 points and FFT convolution for larger arrays to ensure optimal performance. + +[Benchmark](https://github.com/iparsw/differintC/blob/main/BENCHMARK.md) Result for GLI: +| count | differintP | differint | Speedup Factor | +| ----- | ----------- | ------------ | -------------- | +| 1+e2 | 0.7582 ms | 0.0536 ms | x14.14 | +| 1+e3 | 7.391 ms | 0.1832 ms | x40.34 | +| 1+e4 | 128.3709 ms | 0.6801 ms | x188.7 | +| 2+e4 | 15181 ms | 5.3272 ms | x2849 | +- Wiki Pages for + - `CaputoL1point` `CaputoL2point` `CaputoL2Cpoint` `CaputoFromRLpoint` + - `CRONE` `GL` `GLpoint` `GLI` `GL_gpu` + - `RL` `RLpoint` diff --git a/differintP/core.py b/differintP/core.py index e24fb6e..3e08c17 100644 --- a/differintP/core.py +++ b/differintP/core.py @@ -2,11 +2,11 @@ from typing import Callable, cast -from git import Optional import numpy as np from scipy.special import gamma as Gamma - +from scipy.signal import fftconvolve +from numba import njit # CuPy dependency for GPU-accelerated GL_gpu from .gpu_utils import cupy_manager @@ -70,7 +70,7 @@ def functionCheck( f_values = list(map(lambda t: f_name(t), x)) step_size = x[1] - x[0] else: - f_name = cast(np.ndarray | list, f_name) + f_name = cast(np.ndarray | list[float], f_name) num_points = np.size(f_name) f_values = f_name step_size = (domain_end - domain_start) / (num_points - 1) @@ -241,6 +241,17 @@ def GL( return result +@njit +def _GLpoint_loop(alpha: float, f_values: np.ndarray, step: float) -> float: + k = f_values.shape[0] - 1 + acc = 0.0 + c_val = 1.0 + for j in range(k + 1): + acc += c_val * f_values[k - j] + if j < k: + c_val *= (-alpha + j) / (j + 1) + return step ** (-alpha) * acc + def GLpoint( alpha: float, f_name: Callable[[np.ndarray], np.ndarray] | list[float] | np.ndarray, @@ -249,8 +260,8 @@ def GLpoint( num_points: int = 100, ) -> float: """ - Efficiently computes the Grünwald-Letnikov fractional derivative at the endpoint - using vectorized coefficient computation and dot product. + Efficient, robust single-point Grünwald-Letnikov fractional derivative + using a direct recurrence (C++-style) in a JIT-compiled kernel. Parameters ---------- @@ -282,18 +293,8 @@ def GLpoint( raise ValueError("Function array length doesn't match num_points") step = (domain_end - domain_start) / (num_points - 1) - step_power = step ** (-alpha) - k = num_points - 1 - - j = np.arange(0, k + 1) - denom = Gamma(j + 1) * Gamma(alpha - j + 1) - with np.errstate(divide="ignore", invalid="ignore"): - coeffs = Gamma(alpha + 1) / denom - coeffs *= (-1) ** j - coeffs[~np.isfinite(coeffs)] = 0 # set nan/inf to zero - result = np.dot(coeffs, f_values[::-1]) - return step_power * result + return _GLpoint_loop(alpha, f_values, step) ######################################################################################### @@ -382,32 +383,56 @@ def GLI( domain_end: float = 1.0, num_points: int = 100, ) -> np.ndarray: - """Computes the 'improved' GL fractional derivative of a function for an - entire array of function values. The 'improved' definition uses the - 3-point Lagrange interpolation found in: + """ + Computes the 'improved' Grünwald-Letnikov (GL) fractional derivative of a function over an entire domain, + using a hybrid approach for efficiency. - Oldham, K. & Spanier, J. (1974). The Fractional Calculus: Theory - and Applications of Differentiation and Integration to Arbitrary - Order. Academic Press, Inc. + This implementation applies a 3-point Lagrange interpolation (from Oldham & Spanier, 1974) + to the function values before convolution, resulting in a higher-order and more accurate + fractional derivative estimate compared to the standard GL method. - Parameters - ========== - alpha : float - The order of the differintegral to be computed. - f_name : function handle, lambda function, list, or 1d-array of - function values - This is the function that is to be differintegrated. - domain_start : float - The left-endpoint of the function domain. Default value is 0. - domain_end : float - The right-endpoint of the function domain; the point at which the - differintegral is being evaluated. Default value is 1. - num_points : integer - The number of points in the domain. Default value is 100. - - Examples: - >>> GLI_poly = GLI(-0.5, lambda x: x**2 - 1) - >>> GLI_sqrt = GLI(0.5, lambda x: np.sqrt(x), 0., 1., 100) + For best performance, a hybrid approach is used: + - For arrays smaller than 800 points, a direct NumPy convolution is performed. + - For larger arrays, a fast FFT-based convolution is used via `scipy.signal.fftconvolve`. + + Parameters + ---------- + alpha : float + Order of the fractional differintegral to compute. + f_name : Callable, list, or 1d-array + Function, lambda, or sequence of function values to be differintegrated. + If callable, it will be evaluated at `num_points` evenly spaced points. + domain_start : float, optional + Start of the domain. Default is 0.0. + domain_end : float, optional + End of the domain. Default is 1.0. + num_points : int, optional + Number of points to evaluate in the domain. Default is 100. + + Returns + ------- + np.ndarray + Array of 'improved' GL fractional derivative values at each grid point. + + Notes + ----- + This method implements the "improved" Grünwald-Letnikov definition by first applying + quadratic interpolation to each interior function value: + interpolated[i] = a * f[i-1] + b * f[i] + c * f[i+1] + where the coefficients (a, b, c) depend on `alpha` and are chosen to reduce discretization error. + + The resulting interpolated array is then convolved with the GL binomial coefficients, + with the result scaled by the step size raised to `-alpha`. For performance, the function + chooses between direct and FFT convolution based on array size. + + References + ---------- + Oldham, K. & Spanier, J. (1974). The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press. + + Examples + -------- + >>> GLI_poly = GLI(-0.5, lambda x: x**2 - 1) + >>> GLI_sqrt = GLI(0.5, lambda x: np.sqrt(x), 0., 1., 100) """ # Flip the domain limits if they are in the wrong order. @@ -419,22 +444,25 @@ def GLI( f_values, step_size = functionCheck(f_name, domain_start, domain_end, num_points) # Get interpolating values. - IN = GLIinterpolat(0.5) - I = [IN.prv, IN.crr, IN.nxt] + prv = alpha * (alpha - 2) / 8 + crr = (4 - alpha * alpha) / 4 + nxt = alpha * (2 + alpha) / 8 - # Get array of generalized binomial coefficients. - b_coeffs = GLcoeffs(0.5, num_points) + # Build interpolated values (ignoring endpoints for simplicity) + interpolated = np.zeros_like(f_values) + interpolated[1:-1] = prv * f_values[:-2] + crr * f_values[1:-1] + nxt * f_values[2:] # type: ignore - # Calculate the improved GL differintegral using convolution. - GLI = np.zeros(num_points) - for i in range(3, num_points): - F = f_values[:i] - L = len(F) - B = b_coeffs[: (L - 2)] - G = np.convolve(F, B, "valid") - GLI[i] = sum(G * I) + # Precompute coefficients + b_coeffs = GLcoeffs(alpha, num_points - 1) - return GLI * step_size**-alpha + # Single convolution for whole array (valid part only) + if num_points < 800: + result = np.convolve(interpolated, b_coeffs, mode="full")[:num_points] + return result * step_size ** -alpha + #FFT convolution + else: + result = fftconvolve(interpolated, b_coeffs, mode="full")[:num_points] + return result * step_size ** -alpha def CRONE(alpha, f_name): @@ -716,7 +744,6 @@ def RL( ######################################################################################### - def CaputoL1point( alpha: float, f_name: Callable | np.ndarray | list, diff --git a/pyproject.toml b/pyproject.toml index 1d35e97..f56fedc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ dependencies = [ # Add your runtime dependencies here, example: "numpy>=1.19", "scipy>=1.15.3", + "numba>=0.61.2" ] [project.optional-dependencies] diff --git a/tests/test_core.py b/tests/test_core.py index 7e8fb7f..fbd8433 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -2,7 +2,7 @@ import numpy as np # Import from sibling directory. -from differintP.core import * +from differintP.core import * # type: ignore # Define constants to be used in tests. poch_first_argument = 1 @@ -85,17 +85,17 @@ def test_GL_binomial_coefficient_array_size(self): def test_checkValues(self): with self.assertRaises(AssertionError): - checkValues(0.1, 0, 1, 1.1) + checkValues(0.1, 0, 1, 1.1) # type: ignore with self.assertRaises(AssertionError): - checkValues(0.1, 1j, 2, 100) + checkValues(0.1, 1j, 2, 100) # type: ignore with self.assertRaises(AssertionError): - checkValues(0.1, 1, 2j, 100) + checkValues(0.1, 1, 2j, 100) # type: ignore with self.assertRaises(AssertionError): checkValues(0.1, 0, 1, -100) with self.assertRaises(AssertionError): - checkValues(1 + 1j, 0, 1, 100) + checkValues(1 + 1j, 0, 1, 100) # type: ignore checkValues(0.5, 0, 1, 100, support_complex_alpha=True) - checkValues(1 + 1j, 0, 1, 100, support_complex_alpha=True) + checkValues(1 + 1j, 0, 1, 100, support_complex_alpha=True) # type: ignore alpha_vals = np.array([0.1, 0.2]) domain_vals = np.array([0.1, 1, 2.0, -1]) num_vals = np.array([1.0, 100.0])