CoDiPack  2.2.0
A Code Differentiation Package
SciComp TU Kaiserslautern
Loading...
Searching...
No Matches
jacobianReuseTape.hpp
1/*
2 * CoDiPack, a Code Differentiation Package
3 *
4 * Copyright (C) 2015-2024 Chair for Scientific Computing (SciComp), University of Kaiserslautern-Landau
5 * Homepage: http://scicomp.rptu.de
6 * Contact: Prof. Nicolas R. Gauger (codi@scicomp.uni-kl.de)
7 *
8 * Lead developers: Max Sagebaum, Johannes Blühdorn (SciComp, University of Kaiserslautern-Landau)
9 *
10 * This file is part of CoDiPack (http://scicomp.rptu.de/software/codi).
11 *
12 * CoDiPack is free software: you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * CoDiPack is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty
19 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
20 *
21 * See the GNU General Public License for more details.
22 * You should have received a copy of the GNU
23 * General Public License along with CoDiPack.
24 * If not, see <http://www.gnu.org/licenses/>.
25 *
26 * For other licensing options please contact us.
27 *
28 * Authors:
29 * - SciComp, University of Kaiserslautern-Landau:
30 * - Max Sagebaum
31 * - Johannes Blühdorn
32 * - Former members:
33 * - Tim Albring
34 */
35#pragma once
36
37#include <algorithm>
38#include <type_traits>
39
40#include "../config.h"
41#include "../expressions/lhsExpressionInterface.hpp"
42#include "../expressions/logic/compileTimeTraversalLogic.hpp"
43#include "../expressions/logic/traversalLogic.hpp"
44#include "../misc/macros.hpp"
45#include "../traits/expressionTraits.hpp"
46#include "data/chunk.hpp"
47#include "indices/indexManagerInterface.hpp"
48#include "interfaces/editingTapeInterface.hpp"
49#include "interfaces/reverseTapeInterface.hpp"
50#include "jacobianBaseTape.hpp"
51
53namespace codi {
54
62 template<typename T_TapeTypes>
64 : public JacobianBaseTape<T_TapeTypes, JacobianReuseTape<T_TapeTypes>>,
65 public EditingTapeInterface<typename JacobianBaseTape<T_TapeTypes, JacobianReuseTape<T_TapeTypes>>::Position> {
66 public:
67
68 using TapeTypes = CODI_DD(T_TapeTypes,
71
73 friend Base;
74
75 using Real = typename TapeTypes::Real;
76 using Gradient = typename TapeTypes::Gradient;
77 using IndexManager = typename TapeTypes::IndexManager;
78 using Identifier = typename TapeTypes::Identifier;
79 using Position = typename Base::Position;
80 using StatementData = typename TapeTypes::StatementData;
81
82 CODI_STATIC_ASSERT(!IndexManager::IsLinear, "This class requires an index manager with a reuse scheme.");
83
86
87 /*******************************************************************************/
90
92
94 void clearAdjoints(Position const& start, Position const& end,
96 if (AdjointsManagement::Automatic == adjointsManagement) {
97 this->adjoints.beginUse();
98 }
99
100 Identifier adjointsSize = (Identifier)this->adjoints.size();
101
102 auto clearFunc = [&adjointsSize, this](Identifier* index, Config::ArgumentSize* stmtSize) {
103 CODI_UNUSED(stmtSize);
104
105 if (*index < adjointsSize) {
106 this->adjoints[*index] = Gradient();
107 }
108 };
109
110 using StmtPosition = typename StatementData::Position;
111 StmtPosition startStmt = this->llfByteData.template extractPosition<StmtPosition>(start);
112 StmtPosition endStmt = this->llfByteData.template extractPosition<StmtPosition>(end);
113
114 this->statementData.forEachReverse(startStmt, endStmt, clearFunc);
115
116 if (AdjointsManagement::Automatic == adjointsManagement) {
117 this->adjoints.endUse();
118 }
119 }
120
122
123 protected:
124
125 /*******************************************************************************/
128
131 CODI_INLINE void pushStmtData(Identifier const& index, Config::ArgumentSize const& numberOfArguments) {
132 this->statementData.pushData(index, numberOfArguments);
133 }
134
136 template<typename Adjoint>
138 /* data from call */
139 JacobianReuseTape& tape, Adjoint* adjointVector,
140 /* data from low level function byte data vector */
141 size_t& curLLFByteDataPos, size_t const& endLLFByteDataPos, char* dataPtr,
142 /* data from low level function info data vector */
143 size_t& curLLFInfoDataPos, size_t const& endLLFInfoDataPos, Config::LowLevelFunctionToken* const tokenPtr,
144 Config::LowLevelFunctionDataSize* const dataSizePtr,
145 /* data from jacobian vector */
146 size_t& curJacobianPos, size_t const& endJacobianPos, Real const* const rhsJacobians,
147 Identifier const* const rhsIdentifiers,
148 /* data from statement vector */
149 size_t& curStmtPos, size_t const& endStmtPos, Identifier const* const lhsIdentifiers,
150 Config::ArgumentSize const* const numberOfJacobians) {
151 CODI_UNUSED(endJacobianPos, endLLFByteDataPos, endLLFInfoDataPos);
152
153 typename Base::template VectorAccess<Adjoint> vectorAccess(adjointVector);
154
155 while (curStmtPos < endStmtPos) CODI_Likely {
156 Config::ArgumentSize const argsSize = numberOfJacobians[curStmtPos];
157
159 Base::template callLowLevelFunction<LowLevelFunctionEntryCallKind::Forward>(
160 tape, true, curLLFByteDataPos, dataPtr, curLLFInfoDataPos, tokenPtr, dataSizePtr, &vectorAccess);
161 } else CODI_Likely {
162 Adjoint lhsAdjoint = Adjoint();
163 Base::incrementTangents(adjointVector, lhsAdjoint, argsSize, curJacobianPos, rhsJacobians, rhsIdentifiers);
164
165 adjointVector[lhsIdentifiers[curStmtPos]] = lhsAdjoint;
166
168 tape, lhsIdentifiers[curStmtPos], GradientTraits::dim<Adjoint>(),
169 GradientTraits::toArray(lhsAdjoint).data());
170 }
171
172 curStmtPos += 1;
173 }
174 }
175
177 template<typename Adjoint>
179 /* data from call */
180 JacobianReuseTape& tape, Adjoint* adjointVector,
181 /* data from low level function byte data vector */
182 size_t& curLLFByteDataPos, size_t const& endLLFByteDataPos, char* dataPtr,
183 /* data from low level function info data vector */
184 size_t& curLLFInfoDataPos, size_t const& endLLFInfoDataPos, Config::LowLevelFunctionToken* const tokenPtr,
185 Config::LowLevelFunctionDataSize* const dataSizePtr,
186 /* data from jacobianData */
187 size_t& curJacobianPos, size_t const& endJacobianPos, Real const* const rhsJacobians,
188 Identifier const* const rhsIdentifiers,
189 /* data from statementData */
190 size_t& curStmtPos, size_t const& endStmtPos, Identifier const* const lhsIdentifiers,
191 Config::ArgumentSize const* const numberOfJacobians) {
192 CODI_UNUSED(endJacobianPos, endLLFByteDataPos, endLLFInfoDataPos);
193
194 typename Base::template VectorAccess<Adjoint> vectorAccess(adjointVector);
195
196 while (curStmtPos > endStmtPos) CODI_Likely {
197 curStmtPos -= 1;
198
199 Config::ArgumentSize const argsSize = numberOfJacobians[curStmtPos];
200
202 Base::template callLowLevelFunction<LowLevelFunctionEntryCallKind::Reverse>(
203 tape, false, curLLFByteDataPos, dataPtr, curLLFInfoDataPos, tokenPtr, dataSizePtr, &vectorAccess);
204 } else CODI_Likely {
205 Adjoint const lhsAdjoint = adjointVector[lhsIdentifiers[curStmtPos]];
206
208 tape, lhsIdentifiers[curStmtPos], GradientTraits::dim<Adjoint>(),
209 GradientTraits::toArray(lhsAdjoint).data());
210
211 adjointVector[lhsIdentifiers[curStmtPos]] = Adjoint();
212 Base::incrementAdjoints(adjointVector, lhsAdjoint, argsSize, curJacobianPos, rhsJacobians, rhsIdentifiers);
213 }
214 }
215 }
216
218
219 public:
220
221 /*******************************************************************************/
224
228 CODI_INLINE void erase(Position const& start, Position const& end) {
229 JacobianReuseTape emptyTape;
230 erase(start, end, emptyTape);
231 }
232
233 // clang-format off
235 // clang-format on
236 CODI_INLINE void erase(Position const& start, Position const& end, JacobianReuseTape& emptyTape) {
237 // Store the tail after the part to be erased in the helper tape.
238 emptyTape.append(*this, end, this->getPosition());
239 // Reset the tape to before the erased part and re-append the tail. This accounts for external function position
240 // correction.
241 this->resetTo(start);
242 this->append(emptyTape, emptyTape.getZeroPosition(), emptyTape.getPosition());
243 emptyTape.reset();
244 }
245
247 CODI_INLINE void append(JacobianReuseTape& srcTape, Position const& start, Position const& end) {
248 srcTape.llfByteData.evaluateForward(start, end, JacobianReuseTape::internalAppend, this);
249 }
250
252
253 private:
254
255 static CODI_INLINE void internalAppend(
256 /* data from call */
257 JacobianReuseTape* dstTape,
258 /* data from low level function byte data vector */
259 size_t& curLLFByteDataPos, size_t const& endLLFByteDataPos, char* dataPtr,
260 /* data from low level function info data vector */
261 size_t& curLLFInfoDataPos, size_t const& endLLFInfoDataPos, Config::LowLevelFunctionToken* const tokenPtr,
262 Config::LowLevelFunctionDataSize* const dataSizePtr,
263 /* data from jacobianData */
264 size_t& curJacobianPos, size_t const& endJacobianPos, Real const* const rhsJacobians,
265 Identifier const* const rhsIdentifiers,
266 /* data from statementData */
267 size_t& curStmtPos, size_t const& endStmtPos, Identifier const* const lhsIdentifiers,
268 Config::ArgumentSize const* const numberOfJacobians) {
269 CODI_UNUSED(endLLFByteDataPos, endLLFInfoDataPos, endJacobianPos);
270
271 while (curStmtPos < endStmtPos) {
272 Config::ArgumentSize const argsSize = numberOfJacobians[curStmtPos];
274 Config::LowLevelFunctionToken token = tokenPtr[curLLFInfoDataPos];
275 size_t dataSize = dataSizePtr[curLLFInfoDataPos];
276
277 // Create the store on the new tape.
278 ByteDataView dstDataStore = {};
279 dstTape->pushLowLevelFunction(token, dataSize, dstDataStore);
280
281 // Copy the data.
282 dstDataStore.write(&dataPtr[curLLFByteDataPos], dataSize);
283
284 curLLFInfoDataPos += 1;
285 curLLFByteDataPos += dataSize;
286 } else CODI_Likely {
287 // Manual statement push.
288 dstTape->statementData.reserveItems(1);
289 dstTape->jacobianData.reserveItems(numberOfJacobians[curStmtPos]);
290
291 dstTape->pushStmtData(lhsIdentifiers[curStmtPos], numberOfJacobians[curStmtPos]);
292 size_t curJacobianEnd = curJacobianPos + numberOfJacobians[curStmtPos];
293
294 while (curJacobianPos < curJacobianEnd) {
295 dstTape->jacobianData.pushData(rhsJacobians[curJacobianPos], rhsIdentifiers[curJacobianPos]);
296 ++curJacobianPos;
297 }
298 }
299
300 ++curStmtPos;
301 }
302 }
303 };
304}
#define CODI_Unlikely
Declare unlikely evaluation of an execution path.
Definition config.h:399
#define CODI_INLINE
See codi::Config::ForcedInlines.
Definition config.h:457
#define CODI_Likely
Declare likely evaluation of an execution path.
Definition config.h:397
#define CODI_DD(Type, Default)
Abbreviation for CODI_DECLARE_DEFAULT.
Definition macros.hpp:94
#define CODI_STATIC_ASSERT(cond, message)
Static assert in CoDiPack.
Definition macros.hpp:123
#define CODI_T(...)
Abbreviation for CODI_TEMPLATE.
Definition macros.hpp:111
uint16_t LowLevelFunctionDataSize
Size store type for a low level function.
Definition config.h:98
uint16_t LowLevelFunctionToken
Token type for low level functions in the tapes.
Definition config.h:108
size_t constexpr StatementLowLevelFunctionTag
Statement tag for low level functions.
Definition config.h:126
uint8_t ArgumentSize
Type for the number of arguments in statements.
Definition config.h:117
std::array< AtomicTraits::RemoveAtomic< typename TraitsImplementation< Gradient >::Real >, TraitsImplementation< Gradient >::dim > toArray(Gradient const &gradient)
Converts the (possibly multi-component) gradient to an array of Reals.
Definition gradientTraits.hpp:116
CoDiPack - Code Differentiation Package.
Definition codi.hpp:90
void CODI_UNUSED(Args const &...)
Disable unused warnings for an arbitrary number of arguments.
Definition macros.hpp:46
AdjointsManagement
Policies for management of the tape's interal adjoints.
Definition tapeParameters.hpp:98
@ Automatic
Manage internal adjoints automatically, including locking, bounds checking, and resizing.
Represents a concrete lvalue in the CoDiPack expression tree.
Definition activeType.hpp:52
Implementation of VectorAccessInterface for adjoint vectors.
Definition adjointVectorAccess.hpp:59
Definition byteDataView.hpp:51
T * write(T const &data)
Write a single entry of type T.
Definition byteDataView.hpp:114
Data is stored chunk-wise in this DataInterface implementation. If a chunk runs out of space,...
Definition chunkedData.hpp:64
Position getPosition() const
Current position of the tape.
Definition commonTapeImplementation.hpp:629
Position getZeroPosition() const
Initial position of the tape.
Definition commonTapeImplementation.hpp:634
void reset(bool resetAdjoints=true, AdjointsManagement adjointsManagement=AdjointsManagement::Automatic)
Reset the tape to the initial state for a fresh recording.
Definition commonTapeImplementation.hpp:354
LowLevelFunctionByteData llfByteData
Byte data for low level functions.
Definition commonTapeImplementation.hpp:155
void resetTo(Position const &pos, bool resetAdjoints=true, AdjointsManagement adjointsManagement=AdjointsManagement::Automatic)
Reset the tape to the provided position.
Definition commonTapeImplementation.hpp:668
Edit tapes after they have been recorded.
Definition editingTapeInterface.hpp:65
static void notifyStatementEvaluateListeners(Tape &tape, Identifier const &lhsIdentifier, size_t sizeLhsAdjoint, Real const *lhsAdjoint)
Invoke callbacks for StatementEvaluate events.
Definition eventSystem.hpp:712
Indices enable the mapping of primal values to their adjoint counterparts.
Definition indexManagerInterface.hpp:78
Base class for all standard Jacobian tape implementations.
Definition jacobianBaseTape.hpp:125
JacobianData jacobianData
Data stream for argument specific data.
Definition jacobianBaseTape.hpp:171
typename Base::Position Position
See TapeTypesInterface.
Definition jacobianBaseTape.hpp:149
StatementData statementData
Data stream for statement specific data.
Definition jacobianBaseTape.hpp:170
Adjoints adjoints
Evaluation vector for AD.
Definition jacobianBaseTape.hpp:173
void clearAdjoints(AdjointsManagement adjointsManagement=AdjointsManagement::Automatic)
Clear all adjoint values, that is, set them to zero.
Definition jacobianBaseTape.hpp:468
static void incrementAdjoints(Adjoint *adjointVector, Adjoint const &lhsAdjoint, Config::ArgumentSize const &numberOfArguments, size_t &curJacobianPos, Real const *const rhsJacobians, Identifier const *const rhsIdentifiers)
Performs the AD reverse equation for a statement.
Definition jacobianBaseTape.hpp:520
void pushLowLevelFunction(Config::LowLevelFunctionToken token, size_t size, ByteDataView &data)
Push a low level function to the tape.
Definition jacobianBaseTape.hpp:774
static void incrementTangents(Adjoint const *const adjointVector, Adjoint &lhsAdjoint, Config::ArgumentSize const &numberOfArguments, size_t &curJacobianPos, Real const *const rhsJacobians, Identifier const *const rhsIdentifiers)
Performs the AD forward equation for a statement.
Definition jacobianBaseTape.hpp:542
Final implementation for a Jacobian tape with a reuse index management.
Definition jacobianReuseTape.hpp:65
void erase(Position const &start, Position const &end)
Erase a part of the tape. It has to hold start <= end.
Definition jacobianReuseTape.hpp:228
void append(JacobianReuseTape &srcTape, Position const &start, Position const &end)
Definition jacobianReuseTape.hpp:247
typename TapeTypes::Gradient Gradient
See TapeTypesInterface.
Definition jacobianReuseTape.hpp:76
typename TapeTypes::StatementData StatementData
See JacobianTapeTypes.
Definition jacobianReuseTape.hpp:80
void erase(Position const &start, Position const &end, JacobianReuseTape &emptyTape)
Erase a part of the tape. It has to hold start <= end. This variant of erase takes a reference to an ...
Definition jacobianReuseTape.hpp:236
T_TapeTypes TapeTypes
See JacobianReuseTape.
Definition jacobianReuseTape.hpp:68
void clearAdjoints(Position const &start, Position const &end, AdjointsManagement adjointsManagement=AdjointsManagement::Automatic)
Clear all adjoints that would be set in a tape evaluation from start to end. It has to hold start >= ...
Definition jacobianReuseTape.hpp:94
JacobianReuseTape()
Constructor.
Definition jacobianReuseTape.hpp:85
typename TapeTypes::Identifier Identifier
See TapeTypesInterface.
Definition jacobianReuseTape.hpp:78
static void internalEvaluateForward_EvalStatements(JacobianReuseTape &tape, Adjoint *adjointVector, size_t &curLLFByteDataPos, size_t const &endLLFByteDataPos, char *dataPtr, size_t &curLLFInfoDataPos, size_t const &endLLFInfoDataPos, Config::LowLevelFunctionToken *const tokenPtr, Config::LowLevelFunctionDataSize *const dataSizePtr, size_t &curJacobianPos, size_t const &endJacobianPos, Real const *const rhsJacobians, Identifier const *const rhsIdentifiers, size_t &curStmtPos, size_t const &endStmtPos, Identifier const *const lhsIdentifiers, Config::ArgumentSize const *const numberOfJacobians)
Perform a forward evaluation of the tape. Arguments are from the recursive eval methods of the DataIn...
Definition jacobianReuseTape.hpp:137
friend Base
Allow the base class to call protected and private methods.
Definition jacobianReuseTape.hpp:73
static void internalEvaluateReverse_EvalStatements(JacobianReuseTape &tape, Adjoint *adjointVector, size_t &curLLFByteDataPos, size_t const &endLLFByteDataPos, char *dataPtr, size_t &curLLFInfoDataPos, size_t const &endLLFInfoDataPos, Config::LowLevelFunctionToken *const tokenPtr, Config::LowLevelFunctionDataSize *const dataSizePtr, size_t &curJacobianPos, size_t const &endJacobianPos, Real const *const rhsJacobians, Identifier const *const rhsIdentifiers, size_t &curStmtPos, size_t const &endStmtPos, Identifier const *const lhsIdentifiers, Config::ArgumentSize const *const numberOfJacobians)
Perform a reverse evaluation of the tape. Arguments are from the recursive eval methods of the DataIn...
Definition jacobianReuseTape.hpp:178
typename TapeTypes::Real Real
See TapeTypesInterface.
Definition jacobianReuseTape.hpp:75
typename TapeTypes::IndexManager IndexManager
See TapeTypesInterface.
Definition jacobianReuseTape.hpp:77
void pushStmtData(Identifier const &index, Config::ArgumentSize const &numberOfArguments)
Add statement specific data to the data streams.
Definition jacobianReuseTape.hpp:131
typename Base::Position Position
See TapeTypesInterface.
Definition jacobianReuseTape.hpp:79
Type definitions for the Jacobian tapes.
Definition jacobianBaseTape.hpp:75