tboudreaux/libcomposition
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

build(CppAD): brought in CppAD for autodiff

Browse Source 
we need an autodiff library at some point (or we need to roll our own but I do not think that makes sense). CppAD is well tested and header only and easy to include. It is also Liscene compatible with GPL v3.0. Here we bring it in as a dependency
This commit is contained in:
Emily Boudreaux
2025-06-19 14:51:02 -04:00
parent 0f68f2959e
commit 2bca6e447c
368 changed files with 108394 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

661
build-config/cppad/include/cppad/example/atomic_three/mat_mul.hpp Normal file
View File

@@ -0,0 +1,661 @@
# ifndef CPPAD_EXAMPLE_ATOMIC_THREE_MAT_MUL_HPP
# define CPPAD_EXAMPLE_ATOMIC_THREE_MAT_MUL_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-19 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
/*
$begin atomic_three_mat_mul.hpp$$
$spell
Taylor
ty
px
CppAD
jac
hes
nr
nc
afun
mul
$$
$section Matrix Multiply as an Atomic Operation$$
$head See Also$$
$cref atomic_two_eigen_mat_mul.hpp$$
$head Purpose$$
Use scalar $code double$$ operations in an $cref atomic_three$$ operation
that computes the matrix product for $code AD<double$$ operations.
$subhead parameter_x$$
This example demonstrates the use of the
$cref/parameter_x/atomic_three/parameter_x/$$
argument to the $cref atomic_three$$ virtual functions.
$subhead type_x$$
This example also demonstrates the use of the
$cref/type_x/atomic_three/type_x/$$
argument to the $cref atomic_three$$ virtual functions.
$head Matrix Dimensions$$
The first three components of the argument vector $icode ax$$
in the call $icode%afun%(%ax%, %ay%)%$$
are parameters and contain the matrix dimensions.
This enables them to be different for each use of the same atomic
function $icode afun$$.
These dimensions are:
$table
$icode%ax%[0]%$$ $pre $$
$cnext $icode nr_left$$ $pre $$
$cnext number of rows in the left matrix and result matrix
$rend
$icode%ax%[1]%$$ $pre $$
$cnext $icode n_middle$$ $pre $$
$cnext columns in the left matrix and rows in right matrix
$rend
$icode%ax%[2]%$$ $pre $$
$cnext $icode nc_right$$ $pre $$
$cnext number of columns in the right matrix and result matrix
$tend
$head Left Matrix$$
The number of elements in the left matrix is
$codei%
%n_left% = %nr_left% * %n_middle%
%$$
The elements are in
$icode%ax%[3]%$$ through $icode%ax%[2+%n_left%]%$$ in row major order.
$head Right Matrix$$
The number of elements in the right matrix is
$codei%
%n_right% = %n_middle% * %nc_right%
%$$
The elements are in
$icode%ax%[3+%n_left%]%$$ through
$icode%ax%[2+%n_left%+%n_right%]%$$ in row major order.
$head Result Matrix$$
The number of elements in the result matrix is
$codei%
%n_result% = %nr_left% * %nc_right%
%$$
The elements are in
$icode%ay%[0]%$$ through $icode%ay%[%n_result%-1]%$$ in row major order.
$head Start Class Definition$$
$srccode%cpp% */
# include <cppad/cppad.hpp>
namespace { // Begin empty namespace
using CppAD::vector;
//
// matrix result = left * right
class atomic_mat_mul : public CppAD::atomic_three<double> {
/* %$$
$head Constructor$$
$srccode%cpp% */
public:
// ---------------------------------------------------------------------
// constructor
atomic_mat_mul(void) : CppAD::atomic_three<double>("mat_mul")
{ }
private:
/* %$$
$head Left Operand Element Index$$
Index in the Taylor coefficient matrix $icode tx$$ of a left matrix element.
$srccode%cpp% */
size_t left(
size_t i , // left matrix row index
size_t j , // left matrix column index
size_t k , // Taylor coeffocient order
size_t nk , // number of Taylor coefficients in tx
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{ assert( i < nr_left );
assert( j < n_middle );
return (3 + i * n_middle + j) * nk + k;
}
/* %$$
$head Right Operand Element Index$$
Index in the Taylor coefficient matrix $icode tx$$ of a right matrix element.
$srccode%cpp% */
size_t right(
size_t i , // right matrix row index
size_t j , // right matrix column index
size_t k , // Taylor coeffocient order
size_t nk , // number of Taylor coefficients in tx
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{ assert( i < n_middle );
assert( j < nc_right );
size_t offset = 3 + nr_left * n_middle;
return (offset + i * nc_right + j) * nk + k;
}
/* %$$
$head Result Element Index$$
Index in the Taylor coefficient matrix $icode ty$$ of a result matrix element.
$srccode%cpp% */
size_t result(
size_t i , // result matrix row index
size_t j , // result matrix column index
size_t k , // Taylor coeffocient order
size_t nk , // number of Taylor coefficients in ty
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{ assert( i < nr_left );
assert( j < nc_right );
return (i * nc_right + j) * nk + k;
}
/* %$$
$head Forward Matrix Multiply$$
Forward mode multiply Taylor coefficients in $icode tx$$ and sum into
$icode ty$$ (for one pair of left and right orders)
$srccode%cpp% */
void forward_multiply(
size_t k_left , // order for left coefficients
size_t k_right , // order for right coefficients
const vector<double>& tx , // domain space Taylor coefficients
vector<double>& ty , // range space Taylor coefficients
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t nk = tx.size() / nx;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
assert( nk == ty.size() / ny );
# endif
//
size_t k_result = k_left + k_right;
assert( k_result < nk );
//
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ double sum = 0.0;
for(size_t ell = 0; ell < n_middle; ell++)
{ size_t i_left = left(
i, ell, k_left, nk, nr_left, n_middle, nc_right
);
size_t i_right = right(
ell, j, k_right, nk, nr_left, n_middle, nc_right
);
sum += tx[i_left] * tx[i_right];
}
size_t i_result = result(
i, j, k_result, nk, nr_left, n_middle, nc_right
);
ty[i_result] += sum;
}
}
}
/* %$$
$head Reverse Matrix Multiply$$
Reverse mode partials of Taylor coefficients and sum into $icode px$$
(for one pair of left and right orders)
$srccode%cpp% */
void reverse_multiply(
size_t k_left , // order for left coefficients
size_t k_right , // order for right coefficients
const vector<double>& tx , // domain space Taylor coefficients
const vector<double>& ty , // range space Taylor coefficients
vector<double>& px , // partials w.r.t. tx
const vector<double>& py , // partials w.r.t. ty
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t nk = tx.size() / nx;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
assert( nk == ty.size() / ny );
# endif
assert( tx.size() == px.size() );
assert( ty.size() == py.size() );
//
size_t k_result = k_left + k_right;
assert( k_result < nk );
//
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ size_t i_result = result(
i, j, k_result, nk, nr_left, n_middle, nc_right
);
for(size_t ell = 0; ell < n_middle; ell++)
{ size_t i_left = left(
i, ell, k_left, nk, nr_left, n_middle, nc_right
);
size_t i_right = right(
ell, j, k_right, nk, nr_left, n_middle, nc_right
);
// sum += tx[i_left] * tx[i_right];
px[i_left] += tx[i_right] * py[i_result];
px[i_right] += tx[i_left] * py[i_result];
}
}
}
return;
}
/* %$$
$head for_type$$
Routine called by CppAD during $cref/afun(ax, ay)/atomic_three_afun/$$.
$srccode%cpp% */
// calculate type_y
virtual bool for_type(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
vector<CppAD::ad_type_enum>& type_y )
{ assert( parameter_x.size() == type_x.size() );
bool ok = true;
ok &= type_x[0] == CppAD::constant_enum;
ok &= type_x[1] == CppAD::constant_enum;
ok &= type_x[2] == CppAD::constant_enum;
if( ! ok )
return false;
//
size_t nr_left = size_t( parameter_x[0] );
size_t n_middle = size_t( parameter_x[1] );
size_t nc_right = size_t( parameter_x[2] );
//
ok &= type_x.size() == 3 + (nr_left + nc_right) * n_middle;
ok &= type_y.size() == n_middle * nc_right;
if( ! ok )
return false;
//
// commpute type_y
size_t nk = 1; // number of orders
size_t k = 0; // order
for(size_t i = 0; i < nr_left; ++i)
{ for(size_t j = 0; j < nc_right; ++j)
{ // compute type for result[i, j]
CppAD::ad_type_enum type_yij = CppAD::constant_enum;
for(size_t ell = 0; ell < n_middle; ++ell)
{ // index for left(i, ell)
size_t i_left = left(
i, ell, k, nk, nr_left, n_middle, nc_right
);
// indx for right(ell, j)
size_t i_right = right(
ell, j, k, nk, nr_left, n_middle, nc_right
);
// multiplication on left or right by the constant zero
// always results in a constant
bool zero_left = type_x[i_left] == CppAD::constant_enum;
zero_left &= parameter_x[i_left] == 0.0;
bool zero_right = type_x[i_right] == CppAD::constant_enum;
zero_right &= parameter_x[i_right] == 0.0;
if( ! (zero_left | zero_right) )
{ type_yij = std::max(type_yij, type_x[i_left] );
type_yij = std::max(type_yij, type_x[i_right] );
}
}
size_t i_result = result(
i, j, k, nk, nr_left, n_middle, nc_right
);
type_y[i_result] = type_yij;
}
}
return true;
}
/* %$$
$head forward$$
Routine called by CppAD during $cref Forward$$ mode.
$srccode%cpp% */
virtual bool forward(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
size_t need_y ,
size_t q ,
size_t p ,
const vector<double>& tx ,
vector<double>& ty )
{ size_t n_order = p + 1;
size_t nr_left = size_t( tx[ 0 * n_order + 0 ] );
size_t n_middle = size_t( tx[ 1 * n_order + 0 ] );
size_t nc_right = size_t( tx[ 2 * n_order + 0 ] );
# ifndef NDEBUG
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
size_t i, j, ell;
// initialize result as zero
size_t k;
for(i = 0; i < nr_left; i++)
{ for(j = 0; j < nc_right; j++)
{ for(k = q; k <= p; k++)
{ size_t i_result = result(
i, j, k, n_order, nr_left, n_middle, nc_right
);
ty[i_result] = 0.0;
}
}
}
for(k = q; k <= p; k++)
{ // sum the produces that result in order k
for(ell = 0; ell <= k; ell++)
forward_multiply(
ell, k - ell, tx, ty, nr_left, n_middle, nc_right
);
}
// all orders are implemented, so always return true
return true;
}
/* %$$
$head reverse$$
Routine called by CppAD during $cref Reverse$$ mode.
$srccode%cpp% */
virtual bool reverse(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
size_t p ,
const vector<double>& tx ,
const vector<double>& ty ,
vector<double>& px ,
const vector<double>& py )
{ size_t n_order = p + 1;
size_t nr_left = size_t( tx[ 0 * n_order + 0 ] );
size_t n_middle = size_t( tx[ 1 * n_order + 0 ] );
size_t nc_right = size_t( tx[ 2 * n_order + 0 ] );
# ifndef NDEBUG
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
assert( px.size() == tx.size() );
assert( py.size() == ty.size() );
// initialize summation
for(size_t i = 0; i < px.size(); i++)
px[i] = 0.0;
// number of orders to differentiate
size_t k = n_order;
while(k--)
{ // differentiate the produces that result in order k
for(size_t ell = 0; ell <= k; ell++)
reverse_multiply(
ell, k - ell, tx, ty, px, py, nr_left, n_middle, nc_right
);
}
// all orders are implented, so always return true
return true;
}
/* %$$
$head jac_sparsity$$
$srccode%cpp% */
// Jacobian sparsity routine called by CppAD
virtual bool jac_sparsity(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
bool dependency ,
const vector<bool>& select_x ,
const vector<bool>& select_y ,
CppAD::sparse_rc< vector<size_t> >& pattern_out )
{
size_t n = select_x.size();
size_t m = select_y.size();
assert( parameter_x.size() == n );
assert( type_x.size() == n );
//
size_t nr_left = size_t( parameter_x[0] );
size_t n_middle = size_t( parameter_x[1] );
size_t nc_right = size_t( parameter_x[2] );
size_t nk = 1; // only one order
size_t k = 0; // order zero
//
// count number of non-zeros in sparsity pattern
size_t nnz = 0;
for(size_t i = 0; i < nr_left; ++i)
{ for(size_t j = 0; j < nc_right; ++j)
{ size_t i_result = result(
i, j, k, nk, nr_left, n_middle, nc_right
);
if( select_y[i_result] )
{ for(size_t ell = 0; ell < n_middle; ++ell)
{ size_t i_left = left(
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right(
ell, j, k, nk, nr_left, n_middle, nc_right
);
bool zero_left =
type_x[i_left] == CppAD::constant_enum;
zero_left &= parameter_x[i_left] == 0.0;
bool zero_right =
type_x[i_right] == CppAD::constant_enum;
zero_right &= parameter_x[i_right] == 0.0;
if( ! (zero_left | zero_right ) )
{ bool var_left =
type_x[i_left] == CppAD::variable_enum;
bool var_right =
type_x[i_right] == CppAD::variable_enum;
if( select_x[i_left] & var_left )
++nnz;
if( select_x[i_right] & var_right )
++nnz;
}
}
}
}
}
//
// fill in the sparsity pattern
pattern_out.resize(m, n, nnz);
size_t idx = 0;
for(size_t i = 0; i < nr_left; ++i)
{ for(size_t j = 0; j < nc_right; ++j)
{ size_t i_result = result(
i, j, k, nk, nr_left, n_middle, nc_right
);
if( select_y[i_result] )
{ for(size_t ell = 0; ell < n_middle; ++ell)
{ size_t i_left = left(
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right(
ell, j, k, nk, nr_left, n_middle, nc_right
);
bool zero_left =
type_x[i_left] == CppAD::constant_enum;
zero_left &= parameter_x[i_left] == 0.0;
bool zero_right =
type_x[i_right] == CppAD::constant_enum;
zero_right &= parameter_x[i_right] == 0.0;
if( ! (zero_left | zero_right ) )
{ bool var_left =
type_x[i_left] == CppAD::variable_enum;
bool var_right =
type_x[i_right] == CppAD::variable_enum;
if( select_x[i_left] & var_left )
pattern_out.set(idx++, i_result, i_left);
if( select_x[i_right] & var_right )
pattern_out.set(idx++, i_result, i_right);
}
}
}
}
}
assert( idx == nnz );
//
return true;
}
/* %$$
$head hes_sparsity$$
$srccode%cpp% */
// Jacobian sparsity routine called by CppAD
virtual bool hes_sparsity(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
const vector<bool>& select_x ,
const vector<bool>& select_y ,
CppAD::sparse_rc< vector<size_t> >& pattern_out )
{
size_t n = select_x.size();
assert( parameter_x.size() == n );
assert( type_x.size() == n );
//
size_t nr_left = size_t( parameter_x[0] );
size_t n_middle = size_t( parameter_x[1] );
size_t nc_right = size_t( parameter_x[2] );
size_t nk = 1; // only one order
size_t k = 0; // order zero
//
// count number of non-zeros in sparsity pattern
size_t nnz = 0;
for(size_t i = 0; i < nr_left; ++i)
{ for(size_t j = 0; j < nc_right; ++j)
{ size_t i_result = result(
i, j, k, nk, nr_left, n_middle, nc_right
);
if( select_y[i_result] )
{ for(size_t ell = 0; ell < n_middle; ++ell)
{ // i_left depends on i, ell
size_t i_left = left(
i, ell, k, nk, nr_left, n_middle, nc_right
);
// i_right depens on ell, j
size_t i_right = right(
ell, j, k, nk, nr_left, n_middle, nc_right
);
bool var_left = select_x[i_left] &
(type_x[i_left] == CppAD::variable_enum);
bool var_right = select_x[i_right] &
(type_x[i_right] == CppAD::variable_enum);
if( var_left & var_right )
nnz += 2;
}
}
}
}
//
// fill in the sparsity pattern
pattern_out.resize(n, n, nnz);
size_t idx = 0;
for(size_t i = 0; i < nr_left; ++i)
{ for(size_t j = 0; j < nc_right; ++j)
{ size_t i_result = result(
i, j, k, nk, nr_left, n_middle, nc_right
);
if( select_y[i_result] )
{ for(size_t ell = 0; ell < n_middle; ++ell)
{ size_t i_left = left(
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right(
ell, j, k, nk, nr_left, n_middle, nc_right
);
bool var_left = select_x[i_left] &
(type_x[i_left] == CppAD::variable_enum);
bool var_right = select_x[i_right] &
(type_x[i_right] == CppAD::variable_enum);
if( var_left & var_right )
{ // Cannot possibly set the same (i_left, i_right)
// pair twice.
assert( i_left != i_right );
pattern_out.set(idx++, i_left, i_right);
pattern_out.set(idx++, i_right, i_left);
}
}
}
}
}
assert( idx == nnz );
//
return true;
}
/* %$$
$head rev_depend$$
Routine called when a function using $code mat_mul$$ is optimized.
$srccode%cpp% */
// calculate depend_x
virtual bool rev_depend(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
vector<bool>& depend_x ,
const vector<bool>& depend_y )
{ assert( parameter_x.size() == depend_x.size() );
assert( parameter_x.size() == type_x.size() );
bool ok = true;
//
size_t nr_left = size_t( parameter_x[0] );
size_t n_middle = size_t( parameter_x[1] );
size_t nc_right = size_t( parameter_x[2] );
//
ok &= depend_x.size() == 3 + (nr_left + nc_right) * n_middle;
ok &= depend_y.size() == n_middle * nc_right;
if( ! ok )
return false;
//
// initialize depend_x
for(size_t ell = 0; ell < 3; ++ell)
depend_x[ell] = true; // always need these parameters
for(size_t ell = 3; ell < depend_x.size(); ++ell)
depend_x[ell] = false; // initialize as false
//
// commpute depend_x
size_t nk = 1; // number of orders
size_t k = 0; // order
for(size_t i = 0; i < nr_left; ++i)
{ for(size_t j = 0; j < nc_right; ++j)
{ // check depend for result[i, j]
size_t i_result = result(
i, j, k, nk, nr_left, n_middle, nc_right
);
if( depend_y[i_result] )
{ for(size_t ell = 0; ell < n_middle; ++ell)
{ // index for left(i, ell)
size_t i_left = left(
i, ell, k, nk, nr_left, n_middle, nc_right
);
// indx for right(ell, j)
size_t i_right = right(
ell, j, k, nk, nr_left, n_middle, nc_right
);
bool zero_left =
type_x[i_left] == CppAD::constant_enum;
zero_left &= parameter_x[i_left] == 0.0;
bool zero_right =
type_x[i_right] == CppAD::constant_enum;
zero_right &= parameter_x[i_right] == 0.0;
if( ! zero_right )
depend_x[i_left] = true;
if( ! zero_left )
depend_x[i_right] = true;
}
}
}
}
return true;
}
/* %$$
$head End Class Definition$$
$srccode%cpp% */
}; // End of mat_mul class
} // End empty namespace
/* %$$
$comment end nospell$$
$end
*/
# endif

376
build-config/cppad/include/cppad/example/atomic_two/eigen_cholesky.hpp Normal file
View File

@@ -0,0 +1,376 @@
# ifndef CPPAD_EXAMPLE_ATOMIC_TWO_EIGEN_CHOLESKY_HPP
# define CPPAD_EXAMPLE_ATOMIC_TWO_EIGEN_CHOLESKY_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-19 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
/*
$begin atomic_two_eigen_cholesky.hpp$$
$spell
Eigen
Taylor
Cholesky
op
$$
$section atomic_two Eigen Cholesky Factorization Class$$
$head Purpose$$
Construct an atomic operation that computes a lower triangular matrix
$latex L $$ such that $latex L L^\R{T} = A$$
for any positive integer $latex p$$
and symmetric positive definite matrix $latex A \in \B{R}^{p \times p}$$.
$head Start Class Definition$$
$srccode%cpp% */
# include <cppad/cppad.hpp>
# include <Eigen/Dense>
/* %$$
$head Public$$
$subhead Types$$
$srccode%cpp% */
namespace { // BEGIN_EMPTY_NAMESPACE
template <class Base>
class atomic_eigen_cholesky : public CppAD::atomic_base<Base> {
public:
// -----------------------------------------------------------
// type of elements during calculation of derivatives
typedef Base scalar;
// type of elements during taping
typedef CppAD::AD<scalar> ad_scalar;
//
// type of matrix during calculation of derivatives
typedef Eigen::Matrix<
scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
// type of matrix during taping
typedef Eigen::Matrix<
ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
//
// lower triangular scalar matrix
typedef Eigen::TriangularView<matrix, Eigen::Lower> lower_view;
/* %$$
$subhead Constructor$$
$srccode%cpp% */
// constructor
atomic_eigen_cholesky(void) : CppAD::atomic_base<Base>(
"atom_eigen_cholesky" ,
CppAD::atomic_base<Base>::set_sparsity_enum
)
{ }
/* %$$
$subhead op$$
$srccode%cpp% */
// use atomic operation to invert an AD matrix
ad_matrix op(const ad_matrix& arg)
{ size_t nr = size_t( arg.rows() );
size_t ny = ( (nr + 1 ) * nr ) / 2;
size_t nx = 1 + ny;
assert( nr == size_t( arg.cols() ) );
// -------------------------------------------------------------------
// packed version of arg
CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
size_t index = 0;
packed_arg[index++] = ad_scalar( nr );
// lower triangle of symmetric matrix A
for(size_t i = 0; i < nr; i++)
{ for(size_t j = 0; j <= i; j++)
packed_arg[index++] = arg( long(i), long(j) );
}
assert( index == nx );
// -------------------------------------------------------------------
// packed version of result = arg^{-1}.
// This is an atomic_base function call that CppAD uses to
// store the atomic operation on the tape.
CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
(*this)(packed_arg, packed_result);
// -------------------------------------------------------------------
// unpack result matrix L
ad_matrix result = ad_matrix::Zero( long(nr), long(nr) );
index = 0;
for(size_t i = 0; i < nr; i++)
{ for(size_t j = 0; j <= i; j++)
result( long(i), long(j) ) = packed_result[index++];
}
return result;
}
/* %$$
$head Private$$
$subhead Variables$$
$srccode%cpp% */
private:
// -------------------------------------------------------------
// one forward mode vector of matrices for argument and result
CppAD::vector<matrix> f_arg_, f_result_;
// one reverse mode vector of matrices for argument and result
CppAD::vector<matrix> r_arg_, r_result_;
// -------------------------------------------------------------
/* %$$
$subhead forward$$
$srccode%cpp% */
// forward mode routine called by CppAD
virtual bool forward(
// lowest order Taylor coefficient we are evaluating
size_t p ,
// highest order Taylor coefficient we are evaluating
size_t q ,
// which components of x are variables
const CppAD::vector<bool>& vx ,
// which components of y are variables
CppAD::vector<bool>& vy ,
// tx [ j * (q+1) + k ] is x_j^k
const CppAD::vector<scalar>& tx ,
// ty [ i * (q+1) + k ] is y_i^k
CppAD::vector<scalar>& ty
)
{ size_t n_order = q + 1;
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
size_t ny = ((nr + 1) * nr) / 2;
# ifndef NDEBUG
size_t nx = 1 + ny;
# endif
assert( vx.size() == 0 || nx == vx.size() );
assert( vx.size() == 0 || ny == vy.size() );
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
//
// -------------------------------------------------------------------
// make sure f_arg_ and f_result_ are large enough
assert( f_arg_.size() == f_result_.size() );
if( f_arg_.size() < n_order )
{ f_arg_.resize(n_order);
f_result_.resize(n_order);
//
for(size_t k = 0; k < n_order; k++)
{ f_arg_[k].resize( long(nr), long(nr) );
f_result_[k].resize( long(nr), long(nr) );
}
}
// -------------------------------------------------------------------
// unpack tx into f_arg_
for(size_t k = 0; k < n_order; k++)
{ size_t index = 1;
// unpack arg values for this order
for(long i = 0; i < long(nr); i++)
{ for(long j = 0; j <= i; j++)
{ f_arg_[k](i, j) = tx[ index * n_order + k ];
f_arg_[k](j, i) = f_arg_[k](i, j);
index++;
}
}
}
// -------------------------------------------------------------------
// result for each order
// (we could avoid recalculting f_result_[k] for k=0,...,p-1)
//
Eigen::LLT<matrix> cholesky(f_arg_[0]);
f_result_[0] = cholesky.matrixL();
lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
for(size_t k = 1; k < n_order; k++)
{ // initialize sum as A_k
matrix f_sum = f_arg_[k];
// compute A_k - B_k
for(size_t ell = 1; ell < k; ell++)
f_sum -= f_result_[ell] * f_result_[k-ell].transpose();
// compute L_0^{-1} * (A_k - B_k) * L_0^{-T}
matrix temp = L_0.template solve<Eigen::OnTheLeft>(f_sum);
temp = L_0.transpose().template solve<Eigen::OnTheRight>(temp);
// divide the diagonal by 2
for(long i = 0; i < long(nr); i++)
temp(i, i) /= scalar(2.0);
// L_k = L_0 * low[ L_0^{-1} * (A_k - B_k) * L_0^{-T} ]
lower_view view = temp.template triangularView<Eigen::Lower>();
f_result_[k] = f_result_[0] * view;
}
// -------------------------------------------------------------------
// pack result_ into ty
for(size_t k = 0; k < n_order; k++)
{ size_t index = 0;
for(long i = 0; i < long(nr); i++)
{ for(long j = 0; j <= i; j++)
{ ty[ index * n_order + k ] = f_result_[k](i, j);
index++;
}
}
}
// -------------------------------------------------------------------
// check if we are computing vy
if( vx.size() == 0 )
return true;
// ------------------------------------------------------------------
// This is a very dumb algorithm that over estimates which
// elements of the inverse are variables (which is not efficient).
bool var = false;
for(size_t i = 0; i < ny; i++)
var |= vx[1 + i];
for(size_t i = 0; i < ny; i++)
vy[i] = var;
//
return true;
}
/* %$$
$subhead reverse$$
$srccode%cpp% */
// reverse mode routine called by CppAD
virtual bool reverse(
// highest order Taylor coefficient that we are computing derivative of
size_t q ,
// forward mode Taylor coefficients for x variables
const CppAD::vector<double>& tx ,
// forward mode Taylor coefficients for y variables
const CppAD::vector<double>& ty ,
// upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
CppAD::vector<double>& px ,
// derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
const CppAD::vector<double>& py
)
{ size_t n_order = q + 1;
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
# ifndef NDEBUG
size_t ny = ( (nr + 1 ) * nr ) / 2;
size_t nx = 1 + ny;
# endif
//
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
assert( px.size() == tx.size() );
assert( py.size() == ty.size() );
// -------------------------------------------------------------------
// make sure f_arg_ is large enough
assert( f_arg_.size() == f_result_.size() );
// must have previous run forward with order >= n_order
assert( f_arg_.size() >= n_order );
// -------------------------------------------------------------------
// make sure r_arg_, r_result_ are large enough
assert( r_arg_.size() == r_result_.size() );
if( r_arg_.size() < n_order )
{ r_arg_.resize(n_order);
r_result_.resize(n_order);
//
for(size_t k = 0; k < n_order; k++)
{ r_arg_[k].resize( long(nr), long(nr) );
r_result_[k].resize( long(nr), long(nr) );
}
}
// -------------------------------------------------------------------
// unpack tx into f_arg_
for(size_t k = 0; k < n_order; k++)
{ size_t index = 1;
// unpack arg values for this order
for(long i = 0; i < long(nr); i++)
{ for(long j = 0; j <= i; j++)
{ f_arg_[k](i, j) = tx[ index * n_order + k ];
f_arg_[k](j, i) = f_arg_[k](i, j);
index++;
}
}
}
// -------------------------------------------------------------------
// unpack py into r_result_
for(size_t k = 0; k < n_order; k++)
{ r_result_[k] = matrix::Zero( long(nr), long(nr) );
size_t index = 0;
for(long i = 0; i < long(nr); i++)
{ for(long j = 0; j <= i; j++)
{ r_result_[k](i, j) = py[ index * n_order + k ];
index++;
}
}
}
// -------------------------------------------------------------------
// initialize r_arg_ as zero
for(size_t k = 0; k < n_order; k++)
r_arg_[k] = matrix::Zero( long(nr), long(nr) );
// -------------------------------------------------------------------
// matrix reverse mode calculation
lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
//
for(size_t k1 = n_order; k1 > 1; k1--)
{ size_t k = k1 - 1;
//
// L_0^T * bar{L}_k
matrix tmp1 = L_0.transpose() * r_result_[k];
//
//low[ L_0^T * bar{L}_k ]
for(long i = 0; i < long(nr); i++)
tmp1(i, i) /= scalar(2.0);
matrix tmp2 = tmp1.template triangularView<Eigen::Lower>();
//
// L_0^{-T} low[ L_0^T * bar{L}_k ]
tmp1 = L_0.transpose().template solve<Eigen::OnTheLeft>( tmp2 );
//
// M_k = L_0^{-T} * low[ L_0^T * bar{L}_k ]^{T} L_0^{-1}
matrix M_k = L_0.transpose().template
solve<Eigen::OnTheLeft>( tmp1.transpose() );
//
// remove L_k and compute bar{B}_k
matrix barB_k = scalar(0.5) * ( M_k + M_k.transpose() );
r_arg_[k] += barB_k;
barB_k = scalar(-1.0) * barB_k;
//
// 2.0 * lower( bar{B}_k L_k )
matrix temp = scalar(2.0) * barB_k * f_result_[k];
temp = temp.template triangularView<Eigen::Lower>();
//
// remove C_k
r_result_[0] += temp;
//
// remove B_k
for(size_t ell = 1; ell < k; ell++)
{ // bar{L}_ell = 2 * lower( \bar{B}_k * L_{k-ell} )
temp = scalar(2.0) * barB_k * f_result_[k-ell];
r_result_[ell] += temp.template triangularView<Eigen::Lower>();
}
}
// M_0 = L_0^{-T} * low[ L_0^T * bar{L}_0 ]^{T} L_0^{-1}
matrix M_0 = L_0.transpose() * r_result_[0];
for(long i = 0; i < long(nr); i++)
M_0(i, i) /= scalar(2.0);
M_0 = M_0.template triangularView<Eigen::Lower>();
M_0 = L_0.template solve<Eigen::OnTheRight>( M_0 );
M_0 = L_0.transpose().template solve<Eigen::OnTheLeft>( M_0 );
// remove L_0
r_arg_[0] += scalar(0.5) * ( M_0 + M_0.transpose() );
// -------------------------------------------------------------------
// pack r_arg into px
// note that only the lower triangle of barA_k is stored in px
for(size_t k = 0; k < n_order; k++)
{ size_t index = 0;
px[ index * n_order + k ] = 0.0;
index++;
for(long i = 0; i < long(nr); i++)
{ for(long j = 0; j < i; j++)
{ px[ index * n_order + k ] = 2.0 * r_arg_[k](i, j);
index++;
}
px[ index * n_order + k] = r_arg_[k](i, i);
index++;
}
}
// -------------------------------------------------------------------
return true;
}
/* %$$
$head End Class Definition$$
$srccode%cpp% */
}; // End of atomic_eigen_cholesky class
} // END_EMPTY_NAMESPACE
/* %$$
$end
*/
# endif

395
build-config/cppad/include/cppad/example/atomic_two/eigen_mat_inv.hpp Normal file
View File

@@ -0,0 +1,395 @@
# ifndef CPPAD_EXAMPLE_ATOMIC_TWO_EIGEN_MAT_INV_HPP
# define CPPAD_EXAMPLE_ATOMIC_TWO_EIGEN_MAT_INV_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-19 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
/*
$begin atomic_two_eigen_mat_inv.hpp$$
$spell
Eigen
Taylor
$$
$section atomic_two Eigen Matrix Inversion Class$$
$head Purpose$$
Construct an atomic operation that computes the matrix inverse
$latex R = A^{-1}$$
for any positive integer $latex p$$
and invertible matrix $latex A \in \B{R}^{p \times p}$$.
$head Matrix Dimensions$$
This example puts the matrix dimension $latex p$$
in the atomic function arguments,
instead of the $cref/constructor/atomic_two_ctor/$$,
so it can be different for different calls to the atomic function.
$head Theory$$
$subhead Forward$$
The zero order forward mode Taylor coefficient is give by
$latex \[
R_0 = A_0^{-1}
\]$$
For $latex k = 1 , \ldots$$,
the $th k$$ order Taylor coefficient of $latex A R$$ is given by
$latex \[
0 = \sum_{\ell=0}^k A_\ell R_{k-\ell}
\] $$
Solving for $latex R_k$$ in terms of the coefficients
for $latex A$$ and the lower order coefficients for $latex R$$ we have
$latex \[
R_k = - R_0 \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
\] $$
Furthermore, once we have $latex R_k$$ we can compute the sum using
$latex \[
A_0 R_k = - \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
\] $$
$subhead Product of Three Matrices$$
Suppose $latex \bar{E}$$ is the derivative of the
scalar value function $latex s(E)$$ with respect to $latex E$$; i.e.,
$latex \[
\bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} }
\] $$
Also suppose that $latex t$$ is a scalar valued argument and
$latex \[
E(t) = B(t) C(t) D(t)
\] $$
It follows that
$latex \[
E'(t) = B'(t) C(t) D(t) + B(t) C'(t) D(t) + B(t) C(t) D'(t)
\] $$
$latex \[
(s \circ E)'(t)
=
\R{tr} [ \bar{E}^\R{T} E'(t) ]
\] $$
$latex \[
=
\R{tr} [ \bar{E}^\R{T} B'(t) C(t) D(t) ] +
\R{tr} [ \bar{E}^\R{T} B(t) C'(t) D(t) ] +
\R{tr} [ \bar{E}^\R{T} B(t) C(t) D'(t) ]
\] $$
$latex \[
=
\R{tr} [ B(t) D(t) \bar{E}^\R{T} B'(t) ] +
\R{tr} [ D(t) \bar{E}^\R{T} B(t) C'(t) ] +
\R{tr} [ \bar{E}^\R{T} B(t) C(t) D'(t) ]
\] $$
$latex \[
\bar{B} = \bar{E} (C D)^\R{T} \W{,}
\bar{C} = \B{R}^\R{T} \bar{E} D^\R{T} \W{,}
\bar{D} = (B C)^\R{T} \bar{E}
\] $$
$subhead Reverse$$
For $latex k > 0$$, reverse mode
eliminates $latex R_k$$ and expresses the function values
$latex s$$ in terms of the coefficients of $latex A$$
and the lower order coefficients of $latex R$$.
The effect on $latex \bar{R}_0$$
(of eliminating $latex R_k$$) is
$latex \[
\bar{R}_0
= \bar{R}_0 - \bar{R}_k \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)^\R{T}
= \bar{R}_0 + \bar{R}_k ( A_0 R_k )^\R{T}
\] $$
For $latex \ell = 1 , \ldots , k$$,
the effect on $latex \bar{R}_{k-\ell}$$ and $latex A_\ell$$
(of eliminating $latex R_k$$) is
$latex \[
\bar{A}_\ell = \bar{A}_\ell - R_0^\R{T} \bar{R}_k R_{k-\ell}^\R{T}
\] $$
$latex \[
\bar{R}_{k-\ell} = \bar{R}_{k-\ell} - ( R_0 A_\ell )^\R{T} \bar{R}_k
\] $$
We note that
$latex \[
R_0 '(t) A_0 (t) + R_0 (t) A_0 '(t) = 0
\] $$
$latex \[
R_0 '(t) = - R_0 (t) A_0 '(t) R_0 (t)
\] $$
The reverse mode formula that eliminates $latex R_0$$ is
$latex \[
\bar{A}_0
= \bar{A}_0 - R_0^\R{T} \bar{R}_0 R_0^\R{T}
\]$$
$nospell
$head Start Class Definition$$
$srccode%cpp% */
# include <cppad/cppad.hpp>
# include <Eigen/Core>
# include <Eigen/LU>
/* %$$
$head Public$$
$subhead Types$$
$srccode%cpp% */
namespace { // BEGIN_EMPTY_NAMESPACE
template <class Base>
class atomic_eigen_mat_inv : public CppAD::atomic_base<Base> {
public:
// -----------------------------------------------------------
// type of elements during calculation of derivatives
typedef Base scalar;
// type of elements during taping
typedef CppAD::AD<scalar> ad_scalar;
// type of matrix during calculation of derivatives
typedef Eigen::Matrix<
scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
// type of matrix during taping
typedef Eigen::Matrix<
ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
/* %$$
$subhead Constructor$$
$srccode%cpp% */
// constructor
atomic_eigen_mat_inv(void) : CppAD::atomic_base<Base>(
"atom_eigen_mat_inv" ,
CppAD::atomic_base<Base>::set_sparsity_enum
)
{ }
/* %$$
$subhead op$$
$srccode%cpp% */
// use atomic operation to invert an AD matrix
ad_matrix op(const ad_matrix& arg)
{ size_t nr = size_t( arg.rows() );
size_t ny = nr * nr;
size_t nx = 1 + ny;
assert( nr == size_t( arg.cols() ) );
// -------------------------------------------------------------------
// packed version of arg
CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
packed_arg[0] = ad_scalar( nr );
for(size_t i = 0; i < ny; i++)
packed_arg[1 + i] = arg.data()[i];
// -------------------------------------------------------------------
// packed version of result = arg^{-1}.
// This is an atomic_base function call that CppAD uses to
// store the atomic operation on the tape.
CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
(*this)(packed_arg, packed_result);
// -------------------------------------------------------------------
// unpack result matrix
ad_matrix result(nr, nr);
for(size_t i = 0; i < ny; i++)
result.data()[i] = packed_result[i];
return result;
}
/* %$$
$head Private$$
$subhead Variables$$
$srccode%cpp% */
private:
// -------------------------------------------------------------
// one forward mode vector of matrices for argument and result
CppAD::vector<matrix> f_arg_, f_result_;
// one reverse mode vector of matrices for argument and result
CppAD::vector<matrix> r_arg_, r_result_;
// -------------------------------------------------------------
/* %$$
$subhead forward$$
$srccode%cpp% */
// forward mode routine called by CppAD
virtual bool forward(
// lowest order Taylor coefficient we are evaluating
size_t p ,
// highest order Taylor coefficient we are evaluating
size_t q ,
// which components of x are variables
const CppAD::vector<bool>& vx ,
// which components of y are variables
CppAD::vector<bool>& vy ,
// tx [ j * (q+1) + k ] is x_j^k
const CppAD::vector<scalar>& tx ,
// ty [ i * (q+1) + k ] is y_i^k
CppAD::vector<scalar>& ty
)
{ size_t n_order = q + 1;
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
size_t ny = nr * nr;
# ifndef NDEBUG
size_t nx = 1 + ny;
# endif
assert( vx.size() == 0 || nx == vx.size() );
assert( vx.size() == 0 || ny == vy.size() );
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
//
// -------------------------------------------------------------------
// make sure f_arg_ and f_result_ are large enough
assert( f_arg_.size() == f_result_.size() );
if( f_arg_.size() < n_order )
{ f_arg_.resize(n_order);
f_result_.resize(n_order);
//
for(size_t k = 0; k < n_order; k++)
{ f_arg_[k].resize( long(nr), long(nr) );
f_result_[k].resize( long(nr), long(nr) );
}
}
// -------------------------------------------------------------------
// unpack tx into f_arg_
for(size_t k = 0; k < n_order; k++)
{ // unpack arg values for this order
for(size_t i = 0; i < ny; i++)
f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ];
}
// -------------------------------------------------------------------
// result for each order
// (we could avoid recalculting f_result_[k] for k=0,...,p-1)
//
f_result_[0] = f_arg_[0].inverse();
for(size_t k = 1; k < n_order; k++)
{ // initialize sum
matrix f_sum = matrix::Zero( long(nr), long(nr) );
// compute sum
for(size_t ell = 1; ell <= k; ell++)
f_sum -= f_arg_[ell] * f_result_[k-ell];
// result_[k] = arg_[0]^{-1} * sum_
f_result_[k] = f_result_[0] * f_sum;
}
// -------------------------------------------------------------------
// pack result_ into ty
for(size_t k = 0; k < n_order; k++)
{ for(size_t i = 0; i < ny; i++)
ty[ i * n_order + k ] = f_result_[k].data()[i];
}
// -------------------------------------------------------------------
// check if we are computing vy
if( vx.size() == 0 )
return true;
// ------------------------------------------------------------------
// This is a very dumb algorithm that over estimates which
// elements of the inverse are variables (which is not efficient).
bool var = false;
for(size_t i = 0; i < ny; i++)
var |= vx[1 + i];
for(size_t i = 0; i < ny; i++)
vy[i] = var;
return true;
}
/* %$$
$subhead reverse$$
$srccode%cpp% */
// reverse mode routine called by CppAD
virtual bool reverse(
// highest order Taylor coefficient that we are computing derivative of
size_t q ,
// forward mode Taylor coefficients for x variables
const CppAD::vector<double>& tx ,
// forward mode Taylor coefficients for y variables
const CppAD::vector<double>& ty ,
// upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
CppAD::vector<double>& px ,
// derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
const CppAD::vector<double>& py
)
{ size_t n_order = q + 1;
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
size_t ny = nr * nr;
# ifndef NDEBUG
size_t nx = 1 + ny;
# endif
//
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
assert( px.size() == tx.size() );
assert( py.size() == ty.size() );
// -------------------------------------------------------------------
// make sure f_arg_ is large enough
assert( f_arg_.size() == f_result_.size() );
// must have previous run forward with order >= n_order
assert( f_arg_.size() >= n_order );
// -------------------------------------------------------------------
// make sure r_arg_, r_result_ are large enough
assert( r_arg_.size() == r_result_.size() );
if( r_arg_.size() < n_order )
{ r_arg_.resize(n_order);
r_result_.resize(n_order);
//
for(size_t k = 0; k < n_order; k++)
{ r_arg_[k].resize( long(nr), long(nr) );
r_result_[k].resize( long(nr), long(nr) );
}
}
// -------------------------------------------------------------------
// unpack tx into f_arg_
for(size_t k = 0; k < n_order; k++)
{ // unpack arg values for this order
for(size_t i = 0; i < ny; i++)
f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ];
}
// -------------------------------------------------------------------
// unpack py into r_result_
for(size_t k = 0; k < n_order; k++)
{ for(size_t i = 0; i < ny; i++)
r_result_[k].data()[i] = py[ i * n_order + k ];
}
// -------------------------------------------------------------------
// initialize r_arg_ as zero
for(size_t k = 0; k < n_order; k++)
r_arg_[k] = matrix::Zero( long(nr), long(nr) );
// -------------------------------------------------------------------
// matrix reverse mode calculation
//
for(size_t k1 = n_order; k1 > 1; k1--)
{ size_t k = k1 - 1;
// bar{R}_0 = bar{R}_0 + bar{R}_k (A_0 R_k)^T
r_result_[0] +=
r_result_[k] * f_result_[k].transpose() * f_arg_[0].transpose();
//
for(size_t ell = 1; ell <= k; ell++)
{ // bar{A}_l = bar{A}_l - R_0^T bar{R}_k R_{k-l}^T
r_arg_[ell] -= f_result_[0].transpose()
* r_result_[k] * f_result_[k-ell].transpose();
// bar{R}_{k-l} = bar{R}_{k-1} - (R_0 A_l)^T bar{R}_k
r_result_[k-ell] -= f_arg_[ell].transpose()
* f_result_[0].transpose() * r_result_[k];
}
}
r_arg_[0] -=
f_result_[0].transpose() * r_result_[0] * f_result_[0].transpose();
// -------------------------------------------------------------------
// pack r_arg into px
for(size_t k = 0; k < n_order; k++)
{ for(size_t i = 0; i < ny; i++)
px[ (1 + i) * n_order + k ] = r_arg_[k].data()[i];
}
//
return true;
}
/* %$$
$head End Class Definition$$
$srccode%cpp% */
}; // End of atomic_eigen_mat_inv class
} // END_EMPTY_NAMESPACE
/* %$$
$$ $comment end nospell$$
$end
*/
# endif

658
build-config/cppad/include/cppad/example/atomic_two/eigen_mat_mul.hpp Normal file
View File

@@ -0,0 +1,658 @@
# ifndef CPPAD_EXAMPLE_ATOMIC_TWO_EIGEN_MAT_MUL_HPP
# define CPPAD_EXAMPLE_ATOMIC_TWO_EIGEN_MAT_MUL_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-19 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
/*
$begin atomic_two_eigen_mat_mul.hpp$$
$spell
Eigen
Taylor
nr
nc
$$
$section atomic_two Eigen Matrix Multiply Class$$
$head See Also$$
$cref atomic_three_mat_mul.hpp$$
$head Purpose$$
Construct an atomic operation that computes the matrix product,
$latex R = A \times \B{R}$$
for any positive integers $latex r$$, $latex m$$, $latex c$$,
and any $latex A \in \B{R}^{r \times m}$$,
$latex B \in \B{R}^{m \times c}$$.
$head Matrix Dimensions$$
This example puts the matrix dimensions in the atomic function arguments,
instead of the $cref/constructor/atomic_two_ctor/$$, so that they can
be different for different calls to the atomic function.
These dimensions are:
$table
$icode nr_left$$
$cnext number of rows in the left matrix; i.e, $latex r$$ $rend
$icode n_middle$$
$cnext rows in the left matrix and columns in right; i.e, $latex m$$ $rend
$icode nc_right$$
$cnext number of columns in the right matrix; i.e., $latex c$$
$tend
$head Theory$$
$subhead Forward$$
For $latex k = 0 , \ldots $$, the $th k$$ order Taylor coefficient
$latex R_k$$ is given by
$latex \[
R_k = \sum_{\ell = 0}^{k} A_\ell B_{k-\ell}
\] $$
$subhead Product of Two Matrices$$
Suppose $latex \bar{E}$$ is the derivative of the
scalar value function $latex s(E)$$ with respect to $latex E$$; i.e.,
$latex \[
\bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} }
\] $$
Also suppose that $latex t$$ is a scalar valued argument and
$latex \[
E(t) = C(t) D(t)
\] $$
It follows that
$latex \[
E'(t) = C'(t) D(t) + C(t) D'(t)
\] $$
$latex \[
(s \circ E)'(t)
=
\R{tr} [ \bar{E}^\R{T} E'(t) ]
\] $$
$latex \[
=
\R{tr} [ \bar{E}^\R{T} C'(t) D(t) ] +
\R{tr} [ \bar{E}^\R{T} C(t) D'(t) ]
\] $$
$latex \[
=
\R{tr} [ D(t) \bar{E}^\R{T} C'(t) ] +
\R{tr} [ \bar{E}^\R{T} C(t) D'(t) ]
\] $$
$latex \[
\bar{C} = \bar{E} D^\R{T} \W{,}
\bar{D} = C^\R{T} \bar{E}
\] $$
$subhead Reverse$$
Reverse mode eliminates $latex R_k$$ as follows:
for $latex \ell = 0, \ldots , k-1$$,
$latex \[
\bar{A}_\ell = \bar{A}_\ell + \bar{R}_k B_{k-\ell}^\R{T}
\] $$
$latex \[
\bar{B}_{k-\ell} = \bar{B}_{k-\ell} + A_\ell^\R{T} \bar{R}_k
\] $$
$nospell
$head Start Class Definition$$
$srccode%cpp% */
# include <cppad/cppad.hpp>
# include <Eigen/Core>
/* %$$
$head Public$$
$subhead Types$$
$srccode%cpp% */
namespace { // BEGIN_EMPTY_NAMESPACE
template <class Base>
class atomic_eigen_mat_mul : public CppAD::atomic_base<Base> {
public:
// -----------------------------------------------------------
// type of elements during calculation of derivatives
typedef Base scalar;
// type of elements during taping
typedef CppAD::AD<scalar> ad_scalar;
// type of matrix during calculation of derivatives
typedef Eigen::Matrix<
scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
// type of matrix during taping
typedef Eigen::Matrix<
ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
/* %$$
$subhead Constructor$$
$srccode%cpp% */
// constructor
atomic_eigen_mat_mul(void) : CppAD::atomic_base<Base>(
"atom_eigen_mat_mul" ,
CppAD::atomic_base<Base>::set_sparsity_enum
)
{ }
/* %$$
$subhead op$$
$srccode%cpp% */
// use atomic operation to multiply two AD matrices
ad_matrix op(
const ad_matrix& left ,
const ad_matrix& right )
{ size_t nr_left = size_t( left.rows() );
size_t n_middle = size_t( left.cols() );
size_t nc_right = size_t( right.cols() );
assert( n_middle == size_t( right.rows() ) );
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
size_t n_left = nr_left * n_middle;
size_t n_right = n_middle * nc_right;
size_t n_result = nr_left * nc_right;
//
assert( 3 + n_left + n_right == nx );
assert( n_result == ny );
// -----------------------------------------------------------------
// packed version of left and right
CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
//
packed_arg[0] = ad_scalar( nr_left );
packed_arg[1] = ad_scalar( n_middle );
packed_arg[2] = ad_scalar( nc_right );
for(size_t i = 0; i < n_left; i++)
packed_arg[3 + i] = left.data()[i];
for(size_t i = 0; i < n_right; i++)
packed_arg[ 3 + n_left + i ] = right.data()[i];
// ------------------------------------------------------------------
// Packed version of result = left * right.
// This as an atomic_base funciton call that CppAD uses
// to store the atomic operation on the tape.
CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
(*this)(packed_arg, packed_result);
// ------------------------------------------------------------------
// unpack result matrix
ad_matrix result(nr_left, nc_right);
for(size_t i = 0; i < n_result; i++)
result.data()[i] = packed_result[ i ];
//
return result;
}
/* %$$
$head Private$$
$subhead Variables$$
$srccode%cpp% */
private:
// -------------------------------------------------------------
// one forward mode vector of matrices for left, right, and result
CppAD::vector<matrix> f_left_, f_right_, f_result_;
// one reverse mode vector of matrices for left, right, and result
CppAD::vector<matrix> r_left_, r_right_, r_result_;
// -------------------------------------------------------------
/* %$$
$subhead forward$$
$srccode%cpp% */
// forward mode routine called by CppAD
virtual bool forward(
// lowest order Taylor coefficient we are evaluating
size_t p ,
// highest order Taylor coefficient we are evaluating
size_t q ,
// which components of x are variables
const CppAD::vector<bool>& vx ,
// which components of y are variables
CppAD::vector<bool>& vy ,
// tx [ 3 + j * (q+1) + k ] is x_j^k
const CppAD::vector<scalar>& tx ,
// ty [ i * (q+1) + k ] is y_i^k
CppAD::vector<scalar>& ty
)
{ size_t n_order = q + 1;
size_t nr_left = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
size_t n_middle = size_t( CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
size_t nc_right = size_t( CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
# ifndef NDEBUG
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
//
assert( vx.size() == 0 || nx == vx.size() );
assert( vx.size() == 0 || ny == vy.size() );
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
//
size_t n_left = nr_left * n_middle;
size_t n_right = n_middle * nc_right;
size_t n_result = nr_left * nc_right;
assert( 3 + n_left + n_right == nx );
assert( n_result == ny );
//
// -------------------------------------------------------------------
// make sure f_left_, f_right_, and f_result_ are large enough
assert( f_left_.size() == f_right_.size() );
assert( f_left_.size() == f_result_.size() );
if( f_left_.size() < n_order )
{ f_left_.resize(n_order);
f_right_.resize(n_order);
f_result_.resize(n_order);
//
for(size_t k = 0; k < n_order; k++)
{ f_left_[k].resize( long(nr_left), long(n_middle) );
f_right_[k].resize( long(n_middle), long(nc_right) );
f_result_[k].resize( long(nr_left), long(nc_right) );
}
}
// -------------------------------------------------------------------
// unpack tx into f_left and f_right
for(size_t k = 0; k < n_order; k++)
{ // unpack left values for this order
for(size_t i = 0; i < n_left; i++)
f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
//
// unpack right values for this order
for(size_t i = 0; i < n_right; i++)
f_right_[k].data()[i] = tx[ ( 3 + n_left + i) * n_order + k ];
}
// -------------------------------------------------------------------
// result for each order
// (we could avoid recalculting f_result_[k] for k=0,...,p-1)
for(size_t k = 0; k < n_order; k++)
{ // result[k] = sum_ell left[ell] * right[k-ell]
f_result_[k] = matrix::Zero( long(nr_left), long(nc_right) );
for(size_t ell = 0; ell <= k; ell++)
f_result_[k] += f_left_[ell] * f_right_[k-ell];
}
// -------------------------------------------------------------------
// pack result_ into ty
for(size_t k = 0; k < n_order; k++)
{ for(size_t i = 0; i < n_result; i++)
ty[ i * n_order + k ] = f_result_[k].data()[i];
}
// ------------------------------------------------------------------
// check if we are computing vy
if( vx.size() == 0 )
return true;
// ------------------------------------------------------------------
// compute variable information for y; i.e., vy
// (note that the constant zero times a variable is a constant)
scalar zero(0.0);
assert( n_order == 1 );
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ bool var = false;
for(size_t ell = 0; ell < n_middle; ell++)
{ // left information
size_t index = 3 + i * n_middle + ell;
bool var_left = vx[index];
bool nz_left = var_left |
(f_left_[0]( long(i), long(ell) ) != zero);
// right information
index = 3 + n_left + ell * nc_right + j;
bool var_right = vx[index];
bool nz_right = var_right |
(f_right_[0]( long(ell), long(j) ) != zero);
// effect of result
var |= var_left & nz_right;
var |= nz_left & var_right;
}
size_t index = i * nc_right + j;
vy[index] = var;
}
}
return true;
}
/* %$$
$subhead reverse$$
$srccode%cpp% */
// reverse mode routine called by CppAD
virtual bool reverse(
// highest order Taylor coefficient that we are computing derivative of
size_t q ,
// forward mode Taylor coefficients for x variables
const CppAD::vector<double>& tx ,
// forward mode Taylor coefficients for y variables
const CppAD::vector<double>& ty ,
// upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
CppAD::vector<double>& px ,
// derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
const CppAD::vector<double>& py
)
{ size_t n_order = q + 1;
size_t nr_left = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
size_t n_middle = size_t( CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
size_t nc_right = size_t( CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
# ifndef NDEBUG
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
//
assert( nx * n_order == tx.size() );
assert( ny * n_order == ty.size() );
assert( px.size() == tx.size() );
assert( py.size() == ty.size() );
//
size_t n_left = nr_left * n_middle;
size_t n_right = n_middle * nc_right;
size_t n_result = nr_left * nc_right;
assert( 3 + n_left + n_right == nx );
assert( n_result == ny );
// -------------------------------------------------------------------
// make sure f_left_, f_right_ are large enough
assert( f_left_.size() == f_right_.size() );
assert( f_left_.size() == f_result_.size() );
// must have previous run forward with order >= n_order
assert( f_left_.size() >= n_order );
// -------------------------------------------------------------------
// make sure r_left_, r_right_, and r_result_ are large enough
assert( r_left_.size() == r_right_.size() );
assert( r_left_.size() == r_result_.size() );
if( r_left_.size() < n_order )
{ r_left_.resize(n_order);
r_right_.resize(n_order);
r_result_.resize(n_order);
//
for(size_t k = 0; k < n_order; k++)
{ r_left_[k].resize( long(nr_left), long(n_middle) );
r_right_[k].resize( long(n_middle), long(nc_right) );
r_result_[k].resize( long(nr_left), long(nc_right) );
}
}
// -------------------------------------------------------------------
// unpack tx into f_left and f_right
for(size_t k = 0; k < n_order; k++)
{ // unpack left values for this order
for(size_t i = 0; i < n_left; i++)
f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
//
// unpack right values for this order
for(size_t i = 0; i < n_right; i++)
f_right_[k].data()[i] = tx[ (3 + n_left + i) * n_order + k ];
}
// -------------------------------------------------------------------
// unpack py into r_result_
for(size_t k = 0; k < n_order; k++)
{ for(size_t i = 0; i < n_result; i++)
r_result_[k].data()[i] = py[ i * n_order + k ];
}
// -------------------------------------------------------------------
// initialize r_left_ and r_right_ as zero
for(size_t k = 0; k < n_order; k++)
{ r_left_[k] = matrix::Zero( long(nr_left), long(n_middle) );
r_right_[k] = matrix::Zero( long(n_middle), long(nc_right) );
}
// -------------------------------------------------------------------
// matrix reverse mode calculation
for(size_t k1 = n_order; k1 > 0; k1--)
{ size_t k = k1 - 1;
for(size_t ell = 0; ell <= k; ell++)
{ // nr x nm = nr x nc * nc * nm
r_left_[ell] += r_result_[k] * f_right_[k-ell].transpose();
// nm x nc = nm x nr * nr * nc
r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
}
}
// -------------------------------------------------------------------
// pack r_left and r_right int px
for(size_t k = 0; k < n_order; k++)
{ // dimensions are integer constants
px[ 0 * n_order + k ] = 0.0;
px[ 1 * n_order + k ] = 0.0;
px[ 2 * n_order + k ] = 0.0;
//
// pack left values for this order
for(size_t i = 0; i < n_left; i++)
px[ (3 + i) * n_order + k ] = r_left_[k].data()[i];
//
// pack right values for this order
for(size_t i = 0; i < n_right; i++)
px[ (3 + i + n_left) * n_order + k] = r_right_[k].data()[i];
}
//
return true;
}
/* %$$
$subhead for_sparse_jac$$
$srccode%cpp% */
// forward Jacobian sparsity routine called by CppAD
virtual bool for_sparse_jac(
// number of columns in the matrix R
size_t q ,
// sparsity pattern for the matrix R
const CppAD::vector< std::set<size_t> >& r ,
// sparsity pattern for the matrix S = f'(x) * R
CppAD::vector< std::set<size_t> >& s ,
const CppAD::vector<Base>& x )
{
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
# ifndef NDEBUG
size_t nx = 3 + (nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
//
assert( nx == r.size() );
assert( ny == s.size() );
//
size_t n_left = nr_left * n_middle;
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ // pack index for entry (i, j) in result
size_t i_result = i * nc_right + j;
s[i_result].clear();
for(size_t ell = 0; ell < n_middle; ell++)
{ // pack index for entry (i, ell) in left
size_t i_left = 3 + i * n_middle + ell;
// pack index for entry (ell, j) in right
size_t i_right = 3 + n_left + ell * nc_right + j;
// check if result of for this product is alwasy zero
// note that x is nan for commponents that are variables
bool zero = x[i_left] == Base(0.0) || x[i_right] == Base(0);
if( ! zero )
{ s[i_result] =
CppAD::set_union(s[i_result], r[i_left] );
s[i_result] =
CppAD::set_union(s[i_result], r[i_right] );
}
}
}
}
return true;
}
/* %$$
$subhead rev_sparse_jac$$
$srccode%cpp% */
// reverse Jacobian sparsity routine called by CppAD
virtual bool rev_sparse_jac(
// number of columns in the matrix R^T
size_t q ,
// sparsity pattern for the matrix R^T
const CppAD::vector< std::set<size_t> >& rt ,
// sparsoity pattern for the matrix S^T = f'(x)^T * R^T
CppAD::vector< std::set<size_t> >& st ,
const CppAD::vector<Base>& x )
{
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
size_t nx = 3 + (nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
//
assert( nx == st.size() );
assert( ny == rt.size() );
//
// initialize S^T as empty
for(size_t i = 0; i < nx; i++)
st[i].clear();
// sparsity for S(x)^T = f'(x)^T * R^T
size_t n_left = nr_left * n_middle;
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ // pack index for entry (i, j) in result
size_t i_result = i * nc_right + j;
st[i_result].clear();
for(size_t ell = 0; ell < n_middle; ell++)
{ // pack index for entry (i, ell) in left
size_t i_left = 3 + i * n_middle + ell;
// pack index for entry (ell, j) in right
size_t i_right = 3 + n_left + ell * nc_right + j;
//
st[i_left] = CppAD::set_union(st[i_left], rt[i_result]);
st[i_right] = CppAD::set_union(st[i_right], rt[i_result]);
}
}
}
return true;
}
/* %$$
$subhead for_sparse_hes$$
$srccode%cpp% */
virtual bool for_sparse_hes(
// which components of x are variables for this call
const CppAD::vector<bool>& vx,
// sparsity pattern for the diagonal of R
const CppAD::vector<bool>& r ,
// sparsity pattern for the vector S
const CppAD::vector<bool>& s ,
// sparsity patternfor the Hessian H(x)
CppAD::vector< std::set<size_t> >& h ,
const CppAD::vector<Base>& x )
{
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
size_t nx = 3 + (nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
//
assert( vx.size() == nx );
assert( r.size() == nx );
assert( s.size() == ny );
assert( h.size() == nx );
//
// initilize h as empty
for(size_t i = 0; i < nx; i++)
h[i].clear();
//
size_t n_left = nr_left * n_middle;
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ // pack index for entry (i, j) in result
size_t i_result = i * nc_right + j;
if( s[i_result] )
{ for(size_t ell = 0; ell < n_middle; ell++)
{ // pack index for entry (i, ell) in left
size_t i_left = 3 + i * n_middle + ell;
// pack index for entry (ell, j) in right
size_t i_right = 3 + n_left + ell * nc_right + j;
if( r[i_left] & r[i_right] )
{ h[i_left].insert(i_right);
h[i_right].insert(i_left);
}
}
}
}
}
return true;
}
/* %$$
$subhead rev_sparse_hes$$
$srccode%cpp% */
// reverse Hessian sparsity routine called by CppAD
virtual bool rev_sparse_hes(
// which components of x are variables for this call
const CppAD::vector<bool>& vx,
// sparsity pattern for S(x) = g'[f(x)]
const CppAD::vector<bool>& s ,
// sparsity pattern for d/dx g[f(x)] = S(x) * f'(x)
CppAD::vector<bool>& t ,
// number of columns in R, U(x), and V(x)
size_t q ,
// sparsity pattern for R
const CppAD::vector< std::set<size_t> >& r ,
// sparsity pattern for U(x) = g^{(2)} [ f(x) ] * f'(x) * R
const CppAD::vector< std::set<size_t> >& u ,
// sparsity pattern for
// V(x) = f'(x)^T * U(x) + sum_{i=0}^{m-1} S_i(x) f_i^{(2)} (x) * R
CppAD::vector< std::set<size_t> >& v ,
// parameters as integers
const CppAD::vector<Base>& x )
{
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
size_t nx = 3 + (nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
//
assert( vx.size() == nx );
assert( s.size() == ny );
assert( t.size() == nx );
assert( r.size() == nx );
assert( v.size() == nx );
//
// initilaize return sparsity patterns as false
for(size_t j = 0; j < nx; j++)
{ t[j] = false;
v[j].clear();
}
//
size_t n_left = nr_left * n_middle;
for(size_t i = 0; i < nr_left; i++)
{ for(size_t j = 0; j < nc_right; j++)
{ // pack index for entry (i, j) in result
size_t i_result = i * nc_right + j;
for(size_t ell = 0; ell < n_middle; ell++)
{ // pack index for entry (i, ell) in left
size_t i_left = 3 + i * n_middle + ell;
// pack index for entry (ell, j) in right
size_t i_right = 3 + n_left + ell * nc_right + j;
//
// back propagate T(x) = S(x) * f'(x).
t[i_left] |= bool( s[i_result] );
t[i_right] |= bool( s[i_result] );
//
// V(x) = f'(x)^T * U(x) + sum_i S_i(x) * f_i''(x) * R
// U(x) = g''[ f(x) ] * f'(x) * R
// S_i(x) = g_i'[ f(x) ]
//
// back propagate f'(x)^T * U(x)
v[i_left] = CppAD::set_union(v[i_left], u[i_result] );
v[i_right] = CppAD::set_union(v[i_right], u[i_result] );
//
// back propagate S_i(x) * f_i''(x) * R
// (here is where we use vx to check for cross terms)
if( s[i_result] & vx[i_left] & vx[i_right] )
{ v[i_left] = CppAD::set_union(v[i_left], r[i_right] );
v[i_right] = CppAD::set_union(v[i_right], r[i_left] );
}
}
}
}
return true;
}
/* %$$
$head End Class Definition$$
$srccode%cpp% */
}; // End of atomic_eigen_mat_mul class
} // END_EMPTY_NAMESPACE
/* %$$
$$ $comment end nospell$$
$end
*/
# endif

359
build-config/cppad/include/cppad/example/base_adolc.hpp Normal file
View File

@@ -0,0 +1,359 @@
# ifndef CPPAD_EXAMPLE_BASE_ADOLC_HPP
# define CPPAD_EXAMPLE_BASE_ADOLC_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
/*
$begin base_adolc.hpp$$
$spell
stringstream
struct
string
setprecision
str
valgrind
azmul
expm1
atanh
acosh
asinh
erf
erfc
ifndef
define
endif
Rel
codassign
eps
std
abs_geq
fabs
cppad.hpp
undef
Lt
Le
Eq
Ge
Gt
namespace
cassert
condassign
hpp
bool
const
Adolc
adouble
CondExpOp
inline
enum
CppAD
pow
acos
asin
atan
cos
cosh
exp
sqrt
atrig
$$
$section Enable use of AD<Base> where Base is Adolc's adouble Type$$
$head Syntax$$
$codei%# include <cppad/example/base_adolc.hpp>
%$$
$children%
example/general/mul_level_adolc.cpp
%$$
$head Example$$
The file $cref mul_level_adolc.cpp$$ contains an example use of
Adolc's $code adouble$$ type for a CppAD $icode Base$$ type.
The file $cref mul_level_adolc_ode.cpp$$ contains a more realistic
(and complex) example.
$head Include Files$$
This file $code base_adolc.hpp$$ requires $code adouble$$ to be defined.
In addition, it is included before $code <cppad/cppad.hpp>$$,
but it needs to include parts of CppAD that are used by this file.
This is done with the following include commands:
$srccode%cpp% */
# include <adolc/adolc.h>
# include <cppad/base_require.hpp>
/* %$$
$head CondExpOp$$
The type $code adouble$$ supports a conditional assignment function
with the syntax
$codei%
condassign(%a%, %b%, %c%, %d%)
%$$
which evaluates to
$codei%
%a% = (%b% > 0) ? %c% : %d%;
%$$
This enables one to include conditionals in the recording of
$code adouble$$ operations and later evaluation for different
values of the independent variables
(in the same spirit as the CppAD $cref CondExp$$ function).
$srccode%cpp% */
namespace CppAD {
inline adouble CondExpOp(
enum CppAD::CompareOp cop ,
const adouble &left ,
const adouble &right ,
const adouble &trueCase ,
const adouble &falseCase )
{ adouble result;
switch( cop )
{
case CompareLt: // left < right
condassign(result, right - left, trueCase, falseCase);
break;
case CompareLe: // left <= right
condassign(result, left - right, falseCase, trueCase);
break;
case CompareEq: // left == right
condassign(result, left - right, falseCase, trueCase);
condassign(result, right - left, falseCase, result);
break;
case CompareGe: // left >= right
condassign(result, right - left, falseCase, trueCase);
break;
case CompareGt: // left > right
condassign(result, left - right, trueCase, falseCase);
break;
default:
CppAD::ErrorHandler::Call(
true , __LINE__ , __FILE__ ,
"CppAD::CondExp",
"Error: for unknown reason."
);
result = trueCase;
}
return result;
}
}
/* %$$
$head CondExpRel$$
The $cref/CPPAD_COND_EXP_REL/base_cond_exp/CondExpRel/$$ macro invocation
$srccode%cpp% */
namespace CppAD {
CPPAD_COND_EXP_REL(adouble)
}
/* %$$
$head EqualOpSeq$$
The Adolc user interface does not specify a way to determine if
two $code adouble$$ variables correspond to the same operations sequence.
Make $code EqualOpSeq$$ an error if it gets used:
$srccode%cpp% */
namespace CppAD {
inline bool EqualOpSeq(const adouble &x, const adouble &y)
{ CppAD::ErrorHandler::Call(
true , __LINE__ , __FILE__ ,
"CppAD::EqualOpSeq(x, y)",
"Error: adouble does not support EqualOpSeq."
);
return false;
}
}
/* %$$
$head Identical$$
The Adolc user interface does not specify a way to determine if an
$code adouble$$ depends on the independent variables.
To be safe (but slow) return $code false$$ in all the cases below.
$srccode%cpp% */
namespace CppAD {
inline bool IdenticalCon(const adouble &x)
{ return false; }
inline bool IdenticalZero(const adouble &x)
{ return false; }
inline bool IdenticalOne(const adouble &x)
{ return false; }
inline bool IdenticalEqualCon(const adouble &x, const adouble &y)
{ return false; }
}
/* %$$
$head Integer$$
$srccode%cpp% */
inline int Integer(const adouble &x)
{ return static_cast<int>( x.getValue() ); }
/* %$$
$head azmul$$
$srccode%cpp% */
namespace CppAD {
CPPAD_AZMUL( adouble )
}
/* %$$
$head Ordered$$
$srccode%cpp% */
namespace CppAD {
inline bool GreaterThanZero(const adouble &x)
{ return (x > 0); }
inline bool GreaterThanOrZero(const adouble &x)
{ return (x >= 0); }
inline bool LessThanZero(const adouble &x)
{ return (x < 0); }
inline bool LessThanOrZero(const adouble &x)
{ return (x <= 0); }
inline bool abs_geq(const adouble& x, const adouble& y)
{ return fabs(x) >= fabs(y); }
}
/* %$$
$head Unary Standard Math$$
The following $cref/required/base_require/$$ functions
are defined by the Adolc package for the $code adouble$$ base case:
$pre
$$
$code acos$$,
$code acosh$$,
$code asin$$,
$code asinh$$,
$code atan$$,
$code atanh$$,
$code cos$$,
$code cosh$$,
$code erf$$,
$code exp$$,
$code fabs$$,
$code log$$,
$code sin$$,
$code sinh$$,
$code sqrt$$,
$code tan$$.
$head erfc$$
If you provide $code --enable-atrig-erf$$ on the configure command line,
the adolc package supports all the c++11 math functions except
$code erfc$$, $code expm1$$, and $code log1p$$.
For the reason, we make using $code erfc$$ an error:
$srccode%cpp% */
namespace CppAD {
# define CPPAD_BASE_ADOLC_NO_SUPPORT(fun) \
inline adouble fun(const adouble& x) \
{ CPPAD_ASSERT_KNOWN( \
false, \
#fun ": adolc does not support this function" \
); \
return 0.0; \
}
CPPAD_BASE_ADOLC_NO_SUPPORT(erfc)
CPPAD_BASE_ADOLC_NO_SUPPORT(expm1)
CPPAD_BASE_ADOLC_NO_SUPPORT(log1p)
# undef CPPAD_BASE_ADOLC_NO_SUPPORT
}
/* %$$
$head sign$$
This $cref/required/base_require/$$ function is defined using the
$code codassign$$ function so that its $code adouble$$ operation sequence
does not depend on the value of $icode x$$.
$srccode%cpp% */
namespace CppAD {
inline adouble sign(const adouble& x)
{ adouble s_plus, s_minus, half(.5);
// set s_plus to sign(x)/2, except for case x == 0, s_plus = -.5
condassign(s_plus, +x, -half, +half);
// set s_minus to -sign(x)/2, except for case x == 0, s_minus = -.5
condassign(s_minus, -x, -half, +half);
// set s to sign(x)
return s_plus - s_minus;
}
}
/* %$$
$head abs$$
This $cref/required/base_require/$$ function uses the adolc $code fabs$$
function:
$srccode%cpp% */
namespace CppAD {
inline adouble abs(const adouble& x)
{ return fabs(x); }
}
/* %$$
$head pow$$
This $cref/required/base_require/$$ function
is defined by the Adolc package for the $code adouble$$ base case.
$head numeric_limits$$
The following defines the CppAD $cref numeric_limits$$
for the type $code adouble$$:
$srccode%cpp% */
namespace CppAD {
CPPAD_NUMERIC_LIMITS(double, adouble)
}
/* %$$
$head to_string$$
The following defines the CppAD $cref to_string$$ function
for the type $code adouble$$:
$srccode%cpp% */
namespace CppAD {
template <> struct to_string_struct<adouble>
{ std::string operator()(const adouble& x)
{ std::stringstream os;
int n_digits = 1 + std::numeric_limits<double>::digits10;
os << std::setprecision(n_digits);
os << x.value();
return os.str();
}
};
}
/* %$$
$head hash_code$$
It appears that an $code adouble$$ object can have fields
that are not initialized.
This results in a $code valgrind$$ error when these fields are used by the
$cref/default/base_hash/Default/$$ hashing function.
For this reason, the $code adouble$$ class overrides the default definition.
$srccode|cpp| */
namespace CppAD {
inline unsigned short hash_code(const adouble& x)
{ unsigned short code = 0;
double value = x.value();
if( value == 0.0 )
return code;
double log_x = std::log( fabs( value ) );
// assume log( std::numeric_limits<double>::max() ) is near 700
code = static_cast<unsigned short>(
(CPPAD_HASH_TABLE_SIZE / 700 + 1) * log_x
);
code = code % CPPAD_HASH_TABLE_SIZE;
return code;
}
}
/* |$$
Note that after the hash codes match, the
$cref/Identical/base_adolc.hpp/Identical/$$ function will be used
to make sure two values are the same and one can replace the other.
A more sophisticated implementation of the $code Identical$$ function
would detect which $code adouble$$ values depend on the
$code adouble$$ independent variables (and hence can change).
$end
*/
# endif

62
build-config/cppad/include/cppad/example/code_gen_fun.hpp Normal file
View File

@@ -0,0 +1,62 @@
# ifndef CPPAD_EXAMPLE_CODE_GEN_FUN_HPP
# define CPPAD_EXAMPLE_CODE_GEN_FUN_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
-------------------------------------------------------------------------- */
// BEGIN C++
# include <cppad/cg/cppadcg.hpp>
class code_gen_fun {
public:
// type of evaluation for Jacobians (possibly Hessians in the future)
enum evaluation_enum { none_enum, dense_enum, sparse_enum };
private:
// dynamic_lib_
std::unique_ptr< CppAD::cg::DynamicLib<double> > dynamic_lib_;
//
// model_ (contains a reference to dynamic_lib_)
std::unique_ptr< CppAD::cg::GenericModel<double> > model_;
//
public:
// -----------------------------------------------------------------------
// constructors
// -----------------------------------------------------------------------
// fun_name()
code_gen_fun(void);
//
// fun_name( file_name )
code_gen_fun(const std::string& file_name);
//
// fun_name(file_name, cg_fun, eval_jac)
code_gen_fun(
const std::string& file_name ,
CppAD::ADFun< CppAD::cg::CG<double> >& cg_fun ,
evaluation_enum eval_jac = none_enum
);
// -----------------------------------------------------------------------
// operations
// -----------------------------------------------------------------------
// swap(other_fun)
void swap(code_gen_fun& other_fun);
//
// y = fun_name(x)
CppAD::vector<double> operator()(const CppAD::vector<double> & x);
//
// J = fun_name.jacobian(x)
CppAD::vector<double> jacobian(const CppAD::vector<double> & x);
//
// Jrcv = fun_name.sparse_jacobian(x)
CppAD::sparse_rcv< CppAD::vector<size_t>, CppAD::vector<double> >
sparse_jacobian(const CppAD::vector<double>& x);
};
// END C++
# endif

195
build-config/cppad/include/cppad/example/cppad_eigen.hpp Normal file
View File

@@ -0,0 +1,195 @@
# ifndef CPPAD_EXAMPLE_CPPAD_EIGEN_HPP
# define CPPAD_EXAMPLE_CPPAD_EIGEN_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
// cppad.hpp gets included at the end
# define EIGEN_MATRIXBASE_PLUGIN <cppad/example/eigen_plugin.hpp>
# include <Eigen/Core>
/*
$begin cppad_eigen.hpp$$
$spell
impl
typename
Real Real
inline
neg
eps
atan
Num
acos
asin
CppAD
std::numeric
enum
Mul
Eigen
cppad.hpp
namespace
struct
typedef
const
imag
sqrt
exp
cos
$$
$section Enable Use of Eigen Linear Algebra Package with CppAD$$
$head Syntax$$
$codei%# include <cppad/example/cppad_eigen.hpp>
%$$
$children%
include/cppad/example/eigen_plugin.hpp%
example/general/eigen_array.cpp%
example/general/eigen_det.cpp
%$$
$head Purpose$$
Enables the use of the $cref eigen$$
linear algebra package with the type $icode%AD<%Base%>%$$; see
$href%
https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html%
custom scalar types
%$$.
$head Example$$
The files $cref eigen_array.cpp$$ and $cref eigen_det.cpp$$
contain an example and test of this include file.
They return true if they succeed and false otherwise.
$head Include Files$$
The file $code <Eigen/Core>$$ is included before
these definitions and $code <cppad/cppad.hpp>$$ is included after.
$head CppAD Declarations$$
First declare some items that are defined by cppad.hpp:
$srccode%cpp% */
namespace CppAD {
// AD<Base>
template <class Base> class AD;
// numeric_limits<Float>
template <class Float> class numeric_limits;
}
/* %$$
$head Eigen NumTraits$$
Eigen needs the following definitions to work properly
with $codei%AD<%Base%>%$$ scalars:
$srccode%cpp% */
namespace Eigen {
template <class Base> struct NumTraits< CppAD::AD<Base> >
{ // type that corresponds to the real part of an AD<Base> value
typedef CppAD::AD<Base> Real;
// type for AD<Base> operations that result in non-integer values
typedef CppAD::AD<Base> NonInteger;
// type to use for numeric literals such as "2" or "0.5".
typedef CppAD::AD<Base> Literal;
// type for nested value inside an AD<Base> expression tree
typedef CppAD::AD<Base> Nested;
enum {
// does not support complex Base types
IsComplex = 0 ,
// does not support integer Base types
IsInteger = 0 ,
// only support signed Base types
IsSigned = 1 ,
// must initialize an AD<Base> object
RequireInitialization = 1 ,
// computational cost of the corresponding operations
ReadCost = 1 ,
AddCost = 2 ,
MulCost = 2
};
// machine epsilon with type of real part of x
// (use assumption that Base is not complex)
static CppAD::AD<Base> epsilon(void)
{ return CppAD::numeric_limits< CppAD::AD<Base> >::epsilon(); }
// relaxed version of machine epsilon for comparison of different
// operations that should result in the same value
static CppAD::AD<Base> dummy_precision(void)
{ return 100. *
CppAD::numeric_limits< CppAD::AD<Base> >::epsilon();
}
// minimum normalized positive value
static CppAD::AD<Base> lowest(void)
{ return CppAD::numeric_limits< CppAD::AD<Base> >::min(); }
// maximum finite value
static CppAD::AD<Base> highest(void)
{ return CppAD::numeric_limits< CppAD::AD<Base> >::max(); }
// number of decimal digits that can be represented without change.
static int digits10(void)
{ return CppAD::numeric_limits< CppAD::AD<Base> >::digits10; }
};
}
/* %$$
$head CppAD Namespace$$
Eigen also needs the following definitions to work properly
with $codei%AD<%Base%>%$$ scalars:
$srccode%cpp% */
namespace CppAD {
// functions that return references
template <class Base> const AD<Base>& conj(const AD<Base>& x)
{ return x; }
template <class Base> const AD<Base>& real(const AD<Base>& x)
{ return x; }
// functions that return values (note abs is defined by cppad.hpp)
template <class Base> AD<Base> imag(const AD<Base>& x)
{ return CppAD::AD<Base>(0.); }
template <class Base> AD<Base> abs2(const AD<Base>& x)
{ return x * x; }
}
/* %$$
$head eigen_vector$$
The class $code CppAD::eigen_vector$$ is a wrapper for Eigen column vectors
so that they are $cref/simple vectors/SimpleVector/$$.
To be specific, it converts $code Eigen::Index$$ arguments and
return values to $code size_t$$.
$srccode%cpp% */
namespace CppAD {
template <class Scalar>
class eigen_vector : public Eigen::Matrix<Scalar, Eigen::Dynamic, 1> {
private:
// base_class
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> base_class;
public:
// constructor
eigen_vector(size_t n) : base_class( Eigen::Index(n) )
{ }
eigen_vector(void) : base_class()
{ }
// operator[]
Scalar& operator[](size_t i)
{ return base_class::operator[]( Eigen::Index(i) ); }
const Scalar& operator[](size_t i) const
{ return base_class::operator[]( Eigen::Index(i) ); }
// size
size_t size(void) const
{ return size_t( base_class::size() ); }
// resize
void resize(size_t n)
{ base_class::resize( Eigen::Index(n) ); }
};
}
//
# include <cppad/cppad.hpp>
/* %$$
$end
*/
# endif

28
build-config/cppad/include/cppad/example/eigen_plugin.hpp Normal file
View File

@@ -0,0 +1,28 @@
# ifndef CPPAD_EXAMPLE_EIGEN_PLUGIN_HPP
# define CPPAD_EXAMPLE_EIGEN_PLUGIN_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
CppAD is distributed under the terms of the
Eclipse Public License Version 2.0.
This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */
/*$
$begin eigen_plugin.hpp$$
$spell
eigen_plugin.hpp
typedef
$$
$section Source Code for eigen_plugin.hpp$$
$srccode%cpp% */
// Declaration needed, before eigen-3.3.3, so Eigen vector is a simple vector
typedef Scalar value_type;
/* %$$
$end
*/
# endif