ModPiece.h

Go to the documentation of this file.

#ifndef MODPIECE_H
#define MODPIECE_H
 
#include "MUQ/Modeling/WorkPiece.h"
#include "MUQ/Utilities/VariadicMacros.h"
 
#include <Eigen/Core>
#include <vector>
 
namespace muq{
  namespace Modeling{
 
  class ModPiece : public WorkPiece
  {
 
  public:
 
    ModPiece(std::vector<int> const& inputSizes,
             std::vector<int> const& outputSizes);
 
    ModPiece(Eigen::VectorXi const& inputSizes,
             Eigen::VectorXi const& outputSizes);
 
    virtual ~ModPiece() = default;
 
 
    virtual double GetRunTime(const std::string& method="Evaluate") const override;
    virtual void ResetCallTime() override;
 
    virtual unsigned long int GetNumCalls(const std::string& method = "Evaluate") const override;
 
    virtual std::vector<Eigen::VectorXd> const& Evaluate(std::vector<Eigen::VectorXd> const& input);
    virtual std::vector<Eigen::VectorXd> const& Evaluate(ref_vector<Eigen::VectorXd> const& input);
    VARIADIC_TO_REFVECTOR(Evaluate, Eigen::VectorXd, std::vector<Eigen::VectorXd> const&);
 
    virtual Eigen::VectorXd const& Gradient(unsigned int                 const  outputDimWrt,
                                            unsigned int                 const  inputDimWrt,
                                            std::vector<Eigen::VectorXd> const& input,
                                            Eigen::VectorXd              const& sensitivity);
 
    virtual Eigen::VectorXd const& Gradient(unsigned int                const  outputDimWrt,
                                            unsigned int                const  inputDimWrt,
                                            ref_vector<Eigen::VectorXd> const& input,
                                            Eigen::VectorXd             const& sensitivity);
 
    inline Eigen::VectorXd const& Gradient(unsigned int outWrt,
                                           unsigned int inWrt,
                                           Eigen::VectorXd const& last,
                                           Eigen::VectorXd const& sens) {
      ref_vector<Eigen::VectorXd> vec;
      vec.push_back(std::cref(last));
      return Gradient(outWrt, inWrt, vec, sens);
    }
 
    template<typename... Args>
    inline Eigen::VectorXd  const& Gradient(unsigned int wrtOut,
                                            unsigned int wrtIn,
                                            Args const&... args) {
      ref_vector<Eigen::VectorXd> vec;
      return GradientRecurse(wrtOut, wrtIn, vec, args...);
    }
 
 
    inline Eigen::VectorXd const& ApplyHessian(unsigned int outWrt,
                                           unsigned int inWrt1,
                                           unsigned int inWrt2,
                                           Eigen::VectorXd const& last,
                                           Eigen::VectorXd const& sens,
                                           Eigen::VectorXd const& vec) {
      ref_vector<Eigen::VectorXd> invec;
      invec.push_back(std::cref(last));
      return ApplyHessian(outWrt, inWrt1, inWrt2, invec, sens, vec);
    }
 
    template<typename... Args>
    inline Eigen::VectorXd  const& ApplyHessian(unsigned int wrtOut,
                                            unsigned int wrtIn1,
                                            unsigned int wrtIn2,
                                            Args const&... args) {
      ref_vector<Eigen::VectorXd> invec;
      return ApplyHessianRecurse(wrtOut, wrtIn1, wrtIn2, invec, args...);
    }
 
    virtual Eigen::MatrixXd const& Jacobian(unsigned int                 const  outputDimWrt,
                                            unsigned int                 const  inputDimWrt,
                                            std::vector<Eigen::VectorXd> const& input);
 
    virtual Eigen::MatrixXd const& Jacobian(unsigned int                const  outputDimWrt,
                                            unsigned int                const  inputDimWrt,
                                            ref_vector<Eigen::VectorXd> const& input);
 
 
    template<typename... Args>
    inline Eigen::MatrixXd const& Jacobian(unsigned int outWrt, unsigned int inWrt, Args const&... args) {
      ref_vector<Eigen::VectorXd> vec;
      return Jacobian(outWrt, inWrt, vec, args...);
    }
 
    template<typename... Args>
    inline Eigen::MatrixXd JacobianByFD(unsigned int outWrt, unsigned int inWrt, Args const&... args) {
      ref_vector<Eigen::VectorXd> vec;
      return JacobianByFD(outWrt, inWrt, vec, args...);
    }
 
    template<typename... Args>
    inline Eigen::MatrixXd ApplyJacobianByFD(unsigned int outWrt, unsigned int inWrt, Args const&... args) {
      ref_vector<Eigen::VectorXd> vec;
      return ApplyJacobianByFD(outWrt, inWrt, vec, args...);
    }
 
 
    virtual Eigen::VectorXd const& ApplyJacobian(unsigned int                 const  outputDimWrt,
                                                 unsigned int                 const  inputDimWrt,
                                                 std::vector<Eigen::VectorXd> const& input,
                                                 Eigen::VectorXd              const& vec);
 
    virtual Eigen::VectorXd const& ApplyJacobian(unsigned int                const  outputDimWrt,
                                                 unsigned int                const  inputDimWrt,
                                                 ref_vector<Eigen::VectorXd> const& input,
                                                 Eigen::VectorXd             const& vec);
 
    inline Eigen::VectorXd const& ApplyJacobian(unsigned int outWrt, unsigned int inWrt, Eigen::VectorXd const& last, Eigen::VectorXd const& sens) {     \
      ref_vector<Eigen::VectorXd> vec;
      vec.push_back(std::cref(last));                                                                                               \
      return ApplyJacobian(outWrt, inWrt, vec, sens);                                                                              \
    }
    template<typename... Args>
    inline Eigen::VectorXd  const& ApplyJacobian(unsigned int wrtOut, unsigned int wrtIn, Args const&... args) {
      ref_vector<Eigen::VectorXd> vec;
      return ApplyJacobianRecurse(wrtOut, wrtIn, vec, args...);
    }
 
 
    virtual Eigen::VectorXd GradientByFD(unsigned int                 const  outputDimWrt,
                                         unsigned int                 const  inputDimWrt,
                                         std::vector<Eigen::VectorXd> const& input,
                                         Eigen::VectorXd              const& sensitivity);
 
    virtual Eigen::VectorXd GradientByFD(unsigned int                const  outputDimWrt,
                                         unsigned int                const  inputDimWrt,
                                         ref_vector<Eigen::VectorXd> const& input,
                                         Eigen::VectorXd             const& sensitivity);
 
    virtual Eigen::MatrixXd JacobianByFD(unsigned int                 const  outputDimWrt,
                                         unsigned int                 const  inputDimWrt,
                                         std::vector<Eigen::VectorXd> const& input);
 
    virtual Eigen::MatrixXd JacobianByFD(unsigned int                const  outputDimWrt,
                                         unsigned int                const  inputDimWrt,
                                         ref_vector<Eigen::VectorXd> const& input);
 
    virtual Eigen::VectorXd ApplyJacobianByFD(unsigned int                 const  outputDimWrt,
                                              unsigned int                 const  inputDimWrt,
                                              std::vector<Eigen::VectorXd> const& input,
                                              Eigen::VectorXd              const& vec);
 
    virtual Eigen::VectorXd ApplyJacobianByFD(unsigned int                const  outputDimWrt,
                                              unsigned int                const  inputDimWrt,
                                              ref_vector<Eigen::VectorXd> const& input,
                                              Eigen::VectorXd             const& vec);
 
    virtual Eigen::VectorXd const& ApplyHessian(unsigned int const outWrt,
                                         unsigned int const inWrt1,
                                         unsigned int const inWrt2,
                                         std::vector<Eigen::VectorXd> const& input,
                                         Eigen::VectorXd              const& sens,
                                         Eigen::VectorXd              const& vec);
 
    virtual Eigen::VectorXd const& ApplyHessian(unsigned int const outWrt,
                                         unsigned int const inWrt1,
                                         unsigned int const inWrt2,
                                         ref_vector<Eigen::VectorXd> const& input,
                                         Eigen::VectorXd             const& sens,
                                         Eigen::VectorXd             const& vec);
 
   virtual Eigen::VectorXd ApplyHessianByFD(unsigned int                const  outWrt,
                                            unsigned int                const  inWrt1,
                                            unsigned int                const  inWrt2,
                                            std::vector<Eigen::VectorXd> const& input,
                                            Eigen::VectorXd             const& sens,
                                            Eigen::VectorXd             const& vec);
 
   virtual Eigen::VectorXd ApplyHessianByFD(unsigned int                const  outWrt,
                                            unsigned int                const  inWrt1,
                                            unsigned int                const  inWrt2,
                                            ref_vector<Eigen::VectorXd> const& input,
                                            Eigen::VectorXd             const& sens,
                                            Eigen::VectorXd             const& vec);
 
    void EnableCache(){cacheEnabled=true;};
    void DisableCache(){cacheEnabled=false;};
    bool CacheStatus() const{return cacheEnabled;};
 
 
    virtual void SetWarnLevel(unsigned int newLevel){fdWarnLevel=newLevel;};
 
    const Eigen::VectorXi inputSizes;
    const Eigen::VectorXi outputSizes;
 
  protected:
 
    std::vector<Eigen::VectorXd> ToStdVec(ref_vector<Eigen::VectorXd> const& input);
 
    // Should one step caching be used?
    bool cacheEnabled = false;
    std::vector<Eigen::VectorXd> cacheInput;
 
    bool ExistsInCache(ref_vector<Eigen::VectorXd> const& input) const;
 
    // The following variables keep track of how many times the Implemented functions, i.e. EvaluateImpl, GradientImpl,
    // etc... are called.
    unsigned long int numGradCalls    = 0;
    unsigned long int numJacCalls     = 0;
    unsigned long int numJacActCalls  = 0;
    unsigned long int numHessActCalls = 0;
    unsigned long int numGradFDCalls  = 0;
    unsigned long int numJacFDCalls  = 0;
    unsigned long int numJacActFDCalls = 0;
    unsigned long int numHessActFDCalls = 0;
 
    // these variables keep track of the total wall-clock time spent in each of the Implemented functions.  They are in
    // units of milliseconds
    double gradTime    = 0;
    double jacTime     = 0;
    double jacActTime  = 0;
    double hessActTime = 0;
 
    std::vector<Eigen::VectorXd> outputs;
    Eigen::VectorXd gradient;
    Eigen::VectorXd jacobianAction;
    Eigen::MatrixXd jacobian;
    Eigen::VectorXd hessAction;
 
    unsigned int fdWarnLevel = 0;
 
    void CheckInputs(ref_vector<Eigen::VectorXd> const& input, std::string const& funcName);
 
    virtual void EvaluateImpl(ref_vector<boost::any> const& inputs) override;
 
    virtual void EvaluateImpl(ref_vector<Eigen::VectorXd> const& input) = 0;
 
    virtual void GradientImpl(unsigned int                const  outputDimWrt,
                              unsigned int                const  inputDimWrt,
                              ref_vector<Eigen::VectorXd> const& input,
                              Eigen::VectorXd             const& sensitivity);
 
    virtual void JacobianImpl(unsigned int                const  outputDimWrt,
                              unsigned int                const  inputDimWrt,
                              ref_vector<Eigen::VectorXd> const& input);
 
    virtual void ApplyJacobianImpl(unsigned int                const  outputDimWrt,
                                   unsigned int                const  inputDimWrt,
                                   ref_vector<Eigen::VectorXd> const& input,
                                   Eigen::VectorXd             const& vec);
 
    virtual void ApplyHessianImpl(unsigned int                const  outWrt,
                                  unsigned int                const  inWrt1,
                                  unsigned int                const  inWrt2,
                                  ref_vector<Eigen::VectorXd> const& input,
                                  Eigen::VectorXd             const& sens,
                                  Eigen::VectorXd             const& vec);
 
    // virtual void ApplyHessianImpl(int                         const  outputDimWrt,
    //                       int                         const  inputDimWrt1,
    //                       int                         const  inputDimWrt2,
    //                       ref_vector<Eigen::VectorXd> const& input,
    //                       Eigen::VectorXd             const& sensitivity
    //                       Eigen::VectorXd             const& vec);
 
 
    // virtual void ApplyHessianByFD(int                         const  outputDimWrt,
    //                       int                         const  inputDimWrt1,
    //                       int                         const  inputDimWrt2,
    //                       ref_vector<Eigen::VectorXd> const& input,
    //                       Eigen::VectorXd             const& sensitivity
    //                       Eigen::VectorXd             const& vec);
 
  private:
 
    template<typename NextType, typename... Args>
    inline Eigen::VectorXd const& GradientRecurse(unsigned int outWrt,
                                                  unsigned int inWrt,
                                                  ref_vector<Eigen::VectorXd>& vec,
                                                  NextType const& ith,
                                                  Args const&... args) {
 
      static_assert(std::is_same<Eigen::VectorXd,
                    NextType>::value,
                    "In ModPiece::Gradient, cannot cast input to Eigen::VectorXd.");
 
        vec.push_back(std::cref((NextType&)ith));
        return GradientRecurse(outWrt, inWrt, vec, args...);
 
    }
 
    template<typename NextType>
    inline Eigen::VectorXd const& GradientRecurse(unsigned int outWrt,
                                                  unsigned int inWrt,
                                                  ref_vector<Eigen::VectorXd>& vec,
                                                  NextType const& last,
                                                  Eigen::VectorXd const& sens) {
        static_assert(std::is_same<Eigen::VectorXd, NextType>::value, "In ModPiece::Gradient, cannot cast input to Eigen::VectorXd.");
        vec.push_back(std::cref((NextType&)last));
        return Gradient(outWrt, inWrt, vec, sens);
    }
 
    template<typename NextType, typename... Args>
    inline Eigen::VectorXd const& ApplyHessianRecurse(unsigned int outWrt,
                                                      unsigned int inWrt1,
                                                      unsigned int inWrt2,
                                                      ref_vector<Eigen::VectorXd>& invec,
                                                      NextType const& ith,
                                                      Args const&... args) {
 
      static_assert(std::is_same<Eigen::VectorXd,
                  NextType>::value,
                  "In ModPiece::ApplyHessian, cannot cast input to Eigen::VectorXd.");
 
      invec.push_back(std::cref((NextType&)ith));
      return ApplyHessianRecurse(outWrt, inWrt1, inWrt2, invec, args...);
 
    }
 
    template<typename NextType>
    inline Eigen::VectorXd const& ApplyHessianRecurse(unsigned int outWrt,
                                                  unsigned int inWrt1,
                                                  unsigned int inWrt2,
                                                  ref_vector<Eigen::VectorXd>& invec,
                                                  NextType const& last,
                                                  Eigen::VectorXd const& sens,
                                                  Eigen::VectorXd const& vec) {
        static_assert(std::is_same<Eigen::VectorXd, NextType>::value, "In ModPiece::ApplyHessian, cannot cast input to Eigen::VectorXd.");
        invec.push_back(std::cref((NextType&)last));
        return ApplyHessian(outWrt, inWrt1, inWrt2, invec, sens, vec);
    }
 
 
    template<typename... Args>
    inline Eigen::MatrixXd const& Jacobian(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& ith, Args const&... args) {     \
        vec.push_back(std::cref(ith));                                                                                               \
        return Jacobian(outWrt, inWrt, vec, args...);                                                                              \
    }
    inline Eigen::MatrixXd const& Jacobian(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& last) {
      vec.push_back(std::cref(last));
      return Jacobian(outWrt, inWrt, vec);
    }
 
    template<typename NextType, typename... Args>                                                                                               \
    inline Eigen::VectorXd const& ApplyJacobianRecurse(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, NextType const& ith, Args const&... args) {     \
        static_assert(std::is_same<Eigen::VectorXd, NextType>::value, "In ModPiece::Gradient, cannot cast input to Eigen::VectorXd."); \
        vec.push_back(std::cref((NextType&)ith));                                                                                               \
        return ApplyJacobianRecurse(outWrt, inWrt, vec, args...);                                                                              \
    }
 
    template<typename NextType>                                                                                               \
    inline Eigen::VectorXd const& ApplyJacobianRecurse(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, NextType const& last, Eigen::VectorXd const& sens) {     \
        static_assert(std::is_same<Eigen::VectorXd, NextType>::value, "In ModPiece::Gradient, cannot cast input to Eigen::VectorXd."); \
        vec.push_back(std::cref((NextType&)last));                                                                                               \
        return ApplyJacobian(outWrt, inWrt, vec, sens);                                                                              \
    }
 
    template<typename... Args>                                                                                               \
    inline Eigen::MatrixXd JacobianByFD(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& ith, Args const&... args) {     \
        vec.push_back(std::cref(ith));                                                                                               \
        return JacobianByFD(outWrt, inWrt, vec, args...);                                                                              \
    }
    inline Eigen::MatrixXd JacobianByFD(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& last) {
      vec.push_back(std::cref(last));
      return JacobianByFD(outWrt, inWrt, vec);
    }
 
    template<typename NextType, typename... Args>
    inline Eigen::MatrixXd ApplyJacobianByFD(unsigned int outWrt, unsigned int  inWrt, ref_vector<Eigen::VectorXd>& vec, NextType const& ith, Args const&... args) {
      vec.push_back(std::cref(ith));
      return ApplyJacobianByFD(outWrt, inWrt, vec, args...);
    }
    template<typename NextType>
    inline Eigen::MatrixXd ApplyJacobianByFD(unsigned int outWrt, unsigned int inWrt, ref_vector<Eigen::VectorXd>& vec, NextType const& last, Eigen::VectorXd const& sens) {
      vec.push_back(std::cref(last));
      return ApplyJacobianByFD(outWrt, inWrt, vec, sens);
    }
 
  };
 
  }
}
 
 
 
#endif // #ifndef MODPIECE_H