jacobianBaseTape.hpp

/*
 * CoDiPack, a Code Differentiation Package
 *
 * Copyright (C) 2015-2024 Chair for Scientific Computing (SciComp), University of Kaiserslautern-Landau
 * Homepage: http://scicomp.rptu.de
 * Contact:  Prof. Nicolas R. Gauger (codi@scicomp.uni-kl.de)
 *
 * Lead developers: Max Sagebaum, Johannes Blühdorn (SciComp, University of Kaiserslautern-Landau)
 *
 * This file is part of CoDiPack (http://scicomp.rptu.de/software/codi).
 *
 * CoDiPack is free software: you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * CoDiPack is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty
 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 * See the GNU General Public License for more details.
 * You should have received a copy of the GNU
 * General Public License along with CoDiPack.
 * If not, see <http://www.gnu.org/licenses/>.
 *
 * For other licensing options please contact us.
 *
 * Authors:
 *  - SciComp, University of Kaiserslautern-Landau:
 *    - Max Sagebaum
 *    - Johannes Blühdorn
 *    - Former members:
 *      - Tim Albring
 */
#pragma once
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <type_traits>
 
#include "../config.h"
#include "../expressions/lhsExpressionInterface.hpp"
#include "../expressions/logic/compileTimeTraversalLogic.hpp"
#include "../expressions/logic/helpers/forEachLeafLogic.hpp"
#include "../expressions/logic/helpers/jacobianComputationLogic.hpp"
#include "../expressions/logic/traversalLogic.hpp"
#include "../expressions/referenceActiveType.hpp"
#include "../misc/macros.hpp"
#include "../misc/mathUtility.hpp"
#include "../misc/memberStore.hpp"
#include "../traits/adjointVectorTraits.hpp"
#include "../traits/computationTraits.hpp"
#include "../traits/expressionTraits.hpp"
#include "commonTapeImplementation.hpp"
#include "data/chunk.hpp"
#include "data/chunkedData.hpp"
#include "indices/indexManagerInterface.hpp"
#include "io/tapeReaderWriterInterface.hpp"
#include "misc/adjointVectorAccess.hpp"
#include "misc/duplicateJacobianRemover.hpp"
#include "misc/localAdjoints.hpp"
 
namespace codi {
 
  template<typename T_Real, typename T_Gradient, typename T_IndexManager, template<typename, typename> class T_Data,
           template<typename, typename, typename> class T_Adjoints = LocalAdjoints>
  struct JacobianTapeTypes : public TapeTypesInterface {
    public:
 
      using Real = CODI_DD(T_Real, double);                                              
      using Gradient = CODI_DD(T_Gradient, double);                                      
      using IndexManager = CODI_DD(T_IndexManager, CODI_T(IndexManagerInterface<int>));  
      template<typename Chunk, typename Nested>
      using Data = CODI_DD(CODI_T(T_Data<Chunk, Nested>),
                           CODI_T(DataInterface<Nested>));  
 
      using Identifier = typename IndexManager::Index;  
 
      template<typename Impl>
      using Adjoints = CODI_DD(CODI_T(T_Adjoints<Gradient, Identifier, Impl>),
                               CODI_T(InternalAdjointsInterface<Gradient, Identifier, CODI_ANY>));
 
      static bool constexpr IsLinearIndexHandler = IndexManager::IsLinear;  
      static bool constexpr IsStaticIndexHandler =
          IndexManager::NeedsStaticStorage;  
 
      using StatementChunk = typename std::conditional<IsLinearIndexHandler, Chunk1<Config::ArgumentSize>,
                                                       Chunk2<Identifier, Config::ArgumentSize> >::type;
      using StatementData = Data<StatementChunk, IndexManager>;  
 
      using JacobianChunk = Chunk2<Real, Identifier>;           
      using JacobianData = Data<JacobianChunk, StatementData>;  
 
      using NestedData = JacobianData;  
  };
 
  template<typename T_TapeTypes, typename T_Impl>
  struct JacobianBaseTape : public CommonTapeImplementation<T_TapeTypes, T_Impl> {
    public:
 
      using TapeTypes = CODI_DD(
          T_TapeTypes, CODI_T(JacobianTapeTypes<double, double, IndexManagerInterface<int>, DefaultChunkedData>));
      using Impl = CODI_DD(T_Impl, CODI_DEFAULT_TAPE);  
 
      using Base = CommonTapeImplementation<T_TapeTypes, T_Impl>;  
      friend Base;  
 
      using Real = typename TapeTypes::Real;                  
      using Gradient = typename TapeTypes::Gradient;          
      using IndexManager = typename TapeTypes::IndexManager;  
      using Identifier = typename TapeTypes::Identifier;      
 
      using StatementData = typename TapeTypes::StatementData;  
      using JacobianData = typename TapeTypes::JacobianData;    
 
      using Adjoints = typename TapeTypes::template Adjoints<Impl>;  
 
      using PassiveReal = RealTraits::PassiveReal<Real>;  
 
      using NestedPosition = typename JacobianData::Position;  
      using Position = typename Base::Position;                
 
      template<typename AdjointVector>
      using VectorAccess = AdjointVectorAccess<Real, Identifier, AdjointVector>;
 
      static bool constexpr AllowJacobianOptimization = true;  
      static bool constexpr HasPrimalValues = false;           
 
      static bool constexpr LinearIndexHandling =
          TapeTypes::IsLinearIndexHandler;                  
      static bool constexpr RequiresPrimalRestore = false;  
 
    protected:
 
#if CODI_RemoveDuplicateJacobianArguments
      DuplicateJacobianRemover<Real, Identifier> jacobianSorter;  
#endif
 
      MemberStore<IndexManager, Impl, TapeTypes::IsStaticIndexHandler> indexManager;  
      StatementData statementData;  
      JacobianData jacobianData;    
 
      Adjoints adjoints;  
 
    private:
 
      CODI_INLINE Impl const& cast() const {
        return static_cast<Impl const&>(*this);
      }
 
      CODI_INLINE Impl& cast() {
        return static_cast<Impl&>(*this);
      }
 
    protected:
 
      /*******************************************************************************/
 
      template<typename... Args>
      static void internalEvaluateForward_EvalStatements(Args&&... args);
 
      template<typename... Args>
      static void internalEvaluateReverse_EvalStatements(Args&&... args);
 
      void pushStmtData(Identifier const& index, Config::ArgumentSize const& numberOfArguments);
 
 
    public:
 
      JacobianBaseTape()
          : Base(),
#if CODI_RemoveDuplicateJacobianArguments
            jacobianSorter(),
#endif
            indexManager(0),  // Reserve the zero index.
            statementData(Config::ChunkSize),
            jacobianData(std::max(Config::ChunkSize, Config::MaxArgumentSize)),  // Chunk must be large enough to store
                                                                                 // data for all arguments of one
                                                                                 // statement.
            adjoints(1)  // Ensure that adjoint[0] exists, see its use in gradient() const.
 
      {
        statementData.setNested(&indexManager.get());
        jacobianData.setNested(&statementData);
 
        Base::init(&jacobianData);
 
        Base::options.insert(TapeParameters::AdjointSize);
        Base::options.insert(TapeParameters::JacobianSize);
        Base::options.insert(TapeParameters::LargestIdentifier);
        Base::options.insert(TapeParameters::StatementSize);
      }
 
      /*******************************************************************************/
 
      CODI_INLINE Gradient& gradient(Identifier const& identifier,
                                     AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(identifier);
        }
 
        codiAssert(identifier < (Identifier)adjoints.size());
 
        return adjoints[identifier];
      }
 
      CODI_INLINE Gradient const& gradient(
          Identifier const& identifier, AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) const {
        codiAssert(identifier < (Identifier)adjoints.size());
 
        if (AdjointsManagement::Automatic == adjointsManagement && identifier >= (Identifier)adjoints.size()) {
          return adjoints[0];
        } else {
          return adjoints[identifier];
        }
      }
 
      /*******************************************************************************/
 
      template<typename Real>
      CODI_INLINE void initIdentifier(Real& value, Identifier& identifier) {
        CODI_UNUSED(value);
 
        identifier = IndexManager::InactiveIndex;
      }
 
      template<typename Real>
      CODI_INLINE void destroyIdentifier(Real& value, Identifier& identifier) {
        CODI_UNUSED(value);
 
        indexManager.get().template freeIndex<Impl>(identifier);
      }
 
 
    protected:
 
      struct PushJacobianLogic : public JacobianComputationLogic<PushJacobianLogic> {
        public:
          template<typename Node, typename Jacobian, typename DataVector>
          CODI_INLINE void handleJacobianOnActive(Node const& node, Jacobian jacobianExpr, DataVector& dataVector) {
            Real jacobian = ComputationTraits::adjointConversion<Real>(jacobianExpr);
 
            if (CODI_ENABLE_CHECK(Config::CheckZeroIndex, 0 != node.getIdentifier())) {
              if (CODI_ENABLE_CHECK(Config::IgnoreInvalidJacobians, RealTraits::isTotalFinite(jacobian))) {
                if (CODI_ENABLE_CHECK(Config::CheckJacobianIsZero, !RealTraits::isTotalZero(jacobian))) {
                  dataVector.pushData(jacobian, node.getIdentifier());
                }
              }
            }
          }
 
          template<typename Type, typename Jacobian, typename DataVector>
          CODI_INLINE void handleJacobianOnActive(ReferenceActiveType<Type> const& node, Jacobian jacobianExpr,
                                                  DataVector& dataVector) {
            CODI_UNUSED(dataVector);
 
            Real jacobian = ComputationTraits::adjointConversion<Real>(jacobianExpr);
 
            if (CODI_ENABLE_CHECK(Config::IgnoreInvalidJacobians, RealTraits::isTotalFinite(jacobian))) {
              // Do a delayed push for these leaf nodes, accumulate the jacobian in the local member.
              node.jacobian += jacobian;
            }
          }
      };
 
      struct PushDelayedJacobianLogic : public ForEachLeafLogic<PushDelayedJacobianLogic> {
        public:
 
          template<typename Type, typename DataVector>
          CODI_INLINE void handleActive(ReferenceActiveType<Type> const& node, DataVector& dataVector) {
            if (CODI_ENABLE_CHECK(Config::CheckZeroIndex, 0 != node.getIdentifier())) {
              if (CODI_ENABLE_CHECK(Config::CheckJacobianIsZero, !RealTraits::isTotalZero(node.jacobian))) {
                dataVector.pushData(node.jacobian, node.getIdentifier());
 
                // Reset the jacobian here such that it is not pushed multiple times and ready for the next store.
                node.jacobian = Real();
              }
            }
          }
 
          using ForEachLeafLogic<PushDelayedJacobianLogic>::handleActive;
      };
 
      template<typename Rhs>
      CODI_INLINE void pushJacobians(ExpressionInterface<Real, Rhs> const& rhs) {
        PushJacobianLogic pushJacobianLogic;
        PushDelayedJacobianLogic pushDelayedJacobianLogic;
 
#if CODI_RemoveDuplicateJacobianArguments
        auto& insertVector = jacobianSorter;
#else
        auto& insertVector = jacobianData;
#endif
 
        pushJacobianLogic.eval(rhs.cast(), Real(1.0), insertVector);
        pushDelayedJacobianLogic.eval(rhs.cast(), insertVector);
 
#if CODI_RemoveDuplicateJacobianArguments
        jacobianSorter.storeData(jacobianData);
#endif
      }
 
    public:
 
 
      template<typename Lhs, typename Rhs>
      CODI_INLINE void store(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& lhs,
                             ExpressionInterface<Real, Rhs> const& rhs) {
        if (CODI_ENABLE_CHECK(Config::CheckTapeActivity, cast().isActive())) {
          size_t constexpr MaxArgs = ExpressionTraits::NumberOfActiveTypeArguments<Rhs>::value;
 
          codiAssert(MaxArgs < Config::MaxArgumentSize);
 
          statementData.reserveItems(1);
          typename JacobianData::InternalPosHandle jacobianStart = jacobianData.reserveItems(MaxArgs);
 
          pushJacobians(rhs);
 
          size_t numberOfArguments = jacobianData.getPushedDataCount(jacobianStart);
          if (CODI_ENABLE_CHECK(Config::CheckEmptyStatements, 0 != numberOfArguments)) {
            indexManager.get().template assignIndex<Impl>(lhs.cast().getIdentifier());
            cast().pushStmtData(lhs.cast().getIdentifier(), (Config::ArgumentSize)numberOfArguments);
 
            if (Config::StatementEvents) {
              Real* jacobians;
              Identifier* rhsIdentifiers;
              jacobianData.getDataPointers(jacobians, rhsIdentifiers);
              jacobians -= numberOfArguments;
              rhsIdentifiers -= numberOfArguments;
 
              EventSystem<Impl>::notifyStatementStoreOnTapeListeners(cast(), lhs.cast().getIdentifier(),
                                                                     rhs.cast().getValue(), numberOfArguments,
                                                                     rhsIdentifiers, jacobians);
            }
          } else {
            indexManager.get().template freeIndex<Impl>(lhs.cast().getIdentifier());
          }
        } else {
          indexManager.get().template freeIndex<Impl>(lhs.cast().getIdentifier());
        }
 
        lhs.cast().value() = rhs.cast().getValue();
      }
 
      template<typename Lhs, typename Rhs>
      CODI_INLINE void store(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& lhs,
                             LhsExpressionInterface<Real, Gradient, Impl, Rhs> const& rhs) {
        if (CODI_ENABLE_CHECK(Config::CheckTapeActivity, cast().isActive())) {
          if (IndexManager::CopyNeedsStatement || !Config::CopyOptimization) {
            store<Lhs, Rhs>(lhs, static_cast<ExpressionInterface<Real, Rhs> const&>(rhs));
            return;
          } else {
            indexManager.get().template copyIndex<Impl>(lhs.cast().getIdentifier(), rhs.cast().getIdentifier());
          }
        } else {
          indexManager.get().template freeIndex<Impl>(lhs.cast().getIdentifier());
        }
 
        lhs.cast().value() = rhs.cast().getValue();
      }
 
      template<typename Lhs>
      CODI_INLINE void store(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& lhs, Real const& rhs) {
        indexManager.get().template freeIndex<Impl>(lhs.cast().getIdentifier());
 
        lhs.cast().value() = rhs;
      }
 
      /*******************************************************************************
       * Protected helper function for ReverseTapeInterface
       */
 
    protected:
 
      template<typename Lhs>
      CODI_INLINE void internalRegisterInput(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& value,
                                             bool unusedIndex) {
        if (TapeTypes::IsLinearIndexHandler) {
          statementData.reserveItems(1);
        }
 
        if (unusedIndex) {
          indexManager.get().template assignUnusedIndex<Impl>(value.cast().getIdentifier());
        } else {
          indexManager.get().template assignIndex<Impl>(value.cast().getIdentifier());
        }
 
        if (TapeTypes::IsLinearIndexHandler) {
          cast().pushStmtData(value.cast().getIdentifier(), Config::StatementInputTag);
        }
      }
 
    public:
 
      /*******************************************************************************/
 
      template<typename Lhs>
      CODI_INLINE void registerInput(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& value) {
        cast().internalRegisterInput(value, true);
        EventSystem<Impl>::notifyTapeRegisterInputListeners(cast(), value.cast().value(), value.cast().getIdentifier());
      }
 
      CODI_INLINE void clearAdjoints(AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          adjoints.beginUse();
        }
 
        adjoints.zeroAll(indexManager.get().getLargestCreatedIndex());
 
        if (AdjointsManagement::Automatic == adjointsManagement) {
          adjoints.endUse();
        }
      }
 
 
    protected:
 
      CODI_INLINE TapeValues internalGetTapeValues() const {
        std::string name;
        if (TapeTypes::IsLinearIndexHandler) {
          name = "CoDi Tape Statistics ( JacobianLinearTape )";
        } else {
          name = "CoDi Tape Statistics ( JacobianReuseTape )";
        }
        TapeValues values = TapeValues(name);
 
        size_t nAdjoints = indexManager.get().getLargestCreatedIndex();
        double memoryAdjoints = static_cast<double>(nAdjoints) * static_cast<double>(sizeof(Gradient));
 
        bool constexpr globalAdjoints = InternalAdjointVectorTraits::IsGlobal<Adjoints>::value;
        TapeValues::LocalReductionOperation constexpr operation =
            globalAdjoints ? TapeValues::LocalReductionOperation::Max : TapeValues::LocalReductionOperation::Sum;
 
        values.addSection("Adjoint vector");
        values.addUnsignedLongEntry("Number of adjoints", nAdjoints, operation);
        values.addDoubleEntry("Memory allocated", memoryAdjoints, operation, true, true);
 
        values.addSection("Index manager");
        indexManager.get().addToTapeValues(values);
 
        values.addSection("Statement entries");
        statementData.addToTapeValues(values);
        values.addSection("Jacobian entries");
        jacobianData.addToTapeValues(values);
 
        return values;
      }
 
      /******************************************************************************
       * Protected helper function for CustomAdjointVectorEvaluationTapeInterface
       */
 
    protected:
 
      template<typename AdjointVector>
      CODI_INLINE static void incrementAdjoints(AdjointVector& adjointVector,
                                                AdjointVectorTraits::Gradient<AdjointVector> const& lhsAdjoint,
                                                Config::ArgumentSize const& numberOfArguments, size_t& curJacobianPos,
                                                Real const* const rhsJacobians,
                                                Identifier const* const rhsIdentifiers) {
        size_t endJacobianPos = curJacobianPos - numberOfArguments;
 
        if (CODI_ENABLE_CHECK(Config::SkipZeroAdjointEvaluation, !RealTraits::isTotalZero(lhsAdjoint))) CODI_Likely {
          while (endJacobianPos < curJacobianPos) CODI_Likely {
            curJacobianPos -= 1;
            adjointVector[rhsIdentifiers[curJacobianPos]] += rhsJacobians[curJacobianPos] * lhsAdjoint;
          }
        } else CODI_Unlikely {
          curJacobianPos = endJacobianPos;
        }
      }
 
      CODI_WRAP_FUNCTION_TEMPLATE(Wrap_internalEvaluateReverse_EvalStatements,
                                  Impl::template internalEvaluateReverse_EvalStatements);
 
      template<typename AdjointVector>
      CODI_INLINE static void incrementTangents(AdjointVector const& adjointVector,
                                                AdjointVectorTraits::Gradient<AdjointVector>& lhsAdjoint,
                                                Config::ArgumentSize const& numberOfArguments, size_t& curJacobianPos,
                                                Real const* const rhsJacobians,
                                                Identifier const* const rhsIdentifiers) {
        size_t endJacobianPos = curJacobianPos + numberOfArguments;
 
        while (curJacobianPos < endJacobianPos) CODI_Likely {
          lhsAdjoint += rhsJacobians[curJacobianPos] * adjointVector[rhsIdentifiers[curJacobianPos]];
          curJacobianPos += 1;
        }
      }
 
      CODI_WRAP_FUNCTION_TEMPLATE(Wrap_internalEvaluateForward_EvalStatements,
                                  Impl::template internalEvaluateForward_EvalStatements);
 
    public:
 
 
      using Base::evaluate;
 
      template<typename AdjointVector>
      CODI_NO_INLINE void evaluate(Position const& start, Position const& end, AdjointVector&& data) {
        VectorAccess<AdjointVector> adjointWrapper(data);
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(
            cast(), start, end, &adjointWrapper, EventHints::EvaluationKind::Reverse, EventHints::Endpoint::Begin);
 
        Wrap_internalEvaluateReverse_EvalStatements<AdjointVector> evalFunc;
        Base::llfByteData.evaluateReverse(start, end, evalFunc, cast(), std::forward<AdjointVector>(data));
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(cast(), start, end, &adjointWrapper,
                                                       EventHints::EvaluationKind::Reverse, EventHints::Endpoint::End);
      }
 
      template<typename AdjointVector>
      CODI_NO_INLINE void evaluateForward(Position const& start, Position const& end, AdjointVector&& data) {
        VectorAccess<AdjointVector> adjointWrapper(data);
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(
            cast(), start, end, &adjointWrapper, EventHints::EvaluationKind::Forward, EventHints::Endpoint::Begin);
 
        Wrap_internalEvaluateForward_EvalStatements<AdjointVector> evalFunc;
        Base::llfByteData.evaluateForward(start, end, evalFunc, cast(), std::forward<AdjointVector>(data));
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(cast(), start, end, &adjointWrapper,
                                                       EventHints::EvaluationKind::Forward, EventHints::Endpoint::End);
      }
 
      CODI_INLINE decltype(adjoints.data()) getInternalAdjoints() {
        return adjoints.data();
      }
 
      /*******************************************************************************/
 
      CODI_INLINE void swap(Impl& other) {
        // Index manager does not need to be swapped, it is either static or swapped with the vector data.
        // Vectors are swapped recursively in the base class.
 
        adjoints.swap(other.adjoints);
 
        Base::swap(other);
      }
 
      void deleteAdjointVector() {
        adjoints.resize(1);
      }
 
      void resizeAdjointVector() {
        checkAdjointSize(indexManager.get().getLargestCreatedIndex());
      }
 
      void beginUseAdjointVector() {
        adjoints.beginUse();
      }
 
      void endUseAdjointVector() {
        adjoints.endUse();
      }
 
      size_t getParameter(TapeParameters parameter) const {
        switch (parameter) {
          case TapeParameters::AdjointSize:
            return adjoints.size();
            break;
          case TapeParameters::JacobianSize:
            return jacobianData.getDataSize();
            break;
          case TapeParameters::LargestIdentifier:
            return indexManager.get().getLargestCreatedIndex();
            break;
          case TapeParameters::StatementSize:
            return statementData.getDataSize();
          default:
            return Base::getParameter(parameter);
            break;
        }
      }
 
      void setParameter(TapeParameters parameter, size_t value) {
        switch (parameter) {
          case TapeParameters::AdjointSize:
            adjoints.resize(value);
            break;
          case TapeParameters::JacobianSize:
            jacobianData.resize(value);
            break;
          case TapeParameters::LargestIdentifier:
            CODI_EXCEPTION("Tried to set a get only parameter.");
            break;
          case TapeParameters::StatementSize:
            statementData.resize(value);
            break;
          default:
            Base::setParameter(parameter, value);
            break;
        }
      }
 
      VectorAccess<decltype(adjoints.data())>* createVectorAccess() {
        return createVectorAccessCustomAdjoints(adjoints.data());
      }
 
      template<typename AdjointVector>
      VectorAccess<AdjointVector>* createVectorAccessCustomAdjoints(AdjointVector&& data) {
        return new VectorAccess<AdjointVector>(data);
      }
 
      template<typename Adjoint>
      VectorAccess<Adjoint*>* createVectorAccessCustomAdjoints(Adjoint* data) {
        return new VectorAccess<Adjoint*>(data);
      }
 
      void deleteVectorAccess(VectorAccessInterface<Real, Identifier>* access) {
        delete access;
      }
 
      /*******************************************************************************/
 
      template<typename Lhs>
      Real registerExternalFunctionOutput(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& value) {
        cast().internalRegisterInput(value, false);
 
        return Real();
      }
 
      /*******************************************************************************/
 
      using Base::evaluateForward;
 
      void evaluateForward(Position const& start, Position const& end,
                           AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(indexManager.get().getLargestCreatedIndex());
          adjoints.beginUse();
        }
 
        codiAssert(indexManager.get().getLargestCreatedIndex() < (Identifier)adjoints.size());
 
        cast().evaluateForward(start, end, adjoints.data());
 
        if (AdjointsManagement::Automatic == adjointsManagement) {
          adjoints.endUse();
        }
      }
 
      /*******************************************************************************/
 
      void pushJacobianManual(Real const& jacobian, Real const& value, Identifier const& index) {
        CODI_UNUSED(value);
 
        cast().incrementManualPushCounter();
 
        jacobianData.pushData(jacobian, index);
 
        if (Config::StatementEvents) {
          if (this->manualPushCounter == this->manualPushGoal) {
            // emit statement event
            Real* jacobians;
            Identifier* rhsIdentifiers;
            jacobianData.getDataPointers(jacobians, rhsIdentifiers);
            jacobians -= this->manualPushGoal;
            rhsIdentifiers -= this->manualPushGoal;
 
            EventSystem<Impl>::notifyStatementStoreOnTapeListeners(cast(), this->manualPushLhsIdentifier,
                                                                   this->manualPushLhsValue, this->manualPushGoal,
                                                                   rhsIdentifiers, jacobians);
          }
        }
      }
 
      void storeManual(Real const& lhsValue, Identifier& lhsIndex, Config::ArgumentSize const& size) {
        CODI_UNUSED(lhsValue);
        Impl& impl = cast();
 
        codiAssert(size < Config::MaxArgumentSize);
 
        statementData.reserveItems(1);
        jacobianData.reserveItems(size);
 
        indexManager.get().template assignIndex<Impl>(lhsIndex);
        impl.pushStmtData(lhsIndex, (Config::ArgumentSize)size);
 
        impl.initializeManualPushData(lhsValue, lhsIndex, size);
      }
 
      /*******************************************************************************/
 
      void createStatementManual(Real const& lhsValue, Identifier& lhsIndex, Config::ArgumentSize const& size,
                                 Real const* jacobians, Identifier const* rhsIdentifiers) {
        CODI_UNUSED(lhsValue);
        Impl& impl = cast();
 
        statementData.reserveItems(1);
 
        if (Config::StatementInputTag == size && TapeTypes::IsLinearIndexHandler) CODI_Unlikely {
          impl.pushStmtData(lhsIndex, (Config::ArgumentSize)size);
        } else CODI_Likely {
          jacobianData.reserveItems(size);
          impl.pushStmtData(lhsIndex, (Config::ArgumentSize)size);
          impl.initializeManualPushData(lhsValue, lhsIndex, size);
          // Record the rhs of the statement.
          for (size_t rhsCount = 0; rhsCount < size; rhsCount++) {
            cast().incrementManualPushCounter();
            jacobianData.pushData(jacobians[rhsCount], rhsIdentifiers[rhsCount]);
          }
        }
      }
 
      using Base::writeTape;
 
      template<typename Type>
      CODI_NO_INLINE void writeTape(codi::TapeWriterInterface<Type>* writer, Position const& start,
                                    Position const& end) {
        Impl& impl = cast();
        writer->start(impl);
        Base::llfByteData.evaluateForward(start, end, Impl::template internalWriteTape<Type>,  writer);
        writer->finish();
      }
 
      /*******************************************************************************/
 
      IndexManager& getIndexManager() {
        return indexManager.get();
      }
 
      /*******************************************************************************/
 
      CODI_INLINE void pushLowLevelFunction(Config::LowLevelFunctionToken token, size_t size, ByteDataView& data) {
        statementData.reserveItems(1);
 
        Base::internalStoreLowLevelFunction(token, size, data);
 
        Identifier lhsIndex = Identifier();
        if (LinearIndexHandling) {
          indexManager.get().template assignIndex<Impl>(lhsIndex);
        }
        cast().pushStmtData(lhsIndex, Config::StatementLowLevelFunctionTag);
      }
 
      /*******************************************************************************/
 
      CODI_INLINE void evaluate(Position const& start, Position const& end,
                                AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(indexManager.get().getLargestCreatedIndex());
          adjoints.beginUse();
        }
 
        codiAssert(indexManager.get().getLargestCreatedIndex() < (Identifier)adjoints.size());
 
        evaluate(start, end, adjoints.data());
 
        if (AdjointsManagement::Automatic == adjointsManagement) {
          adjoints.endUse();
        }
      }
 
      /*******************************************************************************/
 
      void evaluateKeepState(Position const& start, Position const& end,
                             AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        evaluate(start, end, adjointsManagement);
      }
 
      template<typename AdjointVector>
      void evaluateKeepState(Position const& start, Position const& end, AdjointVector&& data) {
        evaluate(start, end, std::forward<AdjointVector>(data));
      }
 
      void evaluateForwardKeepState(Position const& start, Position const& end,
                                    AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        evaluateForward(start, end, adjointsManagement);
      }
 
      template<typename AdjointVector>
      void evaluateForwardKeepState(Position const& start, Position const& end, AdjointVector&& data) {
        evaluateForward(start, end, std::forward<AdjointVector>(data));
      }
 
      /*******************************************************************************/
 
      using Base::evaluatePrimal;
 
      void evaluatePrimal(Position const& start, Position const& end) {
        CODI_UNUSED(start, end);
 
        CODI_EXCEPTION("Accessing primal evaluation of an Jacobian tape.");
      }
 
      Real& primal(Identifier const& identifier) {
        CODI_UNUSED(identifier);
 
        CODI_EXCEPTION("Accessing primal vector of an Jacobian tape.");
 
        static Real temp;
        return temp;
      }
 
      Real const& primal(Identifier const& identifier) const {
        CODI_UNUSED(identifier);
 
        CODI_EXCEPTION("Accessing primal vector of an Jacobian tape.");
 
        static Real temp;
        return temp;
      }
 
 
    private:
 
      CODI_INLINE void checkAdjointSize(Identifier const& identifier) {
        if (identifier >= (Identifier)adjoints.size()) {
          internalResizeAdjointsVector();
        }
      }
 
      CODI_NO_INLINE void internalResizeAdjointsVector() {
        // overallocate as next multiple of Config::ChunkSize
        adjoints.resize(
            getNextMultiple(indexManager.get().getLargestCreatedIndex() + 1, (Identifier)Config::ChunkSize));
      }
  };
}