0.0.4

.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]

changelog.md

@@ -156,4 +156,16 @@
   - `RL` `RLpoint` `PCsolver` functions Namespaces
 
 - Aditional test and test restructuring
-- Restructuring the Namespaces
+- 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.

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"
 ]

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
+
+

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"

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()