primalValueBaseTape.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 <functional>
#include <type_traits>
#include <utility>
 
#include "../config.h"
#include "../expressions/lhsExpressionInterface.hpp"
#include "../expressions/logic/compileTimeTraversalLogic.hpp"
#include "../expressions/logic/constructStaticContext.hpp"
#include "../expressions/logic/helpers/forEachLeafLogic.hpp"
#include "../expressions/logic/helpers/jacobianComputationLogic.hpp"
#include "../expressions/logic/helpers/mathStatementGenLogic.hpp"
#include "../expressions/logic/traversalLogic.hpp"
#include "../misc/demangleName.hpp"
#include "../misc/macros.hpp"
#include "../misc/mathUtility.hpp"
#include "../misc/memberStore.hpp"
#include "../traits/expressionTraits.hpp"
#include "commonTapeImplementation.hpp"
#include "data/chunk.hpp"
#include "data/chunkedData.hpp"
#include "indices/indexManagerInterface.hpp"
#include "misc/primalAdjointVectorAccess.hpp"
#include "statementEvaluators/statementEvaluatorInterface.hpp"
#include "statementEvaluators/statementEvaluatorTapeInterface.hpp"
namespace codi {
 
  template<typename T_Real, typename T_Gradient, typename T_IndexManager, template<typename> class T_StatementEvaluator,
           template<typename, typename> class T_Data>
  struct PrimalValueTapeTypes : 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>));  
      using StatementEvaluator = CODI_DD(CODI_T(T_StatementEvaluator<Real>),
                                         CODI_T(StatementEvaluatorInterface<double>));  
      template<typename Chunk, typename Nested>
      using Data = CODI_DD(CODI_T(T_Data<Chunk, Nested>),
                           CODI_T(DataInterface<Nested>));  
 
      using Identifier = typename IndexManager::Index;    
      using PassiveReal = RealTraits::PassiveReal<Real>;  
 
      constexpr static bool IsLinearIndexHandler = IndexManager::IsLinear;  
      constexpr static bool IsStaticIndexHandler =
          IndexManager::NeedsStaticStorage;  
 
      using EvalHandle = typename StatementEvaluator::Handle;  
 
      using StatementChunk =
          typename std::conditional<IsLinearIndexHandler, Chunk2<Config::ArgumentSize, EvalHandle>,
                                    Chunk4<Identifier, Config::ArgumentSize, Real, EvalHandle>>::type;
      using StatementData = Data<StatementChunk, IndexManager>;  
 
      using IdentifierChunk = Chunk1<Identifier>;                      
      using RhsIdentifierData = Data<IdentifierChunk, StatementData>;  
 
      using PassiveValueChunk = Chunk1<Real>;                               
      using PassiveValueData = Data<PassiveValueChunk, RhsIdentifierData>;  
 
      using ConstantValueChunk = Chunk1<PassiveReal>;  
      using ConstantValueData = Data<ConstantValueChunk, PassiveValueData>;  
 
      using NestedData = ConstantValueData;  
  };
 
  template<typename T_TapeTypes, typename T_Impl>
  struct PrimalValueBaseTape : public CommonTapeImplementation<T_TapeTypes, T_Impl>,
                               public StatementEvaluatorTapeInterface<typename T_TapeTypes::Real>,
                               public StatementEvaluatorInnerTapeInterface<typename T_TapeTypes::Real> {
    public:
 
      using TapeTypes = CODI_DD(T_TapeTypes,
                                CODI_T(PrimalValueTapeTypes<double, double, IndexManagerInterface<int>,
                                                            StatementEvaluatorInterface, DefaultChunkedData>));
      using Impl = CODI_DD(T_Impl, CODI_T(PrimalValueBaseTape));
 
      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 StatementEvaluator = typename TapeTypes::StatementEvaluator;  
      using Identifier = typename TapeTypes::Identifier;                  
 
      using EvalHandle = typename TapeTypes::EvalHandle;  
 
      using StatementData = typename TapeTypes::StatementData;          
      using RhsIdentifierData = typename TapeTypes::RhsIdentifierData;  
      using PassiveValueData = typename TapeTypes::PassiveValueData;    
      using ConstantValueData = typename TapeTypes::ConstantValueData;  
 
      using PassiveReal = RealTraits::PassiveReal<Real>;  
 
      using NestedPosition = typename ConstantValueData::Position;  
      using Position = typename Base::Position;                     
 
      template<typename AdjointVector>
      using VectorAccess = PrimalAdjointVectorAccess<Real, Identifier, AdjointVector>;
 
      static bool constexpr AllowJacobianOptimization = false;  
      static bool constexpr HasPrimalValues = true;             
      static bool constexpr LinearIndexHandling =
          TapeTypes::IsLinearIndexHandler;  
      static bool constexpr RequiresPrimalRestore =
          !TapeTypes::IsLinearIndexHandler;  
 
    protected:
 
      static EvalHandle const jacobianExpressions[Config::MaxArgumentSize];
 
      MemberStore<IndexManager, Impl, TapeTypes::IsStaticIndexHandler> indexManager;  
      StatementData statementData;          
      RhsIdentifierData rhsIdentifierData;  
      PassiveValueData passiveValueData;    
      ConstantValueData constantValueData;  
 
      std::vector<Gradient> adjoints;  
      std::vector<Real> primals;       
      std::vector<Real> primalsCopy;   
 
    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 internalEvaluatePrimal_EvalStatements(Args&&... args);
 
      template<typename... Args>
      static void internalEvaluateReverse_EvalStatements(Args&&... args);
 
      void internalResetPrimalValues(Position const& pos);
 
      void pushStmtData(Identifier const& index, Config::ArgumentSize const& numberOfPassiveArguments,
                        Real const& oldPrimalValue, EvalHandle evalHandle);
 
 
    public:
 
      PrimalValueBaseTape()
          : Base(),
            indexManager(Config::MaxArgumentSize),  // Reserve first items for passive values.
            statementData(Config::ChunkSize),
            // The following chunks must be large enough to store data for all arguments of one statement.
            rhsIdentifierData(std::max(Config::ChunkSize, Config::MaxArgumentSize)),
            passiveValueData(std::max(Config::ChunkSize, Config::MaxArgumentSize)),
            constantValueData(std::max(Config::ChunkSize, Config::MaxArgumentSize)),
            adjoints(1),  // Ensure that adjoint[0] exists, see its use in gradient() const.
            primals(0),
            primalsCopy(0) {
        checkPrimalSize(true);
 
        statementData.setNested(&indexManager.get());
        rhsIdentifierData.setNested(&statementData);
        passiveValueData.setNested(&rhsIdentifierData);
        constantValueData.setNested(&passiveValueData);
 
        Base::init(&constantValueData);
 
        Base::options.insert(TapeParameters::AdjointSize);
        Base::options.insert(TapeParameters::ConstantValuesSize);
        Base::options.insert(TapeParameters::LargestIdentifier);
        Base::options.insert(TapeParameters::PassiveValuesSize);
        Base::options.insert(TapeParameters::RhsIdentifiersSize);
        Base::options.insert(TapeParameters::PrimalSize);
        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 CountActiveArguments : public ForEachLeafLogic<CountActiveArguments> {
        public:
 
          template<typename Node>
          CODI_INLINE void handleActive(Node const& node, size_t& numberOfActiveArguments) {
            if (CODI_ENABLE_CHECK(Config::CheckZeroIndex, IndexManager::InactiveIndex != node.getIdentifier())) {
              numberOfActiveArguments += 1;
            }
          }
      };
 
      struct PushIdentifierPassiveAndConstant : public ForEachLeafLogic<PushIdentifierPassiveAndConstant> {
        public:
 
          template<typename Node>
          CODI_INLINE void handleActive(Node const& node, RhsIdentifierData& rhsIdentifierData,
                                        PassiveValueData& passiveValueData, ConstantValueData& constantValueData,
                                        size_t& curPassiveArgument) {
            CODI_UNUSED(constantValueData);
 
            Identifier rhsIndex = node.getIdentifier();
            if (CODI_ENABLE_CHECK(Config::CheckZeroIndex, IndexManager::InactiveIndex == rhsIndex)) {
              rhsIndex = curPassiveArgument;
 
              curPassiveArgument += 1;
              passiveValueData.pushData(node.getValue());
            }
 
            rhsIdentifierData.pushData(rhsIndex);
          }
 
          template<typename Node>
          CODI_INLINE void handleConstant(Node const& node, RhsIdentifierData& rhsIdentifierData,
                                          PassiveValueData& passiveValueData, ConstantValueData& constantValueData,
                                          size_t& curPassiveArgument) {
            CODI_UNUSED(rhsIdentifierData, passiveValueData, curPassiveArgument);
 
            using ConversionOperator = typename Node::template ConversionOperator<PassiveReal>;
 
            constantValueData.pushData(ConversionOperator::toDataStore(node.getValue()));
          }
      };
 
      struct JacobianExtractionLogic : public JacobianComputationLogic<JacobianExtractionLogic> {
        private:
          size_t pos;
 
        public:
 
          JacobianExtractionLogic() : pos(0) {}
 
          template<typename Node, typename Jacobian>
          CODI_INLINE void handleJacobianOnActive(Node const& node, Jacobian jacobianExpr, Identifier* rhsIdentifiers,
                                                  Real* jacobians) {
            rhsIdentifiers[pos] = node.getIdentifier();
            jacobians[pos] = ComputationTraits::adjointConversion<Real>(jacobianExpr);
            pos++;
          }
      };
 
    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())) {
          CountActiveArguments countActiveArguments;
          PushIdentifierPassiveAndConstant pushAll;
 
          size_t constexpr MaxActiveArgs = ExpressionTraits::NumberOfActiveTypeArguments<Rhs>::value;
          size_t constexpr MaxConstantArgs = ExpressionTraits::NumberOfConstantTypeArguments<Rhs>::value;
 
          codiAssert(MaxActiveArgs < Config::MaxArgumentSize);
          codiAssert(MaxConstantArgs < Config::MaxArgumentSize);
 
          size_t activeArguments = 0;
          countActiveArguments.eval(rhs.cast(), activeArguments);
 
          if (CODI_ENABLE_CHECK(Config::CheckEmptyStatements, 0 != activeArguments)) {
            statementData.reserveItems(1);
            rhsIdentifierData.reserveItems(MaxActiveArgs);
            passiveValueData.reserveItems(MaxActiveArgs - activeArguments);
            constantValueData.reserveItems(MaxConstantArgs);
 
            size_t passiveArguments = 0;
            pushAll.eval(rhs.cast(), rhsIdentifierData, passiveValueData, constantValueData, passiveArguments);
 
            bool generatedNewIndex = indexManager.get().template assignIndex<Impl>(lhs.cast().getIdentifier());
            checkPrimalSize(generatedNewIndex);
 
            Real& primalEntry = primals[lhs.cast().getIdentifier()];
            cast().pushStmtData(lhs.cast().getIdentifier(), passiveArguments, primalEntry,
                                StatementEvaluator::template createHandle<Impl, Impl, Rhs>());
 
            primalEntry = rhs.cast().getValue();
 
            if (Config::StatementEvents) {
              JacobianExtractionLogic getRhsIdentifiersAndJacobians;
              std::array<Identifier, MaxActiveArgs> rhsIdentifiers;
              std::array<Real, MaxActiveArgs> jacobians;
              getRhsIdentifiersAndJacobians.eval(rhs.cast(), Real(1.0), rhsIdentifiers.data(), jacobians.data());
 
              EventSystem<Impl>::notifyStatementStoreOnTapeListeners(cast(), lhs.cast().getIdentifier(),
                                                                     rhs.cast().getValue(), MaxActiveArgs,
                                                                     rhsIdentifiers.data(), jacobians.data());
            }
          } 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 Real internalRegisterInput(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& value,
                                             bool unusedIndex) {
        if (TapeTypes::IsLinearIndexHandler) {
          statementData.reserveItems(1);
        }
 
        bool generatedNewIndex;
        if (unusedIndex) {
          generatedNewIndex = indexManager.get().template assignUnusedIndex<Impl>(value.cast().getIdentifier());
        } else {
          generatedNewIndex = indexManager.get().template assignIndex<Impl>(value.cast().getIdentifier());
        }
        checkPrimalSize(generatedNewIndex);
 
        Real& primalEntry = primals[value.cast().getIdentifier()];
        if (TapeTypes::IsLinearIndexHandler) {
          cast().pushStmtData(value.cast().getIdentifier(), Config::StatementInputTag, primalEntry,
                              StatementEvaluator::template createHandle<Impl, Impl, Lhs>());
        }
 
        Real oldValue = primalEntry;
        primalEntry = value.cast().value();
 
        return oldValue;
      }
 
    public:
 
 
      template<typename Lhs>
      CODI_INLINE void registerInput(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& value) {
        internalRegisterInput(value, true);
        EventSystem<Impl>::notifyTapeRegisterInputListeners(cast(), value.cast().value(), value.cast().getIdentifier());
      }
 
      CODI_INLINE void clearAdjoints(AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        CODI_UNUSED(adjointsManagement);
 
        size_t maxSize = std::min((size_t)indexManager.get().getLargestCreatedIndex() + 1, adjoints.size());
        for (size_t i = 0; i < maxSize; i += 1) {
          adjoints[i] = Gradient();
        }
      }
 
      CODI_INLINE void reset(bool resetAdjoints = true,
                             AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        for (Real& primal : primals) {
          primal = Real();
        }
 
        Base::reset(resetAdjoints, adjointsManagement);
      }
 
 
    protected:
 
      CODI_INLINE TapeValues internalGetTapeValues() const {
        std::string name;
        if (TapeTypes::IsLinearIndexHandler) {
          name = "CoDi Tape Statistics ( PrimalValueLinearTape )";
        } else {
          name = "CoDi Tape Statistics ( PrimalValueReuseTape )";
        }
        TapeValues values = TapeValues(name);
 
        size_t nAdjoints = indexManager.get().getLargestCreatedIndex();
        double memoryAdjoints = static_cast<double>(nAdjoints) * static_cast<double>(sizeof(Gradient));
 
        size_t nPrimals = indexManager.get().getLargestCreatedIndex();
        double memoryPrimals = static_cast<double>(nPrimals) * static_cast<double>(sizeof(Real));
 
        values.addSection("Adjoint vector");
        values.addUnsignedLongEntry("Number of adjoints", nAdjoints);
        values.addDoubleEntry("Memory allocated", memoryAdjoints, TapeValues::LocalReductionOperation::Sum, true, true);
 
        values.addSection("Primal vector");
        values.addUnsignedLongEntry("Number of primals", nPrimals);
        values.addDoubleEntry("Memory allocated", memoryPrimals, TapeValues::LocalReductionOperation::Sum, true, true);
 
        values.addSection("Index manager");
        indexManager.get().addToTapeValues(values);
 
        values.addSection("Statement entries");
        statementData.addToTapeValues(values);
        values.addSection("Rhs identifiers entries");
        rhsIdentifierData.addToTapeValues(values);
        values.addSection("Passive value entries");
        passiveValueData.addToTapeValues(values);
        values.addSection("Constant value entries");
        constantValueData.addToTapeValues(values);
 
        return values;
      }
 
      /******************************************************************************
       * Protected helper function for CustomAdjointVectorEvaluationTapeInterface
       */
 
      template<typename AdjointVector>
      ADJOINT_VECTOR_TYPE* selectAdjointVector(VectorAccess<AdjointVector>* vectorAccess, AdjointVector data) {
        CODI_UNUSED(vectorAccess, data);
 
#if CODI_VariableAdjointInterfaceInPrimalTapes
        return vectorAccess;
#else
        CODI_STATIC_ASSERT(CODI_T(std::is_same<typename std::remove_reference<AdjointVector>::type, Gradient*>::value),
                           "Please enable 'CODI_VariableAdjointInterfaceInPrimalTapes' in order"
                           " to use custom adjoint vectors in the primal value tapes.");
 
        return data;
#endif
      }
 
      struct IncrementForwardLogic : public JacobianComputationLogic<IncrementForwardLogic> {
        public:
 
          template<typename Node>
          CODI_INLINE void handleJacobianOnActive(Node const& node, Real jacobian, Gradient& lhsTangent,
                                                  ADJOINT_VECTOR_TYPE* adjointVector) {
            CODI_UNUSED(lhsTangent);
 
            if (CODI_ENABLE_CHECK(Config::IgnoreInvalidJacobians, RealTraits::isTotalFinite(jacobian))) CODI_Likely {
#if CODI_VariableAdjointInterfaceInPrimalTapes
              adjointVector->updateTangentWithLhs(node.getIdentifier(), jacobian);
#else
              lhsTangent += jacobian * adjointVector[node.getIdentifier()];
#endif
            }
          }
      };
 
      CODI_WRAP_FUNCTION(Wrap_internalEvaluateForward_EvalStatements, Impl::internalEvaluateForward_EvalStatements);
 
      template<bool copyPrimal, typename AdjointVector>
      CODI_NO_INLINE void internalEvaluateForward(Position const& start, Position const& end, AdjointVector&& data) {
        CODI_STATIC_ASSERT(
            Config::VariableAdjointInterfaceInPrimalTapes ||
                CODI_T(std::is_same<typename std::remove_reference<AdjointVector>::type, Gradient*>::value),
            "Please enable 'CODI_VariableAdjointInterfaceInPrimalTapes' in order"
            " to use custom adjoint vectors in the primal value tapes.");
 
        std::vector<Real> primalsCopy(0);
        Real* primalData = primals.data();
 
        if (copyPrimal) {
          primalsCopy = primals;
          primalData = primalsCopy.data();
        }
 
        VectorAccess<AdjointVector> vectorAccess(data, primalData);
 
        ADJOINT_VECTOR_TYPE* dataVector = selectAdjointVector<AdjointVector>(&vectorAccess, data);
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(
            cast(), start, end, &vectorAccess, EventHints::EvaluationKind::Forward, EventHints::Endpoint::Begin);
 
        Wrap_internalEvaluateForward_EvalStatements evalFunc{};
        Base::llfByteData.evaluateForward(start, end, evalFunc, cast(), primalData, dataVector);
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(cast(), start, end, &vectorAccess,
                                                       EventHints::EvaluationKind::Forward, EventHints::Endpoint::End);
      }
 
      struct IncrementReversalLogic : public JacobianComputationLogic<IncrementReversalLogic> {
        public:
 
          template<typename Node>
          CODI_INLINE void handleJacobianOnActive(Node const& node, Real jacobian, Gradient const& lhsAdjoint,
                                                  ADJOINT_VECTOR_TYPE* adjointVector) {
            CODI_UNUSED(lhsAdjoint);
 
            if (CODI_ENABLE_CHECK(Config::IgnoreInvalidJacobians, RealTraits::isTotalFinite(jacobian))) CODI_Likely {
#if CODI_VariableAdjointInterfaceInPrimalTapes
              adjointVector->updateAdjointWithLhs(node.getIdentifier(), jacobian);
#else
              adjointVector[node.getIdentifier()] += jacobian * lhsAdjoint;
#endif
            }
          }
      };
 
      CODI_WRAP_FUNCTION(Wrap_internalEvaluateReverse_EvalStatements, Impl::internalEvaluateReverse_EvalStatements);
 
      template<bool copyPrimal, typename AdjointVector>
      CODI_INLINE void internalEvaluateReverse(Position const& start, Position const& end, AdjointVector&& data) {
        CODI_STATIC_ASSERT(
            Config::VariableAdjointInterfaceInPrimalTapes ||
                CODI_T(std::is_same<typename std::remove_reference<AdjointVector>::type, Gradient*>::value),
            "Please enable 'CODI_VariableAdjointInterfaceInPrimalTapes' in order"
            " to use custom adjoint vectors in the primal value tapes.");
 
        Real* primalData = primals.data();
 
        if (copyPrimal) {
          primalsCopy = primals;
          primalData = primalsCopy.data();
        }
 
        VectorAccess<AdjointVector> vectorAccess(data, primalData);
 
        ADJOINT_VECTOR_TYPE* dataVector = selectAdjointVector<AdjointVector>(&vectorAccess, data);
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(
            cast(), start, end, &vectorAccess, EventHints::EvaluationKind::Reverse, EventHints::Endpoint::Begin);
 
        Wrap_internalEvaluateReverse_EvalStatements evalFunc;
        Base::llfByteData.evaluateReverse(start, end, evalFunc, cast(), primalData, dataVector);
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(cast(), start, end, &vectorAccess,
                                                       EventHints::EvaluationKind::Reverse, EventHints::Endpoint::End);
      }
 
    public:
 
 
      using Base::evaluate;
 
      template<typename AdjointVector>
      CODI_INLINE void evaluate(Position const& start, Position const& end, AdjointVector&& data) {
        internalEvaluateReverse<!TapeTypes::IsLinearIndexHandler>(start, end, std::forward<AdjointVector>(data));
      }
 
      template<typename AdjointVector>
      CODI_INLINE void evaluateForward(Position const& start, Position const& end, AdjointVector&& data) {
        internalEvaluateForward<!TapeTypes::IsLinearIndexHandler>(start, end, std::forward<AdjointVector>(data));
      }
 
      CODI_INLINE Gradient* 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.
 
        std::swap(adjoints, other.adjoints);
        std::swap(primals, other.primals);
 
        Base::swap(other);
 
        // Ensure that the primals vector of both tapes are sized according to the index manager.
        checkPrimalSize(true);
        other.checkPrimalSize(true);
      }
 
      void deleteAdjointVector() {
        adjoints.resize(1);
      }
 
      void resizeAdjointVector() {
        checkAdjointSize(indexManager.get().getLargestCreatedIndex());
      }
 
      void beginUseAdjointVector() {}
 
      void endUseAdjointVector() {}
 
      size_t getParameter(TapeParameters parameter) const {
        switch (parameter) {
          case TapeParameters::AdjointSize:
            return adjoints.size();
            break;
          case TapeParameters::ConstantValuesSize:
            return constantValueData.getDataSize();
            break;
          case TapeParameters::LargestIdentifier:
            return indexManager.get().getLargestCreatedIndex();
            break;
          case TapeParameters::PassiveValuesSize:
            return passiveValueData.getDataSize();
            break;
          case TapeParameters::RhsIdentifiersSize:
            return rhsIdentifierData.getDataSize();
            break;
          case TapeParameters::PrimalSize:
            return primals.size();
            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::ConstantValuesSize:
            constantValueData.resize(value);
            break;
          case TapeParameters::LargestIdentifier:
            CODI_EXCEPTION("Tried to set a get only option.");
            break;
          case TapeParameters::PassiveValuesSize:
            passiveValueData.resize(value);
            break;
          case TapeParameters::RhsIdentifiersSize:
            rhsIdentifierData.resize(value);
            break;
          case TapeParameters::PrimalSize:
            primals.resize(value);
            break;
          case TapeParameters::StatementSize:
            return statementData.resize(value);
          default:
            Base::setParameter(parameter, value);
            break;
        }
      }
 
      VectorAccess<Gradient*>* createVectorAccess() {
        return createVectorAccessCustomAdjoints(adjoints.data());
      }
 
      template<typename AdjointVector>
      VectorAccess<AdjointVector>* createVectorAccessCustomAdjoints(AdjointVector&& data) {
        return new VectorAccess<AdjointVector>(data, primals.data());
      }
 
      template<typename Adjoint>
      VectorAccess<Adjoint*>* createVectorAccessCustomAdjoints(Adjoint* data) {
        return new VectorAccess<Adjoint*>(data, primals.data());
      }
 
      void deleteVectorAccess(VectorAccessInterface<Real, Identifier>* access) {
        delete access;
      }
 
      /*******************************************************************************/
 
      template<typename Lhs>
      Real registerExternalFunctionOutput(LhsExpressionInterface<Real, Gradient, Impl, Lhs>& value) {
        return internalRegisterInput(value, false);
      }
 
      /*******************************************************************************/
 
      using Base::evaluateForward;
 
      void evaluateForward(Position const& start, Position const& end,
                           AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(indexManager.get().getLargestCreatedIndex());
        }
 
        codiAssert(indexManager.get().getLargestCreatedIndex() < (Identifier)adjoints.size());
 
        cast().evaluateForward(start, end, adjoints.data());
      }
 
      /******************************************************************************
       * Public helper function for ManualStatementPushTapeInterface
       */
 
    public:
 
      template<size_t T_size>
      struct JacobianStatementGenerator {
        public:
 
          static size_t constexpr size = T_size;  
 
          /*******************************************************************************/
 
          template<typename Expr, typename... Args>
          static Real statementEvaluateForward(Args&&... args) {
            CODI_UNUSED(args...);
            CODI_EXCEPTION("Forward evaluation of jacobian statement not possible.");
          }
 
          template<typename Expr, typename... Args>
          static Real statementEvaluatePrimal(Args&&... args) {
            CODI_UNUSED(args...);
            CODI_EXCEPTION("Primal evaluation of jacobian statement not possible.");
          }
 
          template<typename Expr>
          static void statementEvaluateReverse(Real* primalVector, ADJOINT_VECTOR_TYPE* adjointVector,
                                               Gradient lhsAdjoint, Config::ArgumentSize numberOfPassiveArguments,
                                               size_t& curConstantPos, PassiveReal const* const constantValues,
                                               size_t& curPassivePos, Real const* const passiveValues,
                                               size_t& curRhsIdentifiersPos, Identifier const* const rhsIdentifiers) {
            CODI_UNUSED(primalVector, curConstantPos, constantValues);
 
            size_t endPos = curRhsIdentifiersPos - numberOfPassiveArguments;
 
            bool const lhsZero = evalJacobianReverse(adjointVector, lhsAdjoint, curPassivePos, passiveValues,
                                                     curRhsIdentifiersPos, rhsIdentifiers, endPos);
 
            if (Config::SkipZeroAdjointEvaluation && lhsZero) {
              curPassivePos -= numberOfPassiveArguments;
              curRhsIdentifiersPos -= numberOfPassiveArguments;
            }
          }
 
          template<typename Expr, typename... Args>
          static WriteInfo statementGetWriteInformation(Args&&... args) {
            CODI_UNUSED(args...);
 
            WriteInfo writeInfo;
            writeInfo.numberOfActiveArguments = ExpressionTraits::NumberOfActiveTypeArguments<Expr>::value;
            writeInfo.numberOfConstantArguments = ExpressionTraits::NumberOfConstantTypeArguments<Expr>::value;
            writeInfo.stmtExpression = "Impl, typename Impl::template JacobianStatementGenerator<" +
                                       std::to_string(size) + ">, codi::JacobianExpression<" + std::to_string(size) +
                                       ">";
            writeInfo.mathRepresentation = "Jacobian statement";
            return writeInfo;
          }
 
          /*******************************************************************************/
 
          template<typename Expr, typename... Args>
          static Real statementEvaluateForwardInner(Args&&... args) {
            CODI_UNUSED(args...);
            CODI_EXCEPTION("Forward evaluation of jacobian statement not possible.");
 
            return Real();
          }
 
          template<typename Expr, typename... Args>
          static Real statementEvaluatePrimalInner(Args&&... args) {
            CODI_UNUSED(args...);
            CODI_EXCEPTION("Primal evaluation of jacobian statement not possible.");
 
            return Real();
          }
 
          template<typename Expr>
          static void statementEvaluateReverseInner(Real* primalVector, ADJOINT_VECTOR_TYPE* adjointVector,
                                                    Gradient lhsAdjoint, size_t& curConstantPos,
                                                    PassiveReal const* const constantValues,
                                                    size_t& curRhsIdentifiersPos,
                                                    Identifier const* const rhsIdentifiers) {
            CODI_UNUSED(primalVector, curConstantPos, constantValues);
 
            size_t passivePos = size;
            size_t rhsPos = curRhsIdentifiersPos + size;
            size_t endPos = curRhsIdentifiersPos;
 
            evalJacobianReverse(adjointVector, lhsAdjoint, passivePos, primalVector, rhsPos, rhsIdentifiers, endPos);
          }
 
 
        private:
 
          static bool evalJacobianReverse(ADJOINT_VECTOR_TYPE* adjointVector, Gradient lhsAdjoint,
                                          size_t& curPassivePos, Real const* const passiveValues,
                                          size_t& curRhsIdentifiersPos, Identifier const* const rhsIdentifiers,
                                          size_t endRhsIdentifiersPos) {
#if CODI_VariableAdjointInterfaceInPrimalTapes
            CODI_UNUSED(lhsAdjoint);
            bool const lhsZero = adjointVector->isLhsZero();
#else
            bool const lhsZero = RealTraits::isTotalZero(lhsAdjoint);
#endif
 
            if (CODI_ENABLE_CHECK(Config::SkipZeroAdjointEvaluation, !lhsZero)) {
              while (curRhsIdentifiersPos > endRhsIdentifiersPos) {
                curPassivePos -= 1;
                curRhsIdentifiersPos -= 1;
 
                Real const& jacobian = passiveValues[curPassivePos];
#if CODI_VariableAdjointInterfaceInPrimalTapes
                adjointVector->updateAdjointWithLhs(rhsIdentifiers[curRhsIdentifiersPos], jacobian);
#else
                adjointVector[rhsIdentifiers[curRhsIdentifiersPos]] += jacobian * lhsAdjoint;
#endif
              }
            }
 
            return lhsZero;
          }
      };
 
      /*******************************************************************************/
 
      void pushJacobianManual(Real const& jacobian, Real const& value, Identifier const& index) {
        CODI_UNUSED(value);
 
        cast().incrementManualPushCounter();
 
        passiveValueData.pushData(jacobian);
        rhsIdentifierData.pushData(index);
 
        if (Config::StatementEvents) {
          if (this->manualPushCounter == this->manualPushGoal) {
            // emit statement event
            Real* jacobians;
            Identifier* rhsIdentifiers;
            passiveValueData.getDataPointers(jacobians);
            rhsIdentifierData.getDataPointers(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);
 
        codiAssert(size < Config::MaxArgumentSize);
 
        statementData.reserveItems(1);
        rhsIdentifierData.reserveItems(size);
        passiveValueData.reserveItems(size);
 
        indexManager.get().template assignIndex<Impl>(lhsIndex);
        Real& primalEntry = primals[lhsIndex];
        cast().pushStmtData(lhsIndex, size, primalEntry, PrimalValueBaseTape::jacobianExpressions[size]);
 
        primalEntry = lhsValue;
 
        cast().initializeManualPushData(lhsValue, lhsIndex, size);
      }
 
      /*******************************************************************************/
 
      void createStatementManual(Identifier const& lhsIndex, Real const& lhsValue,
                                 Config::ArgumentSize const& nActiveValues, Identifier const* const rhsIdentifiers,
                                 Config::ArgumentSize const& nPassiveValues, Real const* const rhsPrimals,
                                 Config::ArgumentSize const& nConstants, Real const* const rhsConstant,
                                 EvalHandle const& evalHandle) {
        Impl& impl = cast();
        if (Config::StatementLowLevelFunctionTag == nPassiveValues) CODI_Unlikely {
          // TODO.
        } else if (Config::StatementInputTag == nPassiveValues && TapeTypes::IsLinearIndexHandler) CODI_Unlikely {
          statementData.reserveItems(1);
          if (TapeTypes::IsLinearIndexHandler) {
            primals[lhsIndex] = lhsValue;
          }
          impl.pushStmtData(lhsIndex, Config::StatementInputTag, lhsValue, evalHandle);
 
        } else CODI_Likely {
          statementData.reserveItems(1);
          rhsIdentifierData.reserveItems(nActiveValues);
          passiveValueData.reserveItems(nPassiveValues);
          constantValueData.reserveItems(nConstants);
 
          if (TapeTypes::IsLinearIndexHandler) {
            primals[lhsIndex] = lhsValue;
          }
 
          impl.pushStmtData(lhsIndex, nPassiveValues, lhsValue, evalHandle);
          impl.initializeManualPushData(lhsValue, lhsIndex, nPassiveValues);
 
          for (size_t activeCount = 0; activeCount < nActiveValues; activeCount++) {
            rhsIdentifierData.pushData(rhsIdentifiers[activeCount]);
          }
          for (size_t passiveCount = 0; passiveCount < nPassiveValues; passiveCount++) {
            cast().incrementManualPushCounter();
            passiveValueData.pushData(rhsPrimals[passiveCount]);
          }
          for (size_t constantCount = 0; constantCount < nConstants; constantCount++) {
            constantValueData.pushData(rhsConstant[constantCount]);
          }
        }
      }
 
      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>,  primals.data(),
                                          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, Real(), EvalHandle());
      }
 
      /*******************************************************************************/
 
      CODI_INLINE void evaluate(Position const& start, Position const& end,
                                AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(indexManager.get().getLargestCreatedIndex());
        }
 
        codiAssert(indexManager.get().getLargestCreatedIndex() < (Identifier)adjoints.size());
 
        evaluate(start, end, adjoints.data());
      }
 
      CODI_INLINE void resetTo(Position const& pos, bool resetAdjoints = true,
                               AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        cast().internalResetPrimalValues(pos);
 
        Base::resetTo(pos, resetAdjoints, adjointsManagement);
      }
 
      /*******************************************************************************/
 
      void evaluateKeepState(Position const& start, Position const& end,
                             AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(indexManager.get().getLargestCreatedIndex());
        }
 
        codiAssert(indexManager.get().getLargestCreatedIndex() < (Identifier)adjoints.size());
 
        evaluateKeepState(start, end, adjoints.data());
      }
 
      template<typename AdjointVector>
      void evaluateKeepState(Position const& start, Position const& end, AdjointVector&& data) {
        internalEvaluateReverse<false>(start, end, std::forward<AdjointVector>(data));
 
        if (!TapeTypes::IsLinearIndexHandler) {
          evaluatePrimal(end, start);
        }
      }
 
      void evaluateForwardKeepState(Position const& start, Position const& end,
                                    AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        if (AdjointsManagement::Automatic == adjointsManagement) {
          checkAdjointSize(indexManager.get().getLargestCreatedIndex());
        }
 
        codiAssert(indexManager.get().getLargestCreatedIndex() < (Identifier)adjoints.size());
 
        evaluateForwardKeepState(start, end, adjoints.data());
      }
 
      template<typename AdjointVector>
      void evaluateForwardKeepState(Position const& start, Position const& end, AdjointVector&& data) {
        if (!TapeTypes::IsLinearIndexHandler) {
          cast().internalResetPrimalValues(end);
        }
 
        internalEvaluateForward<false>(start, end, std::forward<AdjointVector>(data));
      }
 
    protected:
 
      /******************************************************************************
       * Protected helper function for PrimalEvaluationTapeInterface
       */
 
      CODI_WRAP_FUNCTION(Wrap_internalEvaluatePrimal_EvalStatements, Impl::internalEvaluatePrimal_EvalStatements);
 
    public:
 
      /*******************************************************************************/
 
      using Base::evaluatePrimal;
 
      CODI_NO_INLINE void evaluatePrimal(Position const& start, Position const& end) {
        // TODO: implement primal value only accessor
        PrimalAdjointVectorAccess<Real, Identifier, Gradient*> primalAdjointAccess(adjoints.data(), primals.data());
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(cast(), start, end, &primalAdjointAccess,
                                                       EventHints::EvaluationKind::Primal, EventHints::Endpoint::Begin);
 
        Wrap_internalEvaluatePrimal_EvalStatements evalFunc{};
        Base::llfByteData.evaluateForward(start, end, evalFunc, cast(), primals.data());
 
        EventSystem<Impl>::notifyTapeEvaluateListeners(cast(), start, end, &primalAdjointAccess,
                                                       EventHints::EvaluationKind::Primal, EventHints::Endpoint::End);
      }
 
      Real& primal(Identifier const& identifier) {
        return primals[identifier];
      }
 
      Real const& primal(Identifier const& identifier) const {
        return primals[identifier];
      }
 
      /*******************************************************************************/
 
      template<typename Expr>
      static WriteInfo statementGetWriteInformation(Real* primalVector, Config::ArgumentSize numberOfPassiveArguments,
                                                    size_t& curConstantPos, PassiveReal const* const constantValues,
                                                    size_t& curPassivePos, Real const* const passiveValues,
                                                    size_t& curRhsIdentifiersPos,
                                                    Identifier const* const rhsIdentifiers) {
        if (Config::StatementInputTag != numberOfPassiveArguments) {
          for (Config::ArgumentSize curPos = 0; curPos < numberOfPassiveArguments; curPos += 1) {
            primalVector[curPos] = passiveValues[curPassivePos + curPos];
          }
        }
 
        using Constructor = ConstructStaticContextLogic<Expr, Impl, 0, 0>;
        using StaticRhs = typename Constructor::ResultType;
 
        StaticRhs staticRhs = Constructor::construct(primalVector, &rhsIdentifiers[curRhsIdentifiersPos],
                                                     &constantValues[curConstantPos]);
 
        WriteInfo writeInfo;
        writeInfo.numberOfActiveArguments = ExpressionTraits::NumberOfActiveTypeArguments<Expr>::value;
        writeInfo.numberOfConstantArguments = ExpressionTraits::NumberOfConstantTypeArguments<Expr>::value;
 
        // The Impl template name is added to simplify the writer and reader interfaces EvalHandles. It represents the
        // first and second template args in the createHandle method. For the manual push, the second Impl template is
        // instead replaced with JacobianGenerator<size>.
 
        writeInfo.stmtExpression = "Impl, Impl, ";
        writeInfo.stmtExpression += demangleName<Expr>();
        MathStatementGenLogic<Identifier> mathGen(Config::MaxArgumentSize);
        mathGen.eval(staticRhs, writeInfo.mathRepresentation);
 
        return writeInfo;
      }
 
      template<typename Rhs>
      static Real statementEvaluateForwardInner(Real* primalVector, ADJOINT_VECTOR_TYPE* adjointVector,
                                                Gradient& lhsTangent, size_t& curConstantPos,
                                                PassiveReal const* const constantValues, size_t& curRhsIdentifiersPos,
                                                Identifier const* const rhsIdentifiers) {
        using Constructor = ConstructStaticContextLogic<Rhs, Impl, 0, 0>;
        using StaticRhs = typename Constructor::ResultType;
 
        StaticRhs staticsRhs = Constructor::construct(primalVector, &rhsIdentifiers[curRhsIdentifiersPos],
                                                      &constantValues[curConstantPos]);
 
        IncrementForwardLogic incrementForward;
 
        incrementForward.eval(staticsRhs, Real(1.0), lhsTangent, adjointVector);
        return staticsRhs.getValue();
      }
 
      template<typename Func>
      static Real statementEvaluateForwardFull(Func const& evalInner, size_t const& maxActiveArgs,
                                               size_t const& maxConstantArgs, Real* primalVector,
                                               ADJOINT_VECTOR_TYPE* adjointVector, Gradient& lhsTangent,
                                               Config::ArgumentSize numberOfPassiveArguments, size_t& curConstantPos,
                                               PassiveReal const* const constantValues, size_t& curPassivePos,
                                               Real const* const passiveValues, size_t& curRhsIdentifiersPos,
                                               Identifier const* const rhsIdentifiers) {
        for (Config::ArgumentSize curPos = 0; curPos < numberOfPassiveArguments; curPos += 1) {
          primalVector[curPos] = passiveValues[curPassivePos + curPos];
        }
 
        Real ret = evalInner(primalVector, adjointVector, lhsTangent, curConstantPos, constantValues,
                             curRhsIdentifiersPos, rhsIdentifiers);
 
        // Adapt vector positions.
        curConstantPos += maxConstantArgs;
        curPassivePos += numberOfPassiveArguments;
        curRhsIdentifiersPos += maxActiveArgs;
 
        return ret;
      }
 
      template<typename Rhs>
      static Real statementEvaluatePrimalInner(Real* primalVector, size_t& curConstantPos,
                                               PassiveReal const* const constantValues, size_t& curRhsIdentifiersPos,
                                               Identifier const* const rhsIdentifiers) {
        using Constructor = ConstructStaticContextLogic<Rhs, Impl, 0, 0>;
        using StaticRhs = typename Constructor::ResultType;
 
        StaticRhs staticsRhs = Constructor::construct(primalVector, &rhsIdentifiers[curRhsIdentifiersPos],
                                                      &constantValues[curConstantPos]);
 
        return staticsRhs.getValue();
      }
 
      template<typename Func>
      static Real statementEvaluatePrimalFull(Func const& evalInner, size_t const& maxActiveArgs,
                                              size_t const& maxConstantArgs, Real* primalVector,
                                              Config::ArgumentSize numberOfPassiveArguments, size_t& curConstantPos,
                                              PassiveReal const* const constantValues, size_t& curPassivePos,
                                              Real const* const passiveValues, size_t& curRhsIdentifiersPos,
                                              Identifier const* const rhsIdentifiers) {
        for (Config::ArgumentSize curPos = 0; curPos < numberOfPassiveArguments; curPos += 1) {
          primalVector[curPos] = passiveValues[curPassivePos + curPos];
        }
 
        Real ret = evalInner(primalVector, curConstantPos, constantValues, curRhsIdentifiersPos, rhsIdentifiers);
 
        // Adapt vector positions.
        curConstantPos += maxConstantArgs;
        curPassivePos += numberOfPassiveArguments;
        curRhsIdentifiersPos += maxActiveArgs;
 
        return ret;
      }
 
      template<typename Rhs>
      CODI_INLINE static void statementEvaluateReverseInner(Real* primalVector, ADJOINT_VECTOR_TYPE* adjointVector,
                                                            Gradient lhsAdjoint, size_t& curConstantPos,
                                                            PassiveReal const* const constantValues,
                                                            size_t& curRhsIdentifiersPos,
                                                            Identifier const* const rhsIdentifiers) {
        using Constructor = ConstructStaticContextLogic<Rhs, Impl, 0, 0>;
        using StaticRhs = typename Constructor::ResultType;
 
        StaticRhs staticsRhs = Constructor::construct(primalVector, &rhsIdentifiers[curRhsIdentifiersPos],
                                                      &constantValues[curConstantPos]);
 
        IncrementReversalLogic incrementReverse;
 
        incrementReverse.eval(staticsRhs, Real(1.0), const_cast<Gradient const&>(lhsAdjoint), adjointVector);
      }
 
      template<typename Func>
      CODI_INLINE static void statementEvaluateReverseFull(
          Func const& evalInner, size_t const& maxActiveArgs, size_t const& maxConstantArgs, Real* primalVector,
          ADJOINT_VECTOR_TYPE* adjointVector, Gradient lhsAdjoint, Config::ArgumentSize numberOfPassiveArguments,
          size_t& curConstantPos, PassiveReal const* const constantValues, size_t& curPassivePos,
          Real const* const passiveValues, size_t& curRhsIdentifiersPos, Identifier const* const rhsIdentifiers) {
        // Adapt vector positions.
        curConstantPos -= maxConstantArgs;
        curPassivePos -= numberOfPassiveArguments;
        curRhsIdentifiersPos -= maxActiveArgs;
 
#if CODI_VariableAdjointInterfaceInPrimalTapes
        if (CODI_ENABLE_CHECK(Config::SkipZeroAdjointEvaluation, !adjointVector->isLhsZero())) CODI_Likely {
#else
        if (CODI_ENABLE_CHECK(Config::SkipZeroAdjointEvaluation, !RealTraits::isTotalZero(lhsAdjoint))) CODI_Likely {
#endif
          for (Config::ArgumentSize curPos = 0; curPos < numberOfPassiveArguments; curPos += 1) {
            primalVector[curPos] = passiveValues[curPassivePos + curPos];
          }
 
          evalInner(primalVector, adjointVector, lhsAdjoint, curConstantPos, constantValues, curRhsIdentifiersPos,
                    rhsIdentifiers);
        }
      }
 
      /*******************************************************************************/
 
      template<typename Rhs>
      static Real statementEvaluateForward(Real* primalVector, ADJOINT_VECTOR_TYPE* adjointVector, Gradient& lhsTangent,
                                           Config::ArgumentSize numberOfPassiveArguments, size_t& curConstantPos,
                                           PassiveReal const* const constantValues, size_t& curPassivePos,
                                           Real const* const passiveValues, size_t& curRhsIdentifiersPos,
                                           Identifier const* const rhsIdentifiers) {
        size_t constexpr MaxActiveArgs = ExpressionTraits::NumberOfActiveTypeArguments<Rhs>::value;
        size_t constexpr MaxConstantArgs = ExpressionTraits::NumberOfConstantTypeArguments<Rhs>::value;
 
        return statementEvaluateForwardFull(statementEvaluateForwardInner<Rhs>, MaxActiveArgs, MaxConstantArgs,
                                            primalVector, adjointVector, lhsTangent, numberOfPassiveArguments,
                                            curConstantPos, constantValues, curPassivePos, passiveValues,
                                            curRhsIdentifiersPos, rhsIdentifiers);
      }
 
      template<typename Rhs>
      static Real statementEvaluatePrimal(Real* primalVector, Config::ArgumentSize numberOfPassiveArguments,
                                          size_t& curConstantPos, PassiveReal const* const constantValues,
                                          size_t& curPassivePos, Real const* const passiveValues,
                                          size_t& curRhsIdentifiersPos, Identifier const* const rhsIdentifiers) {
        size_t constexpr MaxActiveArgs = ExpressionTraits::NumberOfActiveTypeArguments<Rhs>::value;
        size_t constexpr MaxConstantArgs = ExpressionTraits::NumberOfConstantTypeArguments<Rhs>::value;
 
        return statementEvaluatePrimalFull(statementEvaluatePrimalInner<Rhs>, MaxActiveArgs, MaxConstantArgs,
                                           primalVector, numberOfPassiveArguments, curConstantPos, constantValues,
                                           curPassivePos, passiveValues, curRhsIdentifiersPos, rhsIdentifiers);
      }
 
      template<typename Rhs>
      CODI_INLINE static void statementEvaluateReverse(Real* primalVector, ADJOINT_VECTOR_TYPE* adjointVector,
                                                       Gradient lhsAdjoint,
                                                       Config::ArgumentSize numberOfPassiveArguments,
                                                       size_t& curConstantPos, PassiveReal const* const constantValues,
                                                       size_t& curPassivePos, Real const* const passiveValues,
                                                       size_t& curRhsIdentifiersPos,
                                                       Identifier const* const rhsIdentifiers) {
        size_t constexpr maxActiveArgs = ExpressionTraits::NumberOfActiveTypeArguments<Rhs>::value;
        size_t constexpr maxConstantArgs = ExpressionTraits::NumberOfConstantTypeArguments<Rhs>::value;
        statementEvaluateReverseFull(statementEvaluateReverseInner<Rhs>, maxActiveArgs, maxConstantArgs, primalVector,
                                     adjointVector, lhsAdjoint, numberOfPassiveArguments, curConstantPos,
                                     constantValues, curPassivePos, passiveValues, curRhsIdentifiersPos,
                                     rhsIdentifiers);
      }
 
    private:
 
      CODI_INLINE void checkAdjointSize(Identifier const& identifier) {
        if (identifier >= (Identifier)adjoints.size()) {
          resizeAdjointsVector();
        }
      }
 
      CODI_INLINE void checkPrimalSize(bool generatedNewIndex) {
        if (generatedNewIndex && indexManager.get().getLargestCreatedIndex() >= (Identifier)primals.size()) {
          resizePrimalVector();
        }
      }
 
      CODI_NO_INLINE void resizeAdjointsVector() {
        // overallocate as next multiple of Config::ChunkSize
        adjoints.resize(getNextMultiple((size_t)indexManager.get().getLargestCreatedIndex() + 1, Config::ChunkSize));
      }
 
      CODI_NO_INLINE void resizePrimalVector() {
        // overallocate as next multiple of Config::ChunkSize
        primals.resize(getNextMultiple((size_t)indexManager.get().getLargestCreatedIndex() + 1, Config::ChunkSize));
      }
  };
 
  template<size_t size>
  struct JacobianExpression {};
 
  template<size_t size>
  struct ExpressionTraits::NumberOfActiveTypeArguments<JacobianExpression<size>> {
      static size_t constexpr value = size;  
  };
 
  template<size_t size>
  struct ExpressionTraits::NumberOfConstantTypeArguments<JacobianExpression<size>> {
      static size_t constexpr value = 0;  
  };
 
#define CREATE_EXPRESSION(size)                                                                                        \
  TapeTypes::StatementEvaluator::template createHandle<Impl, typename Impl::template JacobianStatementGenerator<size>, \
                                                       JacobianExpression<size>>()
 
  template<typename TapeTypes, typename Impl>
  const CODI_DD(typename TapeTypes::EvalHandle,
                CODI_ANY) PrimalValueBaseTape<TapeTypes, Impl>::jacobianExpressions[Config::MaxArgumentSize] = {
      CREATE_EXPRESSION(0),   CREATE_EXPRESSION(1),   CREATE_EXPRESSION(2),   CREATE_EXPRESSION(3),
      CREATE_EXPRESSION(4),   CREATE_EXPRESSION(5),   CREATE_EXPRESSION(6),   CREATE_EXPRESSION(7),
      CREATE_EXPRESSION(8),   CREATE_EXPRESSION(9),   CREATE_EXPRESSION(10),  CREATE_EXPRESSION(11),
      CREATE_EXPRESSION(12),  CREATE_EXPRESSION(13),  CREATE_EXPRESSION(14),  CREATE_EXPRESSION(15),
      CREATE_EXPRESSION(16),  CREATE_EXPRESSION(17),  CREATE_EXPRESSION(18),  CREATE_EXPRESSION(19),
      CREATE_EXPRESSION(20),  CREATE_EXPRESSION(21),  CREATE_EXPRESSION(22),  CREATE_EXPRESSION(23),
      CREATE_EXPRESSION(24),  CREATE_EXPRESSION(25),  CREATE_EXPRESSION(26),  CREATE_EXPRESSION(27),
      CREATE_EXPRESSION(28),  CREATE_EXPRESSION(29),  CREATE_EXPRESSION(30),  CREATE_EXPRESSION(31),
      CREATE_EXPRESSION(32),  CREATE_EXPRESSION(33),  CREATE_EXPRESSION(34),  CREATE_EXPRESSION(35),
      CREATE_EXPRESSION(36),  CREATE_EXPRESSION(37),  CREATE_EXPRESSION(38),  CREATE_EXPRESSION(39),
      CREATE_EXPRESSION(40),  CREATE_EXPRESSION(41),  CREATE_EXPRESSION(42),  CREATE_EXPRESSION(43),
      CREATE_EXPRESSION(44),  CREATE_EXPRESSION(45),  CREATE_EXPRESSION(46),  CREATE_EXPRESSION(47),
      CREATE_EXPRESSION(48),  CREATE_EXPRESSION(49),  CREATE_EXPRESSION(50),  CREATE_EXPRESSION(51),
      CREATE_EXPRESSION(52),  CREATE_EXPRESSION(53),  CREATE_EXPRESSION(54),  CREATE_EXPRESSION(55),
      CREATE_EXPRESSION(56),  CREATE_EXPRESSION(57),  CREATE_EXPRESSION(58),  CREATE_EXPRESSION(59),
      CREATE_EXPRESSION(60),  CREATE_EXPRESSION(61),  CREATE_EXPRESSION(62),  CREATE_EXPRESSION(63),
      CREATE_EXPRESSION(64),  CREATE_EXPRESSION(65),  CREATE_EXPRESSION(66),  CREATE_EXPRESSION(67),
      CREATE_EXPRESSION(68),  CREATE_EXPRESSION(69),  CREATE_EXPRESSION(70),  CREATE_EXPRESSION(71),
      CREATE_EXPRESSION(72),  CREATE_EXPRESSION(73),  CREATE_EXPRESSION(74),  CREATE_EXPRESSION(75),
      CREATE_EXPRESSION(76),  CREATE_EXPRESSION(77),  CREATE_EXPRESSION(78),  CREATE_EXPRESSION(79),
      CREATE_EXPRESSION(80),  CREATE_EXPRESSION(81),  CREATE_EXPRESSION(82),  CREATE_EXPRESSION(83),
      CREATE_EXPRESSION(84),  CREATE_EXPRESSION(85),  CREATE_EXPRESSION(86),  CREATE_EXPRESSION(87),
      CREATE_EXPRESSION(88),  CREATE_EXPRESSION(89),  CREATE_EXPRESSION(90),  CREATE_EXPRESSION(91),
      CREATE_EXPRESSION(92),  CREATE_EXPRESSION(93),  CREATE_EXPRESSION(94),  CREATE_EXPRESSION(95),
      CREATE_EXPRESSION(96),  CREATE_EXPRESSION(97),  CREATE_EXPRESSION(98),  CREATE_EXPRESSION(99),
      CREATE_EXPRESSION(100), CREATE_EXPRESSION(101), CREATE_EXPRESSION(102), CREATE_EXPRESSION(103),
      CREATE_EXPRESSION(104), CREATE_EXPRESSION(105), CREATE_EXPRESSION(106), CREATE_EXPRESSION(107),
      CREATE_EXPRESSION(108), CREATE_EXPRESSION(109), CREATE_EXPRESSION(110), CREATE_EXPRESSION(111),
      CREATE_EXPRESSION(112), CREATE_EXPRESSION(113), CREATE_EXPRESSION(114), CREATE_EXPRESSION(115),
      CREATE_EXPRESSION(116), CREATE_EXPRESSION(117), CREATE_EXPRESSION(118), CREATE_EXPRESSION(119),
      CREATE_EXPRESSION(120), CREATE_EXPRESSION(121), CREATE_EXPRESSION(122), CREATE_EXPRESSION(123),
      CREATE_EXPRESSION(124), CREATE_EXPRESSION(125), CREATE_EXPRESSION(126), CREATE_EXPRESSION(127),
      CREATE_EXPRESSION(128), CREATE_EXPRESSION(129), CREATE_EXPRESSION(130), CREATE_EXPRESSION(131),
      CREATE_EXPRESSION(132), CREATE_EXPRESSION(133), CREATE_EXPRESSION(134), CREATE_EXPRESSION(135),
      CREATE_EXPRESSION(136), CREATE_EXPRESSION(137), CREATE_EXPRESSION(138), CREATE_EXPRESSION(139),
      CREATE_EXPRESSION(140), CREATE_EXPRESSION(141), CREATE_EXPRESSION(142), CREATE_EXPRESSION(143),
      CREATE_EXPRESSION(144), CREATE_EXPRESSION(145), CREATE_EXPRESSION(146), CREATE_EXPRESSION(147),
      CREATE_EXPRESSION(148), CREATE_EXPRESSION(149), CREATE_EXPRESSION(150), CREATE_EXPRESSION(151),
      CREATE_EXPRESSION(152), CREATE_EXPRESSION(153), CREATE_EXPRESSION(154), CREATE_EXPRESSION(155),
      CREATE_EXPRESSION(156), CREATE_EXPRESSION(157), CREATE_EXPRESSION(158), CREATE_EXPRESSION(159),
      CREATE_EXPRESSION(160), CREATE_EXPRESSION(161), CREATE_EXPRESSION(162), CREATE_EXPRESSION(163),
      CREATE_EXPRESSION(164), CREATE_EXPRESSION(165), CREATE_EXPRESSION(166), CREATE_EXPRESSION(167),
      CREATE_EXPRESSION(168), CREATE_EXPRESSION(169), CREATE_EXPRESSION(170), CREATE_EXPRESSION(171),
      CREATE_EXPRESSION(172), CREATE_EXPRESSION(173), CREATE_EXPRESSION(174), CREATE_EXPRESSION(175),
      CREATE_EXPRESSION(176), CREATE_EXPRESSION(177), CREATE_EXPRESSION(178), CREATE_EXPRESSION(179),
      CREATE_EXPRESSION(180), CREATE_EXPRESSION(181), CREATE_EXPRESSION(182), CREATE_EXPRESSION(183),
      CREATE_EXPRESSION(184), CREATE_EXPRESSION(185), CREATE_EXPRESSION(186), CREATE_EXPRESSION(187),
      CREATE_EXPRESSION(188), CREATE_EXPRESSION(189), CREATE_EXPRESSION(190), CREATE_EXPRESSION(191),
      CREATE_EXPRESSION(192), CREATE_EXPRESSION(193), CREATE_EXPRESSION(194), CREATE_EXPRESSION(195),
      CREATE_EXPRESSION(196), CREATE_EXPRESSION(197), CREATE_EXPRESSION(198), CREATE_EXPRESSION(199),
      CREATE_EXPRESSION(200), CREATE_EXPRESSION(201), CREATE_EXPRESSION(202), CREATE_EXPRESSION(203),
      CREATE_EXPRESSION(204), CREATE_EXPRESSION(205), CREATE_EXPRESSION(206), CREATE_EXPRESSION(207),
      CREATE_EXPRESSION(208), CREATE_EXPRESSION(209), CREATE_EXPRESSION(210), CREATE_EXPRESSION(211),
      CREATE_EXPRESSION(212), CREATE_EXPRESSION(213), CREATE_EXPRESSION(214), CREATE_EXPRESSION(215),
      CREATE_EXPRESSION(216), CREATE_EXPRESSION(217), CREATE_EXPRESSION(218), CREATE_EXPRESSION(219),
      CREATE_EXPRESSION(220), CREATE_EXPRESSION(221), CREATE_EXPRESSION(222), CREATE_EXPRESSION(223),
      CREATE_EXPRESSION(224), CREATE_EXPRESSION(225), CREATE_EXPRESSION(226), CREATE_EXPRESSION(227),
      CREATE_EXPRESSION(228), CREATE_EXPRESSION(229), CREATE_EXPRESSION(230), CREATE_EXPRESSION(231),
      CREATE_EXPRESSION(232), CREATE_EXPRESSION(233), CREATE_EXPRESSION(234), CREATE_EXPRESSION(235),
      CREATE_EXPRESSION(236), CREATE_EXPRESSION(237), CREATE_EXPRESSION(238), CREATE_EXPRESSION(239),
      CREATE_EXPRESSION(240), CREATE_EXPRESSION(241), CREATE_EXPRESSION(242), CREATE_EXPRESSION(243),
      CREATE_EXPRESSION(244), CREATE_EXPRESSION(245), CREATE_EXPRESSION(246), CREATE_EXPRESSION(247),
      CREATE_EXPRESSION(248), CREATE_EXPRESSION(249), CREATE_EXPRESSION(250), CREATE_EXPRESSION(251),
      CREATE_EXPRESSION(252)};
 
#undef CREATE_EXPRESSION
}