Prev | Next | atomic_one |
CPPAD_USER_ATOMIC
has been deprecated.
Use atomic_three
instead.
CPPAD_USER_ATOMIC(afun, Tvector, Base,
forward, reverse, for_jac_sparse, rev_jac_sparse, rev_hes_sparse
)
afun(id, ax, ay)
ok = forward(id, k, n, m, vx, vy, tx, ty)
ok = reverse(id, k, n, m, tx, ty, px, py)
ok = for_jac_sparse(id, n, m, q, r, s)
ok = rev_jac_sparse(id, n, m, q, r, s)
ok = rev_hes_sparse(id, n, m, q, r, s, t, u, v)
user_atomic<Base>::clear()
AD<Base>
atomic_base
operations
and letting CppAD do the rest.
In this case, CPPAD_USER_ATOMIC
can be used
add the user code for @(@
f(x)
@)@, and its derivatives,
to the set of
AD<Base>
atomic operations.
Another possible purpose is to reduce the size of the tape.
forward
the routine,
for the case
k = 0
, must be implemented.
Functions with the correct prototype,
that just return false
,
can be used for the other cases
(unless they are required by your calculations).
For example, you need not implement
forward
for the case
k == 2
until you require
forward mode calculation of second derivatives.
CPPAD_USER_ATOMIC(afun, Tvector, Base,
forward, reverse, for_jac_sparse, rev_jac_sparse, rev_hes_sparse
)
defines the
AD<Base>
routine
afun
.
This macro can be placed within a namespace
(not the CppAD
namespace)
but must be outside of any routine.
Tvector
must be a
simple vector template class
.
It determines the type of vectors used as arguments to the routine
afun
.
Base
specifies the
base type
corresponding to
AD<Base>
operation sequences.
Calling the routine
afun
will add the operator defined
by this macro to an
AD<Base>
operation sequence.
ok
has prototype
bool ok
If it is true
, the corresponding evaluation succeeded,
otherwise it failed.
id
has prototype
size_t id
Its value in all other calls is the same as in the corresponding
call to
afun
.
It can be used to store and retrieve extra information about
a specific call to
afun
.
k
has prototype
size_t k
The value
k
is the order of the Taylor coefficient that
we are evaluating (forward
)
or taking the derivative of (reverse
).
n
has prototype
size_t n
It is the size of the vector
ax
in the corresponding call to
afun(id, ax, ay)
; i.e.,
the dimension of the domain space for @(@
y = f(x)
@)@.
m
has prototype
size_t m
It is the size of the vector
ay
in the corresponding call to
afun(id, ax, ay)
; i.e.,
the dimension of the range space for @(@
y = f(x)
@)@.
tx
has prototype
const CppAD::vector<Base>& tx
and
tx.size() >= (k + 1) * n
.
For @(@
j = 0 , \ldots , n-1
@)@ and @(@
\ell = 0 , \ldots , k
@)@,
we use the Taylor coefficient notation
@[@
\begin{array}{rcl}
x_j^\ell & = & tx [ j * ( k + 1 ) + \ell ]
\\
X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^k t^k
\end{array}
@]@
If
tx.size() > (k + 1) * n
,
the other components of
tx
are not specified and should not be used.
Note that superscripts represent an index for @(@
x_j^\ell
@)@
and an exponent for @(@
t^\ell
@)@.
Also note that the Taylor coefficients for @(@
X(t)
@)@ correspond
to the derivatives of @(@
X(t)
@)@ at @(@
t = 0
@)@ in the following way:
@[@
x_j^\ell = \frac{1}{ \ell ! } X_j^{(\ell)} (0)
@]@
ty
has prototype
CppAD::vector<Base>& ty
while in calls to reverse
it has prototype
const CppAD::vector<Base>& ty
For all calls,
tx.size() >= (k + 1) * m
.
For @(@
i = 0 , \ldots , m-1
@)@ and @(@
\ell = 0 , \ldots , k
@)@,
we use the Taylor coefficient notation
@[@
\begin{array}{rcl}
y_i^\ell & = & ty [ i * ( k + 1 ) + \ell ]
\\
Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^k t^k + o ( t^k )
\end{array}
@]@
where @(@
o( t^k ) / t^k \rightarrow 0
@)@ as @(@
t \rightarrow 0
@)@.
If
ty.size() > (k + 1) * m
,
the other components of
ty
are not specified and should not be used.
Note that superscripts represent an index for @(@
y_j^\ell
@)@
and an exponent for @(@
t^\ell
@)@.
Also note that the Taylor coefficients for @(@
Y(t)
@)@ correspond
to the derivatives of @(@
Y(t)
@)@ at @(@
t = 0
@)@ in the following way:
@[@
y_j^\ell = \frac{1}{ \ell ! } Y_j^{(\ell)} (0)
@]@
forward
,
for @(@
i = 0 , \ldots , m-1
@)@, @(@
ty[ i *( k + 1) + k ]
@)@ is an output
and all the other components of
ty
are inputs.
reverse
,
all the components of
ty
are inputs.
afun
,
is the name of the AD function corresponding to this atomic
operation (as it is used in the source code).
CppAD uses the other functions,
where the arguments are vectors with elements of type
Base
,
to implement the function
afun(id, ax, ay)
where the argument are vectors with elements of type
AD<Base>
.
afun
argument
ax
has prototype
const Tvector< AD<Base> >& ax
It is the argument vector @(@
x \in \B{R}^n
@)@
at which the
AD<Base>
version of
@(@
y = f(x)
@)@ is to be evaluated.
The dimension of the domain space for @(@
y = f (x)
@)@
is specified by n
= ax.size()
,
which must be greater than zero.
afun
result
ay
has prototype
Tvector< AD<Base> >& ay
The input values of its elements
are not specified (must not matter).
Upon return, it is the
AD<Base>
version of the
result vector @(@
y = f(x)
@)@.
The dimension of the range space for @(@
y = f (x)
@)@
is specified by m
= ay.size()
,
which must be greater than zero.
afun(id, ax, ay)
must not be in parallel
mode.
In addition, the
atomic_one clear
routine cannot be called while in parallel mode.
forward
is a
user defined function
ok = forward(id, k, n, m, vx, vy, tx, ty)
that computes results during a forward
mode sweep.
For this call, we are given the Taylor coefficients in
tx
form order zero through
k
,
and the Taylor coefficients in
ty
with order less than
k
.
The
forward
routine computes the
k
order Taylor coefficients for @(@
y
@)@ using the definition
@(@
Y(t) = f[ X(t) ]
@)@.
For example, for @(@
i = 0 , \ldots , m-1
@)@,
@[@
\begin{array}{rcl}
y_i^0 & = & Y(0)
= f_i ( x^0 )
\\
y_i^1 & = & Y^{(1)} ( 0 )
= f_i^{(1)} ( x^0 ) X^{(1)} ( 0 )
= f_i^{(1)} ( x^0 ) x^1
\\
y_i^2
& = & \frac{1}{2 !} Y^{(2)} (0)
\\
& = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
+ \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
\\
& = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1
+ f_i^{(1)} ( x^0 ) x^2
\end{array}
@]@
Then, for @(@
i = 0 , \ldots , m-1
@)@, it sets
@[@
ty [ i * (k + 1) + k ] = y_i^k
@]@
The other components of
ty
must be left unchanged.
vx.size() > 0
and
k == 0
,
by calls to
afun
.
It is used,
with
vx.size() = 0
and
k
equal to the order of the derivative begin computed,
by calls to forward
.
forward
argument
vx
has prototype
const CppAD::vector<bool>& vx
The case
vx.size() > 0
occurs
once for each call to
afun
,
during the call,
and before any of the other callbacks corresponding to that call.
Hence such a call can be used to cache information attached to
the corresponding
id
(such as the elements of
vx
).
If
vx.size() > 0
then
k == 0
,
vx.size() >= n
, and
for @(@
j = 0 , \ldots , n-1
@)@,
vx[j]
is true if and only if
ax[j]
is a variable
.
If
vx.size() == 0
,
then
vy.size() == 0
and neither of these vectors
should be used.
forward
argument
vy
has prototype
CppAD::vector<bool>& vy
If
vy.size() == 0
, it should not be used.
Otherwise,
k == 0
and
vy.size() >= m
.
The input values of the elements of
vy
are not specified (must not matter).
Upon return, for @(@
j = 0 , \ldots , m-1
@)@,
vy[i]
is true if and only if
ay[j]
is a variable.
(CppAD uses
vy
to reduce the necessary computations.)
reverse
is a user defined function
ok = reverse(id, k, n, m, tx, ty, px, py)
that computes results during a reverse
mode sweep.
The input value of the vectors
tx
and
ty
contain Taylor coefficient, up to order
k
,
for @(@
X(t)
@)@ and @(@
Y(t)
@)@ respectively.
We use the @(@
\{ x_j^\ell \}
@)@ and @(@
\{ y_i^\ell \}
@)@
to denote these Taylor coefficients where the implicit range indices are
@(@
i = 0 , \ldots , m-1
@)@,
@(@
j = 0 , \ldots , n-1
@)@,
@(@
\ell = 0 , \ldots , k
@)@.
Using the calculations done by forward
,
the Taylor coefficients @(@
\{ y_i^\ell \}
@)@ are a function of the Taylor
coefficients for @(@
\{ x_j^\ell \}
@)@; i.e., given @(@
y = f(x)
@)@
we define the function
@(@
F : \B{R}^{n \times (k+1)} \rightarrow \B{R}^{m \times (k+1)}
@)@ by
@[@
y_i^\ell = F_i^\ell ( \{ x_j^\ell \} )
@]@
We use @(@
G : \B{R}^{m \times (k+1)} \rightarrow \B{R}
@)@
to denote an arbitrary scalar valued function of the Taylor coefficients for
@(@
Y(t)
@)@ and write @(@
z = G( \{ y_i^\ell \} )
@)@.
The reverse
routine
is given the derivative of @(@
z
@)@ with respect to
@(@
\{ y_i^\ell \}
@)@ and computes its derivative with respect
to @(@
\{ x_j^\ell \}
@)@.
k + 1
equal to the order of the derivative being calculated,
by calls to reverse
.
reverse
argument
py
has prototype
const CppAD::vector<Base>& py
and
py.size() >= (k + 1) * m
.
For @(@
i = 0 , \ldots , m-1
@)@ and @(@
\ell = 0 , \ldots , k
@)@,
@[@
py[ i * (k + 1 ) + \ell ] = \partial G / \partial y_i^\ell
@]@
If
py.size() > (k + 1) * m
,
the other components of
py
are not specified and should not be used.
reverse
argument
px
has prototype
CppAD::vector<Base>& px
and
px.size() >= (k + 1) * n
.
The input values of the elements of
px
are not specified (must not matter).
Upon return,
for @(@
j = 0 , \ldots , n-1
@)@ and @(@
p = 0 , \ldots , k
@)@,
@[@
\begin{array}{rcl}
px [ j * (k + 1) + p ] & = & \partial H / \partial x_j^p
\\
& = &
( \partial G / \partial \{ y_i^\ell \} )
( \partial \{ y_i^\ell \} / \partial x_j^p )
\\
& = &
\sum_{i=0}^{m-1} \sum_{\ell=0}^k
( \partial G / \partial y_i^\ell ) ( \partial y_i^\ell / \partial x_j^p )
\\
& = &
\sum_{i=0}^{m-1} \sum_{\ell=p}^k
py[ i * (k + 1 ) + \ell ] ( \partial F_i^\ell / \partial x_j^p )
\end{array}
@]@
Note that we have used the fact that for @(@
\ell < p
@)@,
@(@
\partial F_i^\ell / \partial x_j^p = 0
@)@.
If
px.size() > (k + 1) * n
,
the other components of
px
are not specified and should not be used.
for_jac_sparse
is a user defined function
ok = for_jac_sparse(id, n, m, q, r, s)
that is used to compute results during a forward Jacobian sparsity sweep.
For a fixed @(@
n \times q
@)@ matrix @(@
R
@)@,
the Jacobian of @(@
f( x + R * u)
@)@ with respect to @(@
u \in \B{R}^q
@)@ is
@[@
S(x) = f^{(1)} (x) * R
@]@
Given a sparsity pattern
for @(@
R
@)@,
for_jac_sparse
computes a sparsity pattern for @(@
S(x)
@)@.
for_jac_sparse
argument
q
has prototype
size_t q
It specifies the number of columns in
@(@
R \in \B{R}^{n \times q}
@)@ and the Jacobian
@(@
S(x) \in \B{R}^{m \times q}
@)@.
for_jac_sparse
argument
r
has prototype
const CppAD::vector< std::set<size_t> >& r
and
r.size() >= n
.
For @(@
j = 0 , \ldots , n-1
@)@,
all the elements of
r[j]
are between
zero and
q-1
inclusive.
This specifies a sparsity pattern for the matrix @(@
R
@)@.
for_jac_sparse
return value
s
has prototype
CppAD::vector< std::set<size_t> >& s
and
s.size() >= m
.
The input values of its sets
are not specified (must not matter). Upon return
for @(@
i = 0 , \ldots , m-1
@)@,
all the elements of
s[i]
are between
zero and
q-1
inclusive.
This represents a sparsity pattern for the matrix @(@
S(x)
@)@.
rev_jac_sparse
is a user defined function
ok = rev_jac_sparse(id, n, m, q, r, s)
that is used to compute results during a reverse Jacobian sparsity sweep.
For a fixed @(@
q \times m
@)@ matrix @(@
S
@)@,
the Jacobian of @(@
S * f( x )
@)@ with respect to @(@
x \in \B{R}^n
@)@ is
@[@
R(x) = S * f^{(1)} (x)
@]@
Given a sparsity pattern
for @(@
S
@)@,
rev_jac_sparse
computes a sparsity pattern for @(@
R(x)
@)@.
rev_jac_sparse
argument
q
has prototype
size_t q
It specifies the number of rows in
@(@
S \in \B{R}^{q \times m}
@)@ and the Jacobian
@(@
R(x) \in \B{R}^{q \times n}
@)@.
rev_jac_sparse
argument
s
has prototype
const CppAD::vector< std::set<size_t> >& s
and
s.size() >= m
.
For @(@
i = 0 , \ldots , m-1
@)@,
all the elements of
s[i]
are between zero and
q-1
inclusive.
This specifies a sparsity pattern for the matrix @(@
S^\R{T}
@)@.
rev_jac_sparse
return value
r
has prototype
CppAD::vector< std::set<size_t> >& r
and
r.size() >= n
.
The input values of its sets
are not specified (must not matter).
Upon return for @(@
j = 0 , \ldots , n-1
@)@,
all the elements of
r[j]
are between zero and
q-1
inclusive.
This represents a sparsity pattern for the matrix @(@
R(x)^\R{T}
@)@.
rev_hes_sparse
is a user defined function
ok = rev_hes_sparse(id, n, m, q, r, s, t, u, v)
There is an unspecified scalar valued function
@(@
g : \B{R}^m \rightarrow \B{R}
@)@.
Given a sparsity pattern for @(@
R
@)@
and information about the function @(@
z = g(y)
@)@,
this routine computes the sparsity pattern for
@[@
V(x) = (g \circ f)^{(2)}( x ) R
@]@
rev_hes_sparse
argument
q
has prototype
size_t q
It specifies the number of columns in the sparsity patterns.
rev_hes_sparse
argument
r
has prototype
const CppAD::vector< std::set<size_t> >& r
and
r.size() >= n
.
For @(@
j = 0 , \ldots , n-1
@)@,
all the elements of
r[j]
are between
zero and
q-1
inclusive.
This specifies a sparsity pattern for the matrix @(@
R \in \B{R}^{n \times q}
@)@.
rev_hes_sparse
argument
s
has prototype
const CppAD::vector<bool>& s
and
s.size() >= m
.
This specifies a sparsity pattern for the matrix
@(@
S(x) = g^{(1)} (y) \in \B{R}^{1 \times m}
@)@.
rev_hes_sparse
argument
t
has prototype
CppAD::vector<bool>& t
and
t.size() >= n
.
The input values of its elements
are not specified (must not matter).
Upon return it represents a sparsity pattern for the matrix
@(@
T(x) \in \B{R}^{1 \times n}
@)@ defined by
@[@
T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x)
@]@
rev_hes_sparse
argument
u
has prototype
const CppAD::vector< std::set<size_t> >& u
and
u.size() >= m
.
For @(@
i = 0 , \ldots , m-1
@)@,
all the elements of
u[i]
are between zero and
q-1
inclusive.
This specifies a sparsity pattern
for the matrix @(@
U(x) \in \B{R}^{m \times q}
@)@ defined by
@[@
\begin{array}{rcl}
U(x)
& = &
\partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0}
\\
& = &
\partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0}
\\
& = &
g^{(2)} (y) f^{(1)} (x) R
\end{array}
@]@
rev_hes_sparse
argument
v
has prototype
CppAD::vector< std::set<size_t> >& v
and
v.size() >= n
.
The input values of its elements
are not specified (must not matter).
Upon return, for @(@
j = 0, \ldots , n-1
@)@,
all the elements of
v[j]
are between zero and
q-1
inclusive.
This represents a sparsity pattern for the matrix
@(@
V(x) \in \B{R}^{n \times q}
@)@ defined by
@[@
\begin{array}{rcl}
V(x)
& = &
\partial_u [ \partial_x (g \circ f) ( x + R u ) ]_{u=0}
\\
& = &
\partial_u [ (g \circ f)^{(1)}( x + R u ) ]_{u=0}
\\
& = &
(g \circ f)^{(2)}( x ) R
\\
& = &
f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x) R
+
\sum_{i=1}^m [ g^{(1)} (y) ]_i \; f_i^{(2)} (x) R
\\
& = &
f^{(1)} (x)^\R{T} U(x)
+
\sum_{i=1}^m S(x)_i \; f_i^{(2)} (x) R
\end{array}
@]@
user_atomic<Base>::clear()
makes to work space available
to
for other uses by the same thread.
This should be called when you are done using the
atomic functions for a specific value of
Base
.
clear
routine cannot be called
while in parallel
execution mode.