LinearOperator.h

Go to the documentation of this file.

#ifndef LINEAROPERATOR_H
#define LINEAROPERATOR_H
 
#include "MUQ/Modeling/ModPiece.h"
 
#include <Eigen/Core>
 
#include <memory>
#include <iostream>
#include <exception>
 
#include "MUQ/Utilities/Demangler.h"
 
namespace muq
{
namespace Modeling
{
 
 
class LinearOperatorTypeException: public std::exception
{   
public:
  LinearOperatorTypeException(std::string const& type){
    msg = "Tried creating a linear operator from an unsupported type " + type + ".  Make sure all necessary headers are included and a child of LinearOperator exists for this type.";
  };
 
  virtual const char* what() const throw()
  {
    return msg.c_str();
  }
 
  std::string msg;
};
 
class LinearOperator;
 
template<typename MatrixType>
struct LinearOperatorFactory
{
    static std::shared_ptr<LinearOperator> Create(MatrixType const& A)
    {
        throw LinearOperatorTypeException(muq::Utilities::demangle(typeid(A).name()));
 
        return std::shared_ptr<LinearOperator>();
 
    };
};
 
class LinearOperator : public muq::Modeling::ModPiece{
public:
 
  LinearOperator(int rowsIn, int colsIn, int numInputCols=1);
 
  virtual ~LinearOperator(){};
 
  virtual Eigen::MatrixXd Apply(Eigen::Ref<const Eigen::MatrixXd> const& x) = 0;
 
  virtual Eigen::MatrixXd ApplyTranspose(Eigen::Ref<const Eigen::MatrixXd> const& x) = 0;
 
  virtual void Apply(Eigen::Ref<const Eigen::MatrixXd> const& x, Eigen::Ref<Eigen::MatrixXd> y);
 
  virtual void ApplyTranspose(Eigen::Ref<const Eigen::MatrixXd> const& x, Eigen::Ref<Eigen::MatrixXd> y);
 
  int rows() const { return nrows; }
 
  int cols() const { return ncols; }
 
  virtual Eigen::MatrixXd GetMatrix();
 
  template<typename OtherType>
  static std::shared_ptr<LinearOperator> Create(OtherType const& A)
  {
      return LinearOperatorFactory<OtherType>::Create(A);
  }
 
protected:
 
  const int ncols;
  const int nrows;
 
  // ModPiece overrides
  virtual void EvaluateImpl(muq::Modeling::ref_vector<Eigen::VectorXd> const& input) override;
 
  virtual void GradientImpl(unsigned int                const  outputDimWrt,
                            unsigned int                const  inputDimWrt,
                            muq::Modeling::ref_vector<Eigen::VectorXd> const& input,
                            Eigen::VectorXd             const& sensitivity) override;
 
  virtual void JacobianImpl(unsigned int                const  outputDimWrt,
                            unsigned int                const  inputDimWrt,
                            muq::Modeling::ref_vector<Eigen::VectorXd> const& input) override;
 
  virtual void ApplyJacobianImpl(unsigned int                const  outputDimWrt,
                                 unsigned int                const  inputDimWrt,
                                 muq::Modeling::ref_vector<Eigen::VectorXd> const& input,
                                 Eigen::VectorXd             const& vec) override;
 
  virtual void ApplyHessianImpl(unsigned int const outWrt,
                                unsigned int inWrt1,
                                unsigned int inWrt2,
                                muq::Modeling::ref_vector<Eigen::VectorXd> const& input,
                                Eigen::VectorXd             const& sens,
                                Eigen::VectorXd             const& vec) override;
};
 
 
 
 
} // namespace Modeling
} // namespace MUQ
 
#endif