diff --git a/.gitignore b/.gitignore index 34c869e..9216ac1 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ bentchmark.ipynb +test_gpu.py # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/accuracy_test.ipynb b/accuracy_test.ipynb new file mode 100644 index 0000000..68bdf6b --- /dev/null +++ b/accuracy_test.ipynb @@ -0,0 +1,552 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "900c2a2e", + "metadata": {}, + "source": [ + "# Accuracy Test" + ] + }, + { + "cell_type": "markdown", + "id": "0f2e8118", + "metadata": {}, + "source": [ + "## 📦 **Fractional Derivative Package Comparison**\n", + "\n", + "### **Test Function:**\n", + "\n", + "* **Function:** `sin(x)`\n", + "* **Order α:** 0.5\n", + "* **Domain:** \\[0.0, 1.57]\n", + "* **Points:** 1024\n", + "\n", + "---\n", + "\n", + "### **Absolute Error Comparison**\n", + "\n", + "| Method | differintP\\_err | differint\\_err | differintC\\_err |\n", + "| :-------------- | --------------: | -------------: | --------------: |\n", + "| GL | 0.0871 | 0.0871 | 0.0871 |\n", + "| GLpoint | 0.0871 | 0.0868 | 0.0871 |\n", + "| GLI | 0.0853 | 0.0853 | — |\n", + "| GLpoint\\_direct | 0.0868 | — | — |\n", + "| RL | 0.0873 | 0.0873 | 0.0873 |\n", + "| RLpoint | 0.0873 | 0.0873 | 0.0873 |\n", + "| CaputoL1point | 0.0873 | 0.0873 | — |\n", + "\n", + "*Note:* `—` indicates not available in that package.\n", + "\n", + "---\n", + "\n", + "### **Package Versions**\n", + "\n", + "| Package | Version |\n", + "| :--------- | :------ |\n", + "| differintP | 0.0.3 |\n", + "| differint | 1.0.0 |\n", + "| differintC | 0.0.2.4 |" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "3d1f4239", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from differintP import *\n", + "from differintP.special import GLpoint_direct\n", + "import pandas as pd\n", + "\n", + "import differintP as dp\n", + "\n", + "# Legacy (pip install differint)\n", + "import differint.differint as d\n", + "\n", + "# C-accelerated (if you have it built/installed)\n", + "# pip install differintC\n", + "import differintC as dc" + ] + }, + { + "cell_type": "markdown", + "id": "6024435a", + "metadata": {}, + "source": [ + "## Helper Functions" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "de024fcb", + "metadata": {}, + "outputs": [], + "source": [ + "def accuracy_point(\n", + " gt, # ground truth\n", + " function,\n", + " alpha,\n", + " test_subject,\n", + " start: float = 0.0,\n", + " end: float = 1.0,\n", + " num_points: int = 1024,\n", + "):\n", + " res = test_subject(alpha, function, start, end, num_points)\n", + " abs_e = abs(gt - res)\n", + " r2_e = (gt - res) ** 2\n", + "\n", + " return {\n", + " \"name\" : test_subject.__name__,\n", + " \"abs_e\" : abs_e,\n", + " \"r2_e\" : r2_e,\n", + " \"gt\" : gt,\n", + " \"res\": res,\n", + " }\n", + "\n", + "def accuracy_point_from_array(\n", + " gt, # ground truth\n", + " function,\n", + " alpha,\n", + " test_subject,\n", + " start: float = 0.0,\n", + " end: float = 1.0,\n", + " num_points: int = 1024,\n", + "):\n", + " res = test_subject(alpha, function, start, end, num_points)[-1]\n", + " abs_e = abs(gt - res)\n", + " r2_e = (gt - res) ** 2\n", + "\n", + " return {\n", + " \"name\" : test_subject.__name__,\n", + " \"abs_e\" : abs_e,\n", + " \"r2_e\" : r2_e,\n", + " \"gt\" : gt,\n", + " \"res\": res,\n", + " }\n", + "\n", + "def accuracy_array(\n", + " gt, # ground truth function\n", + " function,\n", + " alpha,\n", + " test_subject,\n", + " start: float = 0.0,\n", + " end: float = 1.0,\n", + " num_points: int = 1024,\n", + "):\n", + " x = np.linspace(start, end, num_points)\n", + " gt_val = gt(x)\n", + " res = test_subject(alpha, function, start, end, num_points)[-1]\n", + "\n", + " abs_e = abs(gt - gt_val)\n", + " avg_abs_e = np.mean(abs_e)\n", + " max_abs_e = np.max(avg_abs_e)\n", + "\n", + " return {\n", + " \"name\" : test_subject.__name__,\n", + " \"avg_abs_e\" : avg_abs_e,\n", + " \"max_abs_e\" : max_abs_e,\n", + " \"gt_val\" : gt_val,\n", + " \"res\": res,\n", + " }\n", + "\n", + "\n", + "## Comparison helpers\n", + "\n", + "def try_array(gt, f, alpha, func, *args, **kwargs):\n", + " if func is None:\n", + " return None\n", + " try:\n", + " res = func(alpha, f, *args, **kwargs)\n", + " val = res[-1]\n", + " abs_e = abs(gt - val)\n", + " return val, abs_e\n", + " except Exception as e:\n", + " return str(e), None\n", + "\n", + "def try_point(gt, f, alpha, func, *args, **kwargs):\n", + " if func is None:\n", + " return None\n", + " try:\n", + " val = func(alpha, f, *args, **kwargs)\n", + " abs_e = abs(gt - val)\n", + " return val, abs_e\n", + " except Exception as e:\n", + " return str(e), None\n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "14fe6d7c", + "metadata": {}, + "source": [ + "## Sin Test" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "c0477010", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ground truth (analytical): sin(1.57 + 0.5 * pi / 2) = 0.7071067811865476\n" + ] + } + ], + "source": [ + "# 1. Setup: Parameters and Ground Truth\n", + "\n", + "# Parameters for the test\n", + "alpha = 0.5\n", + "start = 0.0\n", + "end = np.pi / 2 # Let's pick [0, pi/2] for better nonzero behavior\n", + "num_points = 1024\n", + "\n", + "\n", + "# Test function\n", + "f_name = \"sin\"\n", + "f = np.sin\n", + "gt_expr = f\"sin({end:.2f} + {alpha} * pi / 2)\"\n", + "gt = np.sin(end + alpha * np.pi / 2)\n", + "\n", + "\n", + "\n", + "print(\"Ground truth (analytical):\", gt_expr, \"=\", gt)" + ] + }, + { + "cell_type": "code", + "execution_count": 65, + "id": "8be640e8", + "metadata": {}, + "outputs": [], + "source": [ + "# 2.A Run All Methods on This Test\n", + "\n", + "methods_point = [GLpoint, GLpoint_direct, RLpoint, CaputoL1point]\n", + "methods_array = [GL, GLI, RL]\n", + "\n", + "results = []\n", + "\n", + "# Test point functions\n", + "for method in methods_point:\n", + " result = accuracy_point(\n", + " gt=gt,\n", + " function=f,\n", + " alpha=alpha,\n", + " test_subject=method,\n", + " start=start,\n", + " end=end,\n", + " num_points=num_points,\n", + " )\n", + " results.append(result)\n", + "\n", + "# Test array functions (use last value)\n", + "for method in methods_array:\n", + " result = accuracy_point_from_array(\n", + " gt=gt,\n", + " function=f,\n", + " alpha=alpha,\n", + " test_subject=method,\n", + " start=start,\n", + " end=end,\n", + " num_points=num_points,\n", + " )\n", + " results.append(result)\n", + "df = pd.DataFrame(results)[[\"name\", \"abs_e\", \"r2_e\", \"gt\", \"res\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "id": "17aaa300", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Setup:\n", + "a: 0.5 - start: 0.0 - end: 1.57 - points: 1024 - test: sin\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
nameabs_er2_egtres
0GLpoint0.0870640.0075800.7071070.620043
1GLpoint_direct0.0867610.0075280.7071070.620346
2RLpoint0.0873010.0076210.7071070.619806
3CaputoL1point0.0873010.0076210.7071070.619806
4GL0.0870640.0075800.7071070.620043
5GLI0.0853140.0072780.7071070.621793
6RL0.0873010.0076210.7071070.619806
\n", + "
" + ], + "text/plain": [ + " name abs_e r2_e gt res\n", + "0 GLpoint 0.087064 0.007580 0.707107 0.620043\n", + "1 GLpoint_direct 0.086761 0.007528 0.707107 0.620346\n", + "2 RLpoint 0.087301 0.007621 0.707107 0.619806\n", + "3 CaputoL1point 0.087301 0.007621 0.707107 0.619806\n", + "4 GL 0.087064 0.007580 0.707107 0.620043\n", + "5 GLI 0.085314 0.007278 0.707107 0.621793\n", + "6 RL 0.087301 0.007621 0.707107 0.619806" + ] + }, + "execution_count": 66, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# 3.B Display Results as a Table\n", + "print(\"Setup:\")\n", + "print(f\"a: {alpha} - start: {start} - end: {end:.2f} - points: {num_points} - test: {f_name}\")\n", + "\n", + "df\n", + "# print(df.to_markdown(index=False))" + ] + }, + { + "cell_type": "markdown", + "id": "38e14723", + "metadata": {}, + "source": [ + "### Package Comparison" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "id": "2a1aa976", + "metadata": {}, + "outputs": [], + "source": [ + "# 2.B Intra Package Test\n", + "test_subjects = [\n", + " ('GL', dp.GL, d.GL, dc.GL),\n", + " ('GLpoint', dp.GLpoint, d.GLpoint, dc.GLpoint),\n", + " ('GLI', dp.GLI, d.GLI, None),\n", + " ('GLpoint_direct',GLpoint_direct, None, None),\n", + " ('RL', dp.RL, d.RL, dc.RL),\n", + " ('RLpoint', dp.RLpoint, d.RLpoint, dc.RLpoint),\n", + " ('CaputoL1point', dp.CaputoL1point, d.CaputoL1point, None),\n", + "]\n", + "\n", + "table = []\n", + "for name, f_p, f_l, f_c in test_subjects:\n", + " row = {'Method': name}\n", + " # Choose point or array dispatcher\n", + " dispatcher = try_point if 'point' in name else try_array\n", + " for tag, fn in zip(['differintP', 'differint', 'differintC'], [f_p, f_l, f_c]):\n", + " result = dispatcher(gt, f, alpha, fn, start, end, num_points)\n", + " if result is None:\n", + " row[tag] = ''\n", + " else:\n", + " val, abs_e = result\n", + " row[tag] = f\"{val:.8g}\" if abs_e is not None else val\n", + " row[tag + \"_err\"] = f\"{abs_e:.2e}\" if abs_e is not None else ''\n", + " table.append(row)\n", + "df_intra = pd.DataFrame(table)\n", + "\n", + "err_cols = [col for col in df_intra.columns if col.strip().endswith('_err')]\n", + "if 'Method' in df_intra.columns:\n", + " err_cols = ['Method'] + err_cols\n", + "df_errors = df_intra[err_cols]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ab299c7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Setup:\n", + "a: 0.5 - start: 0.0 - end: 1.57 - points: 1024 - test: sin\n", + "| Method | differintP_err | differint_err | differintC_err |\n", + "|:---------------|-----------------:|----------------:|-----------------:|\n", + "| GL | 0.0871 | 0.0871 | 0.0871 |\n", + "| GLpoint | 0.0871 | 0.0868 | 0.0871 |\n", + "| GLI | 0.0853 | 0.0853 | nan |\n", + "| GLpoint_direct | 0.0868 | nan | nan |\n", + "| RL | 0.0873 | 0.0873 | 0.0873 |\n", + "| RLpoint | 0.0873 | 0.0873 | 0.0873 |\n", + "| CaputoL1point | 0.0873 | 0.0873 | nan |\n" + ] + } + ], + "source": [ + "# 3.B Display Results as a Table\n", + "print(\"Setup:\")\n", + "print(f\"a: {alpha} - start: {start} - end: {end:.2f} - points: {num_points} - test: {f_name}\")\n", + "\n", + "df_errors\n", + "# print(df_errors.to_markdown(index=False))\n" + ] + }, + { + "cell_type": "markdown", + "id": "eb2da7bf", + "metadata": {}, + "source": [ + "## Package Comparison:\n", + "\n", + "| Method | differintP_err | differint_err | differintC_err |\n", + "|:---------------|-----------------:|----------------:|-----------------:|\n", + "| GL | 0.0871 | 0.0871 | 0.0871 |\n", + "| GLpoint | 0.0871 | + 0.0868 | 0.0871 |\n", + "| GLI | 0.0853 | 0.0853 | nan |\n", + "| GLpoint_direct | 0.0868 | nan | nan |\n", + "| RL | 0.0873 | 0.0873 | 0.0873 |\n", + "| RLpoint | 0.0873 | 0.0873 | 0.0873 |\n", + "| CaputoL1point | 0.0873 | 0.0873 | nan |\n", + "\n", + "- Setup:\n", + "- a: 0.5 - start: 0.0 - end: 1.57 - points: 1024 - test: sin\n", + "\n", + "* Versions\n", + "\n", + "| Package | differintP | differint | differintC_err |\n", + "|:---------------|-----------------:|----------------:|-----------------:|\n", + "| Version | 0.0.3 | 1.0.0 | 0.0.2.4 |" + ] + }, + { + "cell_type": "markdown", + "id": "0f0fa328", + "metadata": {}, + "source": [ + "## differintP result\n", + "\n", + "| name | abs_e | r2_e | gt | res |\n", + "|:---------------|----------:|-----------:|---------:|---------:|\n", + "| GLpoint | 0.0870642 | 0.00758018 | 0.707107 | 0.620043 |\n", + "| GLpoint_direct | 0.0867613 | 0.00752752 | 0.707107 | 0.620346 |\n", + "| RLpoint | 0.0873008 | 0.00762142 | 0.707107 | 0.619806 |\n", + "| CaputoL1point | 0.0873008 | 0.00762142 | 0.707107 | 0.619806 |\n", + "| GL | 0.0870642 | 0.00758018 | 0.707107 | 0.620043 |\n", + "| GLI | 0.0853137 | 0.00727842 | 0.707107 | 0.621793 |\n", + "| RL | 0.0873008 | 0.00762142 | 0.707107 | 0.619806 |\n", + "\n", + "v0.0.3 results:\n", + "a: 0.5 - start: 0.0 - end: 1.57 - points: 1024 - test: sin" + ] + }, + { + "cell_type": "markdown", + "id": "b19ed1b2", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/changelog.md b/changelog.md index 62a8a74..4e816a9 100644 --- a/changelog.md +++ b/changelog.md @@ -116,8 +116,10 @@ ### 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. +* **Redesigned `GLpoint`:** + Replaced the previous GLpoint implementation with a highly efficient, single-pass C++-style recurrence. The new method uses a Numba JIT-compiled kernel for substantial speed gains, directly computing the Grünwald-Letnikov fractional derivative at the endpoint in one loop over the function values. + Legacy and alternative endpoint methods (`GLpoint_direct`, `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: @@ -131,23 +133,27 @@ | 1+e7 | 51.2645 ms | 3674.3541 ms | x71.67 | -- **Optimized `GLI` implementation:** +* **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. + * Rewrote core algorithm to use a Numba-accelerated helper for the main rolling-window convolution, matching the original literature. + * Ensured correctness by flipping GL coefficients within the moving-window convolution. + * Replaced class-based coefficient calculation with direct variable computation for improved clarity and efficiency. + * All critical loops are now compiled with Numba for major speed improvements, especially at high grid resolutions. + * Modernized function interface and added robust array conversion and shape checks. [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 | +| 1+e2 | 0.0760 ms | 0.7582 ms | x6.06 | +| 1+e3 | 2.0336 ms | 7.391 ms | x3.63 | +| 1+e4 | 155.7649 ms | 128.3709 ms | x0.8 | +| 5+e4 | 3895.176 ms | 15181.82 ms | x3.89 | - Wiki Pages for - `CaputoL1point` `CaputoL2point` `CaputoL2Cpoint` `CaputoFromRLpoint` - `CRONE` `GL` `GLpoint` `GLI` `GL_gpu` - - `RL` `RLpoint` + - `RL` `RLpoint` `PCsolver` functions Namespaces + +- Aditional test and test restructuring +- Restructuring the Namespaces \ No newline at end of file diff --git a/differintP/__init__.py b/differintP/__init__.py index a5c68e5..95a5b79 100644 --- a/differintP/__init__.py +++ b/differintP/__init__.py @@ -1,15 +1,37 @@ -from .core import * +from .core import ( + GLcoeffs, + GL, + GLpoint, + GL_gpu, + GLI, + CRONE, + RLmatrix, + RLcoeffs, + RLpoint, + RL, + CaputoL1point, + CaputoL2point, + CaputoL2Cpoint, + CaputoFromRLpoint, +) + +from . import functions + __all__ = [ - "GL", "GLcoeffs", - "functionCheck", - "checkValues", - "isInteger", - "GLI", + "GL", "GLpoint", - "RLcoeffs", + "GL_gpu", + "GLI", + "CRONE", "RLmatrix", - "RL", - "poch", + "RLcoeffs", "RLpoint", + "RL", + "CaputoL1point", + "CaputoL2point", + "CaputoL2Cpoint", + "CaputoFromRLpoint", + + "functions" ] diff --git a/differintP/core.py b/differintP/core.py index 0e837eb..03e0b91 100644 --- a/differintP/core.py +++ b/differintP/core.py @@ -1,161 +1,19 @@ from __future__ import print_function -from typing import Callable, cast +from typing import Callable 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 +from .utils import checkValues, functionCheck -def isInteger(n) -> bool: - if n.imag: - return False - if float(n.real).is_integer(): - return True - else: - return False - - -def isPositiveInteger(n) -> bool: - return isInteger(n) and n > 0 - - -def checkValues( - alpha: float, - domain_start: int | float, - domain_end: int | float, - num_points: int, - support_complex_alpha: bool = False, -) -> bool | None: - """Type checking for valid inputs.""" - - assert isPositiveInteger(num_points), ( - "num_points is not an integer: %r" % num_points - ) - - assert isinstance(domain_start, (int, np.integer, float, np.floating)), ( - "domain_start must be integer or float: %r" % domain_start - ) - assert isinstance(domain_end, (int, np.integer, float, np.floating)), ( - "domain_end must be integer or float: %r" % domain_end - ) - if not support_complex_alpha: - assert not isinstance( - alpha, complex - ), "Complex alpha not supported for this algorithm." - - return - - -def functionCheck( - f_name: Callable | list | np.ndarray, - domain_start: float | int, - domain_end: float | int, - num_points: int, -): - """Check if function is callable and assign function values.""" - - # Define the function domain and obtain function values. - if hasattr(f_name, "__call__"): - f_name = cast(Callable, f_name) - # If f_name is callable, call it and save to a list. - x = np.linspace(domain_start, domain_end, num_points) - f_values = list(map(lambda t: f_name(t), x)) - step_size = x[1] - x[0] - else: - 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) - return f_values, step_size - - -def poch(a, n): - """Returns the Pochhammer symbol (a)_n. a can be any complex or real number - except the negative integers and 0. n can be any nonnegative real. - """ - if isPositiveInteger(n): - # Compute the Pochhammer symbol. - n = int(n) - if n == 0: - return 1.0 - else: - poch = 1 - for j in range(n): - poch *= a + j - return poch - - # if a and a + n are both nonpositive integers, we can use another formula... - # see here https://www.mathworks.com/help/symbolic/sym.pochhammer.html - if isPositiveInteger(-1 * a) and isPositiveInteger(-1 * a - n): - sign = -1 if np.abs(n % 2) == 1 else 1 - return sign * Gamma(1 - a) / Gamma(1 - a - n) - return Gamma(a + n) / Gamma(a) - - -def Beta( - x: int | float | np.ndarray | list | complex, - y: int | float | np.ndarray | list | complex, -) -> int | float | np.ndarray | list | complex: - """Beta function using Godfrey's Gamma function.""" - - return Gamma(x) * Gamma(y) / Gamma(x + y) # type: ignore - - -def MittagLeffler( - a: float, - b: float, - x: np.ndarray, - num_terms: int = 50, - ignore_special_cases: bool = False, -) -> np.ndarray: - """Calculate the Mittag-Leffler function by checking for special cases, and trying to - reduce the parameters. If neither of those work, it just brute forces it. - - Parameters - ========== - a : float - The first parameter of the Mittag-Leffler function. - b : float - The second parameter of the Mittag-Leffler function - x : 1D-array of floats (can be len = 1) - The value or values to be evaluated at. - num_terms : int - The number of terms to calculate in the sum. Ignored if - a special case can be used instead. Default value is 100. - ignore_special_cases : bool - Don't use the special cases, use the series definition. - Probably only useful for testing. Default value is False. - """ - # check for quick special cases - if not ignore_special_cases: - if a == 0: - if (np.abs(x) < 1).all(): - return 1 / Gamma(b) * 1 / (1 - x) - return x * np.inf - elif a == 0.5 and b == 1: - # requires calculation of the complementary error function - pass - elif a == 1 and b == 1: - return np.exp(x) - elif a == 2 and b == 1: - return np.cosh(np.sqrt(x)) - elif a == 1 and b == 2: - return (np.exp(x) - 1) / x - elif a == 2 and b == 2: - return np.sinh(np.sqrt(x)) / np.sqrt(x) - # manually calculate with series definition - exponents = np.arange(num_terms) - exp_vals = np.array([x]).T ** exponents - gamma_vals = np.array([Gamma(exponent * a + b) for exponent in exponents]) - return np.sum(exp_vals / gamma_vals, axis=1) ######################################################################################### @@ -252,6 +110,7 @@ def _GLpoint_loop(alpha: float, f_values: np.ndarray, step: float) -> float: 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, @@ -301,6 +160,7 @@ def GLpoint( ######################################## GL Gpu ######################################### ######################################################################################### + def _gpu_GLcoeffs( alpha: float, n: int, @@ -358,22 +218,43 @@ 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. +@njit +def _GLI_core(alpha: float, f_values: np.ndarray, b_coeffs: np.ndarray) -> np.ndarray: """ + Developer function: Efficiently computes the improved Grünwald-Letnikov (GLI) fractional derivative core. + + Given function values and GL coefficients (as 1D NumPy arrays), this routine: + - Applies a 3-point quadratic (Lagrange) interpolation in a moving window. + - For each point, performs a local convolution with a flipped GL coefficient array. + - Combines results with fixed interpolation weights for high-order accuracy. - 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 + Intended for internal use with Numba JIT compilation for speed. + Inputs must be 1D NumPy arrays of matching length. + """ + num_points = f_values.shape[0] + GLI = np.zeros(num_points) + prv = alpha * (alpha - 2) / 8 + crr = (4 - alpha * alpha) / 4 + nxt = alpha * (2 + alpha) / 8 + I = np.array([prv, crr, nxt]) + for i in range(3, num_points): + L = i + F = f_values[:L] + B = b_coeffs[: (L - 2)] + G = np.zeros(3) + B_flip = B[::-1] + for m in range(3): + s = 0.0 + for n in range(len(B)): + s += F[m + n] * B_flip[n] + G[m] = s + GLI[i] = np.sum(G * I) + return GLI def GLI( @@ -385,15 +266,14 @@ def GLI( ) -> np.ndarray: """ Computes the 'improved' Grünwald-Letnikov (GL) fractional derivative of a function over an entire domain, - using a hybrid approach for efficiency. + using the quadratic 3-point Lagrange interpolation method described by Oldham & Spanier (1974). - 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. + This implementation applies a three-point moving-window interpolation to the function values, + followed by a specialized convolution with Grünwald-Letnikov coefficients at each evaluation point. + The result is a higher-order, more accurate fractional derivative estimate compared to the standard GL method, + especially for smooth functions. - 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`. + For performance, a Numba-accelerated helper is used internally. Parameters ---------- @@ -434,37 +314,24 @@ def GLI( >>> 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. 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, step_size = functionCheck(f_name, domain_start, domain_end, num_points) - - f_values = np.asarray(f_values) - - # Get interpolating values. - prv = alpha * (alpha - 2) / 8 - crr = (4 - alpha * alpha) / 4 - nxt = alpha * (2 + alpha) / 8 + # Evaluate function on grid + x = np.linspace(domain_start, domain_end, num_points) + step_size = (domain_end - domain_start) / (num_points - 1) + if callable(f_name): + f_values = np.asarray(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") - # 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 + b_coeffs = GLcoeffs(alpha, num_points) + GLI_vals = _GLI_core(alpha, f_values, b_coeffs) + return GLI_vals * step_size**-alpha - # Precompute coefficients - b_coeffs = GLcoeffs(alpha, num_points - 1) - - # 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): @@ -535,9 +402,10 @@ def CRONEfilter(siz, alpha): ######################################################################################### -########################################### RL ########################################## +######################################### RLmatrix ###################################### ######################################################################################### + def RLmatrix(alpha, N): """Vectorized RL coefficient matrix generation""" # Precompute all required powers @@ -566,6 +434,9 @@ def RLmatrix(alpha, N): # Normalize with Gamma function return coeffMatrix / Gamma(2 - alpha) +######################################################################################### +####################################### RLcoeffs ######################################## +######################################################################################### def RLcoeffs(index_k, index_j, alpha): """Calculates coefficients for the RL differintegral operator. @@ -605,6 +476,9 @@ def RLcoeffs(index_k, index_j, alpha): - 2 * (index_k - index_j) ** (1 - alpha) ) +######################################################################################### +######################################## RLpoint ######################################## +######################################################################################### def RLpoint( alpha: float, @@ -678,6 +552,9 @@ def RLpoint( C = 1 / Gamma(2 - alpha) return C * step_size**-alpha * np.dot(coeffs, f_values) +######################################################################################### +########################################### RL ########################################## +######################################################################################### def RL( alpha: float, @@ -737,10 +614,6 @@ def RL( return result - - - - ######################################################################################### ######################################### Caputo ######################################## ######################################################################################### @@ -1005,125 +878,3 @@ def CaputoFromRLpoint( return C -######################################################################################### -######################################## PCsolver ####################################### -######################################################################################### - -def PCcoeffs(alpha, j, n): - if 1 < alpha: - if j == 0: - return ( - (n + 1) ** alpha * (alpha - n) - + n**alpha * (2 * n - alpha - 1) - - (n - 1) ** (alpha + 1) - ) - elif j == n: - return 2 ** (alpha + 1) - alpha - 3 - return ( - (n - j + 2) ** (alpha + 1) - + 3 * (n - j) ** (alpha + 1) - - 3 * (n - j + 1) ** (alpha + 1) - - (n - j - 1) ** (alpha + 1) - ) - - -def PCsolver( - initial_values, alpha, f_name, domain_start=0, domain_end=1, num_points=100 -): - """Solve an equation of the form D[y(x)]=f(x, y(x)) using the predictor-corrector - method, modified to be compatible with fractional derivatives. - - see Deng, W. (2007) Short memory principle and a predictor-corrector approach for - fractional differential equations. Journal of Computational and Applied - Mathematics. - - test examples from - Baskonus, H.M., Bulut, H. (2015) On the numerical solutions of some fractional - ordinary differential equations by fractional Adams-Bashforth-Moulton method. - De Gruyter. - Weilbeer, M. (2005) Efficient Numerical Methods for Fractional Differential - Equations and their Analytical Background. - - Parameters - ========== - initial_values : float 1d-array - A list of initial values for the IVP. There should be as many IVs - as ceil(alpha). - alpha : float - The order of the differintegral in the equation to be computed. - f_name : function handle or lambda function - This is the function on the right side of the equation, and should - accept two variables; first the independant variable, and second - the equation to be solved. - 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. - - Output - ====== - y_correction : float 1d-array - The calculated solution to the IVP at each of the points - between the left and right endpoint. - - Examples: - >>> f_name = lambda x, y : y - x - 1 - >>> initial_values = [1, 1] - >>> y_solved = PCsolver(initial_values, 1.5, f_name) - >>> theoretical = np.linspace(0, 1, 100) + 1 - >>> np.allclose(y_solved, theoretical) - True - """ - from scipy.special import factorial - - x_points = np.linspace(domain_start, domain_end, num_points) - step_size = x_points[1] - x_points[0] - y_correction = np.zeros(num_points, dtype="complex") - y_prediction = np.zeros(num_points, dtype="complex") - - y_prediction[0] = initial_values[0] - y_correction[0] = initial_values[0] - for x_index in range(num_points - 1): - initial_value_contribution = 0 - if 1 < alpha and alpha < 2: - initial_value_contribution = initial_values[1] * step_size - elif 2 < alpha: - for k in range(1, int(np.ceil(alpha))): - initial_value_contribution += ( - initial_values[k] - / factorial(k) - * (x_points[x_index + 1] ** k - x_points[x_index] ** k) - ) - elif alpha < 1: - raise ValueError("Not yet supported!") - y_prediction[x_index + 1] += initial_value_contribution - y_prediction[x_index + 1] += y_correction[x_index] - y_prediction[x_index + 1] += ( - step_size**alpha - / Gamma(alpha + 1) - * f_name(x_points[x_index], y_correction[x_index]) - ) - subsum = 0 - for j in range(x_index + 1): - subsum += PCcoeffs(alpha, j, x_index) * f_name(x_points[j], y_correction[j]) - y_prediction[x_index + 1] += step_size**alpha / Gamma(alpha + 2) * subsum - - y_correction[x_index + 1] += initial_value_contribution - y_correction[x_index + 1] += y_correction[x_index] - y_correction[x_index + 1] += ( - step_size**alpha - / Gamma(alpha + 2) - * alpha - * f_name(x_points[x_index], y_correction[x_index]) - ) - y_correction[x_index + 1] += ( - step_size**alpha - / Gamma(alpha + 2) - * f_name(x_points[x_index + 1], y_prediction[x_index + 1]) - ) - y_correction[x_index + 1] += step_size**alpha / Gamma(alpha + 2) * subsum - - return y_correction diff --git a/differintP/functions.py b/differintP/functions.py new file mode 100644 index 0000000..94734c4 --- /dev/null +++ b/differintP/functions.py @@ -0,0 +1,87 @@ +import numpy as np + +from scipy.special import gamma + +from .utils import isPositiveInteger + + + +def poch(a, n): + """Returns the Pochhammer symbol (a)_n. a can be any complex or real number + except the negative integers and 0. n can be any nonnegative real. + """ + if isPositiveInteger(n): + # Compute the Pochhammer symbol. + n = int(n) + if n == 0: + return 1.0 + else: + poch = 1 + for j in range(n): + poch *= a + j + return poch + + # if a and a + n are both nonpositive integers, we can use another formula... + # see here https://www.mathworks.com/help/symbolic/sym.pochhammer.html + if isPositiveInteger(-1 * a) and isPositiveInteger(-1 * a - n): + sign = -1 if np.abs(n % 2) == 1 else 1 + return sign * gamma(1 - a) / gamma(1 - a - n) + return gamma(a + n) / gamma(a) + + +def Beta( + x: int | float | np.ndarray | complex, + y: int | float | np.ndarray | complex, +) -> int | float | np.ndarray | complex: + """Beta function using Scipy Gamma function.""" + + return gamma(x) * gamma(y) / gamma(x + y) + + +def MittagLeffler( + a: float, + b: float, + x: np.ndarray, + num_terms: int = 50, + ignore_special_cases: bool = False, +) -> np.ndarray: + """Calculate the Mittag-Leffler function by checking for special cases, and trying to + reduce the parameters. If neither of those work, it just brute forces it. + + Parameters + ========== + a : float + The first parameter of the Mittag-Leffler function. + b : float + The second parameter of the Mittag-Leffler function + x : 1D-array of floats (can be len = 1) + The value or values to be evaluated at. + num_terms : int + The number of terms to calculate in the sum. Ignored if + a special case can be used instead. Default value is 100. + ignore_special_cases : bool + Don't use the special cases, use the series definition. + Probably only useful for testing. Default value is False. + """ + # check for quick special cases + if not ignore_special_cases: + if a == 0: + if (np.abs(x) < 1).all(): + return 1 / gamma(b) * 1 / (1 - x) + return x * np.inf + elif a == 0.5 and b == 1: + # requires calculation of the complementary error function + pass + elif a == 1 and b == 1: + return np.exp(x) + elif a == 2 and b == 1: + return np.cosh(np.sqrt(x)) + elif a == 1 and b == 2: + return (np.exp(x) - 1) / x + elif a == 2 and b == 2: + return np.sinh(np.sqrt(x)) / np.sqrt(x) + # manually calculate with series definition + exponents = np.arange(num_terms) + exp_vals = np.array([x]).T ** exponents + gamma_vals = np.array([gamma(exponent * a + b) for exponent in exponents]) + return np.sum(exp_vals / gamma_vals, axis=1) \ No newline at end of file diff --git a/differintP/solvers.py b/differintP/solvers.py new file mode 100644 index 0000000..90bc02f --- /dev/null +++ b/differintP/solvers.py @@ -0,0 +1,132 @@ +from __future__ import print_function + + +import numpy as np + +from scipy.special import gamma as Gamma + + + +######################################################################################### +######################################## PCsolver ####################################### +######################################################################################### + + +def PCcoeffs(alpha, j, n): + if 1 < alpha: + if j == 0: + return ( + (n + 1) ** alpha * (alpha - n) + + n**alpha * (2 * n - alpha - 1) + - (n - 1) ** (alpha + 1) + ) + elif j == n: + return 2 ** (alpha + 1) - alpha - 3 + return ( + (n - j + 2) ** (alpha + 1) + + 3 * (n - j) ** (alpha + 1) + - 3 * (n - j + 1) ** (alpha + 1) + - (n - j - 1) ** (alpha + 1) + ) + + +def PCsolver( + initial_values, alpha, f_name, domain_start=0, domain_end=1, num_points=100 +): + """Solve an equation of the form D[y(x)]=f(x, y(x)) using the predictor-corrector + method, modified to be compatible with fractional derivatives. + + see Deng, W. (2007) Short memory principle and a predictor-corrector approach for + fractional differential equations. Journal of Computational and Applied + Mathematics. + + test examples from + Baskonus, H.M., Bulut, H. (2015) On the numerical solutions of some fractional + ordinary differential equations by fractional Adams-Bashforth-Moulton method. + De Gruyter. + Weilbeer, M. (2005) Efficient Numerical Methods for Fractional Differential + Equations and their Analytical Background. + + Parameters + ========== + initial_values : float 1d-array + A list of initial values for the IVP. There should be as many IVs + as ceil(alpha). + alpha : float + The order of the differintegral in the equation to be computed. + f_name : function handle or lambda function + This is the function on the right side of the equation, and should + accept two variables; first the independant variable, and second + the equation to be solved. + 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. + + Output + ====== + y_correction : float 1d-array + The calculated solution to the IVP at each of the points + between the left and right endpoint. + + Examples: + >>> f_name = lambda x, y : y - x - 1 + >>> initial_values = [1, 1] + >>> y_solved = PCsolver(initial_values, 1.5, f_name) + >>> theoretical = np.linspace(0, 1, 100) + 1 + >>> np.allclose(y_solved, theoretical) + True + """ + from scipy.special import factorial + + x_points = np.linspace(domain_start, domain_end, num_points) + step_size = x_points[1] - x_points[0] + y_correction = np.zeros(num_points, dtype="complex") + y_prediction = np.zeros(num_points, dtype="complex") + + y_prediction[0] = initial_values[0] + y_correction[0] = initial_values[0] + for x_index in range(num_points - 1): + initial_value_contribution = 0 + if 1 < alpha and alpha < 2: + initial_value_contribution = initial_values[1] * step_size + elif 2 < alpha: + for k in range(1, int(np.ceil(alpha))): + initial_value_contribution += ( + initial_values[k] + / factorial(k) + * (x_points[x_index + 1] ** k - x_points[x_index] ** k) + ) + elif alpha < 1: + raise ValueError("Not yet supported!") + y_prediction[x_index + 1] += initial_value_contribution + y_prediction[x_index + 1] += y_correction[x_index] + y_prediction[x_index + 1] += ( + step_size**alpha + / Gamma(alpha + 1) + * f_name(x_points[x_index], y_correction[x_index]) + ) + subsum = 0 + for j in range(x_index + 1): + subsum += PCcoeffs(alpha, j, x_index) * f_name(x_points[j], y_correction[j]) + y_prediction[x_index + 1] += step_size**alpha / Gamma(alpha + 2) * subsum + + y_correction[x_index + 1] += initial_value_contribution + y_correction[x_index + 1] += y_correction[x_index] + y_correction[x_index + 1] += ( + step_size**alpha + / Gamma(alpha + 2) + * alpha + * f_name(x_points[x_index], y_correction[x_index]) + ) + y_correction[x_index + 1] += ( + step_size**alpha + / Gamma(alpha + 2) + * f_name(x_points[x_index + 1], y_prediction[x_index + 1]) + ) + y_correction[x_index + 1] += step_size**alpha / Gamma(alpha + 2) * subsum + + return y_correction diff --git a/differintP/special.py b/differintP/special.py index 549306c..7d0e191 100644 --- a/differintP/special.py +++ b/differintP/special.py @@ -2,7 +2,8 @@ import numpy as np -from differintP import functionCheck, checkValues, GL +from .utils import functionCheck, checkValues +from .core import GL def GLpoint_via_GL( diff --git a/differintP/utils.py b/differintP/utils.py new file mode 100644 index 0000000..da412bc --- /dev/null +++ b/differintP/utils.py @@ -0,0 +1,68 @@ +from typing import Callable, cast + +import numpy as np + + +def isInteger(n) -> bool: + if n.imag: + return False + if float(n.real).is_integer(): + return True + else: + return False + + +def isPositiveInteger(n) -> bool: + return isInteger(n) and n > 0 + + +def checkValues( + alpha: float, + domain_start: int | float, + domain_end: int | float, + num_points: int, + support_complex_alpha: bool = False, +) -> bool | None: + """Type checking for valid inputs.""" + + assert isPositiveInteger(num_points), ( + "num_points is not an integer: %r" % num_points + ) + + assert isinstance(domain_start, (int, np.integer, float, np.floating)), ( + "domain_start must be integer or float: %r" % domain_start + ) + + assert isinstance(domain_end, (int, np.integer, float, np.floating)), ( + "domain_end must be integer or float: %r" % domain_end + ) + + if not support_complex_alpha: + assert not isinstance( + alpha, complex + ), "Complex alpha not supported for this algorithm." + + return + + +def functionCheck( + f_name: Callable | list | np.ndarray, + domain_start: float | int, + domain_end: float | int, + num_points: int, +): + """Check if function is callable and assign function values.""" + + # Define the function domain and obtain function values. + if hasattr(f_name, "__call__"): + f_name = cast(Callable, f_name) + # If f_name is callable, call it and save to a list. + x = np.linspace(domain_start, domain_end, num_points) + f_values = list(map(lambda t: f_name(t), x)) + step_size = x[1] - x[0] + else: + 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) + return f_values, step_size diff --git a/tests/__init__.py b/tests/__init__.py index e69de29..821d97a 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -0,0 +1 @@ +test_N = 512 \ No newline at end of file diff --git a/tests/test_core.py b/tests/test_core.py index fbd8433..0ecb1b1 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -2,33 +2,29 @@ import numpy as np # Import from sibling directory. -from differintP.core import * # type: ignore +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. -poch_first_argument = 1 -poch_second_argument = 5 -poch_true_answer = 120 size_coefficient_array = 20 -test_N = 512 sqrtpi2 = 0.88622692545275794 truevaluepoly = 0.94031597258 truevaluepoly_caputo = 1.50450555 # 8 / (3 * np.sqrt(np.pi)) truevaluepoly_caputo_higher = 2 / Gamma(1.5) -PC_x_power = np.linspace(0, 1, 100) ** 5.5 - -INTER = GLIinterpolat(1) - -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) -# Get results for checking accuracy. +# 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] @@ -38,189 +34,159 @@ RL_result = RL_r[-1] RL_length = len(RL_r) -# 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 +# --- 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] -class HelperTestCases(unittest.TestCase): - """Tests for helper functions.""" - 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)) +# --- 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) - def test_isPositiveInteger(self): - self.assertTrue(isPositiveInteger(1)) - self.assertFalse(isPositiveInteger(1.1)) - self.assertFalse(isPositiveInteger(-1)) - 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) - 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_GL_binomial_coefficient_array_size(self): - self.assertEqual( - len(GLcoeffs(0.5, size_coefficient_array)) - 1, size_coefficient_array - ) - 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 - ] - - """ Unit tests for gamma function. """ - - def testFiveFactorial(self): - self.assertEqual(Gamma(6), 120) - -# def testNegativePoles(self): -# self.assertEqual(Gamma(-2), np.float64(np.nan)) # The old version was np.inf - - def testRealValue(self): - self.assertEqual(np.round(Gamma(1.25), 12), 0.906402477055) - - def testComplexValue(self): - self.assertEqual(np.round(Gamma(1j), 4), -0.1549 - 0.498j) - - """ 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() - ) +class TestAlgorithmsAccuray(unittest.TestCase): + """Tests for algorithm accuracy.""" - def test_ML_exp(self): - xs = np.arange(10, 0.1) + ####################### + # GLpoint + ####################### + + # sqrt + def test_GLpoint_sqrt_accuracy(self): self.assertTrue( - ( - np.abs(MittagLeffler(1, 1, xs, ignore_special_cases=True) - np.exp(xs)) - <= 1e-3 - ).all() + abs(GLpoint(0.5, lambda x: x**0.5, 0.0, 1.0, 1024) - sqrtpi2) <= 1e-3 ) - def test_ML_geometric(self): - xs = np.arange(1, 0.05) + # x**2 - 1 + def test_GLpoint_accuracy_polynomial(self): self.assertTrue( - ( - np.abs( - MittagLeffler(0, 1, xs, ignore_special_cases=True) - 1 / (1 - xs) - ) - <= 1e-3 - ).all() + 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) -class TestInterpolantCoefficients(unittest.TestCase): - """Test the correctness of the interpolant coefficients.""" - def check_coefficients(self): - self.assertEqual(INTER.prv, -0.125) - self.assertEqual(INTER.crr, 0.75) - self.assertEqual(INTER.nxt, 0.375) + ####################### + # GL + ####################### + # sqrt + def test_GL_accuracy_sqrt(self): + self.assertTrue(abs(GL_result - sqrtpi2) <= 1e-4) -class TestAlgorithms(unittest.TestCase): - """Tests for correct size of algorithm results.""" + # 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 + ) - def test_GLI_result_length(self): - self.assertEqual(GLI_length, test_N) + # 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) - 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)) + ####################### + # GLI + ####################### - """ Tests for algorithm accuracy. """ + # sqrt + def test_GLI_accuracy_sqrt(self): + self.assertTrue(abs(GLI_result - sqrtpi2) <= 1e-4) - def test_GLpoint_sqrt_accuracy(self): + # x**2 - 1 + def test_GLI_accuracy_polynomial(self): self.assertTrue( - abs(GLpoint(0.5, lambda x: x**0.5, 0.0, 1.0, 1024) - sqrtpi2) <= 1e-3 + abs(GLI(0.5, lambda x: x**2 - 1, 0.0, 1.0, 1024)[-1] - truevaluepoly) + <= 6e-3 # low accuracy ) - 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_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 - def test_GL_accuracy_sqrt(self): - self.assertTrue(abs(GL_result - sqrtpi2) <= 1e-4) - def test_GLI_accuracy_sqrt(self): - self.assertTrue(abs(GLI_result - sqrtpi2) <= 1e-4) + ####################### + # 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( @@ -230,6 +196,18 @@ def test_CaputoL1point_accuracy_polynomial(self): <= 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( @@ -239,6 +217,11 @@ def test_CaputoL2point_accuracy_polynomial(self): <= 1e-1 ) + ####################### + # Caputo 2pC + ####################### + + # x**2 (a = 1.5) def test_CaputoL2Cpoint_accuracy_polynomial_higher(self): self.assertTrue( abs( @@ -248,6 +231,7 @@ def test_CaputoL2Cpoint_accuracy_polynomial_higher(self): <= 1e-1 ) + # x**2 - 1 def test_CaputoL2Cpoint_accuracy_polynomial(self): self.assertTrue( abs( @@ -256,36 +240,35 @@ def test_CaputoL2Cpoint_accuracy_polynomial(self): <= 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 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() - ) +class TestAlgorithmsGeneral(unittest.TestCase): + """Tests for correct size of algorithm results.""" - 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_GL_result_length(self): + self.assertEqual(GL_length, test_N) - 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() + 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) diff --git a/tests/test_functions.py b/tests/test_functions.py new file mode 100644 index 0000000..6c60a94 --- /dev/null +++ b/tests/test_functions.py @@ -0,0 +1,73 @@ +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 new file mode 100644 index 0000000..c344577 --- /dev/null +++ b/tests/test_solver.py @@ -0,0 +1,58 @@ +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 new file mode 100644 index 0000000..ab85d04 --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,79 @@ +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