derivativeAccess.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 <utility>
 
#include "../config.h"
#include "../expressions/lhsExpressionInterface.hpp"
#include "../misc/binomial.hpp"
#include "../misc/compileTimeLoop.hpp"
#include "../misc/exceptions.hpp"
#include "../misc/macros.hpp"
#include "../traits/realTraits.hpp"
 
namespace codi {
 
#ifndef DOXYGEN_DISABLE
 
  namespace DerivativeAccessImpl {
    CODI_INLINE size_t constexpr maximumDerivatives(size_t selectionDepth, size_t order) {
      return binomial(selectionDepth, order);
    }
 
    CODI_INLINE size_t constexpr maximumDerivativesPrimalBranch(size_t selectionDepth, size_t order) {
      return binomial(selectionDepth - 1, order);
    }
 
    CODI_INLINE size_t constexpr isPrimalBranch(size_t selectionDepth, size_t order, size_t l) {
      return l < maximumDerivativesPrimalBranch(selectionDepth, order);
    }
 
    template<typename Type, size_t selectionDepth, size_t order, size_t l>
    struct CheckCompileTimeValues {
      public:
        CODI_STATIC_ASSERT(selectionDepth <= RealTraits::MaxDerivativeOrder<Type>(),
                           "Selection depth must not be higher than the maximum derivative order.");
        CODI_STATIC_ASSERT(order <= selectionDepth, "Derivative order must not be higher than the selection depth.");
        CODI_STATIC_ASSERT(l < maximumDerivatives(selectionDepth, order),
                           "Selected derivative must not be greater than the number of available derivatives for that"
                           "order.");
 
        static bool constexpr isValid = true;
    };
 
    template<typename Type, bool constant, size_t selectionDepth, size_t order, size_t l,
             bool primalBranch = isPrimalBranch(selectionDepth, order, l)>
    struct SelectCompileTime;
 
    template<typename T_Type, bool constant, size_t selectionDepth, size_t order, size_t l>
    struct SelectCompileTime<T_Type, constant, selectionDepth, order, l, true> {
      public:
        using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
 
        using Inner = SelectCompileTime<typename Type::Real, constant, selectionDepth - 1, order, l>;
        using ArgType = typename std::conditional<constant, Type const, Type>::type;
        using RType = typename Inner::RType;
 
        CODI_STATIC_ASSERT(CODI_T(CheckCompileTimeValues<Type, selectionDepth, order, l>::isValid),
                           "Checks inside of type.");
 
        static RType& select(ArgType& value) {
          return Inner::select(value.value());
        }
    };
 
    template<typename T_Type, bool constant, size_t selectionDepth, size_t order, size_t l>
    struct SelectCompileTime<T_Type, constant, selectionDepth, order, l, false> {
      public:
        using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
 
        using Inner = SelectCompileTime<typename Type::Gradient, constant, selectionDepth - 1, order - 1,
                                        l - maximumDerivativesPrimalBranch(selectionDepth, order)>;
        using ArgType = typename std::conditional<constant, Type const, Type>::type;
        using RType = typename Inner::RType;
 
        CODI_STATIC_ASSERT(CODI_T(CheckCompileTimeValues<Type, selectionDepth, order, l>::isValid),
                           "Checks inside of type.");
 
        static RType& select(ArgType& value) {
          return Inner::select(value.gradient());
        }
    };
 
    template<typename Type, bool constant>
    struct SelectCompileTime<Type, constant, 0, 0, 0, true> {
      public:
        using ArgType = typename std::conditional<constant, Type const, Type>::type;
        using RType = ArgType;
 
        static RType& select(ArgType& value) {
          return value;
        }
    };
 
    template<typename T_Type, bool constant, size_t T_selectionDepth>
    struct SelectRunTime {
      public:
        using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
        static size_t constexpr selectionDepth = CODI_DD(T_selectionDepth, CODI_UNDEFINED_VALUE);
 
        CODI_STATIC_ASSERT(CODI_T(std::is_same<typename Type::Real, typename Type::Gradient>::value),
                           "CoDiPack type needs to have the same real and gradient value for run time derivative "
                           "selection.");
        CODI_STATIC_ASSERT(selectionDepth <= RealTraits::MaxDerivativeOrder<Type>(),
                           "Selection depth must not be higher than the maximum derivative order");
 
        using Inner = SelectRunTime<typename Type::Real, constant, selectionDepth - 1>;
        using ArgType = typename std::conditional<constant, Type const, Type>::type;
        using RType = typename Inner::RType;
 
        static RType& select(ArgType& v, size_t order, size_t l) {
          size_t const maxDerivativesPrimalBranch = binomial(selectionDepth - 1, order);
          if (l < maxDerivativesPrimalBranch) {
            return Inner::select(v.value(), order, l);
          } else {
            return Inner::select(v.gradient(), order - 1, l - maxDerivativesPrimalBranch);
          }
        }
    };
 
    template<typename T_Type, bool constant>
    struct SelectRunTime<T_Type, constant, 0> {
      public:
        using Type = CODI_DD(T_Type, double);
        using ArgType = typename std::conditional<constant, Type const, Type>::type;
        using RType = ArgType;
 
        static RType& select(ArgType& v, size_t order, size_t l) {
          CODI_UNUSED(order, l);
          return v;
        }
    };
  }
#endif
 
  template<typename T_Type>
  struct DerivativeAccess {
    public:
 
      using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
 
      template<bool constant, size_t selectionDepth>
      using SelectRunTime = DerivativeAccessImpl::SelectRunTime<Type, constant, selectionDepth>;
 
      template<bool constant, size_t selectionDepth, size_t order, size_t l>
      using SelectCompileTime = DerivativeAccessImpl::SelectCompileTime<Type, constant, selectionDepth, order, l>;
 
      template<size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static typename SelectRunTime<true, selectionDepth>::RType const& derivative(Type const& v, size_t order,
                                                                                   size_t l) {
        checkRuntimeSelection<selectionDepth>(order, l);
 
        return SelectRunTime<true, selectionDepth>::select(v, order, l);
      }
 
      template<size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static typename SelectRunTime<false, selectionDepth>::RType& derivative(Type& v, size_t order, size_t l) {
        checkRuntimeSelection<selectionDepth>(order, l);
 
        return SelectRunTime<false, selectionDepth>::select(v, order, l);
      }
 
      template<typename Derivative, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static void setAllDerivatives(Type& v, size_t order, Derivative const& d) {
        size_t const maxDerivatives = binomial(selectionDepth, order);
        for (size_t i = 0; i < maxDerivatives; i += 1) {
          derivative<selectionDepth>(v, order, i) = d;
        }
      }
 
      template<typename Derivative, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static void setAllDerivativesForward(Type& v, size_t order, Derivative const& d) {
        DerivativeAccess<typename Type::Real>::template setAllDerivatives<Derivative, selectionDepth - 1>(v.value(),
                                                                                                          order, d);
      }
 
      template<typename Derivative, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static void setAllDerivativesReverse(Type& v, size_t order, Derivative const& d) {
        DerivativeAccess<typename Type::Gradient>::template setAllDerivatives<Derivative, selectionDepth - 1>(
            v.gradient(), order - 1, d);
      }
 
      template<size_t order, size_t l, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static typename SelectCompileTime<true, selectionDepth, order, l>::RType const& derivative(Type const& v) {
        return SelectCompileTime<true, selectionDepth, order, l>::select(v);
      }
 
      template<size_t order, size_t l, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static typename SelectCompileTime<false, selectionDepth, order, l>::RType& derivative(Type& v) {
        return SelectCompileTime<false, selectionDepth, order, l>::select(v);
      }
 
      template<size_t order, typename Derivative, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static void setAllDerivatives(Type& v, Derivative const& d) {
        CompileTimeLoop<DerivativeAccessImpl::maximumDerivatives(selectionDepth, order)>::eval(
            CallSetDerivative<order, Derivative, selectionDepth>{}, v, d);
      }
 
      template<size_t order, typename Derivative, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static void setAllDerivativesForward(Type& v, Derivative const& d) {
        DerivativeAccess<typename Type::Real>::template setAllDerivatives<order, Derivative, selectionDepth - 1>(
            v.value(), d);
      }
 
      template<size_t order, typename Derivative, size_t selectionDepth = RealTraits::MaxDerivativeOrder<Type>()>
      static void setAllDerivativesReverse(Type& v, Derivative const& d) {
        DerivativeAccess<typename Type::Gradient>::template setAllDerivatives<order - 1, Derivative,
                                                                              selectionDepth - 1>(v.gradient(), d);
      }
 
    private:
 
      template<size_t selectionDepth>
      static void checkRuntimeSelection(size_t order, size_t l) {
        if (order > selectionDepth) {
          CODI_EXCEPTION(
              "The derivative order must be smaller or equal than the maximum possible derivative. order: %d, max "
              "derivative: %d.",
              order, selectionDepth);
        }
 
        size_t numberDerivatives = binomial(selectionDepth, order);
        if (l >= numberDerivatives) {
          CODI_EXCEPTION(
              "The selected derivative must be smaller than the maximum number of derivatives. selected: %d, number "
              "derivatives: %d.",
              l, numberDerivatives);
        }
      }
 
      template<size_t order, typename Derivative, size_t selectionDepth>
      struct CallSetDerivative {
        public:
          template<size_t pos>
          void operator()(std::integral_constant<size_t, pos>, Type& v, Derivative const& d) {
            DerivativeAccess::derivative<order, pos - 1, selectionDepth>(v) = d;
          }
      };
  };
}