forwardEvaluation.hpp

/*
 * CoDiPack, a Code Differentiation Package
 *
 * Copyright (C) 2015-2025 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 "../config.h"
#include "../expressions/aggregate/aggregatedActiveType.hpp"
#include "../expressions/aggregate/arrayAccessExpression.hpp"
#include "../expressions/lhsExpressionInterface.hpp"
#include "../expressions/logic/helpers/jacobianComputationLogic.hpp"
#include "../misc/macros.hpp"
#include "../traits/expressionTraits.hpp"
#include "../traits/realTraits.hpp"
#include "../traits/tapeTraits.hpp"
#include "interfaces/gradientAccessTapeInterface.hpp"
#include "interfaces/internalStatementRecordingTapeInterface.hpp"

namespace codi {

  template<typename T_Real, typename T_Gradient>
  struct ForwardEvaluation : public InternalStatementRecordingTapeInterface<T_Gradient>,
                             public GradientAccessTapeInterface<T_Gradient, T_Gradient> {
    public:
 
      using Real = CODI_DD(T_Real, double);          
      using Gradient = CODI_DD(T_Gradient, double);  
 
      using PassiveReal = RealTraits::PassiveReal<Real>;  
      using Identifier = Gradient;  
 
      /*******************************************************************************/
 
      static bool constexpr AllowJacobianOptimization = true;  

      template<typename Real>
      CODI_INLINE void initIdentifier(Real& value, Identifier& identifier) {
        CODI_UNUSED(value);
        identifier = Identifier();
      }

      template<typename Real>
      CODI_INLINE void destroyIdentifier(Real& value, Identifier& identifier) {
        CODI_UNUSED(value, identifier);
      }

 
    private:
 
      struct LocalReverseLogic : public JacobianComputationLogic<LocalReverseLogic> {
        public:
          template<typename Node>
          CODI_INLINE void handleJacobianOnActive(Node const& node, Real jacobian, Gradient& lhsGradient) {
            if (CODI_ENABLE_CHECK(Config::IgnoreInvalidJacobians, RealTraits::isTotalFinite(jacobian))) {
              lhsGradient += node.gradient() * jacobian;
            }
          }
      };
 
    public:


      template<typename Lhs, typename Rhs>
      CODI_INLINE void store(LhsExpressionInterface<Real, Gradient, ForwardEvaluation, Lhs>& lhs,
                             ExpressionInterface<Real, Rhs> const& rhs) {
        LocalReverseLogic reversal;
 
        Gradient newGradient = Gradient();
        reversal.eval(rhs.cast(), 1.0, newGradient);
 
        lhs.cast().value() = rhs.cast().getValue();
        lhs.cast().gradient() = newGradient;
      }

      template<typename Aggregated, typename Type, typename Lhs, typename Rhs>
      CODI_INLINE void store(AggregatedActiveType<Aggregated, Type, Lhs>& lhs,
                             ExpressionInterface<Aggregated, Rhs> const& rhs) {
        using AggregatedTraits = RealTraits::AggregatedTypeTraits<Aggregated>;
        int constexpr Elements = AggregatedTraits::Elements;
        LocalReverseLogic reversal;
 
        Gradient newGradient[Elements] = {};
        static_for<Elements>([&](auto i) CODI_LAMBDA_INLINE {
          reversal.eval(ArrayAccessExpression<Aggregated, i.value, Rhs>(rhs), Real(1.0), newGradient[i.value]);
        });
 
        Aggregated newValue = rhs.cast().getValue();
        static_for<Elements>([&](auto i) CODI_LAMBDA_INLINE {
          lhs.values[i.value].value() = AggregatedTraits::template arrayAccess<i.value>(newValue);
          lhs.values[i.value].gradient() = newGradient[i.value];
        });
      }

      template<typename Lhs, typename Rhs>
      CODI_INLINE void store(LhsExpressionInterface<Real, Gradient, ForwardEvaluation, Lhs>& lhs,
                             LhsExpressionInterface<Real, Gradient, ForwardEvaluation, Rhs> const& rhs) {
        lhs.cast().value() = rhs.cast().getValue();
        lhs.cast().gradient() = rhs.cast().getGradient();
      }

      template<typename Lhs>
      CODI_INLINE void store(LhsExpressionInterface<Real, Gradient, ForwardEvaluation, Lhs>& lhs, Real const& rhs) {
        lhs.cast().value() = rhs;
        lhs.cast().gradient() = Gradient();
      }

      /*******************************************************************************/

      CODI_INLINE void setGradient(Identifier& identifier, Gradient const& gradient,
                                   AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        CODI_UNUSED(adjointsManagement);
        identifier = gradient;
      }

      CODI_INLINE Gradient const& getGradient(
          Identifier const& identifier, AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) const {
        CODI_UNUSED(adjointsManagement);
        return identifier;
      }

      CODI_INLINE Gradient& gradient(Identifier& identifier,
                                     AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) {
        CODI_UNUSED(adjointsManagement);
        return identifier;
      }

      CODI_INLINE Gradient const& gradient(
          Identifier const& identifier, AdjointsManagement adjointsManagement = AdjointsManagement::Automatic) const {
        CODI_UNUSED(adjointsManagement);
        return identifier;
      }

 
    private:
 
      CODI_INLINE void setGradient(Identifier const& identifier, Gradient const& gradient) {
        CODI_UNUSED(identifier, gradient);
      }
 
      CODI_INLINE Gradient& gradient(Identifier const& identifier) {
        CODI_UNUSED(identifier);
        static Gradient temp = Gradient();
        return temp;
      }
  };

  template<typename T_Type>
  struct RealTraits::IsTotalFinite<T_Type, TapeTraits::EnableIfForwardTape<typename T_Type::Tape>> {
    public:
 
      using Type = CODI_DD(
          T_Type, CODI_T(LhsExpressionInterface<double, double, InternalStatementRecordingTapeInterface<CODI_ANY>,
                                                T_Type>));  

      static CODI_INLINE bool isTotalFinite(Type const& v) {
        using std::isfinite;
        return RealTraits::isTotalFinite(v.getValue()) && RealTraits::isTotalFinite(v.getGradient());
      }
  };

  template<typename T_Type>
  struct RealTraits::IsTotalZero<T_Type, TapeTraits::EnableIfForwardTape<typename T_Type::Tape>> {
    public:
 
      using Type = CODI_DD(
          T_Type, CODI_T(LhsExpressionInterface<double, double, InternalStatementRecordingTapeInterface<CODI_ANY>,
                                                T_Type>));  
      using Real = typename Type::Real;                     

      static CODI_INLINE bool isTotalZero(Type const& v) {
        return Real() == v.getValue() && typename Type::Gradient() == v.getGradient();
      }
  };
}