CoDiPack  2.2.0
A Code Differentiation Package
SciComp TU Kaiserslautern
Loading...
Searching...
No Matches
preaccumulationHelper.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
36#pragma once
37
38#include <vector>
39
40#include "../../config.h"
41#include "../../expressions/lhsExpressionInterface.hpp"
42#include "../../misc/exceptions.hpp"
43#include "../../tapes/interfaces/fullTapeInterface.hpp"
44#include "../../tapes/tagging/tagTapeReverse.hpp"
45#include "../../traits/gradientTraits.hpp"
46#include "../../traits/tapeTraits.hpp"
47#include "../algorithms.hpp"
48#include "../data/jacobian.hpp"
49
51namespace codi {
52
75 template<typename T_Type, typename = void>
77 public:
78
80 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
81
82 using Real = typename Type::Real;
83 using Identifier = typename Type::Identifier;
84 using Gradient = typename Type::Gradient;
85
87 using Tape = CODI_DD(typename Type::Tape, CODI_DEFAULT_TAPE);
88 using Position = typename Tape::Position;
89
90 std::vector<Identifier> inputData;
92 std::vector<Identifier> outputData;
94 std::vector<Type*> outputValues;
96
97 protected:
98
100 std::vector<Gradient> storedAdjoints;
102
103 public:
104
108
110 template<typename... Inputs>
111 void addInput(Inputs const&... inputs) {
112 Tape& tape = Type::getTape();
113
114 if (tape.isActive()) {
115 addInputRecursive(inputs...);
116 }
117 }
118
120 template<typename... Inputs>
121 void start(Inputs const&... inputs) {
122 Tape& tape = Type::getTape();
123
125
126 if (tape.isActive()) {
127 inputData.clear();
128 outputData.clear();
129 outputValues.clear();
130
131 startPos = tape.getPosition();
132
133 addInputRecursive(inputs...);
134 }
135 }
136
138 template<typename... Outputs>
139 void addOutput(Outputs&... outputs) {
140 Tape& tape = Type::getTape();
141
142 if (tape.isActive()) {
143 addOutputRecursive(outputs...);
144 }
145 }
146
148 template<typename... Outputs>
149 void finish(bool const storeAdjoints, Outputs&... outputs) {
150 Tape& tape = Type::getTape();
151
152 if (tape.isActive()) {
153 addOutputRecursive(outputs...);
154
155 if (storeAdjoints) {
156 storeInputAdjoints();
157 }
158
159 tape.setPassive();
160 doPreaccumulation();
161 tape.setActive();
162
163 if (storeAdjoints) {
164 restoreInputAdjoints();
165 }
166 }
167
169 }
170
171 private:
172
173 void addInputLogic(Type const& input) {
174 EventSystem<Tape>::notifyPreaccAddInputListeners(Type::getTape(), input.getValue(), input.getIdentifier());
175 Identifier const& identifier = input.getIdentifier();
176 if (Type::getTape().getPassiveIndex() != identifier) {
177 inputData.push_back(identifier);
178 }
179 }
180
182 void addInputRecursive() {
183 // Terminator implementation.
184 }
185
186 template<typename... Inputs>
187 void addInputRecursive(Type const& input, Inputs const&... r) {
188 addInputLogic(input);
189 addInputRecursive(r...);
190 }
191
192 void addOutputLogic(Type& output) {
193 EventSystem<Tape>::notifyPreaccAddOutputListeners(Type::getTape(), output.value(), output.getIdentifier());
194 Identifier const& identifier = output.getIdentifier();
195 if (Type::getTape().getPassiveIndex() != identifier) {
196 outputData.push_back(identifier);
197 outputValues.push_back(&output);
198 }
199 }
200
202 void addOutputRecursive() {
203 // Terminator implementation.
204 }
205
206 template<typename... Outputs>
207 void addOutputRecursive(Type& output, Outputs&... r) {
208 addOutputLogic(output);
209 addOutputRecursive(r...);
210 }
211
212 void storeInputAdjoints() {
213 Tape& tape = Type::getTape();
214
215 if (storedAdjoints.size() < inputData.size()) {
216 storedAdjoints.resize(inputData.size());
217 }
218
219 for (size_t i = 0; i < inputData.size(); ++i) {
220 Identifier index = inputData[i];
221 Gradient& adjoint = tape.gradient(index);
222 storedAdjoints[i] = adjoint;
223 adjoint = Gradient();
224 }
225 }
226
227 void restoreInputAdjoints() {
228 Tape& tape = Type::getTape();
229
230 for (size_t i = 0; i < inputData.size(); ++i) {
231 Identifier index = inputData[i];
232 tape.gradient(index) = storedAdjoints[i];
233 }
234 }
235
236 void doPreaccumulation() {
237 // Perform the accumulation of the tape part.
238 Tape& tape = Type::getTape();
239
240 Position endPos = tape.getPosition();
241 if (jacobian.getM() != outputData.size() || jacobian.getN() != inputData.size()) {
242 jacobian.resize(outputData.size(), inputData.size());
243 }
244
245 // Manage adjoints manually to reduce the impact of locking on the performance.
246 tape.resizeAdjointVector();
247 tape.beginUseAdjointVector();
248
250 outputData.data(), outputData.size(), jacobian,
252
253 // Store the Jacobian matrix.
254 tape.resetTo(startPos, true, AdjointsManagement::Manual);
255
256 tape.endUseAdjointVector();
257
258 for (size_t curOut = 0; curOut < outputData.size(); ++curOut) {
259 Type& value = *outputValues[curOut];
260 if (0 != jacobian.nonZerosRow(curOut)) {
261 int nonZerosLeft = jacobian.nonZerosRow(curOut);
262 jacobian.nonZerosRow(curOut) = 0;
263
264 // We need to initialize with the output's current identifier such that it is correctly deleted in
265 // storeManual.
266 Identifier lastIdentifier = value.getIdentifier();
267 bool staggeringActive = false;
268 int curIn = 0;
269
270 // Push statements as long as there are nonzeros left.
271 // If there are more than MaxStatementIntValue nonzeros, then we need to stagger the
272 // statement pushes:
273 // e.g. The reverse mode of w = f(u0, ..., u530) which is \bar u_i += df/du_i * \bar w for i = 0 ... 530 is
274 // separated into
275 // Statement 1:
276 // \bar u_i += df/du_i * \bar t_1 for i = 0 ... 253 (254 entries)
277 // Statement 2:
278 // \bar t_1 += \bar w (1 entry)
279 // \bar u_i += df/du_i * \bar t_2 for i = 254 ... 506 (253 entries)
280 // Statement 3:
281 // \bar t_2 += \bar w (1 entry)
282 // \bar u_i += df/du_i * \bar w for i = 507 ... 530 (24 entries)
283 //
284 while (nonZerosLeft > 0) {
285 // Calculate the number of Jacobians for this statement.
286 int jacobiansForStatement = nonZerosLeft;
287 if (jacobiansForStatement > (int)Config::MaxArgumentSize) {
288 jacobiansForStatement = (int)Config::MaxArgumentSize - 1;
289 if (staggeringActive) { // Except in the first round, one Jacobian is reserved for the staggering.
290 jacobiansForStatement -= 1;
291 }
292 }
293 nonZerosLeft -= jacobiansForStatement; // Update nonzeros so that we know if it is the last round.
294
295 Identifier storedIdentifier = lastIdentifier;
296 // storeManual creates a new identifier which is either the identifier of the output w or the temporary
297 // staggering variables t_1, t_2, ...
298 tape.storeManual(value.getValue(), lastIdentifier, jacobiansForStatement + (int)staggeringActive);
299 if (staggeringActive) { // Not the first staggering so push the last output.
300 tape.pushJacobianManual(1.0, 0.0, storedIdentifier);
301 }
302
303 // Push the rest of the Jacobians for the statement.
304 while (jacobiansForStatement > 0) {
305 if (Real() != (Real)jacobian(curOut, curIn)) {
306 tape.pushJacobianManual(jacobian(curOut, curIn), 0.0, inputData[curIn]);
307 jacobiansForStatement -= 1;
308 }
309 curIn += 1;
310 }
311
312 staggeringActive = true;
313 }
314
315 value.getIdentifier() = lastIdentifier; /* now set gradient data for the real output value */
316 } else {
317 // Disable tape index since there is no dependency.
318 tape.destroyIdentifier(value.value(), value.getIdentifier());
319 }
320 }
321 }
322 };
323
324#ifndef DOXYGEN_DISABLE
330 struct PreaccumulationHelperNoOpBase {
331 public:
332
334 template<typename... Inputs>
335 void addInput(Inputs const&... inputs) {
336 CODI_UNUSED(inputs...);
337 // Do nothing.
338 }
339
341 template<typename... Inputs>
342 void start(Inputs const&... inputs) {
343 CODI_UNUSED(inputs...);
344 // Do nothing.
345 }
346
348 template<typename... Outputs>
349 void addOutput(Outputs&... outputs) {
350 CODI_UNUSED(outputs...);
351 // Do nothing.
352 }
353
355 template<typename... Outputs>
356 void finish(bool const storeAdjoints, Outputs&... outputs) {
357 CODI_UNUSED(storeAdjoints, outputs...);
358 // Do nothing.
359 }
360 };
361
363 template<typename Type>
364 struct PreaccumulationHelper<Type, TapeTraits::EnableIfForwardTape<typename Type::Tape>>
365 : public PreaccumulationHelperNoOpBase {};
366
368 template<>
369 struct PreaccumulationHelper<double, void> : public PreaccumulationHelperNoOpBase {};
370
378 template<typename T_Type>
379 struct PreaccumulationHelper<
380 T_Type, typename enable_if_same<typename T_Type::Tape,
381 TagTapeReverse<typename T_Type::Real, typename T_Type::Tape::Tag>>::type> {
382 public:
383
384 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
385
386 using Tape = CODI_DD(typename Type::Tape, CODI_T(TagTapeReverse<Type, int>));
387 using Tag = typename Tape::Tag;
388
389 private:
390
391 std::vector<Type const*> inputLocations;
392 std::vector<Type*> outputLocations;
393 Tag oldTag;
394
395 public:
396
399
401 template<typename... Inputs>
402 void addInput(Inputs const&... inputs) {
403 Tape& tape = getTape();
404
405 if (tape.isActive() && tape.isPreaccumulationHandlingEnabled()) {
406 addInputRecursive(inputs...);
407 }
408 }
409
411 template<typename... Inputs>
412 void start(Inputs const&... inputs) {
413 Tape& tape = getTape();
414
415 if (tape.isActive() && tape.isPreaccumulationHandlingEnabled()) {
416 inputLocations.clear();
417 outputLocations.clear();
418 oldTag = tape.getCurTag();
419 tape.setCurTag(tape.getPreaccumulationHandlingTag());
420
421 addInputRecursive(inputs...);
422 }
423 }
424
426 template<typename... Outputs>
427 void addOutput(Outputs&... outputs) {
428 Tape& tape = getTape();
429
430 if (tape.isActive() && tape.isPreaccumulationHandlingEnabled()) {
431 addOutputRecursive(outputs...);
432 }
433 }
434
436 template<typename... Outputs>
437 void finish(bool const storeAdjoints, Outputs&... outputs) {
438 CODI_UNUSED(storeAdjoints);
439
440 Tape& tape = getTape();
441
442 if (tape.isActive() && tape.isPreaccumulationHandlingEnabled()) {
443 addOutputRecursive(outputs...);
444
445 tape.setCurTag(oldTag);
446 for (Type const* curInput : inputLocations) {
447 tape.setTagOnVariable(*curInput);
448 }
449 for (Type* curOutput : outputLocations) {
450 tape.setTagOnVariable(*curOutput);
451 }
452 }
453 }
454
455 private:
456
458 void addInputRecursive() {
459 // Terminator implementation.
460 }
461
462 template<typename... Inputs>
463 void addInputRecursive(Type const& input, Inputs const&... r) {
464 handleInput(input);
465 addInputRecursive(r...);
466 }
467
468 void handleInput(Type const& input) {
469 inputLocations.push_back(&input);
470 getTape().setTagOnVariable(input);
471 }
472
474 void addOutputRecursive() {
475 // Terminator implementation.
476 }
477
478 template<typename... Outputs>
479 void addOutputRecursive(Type& output, Outputs&... r) {
480 handleOutput(output);
481 addOutputRecursive(r...);
482 }
483
484 void handleOutput(Type& value) {
485 outputLocations.push_back(&value);
486 }
487
488 Tape& getTape() {
489 return Type::getTape();
490 }
491 };
492#endif
493}
#define CODI_INLINE
See codi::Config::ForcedInlines.
Definition config.h:457
#define CODI_DD(Type, Default)
Abbreviation for CODI_DECLARE_DEFAULT.
Definition macros.hpp:94
#define CODI_T(...)
Abbreviation for CODI_TEMPLATE.
Definition macros.hpp:111
size_t constexpr MaxArgumentSize
Maximum number of arguments in a statement.
Definition config.h:120
typename std::enable_if< IsForwardTape< Tape >::value >::type EnableIfForwardTape
Enable if wrapper for IsForwardTape.
Definition tapeTraits.hpp:88
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
std::enable_if< std::is_same< T1, T2 >::value, R > enable_if_same
Enable if abbreviation for "std::is_same".
Definition enableIfHelpers.hpp:51
@ Manual
Do not perform any bounds checking, locking, or resizing.
static void computeJacobian(Tape &tape, Position const &start, Position const &end, Identifier const *input, size_t const inputSize, Identifier const *output, size_t const outputSize, Jac &jac, AdjointsManagement adjointsManagement=AdjointsManagement::Automatic)
Compute the Jacobian with multiple tape sweeps.
Definition algorithms.hpp:141
static void notifyPreaccFinishListeners(Tape &tape)
Invoke callbacks for PreaccFinish events.
Definition eventSystem.hpp:579
static void notifyPreaccAddOutputListeners(Tape &tape, Real &value, Identifier &identifier)
Invoke callbacks for PreaccAddOutput events.
Definition eventSystem.hpp:636
static void notifyPreaccAddInputListeners(Tape &tape, Real const &value, Identifier const &identifier)
Invoke callbacks for PreaccAddInput events.
Definition eventSystem.hpp:607
static void notifyPreaccStartListeners(Tape &tape)
Invoke callbacks for PreaccStart events.
Definition eventSystem.hpp:553
Adds counting of nonzero entries.
Definition jacobian.hpp:129
int & nonZerosRow(size_t const i)
Reference to the number of nonzero entries for the specified row.
Definition jacobian.hpp:167
void resize(size_t const m, size_t const n)
Resize the Jacobian.
Definition jacobian.hpp:161
size_t getM() const
Get size of rows (output variables).
Definition jacobian.hpp:78
size_t getN() const
Get size of columns (input variables).
Definition jacobian.hpp:83
Stores the Jacobian matrix for a code section.
Definition preaccumulationHelper.hpp:76
void finish(bool const storeAdjoints, Outputs &... outputs)
Finish the preaccumulation region and perform the preaccumulation. See addOutput() for outputs.
Definition preaccumulationHelper.hpp:149
std::vector< Gradient > storedAdjoints
If adjoints of inputs should be stored, before the preaccumulation.
Definition preaccumulationHelper.hpp:100
Position startPos
Starting position for the region.
Definition preaccumulationHelper.hpp:99
PreaccumulationHelper()
Constructor.
Definition preaccumulationHelper.hpp:106
typename Type::Tape Tape
See LhsExpressionInterface.
Definition preaccumulationHelper.hpp:87
void addInput(Inputs const &... inputs)
Add multiple additional inputs. Inputs need to be of type Type. Called after start().
Definition preaccumulationHelper.hpp:111
std::vector< Identifier > inputData
Definition preaccumulationHelper.hpp:90
T_Type Type
See PreaccumulationHelper.
Definition preaccumulationHelper.hpp:80
typename Tape::Position Position
See PositionalEvaluationTapeInterface.
Definition preaccumulationHelper.hpp:88
typename Type::Identifier Identifier
See LhsExpressionInterface.
Definition preaccumulationHelper.hpp:83
typename Type::Gradient Gradient
See LhsExpressionInterface.
Definition preaccumulationHelper.hpp:84
std::vector< Type * > outputValues
Definition preaccumulationHelper.hpp:94
std::vector< Identifier > outputData
Definition preaccumulationHelper.hpp:92
typename Type::Real Real
See LhsExpressionInterface.
Definition preaccumulationHelper.hpp:82
void start(Inputs const &... inputs)
Starts a preaccumulation region. Resets the internal state. See addInputs() for inputs.
Definition preaccumulationHelper.hpp:121
void addOutput(Outputs &... outputs)
Add multiple additional outputs. Outputs need to be of type Type. Called before finish().
Definition preaccumulationHelper.hpp:139
JacobianCountNonZerosRow< Real > jacobian
Jacobian for the preaccumulation.
Definition preaccumulationHelper.hpp:101