From d4850bcb840514ac830d9267d40aaf1980629eee Mon Sep 17 00:00:00 2001 From: parsa roshanak Date: Mon, 7 Jul 2025 21:47:55 +0330 Subject: [PATCH 1/2] 0.0.4 c1 --- .github/workflows/python-package.yml | 11 +- changelog.md | 14 +- differintP/__init__.py | 7 + differintP/weyl.py | 113 +++++++++++ {tests => test_pytest}/__init__.py | 0 test_pytest/test_core.py | 168 ++++++++++++++++ test_pytest/test_functions.py | 50 +++++ test_pytest/test_solver.py | 37 ++++ test_pytest/test_utils.py | 73 +++++++ test_pytest/test_weyl.py | 97 ++++++++++ tests/test_core.py | 278 --------------------------- tests/test_functions.py | 73 ------- tests/test_solver.py | 58 ------ tests/test_utils.py | 79 -------- 14 files changed, 566 insertions(+), 492 deletions(-) create mode 100644 differintP/weyl.py rename {tests => test_pytest}/__init__.py (100%) create mode 100644 test_pytest/test_core.py create mode 100644 test_pytest/test_functions.py create mode 100644 test_pytest/test_solver.py create mode 100644 test_pytest/test_utils.py create mode 100644 test_pytest/test_weyl.py delete mode 100644 tests/test_core.py delete mode 100644 tests/test_functions.py delete mode 100644 tests/test_solver.py delete mode 100644 tests/test_utils.py diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 2d4a5f7..3b0df40 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -36,6 +36,11 @@ jobs: flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics - - name: Test with unittest - run: | - python -m unittest discover + # - name: Test with unittest + # run: | + # python -m unittest discover + + # === Added for pytest === + - name: Test with pytest # [pytest] + run: | # [pytest] + pytest # [pytest] \ No newline at end of file diff --git a/changelog.md b/changelog.md index 4e816a9..0cc9c69 100644 --- a/changelog.md +++ b/changelog.md @@ -156,4 +156,16 @@ - `RL` `RLpoint` `PCsolver` functions Namespaces - Aditional test and test restructuring -- Restructuring the Namespaces \ No newline at end of file +- Restructuring the Namespaces + + +### 0.0.4 + +* **Added Weyl and Riesz Fractional Derivative Functions:** + + * **`Weyl`**: Implements the periodic, right-sided fractional derivative using an FFT-based method. For periodic functions, this operator produces phase-shifted analytical derivatives, and is especially well-suited for spectral methods. + * **`Riesz`**: Implements the symmetric (two-sided) fractional derivative, also using FFT. For pure sines/cosines, the Riesz operator multiplies each Fourier mode by $-|k|^\alpha$, resulting in a real-valued, sign-flipped output (not a phase shift). Especially relevant in physics and PDEs. +* **Testing Update:** + + * Both `Weyl` and `Riesz` are covered by **unittest** and **pytest**-style tests for compatibility and reliability. + * The test suite is transitioned from `unittest` to `pytest` for modern, flexible, and more expressive testing. diff --git a/differintP/__init__.py b/differintP/__init__.py index 95a5b79..f264b2f 100644 --- a/differintP/__init__.py +++ b/differintP/__init__.py @@ -15,9 +15,12 @@ CaputoFromRLpoint, ) +from .weyl import Weyl, Riesz + from . import functions __all__ = [ + # Core "GLcoeffs", "GL", "GLpoint", @@ -33,5 +36,9 @@ "CaputoL2Cpoint", "CaputoFromRLpoint", + # weyl + "Weyl", + "Riesz", + "functions" ] diff --git a/differintP/weyl.py b/differintP/weyl.py new file mode 100644 index 0000000..e9e7d92 --- /dev/null +++ b/differintP/weyl.py @@ -0,0 +1,113 @@ +from typing import Callable, cast, Union, List + +import numpy as np + +def Weyl( + alpha: float, + f_name: Union[np.ndarray, List[float], Callable], + domain_start: float = 0.0, + domain_end: float = 2 * np.pi, + num_points: int = 100, +) -> np.ndarray: + """ + Weyl fractional derivative (periodic, Fourier-based). + + Numerically computes the Weyl (right-sided, periodic) fractional derivative of order `alpha` + for a function on a uniform grid, using the FFT. This method is fast and accurate for + periodic functions on [domain_start, domain_end]. + + References: + - Samko, Kilbas, Marichev, *Fractional Integrals and Derivatives* (see Weyl derivative, Ch. 7) + + Parameters + ---------- + alpha : float + Order of the derivative. + f_name : callable or array-like + Function or array of values to differentiate. + domain_start, domain_end : float + Interval (should cover one period for periodic functions). + num_points : int + Number of grid points. + + Returns + ------- + df : np.ndarray + Array of Weyl fractional derivative values at grid points. + """ + + if callable(f_name): + # f_name = cast(Callable, f_name) # type checking + # If f_name is callable, call it and save to a list. + x = np.linspace(domain_start, domain_end, num_points, endpoint=False) + f_values = f_name(x) + else: + f_name = cast(np.ndarray | list[float], f_name) + num_points = np.size(f_name) + f_values = f_name + + # Compute FFT + fhat = np.fft.fft(f_values) # type: ignore + L = domain_end - domain_start + k = np.fft.fftfreq(num_points, d=L / num_points) * 2 * np.pi # Frequency in radians + + # Fractional derivative in Fourier domain + multiplier = np.zeros_like(k, dtype=complex) + multiplier[1:] = (1j * k[1:]) ** alpha # k=0 stays zero + + fhat_new = fhat * multiplier + df = np.fft.ifft(fhat_new) + return df.real if np.all(np.isreal(f_values)) else df # type: ignore + + +def Riesz( + alpha: float, + f_name: Union[np.ndarray, List[float], Callable], + domain_start: float = 0.0, + domain_end: float = 2 * np.pi, + num_points: int = 100, +) -> np.ndarray: + """ + Riesz fractional derivative (symmetric, Fourier-based). + + Numerically computes the Riesz fractional derivative of order `alpha` + for a function on a uniform grid using the FFT. This operator is + symmetric (unlike Weyl) and real for real input. + + Parameters + ---------- + alpha : float + Order of the derivative. + f_name : callable or array-like + Function or array of values to differentiate. + domain_start, domain_end : float + Interval (should cover one period for periodic functions). + num_points : int + Number of grid points. + + Returns + ------- + df : np.ndarray + Array of Riesz fractional derivative values at grid points. + """ + if callable(f_name): + x = np.linspace(domain_start, domain_end, num_points, endpoint=False) + f_values = f_name(x) + else: + f_values = np.asarray(f_name) + num_points = len(f_values) + + # FFT + fhat = np.fft.fft(f_values) + L = domain_end - domain_start + k = np.fft.fftfreq(num_points, d=L / num_points) * 2 * np.pi # radians + + # Riesz symbol: -|k|^alpha (symmetric) + multiplier = -np.abs(k) ** alpha + + # Apply in Fourier domain + fhat_new = fhat * multiplier + df = np.fft.ifft(fhat_new) + return df.real if np.all(np.isreal(f_values)) else df + + diff --git a/tests/__init__.py b/test_pytest/__init__.py similarity index 100% rename from tests/__init__.py rename to test_pytest/__init__.py diff --git a/test_pytest/test_core.py b/test_pytest/test_core.py new file mode 100644 index 0000000..7f8109c --- /dev/null +++ b/test_pytest/test_core.py @@ -0,0 +1,168 @@ +import numpy as np + +import sys +import os + +# Add the parent directory (containing differintP) to the import path +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) + +from differintP.core import * # type: ignore +from differintP.functions import MittagLeffler +from .__init__ import test_N + +# Define constants to be used in tests. +size_coefficient_array = 20 +sqrtpi2 = 0.88622692545275794 +truevaluepoly = 0.94031597258 +truevaluepoly_caputo = 1.50450555 # 8 / (3 * np.sqrt(np.pi)) +truevaluepoly_caputo_higher = 2 / Gamma(1.5) + +# Get SQRT results for checking accuracy. +GL_r = GL(0.5, lambda x: np.sqrt(x), 0, 1, test_N) +GL_result = GL_r[-1] +GL_length = len(GL_r) + +GLI_r = GLI(0.5, lambda x: np.sqrt(x), 0, 1, test_N) +GLI_result = GLI_r[-1] +GLI_length = len(GLI_r) + +RL_r = RL(0.5, lambda x: np.sqrt(x), 0, 1, test_N) +RL_result = RL_r[-1] +RL_length = len(RL_r) + +# --- Exponential function --- +alpha_exp = 0.5 +groundtruth_exp = MittagLeffler(1, 1 - 0.5, np.array([1]))[0] + +# --- Sine function --- +alpha_sin = 0.5 +groundtruth_sin = np.sin(1 + alpha_sin * np.pi / 2) + +####################### +# GLpoint +####################### + +def test_GLpoint_sqrt_accuracy(): + assert abs(GLpoint(0.5, lambda x: x**0.5, 0.0, 1.0, 1024) - sqrtpi2) <= 1e-3 + +def test_GLpoint_accuracy_polynomial(): + assert abs(GLpoint(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024) - truevaluepoly) <= 1e-3 + +def test_GLpoint_accuracy_exp(): + val = GLpoint(alpha_exp, np.exp, 0, 1, 1024) + assert abs(val - groundtruth_exp) < 1e-3 + +####################### +# GL +####################### + +def test_GL_accuracy_sqrt(): + assert abs(GL_result - sqrtpi2) <= 1e-4 + +def test_GL_accuracy_polynomial(): + assert abs(GL(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) <= 1e-3 + +def test_GL_accuracy_exp(): + val = GL(alpha_exp, np.exp, 0, 1, 1024)[-1] + assert abs(val - groundtruth_exp) < 1e-3 + +####################### +# GLI +####################### + +def test_GLI_accuracy_sqrt(): + assert abs(GLI_result - sqrtpi2) <= 1e-4 + +def test_GLI_accuracy_polynomial(): + assert abs(GLI(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) <= 6e-3 # low accuracy + +def test_GLI_accuracy_exp(): + val = GLI(alpha_exp, np.exp, 0, 1, 1024)[-1] + assert abs(val - groundtruth_exp) < 5e-3 # low accuracy + +####################### +# RLpoint +####################### + +def test_RLpoint_sqrt_accuracy(): + assert abs(RLpoint(0.5, lambda x: x**0.5, 0.0, 1.0, 1024) - sqrtpi2) <= 1e-3 + +def test_RLpoint_accuracy_polynomial(): + assert abs(RLpoint(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024) - truevaluepoly) <= 1e-2 + +def test_RLpoint_accuracy_exp(): + val = RLpoint(alpha_exp, np.exp, 0, 1, 1024) + assert abs(val - groundtruth_exp) < 1e-3 + +####################### +# RL +####################### + +def test_RL_accuracy_sqrt(): + assert abs(RL_result - sqrtpi2) <= 1e-4 + +def test_RL_accuracy_polynomial(): + assert abs(RL(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) <= 1e-3 + +def test_RL_accuracy_exp(): + val = RL(alpha_exp, np.exp, 0, 1, 1024)[-1] + assert abs(val - groundtruth_exp) < 1e-3 + +####################### +# Caputo 1p +####################### + +def test_CaputoL1point_accuracy_sqrt(): + assert abs(CaputoL1point(0.5, lambda x: x**0.5, 0, 1.0, 1024) - sqrtpi2) <= 1e-2 + +def test_CaputoL1point_accuracy_polynomial(): + assert abs(CaputoL1point(0.5, lambda x: x**2 - 1, 0, 1.0, 1024) - truevaluepoly_caputo) <= 1e-3 + +def test_CaputoL1point_accuracy_exp(): + val = CaputoL1point(alpha_exp, np.exp, 0, 1, 1024) + assert abs(val - groundtruth_exp) < 0.6 # really bad accuracy + +####################### +# Caputo 2p +####################### + +def test_CaputoL2point_accuracy_polynomial(): + assert abs(CaputoL2point(1.5, lambda x: x**2, 0, 1.0, 1024) - truevaluepoly_caputo_higher) <= 1e-1 + +####################### +# Caputo 2pC +####################### + +def test_CaputoL2Cpoint_accuracy_polynomial_higher(): + assert abs(CaputoL2Cpoint(1.5, lambda x: x**2, 0, 1.0, 1024) - truevaluepoly_caputo_higher) <= 1e-1 + +def test_CaputoL2Cpoint_accuracy_polynomial(): + assert abs(CaputoL2Cpoint(0.5, lambda x: x**2, 0, 1.0, 1024) - truevaluepoly_caputo) <= 1e-3 + +def test_CaputoL2Cpoint_accuracy_exp(): + val = CaputoL2Cpoint(alpha_exp, np.exp, 0, 1, 1024) + assert abs(val - groundtruth_exp) < 2 # unusable + +####################### +# General algorithm output checks +####################### + +def test_GL_result_length(): + assert GL_length == test_N + +def test_GLI_result_length(): + assert GLI_length == test_N + +def test_RL_result_length(): + assert RL_length == test_N + +def test_RL_matrix_shape(): + assert np.shape(RLmatrix(0.4, test_N)) == (test_N, test_N) + +def test_GL_binomial_coefficient_array_size(): + assert len(GLcoeffs(0.5, size_coefficient_array)) - 1 == size_coefficient_array + +# Optional: Run doctest if called directly +if __name__ == "__main__": + import doctest + doctest.testmod() diff --git a/test_pytest/test_functions.py b/test_pytest/test_functions.py new file mode 100644 index 0000000..71de34a --- /dev/null +++ b/test_pytest/test_functions.py @@ -0,0 +1,50 @@ +import numpy as np + +import sys +import os + +# Add the parent directory (containing differintP) to the import path +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) + +from differintP.functions import * # type: ignore + +# Define constants to be used in tests. +poch_first_argument = 1 +poch_second_argument = 5 +poch_true_answer = 120 + +def test_pochhammer(): + assert poch(poch_first_argument, poch_second_argument) == poch_true_answer + assert poch(-1, 3) == 0 + # assert poch(-1.5, 0.5) == np.inf # commented as in original + assert np.round(poch(1j, 1), 3) == 0.000 + 1.000j + assert poch(-10, 2) == 90 + +def test_ML_cosh_root(): + xs = np.arange(0.1, 10, 0.1) + assert np.all( + np.abs( + MittagLeffler(2, 1, xs, ignore_special_cases=True) - np.cosh(np.sqrt(xs)) + ) <= 1e-3 + ) + +def test_ML_exp(): + xs = np.arange(0.1, 10, 0.1) + assert np.all( + np.abs( + MittagLeffler(1, 1, xs, ignore_special_cases=True) - np.exp(xs) + ) <= 1e-3 + ) + +def test_ML_geometric(): + xs = np.arange(0.05, 1, 0.05) + assert np.all( + np.abs( + MittagLeffler(0, 1, xs, ignore_special_cases=True) - 1 / (1 - xs) + ) <= 1e-3 + ) + +# Optional: Run doctest if called directly +if __name__ == "__main__": + import doctest + doctest.testmod() diff --git a/test_pytest/test_solver.py b/test_pytest/test_solver.py new file mode 100644 index 0000000..0544155 --- /dev/null +++ b/test_pytest/test_solver.py @@ -0,0 +1,37 @@ +import numpy as np + +import sys +import os + +# Add the parent directory (containing differintP) to the import path +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + +from differintP.core import Gamma +from differintP.solvers import PCsolver +from differintP.functions import MittagLeffler + +PC_x_power = np.linspace(0, 1, 100) ** 5.5 + +# Get FODE function for solving. +PC_func_power = lambda x, y: 1 / 24 * Gamma(5 + 1.5) * x**4 + x ** (8 + 2 * 1.5) - y**2 +PC_func_ML = lambda x, y: y + +def test_PC_solution_three_halves(): + result = PCsolver([0, 0], 1.5, PC_func_power, 0, 1, 100) + assert np.all(np.abs(result - PC_x_power) <= 1e-2) + +def test_PC_solution_ML(): + xs = np.linspace(0, 1, 100) + ML_alpha = MittagLeffler(5.5, 1, xs**5.5) + result = PCsolver([1, 0, 0, 0, 0, 0], 5.5, PC_func_ML) + assert np.all(np.abs(result - ML_alpha) <= 1e-2) + +def test_PC_solution_linear(): + xs = np.linspace(0, 1, 100) + result = PCsolver([1, 1], 1.5, lambda x, y: y - x - 1) + assert np.all(np.abs(result - (xs + 1)) <= 1e-2) + +# Optional: Run doctest (only if you want this to run on manual invocation) +if __name__ == "__main__": + import doctest + doctest.testmod() diff --git a/test_pytest/test_utils.py b/test_pytest/test_utils.py new file mode 100644 index 0000000..77b36c0 --- /dev/null +++ b/test_pytest/test_utils.py @@ -0,0 +1,73 @@ +import numpy as np +import pytest + +import sys +import os + +# Add the parent directory (containing differintP) to the import path +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + +from differintP.utils import * # type: ignore +from .__init__ import test_N + +stepsize = 1 / (test_N - 1) + +# Testing if callable functions and arrays of function values will work. +checked_function1, test_stepsize1 = functionCheck( + lambda x: 2 * np.exp(3 * x) * x - x**2 + x - 5, 0, 1, test_N +) +checked_function2, test_stepsize2 = functionCheck(np.ones(test_N), 0, 1, test_N) + +def test_isInteger(): + assert isInteger(1) + assert isInteger(1.0) + assert isInteger(1 + 0j) + assert not isInteger(1.1) + assert not isInteger(1.1 + 0j) + assert not isInteger(1 + 1j) + +def test_isPositiveInteger(): + assert isPositiveInteger(1) + assert not isPositiveInteger(1.1) + assert not isPositiveInteger(-1) + +def test_functionCheck(): + assert len(checked_function1) == test_N + assert len(checked_function2) == test_N + + # Make sure it treats defined functions and arrays of function values the same. + assert len(checked_function1) == len(checked_function2) + assert test_stepsize1 == stepsize + assert test_stepsize2 == stepsize + assert test_stepsize1 == test_stepsize2 + +def test_checkValues(): + with pytest.raises(AssertionError): + checkValues(0.1, 0, 1, 1.1) # type: ignore + with pytest.raises(AssertionError): + checkValues(0.1, 1j, 2, 100) # type: ignore + with pytest.raises(AssertionError): + checkValues(0.1, 1, 2j, 100) # type: ignore + with pytest.raises(AssertionError): + checkValues(0.1, 0, 1, -100) + with pytest.raises(AssertionError): + checkValues(1 + 1j, 0, 1, 100) # type: ignore + # These should not raise + checkValues(0.5, 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]) + [ + [ + [ + [ + checkValues(alpha, domain_start, domain_end, num_points) + for alpha in alpha_vals + ] + for domain_start in domain_vals + ] + for domain_end in domain_vals + ] + for num_points in num_vals + ] diff --git a/test_pytest/test_weyl.py b/test_pytest/test_weyl.py new file mode 100644 index 0000000..cf4f1ea --- /dev/null +++ b/test_pytest/test_weyl.py @@ -0,0 +1,97 @@ +import pytest + +import numpy as np + +import sys +import os + +# Add the parent directory (containing differintP) to the import path +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) + +from differintP.weyl import Weyl, Riesz + +class TestWeylAccuray: + ####################### + # Weyl + ####################### + + def test_Weyl_sin_first_derivative(self): + """Weyl(1, sin(x)) should approximate cos(x) for periodic domain""" + alpha = 1 + domain_start = 0.0 + domain_end = 2 * np.pi + num_points = 128 + x = np.linspace(domain_start, domain_end, num_points, endpoint=False) + gt = np.cos(x) + result = Weyl(alpha, np.sin, domain_start, domain_end, num_points) + max_error = np.max(np.abs(result - gt)) + avg_error = np.mean(np.abs(result - gt)) + assert max_error < 2e-14 # Great Performance + assert avg_error < 5e-15 # Great Performance + + def test_Weyl_sin_half_derivative(self): + """Weyl(0.5, sin(x)) should match analytical fractional derivative""" + alpha = 0.5 + domain_start = 0.0 + domain_end = 2 * np.pi + num_points = 128 + x = np.linspace(domain_start, domain_end, num_points, endpoint=False) + gt = np.sin(x + alpha * np.pi / 2) + result = Weyl(alpha, np.sin, domain_start, domain_end, num_points) + max_error = np.max(np.abs(result - gt)) + avg_error = np.mean(np.abs(result - gt)) + assert max_error < 1e-14 # Great Performance + assert avg_error < 1e-15 # Great Performance + + def test_Weyl_sin_half_derivative_middle_third(self): + """Weyl(0.5, sin(x)) should match analytical fractional derivative (middle third only)""" + alpha = 0.5 + domain_start = 0.0 + domain_end = 2 * np.pi + num_points = 1024 + x = np.linspace(domain_start, domain_end, num_points, endpoint=False) + gt = np.sin(x + alpha * np.pi / 2) + result = Weyl(alpha, np.sin, domain_start, domain_end, num_points) + + # Indices for the middle third + start_idx = num_points // 3 + end_idx = 2 * num_points // 3 + + mid_result = result[start_idx:end_idx] + mid_gt = gt[start_idx:end_idx] + + max_error = np.max(np.abs(mid_result - mid_gt)) + avg_error = np.mean(np.abs(mid_result - mid_gt)) + + assert max_error < 1e-14 # 6.99e-15 + assert avg_error < 1e-14 # 1.97e-15 + + + +####################### +# Riesz +####################### +@pytest.mark.parametrize( + "alpha,max_error_thresh,avg_error_thresh", + [ + (0.1, 1e-14, 1e-15), # Riesz(0.5, sin(x)) -> -sin(x) + (0.3, 1e-14, 1e-15), # Riesz(0.5, sin(x)) -> -sin(x) + (0.5, 1e-14, 1e-15), # Riesz(0.5, sin(x)) -> -sin(x) + (0.7, 2e-14, 5e-15), # Riesz(0.5, sin(x)) -> -sin(x) + (1.0, 2e-14, 5e-15), # Riesz(1.0, sin(x)) -> -sin(x) + (1.5, 1e-13, 5e-14), # Riesz(1.5, sin(x)) -> -sin(x) + ] +) +def test_Riesz_sin_derivative(alpha, max_error_thresh, avg_error_thresh): + """Riesz(alpha, sin(x)) should approximate -sin(x) for periodic domain for various alphas.""" + domain_start = 0.0 + domain_end = 2 * np.pi + num_points = 128 + x = np.linspace(domain_start, domain_end, num_points, endpoint=False) + gt = -np.sin(x) + # import Riesz however you do in your test file, e.g.: + result = Riesz(alpha, np.sin, domain_start, domain_end, num_points) # type: ignore + max_error = np.max(np.abs(result - gt)) + avg_error = np.mean(np.abs(result - gt)) + assert max_error < max_error_thresh + assert avg_error < avg_error_thresh \ No newline at end of file diff --git a/tests/test_core.py b/tests/test_core.py deleted file mode 100644 index 0ecb1b1..0000000 --- a/tests/test_core.py +++ /dev/null @@ -1,278 +0,0 @@ -import unittest -import numpy as np - -# Import from sibling directory. -import sys -import os - -# Add the parent directory (containing differintP) to the import path -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) - -from differintP.core import * -from differintP.functions import MittagLeffler - -from .__init__ import test_N - -# Define constants to be used in tests. -size_coefficient_array = 20 -sqrtpi2 = 0.88622692545275794 -truevaluepoly = 0.94031597258 -truevaluepoly_caputo = 1.50450555 # 8 / (3 * np.sqrt(np.pi)) -truevaluepoly_caputo_higher = 2 / Gamma(1.5) - - -# Get SQRT results for checking accuracy. -GL_r = GL(0.5, lambda x: np.sqrt(x), 0, 1, test_N) -GL_result = GL_r[-1] -GL_length = len(GL_r) - -GLI_r = GLI(0.5, lambda x: np.sqrt(x), 0, 1, test_N) -GLI_result = GLI_r[-1] -GLI_length = len(GLI_r) - -RL_r = RL(0.5, lambda x: np.sqrt(x), 0, 1, test_N) -RL_result = RL_r[-1] -RL_length = len(RL_r) - - -# --- Exponential function --- -# For f(x) = exp(x), the fractional derivative of order alpha at x = 1: -# D^{alpha} e^{x} |_{x=1} = exp(1) / Gamma(1 - alpha) -alpha_exp = 0.5 -groundtruth_exp_expr = "exp(1) * sum_{k=0}^∞ 1^k / Gamma(k + 1 - 0.5)" -groundtruth_exp = MittagLeffler(1, 1 - 0.5, np.array([1]))[0] - - -# --- Sine function --- -# D^{alpha} sin(x) = sin(x + alpha * pi/2) -# For alpha=0.5, x=1: -# Analytical: sin(1 + 0.5 * pi/2) -alpha_sin = 0.5 -groundtruth_sin_expr = "sin(1 + alpha_sin * pi/2)" -groundtruth_sin = np.sin(1 + alpha_sin * np.pi / 2) - - - - - - -class TestAlgorithmsAccuray(unittest.TestCase): - """Tests for algorithm accuracy.""" - - ####################### - # GLpoint - ####################### - - # sqrt - def test_GLpoint_sqrt_accuracy(self): - self.assertTrue( - abs(GLpoint(0.5, lambda x: x**0.5, 0.0, 1.0, 1024) - sqrtpi2) <= 1e-3 - ) - - # x**2 - 1 - def test_GLpoint_accuracy_polynomial(self): - self.assertTrue( - abs(GLpoint(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024) - truevaluepoly) - <= 1e-3 - ) - - # exp - def test_GLpoint_accuracy_exp(self): - """Test GLpoint on f(x) = exp(x), alpha=0.5. Analytical: exp(1)/Gamma(0.5)""" - val = GLpoint(alpha_exp, np.exp, 0, 1, 1024) - self.assertTrue(abs(val - groundtruth_exp) < 1e-3) - - - ####################### - # GL - ####################### - - # sqrt - def test_GL_accuracy_sqrt(self): - self.assertTrue(abs(GL_result - sqrtpi2) <= 1e-4) - - # x**2 - 1 - def test_GL_accuracy_polynomial(self): - self.assertTrue( - abs(GL(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) - <= 1e-3 - ) - - # exp - def test_GL_accuracy_exp(self): - """Test GL on f(x) = exp(x), alpha=0.5. Analytical: exp(1)/Gamma(0.5)""" - val = GL(alpha_exp, np.exp, 0, 1, 1024)[-1] - # print(f"exp: numeric={val}, expected={groundtruth_exp}") - self.assertTrue(abs(val - groundtruth_exp) < 1e-3) - - - ####################### - # GLI - ####################### - - # sqrt - def test_GLI_accuracy_sqrt(self): - self.assertTrue(abs(GLI_result - sqrtpi2) <= 1e-4) - - # x**2 - 1 - def test_GLI_accuracy_polynomial(self): - self.assertTrue( - abs(GLI(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) - <= 6e-3 # low accuracy - ) - - # exp - def test_GLI_accuracy_exp(self): - """Test GLI on f(x) = exp(x), alpha=0.5. Analytical: Gamma(0.5)""" - val = GLI(alpha_exp, np.exp, 0, 1, 1024)[-1] - # print(f"exp: numeric={val}, expected={groundtruth_exp}") - self.assertTrue(abs(val - groundtruth_exp) < 5e-3) # low accuracy - # Abs error: 0.00481 - - - ####################### - # RLpoint - ####################### - - # sqrt - def test_RLpoint_sqrt_accuracy(self): - self.assertTrue( - abs(RLpoint(0.5, lambda x: x**0.5, 0.0, 1.0, 1024) - sqrtpi2) <= 1e-3 - ) - - # poly x**2 - 1 - def test_RLpoint_accuracy_polynomial(self): - self.assertTrue( - abs(RLpoint(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024) - truevaluepoly) - <= 1e-2 - ) - - # exp - def test_RLpoint_accuracy_exp(self): - """Test RLpoint on f(x) = exp(x), alpha=0.5. Analytical: exp(1)/Gamma(0.5)""" - val = RLpoint(alpha_exp, np.exp, 0, 1, 1024) - self.assertTrue(abs(val - groundtruth_exp) < 1e-3) - - - ####################### - # RL - ####################### - - # sqrt - def test_RL_accuracy_sqrt(self): - self.assertTrue(abs(RL_result - sqrtpi2) <= 1e-4) - - # x**2 - 1 - def test_RL_accuracy_polynomial(self): - self.assertTrue( - abs(RL(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) - <= 1e-3 - ) - - # exp - def test_RL_accuracy_exp(self): - """Test RL on f(x) = exp(x), alpha=0.5. Analytical: exp(1)/Gamma(0.5)""" - val = RL(alpha_exp, np.exp, 0, 1, 1024)[-1] - # print(f"exp: numeric={val}, expected={groundtruth_exp}") - self.assertTrue(abs(val - groundtruth_exp) < 1e-3) - - ####################### - # Caputo 1p - ####################### - - # sqrt - def test_CaputoL1point_accuracy_sqrt(self): - self.assertTrue( - abs(CaputoL1point(0.5, lambda x: x**0.5, 0, 1.0, 1024) - sqrtpi2) <= 1e-2 - ) - - # x**2 - 1 - def test_CaputoL1point_accuracy_polynomial(self): - self.assertTrue( - abs( - CaputoL1point(0.5, lambda x: x**2 - 1, 0, 1.0, 1024) - - truevaluepoly_caputo - ) - <= 1e-3 - ) - - # exp - def test_CaputoL1point_accuracy_exp(self): - """Test RLpoint on f(x) = exp(x), alpha=0.5. Analytical: exp(1)/Gamma(0.5)""" - val = CaputoL1point(alpha_exp, np.exp, 0, 1, 1024) - self.assertTrue(abs(val - groundtruth_exp) < 0.6) # really bad accuracy - - - ####################### - # Caputo 2p - ####################### - - # x**2 - 1 - def test_CaputoL2point_accuracy_polynomial(self): - self.assertTrue( - abs( - CaputoL2point(1.5, lambda x: x**2, 0, 1.0, 1024) - - truevaluepoly_caputo_higher - ) - <= 1e-1 - ) - - ####################### - # Caputo 2pC - ####################### - - # x**2 (a = 1.5) - def test_CaputoL2Cpoint_accuracy_polynomial_higher(self): - self.assertTrue( - abs( - CaputoL2Cpoint(1.5, lambda x: x**2, 0, 1.0, 1024) - - truevaluepoly_caputo_higher - ) - <= 1e-1 - ) - - # x**2 - 1 - def test_CaputoL2Cpoint_accuracy_polynomial(self): - self.assertTrue( - abs( - CaputoL2Cpoint(0.5, lambda x: x**2, 0, 1.0, 1024) - truevaluepoly_caputo - ) - <= 1e-3 - ) - - # exp - def test_CaputoL2Cpoint_accuracy_exp(self): - """Test RLpoint on f(x) = exp(x), alpha=0.5. Analytical: exp(1)/Gamma(0.5)""" - val = CaputoL2Cpoint(alpha_exp, np.exp, 0, 1, 1024) - self.assertTrue(abs(val - groundtruth_exp) < 2) # unusable - - -class TestAlgorithmsGeneral(unittest.TestCase): - """Tests for correct size of algorithm results.""" - - def test_GL_result_length(self): - self.assertEqual(GL_length, test_N) - - def test_GLI_result_length(self): - self.assertEqual(GLI_length, test_N) - - def test_RL_result_length(self): - self.assertEqual(RL_length, test_N) - - def test_RL_matrix_shape(self): - self.assertTrue(np.shape(RLmatrix(0.4, test_N)) == (test_N, test_N)) - - def test_GL_binomial_coefficient_array_size(self): - self.assertEqual( - len(GLcoeffs(0.5, size_coefficient_array)) - 1, size_coefficient_array - ) - - - -if __name__ == "__main__": - unittest.main(argv=["first-arg-is-ignored"], exit=False) - - # Ensure all docstring examples work. - import doctest - - doctest.testmod() diff --git a/tests/test_functions.py b/tests/test_functions.py deleted file mode 100644 index 6c60a94..0000000 --- a/tests/test_functions.py +++ /dev/null @@ -1,73 +0,0 @@ -import unittest -import numpy as np - -# Import from sibling directory. -import sys -import os - -# Add the parent directory (containing differintP) to the import path -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) - -from differintP.functions import * - -# Define constants to be used in tests. -poch_first_argument = 1 -poch_second_argument = 5 -poch_true_answer = 120 - - -class HelperTestCases(unittest.TestCase): - """Tests for helper functions.""" - - def test_pochhammer(self): - self.assertEqual( - poch(poch_first_argument, poch_second_argument), poch_true_answer - ) - self.assertEqual(poch(-1, 3), 0) - # self.assertEqual(poch(-1.5, 0.5), np.inf) - self.assertEqual(np.round(poch(1j, 1), 3), 0.000 + 1.000j) - self.assertEqual(poch(-10, 2), 90) - - - """ Unit tests for Mittag-Leffler function. """ - - def test_ML_cosh_root(self): - xs = np.arange(10, 0.1) - self.assertTrue( - ( - np.abs( - MittagLeffler(2, 1, xs, ignore_special_cases=True) - - np.cosh(np.sqrt(xs)) - ) - <= 1e-3 - ).all() - ) - - def test_ML_exp(self): - xs = np.arange(10, 0.1) - self.assertTrue( - ( - np.abs(MittagLeffler(1, 1, xs, ignore_special_cases=True) - np.exp(xs)) - <= 1e-3 - ).all() - ) - - def test_ML_geometric(self): - xs = np.arange(1, 0.05) - self.assertTrue( - ( - np.abs( - MittagLeffler(0, 1, xs, ignore_special_cases=True) - 1 / (1 - xs) - ) - <= 1e-3 - ).all() - ) - - -if __name__ == "__main__": - unittest.main(argv=["first-arg-is-ignored"], exit=False) - - # Ensure all docstring examples work. - import doctest - - doctest.testmod() diff --git a/tests/test_solver.py b/tests/test_solver.py deleted file mode 100644 index c344577..0000000 --- a/tests/test_solver.py +++ /dev/null @@ -1,58 +0,0 @@ -import unittest -import numpy as np - -# Import from sibling directory. -import sys -import os - -# Add the parent directory (containing differintP) to the import path -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) - -from differintP.core import Gamma -from differintP.solvers import PCsolver -from differintP.functions import MittagLeffler - - -PC_x_power = np.linspace(0, 1, 100) ** 5.5 - -# Get FODE function for solving. -PC_func_power = lambda x, y: 1 / 24 * Gamma(5 + 1.5) * x**4 + x ** (8 + 2 * 1.5) - y**2 -PC_func_ML = lambda x, y: y - - -class TestSolvers(unittest.TestCase): - """Tests for the correct solution to the equations.""" - - def test_PC_solution_three_halves(self): - self.assertTrue( - ( - np.abs(PCsolver([0, 0], 1.5, PC_func_power, 0, 1, 100) - PC_x_power) - <= 1e-2 - ).all() - ) - - def test_PC_solution_ML(self): - xs = np.linspace(0, 1, 100) - ML_alpha = MittagLeffler(5.5, 1, xs**5.5) - self.assertTrue( - ( - np.abs(PCsolver([1, 0, 0, 0, 0, 0], 5.5, PC_func_ML) - ML_alpha) <= 1e-2 - ).all() - ) - - def test_PC_solution_linear(self): - xs = np.linspace(0, 1, 100) - self.assertTrue( - ( - np.abs(PCsolver([1, 1], 1.5, lambda x, y: y - x - 1) - (xs + 1)) <= 1e-2 - ).all() - ) - - -if __name__ == "__main__": - unittest.main(argv=["first-arg-is-ignored"], exit=False) - - # Ensure all docstring examples work. - import doctest - - doctest.testmod() diff --git a/tests/test_utils.py b/tests/test_utils.py deleted file mode 100644 index ab85d04..0000000 --- a/tests/test_utils.py +++ /dev/null @@ -1,79 +0,0 @@ -import unittest -import numpy as np - -# Import from sibling directory. -import sys -import os - -# Add the parent directory (containing differintP) to the import path -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) - -from differintP.utils import * # type: ignore -from .__init__ import test_N - - - - -stepsize = 1 / (test_N - 1) - -# Testing if callable functions and arrays of function values will work. -checked_function1, test_stepsize1 = functionCheck( - lambda x: 2 * np.exp(3 * x) * x - x**2 + x - 5, 0, 1, test_N -) -checked_function2, test_stepsize2 = functionCheck(np.ones(test_N), 0, 1, test_N) - - -class UtilsTestCases(unittest.TestCase): - def test_isInteger(self): - self.assertTrue(isInteger(1)) - self.assertTrue(isInteger(1.0)) - self.assertTrue(isInteger(1 + 0j)) - self.assertFalse(isInteger(1.1)) - self.assertFalse(isInteger(1.1 + 0j)) - self.assertFalse(isInteger(1 + 1j)) - - def test_isPositiveInteger(self): - self.assertTrue(isPositiveInteger(1)) - self.assertFalse(isPositiveInteger(1.1)) - self.assertFalse(isPositiveInteger(-1)) - - def test_functionCheck(self): - self.assertEqual(len(checked_function1), test_N) - self.assertEqual(len(checked_function2), test_N) - - # Make sure it treats defined functions and arrays of function values the same. - self.assertEqual(len(checked_function1), len(checked_function2)) - self.assertEqual(test_stepsize1, stepsize) - self.assertEqual(test_stepsize2, stepsize) - self.assertEqual(test_stepsize1, test_stepsize2) - - - def test_checkValues(self): - with self.assertRaises(AssertionError): - checkValues(0.1, 0, 1, 1.1) # type: ignore - with self.assertRaises(AssertionError): - checkValues(0.1, 1j, 2, 100) # type: ignore - with self.assertRaises(AssertionError): - 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) # type: ignore - checkValues(0.5, 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]) - [ - [ - [ - [ - checkValues(alpha, domain_start, domain_end, num_points) - for alpha in alpha_vals - ] - for domain_start in domain_vals - ] - for domain_end in domain_vals - ] - for num_points in num_vals - ] \ No newline at end of file From c3ee27a70b812d3d1424cb9f3e5a960d7343becd Mon Sep 17 00:00:00 2001 From: parsa roshanak Date: Mon, 7 Jul 2025 21:51:44 +0330 Subject: [PATCH 2/2] 0.0.4 c2 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index bb75657..568534e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "differintP" -version = "0.0.3" +version = "0.0.4" description = "Modern, pure Python fractional calculus library (fork of differint)" readme = "README.md" requires-python = ">=3.10"