CoDiPack  2.2.0
A Code Differentiation Package
SciComp TU Kaiserslautern
Loading...
Searching...
No Matches
tapeHelper.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 <vector>
38
39#include "../../config.h"
40#include "../../expressions/lhsExpressionInterface.hpp"
41#include "../../traits/realTraits.hpp"
42#include "../../traits/tapeTraits.hpp"
43#include "../algorithms.hpp"
44#include "../data/hessian.hpp"
45#include "../data/jacobian.hpp"
46
48namespace codi {
49
102 template<typename T_Type, typename T_Impl>
104 public:
105
107 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
108 using Impl = CODI_DD(T_Impl, TapeHelperBase);
109
110 using Real = typename Type::Real;
111 using Identifier = typename Type::Identifier;
112 using Gradient = typename Type::Gradient;
113
115
118
119 protected:
120
121 using Tape = typename Type::Tape;
122
124
125 std::vector<Identifier> inputValues;
126 std::vector<Identifier> outputValues;
127
129
130 public:
131
134
136 virtual ~TapeHelperBase() {}
137
149
161
171 JacobianType* jacPointer = new JacobianType(getOutputSize(), getInputSize());
172
173 return *jacPointer;
174 }
175
185 HessianType* hesPointer = new HessianType(getOutputSize(), getInputSize());
186
187 return *hesPointer;
188 }
189
199 return new Real[getInputSize()];
200 }
201
211 return new Real[getOutputSize()];
212 }
213
216 delete[] vec;
217
218 vec = nullptr;
219 }
220
223 JacobianType* jacPointer = &jac;
224
225 delete jacPointer;
226 }
227
230 HessianType* hesPointer = &hes;
231
232 delete hesPointer;
233 }
234
237 delete[] vec;
238
239 vec = nullptr;
240 }
241
244 size_t getInputSize() {
245 return inputValues.size();
246 }
247
250 size_t getOutputSize() {
251 return outputValues.size();
252 }
253
265 void registerInput(Type& value) {
266 tape.registerInput(value);
267 inputValues.push_back(value.getIdentifier());
268 }
269
281 void registerOutput(Type& value) {
282 tape.registerOutput(value);
283 outputValues.push_back(value.getIdentifier());
284 }
285
288 tape.reset();
289 inputValues.clear();
290 outputValues.clear();
291
292 tape.setActive();
293 }
294
297 tape.setPassive();
298 }
299
311 virtual void evalPrimal(Real const* x, Real* y = nullptr) = 0;
312
320 CODI_INLINE void evalForward(Gradient const* x_d, Gradient* y_d) {
322
323 for (size_t j = 0; j < inputValues.size(); j += 1) {
324 tape.setGradient(inputValues[j], x_d[j]);
325 }
326
327 tape.evaluateForward();
328
329 for (size_t i = 0; i < outputValues.size(); i += 1) {
330 y_d[i] = tape.getGradient(outputValues[i]);
331 tape.setGradient(outputValues[i], Gradient());
332 }
333 }
334
350 CODI_INLINE void evalForwardAt(Real const* x, Gradient const* x_d, Gradient* y_d, Real* y = nullptr) {
351 evalPrimal(x, y);
352
353 evalForward(x_d, y_d);
354 }
355
363 CODI_INLINE void evalReverse(Gradient const* y_b, Gradient* x_b) {
365
366 for (size_t i = 0; i < outputValues.size(); i += 1) {
367 tape.setGradient(outputValues[i], y_b[i]);
368 }
369
370 tape.evaluate();
371
372 for (size_t j = 0; j < inputValues.size(); j += 1) {
373 x_b[j] = tape.getGradient(inputValues[j]);
374 tape.setGradient(inputValues[j], Gradient());
375 }
376
378 tape.clearAdjoints();
379 }
380 }
381
397 CODI_INLINE void evalReverseAt(Real const* x, Gradient const* y_b, Gradient* x_b, Real* y = nullptr) {
398 evalPrimal(x, y);
399
400 evalReverse(y_b, x_b);
401 }
402
414
415 evalJacobianGen(wrapper);
416 }
417
431 CODI_INLINE void evalJacobianAt(Real const* x, JacobianType& jac, Real* y = nullptr) {
432 evalPrimal(x, y);
433
434 evalJacobian(jac);
435 }
436
444 template<typename Jac>
446 using Algo = Algorithms<Type>;
447 typename Algo::EvaluationType evalType = Algo::getEvaluationChoice(inputValues.size(), outputValues.size());
448
449 if (Algo::EvaluationType::Forward == evalType) {
451 } else if (Algo::EvaluationType::Reverse == evalType) {
453 } else {
454 CODI_EXCEPTION("Evaluation type not implemented.");
455 }
456
457 Algorithms<Type>::template computeJacobian<Jac, false>(tape, tape.getZeroPosition(), tape.getPosition(),
458 inputValues.data(), inputValues.size(),
459 outputValues.data(), outputValues.size(), jac);
460 }
461
475 template<typename Jac = DummyJacobian>
477
495 template<typename Jac = DummyJacobian>
496 CODI_INLINE void evalHessianAt(Real const* x, HessianType& hes, Real* y = nullptr,
498 evalPrimal(x, y);
499
500 cast().evalHessian(hes, jac);
501 }
502
503 protected:
504
507 return static_cast<Impl&>(*this);
508 }
509
512 wasForwardEvaluated = true;
513
514 // No cleanup to do.
515 }
516
520 // Forward evaluation leaves the adjoint vector dirty.
521
522 tape.clearAdjoints();
523 }
524
525 wasForwardEvaluated = false;
526 }
527 };
528
529 // clang-format off
535 template<typename T_Type>
536 struct TapeHelperNoImpl : public TapeHelperBase<T_Type, TapeHelperNoImpl<T_Type>> {
537 public:
538
539 CODI_STATIC_ASSERT(false && std::is_void<T_Type>::value, "Tape helper not implemented for this tape.");
540
541 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
542 using Real = typename Type::Real;
543
545
547 virtual void evalPrimal(Real const* x, Real* y = nullptr) CODI_DD(= 0;, {})
548
550 template<typename Jac = DummyJacobian>
552 };
553 // clang-format on
554
560 template<typename T_Type>
561 struct TapeHelperJacobi : public TapeHelperBase<T_Type, TapeHelperJacobi<T_Type>> {
562 public:
563
564 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
565 using Real = typename Type::Real;
566
568
570 virtual void evalPrimal(Real const* x, Real* y = nullptr) {
571 CODI_UNUSED(x, y);
572
573 CODI_EXCEPTION(
574 "No primal evaluation for Jacobian tapes. "
575 "Please use codi::RealReversePrimal or codi::RealReversePrimalIndex types for this kind of functionality.");
576 }
577
579 template<typename Jac = DummyJacobian>
581 CODI_UNUSED(hes, jac);
582
583 CODI_EXCEPTION(
584 "No direct Hessian evaluation for Jacobian tapes. "
585 "Please use codi::RealReversePrimal or codi::RealReversePrimalIndex types for this kind of functionality "
586 "or the EvaluationHelper class.");
587 }
588 };
589
595 template<typename T_Type>
596 struct TapeHelperPrimal : public TapeHelperBase<T_Type, TapeHelperPrimal<T_Type>> {
597 public:
598
599 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
600 using Real = typename Type::Real;
601
603
605 virtual void evalPrimal(Real const* x, Real* y = nullptr) {
606 for (size_t j = 0; j < this->inputValues.size(); j += 1) {
607 this->tape.primal(this->inputValues[j]) = x[j];
608 }
609
610 this->tape.evaluatePrimal();
611
612 if (nullptr != y) {
613 for (size_t i = 0; i < this->outputValues.size(); i += 1) {
614 y[i] = this->tape.primal(this->outputValues[i]);
615 }
616 }
617 }
618
620 template<typename Jac = DummyJacobian>
622 using Algo = Algorithms<Type>;
623 typename Algo::EvaluationType evalType =
624 Algo::getEvaluationChoice(this->inputValues.size(), this->outputValues.size());
625
626 if (Algo::EvaluationType::Forward == evalType) {
628 } else if (Algo::EvaluationType::Reverse == evalType) {
630 } else {
631 CODI_EXCEPTION("Evaluation type not implemented.");
632 }
633
635 this->tape, this->tape.getZeroPosition(), this->tape.getPosition(), this->inputValues.data(),
636 this->inputValues.size(), this->outputValues.data(), this->outputValues.size(), hes, jac);
637 }
638 };
639
641 template<typename Type, typename = void>
642 struct TapeHelper : public TapeHelperNoImpl<Type> {};
643
644#ifndef DOXYGEN_DISABLE
646 template<typename Type>
647 struct TapeHelper<Type, TapeTraits::EnableIfJacobianTape<typename Type::Tape>> : public TapeHelperJacobi<Type> {};
648
650 template<typename Type>
651 struct TapeHelper<Type, TapeTraits::EnableIfPrimalValueTape<typename Type::Tape>> : public TapeHelperPrimal<Type> {};
652#endif
653
654}
#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_STATIC_ASSERT(cond, message)
Static assert in CoDiPack.
Definition macros.hpp:123
bool constexpr ReversalZeroesAdjoints
With a linear index management, control if adjoints are set to zero during reversal.
Definition config.h:289
typename TraitsImplementation< Type >::PassiveReal PassiveReal
The original computation type, that was used in the application.
Definition realTraits.hpp:117
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
Basic algorithms for tape evaluation in CoDiPack.
Definition algorithms.hpp:68
static void computeHessianPrimalValueTape(Tape &tape, Position const &start, Position const &end, Identifier const *input, size_t const inputSize, Identifier const *output, size_t const outputSize, Hes &hes, Jac &jac=StaticDummy< DummyJacobian >::dummy)
Compute the Hessian with multiple tape sweeps.
Definition algorithms.hpp:244
Default implementation of the Hessian interface.
Definition hessian.hpp:59
Wrapper for JacboianInterfaces that requires a passive value conversion.
Definition jacobian.hpp:186
Default implementation of the Jacobian interface.
Definition jacobian.hpp:60
Static dummy objects for e.g. default reference arguments.
Definition staticDummy.hpp:42
A helper class that allows for a simpler access and management of a CoDiPack tape.
Definition tapeHelper.hpp:103
Gradient * createGradientVectorInput()
Create a gradient vector that can hold the tangent/adjoint of the input variables.
Definition tapeHelper.hpp:146
void evalForward(Gradient const *x_d, Gradient *y_d)
Perform a forward (tangent) evaluation of the recorded tape.
Definition tapeHelper.hpp:320
virtual ~TapeHelperBase()
Destructor.
Definition tapeHelper.hpp:136
void evalHessian(HessianType &hes, Jac &jac=StaticDummy< DummyJacobian >::dummy)
Evaluates the full Hessian of the recorded tape.
JacobianType & createJacobian()
Create a Jacobian that can hold the Jacobian of the recorded tape.
Definition tapeHelper.hpp:170
void deleteGradientVector(Gradient *vec)
Delete a gradient vector that was created with createGradientVectorInput or createGradientVectorOutpu...
Definition tapeHelper.hpp:215
Tape & tape
Reference to the global tape.
Definition tapeHelper.hpp:123
void changeStateToForwardEvaluation()
Change state.
Definition tapeHelper.hpp:511
void registerOutput(Type &value)
Add an output variable to the tape.
Definition tapeHelper.hpp:281
void evalJacobianAt(Real const *x, JacobianType &jac, Real *y=nullptr)
Re-evaluate the tape with new input variables and compute the full Jacobian at the new inputs.
Definition tapeHelper.hpp:431
void deletePrimalVector(Real *vec)
Delete a primal vector that was created with createPrimalVectorInput or createPrimalVectorOutput.
Definition tapeHelper.hpp:236
size_t getOutputSize()
Definition tapeHelper.hpp:250
void deleteHessian(HessianType &hes)
Delete the Hessian that was created with createHessian function.
Definition tapeHelper.hpp:229
void evalReverseAt(Real const *x, Gradient const *y_b, Gradient *x_b, Real *y=nullptr)
Re-evaluate the tape with new input variables and compute the AD forward mode.
Definition tapeHelper.hpp:397
void changeStateToReverseEvaluation()
Change state and clear the adjoints.
Definition tapeHelper.hpp:518
void evalForwardAt(Real const *x, Gradient const *x_d, Gradient *y_d, Real *y=nullptr)
Re-evaluate the tape with new input variables and compute the AD forward mode.
Definition tapeHelper.hpp:350
void evalJacobian(JacobianType &jac)
Evaluates the full Jacobian of the recorded tape.
Definition tapeHelper.hpp:412
void stopRecording()
Stop the recording process.
Definition tapeHelper.hpp:296
void evalJacobianGen(Jac &jac)
Evaluates the full Jacobian of the recorded tape with a custom Jacobian type chosen by the user.
Definition tapeHelper.hpp:445
Jacobian< PassiveReal > JacobianType
Type of the Jacobian.
Definition tapeHelper.hpp:116
T_Type Type
See TapeHelperBase.
Definition tapeHelper.hpp:107
typename Type::Tape Tape
Underlying tape type.
Definition tapeHelper.hpp:121
void registerInput(Type &value)
Add an input variable to the tape.
Definition tapeHelper.hpp:265
Real * createPrimalVectorInput()
Create a primal vector that can hold the primal seeding of the input variables.
Definition tapeHelper.hpp:198
Impl & cast()
Cast to the implementing class.
Definition tapeHelper.hpp:506
std::vector< Identifier > inputValues
Input value identifiers.
Definition tapeHelper.hpp:125
Gradient * createGradientVectorOutput()
Create a gradient vector that can hold the tangent/adjoint of the output variables.
Definition tapeHelper.hpp:158
typename Type::Real Real
See LhsExpressionInterface.
Definition tapeHelper.hpp:110
T_Impl Impl
See TapeHelperBase.
Definition tapeHelper.hpp:108
Hessian< PassiveReal > HessianType
Type of the Hessian.
Definition tapeHelper.hpp:117
typename Type::Identifier Identifier
See LhsExpressionInterface.
Definition tapeHelper.hpp:111
void evalHessianAt(Real const *x, HessianType &hes, Real *y=nullptr, Jac &jac=StaticDummy< DummyJacobian >::dummy)
Re-evaluate the tape with new input variables and compute the full Hessian at the new inputs.
Definition tapeHelper.hpp:496
typename Type::Gradient Gradient
See LhsExpressionInterface.
Definition tapeHelper.hpp:112
void startRecording()
Start the recording process. Deletes the old tape.
Definition tapeHelper.hpp:287
size_t getInputSize()
Definition tapeHelper.hpp:244
typename RealTraits::PassiveReal< Real > PassiveReal
Passive base of the CoDiPack type.
Definition tapeHelper.hpp:114
Real * createPrimalVectorOutput()
Create a primal vector that can hold the primal result of the output variables.
Definition tapeHelper.hpp:210
void evalReverse(Gradient const *y_b, Gradient *x_b)
Perform a reverse (adjoint) evaluation of the recorded tape.
Definition tapeHelper.hpp:363
std::vector< Identifier > outputValues
Input value identifiers.
Definition tapeHelper.hpp:126
virtual void evalPrimal(Real const *x, Real *y=nullptr)=0
Perform a primal re-evaluation of the tape.
bool wasForwardEvaluated
State of the last evaluation.
Definition tapeHelper.hpp:128
TapeHelperBase()
Constructor.
Definition tapeHelper.hpp:133
HessianType & createHessian()
Create a Hessian that can hold the Hessian of the recorded tape.
Definition tapeHelper.hpp:184
void deleteJacobian(JacobianType &jac)
Delete the Jacobian that was created with createJacobian function.
Definition tapeHelper.hpp:222
Definition tapeHelper.hpp:561
virtual void evalPrimal(Real const *x, Real *y=nullptr)
Throws an exception since primal evaluations are not support by Jacobian tapes.
Definition tapeHelper.hpp:570
void evalHessian(typename Base::HessianType &hes, Jac &jac=StaticDummy< DummyJacobian >::dummy)
Throws an exception since primal evaluations are not supported by Jacobian tapes.
Definition tapeHelper.hpp:580
typename Type::Real Real
See TapeHelperBase.
Definition tapeHelper.hpp:565
T_Type Type
See TapeHelperBase.
Definition tapeHelper.hpp:564
Definition tapeHelper.hpp:536
virtual void evalPrimal(Real const *x, Real *y=nullptr)=0
Missing implementation will yield linker errors.
typename Type::Real Real
See TapeHelperBase.
Definition tapeHelper.hpp:542
T_Type Type
See TapeHelperBase.
Definition tapeHelper.hpp:541
void evalHessian(typename Base::HessianType &hes, Jac &jac=StaticDummy< DummyJacobian >::dummy)
Missing implementation will yield linker errors.
Definition tapeHelper.hpp:596
virtual void evalPrimal(Real const *x, Real *y=nullptr)
Perform a primal re-evaluation of the tape.
Definition tapeHelper.hpp:605
T_Type Type
See TapeHelperBase.
Definition tapeHelper.hpp:599
typename Type::Real Real
See TapeHelperBase.
Definition tapeHelper.hpp:600
void evalHessian(typename Base::HessianType &hes, Jac &jac=StaticDummy< DummyJacobian >::dummy)
Evaluates the full Hessian of the recorded tape.
Definition tapeHelper.hpp:621
See TapeHelperBase.
Definition tapeHelper.hpp:642