diff --git a/.github/workflows/deploy-docs.yml b/.github/workflows/deploy-docs.yml new file mode 100644 index 0000000..6c3de1a --- /dev/null +++ b/.github/workflows/deploy-docs.yml @@ -0,0 +1,35 @@ +name: Build & Publish Docs with Sphinx + +on: + push: + branches: + - main + +permissions: + contents: write + +jobs: + release: + name: Build + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - name: Set up Python 3.10 + uses: actions/setup-python@v2 + with: + python-version: "3.10" + - name: Install dependencies + run: | + run: python -m pip install .[docs] --upgrade pip + - name: Generate API docs & Build sphinx documentation + run: | + cd docs + make clean + make html + cd .. + - name: Deploy 🚀 + uses: JamesIves/github-pages-deploy-action@v4 + with: + branch: gh-pages # The branch the action should deploy to. + folder: docs/_build/html # The folder the action should deploy. + ssk-key: ${{ secrets.DEPLOY_KEY }} \ No newline at end of file diff --git a/.github/workflows/python-build.yml b/.github/workflows/python-build.yml index 855040b..4950c65 100644 --- a/.github/workflows/python-build.yml +++ b/.github/workflows/python-build.yml @@ -7,12 +7,12 @@ jobs: steps: - uses: actions/checkout@v2 - - name: Set up Python 3.8 + - name: Set up Python 3.10 uses: actions/setup-python@v2 with: - python-version: 3.8 + python-version: "3.10" - name: Install dependencies - run: python -m pip install .[dev] --upgrade pip + run: python -m pip install .[test] --upgrade pip - name: Test with pytest run: pytest - name: Build package distribution @@ -22,7 +22,7 @@ jobs: python -m build --sdist --wheel --outdir dist/ . - name: Publish package distribution to PyPI if: startsWith(github.ref, 'refs/tags') - uses: pypa/gh-action-pypi-publish@master + uses: pypa/gh-action-pypi-publish@release/v1 with: skip_existing: true user: __token__ diff --git a/CODE_OF_CONDUCT.rst b/CODE_OF_CONDUCT.rst new file mode 100644 index 0000000..80e699d --- /dev/null +++ b/CODE_OF_CONDUCT.rst @@ -0,0 +1,113 @@ +Contributor Covenant Code of Conduct +==================================== + +Our Pledge +---------- + +We as members, contributors, and leaders pledge to make participation in our community a harassment-free experience for +everyone, regardless of age, body size, visible or invisible disability, ethnicity, sex characteristics, gender identity +and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, religion, +or sexual identity and orientation. + +We pledge to act and interact in ways that contribute to an open, welcoming, diverse, inclusive, and healthy community. + +Our Standards +------------- + +Examples of behavior that contributes to a positive environment for our community include: + +* Demonstrating empathy and kindness toward other people +* Being respectful of differing opinions, viewpoints, and experiences +* Giving and gracefully accepting constructive feedback +* Accepting responsibility and apologizing to those affected by our mistakes, and learning from the experience +* Focusing on what is best not just for us as individuals, but for the overall community + +Examples of unacceptable behavior include: + +* The use of sexualized language or imagery, and sexual attention or advances of any kind +* Trolling, insulting or derogatory comments, and personal or political attacks +* Public or private harassment +* Publishing others' private information, such as a physical or email address, without their explicit permission +* Other conduct which could reasonably be considered inappropriate in a professional setting + +Enforcement Responsibilities +---------------------------- + +Community leaders are responsible for clarifying and enforcing our standards of acceptable behavior and will take +appropriate and fair corrective action in response to any behavior that they deem inappropriate, threatening, offensive, +or harmful. + +Community leaders have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, +issues, and other contributions that are not aligned to this Code of Conduct, and will communicate reasons for +moderation decisions when appropriate. + +Scope +----- + +This Code of Conduct applies within all community spaces, and also applies when an individual is officially representing +the community in public spaces. Examples of representing our community include using an official e-mail address, posting +via an official social media account, or acting as an appointed representative at an online or offline event. + +Enforcement +----------- + +Instances of abusive, harassing, or otherwise unacceptable behavior may be reported to the community leaders responsible +for enforcement at `zhengp@uw.edu `_. +All complaints will be reviewed and investigated promptly and fairly. + +All community leaders are obligated to respect the privacy and security of the reporter of any incident. + +Enforcement Guidelines +---------------------- + +Community leaders will follow these Community Impact Guidelines in determining the consequences for any action they deem +in violation of this Code of Conduct: + +1. Correction +~~~~~~~~~~~~~ + +**Community Impact**: Use of inappropriate language or other behavior deemed unprofessional or unwelcome in the +community. + +**Consequence**: A private, written warning from community leaders, providing clarity around the nature of the violation +and an explanation of why the behavior was inappropriate. A public apology may be requested. + +2. Warning +~~~~~~~~~~ + +**Community Impact**: A violation through a single incident or series of actions. + +**Consequence**: A warning with consequences for continued behavior. No interaction with the people involved, including +unsolicited interaction with those enforcing the Code of Conduct, for a specified period of time. This includes avoiding +interactions in community spaces as well as external channels like social media. Violating these terms may lead to a +temporary or permanent ban. + +3. Temporary Ban +~~~~~~~~~~~~~~~~ + +**Community Impact**: A serious violation of community standards, including sustained inappropriate behavior. + +**Consequence**: A temporary ban from any sort of interaction or public communication with the community for a specified +period of time. No public or private interaction with the people involved, including unsolicited interaction with those +enforcing the Code of Conduct, is allowed during this period. Violating these terms may lead to a permanent ban. + +4. Permanent Ban +~~~~~~~~~~~~~~~~ + +**Community Impact**: Demonstrating a pattern of violation of community standards, including sustained inappropriate +behavior, harassment of an individual, or aggression toward or disparagement of classes of individuals. + +**Consequence**: A permanent ban from any sort of public interaction within the community. + +Attribution +----------- + +This Code of Conduct is adapted from the `Contributor Covenant `_, +version 2.0, available at ``_. + +Community Impact Guidelines were inspired by +`Mozilla's code of conduct enforcement ladder `_. + +For answers to common questions about this code of conduct, see the FAQ at +``_. Translations are available +at ``_. diff --git a/LICENSE b/LICENSE index 931f1bf..ba43158 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ BSD 2-Clause License -Copyright (c) 2020, Peng Zheng +Copyright (c) 2023, IHME Math Sciences All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/Makefile b/Makefile deleted file mode 100644 index de405f1..0000000 --- a/Makefile +++ /dev/null @@ -1,15 +0,0 @@ -# makefile for xspline package -.PHONY: clean, tests - -clean: - find . -name "*.so*" | xargs rm -rf - find . -name "*.pyc" | xargs rm -rf - find . -name "__pycache__" | xargs rm -rf - find . -name "build" | xargs rm -rf - find . -name "dist" | xargs rm -rf - find . -name "MANIFEST" | xargs rm -rf - find . -name "*.egg-info" | xargs rm -rf - find . -name ".pytest_cache" | xargs rm -rf - -uninstall: - find $(CONDA_PREFIX)/lib/ -name "*xspline*" | xargs rm -rf diff --git a/README.rst b/README.rst index 5471528..a9571ed 100644 --- a/README.rst +++ b/README.rst @@ -1,10 +1,124 @@ -.. image:: https://github.com/zhengp0/xspline/workflows/python-build/badge.svg +.. image:: https://img.shields.io/pypi/l/xspline + :target: https://github.com/zhengp0/xspline/LICENSE + +.. image:: https://img.shields.io/pypi/v/xspline + :target: https://pypi.org/project/xspline + +.. image:: https://img.shields.io/github/actions/workflow/status/zhengp0/xspline/python-build.yml?branch=main :target: https://github.com/zhengp0/xspline/actions -.. image:: https://badge.fury.io/py/xspline.svg - :target: https://badge.fury.io/py/xspline +.. image:: https://img.shields.io/badge/docs-here-green + :target: https://zhengp0.github.io/xspline + XSpline ======= Advanced spline package that provides b-spline bases, their derivatives and integrals. + + +Installation +------------ + +XSpline requires python 3.10 or higher. XSpline only depends on ``numpy>=1.25.1``. +It can be installed via + +.. code:: bash + + pip install xspline>=0.1.0 + +For developers, you can clone the repository and install the package in the +development mode. + +.. code:: + + git clone https://github.com/zhengp0/xspline.git + cd xspline + pip install -e ".[test,docs]" + + +Usage +----- + +You can use XSpline as a univariate function or use it to get design matrix. + +.. code:: python + + import numpy as np + import matplotlib.pyplot as plt + from xspline import XSpline + + spline = XSpline(knots=[0, 0.25, 0.5, 0.75, 1], degree=3) + x = np.arange(0, 1.01, 0.01) + + +One is to use XSpline as a univariate function. In this case, user must provide +coefficients for the spline bases. + +.. code:: python + + np.random.seed(123) + spline.coef = np.random.randn(len(spline)) + y, design_mat = spline(x), spline.get_design_mat(x) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, y) + ax[1].plot(x, design_mat) + +.. image:: docs/images/readme_usage_0.png + +XSpline can be used to obtain derivatives. + +.. code:: python + + dy, ddesign_mat = spline(x, order=1), spline.get_design_mat(x, order=1) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, dy) + ax[1].plot(x, ddesign_mat) + +.. image:: docs/images/readme_usage_1.png + +XSpline can be used to obtain definite integrals. + +.. code:: python + + iy, idesign_mat = spline(x, order=-1), spline.get_design_mat(x, order=-1) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, iy) + ax[1].plot(x, idesign_mat) + +.. image:: docs/images/readme_usage_2.png + +XSpline can extrapolate with different polynomial options + +.. code:: python + + np.random.seed(123) + # constant extrapolation one the left and linear extrapolation on the right + spline = XSpline( + knots=[0, 0.25, 0.5, 0.75, 1], + degree=3, + ldegree=0, + rdegree=1, + coef=np.random.randn(len(spline)), + ) + x = np.arange(-0.5, 1.51, 0.01) + y, design_mat = spline(x), spline.get_design_mat(x) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, y) + ax[1].plot(x, design_mat) + for i in range(len(ax)): + ax[i].vlines( + [0, 1], + ymin=0, + ymax=1, + transform=ax[i].get_xaxis_transform(), + linestyle="--", + linewidth=1, + color="grey", + ) + +.. image:: docs/images/readme_usage_3.png diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..41c270b --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file diff --git a/docs/_static/css/custom.css b/docs/_static/css/custom.css new file mode 100644 index 0000000..1759363 --- /dev/null +++ b/docs/_static/css/custom.css @@ -0,0 +1,20 @@ +/* reduce the size of the main text */ + +p { + font-size: 0.95rem; +} + +h1, +h2, +h3, +h4, +h5, +h6 { + font-weight: 500; +} + +.sidebar-brand-text { + font-size: 1rem; + font-weight: 500; + margin: auto; +} \ No newline at end of file diff --git a/docs/_static/logo/logo-dark.png b/docs/_static/logo/logo-dark.png new file mode 100644 index 0000000..b321abd Binary files /dev/null and b/docs/_static/logo/logo-dark.png differ diff --git a/docs/_static/logo/logo-light.png b/docs/_static/logo/logo-light.png new file mode 100644 index 0000000..1786e8f Binary files /dev/null and b/docs/_static/logo/logo-light.png differ diff --git a/docs/api_reference/index.rst b/docs/api_reference/index.rst new file mode 100644 index 0000000..f1528f3 --- /dev/null +++ b/docs/api_reference/index.rst @@ -0,0 +1,8 @@ +API Reference +============= + +.. toctree:: + :maxdepth: 2 + :glob: + + * \ No newline at end of file diff --git a/docs/api_reference/xspline.bspl.rst b/docs/api_reference/xspline.bspl.rst new file mode 100644 index 0000000..90c4745 --- /dev/null +++ b/docs/api_reference/xspline.bspl.rst @@ -0,0 +1,7 @@ +xspline.bspl +============ + +.. automodule:: xspline.bspl + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/api_reference/xspline.indi.rst b/docs/api_reference/xspline.indi.rst new file mode 100644 index 0000000..356fe9e --- /dev/null +++ b/docs/api_reference/xspline.indi.rst @@ -0,0 +1,7 @@ +xspline.indi +============= + +.. automodule:: xspline.indi + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/api_reference/xspline.poly.rst b/docs/api_reference/xspline.poly.rst new file mode 100644 index 0000000..ba2e359 --- /dev/null +++ b/docs/api_reference/xspline.poly.rst @@ -0,0 +1,7 @@ +xspline.poly +============= + +.. automodule:: xspline.poly + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/api_reference/xspline.xfunction.rst b/docs/api_reference/xspline.xfunction.rst new file mode 100644 index 0000000..347c4ea --- /dev/null +++ b/docs/api_reference/xspline.xfunction.rst @@ -0,0 +1,7 @@ +xspline.xfunction +================= + +.. automodule:: xspline.xfunction + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/api_reference/xspline.xspl.rst b/docs/api_reference/xspline.xspl.rst new file mode 100644 index 0000000..fa9d091 --- /dev/null +++ b/docs/api_reference/xspline.xspl.rst @@ -0,0 +1,7 @@ +xspline.xspl +============= + +.. automodule:: xspline.xspl + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..c1dd8cd --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,91 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import tomllib + +with open("../pyproject.toml", "rb") as f: + meta = tomllib.load(f)["tool"]["sphinx"] + + +# -- Project information ----------------------------------------------------- + +project = meta["project"] +author = meta["author"] +copyright = meta["copyright"] + + +# The full version, including alpha/beta/rc tags +version = meta["version"] + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named "sphinx.ext.*") or your custom +# ones. +extensions = [ + "sphinx.ext.autodoc", + "sphinx.ext.autosectionlabel", + "sphinx.ext.extlinks", + "sphinx.ext.intersphinx", + "sphinx.ext.mathjax", + "sphinx.ext.todo", + "sphinx.ext.viewcode", + "sphinx.ext.napoleon", +] +autodoc_typehints = "description" +autodoc_member_order = "bysource" +autodoc_type_aliases = { + "ArrayLike": "ArrayLike", + "NDArray": "NDArray", + "DataFrame": "DataFrame", +} + +# Add any paths that contain templates here, relative to this directory. +templates_path = ["_templates"] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = "furo" + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ["_static"] +html_css_files = ["css/custom.css"] +html_title = f"{project} {version}" +html_theme_options = { + "sidebar_hide_name": False, + "light_logo": "logo/logo-light.png", + "dark_logo": "logo/logo-dark.png", + "light_css_variables": { + "color-brand-primary": "#008080", + "color-brand-content": "#008080", + "color-problematic": "#BF5844", + "color-background-secondary": "#F8F8F8", + }, + "dark_css_variables": { + "color-brand-primary": "#6FD8D1", + "color-brand-content": "#6FD8D1", + "color-problematic": "#FA9F50", + "color-background-secondary": "#202020", + }, +} diff --git a/docs/images/readme_usage_0.png b/docs/images/readme_usage_0.png new file mode 100644 index 0000000..02bbdf1 Binary files /dev/null and b/docs/images/readme_usage_0.png differ diff --git a/docs/images/readme_usage_1.png b/docs/images/readme_usage_1.png new file mode 100644 index 0000000..922d7a0 Binary files /dev/null and b/docs/images/readme_usage_1.png differ diff --git a/docs/images/readme_usage_2.png b/docs/images/readme_usage_2.png new file mode 100644 index 0000000..d8b7a2c Binary files /dev/null and b/docs/images/readme_usage_2.png differ diff --git a/docs/images/readme_usage_3.png b/docs/images/readme_usage_3.png new file mode 100644 index 0000000..6809ce0 Binary files /dev/null and b/docs/images/readme_usage_3.png differ diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..a44d138 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,26 @@ +XSpline +======= + +Advanced spline package that provides b-spline bases, their derivatives and integrals. + + +To get started please check section :ref:`User's Guide`. +And the full reference is provided in section :ref:`Interface Reference`. + +User's guide +------------ + +.. toctree:: + :maxdepth: 2 + + installation + quickstart + + +Interface reference +------------------- + +.. toctree:: + :maxdepth: 2 + + api_reference/index diff --git a/docs/installation.rst b/docs/installation.rst new file mode 100644 index 0000000..ef72c54 --- /dev/null +++ b/docs/installation.rst @@ -0,0 +1,20 @@ +============ +Installation +============ + + +XSpline requires python 3.10 or higher. XSpline only depends on ``numpy>=1.25.1``. +It can be installed via + +.. code:: bash + + pip install xspline>=0.1.0 + +For developers, you can clone the repository and install the package in the +development mode. + +.. code:: + + git clone https://github.com/zhengp0/xspline.git + cd xspline + pip install -e ".[test,docs]" \ No newline at end of file diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..32bb245 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/quickstart.rst b/docs/quickstart.rst new file mode 100644 index 0000000..f76a730 --- /dev/null +++ b/docs/quickstart.rst @@ -0,0 +1,89 @@ +========== +Quickstart +========== + +Usage +----- + +You can use XSpline as a univariate function or use it to get design matrix. + +.. code:: python + + import numpy as np + import matplotlib.pyplot as plt + from xspline import XSpline + + spline = XSpline(knots=[0, 0.25, 0.5, 0.75, 1], degree=3) + x = np.arange(0, 1.01, 0.01) + + +One is to use XSpline as a univariate function. In this case, user must provide +coefficients for the spline bases. + +.. code:: python + + np.random.seed(123) + spline.coef = np.random.randn(len(spline)) + y, design_mat = spline(x), spline.get_design_mat(x) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, y) + ax[1].plot(x, design_mat) + +.. image:: images/readme_usage_0.png + +XSpline can be used to obtain derivatives. + +.. code:: python + + dy, ddesign_mat = spline(x, order=1), spline.get_design_mat(x, order=1) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, dy) + ax[1].plot(x, ddesign_mat) + +.. image:: images/readme_usage_1.png + +XSpline can be used to obtain definite integrals. + +.. code:: python + + iy, idesign_mat = spline(x, order=-1), spline.get_design_mat(x, order=-1) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, iy) + ax[1].plot(x, idesign_mat) + +.. image:: images/readme_usage_2.png + +XSpline can extrapolate with different polynomial options + +.. code:: python + + np.random.seed(123) + # constant extrapolation one the left and linear extrapolation on the right + spline = XSpline( + knots=[0, 0.25, 0.5, 0.75, 1], + degree=3, + ldegree=0, + rdegree=1, + coef=np.random.randn(len(spline)), + ) + x = np.arange(-0.5, 1.51, 0.01) + y, design_mat = spline(x), spline.get_design_mat(x) + + fig, ax = plt.subplots(1, 2, figsize=(10, 3)) + ax[0].plot(x, y) + ax[1].plot(x, design_mat) + for i in range(len(ax)): + ax[i].vlines( + [0, 1], + ymin=0, + ymax=1, + transform=ax[i].get_xaxis_transform(), + linestyle="--", + linewidth=1, + color="grey", + ) + +.. image:: images/readme_usage_3.png \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..97c4cb7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,28 @@ +[build-system] +requires = ["setuptools"] +build-backend = "setuptools.build_meta" + +[project] +name = "xspline" +version = "0.1.0" +description = "Advanced spline package that provides b-spline bases, their derivatives and integrals" +readme = "REDME.rst" +requires-python = ">=3.10" +license = { file = "LICENSE" } +authors = [ + { name = "IHME Math Sciences", email = "ihme.math.sciences@gmail.com" }, +] +dependencies = ["numpy>=1.25.1"] + +[project.optional-dependencies] +test = ["pytest"] +docs = ["sphinx", "sphinx-autodoc-typehints", "furo"] + +[project.urls] +github = "https://github.com/zhengp0/xspline" + +[tool.sphinx] +project = "xspline" +author = "IHME Math Sciences" +copyright = "2023, IHME Math Sciences" +version = "0.1.0" \ No newline at end of file diff --git a/setup.py b/setup.py deleted file mode 100644 index 78ed98f..0000000 --- a/setup.py +++ /dev/null @@ -1,49 +0,0 @@ -import sys -from pathlib import Path -from setuptools import setup, find_packages - - -if __name__ == '__main__': - base_dir = Path(__file__).parent - src_dir = base_dir/'src'/'xspline' - - sys.path.insert(0, src_dir.as_posix()) - import __about__ as about - - with (base_dir/'README.rst').open() as f: - long_description = f.read() - - install_requirements = [ - 'numpy' - ] - - test_requirements = [ - 'pytest', - 'pytest-mock' - ] - - doc_requirements = [] - - setup(name=about.__title__, - version=about.__version__, - - description=about.__summary__, - long_description=long_description, - license=about.__license__, - url=about.__uri__, - - author=about.__author__, - author_email=about.__email__, - - package_dir={'': 'src'}, - packages=find_packages(where='src'), - include_package_data=True, - - install_requires=install_requirements, - tests_require=test_requirements, - extras_require={ - 'docs': doc_requirements, - 'test': test_requirements, - 'dev': doc_requirements + test_requirements - }, - zip_safe=False,) diff --git a/src/xspline/__about__.py b/src/xspline/__about__.py deleted file mode 100644 index 6b8dac6..0000000 --- a/src/xspline/__about__.py +++ /dev/null @@ -1,16 +0,0 @@ -__all__ = [ - "__title__", "__summary__", "__uri__", "__version__", "__author__", - "__email__", "__license__", "__copyright__", -] - -__title__ = "xspline" -__summary__ = "xspline: Advanced spline tools" -__uri__ = "https://github.com/zhengp0/xspline" - -__version__ = "0.0.7" - -__author__ = "Peng Zheng" -__email__ = "zhengp@uw.edu" - -__license__ = "BSD 2-Clause License" -__copyright__ = f"Copyright 2020 {__author__}" diff --git a/src/xspline/__init__.py b/src/xspline/__init__.py index e110797..8bade48 100644 --- a/src/xspline/__init__.py +++ b/src/xspline/__init__.py @@ -1,9 +1,15 @@ -# -*- coding: utf-8 -*- -""" - xspline - ~~~~~~~ +from .bspl import Bspl +from .indi import Indi +from .poly import Poly +from .xfunction import BasisXFunction, BundleXFunction, XFunction +from .xspl import XSpline - xspline package -""" -from .core import * -from . import utils +__all__ = [ + "Bspl", + "Indi", + "Poly", + "XSpline", + "BasisXFunction", + "BundleXFunction", + "XFunction", +] diff --git a/src/xspline/bspl.py b/src/xspline/bspl.py new file mode 100644 index 0000000..617a971 --- /dev/null +++ b/src/xspline/bspl.py @@ -0,0 +1,231 @@ +import numpy as np +from numpy.typing import NDArray + +from xspline.indi import indi_int +from xspline.typing import BsplParams, RawFunction +from xspline.xfunction import BundleXFunction + + +def cache_bspl(function: RawFunction) -> RawFunction: + """Cache implementation for bspline basis functions, to avoid repetitively + evaluate functions. + + Parameters + ---------- + function + Raw value, derivative and definite integral functions. + + Returns + ------- + describe + Cached version of the raw functions. + + """ + cache = {} + + def wrapper_function(*args, **kwargs) -> NDArray: + key = tuple(tuple(x.ravel()) if isinstance(x, np.ndarray) else x for x in args) + if key in cache: + return cache[key] + result = function(*args, **kwargs) + cache[key] = result + return result + + def cache_clear(): + cache.clear() + + wrapper_function.cache_clear = cache_clear + return wrapper_function + + +@cache_bspl +def bspl_val(params: BsplParams, x: NDArray) -> NDArray: + """Value of the bspline function. + + Parameters + ---------- + params + Bspline function parameters as a tuple including, knots, degree and the + index of the spline basis. + x + Data points. + + Returns + ------- + describe + Function value of the bspline function. + + """ + # knots, degree, and index + t, k, i = params + + if k == 0: + val = np.zeros(x.shape, dtype=x.dtype) + val[(x >= t[i]) & (x < t[i + 1])] = 1.0 + if i == len(t) - 2: + val[x == t[i + 1]] = 1.0 + else: + ii = np.maximum(np.minimum([i, i + 1, i + k, i + k + 1], len(t) - 1), 0) + + val0 = np.zeros(x.shape, dtype=x.dtype) + val1 = np.zeros(x.shape, dtype=x.dtype) + + if t[ii[0]] != t[ii[2]]: + n0 = bspl_val((t, k - 1, i), x) + val0 = (x - t[ii[0]]) * n0 / (t[ii[2]] - t[ii[0]]) + if t[ii[1]] != t[ii[3]]: + n1 = bspl_val((t, k - 1, i + 1), x) + val1 = (t[ii[3]] - x) * n1 / (t[ii[3]] - t[ii[1]]) + + val = val0 + val1 + + return val + + +@cache_bspl +def bspl_der(params: BsplParams, x: NDArray, order: int) -> NDArray: + """Derivative of the bspline function. + + Parameters + ---------- + params + Bspline function parameters as a tuple including, knots, degree and the + index of the spline basis. + x + Data points. + + Returns + ------- + describe + Derivative of the bspline function. + + """ + # knots, degree, and index + t, k, i = params + + if order == 0: + return bspl_val(params, x) + + if order > k: + return np.zeros(x.shape, dtype=x.dtype) + + ii = np.maximum(np.minimum([i, i + 1, i + k, i + k + 1], len(t) - 1), 0) + + val0 = np.zeros(x.shape, dtype=x.dtype) + val1 = np.zeros(x.shape, dtype=x.dtype) + + if t[ii[0]] != t[ii[2]]: + n0 = bspl_der((t, k - 1, i), x, order - 1) + val0 = k * n0 / (t[ii[2]] - t[ii[0]]) + if t[ii[1]] != t[ii[3]]: + n1 = bspl_der((t, k - 1, i + 1), x, order - 1) + val1 = k * n1 / (t[ii[3]] - t[ii[1]]) + + val = val0 - val1 + + return val + + +@cache_bspl +def bspl_int(params: BsplParams, x: NDArray, order: int) -> NDArray: + """Definite integral of the bspline function. + + Parameters + ---------- + params + Bspline function parameters as a tuple including, knots, degree and the + index of the spline basis. + x + Data points. + + Returns + ------- + describe + Definite integral of the bspline function. + + """ + # knots, degree, and index + t, k, i = params + + if order == 0: + return bspl_val(params, x) + + ii = np.maximum(np.minimum([i, i + 1, i + k, i + k + 1], len(t) - 1), 0) + if k == 0: + val = np.zeros(x.shape, dtype=x.dtype) + if t[ii[0]] != t[ii[1]]: + val = indi_int(((t[ii[0]], True), (t[ii[1]], True)), x, order) + else: + val0 = np.zeros(x.shape, dtype=x.dtype) + val1 = np.zeros(x.shape, dtype=x.dtype) + + if t[ii[0]] != t[ii[2]]: + val0 = ( + (x - t[ii[0]]) * bspl_int((t, k - 1, i), x, order) + + order * bspl_int((t, k - 1, i), x, order - 1) + ) / (t[ii[2]] - t[ii[0]]) + if t[ii[1]] != t[ii[3]]: + val1 = ( + (t[ii[3]] - x) * bspl_int((t, k - 1, i + 1), x, order) + - order * bspl_int((t, k - 1, i + 1), x, order - 1) + ) / (t[ii[3]] - t[ii[1]]) + + val = val0 + val1 + + return val + + +def clear_bspl_cache() -> None: + """Clear all cache of the value, derivative and definite integral for + bspline function. + + """ + bspl_val.cache_clear() + bspl_der.cache_clear() + bspl_int.cache_clear() + + +class Bspl(BundleXFunction): + """Basis spline function. + + Parameters + ---------- + params + This is a tuple that contains knots, degree and index of the basis + function. + + Example + ------- + >>> bspl = Bspl(((0.0, 1.0), 1, 0)) # knots=(0.0, 1.0), degree=2, index=0 + >>> bspl([0.0, 1.0]) + array([0., 1.]) + >>> bspl([0.0, 1.0], order=1) + array([1., 1.]) + >>> bspl([0.0, 1.0], order=2) + array([0., 0.]) + >>> bspl([0.0, 1.0], order=-1) + array([0. , 0.5]) + + """ + + def __init__(self, params: BsplParams) -> None: + super().__init__(params, bspl_val, bspl_der, bspl_int) + + +def get_bspl_funs(knots: tuple[float, ...], degree: int) -> tuple[Bspl]: + """Create the bspline basis functions give knots and degree. + + Parameters + ---------- + knots + Bspline knots. + degree + Bspline degree. + + Returns + ------- + describe + A full set of bspline functions. + + """ + return tuple(Bspl((knots, degree, i)) for i in range(-degree, len(knots) - 1)) diff --git a/src/xspline/core.py b/src/xspline/core.py deleted file mode 100644 index 4189aab..0000000 --- a/src/xspline/core.py +++ /dev/null @@ -1,995 +0,0 @@ -# -*- coding: utf-8 -*- -""" - core - ~~~~ - - core module contains main functions and classes. -""" -import numpy as np -from . import utils - - -def bspline_domain(knots, degree, idx, l_extra=False, r_extra=False): - r"""Compute the support for the spline basis, knots degree and the index of - the basis. - - Args: - knots (numpy.ndarray): - 1D array that stores the knots of the splines. - - degree (int): - A non-negative integer that indicates the degree of the polynomial. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - numpy.ndarray: - 1D array with two elements represents that left and right end of the - support of the spline basis. - """ - num_knots = knots.size - num_intervals = num_knots - 1 - num_splines = num_intervals + degree - - if idx == -1: - idx = num_splines - 1 - - lb = knots[max(idx - degree, 0)] - ub = knots[min(idx + 1, num_intervals)] - - if idx == 0 and l_extra: - lb = -np.inf - if idx == num_splines - 1 and r_extra: - ub = np.inf - - return np.array([lb, ub]) - - -def bspline_fun(x, knots, degree, idx, l_extra=False, r_extra=False): - r"""Compute the spline basis. - - Args: - x (float | numpy.ndarray): - Scalar or numpy array that store the independent variables. - - knots (numpy.ndarray): - 1D array that stores the knots of the splines. - - degree (int): - A non-negative integer that indicates the degree of the polynomial. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - float | numpy.ndarray: - Function values of the corresponding spline bases. - """ - num_knots = knots.size - num_intervals = num_knots - 1 - num_splines = num_intervals + degree - - if idx == -1: - idx = num_splines - 1 - - b = bspline_domain(knots, degree, idx, l_extra=l_extra, r_extra=r_extra) - - if degree == 0: - f = utils.indicator_f(x, b, r_close=(idx == num_splines - 1)) - return f - - if idx == 0: - b_effect = bspline_domain(knots, degree, idx) - y = utils.indicator_f(x, b) - z = utils.linear_rf(x, b_effect) - return y*(z**degree) - - if idx == num_splines - 1: - b_effect = bspline_domain(knots, degree, idx) - y = utils.indicator_f(x, b, r_close=True) - z = utils.linear_lf(x, b_effect) - return y*(z**degree) - - lf = bspline_fun(x, knots, degree - 1, idx - 1, - l_extra=l_extra, r_extra=r_extra) - lf *= utils.linear_lf(x, bspline_domain(knots, degree - 1, idx - 1)) - - rf = bspline_fun(x, knots, degree - 1, idx, - l_extra=l_extra, r_extra=r_extra) - rf *= utils.linear_rf(x, bspline_domain(knots, degree - 1, idx)) - - return lf + rf - - -def bspline_dfun(x, knots, degree, order, idx, l_extra=False, r_extra=False): - r"""Compute the derivative of the spline basis. - - Args: - x (float | numpy.ndarray): - Scalar or numpy array that store the independent variables. - - knots (numpy.ndarray): - 1D array that stores the knots of the splines. - - degree (int): - A non-negative integer that indicates the degree of the polynomial. - - order (int): - A non-negative integer that indicates the order of differentiation. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - float | numpy.ndarray: - Derivative values of the corresponding spline bases. - """ - num_knots = knots.size - num_intervals = num_knots - 1 - num_splines = num_intervals + degree - - if idx == -1: - idx = num_splines - 1 - - if order == 0: - return bspline_fun(x, knots, degree, idx, - l_extra=l_extra, r_extra=r_extra) - - if order > degree: - if np.isscalar(x): - return 0.0 - else: - return np.zeros(len(x)) - - if idx == 0: - rdf = 0.0 - else: - b = bspline_domain(knots, degree - 1, idx - 1) - d = b[1] - b[0] - f = (x - b[0])/d - rdf = f*bspline_dfun(x, knots, degree - 1, order, idx - 1, - l_extra=l_extra, r_extra=r_extra) - rdf += order*bspline_dfun(x, knots, degree - 1, order - 1, idx - 1, - l_extra=l_extra, r_extra=r_extra)/d - - if idx == num_splines - 1: - ldf = 0.0 - else: - b = bspline_domain(knots, degree - 1, idx) - d = b[0] - b[1] - f = (x - b[1])/d - ldf = f*bspline_dfun(x, knots, degree - 1, order, idx, - l_extra=l_extra, r_extra=r_extra) - ldf += order*bspline_dfun(x, knots, degree - 1, order - 1, idx, - l_extra=l_extra, r_extra=r_extra)/d - - return ldf + rdf - - -def bspline_ifun(a, x, knots, degree, order, idx, l_extra=False, r_extra=False): - r"""Compute the integral of the spline basis. - - Args: - a (float | numpy.ndarray): - Scalar or numpy array that store the starting point of the integration. - - x (float | numpy.ndarray): - Scalar or numpy array that store the ending point of the integration. - - knots (numpy.ndarray): - 1D array that stores the knots of the splines. - - degree (int): - A non-negative integer that indicates the degree of the polynomial. - - order (int): - A non-negative integer that indicates the order of integration. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - float | numpy.ndarray: - Integral values of the corresponding spline bases. - """ - num_knots = knots.size - num_intervals = num_knots - 1 - num_splines = num_intervals + degree - - if idx == -1: - idx = num_splines - 1 - - if order == 0: - return bspline_fun(x, knots, degree, idx, - l_extra=l_extra, r_extra=r_extra) - - if degree == 0: - b = bspline_domain(knots, degree, idx, l_extra=l_extra, r_extra=r_extra) - return utils.indicator_if(a, x, order, b) - - if idx == 0: - rif = 0.0 - else: - b = bspline_domain(knots, degree - 1, idx - 1) - d = b[1] - b[0] - f = (x - b[0]) / d - rif = f*bspline_ifun(a, x, knots, degree - 1, order, idx - 1, - l_extra=l_extra, r_extra=r_extra) - rif -= order*bspline_ifun(a, x, knots, degree - 1, order + 1, idx - 1, - l_extra=l_extra, r_extra=r_extra)/d - - if idx == num_splines - 1: - lif = 0.0 - else: - b = bspline_domain(knots, degree - 1, idx) - d = b[0] - b[1] - f = (x - b[1]) / d - lif = f*bspline_ifun(a, x, knots, degree - 1, order, idx, - l_extra=l_extra, r_extra=r_extra) - lif -= order*bspline_ifun(a, x, knots, degree - 1, order + 1, idx, - l_extra=l_extra, r_extra=r_extra)/d - - return lif + rif - - -class XSpline: - """XSpline main class of the package. - """ - - def __init__(self, - knots, - degree, - l_linear=False, - r_linear=False, - include_first_basis: bool = True): - r"""Constructor of the XSpline class. - - knots (numpy.ndarray): - 1D numpy array that store the knots, must including that boundary knots. - - degree (int): - A non-negative integer that indicates the degree of the spline. - - l_linear (bool, optional): - A bool variable, that if using the linear tail at left end. - - r_linear (bool, optional): - A bool variable, that if using the linear tail at right end. - """ - # pre-process the knots vector - knots = list(set(knots)) - knots = np.sort(np.array(knots)) - - self.knots = knots - self.degree = degree - self.l_linear = l_linear - self.r_linear = r_linear - self.basis_start = int(not include_first_basis) - - # dimensions - self.num_knots = knots.size - self.num_intervals = knots.size - 1 - - # check inputs - int_l_linear = int(l_linear) - int_r_linear = int(r_linear) - assert self.num_intervals >= 1 + int_l_linear + int_r_linear - assert isinstance(self.degree, int) and self.degree >= 0 - - # create inner knots - self.inner_knots = self.knots[int_l_linear: - self.num_knots - int_r_linear] - self.lb = self.knots[0] - self.ub = self.knots[-1] - self.inner_lb = self.inner_knots[0] - self.inner_ub = self.inner_knots[-1] - - self.num_spline_bases = self.inner_knots.size - 1 + self.degree - self.basis_start - - def domain(self, idx, l_extra=False, r_extra=False): - """Return the support of the XSpline. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - numpy.ndarray: - 1D array with two elements represents that left and right end of the - support of the spline basis. - """ - inner_domain = bspline_domain(self.inner_knots, - self.degree, - idx, - l_extra=l_extra, - r_extra=r_extra) - lb = inner_domain[0] - ub = inner_domain[1] - - lb = self.lb if inner_domain[0] == self.inner_lb else lb - ub = self.ub if inner_domain[1] == self.inner_ub else ub - - return np.array([lb, ub]) - - def fun(self, x, idx, l_extra=False, r_extra=False): - r"""Compute the spline basis. - - Args: - x (float | numpy.ndarray): - Scalar or numpy array that store the independent variables. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - float | numpy.ndarray: - Function values of the corresponding spline bases. - """ - if not self.l_linear and not self.r_linear: - return bspline_fun(x, - self.inner_knots, - self.degree, - idx, - l_extra=l_extra, - r_extra=r_extra) - - x_is_scalar = np.isscalar(x) - if x_is_scalar: - x = np.array([x]) - - f = np.zeros(x.size) - m_idx = np.array([True] * x.size) - - if self.l_linear: - l_idx = (x < self.inner_lb) & ((x >= self.lb) | l_extra) - m_idx &= (x >= self.inner_lb) - - inner_lb_yun = bspline_fun(self.inner_lb, - self.inner_knots, - self.degree, - idx) - inner_lb_dfun = bspline_dfun(self.inner_lb, - self.inner_knots, - self.degree, - 1, idx) - - f[l_idx] = inner_lb_yun + inner_lb_dfun * (x[l_idx] - self.inner_lb) - - if self.r_linear: - u_idx = (x > self.inner_ub) & ((x <= self.ub) | r_extra) - m_idx &= (x <= self.inner_ub) - - inner_ub_yun = bspline_fun(self.inner_ub, - self.inner_knots, - self.degree, - idx) - inner_ub_dfun = bspline_dfun(self.inner_ub, - self.inner_knots, - self.degree, - 1, idx) - - f[u_idx] = inner_ub_yun + inner_ub_dfun * (x[u_idx] - self.inner_ub) - - f[m_idx] = bspline_fun(x[m_idx], - self.inner_knots, - self.degree, - idx, - l_extra=l_extra, - r_extra=r_extra) - - if x_is_scalar: - return f[0] - else: - return f - - def dfun(self, x, order, idx, l_extra=False, r_extra=False): - r"""Compute the derivative of the spline basis. - - Args: - x (float | numpy.ndarray): - Scalar or numpy array that store the independent variables. - - order (int): - A non-negative integer that indicates the order of differentiation. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - float | numpy.ndarray: - Derivative values of the corresponding spline bases. - """ - if order == 0: - return self.fun(x, idx, l_extra=l_extra, r_extra=r_extra) - - if (not self.l_linear) and (not self.r_linear): - return bspline_dfun(x, - self.knots, - self.degree, - order, - idx, - l_extra=l_extra, - r_extra=r_extra) - - x_is_scalar = np.isscalar(x) - if x_is_scalar: - x = np.array([x]) - - dy = np.zeros(x.size) - m_idx = np.array([True] * x.size) - - if self.l_linear: - l_idx = (x < self.inner_lb) & ((x >= self.lb) | l_extra) - m_idx &= (x >= self.inner_lb) - - if order == 1: - inner_lb_dy = bspline_dfun(self.inner_lb, - self.inner_knots, - self.degree, - order, idx) - dy[l_idx] = np.repeat(inner_lb_dy, np.sum(l_idx)) - - if self.r_linear: - u_idx = (x > self.inner_ub) & ((x <= self.ub) | r_extra) - m_idx &= (x <= self.inner_ub) - - if order == 1: - inner_ub_dy = bspline_dfun(self.inner_ub, - self.inner_knots, - self.degree, - order, idx) - dy[u_idx] = np.repeat(inner_ub_dy, np.sum(u_idx)) - - dy[m_idx] = bspline_dfun(x[m_idx], - self.inner_knots, - self.degree, - order, - idx, - l_extra=l_extra, - r_extra=r_extra) - - if x_is_scalar: - return dy[0] - else: - return dy - - def ifun(self, a, x, order, idx, l_extra=False, r_extra=False): - r"""Compute the integral of the spline basis. - - Args: - a (float | numpy.ndarray): - Scalar or numpy array that store the starting point of the - integration. - - x (float | numpy.ndarray): - Scalar or numpy array that store the ending point of the - integration. - - order (int): - A non-negative integer that indicates the order of integration. - - idx (int): - A non-negative integer that indicates the index in the spline bases - list. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - float | numpy.ndarray: - Integral values of the corresponding spline bases. - """ - if order == 0: - return self.fun(x, idx, l_extra=l_extra, r_extra=r_extra) - - if (not self.l_linear) and (not self.r_linear): - return bspline_ifun(a, x, - self.knots, - self.degree, - order, - idx, - l_extra=l_extra, - r_extra=r_extra) - # verify the inputs - assert np.all(a <= x) - - # function and derivative values at inner lb and inner rb - inner_lb_y = bspline_fun(self.inner_lb, - self.inner_knots, - self.degree, - idx) - inner_ub_y = bspline_fun(self.inner_ub, - self.inner_knots, - self.degree, - idx) - inner_lb_dy = bspline_dfun(self.inner_lb, - self.inner_knots, - self.degree, - 1, idx) - inner_ub_dy = bspline_dfun(self.inner_ub, - self.inner_knots, - self.degree, - 1, idx) - - # there are in total 5 pieces functions - def l_piece(a, x, order): - return utils.linear_if(a, x, order, - self.inner_lb, inner_lb_y, inner_lb_dy) - - def m_piece(a, x, order): - return bspline_ifun(a, x, - self.inner_knots, - self.degree, - order, idx, - l_extra=l_extra, r_extra=r_extra) - - def r_piece(a, x, order): - return utils.linear_if(a, x, order, - self.inner_ub, inner_ub_y, inner_ub_dy) - - def zero_piece(a, x, order): - if np.isscalar(a) and np.isscalar(x): - return 0.0 - elif np.isscalar(a): - return np.zeros(x.size) - else: - return np.zeros(a.size) - - funcs = [] - knots = [] - if not l_extra: - funcs.append(zero_piece) - if self.l_linear: - funcs.append(l_piece) - funcs.append(m_piece) - if self.r_linear: - funcs.append(r_piece) - if not r_extra: - funcs.append(zero_piece) - - if not l_extra: - knots.append(self.lb) - knots.append(self.inner_lb) - if self.l_linear: - knots.append(self.inner_lb) - if self.r_linear: - knots.append(self.inner_ub) - if not r_extra: - knots.append(self.inner_ub) - knots.append(self.ub) - - knots = np.sort(list(set(knots))) - - return utils.pieces_if(a, x, order, funcs, knots) - - def design_mat(self, x, l_extra=False, r_extra=False): - r"""Compute the design matrix of spline basis. - - Args: - x (float | numpy.ndarray): - Scalar or numpy array that store the independent variables. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - numpy.ndarray: - Return design matrix. - """ - mat = np.vstack([ - self.fun(x, idx, l_extra=l_extra, r_extra=r_extra) - for idx in range(self.basis_start, self.num_spline_bases + self.basis_start) - ]).T - return mat - - def design_dmat(self, x, order, l_extra=False, r_extra=False): - r"""Compute the design matrix of spline basis derivatives. - - Args: - x (float | numpy.ndarray): - Scalar or numpy array that store the independent variables. - - order (int): - A non-negative integer that indicates the order of differentiation. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - numpy.ndarray: - Return design matrix. - """ - dmat = np.vstack([ - self.dfun(x, order, idx, l_extra=l_extra, r_extra=r_extra) - for idx in range(self.basis_start, self.num_spline_bases + self.basis_start) - ]).T - return dmat - - def design_imat(self, a, x, order, l_extra=False, r_extra=False): - r"""Compute the design matrix of the integrals of the spline bases. - - Args: - a (float | numpy.ndarray): - Scalar or numpy array that store the starting point of the - integration. - - x (float | numpy.ndarray): - Scalar or numpy array that store the ending point of the - integration. - - order (int): - A non-negative integer that indicates the order of integration. - - l_extra (bool, optional): - A optional bool variable indicates that if extrapolate at left end. - Default to be False. - - r_extra (bool, optional): - A optional bool variable indicates that if extrapolate at right end. - Default to be False. - - Returns: - numpy.ndarray: - Return design matrix. - """ - imat = np.vstack([ - self.ifun(a, x, order, idx, l_extra=l_extra, r_extra=r_extra) - for idx in range(self.basis_start, self.num_spline_bases + self.basis_start) - ]).T - return imat - - def last_dmat(self): - """Compute highest order of derivative in domain. - - Returns: - numpy.ndarray: - 1D array that contains highest order of derivative for intervals. - """ - # compute the last dmat for the inner domain - dmat = self.design_dmat(self.inner_knots[:-1], self.degree) - - if self.l_linear: - dmat = np.vstack((self.design_dmat(np.array([self.inner_lb]), 1), - dmat)) - - if self.r_linear: - dmat = np.vstack((dmat, - self.design_dmat(np.array([self.inner_ub]), 1))) - - return dmat - - -class NDXSpline: - """Multi-dimensional xspline. - """ - - def __init__(self, ndim, knots_list, degree_list, - l_linear_list=None, - r_linear_list=None, - include_first_basis_list=True): - """Constructor of ndXSpline class - - Args: - ndim (int): - Number of dimension. - knots_list (list{numpy.ndarray}): - List of knots for every dimension. - degree_list (list{int}): - List of degree for every dimension. - l_linear_list (list{bool} | None, optional): - List of indicator of if have left linear tail for each - dimension. - r_linear_list (list{bool} | None, optional): - List of indicator of if have right linear tail for each - dimension. - """ - self.ndim = ndim - self.knots_list = knots_list - self.degree_list = degree_list - self.l_linear_list = utils.option_to_list(l_linear_list, self.ndim) - self.r_linear_list = utils.option_to_list(r_linear_list, self.ndim) - self.include_first_basis_list = utils.option_to_list(include_first_basis_list, self.ndim) - - self.spline_list = [ - XSpline(self.knots_list[i], self.degree_list[i], - l_linear=self.l_linear_list[i], - r_linear=self.r_linear_list[i], - include_first_basis=self.include_first_basis_list[i]) - for i in range(self.ndim) - ] - - self.num_knots_list = np.array([ - spline.num_knots for spline in self.spline_list]) - self.num_intervals_list = np.array([ - spline.num_intervals for spline in self.spline_list]) - self.num_spline_bases_list = np.array([ - spline.num_spline_bases for spline in self.spline_list]) - - self.num_knots = self.num_knots_list.prod() - self.num_intervals = self.num_intervals_list.prod() - self.num_spline_bases = self.num_spline_bases_list.prod() - - def design_mat(self, x_list, - is_grid=True, - l_extra_list=None, - r_extra_list=None): - """Design matrix of the spline basis. - - Args: - x_list (list{numpy.ndarray}): - A list of coordinates for each dimension, they should have the - same dimension or come in matrix form. - is_grid (bool, optional): - If `True` treat the coordinates from `x_list` as the grid points - and compute the mesh grid from it, otherwise, treat each group - of the coordinates independent. - l_extra_list (list{bool} | None, optional): - Indicators of if extrapolate in the left side for each - dimension. - r_extra_list (list{bool} | None, optional): - Indicators of if extrapolate in the right side for each - dimension. - - Returns: - numpy.ndarray: - Design matrix. - """ - l_extra_list = utils.option_to_list(l_extra_list, self.ndim) - r_extra_list = utils.option_to_list(r_extra_list, self.ndim) - - assert len(x_list) == self.ndim - assert len(l_extra_list) == self.ndim - assert len(r_extra_list) == self.ndim - - mat_list = [spline.design_mat(x_list[i], - l_extra=l_extra_list[i], - r_extra=r_extra_list[i]) - for i, spline in enumerate(self.spline_list)] - - if is_grid: - mat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [mat_list[j][:, index_list[j]] - for j in range(self.ndim)] - mat.append(utils.outer_flatten(*bases_list)) - else: - num_points = x_list[0].size - assert np.all([x_list[i].size == num_points - for i in range(self.ndim)]) - mat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [mat_list[j][:, index_list[j]] - for j in range(self.ndim)] - mat.append(np.prod(bases_list, axis=0)) - - return np.ascontiguousarray(np.vstack(mat).T) - - def design_dmat(self, x_list, n_list, - is_grid=True, - l_extra_list=None, - r_extra_list=None): - """Design matrix of the derivatives of spline basis. - - Args: - x_list (list{numpy.ndarray}): - A list of coordinates for each dimension, they should have the - same dimension or come in matrix form. - n_list (list{int}): - A list of integers indicates the order of differentiation for - each dimension. - is_grid (bool, optional): - If `True` treat the coordinates from `x_list` as the grid points - and compute the mesh grid from it, otherwise, treat each group - of the coordinates independent. - l_extra_list (list{bool} | None, optional): - Indicators of if extrapolate in the left side for each - dimension. - r_extra_list (list{bool} | None, optional): - Indicators of if extrapolate in the right side for each - dimension. - - Returns: - numpy.ndarray: - Differentiation design matrix. - """ - l_extra_list = utils.option_to_list(l_extra_list, self.ndim) - r_extra_list = utils.option_to_list(r_extra_list, self.ndim) - - assert len(x_list) == self.ndim - assert len(n_list) == self.ndim - assert len(l_extra_list) == self.ndim - assert len(r_extra_list) == self.ndim - - dmat_list = [spline.design_dmat(x_list[i], n_list[i], - l_extra=l_extra_list[i], - r_extra=r_extra_list[i]) - for i, spline in enumerate(self.spline_list)] - - if is_grid: - dmat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [dmat_list[j][:, index_list[j]] - for j in range(self.ndim)] - dmat.append(utils.outer_flatten(*bases_list)) - else: - num_points = x_list[0].size - assert np.all([x_list[i].size == num_points - for i in range(self.ndim)]) - dmat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [dmat_list[j][:, index_list[j]] - for j in range(self.ndim)] - dmat.append(np.prod(bases_list, axis=0)) - - return np.ascontiguousarray(np.vstack(dmat).T) - - def design_imat(self, a_list, x_list, n_list, - is_grid=True, - l_extra_list=None, - r_extra_list=None): - """Design matrix of the spline basis. - - Args: - a_list (list{numpy.ndarray}): - Start of integration of coordinates for each dimension. - x_list (list{numpy.ndarray}): - A list of coordinates for each dimension, they should have the - same dimension or come in matrix form. - n_list (list{int}): - A list of integers indicates the order of integration for - each dimension. - is_grid (bool, optional): - If `True` treat the coordinates from `x_list` as the grid points - and compute the mesh grid from it, otherwise, treat each group - of the coordinates independent. - l_extra_list (list{bool} | None, optional): - Indicators of if extrapolate in the left side for each - dimension. - r_extra_list (list{bool} | None, optional): - Indicators of if extrapolate in the right side for each - dimension. - - Returns: - numpy.ndarray: - Integration design matrix. - """ - l_extra_list = utils.option_to_list(l_extra_list, self.ndim) - r_extra_list = utils.option_to_list(r_extra_list, self.ndim) - - assert len(a_list) == self.ndim - assert len(x_list) == self.ndim - assert len(n_list) == self.ndim - assert len(l_extra_list) == self.ndim - assert len(r_extra_list) == self.ndim - - imat_list = [spline.design_imat(a_list[i], x_list[i], n_list[i], - l_extra=l_extra_list[i], - r_extra=r_extra_list[i]) - for i, spline in enumerate(self.spline_list)] - - if is_grid: - imat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [imat_list[j][:, index_list[j]] - for j in range(self.ndim)] - imat.append(utils.outer_flatten(*bases_list)) - else: - num_points = x_list[0].size - assert np.all([x_list[i].size == num_points - for i in range(self.ndim)]) - imat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [imat_list[j][:, index_list[j]] - for j in range(self.ndim)] - imat.append(np.prod(bases_list, axis=0)) - - return np.ascontiguousarray(np.vstack(imat).T) - - def last_dmat(self): - """Highest order of derivative matrix. - - Returns: - numpy.ndarray: - Design matrix contain the highest order of derivative. - """ - mat_list = [spline.last_dmat() for spline in self.spline_list] - - mat = [] - for i in range(self.num_spline_bases): - index_list = utils.order_to_index(i, self.num_spline_bases_list) - bases_list = [mat_list[j][:, index_list[j]] - for j in range(self.ndim)] - mat.append(utils.outer_flatten(*bases_list)) - - return np.ascontiguousarray(np.vstack(mat).T) - - -# TODO: -# 1. bspline function pass in too many default every time -# 2. name of f, df and if -# 3. the way to deal with the scalar vs array. -# 4. keep the naming scheme consistent. diff --git a/src/xspline/indi.py b/src/xspline/indi.py new file mode 100644 index 0000000..c391b7f --- /dev/null +++ b/src/xspline/indi.py @@ -0,0 +1,116 @@ +from math import factorial + +import numpy as np + +from xspline.typing import IndiParams, NDArray +from xspline.xfunction import BundleXFunction + + +def indi_val(params: IndiParams, x: NDArray) -> NDArray: + """Indicator value function, + + Parameters + ---------- + params + Indicator parameters as a tuple consists of lower and upper bound of + the interval corresponding to the indicator function. + x + Data points. + + Returns + ------- + describe + Indicator function value. + + """ + # lower and upper bounds + lb, ub = params + + val = np.zeros(x.size, dtype=x.dtype) + ind = (x > lb[0]) & (x < ub[0]) + if lb[1]: + ind = ind | (x == lb[0]) + if ub[1]: + ind = ind | (x == ub[0]) + val[ind] = 1.0 + return val + + +def indi_der(params: IndiParams, x: NDArray, order: int) -> NDArray: + """Indicator derivative function. Since indicator function is a piecewise + constant function, its derivative will always be zero. + + Parameters + ---------- + params + Indicator parameters as a tuple consists of lower and upper bound of + the interval corresponding to the indicator function. + x + Data points. + + Returns + ------- + describe + Indicator deviative value. + + """ + return np.zeros(x.size, dtype=x.dtype) + + +def indi_int(params: IndiParams, x: NDArray, order: int) -> NDArray: + """Indicator definite integral function. It is a piecewise polynomial + function. + + Parameters + ---------- + params + Indicator parameters as a tuple consists of lower and upper bound of + the interval corresponding to the indicator function. + x + Data points. + + Returns + ------- + describe + Indicator definite integral value. + + """ + # lower and upper bounds + lb, ub = params + + val = np.zeros(x.size, dtype=x.dtype) + ind0 = (x >= lb[0]) & (x <= ub[0]) + ind1 = x > ub[0] + val[ind0] = (x[ind0] - lb[0]) ** (-order) / factorial(-order) + for i in range(-order): + val[ind1] += ((ub[0] - lb[0]) ** (-order - i) / factorial(-order - i)) * ( + (x[ind1] - ub[0]) ** i / factorial(i) + ) + return val + + +class Indi(BundleXFunction): + """Indicator function. + + Parameters + ---------- + params + This is a tuple contains the lower and upper bounds of the indicator + function. For each bound it consists of a number for the location of the + bound and a boolean for the inclusion of the bound. For example, if we + pass in `((0.0, True), (1.0, False))`, this represents interval [0, 1). + + Example + ------- + >>> indi = Indi(((0.0, True), (1.0, False))) + >>> indi([-1.0, 0.0, 1.0]) + array([0., 1., 0.]) + >>> indi([-1.0, 0.0, 1.0], order=1) + array([0., 0., 0.]) + >>> indi([-1.0, 0.0, 1.0], order=-1) + array([0., 0., 1.]) + + """ + + def __init__(self, params: IndiParams) -> None: + super().__init__(params, indi_val, indi_der, indi_int) diff --git a/src/xspline/poly.py b/src/xspline/poly.py new file mode 100644 index 0000000..ec250d4 --- /dev/null +++ b/src/xspline/poly.py @@ -0,0 +1,144 @@ +from math import factorial + +import numpy as np + +from xspline.typing import PolyParams, NDArray +from xspline.xfunction import BundleXFunction, XFunction + + +def poly_val(params: PolyParams, x: NDArray) -> NDArray: + """Polynomial value function. + + Parameters + ---------- + params + Polynomial coefficients. + x + Data points. + + Returns + ------- + describe + Polynomial function values. + + """ + return np.polyval(params, x) + + +def poly_der(params: PolyParams, x: NDArray, order: int) -> NDArray: + """Polynomial derivative function. + + Parameters + ---------- + params + Polynomial coefficients. + x + Data points. + order + Order of differentiation. + + Returns + ------- + describe + Polynomial derivative values. + + """ + return np.polyval(np.polyder(params, order), x) + + +def poly_int(params: PolyParams, x: NDArray, order: int) -> NDArray: + """Polynomial definite integral function. + + Parameters + ---------- + params + Polynomial coefficients. + x + Data points. + order + Order of integration. Here we use negative integer. + + Returns + ------- + describe + Polynomial definite integral values. + + """ + + return np.polyval(np.polyint(params, -order), x) + + +class Poly(BundleXFunction): + """Polynomial function. A simple wrapper for the numpy poly functions. + + Parameters + ---------- + params + This a tuple contains coefficients for the terms in polynomial. + + Example + ------- + >>> poly = Poly((1.0, 0.0)) + >>> poly([0.0, 1.0]) + array([0.0, 1.0]) + >>> poly([0.0, 1.0], order=1) + array([1.0, 1.0]) + >>> poly([0.0, 1.0], order=2) + array([0.0, 0.0]) + >>> poly([0.0, 1.0]], order=-1) + array([0.0, 0.5]) + + """ + + def __init__(self, params: PolyParams) -> None: + super().__init__(params, poly_val, poly_der, poly_int) + + +def get_poly_params(fun: XFunction, x: float, degree: int) -> tuple[float, ...]: + """Solve polynomial (taylor) coefficients provided the ``XFunction``. + + Parameters + ---------- + fun + Provided ``XFunction`` to be approximated. + x + The point where we want to approximate ``XFunction`` by the polynomial. + degree + Degree of the approximation polynomial. + + Returns + ------- + describe + The approximation polynomial coefficients. + + """ + if degree < 0: + return (0.0,) + rhs = np.array([fun(x, order=i) for i in range(degree, -1, -1)]) + mat = np.zeros((degree + 1, degree + 1)) + for i in range(degree + 1): + for j in range(i + 1): + mat[i, j] = factorial(degree - j) / factorial(i - j) * x ** (i - j) + return tuple(np.linalg.solve(mat, rhs)) + + +def get_poly_fun(fun: XFunction, x: float, degree: int) -> Poly: + """Get the approximation polynomial function. + + Parameters + ---------- + fun + Provided ``XFunction`` to be approximated. + x + The point where we want to approximate ``XFunction`` by the polynomial. + degree + Degree of the approximation polynomial. + + Returns + ------- + describe + Instance of the ``Poly`` class to approximate provided ``XFunction``. + + """ + params = get_poly_params(fun, x, degree) + return Poly(params) diff --git a/src/xspline/typing.py b/src/xspline/typing.py new file mode 100644 index 0000000..7ba201f --- /dev/null +++ b/src/xspline/typing.py @@ -0,0 +1,17 @@ +from typing import Callable + +from numpy.typing import NDArray + +VFunction = Callable[[NDArray], NDArray] +DFunction = Callable[[NDArray, int], NDArray] +IFunction = Callable[[NDArray, int], NDArray] + +RawVFunction = Callable[[tuple, NDArray], NDArray] +RawDFunction = Callable[[tuple, NDArray, int], NDArray] +RawIFunction = Callable[[tuple, NDArray, int], NDArray] +RawFunction = RawVFunction | RawDFunction | RawIFunction + +BoundaryPoint = tuple[float, bool] +IndiParams = tuple[BoundaryPoint, BoundaryPoint] +PolyParams = tuple[float, ...] +BsplParams = tuple[tuple[float, ...], int, int] diff --git a/src/xspline/utils.py b/src/xspline/utils.py deleted file mode 100644 index 2253912..0000000 --- a/src/xspline/utils.py +++ /dev/null @@ -1,449 +0,0 @@ -# utility functions for the bsplinex class -import numpy as np - - -def indicator_f(x, b, l_close=True, r_close=False): - r"""Indicator function for provided interval. - - Args: - x (float | numpy.ndarray): - Float scalar or numpy array store the variable(s). - - b (numpy.ndarray): - 1D array with 2 elements represent the left and right end of the - interval. - - l_close (bool | True, optional): - Bool variable indicate that if include the left end of the interval. - - r_close (bool | False, optional): - Bool variable indicate that if include the right end of the interval. - - Returns: - float | numpy.ndarray: - Return function value(s) at ``x``. The result has the same shape with - ``x``. 0 in the result indicate corresponding element in ``x`` is not in - the interval and 1 means the it is in the interval. - """ - if l_close: - lb = (x >= b[0]) - else: - lb = (x > b[0]) - - if r_close: - rb = (x <= b[1]) - else: - rb = (x < b[1]) - - if np.isscalar(x): - return float(lb & rb) - else: - return (lb & rb).astype(np.double) - - -def linear_f(x, z, fz, dfz): - r"""Linear function construct by reference point(s). - - Args: - x (float | numpy.ndarray): - Float scalar or numpy array store the variable(s). - - z (float | numpy.ndarray): - Float scalar or numpy array store the reference point(s). When ``x`` is - ``np.ndarray``, ``z`` has to have the same shape with ``x``. - - fz (float | numpy.ndarray): - Float scalar or numpy array store the function value(s) at ``z``. Same - requirements with ``z``. - - dfz (float | numpy.ndarray): - Float scalar or numpy array store the function derivative(s) at ``z``. - Same requirements with ``z``. - - Returns: - float | numpy.ndarray: - Return function value(s) at ``x`` with linear function constructed at - base point(s) ``z``. The result has the same shape with ``x`` when - ``z`` is scalar or same shape with ``z`` when ``x`` is scalar. - """ - return fz + dfz*(x - z) - - -def linear_lf(x, b): - r"""Linear function constructed by linearly interpolate 0 and 1 from left - end point to right end point. - - Args: - x (float | numpy.ndarray): - Float scalar or numpy array that store the variable(s). - - b (numpy.ndarray): - 1D array with 2 elements represent the left and right end of the - interval. - - Returns: - float | numpy.ndarray: - Return function value(s) at ``x``. The result has the same shape with - ``x``. - """ - return (x - b[0])/(b[1] - b[0]) - - -def linear_rf(x, b): - r"""Linear function constructed by linearly interpolate 0 and 1 from right - end point to left end point. - - Args: - x (float | numpy.ndarray): - Float scalar or numpy array that store the variable(s). - - b (numpy.ndarray): - 1D array with 2 elements represent the left and right end of the - interval. - - Returns: - float | numpy.ndarray: - Return function value(s) at ``x``. The result has the same shape with - ``x``. - """ - return (x - b[1])/(b[0] - b[1]) - - -def constant_if(a, x, order, c): - r"""Integration of constant function. - - Args: - a (float | numpy.ndarray): - Starting point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - x (float | numpy.ndarray): - Ending point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - order (int): - Non-negative integer number indicate the order of integration. In other - words, how many time(s) we integrate. - - c (float): - Constant function value. - - Returns: - float | numpy.ndarray: - Integration value(s) of the constant function. - """ - # determine the result size - a_is_ndarray = isinstance(a, np.ndarray) - x_is_ndarray = isinstance(x, np.ndarray) - - if a_is_ndarray and x_is_ndarray: - assert a.size == x.size - - result_is_ndarray = a_is_ndarray or x_is_ndarray - if a_is_ndarray: - result_size = a.size - elif x_is_ndarray: - result_size = x.size - else: - result_size = 1 - - # special case when c is 0 - if c == 0.0: - if result_is_ndarray: - return np.zeros(result_size) - else: - return 0.0 - - return c*(x - a)**order/np.math.factorial(order) - - -def linear_if(a, x, order, z, fz, dfz): - r"""Integrate linear function constructed by the reference point(s). - - Args: - a (float | numpy.ndarray): - Starting point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - x (float | numpy.ndarray): - Ending point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - order (int): - Non-negative integer number indicate the order of integration. In other - words, how many time(s) we integrate. - - z (float | numpy.ndarray): - Float scalar or numpy array store the reference point(s). When ``x`` is - ``np.ndarray``, ``z`` has to have the same shape with ``x``. - - fz (float | numpy.ndarray): - Float scalar or numpy array store the function value(s) at ``z``. Same - requirements with ``z``. - - dfz (float | numpy.ndarray): - Float scalar or numpy array store the function derivative(s) at ``z``. - Same requirements with ``z``. - - Returns: - float | numpy.ndarray: - Integration value(s) of the constant function. - """ - fa = fz + dfz*(a - z) - dfa = dfz - - return dfa*(x - a)**(order + 1)/np.math.factorial(order + 1) + \ - fa*(x - a)**order/np.math.factorial(order) - - -def integrate_across_pieces(a, x, order, funcs, knots): - r"""Integrate Across piecewise functions. - - Args: - a (float | numpy.ndarray): - Starting point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. ``a`` has to be less - than ``knots[0]``. - - x (float | numpy.ndarray): - Ending point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. ``x`` has to be greater - than ``knots[-1]``. - - order (int): - Non-negative integer number indicate the order of integration. In other - words, how many time(s) we integrate. - - funcs (list): - List of functions defined on pieces of domain. - - knots (numpy.ndarray): - 1D numpy array contains all the breaking points for the piecewise - functions. - - Returns: - float | numpy.ndarray: - Integration value(s) of the constant function. - """ - assert len(funcs) == len(knots) + 1 - if len(funcs) == 1: - return funcs[0](a, x, order) - else: - assert np.all(a < knots[0]) and np.all(x > knots[-1]), f"{a} should be in [{knots[0]}, {knots[-1]}]." - - if np.isscalar(a): - b = knots[0] - else: - b = np.repeat(knots[0], a.size) - - val = integrate_across_pieces(b, x, order, funcs[1:], knots[1:]) - - for j in range(order): - val += funcs[0](a, b, order - j)*(x - b)**j / np.math.factorial(j) - - return val - - -def pieces_if(a, x, order, funcs, knots): - r"""Integrate pieces of the functions. - - Args: - a (float | numpy.ndarray): - Starting point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - x (float | numpy.ndarray): - Ending point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - order (int): - Non-negative integer number indicate the order of integration. In other - words, how many time(s) we integrate. - - funcs (list): - List of functions defined on pieces of domain. - - knots (numpy.ndarray): - 1D numpy array contains all the breaking points for the piecewise - functions. - - Returns: - float | numpy.ndarray: - Integration value(s) of the constant function. - """ - # verify the input - assert len(funcs) == len(knots) + 1 - if np.isscalar(a) and not np.isscalar(x): - a = np.repeat(a, x.size) - if np.isscalar(x) and not np.isscalar(a): - x = np.repeat(x, a.size) - if np.isscalar(a) and np.isscalar(x): - a = np.array([a]) - x = np.array([x]) - result_is_scalar = True - else: - result_is_scalar = False - - assert a.size == x.size - assert np.all(a <= x) - assert len(funcs) == len(knots) + 1 - - num_knots = len(knots) - - # different cases - a_ind = [a < knots[0]] +\ - [(a >= knots[i]) & (a < knots[i + 1]) - for i in range(num_knots - 1)] +\ - [a >= knots[-1]] - - x_ind = [x <= knots[0]] +\ - [(x > knots[i]) & (x <= knots[i + 1]) - for i in range(num_knots - 1)] +\ - [x > knots[-1]] - - int_f = np.zeros(a.size) - for ia in range(len(funcs)): - for ix in range(ia, len(funcs)): - case_id = a_ind[ia] & x_ind[ix] - if np.any(case_id): - int_f[case_id] = integrate_across_pieces(a[case_id], - x[case_id], - order, - funcs[ia:ix + 1], - knots[ia:ix]) - - if result_is_scalar: - return int_f[0] - else: - return int_f - - -def indicator_if(a, x, order, b, l_close=True, r_close=False): - r"""Integrate indicator function. - - Args: - a (float | numpy.ndarray): - Starting point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - x (float | numpy.ndarray): - Ending point(s) of the integration. If both ``a`` and ``x`` are - ``np.ndarray``, they should have the same shape. - - order (int): - Non-negative integer number indicate the order of integration. In other - words, how many time(s) we integrate. - - b (numpy.ndarray): - 1D array with 2 elements represent the left and right end of the - interval. - - l_close (bool | True, optional): - Bool variable indicate that if include the left end of the interval. - - r_close (bool | False, optional): - Bool variable indicate that if include the right end of the interval. - - Returns: - float | numpy.ndarray: - Integration value(s) of the constant function. - """ - if order == 0: - return indicator_f(x, b, l_close=l_close, r_close=r_close) - else: - return pieces_if(a, x, order, - [lambda *params: constant_if(*params, 0.0), - lambda *params: constant_if(*params, 1.0), - lambda *params: constant_if(*params, 0.0)], b) - - -def seq_diff_mat(size): - r"""Compute sequencial difference matrix. - - Args: - size (int): - Positive integer that indicate the number of columns of the matrix. - Number of rows will be 1 less than number of columns. - - Returns: - ndarray: - A matrix consist of one and minus one. - """ - assert isinstance(size, int) and size >= 2 - - mat = np.zeros((size - 1, size)) - id_d0 = np.diag_indices(size - 1) - id_d1 = (id_d0[0], id_d0[1] + 1) - - mat[id_d0] = -1.0 - mat[id_d1] = 1.0 - - return mat - - -def order_to_index(order, shape): - r"""Compute the index of the element in a high dimensional array provided - the order of the element. - - Args: - order: int - Non-negative integer present the order of the element. - - shape: tuple - Shape tuple of the high dimensional array. - - Returns: - tuple: - The index element in the array. - """ - assert hasattr(shape, '__iter__') - assert isinstance(order, int) - assert 0 <= order < np.prod(shape) - - ndim = len(shape) - n = np.cumprod(np.insert(shape[::-1], 0, 1))[-2::-1] - - index = [] - for j in range(ndim): - quotient = order // n[j] - index.append(quotient) - order -= quotient*n[j] - - return tuple(index) - - -def option_to_list(opt, size): - r"""Convert default option to list of options. - - Args: - opt (bool | None): - A single option in the form of bool or None. - - size (int): - Positive integer indicate the size of the option list. - - Returns: - list: - Option list consist of bool elements. - """ - assert isinstance(size, int) - assert size > 0 - if not opt: - return [False]*size - else: - return [True]*size - - -def outer_flatten(*args): - r"""Outer product of multiple vectors and then flatten the result. - - Args: - args (list | tuple): - A list or tuple of 1D numpy arrays. - - Return: - numpy.ndarray: - 1D numpy array that store the flattened outer product. - """ - result = np.prod(np.ix_(*args)) - return result.reshape(result.size,) diff --git a/src/xspline/xfunction.py b/src/xspline/xfunction.py new file mode 100644 index 0000000..3818b84 --- /dev/null +++ b/src/xspline/xfunction.py @@ -0,0 +1,258 @@ +from __future__ import annotations +from functools import partial +from math import factorial +from operator import attrgetter + +import numpy as np + +from xspline.typing import ( + BoundaryPoint, + Callable, + RawDFunction, + RawIFunction, + RawVFunction, + NDArray, +) + + +class XFunction: + """Function interface that provide easy access to function value, + derivatives and definite integrals. + + Parameters + ---------- + fun + Function implementation. + + TODO: describe the interface of the function implementation + + """ + + def __init__(self, fun: Callable) -> None: + self._fun = fun + + def _check_args(self, x: NDArray, order: int) -> tuple[NDArray, int, bool]: + x, order = np.asarray(x, dtype=float), int(order) + if (x.ndim not in [0, 1, 2]) or (x.ndim == 2 and len(x) != 2): + raise ValueError( + "please provide a scalar, an 1d array, or a 2d " "array with two rows" + ) + + # reshape array + isscalar = x.ndim == 0 + if isscalar: + x = x.ravel() + if order >= 0 and x.ndim == 2: + raise ValueError( + "please provide an 1d array for function value " + "defivative computation" + ) + if order < 0 and x.ndim == 1: + x = np.vstack([np.repeat(x.min(), x.size), x]) + + # check interval bounds + if order < 0 and (x[0] > x[1]).any(): + raise ValueError("to integrate, `x` must satisfy `x[0] <= x[1]`") + + return x, order, isscalar + + def __call__(self, x: NDArray, order: int = 0) -> NDArray: + """Function returns function values, derivatives and definite integrals. + + Parameters + ---------- + x + Data points where the function is evaluated. If `order < 0` and `x` + is a 2d array with two rows, the rows will be treated as the + starting and ending points for definite interval. If `order < 0` and + `x` is a 1d array, function will use the smallest number in `x` as + the starting point of the definite interval. + order + Order of differentiation or integration. When `order = 0`, function + value will be returned. When `order > 0` function derviative will + be returned. When `order < 0`, function integral will be returned. + Default is `0`. + + Returns + ------- + describe + Return function values, derivatives or definite integrals. + + Raises + ------ + AttributeError + Raised when the function implementation is not provided. + ValueError + Raised when `x` is not a scalar, 1d array or 2d array with two rows. + ValueErorr + Raised when `order >= 0` and `x` is a 2d array. + ValueError + Raised when `order < 0` and `any(x[0] > x[1])`. + + """ + x, order, isscalar = self._check_args(x, order) + if x.size == 0: + return np.empty(shape=x.shape, dtype=x.dtype) + result = self._fun(x, order) + if isscalar: + return result[0] + return result + + def append(self, other: XFunction, sep: BoundaryPoint) -> XFunction: + """Splice with another instance of ``XFunction`` to create a new + ``XFunction``. + + Parameters + ---------- + other + Another ``XFunction`` after the current function. + sep + The boundary point to separate two functions, before is the current + function and after the the ``other`` function. + + """ + lfun, rfun = self._fun, other._fun + + def fun(x: NDArray, order: int = 0) -> NDArray: + left = x <= sep[0] if sep[1] else x < sep[0] + + if order >= 0: + return np.where(left, lfun(x, order), rfun(x, order)) + + lboth, rboth = left.all(axis=0), (~left).all(axis=0) + landr = (~lboth) & (~rboth) + + result = np.zeros(x.shape[1], dtype=x.dtype) + result[lboth] = lfun(x[:, lboth], order) + result[rboth] = rfun(x[:, rboth], order) + + if landr.any(): + lx = np.insert(x[np.ix_([0], landr)], 1, sep[0], axis=0) + rx = np.insert(x[np.ix_([1], landr)], 0, sep[0], axis=0) + dx = x[1][landr] - sep[0] + + for i in range(1, -order): + result[landr] += lfun(lx, order + i) * (dx**i / factorial(i)) + result[landr] += lfun(lx, order) + rfun(rx, order) + return result + + return XFunction(fun) + + +class BundleXFunction(XFunction): + """This is one implementation of the ``XFunction``, it takes the value, + derivative and definite integral function and bundle them together as a + ``XFunction``. + + Parameters + ---------- + params + This is the parameters that is needed for the value, derivatives and + the definitely integral function. + val_fun + Value function. + der_fun + Derviative function. + int_fun + Defintie integral function. + + """ + + def __init__( + self, + params: tuple, + val_fun: RawVFunction, + der_fun: RawDFunction, + int_fun: RawIFunction, + ) -> None: + self.params = params + self.val_fun = partial(val_fun, params) + self.der_fun = partial(der_fun, params) + self.int_fun = partial(int_fun, params) + + def fun(x: NDArray, order: int = 0) -> NDArray: + if order == 0: + return self.val_fun(x) + if order > 0: + return self.der_fun(x, order) + dx = np.diff(x, axis=0)[0] + val = self.int_fun(x[1], order) + for i in range(-order): + val -= self.int_fun(x[0], order + i) * (dx**i / factorial(i)) + return val + + super().__init__(fun) + + +class BasisXFunction(XFunction): + """This is one implementation of ``XFunction`` by taking in a set of + instances of ``XFunction`` as basis functions. And the linear combination + coefficients to provide function value, derivative and definite integral. + + Parameters + ---------- + basis_funs + A set of instances of ``XFunction`` as basis functions. + coef + Coefficients for the linearly combine the basis functions. + + """ + + coef = property(attrgetter("_coef")) + + def __init__( + self, basis_funs: tuple[XFunction, ...], coef: NDArray | None = None + ) -> None: + if not all(isinstance(fun, XFunction) for fun in basis_funs): + raise TypeError("basis functions must all be instances of 'XFunction'") + self.basis_funs = tuple(basis_funs) + self.coef = coef + + def fun(x: NDArray, order: int = 0) -> NDArray: + if self.coef is None: + raise ValueError( + "please provide the coefficients for the basis functions" + ) + design_mat = self.get_design_mat(x, order=order, check_args=False) + return design_mat.dot(self.coef) + + super().__init__(fun) + + @coef.setter + def coef(self, coef: NDArray | None) -> None: + if coef is not None: + coef = np.asarray(coef, dtype=float).ravel() + if coef.size != len(self): + raise ValueError( + "number of coeffcients does not match number of basis functions" + ) + self._coef = coef + + def get_design_mat( + self, x: NDArray, order: int = 0, check_args: bool = True + ) -> NDArray: + """Provide design matrix from the set of basis functions. + + Parameters + ---------- + x + Data points + order + Order of differentiation/integration. + check_args + If ``True`` it will check and parse the arguments. + + Returns + ------- + describe + Design matrix with dimention number of data points by number of + basis functions. + + """ + if check_args: + x, order, _ = self._check_args(x, order) + return np.vstack([xfun._fun(x, order) for xfun in self.basis_funs]).T + + def __len__(self) -> int: + """Number of basis functions.""" + return len(self.basis_funs) diff --git a/src/xspline/xspl.py b/src/xspline/xspl.py new file mode 100644 index 0000000..4efa437 --- /dev/null +++ b/src/xspline/xspl.py @@ -0,0 +1,79 @@ +from typing import Optional + +from xspline.bspl import clear_bspl_cache, get_bspl_funs +from xspline.poly import get_poly_fun +from xspline.xfunction import BasisXFunction +from xspline.typing import NDArray + + +class XSpline(BasisXFunction): + """Main class for xspline functions. + + Parameters + ---------- + knots + Knots of the spline. + degree + Degree of the spline. + ldegree + Left extrapolation polynomial degree. + rdegree + Right extrapolation polynomial degree. + coef + The coefficients for linear combining the spline basis. + + """ + + def __init__( + self, + knots: tuple[float, ...], + degree: int, + ldegree: Optional[int] = None, + rdegree: Optional[int] = None, + coef: Optional[NDArray] = None, + ) -> None: + # validate inputs + knots, degree = tuple(sorted(map(float, knots))), int(degree) + if len(set(knots)) < 2: + raise ValueError("please provide at least provide 2 distinct knots") + if degree < 0: + raise ValueError("degree must be nonnegative") + ldegree = min(int(degree if ldegree is None else ldegree), degree) + rdegree = min(int(degree if rdegree is None else rdegree), degree) + + # create basis functions + mfuns = get_bspl_funs(knots, degree) + lfuns = tuple(get_poly_fun(fun, knots[0], ldegree) for fun in mfuns) + rfuns = tuple(get_poly_fun(fun, knots[-1], rdegree) for fun in mfuns) + funs = tuple( + lfun.append(mfun, (knots[0], False)).append(rfun, (knots[-1], True)) + for lfun, mfun, rfun in zip(lfuns, mfuns, rfuns) + ) + + self.knots, self.degree = knots, degree + self.ldegree, self.rdegree = ldegree, rdegree + super().__init__(funs, coef=coef) + + def get_design_mat( + self, x: NDArray, order: int = 0, check_args: bool = True + ) -> NDArray: + """Create design matrix from spline basis functions. + + Parameters + ---------- + x + Data points. + order + Order of differentiation/integration. + check_args + If ``True``, it will automatically check and parse the arguments. + + Returns + ------- + describe + Design matrix from spline basis functions. + + """ + design_mat = super().get_design_mat(x, order, check_args) + clear_bspl_cache() + return design_mat diff --git a/tests/test_basis_xfunction.py b/tests/test_basis_xfunction.py new file mode 100644 index 0000000..0197c34 --- /dev/null +++ b/tests/test_basis_xfunction.py @@ -0,0 +1,43 @@ +import numpy as np +import pytest +from xspline.indi import Indi +from xspline.xfunction import BasisXFunction + + +@pytest.fixture +def xfun(): + fun0 = Indi(((0.0, True), (1.0, True))) + fun1 = Indi(((1.5, True), (2.5, True))) + + return BasisXFunction([fun0, fun1]) + + +def test_len(xfun): + assert len(xfun) == 2 + + +def test_coef_error(xfun): + coef = [1, 2, 3] + with pytest.raises(ValueError): + xfun.coef = coef + + +@pytest.mark.parametrize("order", [-1, 0, 1]) +def test_get_design_mat(xfun, order): + x = np.linspace(-0.5, 3.0, 101) + design_mat = xfun.get_design_mat(x, order=order) + assert design_mat.shape == (x.size, len(xfun)) + + +@pytest.mark.parametrize("order", [-1, 0, 1]) +def test_fun(xfun, order): + xfun.coef = [1, 1] + x = np.linspace(-0.5, 3.0, 101) + result = xfun(x, order=order) + assert result.shape == x.shape + + +def test_fun_error(xfun): + x = np.linspace(-0.5, 3.0, 101) + with pytest.raises(ValueError): + xfun(x) diff --git a/tests/test_bspl.py b/tests/test_bspl.py new file mode 100644 index 0000000..829e2f8 --- /dev/null +++ b/tests/test_bspl.py @@ -0,0 +1,75 @@ +import numpy as np +import pytest +from xspline.bspl import Bspl, clear_bspl_cache + +knots = (0.0, 0.5, 1.0) +degree = 2 +x = np.linspace(knots[0], knots[1], 101) +poly_params = [ + [(4.0, -4.0, 1.0), (0.0, 0.0, 0.0)], + [(-6.0, 4.0, 0.0), (2.0, -4.0, 2.0)], + [(2.0, 0.0, 0.0), (-6.0, 8.0, -2.0)], + [(0.0, 0.0, 0.0), (4.0, -4.0, 1.0)], +] +indices = list(range(-degree, len(knots) - 1)) + + +@pytest.mark.parametrize("index", indices) +def test_bspl_val(index): + bspl = Bspl((knots, degree, index)) + sub_poly_params = poly_params[index + degree] + + my_val = bspl(x) + tr_val = np.where( + x < knots[1], + np.polyval(sub_poly_params[0], x), + np.polyval(sub_poly_params[1], x) + ) + assert np.allclose(my_val, tr_val) + + clear_bspl_cache() + + +@pytest.mark.parametrize("index", indices) +def test_bspl_der(index): + bspl = Bspl((knots, degree, index)) + sub_poly_params = list( + map(lambda c: np.polyder(c, 1), poly_params[index + degree]) + ) + + my_val = bspl(x, order=1) + tr_val = np.where( + x < knots[1], + np.polyval(sub_poly_params[0], x), + np.polyval(sub_poly_params[1], x) + ) + assert np.allclose(my_val, tr_val) + + clear_bspl_cache() + + +@pytest.mark.parametrize("index", indices) +def test_bspl_int(index): + bspl = Bspl((knots, degree, index)) + sub_poly_params = list( + map(lambda c: np.polyint(c, 1), poly_params[index + degree]) + ) + offsets = [ + -np.polyval(sub_poly_params[0], 0.25), + -np.polyval(sub_poly_params[1], knots[1]) + + ( + np.polyval(sub_poly_params[0], knots[1]) - + np.polyval(sub_poly_params[0], 0.25) + ) + ] + a = 0.25 + x = np.vstack([np.repeat(a, 76), np.linspace(a, knots[2], 76)]) + my_val = bspl(x, order=-1) + tr_val = np.where( + x[1] < knots[1], + np.polyval(sub_poly_params[0], x[1]) + offsets[0], + np.polyval(sub_poly_params[1], x[1]) + offsets[1], + ) + assert np.allclose(my_val, tr_val) + + clear_bspl_cache() diff --git a/tests/test_core.py b/tests/test_core.py deleted file mode 100644 index 5332c3d..0000000 --- a/tests/test_core.py +++ /dev/null @@ -1,204 +0,0 @@ -# -*- coding: utf-8 -*- -""" - test - ~~~~~~~~~ - - unit tests for xspline.core -""" -import numpy as np -import pytest -from xspline import core - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("degree", [0, 1, 2]) -@pytest.mark.parametrize("idx", [0, 1, -1]) -@pytest.mark.parametrize("l_extra", [False, True]) -@pytest.mark.parametrize("r_extra", [False, True]) -def test_bspline_domain(knots, degree, idx, l_extra, r_extra): - my_domain = core.bspline_domain(knots, degree, idx, - l_extra=l_extra, - r_extra=r_extra) - if idx == 0: - tr_domain = knots[:2].copy() - if l_extra: - tr_domain[0] = -np.inf - elif idx == -1: - tr_domain = knots[-2:].copy() - if r_extra: - tr_domain[1] = np.inf - elif idx == 1: - if degree == 0: - tr_domain = np.array([knots[1], knots[2]]) - else: - tr_domain = np.array([knots[0], knots[2]]) - - assert tr_domain[0] == my_domain[0] - assert tr_domain[1] == my_domain[1] - - -@pytest.mark.parametrize("knots", [np.arange(5).astype(float)]) -@pytest.mark.parametrize("degree", [0, 1]) -def test_bspline_domain_l_extra(knots, degree): - idx = 0 - result = np.array([-np.inf, knots[1]]) - my_result = core.bspline_domain(knots, degree, idx, l_extra=True) - assert np.isneginf(my_result[0]) - assert np.abs(my_result[1] - result[1]) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.arange(5).astype(float)]) -@pytest.mark.parametrize("degree", [0, 1]) -def test_bspline_domain_r_extra(knots, degree): - idx = knots.size + degree - 2 - result = np.array([knots[-2], np.inf]) - my_result = core.bspline_domain(knots, degree, idx, r_extra=True) - assert np.isposinf(my_result[1]) - assert np.abs(my_result[0] - result[0]) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("degree", [1]) -@pytest.mark.parametrize("idx", [0, 1]) -@pytest.mark.parametrize("x", [np.linspace(0.0, 1.0, 101)]) -def test_bspline_fun(x, knots, degree, idx): - my_y = core.bspline_fun(x, knots, degree, idx) - if idx == 0: - tr_y = np.maximum((knots[1] - x)/knots[1], 0.0) - else: - tr_y = np.zeros(x.size) - idx1 = (x >= knots[0]) & (x < knots[1]) - idx2 = (x >= knots[1]) & (x < knots[2]) - tr_y[idx1] = x[idx1]/knots[1] - tr_y[idx2] = (knots[2] - x[idx2])/(knots[2] - knots[1]) - - assert np.linalg.norm(tr_y - my_y) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("x", [np.linspace(-1.0, -0.1, 10)]) -@pytest.mark.parametrize("l_extra", [True, False]) -def test_bspline_fun_l_extra(x, knots, l_extra): - degree = 0 - idx = 0 - my_y = core.bspline_fun(x, knots, degree, idx, l_extra=l_extra) - if l_extra: - tr_y = 1.0 - else: - tr_y = 0.0 - assert np.linalg.norm(my_y - tr_y) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("x", [np.linspace(1.1, 2.0, 10)]) -@pytest.mark.parametrize("r_extra", [True, False]) -def test_bspline_fun_r_extra(x, knots, r_extra): - degree = 0 - idx = -1 - my_y = core.bspline_fun(x, knots, degree, idx, r_extra=r_extra) - if r_extra: - tr_y = 1.0 - else: - tr_y = 0.0 - assert np.linalg.norm(my_y - tr_y) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("degree", [1]) -@pytest.mark.parametrize("order", [1]) -@pytest.mark.parametrize("idx", [0, 1]) -@pytest.mark.parametrize("x", [np.linspace(0.0, 1.0, 101)]) -def test_bspline_dfun(x, knots, degree, order, idx): - my_dy = core.bspline_dfun(x, knots, degree, order, idx) - tr_dy = np.zeros(x.size) - idx1 = (x >= knots[0]) & (x < knots[1]) - idx2 = (x >= knots[1]) & (x < knots[2]) - if idx == 0: - tr_dy[idx1] = -1.0/knots[1] - else: - tr_dy[idx1] = 1.0/knots[1] - tr_dy[idx2] = -1.0/(knots[2] - knots[1]) - - assert np.linalg.norm(tr_dy - my_dy) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("x", [np.linspace(-1.0, -0.1, 10)]) -@pytest.mark.parametrize("l_extra", [True, False]) -def test_bspline_dfun_l_extra(x, knots, l_extra): - degree = 1 - idx = 0 - order = 1 - my_dy = core.bspline_dfun(x, knots, degree, order, idx, - l_extra=l_extra) - if l_extra: - tr_dy = -1.0/knots[1] - else: - tr_dy = 0.0 - - assert np.linalg.norm(my_dy - tr_dy) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("x", [np.linspace(1.1, 2.0, 10)]) -@pytest.mark.parametrize("r_extra", [True, False]) -def test_bspline_dfun_r_extra(x, knots, r_extra): - degree = 1 - idx = -1 - order = 1 - my_dy = core.bspline_dfun(x, knots, degree, order, idx, - r_extra=r_extra) - if r_extra: - tr_dy = 1.0/(knots[-1] - knots[-2]) - else: - tr_dy = 0.0 - assert np.linalg.norm(my_dy - tr_dy) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("degree", [1]) -@pytest.mark.parametrize("order", [1]) -@pytest.mark.parametrize("idx", [0]) -@pytest.mark.parametrize("x", [np.linspace(0.0, 1.0, 101)]) -def test_bspline_ifun(x, knots, degree, order, idx): - my_iy = core.bspline_ifun(knots[0], x, knots, degree, order, idx) - tr_iy = np.zeros(x.size) - idx1 = (x >= knots[0]) & (x <= knots[1]) - - tr_iy[idx1] = x[idx1] - 0.5/knots[1]*x[idx1]**2 - tr_iy[~idx1] = tr_iy[idx1][-1] - - assert np.linalg.norm(tr_iy - my_iy) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("x", [np.linspace(-1.0, -0.1, 10)]) -@pytest.mark.parametrize("l_extra", [True, False]) -def test_bspline_ifun_l_extra(x, knots, l_extra): - degree = 0 - idx = 0 - order = 1 - my_iy = core.bspline_ifun(x[0], x, knots, degree, order, idx, - l_extra=l_extra) - if l_extra: - tr_iy = x - x[0] - else: - tr_iy = 0.0 - - assert np.linalg.norm(my_iy - tr_iy) < 1e-10 - - -@pytest.mark.parametrize("knots", [np.linspace(0.0, 1.0, 5)]) -@pytest.mark.parametrize("x", [np.linspace(1.1, 2.0, 10)]) -@pytest.mark.parametrize("r_extra", [True, False]) -def test_bspline_ifun_r_extra(x, knots, r_extra): - degree = 0 - idx = -1 - order = 1 - my_iy = core.bspline_ifun(x[0], x, knots, degree, order, idx, - r_extra=r_extra) - if r_extra: - tr_iy = x - x[0] - else: - tr_iy = 0.0 - assert np.linalg.norm(my_iy - tr_iy) < 1e-10 diff --git a/tests/test_indi.py b/tests/test_indi.py new file mode 100644 index 0000000..dc248a4 --- /dev/null +++ b/tests/test_indi.py @@ -0,0 +1,44 @@ +import numpy as np +import pytest +from numpy.typing import NDArray +from xspline.indi import Indi + +# define parameters +params = [((1.0, True), (2.0, False))] +x = [np.linspace(1, 2, 101)] +order = [-2, -1, 0, 1, 2] + + +def truth(x: NDArray, order: int) -> NDArray: + result = np.zeros(x.size) + index0 = (x >= 1.0) & (x < 2) + index1 = (x >= 2.0) + if order == 0: + result[index0] = 1.0 + return result + if order == 1: + return result + if order == 2: + return result + if order == -1: + result[index0] = x[index0] - 1.0 + result[index1] = 1.0 + return result + if order == -2: + result[index0] = 0.5*x[index0]**2 - x[index0] + 0.5 + result[index1] = x[index1] - 1.5 + return result + raise ValueError("not support value for 'order'") + + +@pytest.mark.parametrize("params", params) +@pytest.mark.parametrize("x", x) +@pytest.mark.parametrize("order", order) +def test_indi(params, x, order): + indi = Indi(params) + z = x + if order < 0: + z = np.vstack([np.ones(x.size), x]) + assert np.allclose( + indi(z, order), truth(x, order) + ) diff --git a/tests/test_ndxspline.py b/tests/test_ndxspline.py deleted file mode 100644 index f9e859e..0000000 --- a/tests/test_ndxspline.py +++ /dev/null @@ -1,115 +0,0 @@ -# -*- coding: utf-8 -*- -""" - test_xspline - ~~~~~~~~~~~~ - - Test XSpline class. -""" -import numpy as np -import pytest -from xspline import NDXSpline - -@pytest.fixture -def spline(): - knots = [ - np.linspace(0.0, 1.0, 3), - np.linspace(0.0, 1.0, 2) - ] - degree = [2, 3] - return NDXSpline(2, knots, degree) - - -def test_ndim(spline): - assert spline.ndim == 2 - - -def test_num_knots(spline): - assert all([spline.num_knots_list[i] == spline.knots_list[i].size - for i in range(spline.ndim)]) - assert spline.num_knots == np.prod(spline.num_knots_list) - - -def test_num_intervals(spline): - assert all([spline.num_intervals_list[i] == spline.knots_list[i].size - 1 - for i in range(spline.ndim)]) - assert spline.num_intervals == np.prod(spline.num_intervals_list) - - -def test_num_spline_bases(spline): - assert all([spline.num_spline_bases_list[i] == - spline.spline_list[i].num_spline_bases - for i in range(spline.ndim)]) - assert spline.num_spline_bases == np.prod(spline.num_spline_bases_list) - - -@pytest.mark.parametrize('l_extra_list', [None, - [False, False], - [True, False], - [False, True], - [True, True]]) -@pytest.mark.parametrize('r_extra_list', [None, - [False, False], - [True, False], - [False, True], - [True, True]]) -@pytest.mark.parametrize('is_grid', [True, False]) -def test_design_mat(spline, is_grid, l_extra_list, r_extra_list): - x = np.linspace(0.0, 1.0, 11) - y = np.linspace(0.0, 1.0, 11) - mat = spline.design_mat([x, y], - is_grid=is_grid, - l_extra_list=l_extra_list, - r_extra_list=r_extra_list) - if is_grid: - assert mat.shape == (x.size*y.size, spline.num_spline_bases) - else: - assert mat.shape == (x.size, spline.num_spline_bases) - - -@pytest.mark.parametrize('l_extra_list', [None, - [False, False], - [True, False], - [False, True], - [True, True]]) -@pytest.mark.parametrize('r_extra_list', [None, - [False, False], - [True, False], - [False, True], - [True, True]]) -@pytest.mark.parametrize('is_grid', [True, False]) -def test_design_dmat(spline, is_grid, l_extra_list, r_extra_list): - x = np.linspace(0.0, 1.0, 11) - y = np.linspace(0.0, 1.0, 11) - dmat = spline.design_dmat([x, y], [1, 1], - is_grid=is_grid, - l_extra_list=l_extra_list, - r_extra_list=r_extra_list) - if is_grid: - assert dmat.shape == (x.size*y.size, spline.num_spline_bases) - else: - assert dmat.shape == (x.size, spline.num_spline_bases) - - -@pytest.mark.parametrize('l_extra_list', [None, - [False, False], - [True, False], - [False, True], - [True, True]]) -@pytest.mark.parametrize('r_extra_list', [None, - [False, False], - [True, False], - [False, True], - [True, True]]) -@pytest.mark.parametrize('is_grid', [True, False]) -def test_design_imat(spline, is_grid, l_extra_list, r_extra_list): - x = np.linspace(0.0, 1.0, 11) - y = np.linspace(0.0, 1.0, 11) - imat = spline.design_imat([x, y], [x, y], [1, 1], - is_grid=is_grid, - l_extra_list=l_extra_list, - r_extra_list=r_extra_list) - assert np.allclose(imat, 0.0) - if is_grid: - assert imat.shape == (x.size*y.size, spline.num_spline_bases) - else: - assert imat.shape == (x.size, spline.num_spline_bases) \ No newline at end of file diff --git a/tests/test_poly.py b/tests/test_poly.py new file mode 100644 index 0000000..a13ae16 --- /dev/null +++ b/tests/test_poly.py @@ -0,0 +1,36 @@ +import numpy as np +import pytest +from numpy.typing import NDArray +from xspline.poly import Poly + +# define the parameters +params = [(12, -6, 2)] +x = [np.linspace(1.0, 2, 101)] +order = [-2, -1, 0, 1, 2] + + +def truth(x: NDArray, order: int) -> NDArray: + if order == 0: + return 12*x**2 - 6*x + 2 + if order == 1: + return 24*x - 6 + if order == 2: + return np.repeat(24, x.size) + if order == -1: + return 4*x**3 - 3*x**2 + 2*x - 3 + if order == -2: + return x**4 - x**3 + x**2 - 3*x + 2 + raise ValueError("not support value for 'order'") + + +@pytest.mark.parametrize("params", params) +@pytest.mark.parametrize("x", x) +@pytest.mark.parametrize("order", order) +def test_poly(params, x, order): + poly = Poly(params) + z = x + if order < 0: + z = np.vstack([np.ones(x.size), x]) + assert np.allclose( + poly(z, order), truth(x, order) + ) diff --git a/tests/test_utils.py b/tests/test_utils.py deleted file mode 100644 index 8ef7eb7..0000000 --- a/tests/test_utils.py +++ /dev/null @@ -1,176 +0,0 @@ -# -*- coding: utf-8 -*- -""" - test_utils - ~~~~~~~~~~ - - unit tests for xspline.utils -""" -import numpy as np -import pytest -from xspline import utils - - -@pytest.mark.parametrize(("x", "l_close", "r_close", "result"), - [(0.5, False, False, 1.0), - (2.0, False, False, 0.0), - (1.0, False, False, 0.0), - (1.0, False, True, 1.0), - (0.0, False, False, 0.0), - (0.0, True, False, 1.0), - (np.repeat(0.5, 3), False, False, np.ones(3)), - (np.repeat(2.0, 3), False, False, np.zeros(3)), - (np.repeat(1.0, 3), False, False, np.zeros(3)), - (np.repeat(1.0, 3), False, True, np.ones(3)), - (np.repeat(0.0, 3), False, False, np.zeros(3)), - (np.repeat(0.0, 3), True, False, np.ones(3))]) -def test_utils_indicator_f(x, l_close, r_close, result): - b = np.array([0.0, 1.0]) - my_result = utils.indicator_f(x, b, l_close=l_close, r_close=r_close) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize(("x", "z", "fz", "dfz"), - [(1.0, 0.0, -1.0, 2.0), - (np.ones(5), 0.0, -1.0, 2.0), - (np.ones(5), - np.zeros(5), np.repeat(-1.0, 5), np.repeat(2.0, 5)), - (1.0, - np.zeros(5), np.repeat(-1.0, 5), np.repeat(2.0, 5))]) -def test_utils_linear_f(x, z, fz, dfz): - result = fz + dfz*(x - z) - my_result = utils.linear_f(x, z, fz, dfz) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize(("x", "b", "result"), - [(0.0, np.array([0.0, 1.0]), 0.0), - (1.0, np.array([0.0, 1.0]), 1.0), - (np.repeat(0.5, 5), np.array([0.0, 1.0]), - np.repeat(0.5, 5))]) -def test_utils_linear_lf(x, b, result): - my_result = utils.linear_lf(x, b) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize(("x", "b", "result"), - [(0.0, np.array([0.0, 1.0]), 1.0), - (1.0, np.array([0.0, 1.0]), 0.0), - (np.repeat(0.5, 5), np.array([0.0, 1.0]), - np.repeat(0.5, 5))]) -def test_utils_linear_rf(x, b, result): - my_result = utils.linear_rf(x, b) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize("a", [0.0, np.zeros(5)]) -@pytest.mark.parametrize("x", [1.0, np.ones(5)]) -@pytest.mark.parametrize("order", [0, 1]) -def test_utils_constant_if_0(a, x, order): - my_result = utils.constant_if(a, x, order, 0.0) - assert np.linalg.norm(my_result) < 1e-10 - - -@pytest.mark.parametrize("a", [0.0, np.zeros(5)]) -@pytest.mark.parametrize("x", [1.0, np.ones(5)]) -@pytest.mark.parametrize(("order", "result"), - [(1, 1.0), (2, 0.5)]) -def test_utils_constant_if_1(a, x, order, result): - my_result = utils.constant_if(a, x, order, 1.0) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize(("z", "fz", "dfz"), - [(1.0, 0.0, 2.0)]) -@pytest.mark.parametrize("a", [0.0, np.zeros(5)]) -@pytest.mark.parametrize("x", [1.0, np.ones(5)]) -@pytest.mark.parametrize(("order", "result"), - [(0, 0.0), (1, -1.0), (2, -2.0/3.0)]) -def test_utils_linear_if(a, x, order, z, fz, dfz, result): - my_result = utils.linear_if(a, x, order, z, fz, dfz) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize("a", [0.0, np.zeros(5)]) -@pytest.mark.parametrize("x", [1.0, np.ones(5)]) -@pytest.mark.parametrize("knots", [np.array([0.5])]) -@pytest.mark.parametrize(("order", "result"), - [(0, 2.0), (1, 1.5), (2, 0.625), (3, 3.0/16.0)]) -def test_utils_integrate_across_pieces(a, x, order, knots, result): - my_result = utils.integrate_across_pieces( - a, x, order, - [ - lambda *params: utils.constant_if(*params, 1.0), - lambda *params: utils.constant_if(*params, 2.0), - ], - knots - ) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize("a", [0.0, np.zeros(5)]) -@pytest.mark.parametrize("x", [1.0, np.ones(5)]) -@pytest.mark.parametrize("knots", [np.array([0.0, 0.5, 1.0])]) -@pytest.mark.parametrize(("order", "result"), - [(0, 2.0), (1, 1.5), (2, 0.625), (3, 3.0/16.0)]) -def test_utils_pieces_if(a, x, order, knots, result): - my_result = utils.pieces_if( - a, x, order, - [ - lambda *params: utils.constant_if(*params, 0.0), - lambda *params: utils.constant_if(*params, 1.0), - lambda *params: utils.constant_if(*params, 2.0), - lambda *params: utils.constant_if(*params, 0.0), - ], - knots - ) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize("a", [0.0, np.zeros(5)]) -@pytest.mark.parametrize("x", [1.0, np.ones(5)]) -@pytest.mark.parametrize("b", [np.array([0.5, 1.0])]) -@pytest.mark.parametrize(("order", "result"), - [(0, 0.0), (1, 0.5), (2, 0.125)]) -def test_utils_indicator_if(a, x, order, b, result): - my_result = utils.indicator_if(a, x, order, b) - assert np.linalg.norm(my_result - result) < 1e-10 - - -@pytest.mark.parametrize("size", [5, 10]) -def test_utils_seq_diff_mat(size): - x = np.random.randn(size) - y = np.array([x[i+1] - x[i] for i in range(size - 1)]) - - mat = utils.seq_diff_mat(size) - my_y = mat.dot(x) - - assert np.linalg.norm(my_y - y) < 1e-10 - - -@pytest.mark.parametrize("order", [0, 2, 4]) -@pytest.mark.parametrize("shape", [(5,), (2, 5), (1, 2, 3)]) -def test_utils_order_to_index(order, shape): - x = np.random.randn(*shape) - z = x.reshape(x.size, 1) - index = utils.order_to_index(order, shape) - assert z[order] == x[index] - - -@pytest.mark.parametrize("option", [True, False, None]) -@pytest.mark.parametrize("size", [2, 3]) -def test_utils_option_to_list(option, size): - result = not option - my_result = utils.option_to_list(option, size) - assert len(my_result) == size - assert all([~(my_result[i] ^ result) for i in range(size)]) - - -@pytest.mark.parametrize("args", - [(np.arange(2), np.arange(3)), - (np.arange(2), np.arange(3), np.arange(4))]) -def test_utils_outer_flatten(args): - result = np.prod(np.ix_(*args)) - result = result.reshape(result.size,) - my_result = utils.outer_flatten(*args) - - assert np.linalg.norm(my_result - result) < 1e-10 diff --git a/tests/test_xfunction.py b/tests/test_xfunction.py new file mode 100644 index 0000000..f9e1dc8 --- /dev/null +++ b/tests/test_xfunction.py @@ -0,0 +1,111 @@ +from math import factorial + +import numpy as np +import pytest +from numpy.typing import NDArray +from xspline.xfunction import XFunction + + +@pytest.mark.parametrize("order", [-2, -1, 0, 1, 2]) +def test_append(order): + # test append function + sep = (1.0, False) + + @XFunction + def fun0(x: NDArray, order: int = 0) -> NDArray: + return np.zeros(x.shape[-1], dtype=x.dtype) + + @XFunction + def fun1(x: NDArray, order: int = 0) -> NDArray: + if order == 0: + return np.ones(x.size, dtype=x.dtype) + if order > 0: + return np.zeros(x.size, dtype=x.dtype) + dx = x[1] - x[0] + return dx**(-order) / factorial(-order) + + my_fun = fun0.append(fun1, sep) + + def tr_fun(x: NDArray, order: int = 0) -> NDArray: + left = x <= sep[0] if sep[1] else x < sep[0] + if order == 0: + return np.where(left, 0.0, 1.0) + if order > 0: + return np.zeros(x.size, dtype=x.dtype) + z = np.maximum(0.0, x - sep[0]) + return z**(-order) / factorial(-order) + + x = np.linspace(0.0, 2.0, 101) + my_val = my_fun(x, order=order) + tr_val = tr_fun(x, order=order) + assert np.allclose(my_val, tr_val) + + +@pytest.fixture +def xfun(): + @XFunction + def fun(x: NDArray, order: int = 0) -> NDArray: + if order == 0: + return np.ones(x.size, dtype=x.dtype) + if order > 0: + return np.zeros(x.size, dtype=x.dtype) + dx = x[1] - x[0] + return dx**(-order) / factorial(-order) + return fun + + +def test_check_args_0(xfun): + """check value error, when x dimension is greater than 2 + """ + x = np.ones((2, 2, 2)) + with pytest.raises(ValueError): + xfun(x) + + +def test_check_args_1(xfun): + """check value error, when x dimension is 2 but x have more rows than 2 + """ + x = np.ones((3, 2)) + with pytest.raises(ValueError): + xfun(x) + + +@pytest.mark.parametrize("order", [0, 1]) +def test_check_args_2(xfun, order): + """check collapes behavior when x dimension is 2, and compute val or der + """ + x = np.ones((2, 3)) + with pytest.raises(ValueError): + xfun._check_args(x, order=order) + + +def test_check_args_3(xfun): + """check extend behavior when x dimension is 1, and compute integral + """ + x = np.ones(3) + x, _, _ = xfun._check_args(x, order=-1) + assert np.allclose(x, np.ones((2, 3))) + + +def test_check_args_4(xfun): + """check when x is a scalar + """ + x = 1.0 + x, _, isscalar = xfun._check_args(x, 0) + assert isscalar and x.shape == (1,) and np.allclose(x, 1.0) + + +def test_scalar_input_output(xfun): + """check when input is a scalar, output is also a scalar + """ + x = 1.0 + y = xfun(x) + assert np.isscalar(y) + + +def test_empty_input_output(xfun): + """check when input is an empty array, output is also an empty array. + """ + x = np.array([], dtype=float) + y = xfun(x) + assert y.size == 0 and y.shape == (0,) diff --git a/tests/test_xspl.py b/tests/test_xspl.py new file mode 100644 index 0000000..7d75f9d --- /dev/null +++ b/tests/test_xspl.py @@ -0,0 +1,35 @@ +import pytest +from xspline.xspl import XSpline + + +def test_knots_error(): + knots = (1, 1, 1) + degree = 2 + with pytest.raises(ValueError): + XSpline(knots, degree) + + +@pytest.mark.parametrize("degree", [-1, -1.0]) +def test_degree_error(degree): + knots = (0.0, 1.0) + with pytest.raises(ValueError): + XSpline(knots, degree) + + +@pytest.mark.parametrize(("degree", "sdegree", "result"), + [(2, None, 2), + (1, None, 1), + (2, 1, 1), + (2, 3, 2), + (2, -1, -1)]) +def test_ldegree_rdegree(degree, sdegree, result): + knots = (0.0, 1.0) + xspline = XSpline(knots, degree, ldegree=sdegree, rdegree=sdegree) + assert result == xspline.ldegree and result == xspline.rdegree + + +def test_len(): + knots = (0.0, 1.0) + degree = 3 + xspline = XSpline(knots, degree) + assert len(xspline) == degree + len(knots) - 1 diff --git a/tests/test_xspline.py b/tests/test_xspline.py deleted file mode 100644 index 54ae8af..0000000 --- a/tests/test_xspline.py +++ /dev/null @@ -1,185 +0,0 @@ -# -*- coding: utf-8 -*- -""" - test_xspline - ~~~~~~~~~~~~ - - Test XSpline class. -""" -import numpy as np -import pytest -from xspline.core import XSpline - - -@pytest.fixture -def knots(): - return np.linspace(0.0, 1.0, 5) - - -@pytest.fixture -def degree(): - return 1 - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("idx", [0, -1]) -def test_domain(knots, degree, idx, l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - my_domain = xs.domain(idx, l_extra=l_extra, r_extra=r_extra) - if idx == 0: - lb = -np.inf if l_extra else xs.knots[0] - else: - lb = xs.inner_knots[-2] - if idx == -1: - ub = np.inf if r_extra else xs.knots[-1] - else: - ub = xs.inner_knots[1] - - assert my_domain[0] == lb - assert my_domain[1] == ub - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("idx", [0, -1]) -def test_fun(knots, degree, idx, l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - if idx == 0: - x = np.linspace(xs.knots[0] - 1.0, xs.knots[1], 101) - else: - x = np.linspace(xs.knots[-2], xs.knots[-1], 101) - my_y = xs.fun(x, idx, l_extra=l_extra, r_extra=r_extra) - - if idx == 0: - tr_y = (xs.inner_knots[1] - x) / \ - (xs.inner_knots[1] - xs.inner_knots[0]) - if not l_extra: - tr_y[x < xs.knots[0]] = 0.0 - else: - tr_y = (x - xs.inner_knots[-2]) / \ - (xs.inner_knots[-1] - xs.inner_knots[-2]) - if not r_extra: - tr_y[x > xs.knots[-1]] = 0.0 - - assert np.linalg.norm(my_y - tr_y) < 1e-10 - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("order", [1]) -@pytest.mark.parametrize("idx", [0, -1]) -def test_dfun(knots, degree, order, idx, l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - if idx == 0: - x = np.linspace(xs.knots[0] - 1.0, xs.knots[1] - 1e-10, 101) - else: - x = np.linspace(xs.knots[-2], xs.knots[-1] - 1e-10, 101) - my_dy = xs.dfun(x, order, idx, l_extra=l_extra, r_extra=r_extra) - - if idx == 0: - tr_dy = np.repeat(1.0/(xs.inner_knots[0] - xs.inner_knots[1]), x.size) - if not l_extra: - tr_dy[x < xs.knots[0]] = 0.0 - else: - tr_dy = np.repeat(1.0/(xs.inner_knots[-1] - xs.inner_knots[-2]), x.size) - if not r_extra: - tr_dy[x > xs.knots[-1]] = 0.0 - - assert np.linalg.norm(my_dy - tr_dy) < 1e-10 - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("order", [1]) -@pytest.mark.parametrize("idx", [0, -1]) -def test_ifun(knots, degree, order, idx, l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - if idx == 0: - x = np.linspace(xs.knots[0] - 1.0, xs.knots[1], 101) - else: - x = np.linspace(xs.knots[-2], xs.knots[-1], 101) - my_iy = xs.ifun(x[0], x, order, idx, l_extra=l_extra, r_extra=r_extra) - - domain = xs.domain(idx, l_extra=l_extra, r_extra=r_extra) - lb = domain[0] - ub = domain[1] - v_idx = (x >= lb) & (x <= ub) - tr_iy = np.zeros(x.size) - - if idx == 0: - def ifun(a): - return 0.5*(xs.inner_knots[1] - a)**2 / \ - (xs.inner_knots[0] - xs.inner_knots[1]) - else: - def ifun(a): - return 0.5*(a - xs.inner_knots[-2])**2 / \ - (xs.inner_knots[-1] - xs.inner_knots[-2]) - - tr_iy[v_idx] = ifun(x[v_idx]) - ifun(lb) - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("x", [np.linspace(0.0, 1.0, 100), - np.linspace(-1.0, 0.5, 100), - np.linspace(0.5, 2.0, 100)]) -def test_design_mat(knots, degree, x, l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - mat = xs.design_mat(x, l_extra=l_extra, r_extra=r_extra) - - assert mat.shape == (x.size, xs.num_spline_bases) - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("order", [1]) -@pytest.mark.parametrize("x", [np.linspace(0.0, 1.0, 100), - np.linspace(-1.0, 0.5, 100), - np.linspace(0.5, 2.0, 100)]) -def test_design_dmat(knots, degree, order, x, - l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - dmat = xs.design_dmat(x, order, l_extra=l_extra, r_extra=r_extra) - - assert dmat.shape == (x.size, xs.num_spline_bases) - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -@pytest.mark.parametrize("l_extra", [True, False]) -@pytest.mark.parametrize("r_extra", [True, False]) -@pytest.mark.parametrize("order", [1]) -@pytest.mark.parametrize("x", [np.linspace(0.0, 1.0, 100), - np.linspace(-1.0, 0.5, 100), - np.linspace(0.5, 2.0, 100)]) -def test_design_imat(knots, degree, order, x, - l_linear, r_linear, l_extra, r_extra): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - imat = xs.design_imat(x[0], x, order, l_extra=l_extra, r_extra=r_extra) - - assert imat.shape == (x.size, xs.num_spline_bases) - - -@pytest.mark.parametrize("l_linear", [True, False]) -@pytest.mark.parametrize("r_linear", [True, False]) -def test_last_dmat(knots, degree, l_linear, r_linear): - xs = XSpline(knots, degree, l_linear=l_linear, r_linear=r_linear) - my_last_dmat = xs.last_dmat() - - x = np.array([0.5*(xs.knots[i] + xs.knots[i+1]) - for i in range(xs.knots.size - 1)]) - tr_last_dmat = xs.design_dmat(x, 1) - - assert np.linalg.norm(my_last_dmat - tr_last_dmat) < 1e-10 \ No newline at end of file