import 'autodiff' third party library

this is version v1.0.3
This commit is contained in:
Arne Morten Kvarving 2023-09-20 07:05:03 +02:00
parent a8a4cd780e
commit c9c13d888c
24 changed files with 6697 additions and 0 deletions

107
3rdparty/autodiff/CHANGELOG.md vendored Normal file
View File

@ -0,0 +1,107 @@
# Change Log
This file documents important changes to this project, which uses [Semantic Versioning](http://semver.org/).
#[v0.5](https://github.com/autodiff/autodiff/releases/tag/v0.5.0) (June 17, 2019)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v0.5.0...v0.4.2)
This release introduces a **BREAKING CHANGE**!
From now on, methods `derivative`, `gradient`, and `jacobian` require the
use of auxiliary functions `wrt` and the newly introduced one `at`.
Examples:
~~~c++
// f = f(x)
double dudx = derivative(f, wrt(x), at(x));
~~~
~~~c++
// f = f(x, y, z)
double dudx = derivative(f, wrt(x), at(x, y, z));
double dudy = derivative(f, wrt(y), at(x, y, z));
double dudz = derivative(f, wrt(z), at(x, y, z));
~~~
~~~c++
// f = f(x), scalar function, where x is an Eigen vector
VectorXd g = gradient(f, wrt(x), at(x));
// Compuring gradient with respect to only some variables
VectorXd gpartial = gradient(f, wrt(x.tail(5)), at(x));
~~~
~~~c++
// F = F(x), vector function, where x is an Eigen vector
MatrixXd J = jacobian(f, wrt(x), at(x));
// F = F(x, p), vector function with params, where x and p are Eigen vectors
MatrixXd Jx = jacobian(f, wrt(x), at(x, p));
MatrixXd Jp = jacobian(f, wrt(p), at(x, p));
// Compuring Jacobian with respect to only some variables
MatrixXd Jpartial = jacobian(f, wrt(x.tail(5)), at(x));
~~~
This release also permits one to retrieve the evaluated value of function during
a call to the methods `derivative`, `gradient`, and `jacobian`:
~~~c++
// f = f(x)
dual u;
double dudx = derivative(f, wrt(x), at(x), u);
~~~
~~~c++
// f = f(x), scalar function, where x is an Eigen vector
dual u;
VectorXd g = gradient(f, wrt(x), at(x), u);
~~~
~~~c++
// F = F(x), vector function, where x is an Eigen vector
VectorXdual F;
MatrixXd J = jacobian(f, wrt(x), at(x), F);
~~~
#[v0.4.2](https://github.com/autodiff/autodiff/releases/tag/v0.4.2) (Mar 28, 2019)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v0.4.2...v0.4.1)
This is to force conda-forge to produce a new version (now 0.4.2) since the last
one (0.4.1) did not work.
#[v0.4.1](https://github.com/autodiff/autodiff/releases/tag/v0.4.1) (Mar 26, 2019)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v0.4.1...v0.4.0)
This release fixes a bug in the computation of Jacobian matrices when the input
and output vectors in a vector-valued function have different dimensions (see
issue [#24](https://github.com/autodiff/autodiff/issues/24)).
#[v0.4.0](https://github.com/autodiff/autodiff/releases/tag/v0.4.0) (Feb 20, 2019)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v0.4.0...v0.3.0)
This release contains changes that enable autodiff to be successfully compiled
in Linux, macOS, and Windows. Compilers tested were GCC 7, Clang 9, and Visual
Studio 2017. Compilers should support C++17.
#[v0.3.0](https://github.com/autodiff/autodiff/releases/tag/v0.3.0) (Feb 5, 2019)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v0.3.0...v0.2.0)
This release improves the forward mode algorithm to compute derivatives of any
order. It also introduces a proper website containing a more detailed
documentation of autodiff library: https://autodiff.github.io
#[v0.2.0](https://github.com/autodiff/autodiff/releases/tag/v0.2.0) (Jul 26, 2018)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v0.2.0...v0.1.0)
This release permits higher order derivatives to be computed with
`autodiff::gradx` function and it also enables the use of `autodiff::var` type
with Eigen vector and matrix types.
#[v0.1.0](https://github.com/autodiff/autodiff/releases/tag/v0.1.0) (Jul 19, 2018)
[Full Changelog](https://github.com/autodiff/autodiff/compare/v3.6.0...v3.6.1)
This is the first release of autodiff. Please note breaking changes might be
introduced, but not something that would take you more than a few minutes to
correct.

62
3rdparty/autodiff/CMakeLists.txt vendored Normal file
View File

@ -0,0 +1,62 @@
# The minimum required cmake version
cmake_minimum_required(VERSION 3.16)
# Add cmake modules of this project to the module path
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
# Use ccache to speed up repeated compilations
include(CCache)
# Name and details of the project
project(autodiff VERSION 1.0.3 LANGUAGES CXX)
# Enable parallel build if MSVC is used
add_compile_options($<$<CXX_COMPILER_ID:MSVC>:/MP>)
# Include the cmake variables with values for installation directories
include(GNUInstallDirs)
# Ensure proper configuration if in a conda environment
include(CondaAware)
# Define build options
option(AUTODIFF_BUILD_TESTS "Enable the compilation of the test files." ON)
option(AUTODIFF_BUILD_PYTHON "Enable the compilation of the python bindings." ON)
option(AUTODIFF_BUILD_EXAMPLES "Enable the compilation of the example files." ON)
option(AUTODIFF_BUILD_DOCS "Enable the build of the documentation and website." ON)
# Find eigen library
find_package(Eigen3 REQUIRED)
# autodiff requires a c++17 enabled compiler
set(CMAKE_CXX_STANDARD 17) # ensure cmake instructs compiler to use C++17
set(CMAKE_CXX_STANDARD_REQUIRED ON) # ensure the C++ standard given before is actually used
set(CMAKE_CXX_EXTENSIONS OFF) # avoid compile flags of the type -std=gnu++1z added by cmake
# Build the library (actually, just provide details about the library target, specify required compile options, etc. because autodiff is header-only)
add_subdirectory(autodiff)
# Build the tests
if(AUTODIFF_BUILD_TESTS)
add_subdirectory(tests)
endif()
# Build the python wrappers
if(AUTODIFF_BUILD_PYTHON)
set(PYBIND11_CPP_STANDARD -std=c++17) # Ensure pybind11 uses C++17 standard
find_package(pybind11 REQUIRED)
add_subdirectory(python)
endif()
# Build the examples
if(AUTODIFF_BUILD_EXAMPLES)
add_subdirectory(examples)
endif()
# Build the docs
if(AUTODIFF_BUILD_DOCS)
add_subdirectory(docs)
endif()
# Install the cmake config files that permit users to use find_package(autodiff)
include(autodiffInstallCMakeConfigFiles)

21
3rdparty/autodiff/LICENSE vendored Normal file
View File

@ -0,0 +1,21 @@
MIT License
Copyright (c) 20182021 Allan Leal
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

179
3rdparty/autodiff/README.md vendored Normal file
View File

@ -0,0 +1,179 @@
<a href="https://autodiff.github.io" target="_blank">
<img src='art/autodiff-header.svg' width='100%'>
</a>
[![Gitter chat](https://badges.gitter.im/autodiff/gitter.png)](https://gitter.im/autodiff/community)
![Linux build status](https://github.com/autodiff/autodiff/workflows/linux/badge.svg?branch=master)
![macOS build status](https://github.com/autodiff/autodiff/workflows/osx/badge.svg?branch=master)
![Windows build status](https://github.com/autodiff/autodiff/workflows/windows/badge.svg?branch=master)
# Overview
**autodiff** is a C++17 library that uses modern and advanced programming
techniques to enable automatic computation of derivatives in an efficient, easy,
and intuitive way.
We welcome you to use **autodiff** and recommend us any improvements you think
it is necessary. You may want to do so by chatting with us on our [Gitter
Community Channel][gitter] and/or by making proposals by creating a [GitHub
issue][issues].
## Demonstration
Consider the following function *f(x, y, z)*:
```c++
double f(double x, double y, double z)
{
return (x + y + z) * exp(x * y * z);
}
```
which we use use to evaluate the variable *u = f(x, y, z)*:
```c++
double x = 1.0;
double y = 2.0;
double z = 3.0;
double u = f(x, y, z);
```
How can we minimally transform this code so that not only *u*, but also its
derivatives *∂u/∂x*, *∂u/∂y*, and *∂u/∂z*, can be computed?
The next two sections present how this can be achieved using two automatic
differentiation algorithms implemented in **autodiff**: **forward mode** and
**reverse mode**.
### Forward mode
In a *forward mode automatic differentiation* algorithm, both output variables
and one or more of their derivatives are computed together. For example, the
function evaluation *f(x, y, z)* can be transformed in a way that it will not
only produce the value of *u*, *the output variable*, but also one or more of
its derivatives *(∂u/∂x, ∂u/∂y, ∂u/∂z)* with respect to the *input
variables* *(x, y, z)*.
Enabling forward automatic differentiation for the calculation of derivatives
using **autodiff** is relatively simple. For our previous function *f*, we only
need to replace the floating-point type `double` with `autodiff::dual` for both
input and output variables:
```c++
dual f(const dual& x, const dual& y, const dual& z)
{
return (x + y + z) * exp(x * y * z);
}
```
We can now compute the derivatives *∂u/∂x*, *∂u/∂y*, and *∂u/∂z* as follows:
```c++
dual x = 1.0;
dual y = 2.0;
dual z = 3.0;
dual u = f(x, y, z);
double dudx = derivative(f, wrt(x), at(x, y, z));
double dudy = derivative(f, wrt(y), at(x, y, z));
double dudz = derivative(f, wrt(z), at(x, y, z));
```
The auxiliary function `autodiff::wrt`, an acronym for **with respect to**,
is used to indicate which input variable *(x, y, z)* is the selected one to
compute the partial derivative of *f*. The auxiliary function `autodiff::at`
is used to indicate where (at which values of its parameters) the derivative
of *f* is evaluated.
### Reverse mode
In a *reverse mode automatic differentiation* algorithm, the output variable of
a function is evaluated first. During this function evaluation, all
mathematical operations between the input variables are *"recorded"* in an
*expression tree*. By traversing this tree from top-level (output variable as
the root node) to bottom-level (input variables as the leaf nodes), it is
possible to compute the contribution of each branch on the derivatives of the
output variable with respect to input variables.
<img
src='art/expression-tree-diagram.svg'
style='max-width:100%;'
title='Expression tree diagram.'>
Thus, a single pass in a reverse mode calculation **computes all derivatives**,
in contrast with forward mode, which requires one pass for each input variable.
Note, however, that it is possible to change the behavior of a forward pass so
that many (perhaps even all) derivatives of an output variable are computed
simultaneously (e.g., in a single forward pass, *∂u/∂x*, *∂u/∂y*, and *∂u/∂z*
are evaluated together with *u*, in contrast with three forward passes, each
one computing the individual derivatives).
Similar as before, we can use **autodiff** to enable reverse automatic
differentiation for our function *f* by simply replacing type `double` with
`autodiff::var` as follows:
```c++
var f(var x, var y, var z)
{
return (x + y + z) * exp(x * y * z);
}
```
The code below demonstrates how the derivatives *∂u/∂x*, *∂u/∂y*, and *∂u/∂z*
can be calculated:
```c++
var x = 1.0;
var y = 2.0;
var z = 3.0;
var u = f(x, y, z);
Derivatives dud = derivatives(u);
double dudx = dud(x);
double dudy = dud(y);
double dudz = dud(z);
```
The function `autodiff::derivatives` will traverse the expression tree stored
in variable `u` and compute all its derivatives with respect to the input
variables *(x, y, z)*, which are then stored in the object `dud`. The
derivative of `u` with respect to input variable `x` (i.e., *∂u/∂x*) can then
be extracted from `dud` using `dud(x)`. The operations `dud(x)`, `dud(y)`,
`dud(z)` involve no computations! Just extraction of derivatives previously
computed with a call to function `autodiff::derivatives`.
## Documentation
Check the documentation website for more details:
<a href="https://autodiff.github.io" target="_blank">
<img src='art/autodiff.github.io.svg' width='100%'>
</a>
# License
MIT License
Copyright (c) 20182023 Allan Leal
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
[gitter]: https://gitter.im/autodiff/community
[issues]: https://github.com/autodiff/autodiff/issues/new

View File

@ -0,0 +1,23 @@
# Create the autodiff interface library
add_library(autodiff INTERFACE)
# Add an alias autodiff::autodiff to the target library autodiff
add_library(autodiff::autodiff ALIAS autodiff)
# Set autodiff compilation features to be propagated to client code.
target_compile_features(autodiff INTERFACE cxx_std_17)
# Add the include paths to the Reaktoro target
target_include_directories(autodiff
INTERFACE
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}>
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
)
# Install autodiff interface library
install(TARGETS autodiff EXPORT autodiffTargets)
# Install autodiff header files
install(DIRECTORY ${PROJECT_SOURCE_DIR}/autodiff
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT headers
PATTERN "CMakeLists.txt" EXCLUDE)

View File

@ -0,0 +1,113 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// C++ includes
#include <cstddef>
namespace autodiff {
namespace detail {
/// The array containing pre calculated binomial coefficients up to order 50.
constexpr double binomialcoeffs_data[] = {
1,
1, 1,
1, 2, 1,
1, 3, 3, 1,
1, 4, 6, 4, 1,
1, 5, 10, 10, 5, 1,
1, 6, 15, 20, 15, 6, 1,
1, 7, 21, 35, 35, 21, 7, 1,
1, 8, 28, 56, 70, 56, 28, 8, 1,
1, 9, 36, 84, 126, 126, 84, 36, 9, 1,
1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1,
1, 11, 55, 165, 330, 462, 462, 330, 165, 55, 11, 1,
1, 12, 66, 220, 495, 792, 924, 792, 495, 220, 66, 12, 1,
1, 13, 78, 286, 715, 1287, 1716, 1716, 1287, 715, 286, 78, 13, 1,
1, 14, 91, 364, 1001, 2002, 3003, 3432, 3003, 2002, 1001, 364, 91, 14, 1,
1, 15, 105, 455, 1365, 3003, 5005, 6435, 6435, 5005, 3003, 1365, 455, 105, 15, 1,
1, 16, 120, 560, 1820, 4368, 8008, 11440, 12870, 11440, 8008, 4368, 1820, 560, 120, 16, 1,
1, 17, 136, 680, 2380, 6188, 12376, 19448, 24310, 24310, 19448, 12376, 6188, 2380, 680, 136, 17, 1,
1, 18, 153, 816, 3060, 8568, 18564, 31824, 43758, 48620, 43758, 31824, 18564, 8568, 3060, 816, 153, 18, 1,
1, 19, 171, 969, 3876, 11628, 27132, 50388, 75582, 92378, 92378, 75582, 50388, 27132, 11628, 3876, 969, 171, 19, 1,
1, 20, 190, 1140, 4845, 15504, 38760, 77520, 125970, 167960, 184756, 167960, 125970, 77520, 38760, 15504, 4845, 1140, 190, 20, 1,
1, 21, 210, 1330, 5985, 20349, 54264, 116280, 203490, 293930, 352716, 352716, 293930, 203490, 116280, 54264, 20349, 5985, 1330, 210, 21, 1,
1, 22, 231, 1540, 7315, 26334, 74613, 170544, 319770, 497420, 646646, 705432, 646646, 497420, 319770, 170544, 74613, 26334, 7315, 1540, 231, 22, 1,
1, 23, 253, 1771, 8855, 33649, 100947, 245157, 490314, 817190, 1144066, 1352078, 1352078, 1144066, 817190, 490314, 245157, 100947, 33649, 8855, 1771, 253, 23, 1,
1, 24, 276, 2024, 10626, 42504, 134596, 346104, 735471, 1307504, 1961256, 2496144, 2704156, 2496144, 1961256, 1307504, 735471, 346104, 134596, 42504, 10626, 2024, 276, 24, 1,
1, 25, 300, 2300, 12650, 53130, 177100, 480700, 1081575, 2042975, 3268760, 4457400, 5200300, 5200300, 4457400, 3268760, 2042975, 1081575, 480700, 177100, 53130, 12650, 2300, 300, 25, 1,
1, 26, 325, 2600, 14950, 65780, 230230, 657800, 1562275, 3124550, 5311735, 7726160, 9657700, 10400600, 9657700, 7726160, 5311735, 3124550, 1562275, 657800, 230230, 65780, 14950, 2600, 325, 26, 1,
1, 27, 351, 2925, 17550, 80730, 296010, 888030, 2220075, 4686825, 8436285, 13037895, 17383860, 20058300, 20058300, 17383860, 13037895, 8436285, 4686825, 2220075, 888030, 296010, 80730, 17550, 2925, 351, 27, 1,
1, 28, 378, 3276, 20475, 98280, 376740, 1184040, 3108105, 6906900, 13123110, 21474180, 30421755, 37442160, 40116600, 37442160, 30421755, 21474180, 13123110, 6906900, 3108105, 1184040, 376740, 98280, 20475, 3276, 378, 28, 1,
1, 29, 406, 3654, 23751, 118755, 475020, 1560780, 4292145, 10015005, 20030010, 34597290, 51895935, 67863915, 77558760, 77558760, 67863915, 51895935, 34597290, 20030010, 10015005, 4292145, 1560780, 475020, 118755, 23751, 3654, 406, 29, 1,
1, 30, 435, 4060, 27405, 142506, 593775, 2035800, 5852925, 14307150, 30045015, 54627300, 86493225, 119759850, 145422675, 155117520, 145422675, 119759850, 86493225, 54627300, 30045015, 14307150, 5852925, 2035800, 593775, 142506, 27405, 4060, 435, 30, 1,
1, 31, 465, 4495, 31465, 169911, 736281, 2629575, 7888725, 20160075, 44352165, 84672315, 141120525, 206253075, 265182525, 300540195, 300540195, 265182525, 206253075, 141120525, 84672315, 44352165, 20160075, 7888725, 2629575, 736281, 169911, 31465, 4495, 465, 31, 1,
1, 32, 496, 4960, 35960, 201376, 906192, 3365856, 10518300, 28048800, 64512240, 129024480, 225792840, 347373600, 471435600, 565722720, 601080390, 565722720, 471435600, 347373600, 225792840, 129024480, 64512240, 28048800, 10518300, 3365856, 906192, 201376, 35960, 4960, 496, 32, 1,
1, 33, 528, 5456, 40920, 237336, 1107568, 4272048, 13884156, 38567100, 92561040, 193536720, 354817320, 573166440, 818809200, 1037158320, 1166803110, 1166803110, 1037158320, 818809200, 573166440, 354817320, 193536720, 92561040, 38567100, 13884156, 4272048, 1107568, 237336, 40920, 5456, 528, 33, 1,
1, 34, 561, 5984, 46376, 278256, 1344904, 5379616, 18156204, 52451256, 131128140, 286097760, 548354040, 927983760, 1391975640, 1855967520, 2203961430, 2333606220, 2203961430, 1855967520, 1391975640, 927983760, 548354040, 286097760, 131128140, 52451256, 18156204, 5379616, 1344904, 278256, 46376, 5984, 561, 34, 1,
1, 35, 595, 6545, 52360, 324632, 1623160, 6724520, 23535820, 70607460, 183579396, 417225900, 834451800, 1476337800, 2319959400, 3247943160, 4059928950, 4537567650, 4537567650, 4059928950, 3247943160, 2319959400, 1476337800, 834451800, 417225900, 183579396, 70607460, 23535820, 6724520, 1623160, 324632, 52360, 6545, 595, 35, 1,
1, 36, 630, 7140, 58905, 376992, 1947792, 8347680, 30260340, 94143280, 254186856, 600805296, 1251677700, 2310789600, 3796297200, 5567902560, 7307872110, 8597496600, 9075135300, 8597496600, 7307872110, 5567902560, 3796297200, 2310789600, 1251677700, 600805296, 254186856, 94143280, 30260340, 8347680, 1947792, 376992, 58905, 7140, 630, 36, 1,
1, 37, 666, 7770, 66045, 435897, 2324784, 10295472, 38608020, 124403620, 348330136, 854992152, 1852482996, 3562467300, 6107086800, 9364199760, 12875774670, 15905368710, 17672631900, 17672631900, 15905368710, 12875774670, 9364199760, 6107086800, 3562467300, 1852482996, 854992152, 348330136, 124403620, 38608020, 10295472, 2324784, 435897, 66045, 7770, 666, 37, 1,
1, 38, 703, 8436, 73815, 501942, 2760681, 12620256, 48903492, 163011640, 472733756, 1203322288, 2707475148, 5414950296, 9669554100, 15471286560, 22239974430, 28781143380, 33578000610, 35345263800, 33578000610, 28781143380, 22239974430, 15471286560, 9669554100, 5414950296, 2707475148, 1203322288, 472733756, 163011640, 48903492, 12620256, 2760681, 501942, 73815, 8436, 703, 38, 1,
1, 39, 741, 9139, 82251, 575757, 3262623, 15380937, 61523748, 211915132, 635745396, 1676056044, 3910797436, 8122425444, 15084504396, 25140840660, 37711260990, 51021117810, 62359143990, 68923264410, 68923264410, 62359143990, 51021117810, 37711260990, 25140840660, 15084504396, 8122425444, 3910797436, 1676056044, 635745396, 211915132, 61523748, 15380937, 3262623, 575757, 82251, 9139, 741, 39, 1,
1, 40, 780, 9880, 91390, 658008, 3838380, 18643560, 76904685, 273438880, 847660528, 2311801440, 5586853480, 12033222880, 23206929840, 40225345056, 62852101650, 88732378800, 113380261800, 131282408400, 137846528820, 131282408400, 113380261800, 88732378800, 62852101650, 40225345056, 23206929840, 12033222880, 5586853480, 2311801440, 847660528, 273438880, 76904685, 18643560, 3838380, 658008, 91390, 9880, 780, 40, 1,
1, 41, 820, 10660, 101270, 749398, 4496388, 22481940, 95548245, 350343565, 1121099408, 3159461968, 7898654920, 17620076360, 35240152720, 63432274896, 103077446706, 151584480450, 202112640600, 244662670200, 269128937220, 269128937220, 244662670200, 202112640600, 151584480450, 103077446706, 63432274896, 35240152720, 17620076360, 7898654920, 3159461968, 1121099408, 350343565, 95548245, 22481940, 4496388, 749398, 101270, 10660, 820, 41, 1,
1, 42, 861, 11480, 111930, 850668, 5245786, 26978328, 118030185, 445891810, 1471442973, 4280561376, 11058116888, 25518731280, 52860229080, 98672427616, 166509721602, 254661927156, 353697121050, 446775310800, 513791607420, 538257874440, 513791607420, 446775310800, 353697121050, 254661927156, 166509721602, 98672427616, 52860229080, 25518731280, 11058116888, 4280561376, 1471442973, 445891810, 118030185, 26978328, 5245786, 850668, 111930, 11480, 861, 42, 1,
1, 43, 903, 12341, 123410, 962598, 6096454, 32224114, 145008513, 563921995, 1917334783, 5752004349, 15338678264, 36576848168, 78378960360, 151532656696, 265182149218, 421171648758, 608359048206, 800472431850, 960566918220, 1052049481860, 1052049481860, 960566918220, 800472431850, 608359048206, 421171648758, 265182149218, 151532656696, 78378960360, 36576848168, 15338678264, 5752004349, 1917334783, 563921995, 145008513, 32224114, 6096454, 962598, 123410, 12341, 903, 43, 1,
1, 44, 946, 13244, 135751, 1086008, 7059052, 38320568, 177232627, 708930508, 2481256778, 7669339132, 21090682613, 51915526432, 114955808528, 229911617056, 416714805914, 686353797976, 1029530696964, 1408831480056, 1761039350070, 2012616400080, 2104098963720, 2012616400080, 1761039350070, 1408831480056, 1029530696964, 686353797976, 416714805914, 229911617056, 114955808528, 51915526432, 21090682613, 7669339132, 2481256778, 708930508, 177232627, 38320568, 7059052, 1086008, 135751, 13244, 946, 44, 1,
1, 45, 990, 14190, 148995, 1221759, 8145060, 45379620, 215553195, 886163135, 3190187286, 10150595910, 28760021745, 73006209045, 166871334960, 344867425584, 646626422970, 1103068603890, 1715884494940, 2438362177020, 3169870830126, 3773655750150, 4116715363800, 4116715363800, 3773655750150, 3169870830126, 2438362177020, 1715884494940, 1103068603890, 646626422970, 344867425584, 166871334960, 73006209045, 28760021745, 10150595910, 3190187286, 886163135, 215553195, 45379620, 8145060, 1221759, 148995, 14190, 990, 45, 1,
1, 46, 1035, 15180, 163185, 1370754, 9366819, 53524680, 260932815, 1101716330, 4076350421, 13340783196, 38910617655, 101766230790, 239877544005, 511738760544, 991493848554, 1749695026860, 2818953098830, 4154246671960, 5608233007146, 6943526580276, 7890371113950, 8233430727600, 7890371113950, 6943526580276, 5608233007146, 4154246671960, 2818953098830, 1749695026860, 991493848554, 511738760544, 239877544005, 101766230790, 38910617655, 13340783196, 4076350421, 1101716330, 260932815, 53524680, 9366819, 1370754, 163185, 15180, 1035, 46, 1,
1, 47, 1081, 16215, 178365, 1533939, 10737573, 62891499, 314457495, 1362649145, 5178066751, 17417133617, 52251400851, 140676848445, 341643774795, 751616304549, 1503232609098, 2741188875414, 4568648125690, 6973199770790, 9762479679106, 12551759587422, 14833897694226, 16123801841550, 16123801841550, 14833897694226, 12551759587422, 9762479679106, 6973199770790, 4568648125690, 2741188875414, 1503232609098, 751616304549, 341643774795, 140676848445, 52251400851, 17417133617, 5178066751, 1362649145, 314457495, 62891499, 10737573, 1533939, 178365, 16215, 1081, 47, 1,
1, 48, 1128, 17296, 194580, 1712304, 12271512, 73629072, 377348994, 1677106640, 6540715896, 22595200368, 69668534468, 192928249296, 482320623240, 1093260079344, 2254848913647, 4244421484512, 7309837001104, 11541847896480, 16735679449896, 22314239266528, 27385657281648, 30957699535776, 32247603683100, 30957699535776, 27385657281648, 22314239266528, 16735679449896, 11541847896480, 7309837001104, 4244421484512, 2254848913647, 1093260079344, 482320623240, 192928249296, 69668534468, 22595200368, 6540715896, 1677106640, 377348994, 73629072, 12271512, 1712304, 194580, 17296, 1128, 48, 1,
1, 49, 1176, 18424, 211876, 1906884, 13983816, 85900584, 450978066, 2054455634, 8217822536, 29135916264, 92263734836, 262596783764, 675248872536, 1575580702584, 3348108992991, 6499270398159, 11554258485616, 18851684897584, 28277527346376, 39049918716424, 49699896548176, 58343356817424, 63205303218876, 63205303218876, 58343356817424, 49699896548176, 39049918716424, 28277527346376, 18851684897584, 11554258485616, 6499270398159, 3348108992991, 1575580702584, 675248872536, 262596783764, 92263734836, 29135916264, 8217822536, 2054455634, 450978066, 85900584, 13983816, 1906884, 211876, 18424, 1176, 49, 1,
1, 50, 1225, 19600, 230300, 2118760, 15890700, 99884400, 536878650, 2505433700, 10272278170, 37353738800, 121399651100, 354860518600, 937845656300, 2250829575120, 4923689695575, 9847379391150, 18053528883775, 30405943383200, 47129212243960, 67327446062800, 88749815264600, 108043253365600, 121548660036300, 126410606437752, 121548660036300, 108043253365600, 88749815264600, 67327446062800, 47129212243960, 30405943383200, 18053528883775, 9847379391150, 4923689695575, 2250829575120, 937845656300, 354860518600, 121399651100, 37353738800, 10272278170, 2505433700, 536878650, 99884400, 15890700, 2118760, 230300, 19600, 1225, 50, 1
};
/// The array containing the offsets for each row in the array binomialcoeffs_data.
constexpr size_t binomialcoeffs_offsets[] = { 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171, 190, 210, 231, 253, 276, 300, 325, 351, 378, 406, 435, 465, 496, 528, 561, 595, 630, 666, 703, 741, 780, 820, 861, 903, 946, 990, 1035, 1081, 1128, 1176, 1225, 1275, 1326 };
/// The maximum order for binomial coefficients supported currently (50).
constexpr size_t binomialcoeffs_nmax = 50;
/// Return the binomial coefficient C(i,j) at compile time.
template<size_t i, size_t j>
struct AuxBinomialCoefficient
{
static_assert(i <= binomialcoeffs_nmax, "Violation of maximum order for binomial coefficient retrieval.");
static_assert(j <= i, "Violation of j <= i condition for retrieving binomial coefficient C(i,j).");
constexpr static double value = binomialcoeffs_data[binomialcoeffs_offsets[i] + j];
};
/// The binomial coefficient C(i,j) at compile time.
template<size_t i, size_t j>
constexpr double BinomialCoefficient = AuxBinomialCoefficient<i, j>::value;
} // namespace detail
} // namespace autodiff

View File

@ -0,0 +1,106 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// autodiff includes
#include <autodiff/common/meta.hpp>
namespace autodiff {
namespace detail {
//==========================================================================================================================================================
// The code below was taken from: https://stackoverflow.com/questions/87372/check-if-a-class-has-a-member-function-of-a-given-signature/16867422#16867422
// It implements checkers to determine if a type has a member variable, function, etc.
//==========================================================================================================================================================
template <typename... Args> struct ambiguate : public Args... {};
template<typename A, typename = void>
struct got_type : std::false_type {};
template<typename A>
struct got_type<A> : std::true_type { typedef A type; };
template<typename T, T>
struct sig_check : std::true_type {};
template<typename Alias, typename AmbiguitySeed>
struct has_member {
template<typename C> static char ((&f(decltype(&C::value))))[1];
template<typename C> static char ((&f(...)))[2];
//Make sure the member name is consistently spelled the same.
static_assert(
(sizeof(f<AmbiguitySeed>(0)) == 1)
, "Member name specified in AmbiguitySeed is different from member name specified in Alias, or wrong Alias/AmbiguitySeed has been specified."
);
static bool const value = sizeof(f<Alias>(0)) == 2;
};
//Check for any member with given name, whether var, func, class, union, enum.
#define CREATE_MEMBER_CHECK(member) \
\
template<typename T, typename = std::true_type> \
struct Alias_##member; \
\
template<typename T> \
struct Alias_##member < \
T, std::integral_constant<bool, got_type<decltype(&T::member)>::value> \
> { static const decltype(&T::member) value; }; \
\
struct AmbiguitySeed_##member { char member; }; \
\
template<typename T> \
struct has_member_##member { \
static const bool value \
= has_member< \
Alias_##member<ambiguate<T, AmbiguitySeed_##member>> \
, Alias_##member<AmbiguitySeed_##member> \
>::value \
; \
}
// Create type trait struct `has_member_size`.
CREATE_MEMBER_CHECK(size);
/// Boolean constant that is true if type T implements `size` method.
template<typename T>
constexpr bool hasSize = has_member_size<PlainType<T>>::value;
// Create type trait struct `has_operator_bracket`.
template<class, typename T> struct has_operator_bracket_impl : std::false_type {};
template<typename T> struct has_operator_bracket_impl<decltype( void(std::declval<T>().operator [](0)) ), T> : std::true_type {};
/// Boolean type that is true if type T implements `operator[](int)` method.
template<typename T> struct has_operator_bracket : has_operator_bracket_impl<void, T> {};
} // namespace detail
} // namespace autodiff

View File

@ -0,0 +1,155 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// Eigen includes
#include <Eigen/Core>
// autodiff includes
#include <autodiff/common/vectortraits.hpp>
//=====================================================================================================================
//
// EIGEN MACROS FOR CREATING NEW TYPE ALIASES
//
//=====================================================================================================================
#define AUTODIFF_DEFINE_EIGEN_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix) \
using Array##SizeSuffix##SizeSuffix##TypeSuffix = Eigen::Array<Type, Size, Size>; \
using Array##SizeSuffix##TypeSuffix = Eigen::Array<Type, Size, 1>; \
using Matrix##SizeSuffix##TypeSuffix = Eigen::Matrix<Type, Size, Size, 0, Size, Size>; \
using Vector##SizeSuffix##TypeSuffix = Eigen::Matrix<Type, Size, 1, 0, Size, 1>; \
using RowVector##SizeSuffix##TypeSuffix = Eigen::Matrix<Type, 1, Size, 1, 1, Size>;
#define AUTODIFF_DEFINE_EIGEN_FIXED_TYPEDEFS(Type, TypeSuffix, Size) \
using Array##Size##X##TypeSuffix = Eigen::Array<Type, Size, -1>; \
using Array##X##Size##TypeSuffix = Eigen::Array<Type, -1, Size>; \
using Matrix##Size##X##TypeSuffix = Eigen::Matrix<Type, Size, -1, 0, Size, -1>; \
using Matrix##X##Size##TypeSuffix = Eigen::Matrix<Type, -1, Size, 0, -1, Size>;
#define AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(Type, TypeSuffix) \
AUTODIFF_DEFINE_EIGEN_TYPEDEFS(Type, TypeSuffix, 2, 2) \
AUTODIFF_DEFINE_EIGEN_TYPEDEFS(Type, TypeSuffix, 3, 3) \
AUTODIFF_DEFINE_EIGEN_TYPEDEFS(Type, TypeSuffix, 4, 4) \
AUTODIFF_DEFINE_EIGEN_TYPEDEFS(Type, TypeSuffix, -1, X) \
AUTODIFF_DEFINE_EIGEN_FIXED_TYPEDEFS(Type, TypeSuffix, 2) \
AUTODIFF_DEFINE_EIGEN_FIXED_TYPEDEFS(Type, TypeSuffix, 3) \
AUTODIFF_DEFINE_EIGEN_FIXED_TYPEDEFS(Type, TypeSuffix, 4)
namespace autodiff {
namespace detail {
//=====================================================================================================================
//
// DEFINE VECTOR TRAITS FOR EIGEN TYPES
//
//=====================================================================================================================
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct VectorTraits<Eigen::Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
using ValueType = Scalar;
template<typename NewValueType>
using ReplaceValueType = Eigen::Matrix<NewValueType, Rows, Cols, Options, MaxRows, MaxCols>;
};
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct VectorTraits<Eigen::Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols>>
{
using ValueType = Scalar;
template<typename NewValueType>
using ReplaceValueType = Eigen::Array<NewValueType, Rows, Cols, Options, MaxRows, MaxCols>;
};
template<typename VectorType, int Size>
struct VectorTraits<Eigen::VectorBlock<VectorType, Size>>
{
using ValueType = typename PlainType<VectorType>::Scalar;
template<typename NewValueType>
using ReplaceValueType = VectorReplaceValueType<VectorType, NewValueType>;
};
#if EIGEN_VERSION_AT_LEAST(3, 3, 90)
template<typename VectorType, typename IndicesType>
struct VectorTraits<Eigen::IndexedView<VectorType, IndicesType, Eigen::internal::SingleRange>>
{
using ValueType = typename PlainType<VectorType>::Scalar;
template<typename NewValueType>
using ReplaceValueType = VectorReplaceValueType<VectorType, NewValueType>;
};
template<typename VectorType, typename IndicesType>
struct VectorTraits<Eigen::IndexedView<VectorType, Eigen::internal::SingleRange, IndicesType>>
{
using ValueType = typename PlainType<VectorType>::Scalar;
template<typename NewValueType>
using ReplaceValueType = VectorReplaceValueType<VectorType, NewValueType>;
};
#endif
template<typename MatrixType>
struct VectorTraits<Eigen::Ref<MatrixType>>
{
using ValueType = VectorValueType<MatrixType>;
template<typename NewValueType>
using ReplaceValueType = VectorReplaceValueType<MatrixType, NewValueType>;
};
template<typename VectorType, int MapOptions, typename StrideType>
struct VectorTraits<Eigen::Map<VectorType, MapOptions, StrideType>>
{
using ValueType = VectorValueType<VectorType>;
template<typename NewValueType>
using ReplaceValueType = Eigen::Map<VectorReplaceValueType<VectorType, NewValueType>, MapOptions, StrideType>;
};
//=====================================================================================================================
//
// AUXILIARY TEMPLATE TYPE ALIASES
//
//=====================================================================================================================
template<typename Scalar>
using VectorX = Eigen::Matrix<Scalar, Eigen::Dynamic, 1, 0, Eigen::Dynamic, 1>;
template<typename Scalar>
using MatrixX = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, 0, Eigen::Dynamic, Eigen::Dynamic>;
} // namespace detail
} // namespace autodiff

View File

@ -0,0 +1,186 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// C++ includes
#include <cstddef>
#include <tuple>
#include <type_traits>
namespace autodiff {
namespace detail {
template<bool value>
using EnableIf = std::enable_if_t<value>;
template<bool value>
using Requires = std::enable_if_t<value, bool>;
template<typename T>
using PlainType = std::remove_cv_t<std::remove_reference_t<T>>;
template<bool Cond, typename WhenTrue, typename WhenFalse>
using ConditionalType = std::conditional_t<Cond, WhenTrue, WhenFalse>;
template<typename A, typename B>
using CommonType = std::common_type_t<A, B>;
template<typename Fun, typename... Args>
using ReturnType = std::invoke_result_t<Fun, Args...>;
template<typename T>
constexpr bool isConst = std::is_const_v<std::remove_reference_t<T>>;
template<typename T, typename U>
constexpr bool isConvertible = std::is_convertible<PlainType<T>, U>::value;
template<typename A, typename B>
constexpr bool isSame = std::is_same_v<A, B>;
template<typename Tuple>
constexpr auto TupleSize = std::tuple_size_v<std::decay_t<Tuple>>;
template<typename Tuple>
constexpr auto TupleHead(Tuple&& tuple)
{
return std::get<0>(std::forward<Tuple>(tuple));
}
template<typename Tuple>
constexpr auto TupleTail(Tuple&& tuple)
{
auto g = [&](auto&&, auto&&... args) constexpr {
return std::forward_as_tuple(args...);
};
return std::apply(g, std::forward<Tuple>(tuple));
}
template<size_t i>
struct Index
{
constexpr static size_t index = i;
constexpr operator size_t() const { return index; }
constexpr operator size_t() { return index; }
};
template<size_t i, size_t ibegin, size_t iend, typename Function>
constexpr auto AuxFor(Function&& f)
{
if constexpr (i < iend) {
f(Index<i>{});
AuxFor<i + 1, ibegin, iend>(std::forward<Function>(f));
}
}
template<size_t ibegin, size_t iend, typename Function>
constexpr auto For(Function&& f)
{
AuxFor<ibegin, ibegin, iend>(std::forward<Function>(f));
}
template<size_t iend, typename Function>
constexpr auto For(Function&& f)
{
For<0, iend>(std::forward<Function>(f));
}
template<size_t i, size_t ibegin, size_t iend, typename Function>
constexpr auto AuxReverseFor(Function&& f)
{
if constexpr (i < iend)
{
AuxReverseFor<i + 1, ibegin, iend>(std::forward<Function>(f));
f(Index<i>{});
}
}
template<size_t ibegin, size_t iend, typename Function>
constexpr auto ReverseFor(Function&& f)
{
AuxReverseFor<ibegin, ibegin, iend>(std::forward<Function>(f));
}
template<size_t iend, typename Function>
constexpr auto ReverseFor(Function&& f)
{
ReverseFor<0, iend>(std::forward<Function>(f));
}
template<typename Tuple, typename Function>
constexpr auto ForEach(Tuple&& tuple, Function&& f)
{
constexpr auto N = TupleSize<Tuple>;
For<N>([&](auto i) constexpr {
f(std::get<i>(tuple));
});
//------------------------------------------------------------
// ALTERNATIVE IMPLEMENTATION POSSIBLY USEFUL TO KEEP IT HERE
// auto g = [&](auto&&... args) constexpr {
// ( f(std::forward<decltype(args)>(args)), ...);
// };
// std::apply(g, std::forward<Tuple>(tuple));
//------------------------------------------------------------
}
template<typename Tuple1, typename Tuple2, typename Function>
constexpr auto ForEach(Tuple1&& tuple1, Tuple2&& tuple2, Function&& f)
{
constexpr auto N1 = TupleSize<Tuple1>;
constexpr auto N2 = TupleSize<Tuple2>;
static_assert(N1 == N2);
For<N1>([&](auto i) constexpr {
f(std::get<i>(tuple1), std::get<i>(tuple2));
});
}
template<size_t ibegin, size_t iend, typename Function>
constexpr auto Sum(Function&& f)
{
using ResultType = std::invoke_result_t<Function, Index<ibegin>>;
ResultType res = {};
For<ibegin, iend>([&](auto i) constexpr {
res += f(Index<i>{});
});
return res;
}
template<typename Tuple, typename Function>
constexpr auto Reduce(Tuple&& tuple, Function&& f)
{
using ResultType = std::invoke_result_t<Function, decltype(std::get<0>(tuple))>;
ResultType res = {};
ForEach(tuple, [&](auto&& item) constexpr {
res += f(item);
});
return res;
}
} // namespace detail
} // namespace autodiff

View File

@ -0,0 +1,73 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// C++ includes
#include <autodiff/common/meta.hpp>
namespace autodiff {
namespace detail {
/// A trait class used to specify whether a type is arithmetic.
template<typename T>
struct ArithmeticTraits
{
static constexpr bool isArithmetic = std::is_arithmetic_v<T>;
};
/// A compile-time constant that indicates whether a type is arithmetic.
template<typename T>
constexpr bool isArithmetic = ArithmeticTraits<PlainType<T>>::isArithmetic;
/// An auxiliary template type to indicate NumberTraits has not been defined for a type.
template<typename T>
struct NumericTypeInfoNotDefinedFor { using type = T; };
/// A trait class used to specify whether a type is an autodiff number.
template<typename T>
struct NumberTraits
{
/// The underlying floating point type of the autodiff number type.
using NumericType = std::conditional_t<isArithmetic<T>, T, NumericTypeInfoNotDefinedFor<T>>;
/// The order of the autodiff number type.
static constexpr auto Order = 0;
};
/// A template alias to get the underlying floating point type of an autodiff number.
template<typename T>
using NumericType = typename NumberTraits<PlainType<T>>::NumericType;
/// A compile-time constant with the order of an autodiff number.
template<typename T>
constexpr auto Order = NumberTraits<PlainType<T>>::Order;
} // namespace detail
} // namespace autodiff

View File

@ -0,0 +1,84 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// C++ includes
#include <vector>
// autodiff includes
#include <autodiff/common/meta.hpp>
namespace autodiff {
namespace detail {
/// An auxiliary template type to indicate VectorTraits has not been defined for a type.
template<typename V>
struct VectorTraitsNotDefinedFor {};
/// An auxiliary template type to indicate VectorTraits::ReplaceValueType is not supported for a type.
template<typename V>
struct VectorReplaceValueTypeNotSupportedFor {};
/// A vector traits to be defined for each autodiff number.
template<typename V, class Enable = void>
struct VectorTraits
{
/// The value type of each entry in the vector.
using ValueType = VectorTraitsNotDefinedFor<V>;
/// The template alias to replace the value type of a vector type with another value type.
using ReplaceValueType = VectorReplaceValueTypeNotSupportedFor<V>;
};
/// A template alias used to get the type of the values in a vector type.
template<typename V>
using VectorValueType = typename VectorTraits<PlainType<V>>::ValueType;
/// A template alias used to get the type of a vector that is equivalent to another but with a different value type.
template<typename V, typename NewValueType>
using VectorReplaceValueType = typename VectorTraits<PlainType<V>>::template ReplaceValueType<NewValueType>;
/// A compile-time constant that indicates with a type is a vector type.
template<typename V>
constexpr bool isVector = !std::is_same_v<VectorValueType<PlainType<V>>, VectorTraitsNotDefinedFor<PlainType<V>>>;
/// Implementation of VectorTraits for std::vector.
template<typename T, template<class> typename Allocator>
struct VectorTraits<std::vector<T, Allocator<T>>>
{
using ValueType = T;
template<typename NewValueType>
using ReplaceValueType = std::vector<NewValueType, Allocator<NewValueType>>;
};
} // namespace detail
} // namespace autodiff

View File

@ -0,0 +1,35 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// autodiff includes
#include <autodiff/forward/dual/dual.hpp>
#include <autodiff/forward/utils/derivative.hpp>
#include <autodiff/forward/utils/taylorseries.hpp>

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,129 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// Eigen includes
#include <Eigen/Core>
// autodiff includes
#include <autodiff/forward/dual.hpp>
#include <autodiff/forward/utils/gradient.hpp>
#include <autodiff/common/eigen.hpp>
//------------------------------------------------------------------------------
// SUPPORT FOR EIGEN MATRICES AND VECTORS OF DUAL
//------------------------------------------------------------------------------
namespace Eigen {
using namespace autodiff;
using namespace autodiff::detail;
template<typename T>
struct NumTraits;
template<typename T, typename G>
struct NumTraits<autodiff::Dual<T, G>> : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
typedef autodiff::Dual<T, G> Real;
typedef autodiff::Dual<T, G> NonInteger;
typedef autodiff::Dual<T, G> Nested;
enum
{
IsComplex = 0,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1,
ReadCost = 1,
AddCost = 3,
MulCost = 3
};
};
template<typename T, typename G, typename BinOp>
struct ScalarBinaryOpTraits<Dual<T, G>, NumericType<T>, BinOp>
{
typedef DualType<Dual<T, G>> ReturnType;
};
template<typename Op, typename R, typename BinOp>
struct ScalarBinaryOpTraits<UnaryExpr<Op, R>, NumericType<UnaryExpr<Op, R>>, BinOp>
{
typedef DualType<UnaryExpr<Op, R>> ReturnType;
};
template<typename Op, typename L, typename R, typename BinOp>
struct ScalarBinaryOpTraits<BinaryExpr<Op, L, R>, NumericType<BinaryExpr<Op, L, R>>, BinOp>
{
typedef DualType<BinaryExpr<Op, L, R>> ReturnType;
};
template<typename Op, typename L, typename C, typename R, typename BinOp>
struct ScalarBinaryOpTraits<TernaryExpr<Op, L, C, R>, NumericType<TernaryExpr<Op, L, C, R>>, BinOp>
{
typedef DualType<TernaryExpr<Op, L, C, R>> ReturnType;
};
template<typename T, typename G, typename BinOp>
struct ScalarBinaryOpTraits<NumericType<T>, Dual<T, G>, BinOp>
{
typedef DualType<Dual<T, G>> ReturnType;
};
template<typename Op, typename R, typename BinOp>
struct ScalarBinaryOpTraits<NumericType<UnaryExpr<Op, R>>, UnaryExpr<Op, R>, BinOp>
{
typedef DualType<UnaryExpr<Op, R>> ReturnType;
};
template<typename Op, typename L, typename R, typename BinOp>
struct ScalarBinaryOpTraits<NumericType<BinaryExpr<Op, L, R>>, BinaryExpr<Op, L, R>, BinOp>
{
typedef DualType<BinaryExpr<Op, L, R>> ReturnType;
};
template<typename Op, typename L, typename C, typename R, typename BinOp>
struct ScalarBinaryOpTraits<NumericType<TernaryExpr<Op, L, C, R>>, TernaryExpr<Op, L, C, R>, BinOp>
{
typedef DualType<TernaryExpr<Op, L, C, R>> ReturnType;
};
} // namespace Eigen
namespace autodiff {
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(dual0th, dual0th);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(dual1st, dual1st);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(dual2nd, dual2nd);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(dual3rd, dual3rd);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(dual4th, dual4th);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(dual, dual)
} // namespace autodiff

View File

@ -0,0 +1,35 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// autodiff includes
#include <autodiff/forward/real/real.hpp>
#include <autodiff/forward/utils/derivative.hpp>
#include <autodiff/forward/utils/taylorseries.hpp>

View File

@ -0,0 +1,90 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// Eigen includes
#include <Eigen/Core>
// autodiff includes
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/utils/gradient.hpp>
#include <autodiff/common/eigen.hpp>
//------------------------------------------------------------------------------
// SUPPORT FOR EIGEN MATRICES AND VECTORS OF REAL
//------------------------------------------------------------------------------
namespace Eigen {
template<typename T>
struct NumTraits;
template<size_t N, typename T>
struct NumTraits<autodiff::Real<N, T>> : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
typedef autodiff::Real<N, T> Real;
typedef autodiff::Real<N, T> NonInteger;
typedef autodiff::Real<N, T> Nested;
enum
{
IsComplex = 0,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1,
ReadCost = 1,
AddCost = 3,
MulCost = 3
};
};
template<size_t N, typename T, typename BinOp>
struct ScalarBinaryOpTraits<autodiff::Real<N, T>, T, BinOp>
{
typedef autodiff::Real<N, T> ReturnType;
};
template<size_t N, typename T, typename BinOp>
struct ScalarBinaryOpTraits<T, autodiff::Real<N, T>, BinOp>
{
typedef autodiff::Real<N, T> ReturnType;
};
} // namespace Eigen
namespace autodiff {
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(real0th, real0th);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(real1st, real1st);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(real2nd, real2nd);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(real3rd, real3rd);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(real4th, real4th);
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(real, real)
} // namespace autodiff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,268 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
// C++ includes
#include <cassert>
#include <cstddef>
// autodiff includes
#include <autodiff/common/meta.hpp>
#include <autodiff/common/vectortraits.hpp>
#pragma once
namespace autodiff {
namespace detail {
template<typename... Args>
struct At
{
std::tuple<Args...> args;
};
template<typename... Args>
struct Wrt
{
std::tuple<Args...> args;
};
template<typename... Args>
struct Along
{
std::tuple<Args...> args;
};
/// The keyword used to denote the variables *with respect to* the derivative is calculated.
template<typename... Args>
auto wrt(Args&&... args)
{
return Wrt<Args&&...>{ std::forward_as_tuple(std::forward<Args>(args)...) };
}
/// The keyword used to denote the variables *at* which the derivatives are calculated.
template<typename... Args>
auto at(Args&&... args)
{
return At<Args&&...>{ std::forward_as_tuple(std::forward<Args>(args)...) };
}
/// The keyword used to denote the direction vector *along* which the derivatives are calculated.
template<typename... Args>
auto along(Args&&... args)
{
return Along<Args&&...>{ std::forward_as_tuple(std::forward<Args>(args)...) };
}
/// Seed each dual number in the **wrt** list using its position as the derivative order to be seeded.
/// Using `seed(wrt(x, y, z), 1)` will set the 1st order derivative of `x`, the
/// 2nd order derivative of `y`, and the 3rd order derivative of `z` to 1. If
/// these dual numbers have order greater than 3, then the last dual number will
/// be used for the remaining higher-order derivatives. For example, if these
/// numbers are 5th order, than the 4th and 5th order derivatives of `z` will be
/// set to 1 as well. In this example, `wrt(x, y, z)` is equivalent to `wrt(x,
/// y, z, z, z)`. This automatic seeding permits derivatives `fx`, `fxy`,
/// `fxyz`, `fxyzz`, and `fxyzzz` to be computed in a more convenient way.
template<typename Var, typename... Vars, typename T>
auto seed(const Wrt<Var&, Vars&...>& wrt, T&& seedval)
{
constexpr static auto N = Order<Var>;
constexpr static auto size = 1 + sizeof...(Vars);
static_assert(size <= N, "It is not possible to compute higher-order derivatives with order greater than that of the autodiff number (e.g., using wrt(x, x, y, z) will fail if the autodiff numbers in use have order below 4).");
For<N>([&](auto i) constexpr {
if constexpr (i.index < size)
seed<i.index + 1>(std::get<i>(wrt.args), seedval);
else
seed<i.index + 1>(std::get<size - 1>(wrt.args), seedval); // use the last variable in the wrt list as the variable for which the remaining higher-order derivatives are calculated (e.g., derivatives(f, wrt(x), at(x)) will produce [f0, fx, fxx, fxxx, fxxxx] when x is a 4th order dual number).
});
}
template<typename... Vars>
auto seed(const Wrt<Vars...>& wrt)
{
seed(wrt, 1.0);
}
template<typename... Vars>
auto unseed(const Wrt<Vars...>& wrt)
{
seed(wrt, 0.0);
}
template<typename... Args, typename... Vecs>
auto seed(const At<Args...>& at, const Along<Vecs...>& along)
{
static_assert(sizeof...(Args) == sizeof...(Vecs));
ForEach(at.args, along.args, [&](auto& arg, auto&& dir) constexpr {
if constexpr (isVector<decltype(arg)>) {
static_assert(isVector<decltype(dir)>);
assert(arg.size() == dir.size());
for(auto i = 0; i < dir.size(); ++i)
seed<1>(arg[i], dir[i]);
}
else seed<1>(arg, dir);
});
}
template<typename... Args>
auto unseed(const At<Args...>& at)
{
ForEach(at.args, [&](auto& arg) constexpr {
if constexpr (isVector<decltype(arg)>) {
for(auto i = 0; i < arg.size(); ++i)
seed<1>(arg[i], 0.0);
}
else seed<1>(arg, 0.0);
});
}
template<size_t order = 1, typename T, Requires<Order<T>> = true>
auto seed(T& x)
{
seed<order>(x, 1.0);
}
template<size_t order = 1, typename T, Requires<Order<T>> = true>
auto unseed(T& x)
{
seed<order>(x, 0.0);
}
template<typename Fun, typename... Args, typename... Vars>
auto eval(const Fun& f, const At<Args...>& at, const Wrt<Vars...>& wrt)
{
seed(wrt);
auto u = std::apply(f, at.args);
unseed(wrt);
return u;
}
template<typename Fun, typename... Args, typename... Vecs>
auto eval(const Fun& f, const At<Args...>& at, const Along<Vecs...>& along)
{
seed(at, along);
auto u = std::apply(f, at.args);
unseed(at);
return u;
}
/// Extract the derivative of given order from a vector of dual/real numbers.
template<size_t order = 1, typename Vec, Requires<isVector<Vec>> = true>
auto derivative(const Vec& u)
{
size_t len = u.size(); // the length of the vector containing dual/real numbers
using NumType = decltype(u[0]); // get the type of the dual/real number
using T = NumericType<NumType>; // get the numeric/floating point type of the dual/real number
using Res = VectorReplaceValueType<Vec, T>; // get the type of the vector containing numeric values instead of dual/real numbers (e.g., vector<real> becomes vector<double>, VectorXdual becomes VectorXd, etc.)
Res res(len); // create an array to store the derivatives stored inside the dual/real number
for(auto i = 0; i < len; ++i)
res[i] = derivative<order>(u[i]); // get the derivative of given order from i-th dual/real number
return res;
}
/// Alias method to `derivative<order>(x)` where `x` is either a dual/real number or vector/array of such numbers.
template<size_t order = 1, typename T>
auto grad(const T& x)
{
return derivative<order>(x);
}
/// Unpack the derivatives from the result of an @ref eval call into an array.
template<typename Result>
auto derivatives(const Result& result)
{
if constexpr (isVector<Result>) // check if the argument is a vector container of dual/real numbers
{
size_t len = result.size(); // the length of the vector containing dual/real numbers
using NumType = decltype(result[0]); // get the type of the dual/real number
using T = NumericType<NumType>; // get the numeric/floating point type of the dual/real number
using Vec = VectorReplaceValueType<Result, T>; // get the type of the vector containing numeric values instead of dual/real numbers (e.g., vector<real> becomes vector<double>, VectorXdual becomes VectorXd, etc.)
constexpr auto N = Order<NumType>; // the order of the dual/real number
std::array<Vec, N + 1> values; // create an array to store the derivatives stored inside the dual/real number
For<N + 1>([&](auto i) constexpr {
values[i].resize(len);
for(auto j = 0; j < len; ++j)
values[i][j] = derivative<i>(result[j]); // get the ith derivative of the jth dual/real number
});
return values;
}
else // result is then just a dual/real number
{
using T = NumericType<Result>; // get the numeric/floating point type of the dual/real result number
constexpr auto N = Order<Result>; // the order of the dual/real result number
std::array<T, N + 1> values; // create an array to store the derivatives stored inside the dual/real number
For<N + 1>([&](auto i) constexpr {
values[i] = derivative<i>(result);
});
return values;
}
}
template<typename Fun, typename... Vars, typename... Args>
auto derivatives(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at)
{
return derivatives(eval(f, at, wrt));
}
template<size_t order=1, typename Fun, typename... Vars, typename... Args, typename Result>
auto derivative(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, Result& u)
{
u = derivatives(f, wrt, at);
return derivative<order>(u);
}
template<size_t order=1, typename Fun, typename... Vars, typename... Args>
auto derivative(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at)
{
auto u = eval(f, at, wrt);
return derivative<order>(u);
}
template<typename Fun, typename... Vecs, typename... Args>
auto derivatives(const Fun& f, const Along<Vecs...>& along, const At<Args...>& at)
{
return derivatives(eval(f, at, along));
}
} // namespace detail
using detail::derivatives;
using detail::derivative;
using detail::grad;
using detail::along;
using detail::wrt;
using detail::at;
using detail::seed;
using detail::unseed;
using detail::Along;
using detail::At;
using detail::Wrt;
} // namespace autodiff

View File

@ -0,0 +1,228 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// autodiff includes
#include <autodiff/common/eigen.hpp>
#include <autodiff/common/meta.hpp>
#include <autodiff/common/classtraits.hpp>
#include <autodiff/forward/utils/derivative.hpp>
namespace autodiff {
namespace detail {
/// Return the length of an item in a `wrt(...)` list.
template<typename Item>
auto wrt_item_length(const Item& item) -> size_t
{
if constexpr (isVector<Item>)
return item.size(); // if item is a vector, return its size
else return 1; // if not a vector, say, a number, return 1 for its length
}
/// Return the sum of lengths of all itens in a `wrt(...)` list.
template<typename... Vars>
auto wrt_total_length(const Wrt<Vars...>& wrt) -> size_t
{
return Reduce(wrt.args, [&](auto&& item) constexpr {
return wrt_item_length(item);
});
}
// Loop through each variable in a wrt list and apply a function f(i, x) that
// accepts an index i and the variable x[i], where i is the global index of the
// variable in the list.
template<typename Function, typename... Vars>
constexpr auto ForEachWrtVar(const Wrt<Vars...>& wrt, Function&& f)
{
auto i = 0; // the current index of the variable in the wrt list
ForEach(wrt.args, [&](auto& item) constexpr
{
using T = decltype(item);
static_assert(isVector<T> || Order<T> > 0, "Expecting a wrt list with either vectors or individual autodiff numbers.");
if constexpr (isVector<T>) {
for(auto j = 0; j < item.size(); ++j)
// call given f with current index and variable from item (a vector)
if constexpr (detail::has_operator_bracket<T>()) {
f(i++, item[j]);
} else {
f(i++, item(j));
}
}
else f(i++, item); // call given f with current index and variable from item (a number, not a vector)
});
}
/// Return the gradient of scalar function *f* with respect to some or all variables *x*.
template<typename Fun, typename... Vars, typename... Args, typename Y, typename G>
void gradient(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, Y& u, G& g)
{
static_assert(sizeof...(Vars) >= 1);
static_assert(sizeof...(Args) >= 1);
const size_t n = wrt_total_length(wrt);
g.resize(n);
if(n == 0) return;
ForEachWrtVar(wrt, [&](auto&& i, auto&& xi) constexpr
{
static_assert(!isConst<decltype(xi)>, "Expecting non-const autodiff numbers in wrt list because these need to be seeded, and thus altered!");
u = eval(f, at, detail::wrt(xi)); // evaluate u with xi seeded so that du/dxi is also computed
g[i] = derivative<1>(u);
});
}
/// Return the gradient of scalar function *f* with respect to some or all variables *x*.
template<typename Fun, typename... Vars, typename... Args, typename Y>
auto gradient(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, Y& u)
{
using T = NumericType<decltype(u)>; // the underlying numeric floating point type in the autodiff number u
using Vec = VectorX<T>; // the gradient vector type with floating point values (not autodiff numbers!)
Vec g;
gradient(f, wrt, at, u, g);
return g;
}
/// Return the gradient of scalar function *f* with respect to some or all variables *x*.
template<typename Fun, typename... Vars, typename... Args>
auto gradient(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at)
{
ReturnType<Fun, Args...> u;
return gradient(f, wrt, at, u);
}
/// Return the Jacobian matrix of a function *f* with respect to some or all variables.
template<typename Fun, typename... Vars, typename... Args, typename Y, typename Jac>
void jacobian(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, Y& F, Jac& J)
{
static_assert(sizeof...(Vars) >= 1);
static_assert(sizeof...(Args) >= 1);
size_t n = wrt_total_length(wrt); /// using const size_t produces an error in GCC 7.3 because of the capture in the constexpr lambda in the ForEach block
size_t m = 0;
ForEachWrtVar(wrt, [&](auto&& i, auto&& xi) constexpr {
static_assert(!isConst<decltype(xi)>, "Expecting non-const autodiff numbers in wrt list because these need to be seeded, and thus altered!");
F = eval(f, at, detail::wrt(xi)); // evaluate F with xi seeded so that dF/dxi is also computed
if(m == 0) { m = F.size(); J.resize(m, n); };
for(size_t row = 0; row < m; ++row)
J(row, i) = derivative<1>(F[row]);
});
}
/// Return the Jacobian matrix of a function *f* with respect to some or all variables.
template<typename Fun, typename... Vars, typename... Args, typename Y>
auto jacobian(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, Y& F)
{
using U = VectorValueType<decltype(F)>; // the type of the autodiff numbers in vector F
using T = NumericType<U>; // the underlying numeric floating point type in the autodiff number U
using Mat = MatrixX<T>; // the jacobian matrix type with floating point values (not autodiff numbers!)
Mat J;
jacobian(f, wrt, at, F, J);
return J;
}
/// Return the Jacobian matrix of a function *f* with respect to some or all variables.
template<typename Fun, typename... Vars, typename... Args>
auto jacobian(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at)
{
using Y = ReturnType<Fun, Args...>;
static_assert(!std::is_same_v<Y, void>,
"In jacobian(f, wrt(x), at(x)), the type of x "
"might not be the same as in the definition of f. "
"For example, x is Eigen::VectorXdual but the "
"definition of f uses Eigen::Ref<const Eigen::VectorXdual>.");
Y F;
return jacobian(f, wrt, at, F);
}
/// Return the hessian matrix of scalar function *f* with respect to some or all variables *x*.
template<typename Fun, typename... Vars, typename... Args, typename U, typename G, typename H>
void hessian(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, U& u, G& g, H& h)
{
static_assert(sizeof...(Vars) >= 1);
static_assert(sizeof...(Args) >= 1);
size_t n = wrt_total_length(wrt);
g.resize(n);
h.resize(n, n);
ForEachWrtVar(wrt, [&](auto&& i, auto&& xi) constexpr {
ForEachWrtVar(wrt, [&](auto&& j, auto&& xj) constexpr
{
static_assert(!isConst<decltype(xi)> && !isConst<decltype(xj)>, "Expecting non-const autodiff numbers in wrt list because these need to be seeded, and thus altered!");
if(j >= i) { // this take advantage of the fact the Hessian matrix is symmetric
u = eval(f, at, detail::wrt(xi, xj)); // evaluate u with xi and xj seeded to produce u0, du/dxi, d2u/dxidxj
g[i] = derivative<1>(u); // get du/dxi from u
h(i, j) = h(j, i) = derivative<2>(u); // get d2u/dxidxj from u
}
});
});
}
/// Return the hessian matrix of scalar function *f* with respect to some or all variables *x*.
template<typename Fun, typename... Vars, typename... Args, typename U, typename G>
auto hessian(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at, U& u, G& g)
{
using T = NumericType<decltype(u)>; // the underlying numeric floating point type in the autodiff number u
using Mat = MatrixX<T>; // the Hessian matrix type with floating point values (not autodiff numbers!)
Mat H;
hessian(f, wrt, at, u, g, H);
return H;
}
/// Return the hessian matrix of scalar function *f* with respect to some or all variables *x*.
template<typename Fun, typename... Vars, typename... Args>
auto hessian(const Fun& f, const Wrt<Vars...>& wrt, const At<Args...>& at)
{
using U = ReturnType<Fun, Args...>;
using T = NumericType<U>;
using Vec = VectorX<T>;
U u;
Vec g;
return hessian(f, wrt, at, u, g);
}
} // namespace detail
using detail::gradient;
using detail::jacobian;
using detail::hessian;
} // namespace autodiff

View File

@ -0,0 +1,97 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// C++ includes
#include <array>
// autodiff includes
#include <autodiff/forward/utils/derivative.hpp>
#include <autodiff/common/meta.hpp>
#include <autodiff/common/vectortraits.hpp>
namespace autodiff {
namespace detail {
/// Represents a Taylor series along a direction for either a scalar or vector function.
/// @see taylorseries
template<size_t N, typename V>
class TaylorSeries
{
public:
/// The numeric floating point type of the derivatives, which can be a vector of values or just one.
using T = std::conditional_t<isVector<V>, VectorValueType<V>, V>;
/// Construct a default TaylorSeries object.
TaylorSeries() = default;
/// Construct a TaylorSeries object with given directional derivatives.
explicit TaylorSeries(const std::array<V, N + 1>& derivatives)
: _derivatives(derivatives)
{}
/// Evaluate the Taylor series object with given directional derivatives.
auto operator()(const T& t)
{
auto res = _derivatives[0];
auto factor = t;
For<1, N + 1>([&](auto&& i) constexpr {
res += factor * _derivatives[i];
factor *= t / static_cast<T>(i + 1);
});
return res;
}
/// Return the directional derivatives of this TaylorSeries.
auto derivatives()
{
return _derivatives;
}
private:
/// The directional derivatives of the function up to Nth order.
std::array<V, N + 1> _derivatives;
};
/// Return a TaylorSeries of a scalar or vector function *f* along a direction *v* at *x*.
template<typename Fun, typename...Vecs, typename... Args>
auto taylorseries(const Fun& f, const Along<Vecs...>& along, const At<Args...>& at)
{
auto data = derivatives(f, along, at);
constexpr auto N = data.size() - 1;
using V = typename decltype(data)::value_type;
return TaylorSeries<N, V>(data);
}
} // namespace detail
using detail::taylorseries;
} // namespace autodiff

View File

@ -0,0 +1,89 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// pybind11 includes
#include "pybind11.hxx"
// autodiff includes
#include <autodiff/forward/dual.hpp>
#include <autodiff/forward/dual/eigen.hpp>
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
#define PYBIND11_MAKE_OPAQUE_EIGEN_TYPES(scalar) \
PYBIND11_MAKE_OPAQUE(autodiff::VectorX##scalar##0th); \
PYBIND11_MAKE_OPAQUE(autodiff::VectorX##scalar##1st); \
PYBIND11_MAKE_OPAQUE(autodiff::VectorX##scalar##2nd); \
PYBIND11_MAKE_OPAQUE(autodiff::VectorX##scalar##3rd); \
PYBIND11_MAKE_OPAQUE(autodiff::VectorX##scalar##4th); \
PYBIND11_MAKE_OPAQUE(autodiff::MatrixX##scalar##0th); \
PYBIND11_MAKE_OPAQUE(autodiff::MatrixX##scalar##1st); \
PYBIND11_MAKE_OPAQUE(autodiff::MatrixX##scalar##2nd); \
PYBIND11_MAKE_OPAQUE(autodiff::MatrixX##scalar##3rd); \
PYBIND11_MAKE_OPAQUE(autodiff::MatrixX##scalar##4th); \
PYBIND11_MAKE_OPAQUE(autodiff::ArrayX##scalar##0th); \
PYBIND11_MAKE_OPAQUE(autodiff::ArrayX##scalar##1st); \
PYBIND11_MAKE_OPAQUE(autodiff::ArrayX##scalar##2nd); \
PYBIND11_MAKE_OPAQUE(autodiff::ArrayX##scalar##3rd); \
PYBIND11_MAKE_OPAQUE(autodiff::ArrayX##scalar##4th); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::VectorX##scalar##0th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::VectorX##scalar##1st>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::VectorX##scalar##2nd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::VectorX##scalar##3rd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::VectorX##scalar##4th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::MatrixX##scalar##0th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::MatrixX##scalar##1st>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::MatrixX##scalar##2nd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::MatrixX##scalar##3rd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::MatrixX##scalar##4th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::ArrayX##scalar##0th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::ArrayX##scalar##1st>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::ArrayX##scalar##2nd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::ArrayX##scalar##3rd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<autodiff::ArrayX##scalar##4th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::VectorX##scalar##0th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::VectorX##scalar##1st>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::VectorX##scalar##2nd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::VectorX##scalar##3rd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::VectorX##scalar##4th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::MatrixX##scalar##0th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::MatrixX##scalar##1st>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::MatrixX##scalar##2nd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::MatrixX##scalar##3rd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::MatrixX##scalar##4th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::ArrayX##scalar##0th>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::ArrayX##scalar##1st>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::ArrayX##scalar##2nd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::ArrayX##scalar##3rd>); \
PYBIND11_MAKE_OPAQUE(Eigen::Ref<const autodiff::ArrayX##scalar##4th>);
PYBIND11_MAKE_OPAQUE_EIGEN_TYPES(real);
PYBIND11_MAKE_OPAQUE_EIGEN_TYPES(dual);

View File

@ -0,0 +1,33 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// autodiff includes
#include <autodiff/reverse/var/var.hpp>

View File

@ -0,0 +1,223 @@
// _ _
// _ _|_ _ _|o_|__|_
// (_||_||_(_)(_|| | |
//
// automatic differentiation made easier in C++
// https://github.com/autodiff/autodiff
//
// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
//
// Copyright (c) 2018-2022 Allan Leal
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
// Eigen includes
#include <Eigen/Core>
// autodiff includes
#include <autodiff/common/eigen.hpp>
#include <autodiff/common/meta.hpp>
#include <autodiff/reverse/var/var.hpp>
//------------------------------------------------------------------------------
// SUPPORT FOR EIGEN MATRICES AND VECTORS OF VAR
//------------------------------------------------------------------------------
namespace Eigen {
template<typename T>
struct NumTraits;
template<typename T>
struct NumTraits<autodiff::Variable<T>> : NumTraits<T> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
typedef autodiff::Variable<T> Real;
typedef autodiff::Variable<T> NonInteger;
typedef autodiff::Variable<T> Nested;
enum
{
IsComplex = 0,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1,
ReadCost = 1,
AddCost = 3,
MulCost = 3
};
};
template<typename T, typename BinOp>
struct ScalarBinaryOpTraits<autodiff::Variable<T>, T, BinOp>
{
typedef autodiff::Variable<T> ReturnType;
};
template<typename T, typename BinOp>
struct ScalarBinaryOpTraits<T, autodiff::Variable<T>, BinOp>
{
typedef autodiff::Variable<T> ReturnType;
};
template<typename T>
struct NumTraits<autodiff::detail::ExprPtr<T>> : NumTraits<T> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
typedef autodiff::Variable<T> Real;
typedef autodiff::Variable<T> NonInteger;
typedef autodiff::Variable<T> Nested;
enum
{
IsComplex = 0,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1,
ReadCost = 1,
AddCost = 3,
MulCost = 3
};
};
template<typename T, typename BinOp>
struct ScalarBinaryOpTraits<autodiff::detail::ExprPtr<T>, T, BinOp>
{
typedef autodiff::Variable<T> ReturnType;
};
template<typename T, typename BinOp>
struct ScalarBinaryOpTraits<T, autodiff::detail::ExprPtr<T>, BinOp>
{
typedef autodiff::Variable<T> ReturnType;
};
} // namespace Eigen
namespace autodiff {
namespace detail {
template<typename T, int Rows, int MaxRows>
using Vec = Eigen::Matrix<T, Rows, 1, 0, MaxRows, 1>;
template<typename T, int Rows, int Cols, int MaxRows, int MaxCols>
using Mat = Eigen::Matrix<T, Rows, Cols, 0, MaxRows, MaxCols>;
/// Return the gradient vector of variable y with respect to variables x.
template<typename T, typename X>
auto gradient(const Variable<T>& y, Eigen::DenseBase<X>& x)
{
using U = VariableValueType<T>;
using ScalarX = typename X::Scalar;
static_assert(isVariable<ScalarX>, "Argument x is not a vector with Variable<T> (aka var) objects..");
constexpr auto isVec = X::IsVectorAtCompileTime;
static_assert(isVec, "Argument x is not a vector.");
constexpr auto Rows = X::RowsAtCompileTime;
constexpr auto MaxRows = X::MaxRowsAtCompileTime;
const auto n = x.size();
using Gradient = Vec<U, Rows, MaxRows>;
Gradient g = Gradient::Zero(n);
for(auto i = 0; i < n; ++i)
x[i].expr->bind_value(&g[i]);
y.expr->propagate(1.0);
for(auto i = 0; i < n; ++i)
x[i].expr->bind_value(nullptr);
return g;
}
/// Return the Hessian matrix of variable y with respect to variables x.
template<typename T, typename X, typename GradientVec>
auto hessian(const Variable<T>& y, Eigen::DenseBase<X>& x, GradientVec& g)
{
using U = VariableValueType<T>;
using ScalarX = typename X::Scalar;
static_assert(isVariable<ScalarX>, "Argument x is not a vector with Variable<T> (aka var) objects.");
using ScalarG = typename GradientVec::Scalar;
static_assert(std::is_same_v<U, ScalarG>, "Argument g does not have the same arithmetic type as y.");
constexpr auto Rows = X::RowsAtCompileTime;
constexpr auto MaxRows = X::MaxRowsAtCompileTime;
const auto n = x.size();
// Form a vector containing gradient expressions for each variable
using ExpressionGradient = Vec<ScalarX, Rows, MaxRows>;
ExpressionGradient G(n);
for(auto k = 0; k < n; ++k)
x[k].expr->bind_expr(&G(k).expr);
/* Build a full gradient expression in DFS tree traversal, updating
* gradient expressions when encountering variables
*/
y.expr->propagatex(constant<T>(1.0));
for(auto k = 0; k < n; ++k) {
x[k].expr->bind_expr(nullptr);
}
// Read the gradient value from gradient expressions' cached values
g.resize(n);
for(auto i = 0; i < n; ++i)
g[i] = val(G[i]);
// Form a numeric hessian using the gradient expressions
using Hessian = Mat<U, Rows, Rows, MaxRows, MaxRows>;
Hessian H = Hessian::Zero(n, n);
for(auto i = 0; i < n; ++i)
{
for(auto k = 0; k < n; ++k)
x[k].expr->bind_value(&H(i, k));
// Propagate a second derivative value calculation down the gradient expression tree for variable i
G[i].expr->propagate(1.0);
for(auto k = 0; k < n; ++k)
x[k].expr->bind_value(nullptr);
}
return H;
}
/// Return the Hessian matrix of variable y with respect to variables x.
template<typename T, typename X>
auto hessian(const Variable<T>& y, Eigen::DenseBase<X>& x)
{
using U = VariableValueType<T>;
constexpr auto Rows = X::RowsAtCompileTime;
constexpr auto MaxRows = X::MaxRowsAtCompileTime;
Vec<U, Rows, MaxRows> g;
return hessian(y, x, g);
}
} // namespace detail
AUTODIFF_DEFINE_EIGEN_TYPEDEFS_ALL_SIZES(autodiff::var, var)
using detail::gradient;
using detail::hessian;
} // namespace autodiff

File diff suppressed because it is too large Load Diff