HyperJet — Algorithmic Differentiation with Hyper-Dual numbers for Python and C++
+HyperJet — Algorithmic Differentiation with Hyper-Dual Numbers for C++ and Python
--- -A header-only library for algorithmic differentiation with hyper-dual numbers. Written in C++17 with an extensive Python interface. +A header-only C++23 library for algorithmic differentiation with hyper-dual numbers. Supports first- and second-order derivatives with both dense indexed variables (`DDScalar`) and sparse named variables (`SScalar`). Includes an extensive Python interface via [pybind11](https://github.com/pybind/pybind11). -[](https://pypi.org/project/hyperjet) [](https://zenodo.org/badge/latestdoi/165487832) [](https://github.com/oberbichler/HyperJet/actions)    +[](https://pypi.org/project/hyperjet) + + +[](https://github.com/oberbichler/HyperJet/actions/workflows/test-python.yml) +[](https://github.com/oberbichler/HyperJet/actions/workflows/test-cpp.yml) +[](https://github.com/oberbichler/HyperJet/blob/main/LICENSE) +[](https://zenodo.org/badge/latestdoi/165487832) ## Installation @@ -14,6 +20,8 @@ A header-only library for algorithmic differentiation with hyper-dual numbers. W pip install hyperjet ``` +**Requirements:** Python ≥ 3.11, NumPy ≥ 2.0 + ## Quickstart Import the module: @@ -22,47 +30,43 @@ Import the module: import hyperjet as hj ``` -Create a set of variables e.g. `x=3` and `y=6`: +Create a set of hyper-dual variables, e.g. `x=3` and `y=6`: ```python x, y = hj.variables([3, 6]) ``` -`x` and `y` are hyper-dual numbers. This is indicated by the postfix `hj`: +`x` and `y` are second-order hyper-dual numbers, indicated by the `hj` postfix: ```python x >>> 3hj ``` -Get the value as a simple `float`: +Get the value as a plain `float`: ```python x.f >>> 3 ``` -The hyper-dual number stores the derivatives as a numpy array. - -Get the first order derivatives (Gradient) of a hyper-dual number: +Get the first-order derivatives (gradient) as a NumPy array: ```python -x.g # = [dx/dx, dx/dy] +x.g # [dx/dx, dx/dy] >>> array([1., 0.]) ``` -Get the second order derivatives (Hessian matrix): +Get the second-order derivatives (Hessian matrix): ```python -x.hm() # = [[d^2 x/ dx**2 , d^2 x/(dx*dy)], - # [d^2 x/(dx*dy), d^2 x/ dy**2 ]] +x.hm() # [[d²x/dx², d²x/(dx·dy)], + # [d²x/(dx·dy), d²x/dy²]] >>> array([[0., 0.], [0., 0.]]) ``` -For a simple variable these derivatives are trivial. - -Now do some computations: +For a single variable these derivatives are trivial. Now do a computation: ```python f = (x * y) / (x - y) @@ -70,27 +74,59 @@ f >>> -6hj ``` -The result is again a hyper-dual number. - -Get the first order derivatives of `f` with respect to `x` and `y`: +The result is again a hyper-dual number carrying all derivatives: ```python -f.g # = [df/dx, df/dy] +f.g # [df/dx, df/dy] >>> array([-4., 1.]) + +f.hm() # [[d²f/dx², d²f/(dx·dy)], + # [d²f/(dx·dy), d²f/dy²]] +>>> array([[-2.66666667, 1.33333333], + [ 1.33333333, -0.66666667]]) ``` -Get the second order derivatives of `f`: +## Types + +HyperJet provides two families of scalar types: + +### `DDScalar` — Dense dual numbers with indexed variables + +Stores derivatives in dense arrays. Supports first-order (gradient only) and second-order (gradient + Hessian) derivatives. Available in **static** variants with a compile-time-fixed number of variables, and a **dynamic** variant for arbitrary sizes. + +| Python type | Order | Variables | C++ type | +|-------------|-------|-----------|----------| +| `DScalar` | 1 | dynamic | `DDScalar<1, double>` | +| `DDScalar` | 2 | dynamic | `DDScalar<2, double>` | +| `D3Scalar` | 1 | 3 (static) | `DDScalar<1, double, 3>` | +| `DD3Scalar` | 2 | 3 (static) | `DDScalar<2, double, 3>` | + +Static variants (`D1Scalar`–`D15Scalar`, `DD1Scalar`–`DD15Scalar`) avoid heap allocation and enable better compiler optimization. The dynamic variants (`DScalar`, `DDScalar`) accept any number of variables at runtime. + +The convenience function `hj.variables(values, order=2)` automatically selects the appropriate static type when the number of variables is ≤ 15, and falls back to the dynamic variant otherwise. + +### `SScalar` — Sparse dual numbers with named variables + +Stores first-order derivatives in a sparse map keyed by variable name (string). Useful when variables are identified by name rather than index, or when only a small subset of derivatives is non-zero. ```python -f.hm() # = [[d^2 f/ dx**2 , d^2 f/(dx*dy)], - # [d^2 f/(dx*dy), d^2 f/ dy**2 ]] ->>> array([[-2.66666667, 1.33333333], - [ 1.33333333, -0.66666667]]) +x = hj.SScalar.variable("x", 3.0) +y = hj.SScalar.variable("y", 6.0) + +f = (x * y) / (x - y) +f.f # value +>>> -6.0 +f.d("x") # df/dx +>>> -4.0 +f.d("y") # df/dy +>>> 1.0 ``` -You can use numpy to perform vector and matrix operations. +## NumPy Integration + +HyperJet scalars work with NumPy for vector and matrix operations. -Compute the nomalized cross product of two vectors `u = [1, 2, 2]` and `v = [4, 1, -1]` with hyper-dual numbers: +Compute the normalized cross product of `u = [1, 2, 2]` and `v = [4, 1, -1]` with full second-order derivatives: ```python import numpy as np @@ -98,8 +134,8 @@ import numpy as np variables = hj.DDScalar.variables([1, 2, 2, 4, 1, -1]) -u = np.array(variables[:3]) # = [1hj, 2hj, 2hj] -v = np.array(variables[3:]) # = [4hj, 1hj, -1hj] +u = np.array(variables[:3]) # [1hj, 2hj, 2hj] +v = np.array(variables[3:]) # [4hj, 1hj, -1hj] normal = np.cross(u, v) normal /= np.linalg.norm(normal) @@ -107,9 +143,7 @@ normal >>> array([-0.331042hj, 0.744845hj, -0.579324hj], dtype=object) ``` -The result is a three-dimensional numpy array containing hyper-dual numbers. - -Get the value and derivatives of the x-component: +The result is a three-dimensional NumPy array of hyper-dual numbers. Extract value and derivatives from any component: ```python normal[0].f @@ -127,6 +161,65 @@ normal[0].hm() [-0.02335746, 0.04858632, -0.03690759, -0.01546811, -0.02868433, 0.03641839]]) ``` +## C++ Usage + +HyperJet is a single header-only library. Add `include/` to your include path and use C++23: + +```cpp +#include