linear_transform.h

Go to the documentation of this file.

#ifndef SOPT_OPERATORS_H
#define SOPT_OPERATORS_H
 
#include "sopt/config.h"
#include <array>
#include <memory>
#include <type_traits>
#include <utility>  // for std::move<>
#include <Eigen/Core>
#include "sopt/logging.h"
#include "sopt/maths.h"
#include "sopt/types.h"
#include "sopt/wrapper.h"
 
namespace sopt {
 
namespace details {
template <typename EIGEN>
class MatrixToLinearTransform;
template <typename EIGEN>
class MatrixAdjointToLinearTransform;
}  // namespace details
 
template <typename VECTOR>
class LinearTransform : public details::WrapFunction<VECTOR> {
 public:
  using t_Function = OperatorFunction<VECTOR>;
 
  LinearTransform(t_Function const &direct, t_Function const &indirect,
                  std::array<t_int, 3> sizes = {{1, 1, 0}})
      : LinearTransform(
            direct, sizes, indirect,
            {{sizes[1], sizes[0], sizes[0] == 0 ? 0 : -(sizes[2] * sizes[1]) / sizes[0]}}) {
    assert(sizes[0] != 0);
  }
  LinearTransform(t_Function const &direct, std::array<t_int, 3> dsizes, t_Function const &indirect,
                  std::array<t_int, 3> isizes)
      : LinearTransform(details::wrap(direct, dsizes), details::wrap(indirect, isizes)) {}
  LinearTransform(details::WrapFunction<VECTOR> const &direct,
                  details::WrapFunction<VECTOR> const &indirect)
      : details::WrapFunction<VECTOR>(direct), indirect_(indirect) {}
  LinearTransform(LinearTransform const &c)
      : details::WrapFunction<VECTOR>(c), indirect_(c.indirect_), norm_(c.norm_), sq_norm_(c.sq_norm_) {}
  LinearTransform(LinearTransform &&c)
      : details::WrapFunction<VECTOR>(std::move(c)), indirect_(std::move(c.indirect_)), norm_(c.norm_), sq_norm_(c.sq_norm_) {}
  void operator=(LinearTransform const &c) {
    details::WrapFunction<VECTOR>::operator=(c);
    indirect_ = c.indirect_;
    norm_ = c.norm_;
    sq_norm_ = c.sq_norm_;
  }
  void operator=(LinearTransform &&c) {
    details::WrapFunction<VECTOR>::operator=(std::move(c));
    indirect_ = std::move(c.indirect_);
    norm_ = c.norm_;
    sq_norm_ = c.sq_norm_;
  }
 
  LinearTransform<VECTOR> adjoint() const {
    return LinearTransform<VECTOR>(indirect_, static_cast<details::WrapFunction<VECTOR> const &>(*this));
  }
 
  using details::WrapFunction<VECTOR>::operator*;
  using details::WrapFunction<VECTOR>::sizes;
  using details::WrapFunction<VECTOR>::rows;
 
  void set_norm(t_real n)
  {
    norm_ = n;
    sq_norm_ = n*n;
  }
 
  sopt::t_real norm() const
  {
    return norm_;
  }
 
  sopt::t_real sq_norm() const
  {
    return sq_norm_;
  }
 
 private:
  details::WrapFunction<VECTOR> indirect_;
 
  sopt::t_real norm_ = 1.0;
  sopt::t_real sq_norm_ = 1.0;
};
 
template <typename VECTOR>
LinearTransform<VECTOR> linear_transform(OperatorFunction<VECTOR> const &direct,
                                         OperatorFunction<VECTOR> const &indirect,
                                         std::array<t_int, 3> const &sizes = {{1, 1, 0}}) {
  return {direct, indirect, sizes};
}
template <typename VECTOR>
LinearTransform<VECTOR> linear_transform(OperatorFunction<VECTOR> const &direct,
                                         std::array<t_int, 3> const &dsizes,
                                         OperatorFunction<VECTOR> const &indirect,
                                         std::array<t_int, 3> const &isizes) {
  return {direct, dsizes, indirect, isizes};
}
 
template <typename VECTOR>
LinearTransform<VECTOR> &linear_transform(LinearTransform<VECTOR> &passthrough) {
  return passthrough;
}
template <typename VECTOR>
LinearTransform<VECTOR> linear_transform(details::WrapFunction<VECTOR> const &direct,
                                         details::WrapFunction<VECTOR> const &adjoint) {
  return {direct, adjoint};
}
 
namespace details {
 
template <typename EIGEN>
class MatrixToLinearTransform {
  using Raw = typename std::remove_const<typename std::remove_reference<EIGEN>::type>::type;
  using PlainMatrix = typename Raw::PlainObject;
 
 public:
  using PlainObject =
      typename std::conditional<std::is_base_of<Eigen::MatrixBase<PlainMatrix>, PlainMatrix>::value,
                                Vector<typename PlainMatrix::Scalar>,
                                Array<typename PlainMatrix::Scalar>>::type;
  template <typename T0>
  MatrixToLinearTransform(Eigen::MatrixBase<T0> const &A) : matrix(std::make_shared<EIGEN>(A)) {}
  MatrixToLinearTransform(std::shared_ptr<EIGEN> const &x) : matrix(x){}
 
  void operator()(PlainObject &out, PlainObject const &x) const {
#ifndef NDEBUG
    if ((*matrix).cols() != x.size())
      SOPT_THROW("Input vector and matrix do not match: ")
          << out.cols() << " columns for " << x.size() << " elements.";
#endif
    out = (*matrix) * x;
  }
  MatrixAdjointToLinearTransform<EIGEN> adjoint() const {
        return MatrixAdjointToLinearTransform<EIGEN>(matrix);
  }
 
 private:
  std::shared_ptr<EIGEN> matrix;
};
 
template <typename EIGEN>
class MatrixAdjointToLinearTransform {
 public:
  using PlainObject = typename MatrixToLinearTransform<EIGEN>::PlainObject;
  template <typename T0>
  MatrixAdjointToLinearTransform(Eigen::MatrixBase<T0> const &A)
      : matrix(std::make_shared<EIGEN>(A)) {}
  MatrixAdjointToLinearTransform(std::shared_ptr<EIGEN> const &x) : matrix(x){}
 
  void operator()(PlainObject &out, PlainObject const &x) const {
#ifndef NDEBUG
    if ((*matrix).rows() != x.size())
      SOPT_THROW("Input vector and matrix adjoint do not match: ")
          << out.cols() << " rows for " << x.size() << " elements.";
#endif
    out = matrix->adjoint() * x;
  }
  MatrixToLinearTransform<EIGEN> adjoint() const { return MatrixToLinearTransform<EIGEN>(matrix); }
 
 private:
  std::shared_ptr<EIGEN> matrix;
};
}  // namespace details
 
template <typename DERIVED>
LinearTransform<Vector<typename DERIVED::Scalar>> linear_transform(
    Eigen::MatrixBase<DERIVED> const &A) {
  details::MatrixToLinearTransform<Matrix<typename DERIVED::Scalar>> const matrix(A);
  if (A.rows() == A.cols())
    return {matrix, matrix.adjoint()};
  else {
    t_int const gcd = details::gcd(A.cols(), A.rows());
    t_int const a = A.cols() / gcd;
    t_int const b = A.rows() / gcd;
    return {matrix, matrix.adjoint(), {{b, a, 0}}};
  }
}
 
template <typename SCALAR>
LinearTransform<Vector<SCALAR>> linear_transform_identity() {
  return {[](Vector<SCALAR> &out, Vector<SCALAR> const &in) { out = in; },
          [](Vector<SCALAR> &out, Vector<SCALAR> const &in) { out = in; }};
}
}  // namespace sopt
#endif