analytic_function_operator.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
#include <kernel/base_header.hpp>
#include <kernel/analytic/function.hpp>
 
#include <deque>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename DT_, typename IT_, typename Function_>
    class AnalyticFunctionOperator
    {
      public:
        static_assert(std::is_same<typename Function_::ImageType, Analytic::Image::Scalar>::value,
        "AnalyticFunctionOperator is implemented for scalar functions only");
 
        typedef DT_ DataType;
        typedef IT_ IndexType;
 
        typedef Function_ FunctionType;
 
        static constexpr bool can_hess = FunctionType::can_hess;
        static constexpr bool can_diag_hess = can_hess;
 
        typedef Analytic::EvalTraits<DT_, Function_> FuncEvalTraits;
        typedef typename Function_::template Evaluator<FuncEvalTraits> EvalType;
        typedef typename EvalType::PointType PointType;
 
        static constexpr int dim = PointType::n;
 
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, dim> VectorTypeR;
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, dim> VectorTypeL;
        typedef VectorTypeR GradientType;
        typedef typename EvalType::HessianType HessianType;
 
      private:
        Function_& _function;
        EvalType _func_eval;
 
        PointType _my_state;
        typename EvalType::GradientType _my_grad;
        typename EvalType::HessianType _my_hess;
 
        Index _num_func_evals;
        Index _num_grad_evals;
        Index _num_hess_evals;
 
      public:
        explicit AnalyticFunctionOperator(Function_& func) :
          _function(func),
          _func_eval(func),
          _my_state(DT_(0)),
          _my_grad(DT_(0)),
          _my_hess(DT_(0)),
          _num_func_evals(Index(0)),
          _num_grad_evals(Index(0)),
          _num_hess_evals(Index(0))
          {
          }
 
        static String name()
        {
          return "AnalyticFunctionOperator";
        }
 
        VectorTypeL create_vector_l() const
        {
          return VectorTypeL(IT_(1));
        }
 
        VectorTypeR create_vector_r() const
        {
          return VectorTypeR(IT_(1));
        }
 
        Index columns()
        {
          return Index(dim);
        }
 
        Index rows()
        {
          return Index(dim);
        }
 
        Index get_num_func_evals() const
        {
          return _num_func_evals;
        }
 
        Index get_num_grad_evals() const
        {
          return _num_grad_evals;
        }
 
        Index get_num_hess_evals() const
        {
          return _num_hess_evals;
        }
 
        void reset_num_evals()
        {
          _num_func_evals = Index(0);
          _num_grad_evals = Index(0);
          _num_hess_evals = Index(0);
        }
 
        template<typename FilterType_>
        void prepare(const VectorTypeR& vec_state, FilterType_ DOXY(filter))
        {
          _my_state = vec_state(0);
        }
 
        void eval_fval_grad(typename EvalType::ValueType& fval, GradientType& grad)
        {
          ++_num_func_evals;
          ++_num_grad_evals;
          fval = _func_eval.value(_my_state);
          _my_grad = _func_eval.gradient(_my_state);
          grad(0, _my_grad);
        }
 
        void apply_hess(VectorTypeL& vec_out, const VectorTypeR& vec_in)
        {
          ++_num_hess_evals;
          _my_hess = _func_eval.hessian(_my_state);
          vec_out(0, _my_hess*vec_in(0));
        }
 
        void compute_hess(HessianType& mat_out)
        {
          ++_num_hess_evals;
          mat_out = _func_eval.hessian(_my_state);
        }
 
        void compute_approx_hess(HessianType& mat_out)
        {
          ++_num_hess_evals;
          mat_out = _func_eval.hessian(_my_state);
          for(int i(1); i < mat_out.m; ++i)
          {
            for(int j(0); j < i; ++j)
            {
              mat_out(i,j) = DT_(0);
              mat_out(j,i) = DT_(0);
            }
          }
        }
    };
  } // namespace Solver
} // namespace FEAT