customAdjointVectorHelper.hpp

/*
 * CoDiPack, a Code Differentiation Package
 *
 * Copyright (C) 2015-2024 Chair for Scientific Computing (SciComp), University of Kaiserslautern-Landau
 * Homepage: http://www.scicomp.uni-kl.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://www.scicomp.uni-kl.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/lhsExpressionInterface.hpp"
#include "../../misc/macros.hpp"
#include "../../tapes/interfaces/fullTapeInterface.hpp"
#include "../../tapes/misc/adjointVectorAccess.hpp"
#include "../../tapes/misc/vectorAccessInterface.hpp"
#include "../../traits/tapeTraits.hpp"
 
namespace codi {
 
  template<typename T_Type>
  struct CustomAdjointVectorInterface {
    public:
 
      using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
 
      using Real = typename Type::Real;              
      using Identifier = typename Type::Identifier;  
 
      using Tape = typename Type::Tape;
      using Position = typename Tape::Position;  
 
    protected:
 
      Tape& tape;  
 
    public:
 
      CustomAdjointVectorInterface() : tape(Type::getTape()) {}
 
      virtual ~CustomAdjointVectorInterface() {}
 
      /*******************************************************************************/
 
      virtual void clearAdjoints() = 0;
 
      virtual void deleteAdjointVector() = 0;
 
      virtual void evaluate(Position const& start, Position const& end) = 0;
 
      // clang-format off
      // clang-format on
      virtual void evaluateForward(Position const& start, Position const& end) = 0;
 
      virtual VectorAccessInterface<Real, Identifier>* getVectorInterface() = 0;
 
      /*******************************************************************************/
 
      void evaluate() {
        evaluate(tape.getPosition(), tape.getZeroPosition());
      }
 
      void evaluateForward() {
        evaluate(tape.getPosition(), tape.getZeroPosition());
      }
 
      void setTape(Tape& tape) {
        this->tape = tape;
      }
 
  };
 
  template<typename T_Type, typename T_Gradient>
  struct CustomAdjointVectorHelper : public CustomAdjointVectorInterface<T_Type> {
    public:
 
      using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
      using Gradient = CODI_DD(T_Gradient, double);  
 
      using Base = CustomAdjointVectorInterface<Type>;  
 
      using Real = typename Type::Real;              
      using Identifier = typename Type::Identifier;  
 
      using Tape = CODI_DD(typename Type::Tape, CODI_DEFAULT_TAPE);
      using Position = typename Tape::Position;  
 
    protected:
      std::vector<Gradient> adjointVector;  
 
      Gradient zeroValue;             
      Gradient const constZeroValue;  
 
      AdjointVectorAccess<Real, Identifier, Gradient>* adjointInterface;  
 
    public:
 
      CustomAdjointVectorHelper()
          : Base(), adjointVector(0), zeroValue(), constZeroValue(), adjointInterface(nullptr) {}
 
      ~CustomAdjointVectorHelper() {
        if (nullptr != adjointInterface) {
          delete adjointInterface;
        }
      }
 
      /*******************************************************************************/
 
      void clearAdjoints() {
        for (size_t i = 0; i < adjointVector.size(); i += 1) {
          adjointVector[i] = Gradient();
        }
      }
 
      void deleteAdjointVector() {
        adjointVector.resize(0);
        adjointVector.shrink_to_fit();
      }
 
      void evaluate(Position const& start, Position const& end) {
        checkAdjointVectorSize();
 
        Base::tape.evaluate(start, end, adjointVector.data());
      }
      using Base::evaluate;
 
      void evaluateForward(Position const& start, Position const& end) {
        checkAdjointVectorSize();
 
        Base::tape.evaluateForward(start, end, adjointVector.data());
      }
      using Base::evaluateForward;
 
      VectorAccessInterface<Real, Identifier>* getVectorInterface() {
        if (nullptr != adjointInterface) {
          delete adjointInterface;
        }
 
        checkAdjointVectorSize();
        adjointInterface = new AdjointVectorAccess<Real, Identifier, Gradient>(adjointVector.data());
        return adjointInterface;
      }
 
      /*******************************************************************************/
 
      Gradient const& getGradient(Identifier const& identifier) const {
        return gradient(identifier);
      }
 
      Gradient& gradientUnchecked(Identifier const& identifier) {
        return adjointVector[identifier];
      }
 
      Gradient const& gradientUnchecked(Identifier const& identifier) const {
        return adjointVector[identifier];
      }
 
      Gradient& gradient(Identifier const& identifier) {
        checkAdjointVectorSize();
 
        if (0 != identifier && identifier < (Identifier)adjointVector.size()) {
          return adjointVector[identifier];
        } else {
          zeroValue = Gradient();
          return zeroValue;
        }
      }
 
      Gradient const& gradient(Identifier const& identifier) const {
        if (0 != identifier && identifier < (Identifier)adjointVector.size()) {
          return adjointVector[identifier];
        } else {
          return constZeroValue;
        }
      }
 
      void setGradient(Identifier& identifier, Gradient const& gradientValue) {
        gradient(identifier) = gradientValue;
      }
 
 
    private:
 
      void checkAdjointVectorSize() {
        if (adjointVector.size() <= Base::tape.getParameter(TapeParameters::LargestIdentifier)) {
          adjointVector.resize(Base::tape.getParameter(TapeParameters::LargestIdentifier) + 1);
        }
      }
  };
}