eigenutil.h

Go to the documentation of this file.

/* This file is part of the Tomographer project, which is distributed under the
 * terms of the MIT license.
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
 * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef TOMOGRAPHER_TOOLS_EIGENUTIL_H
#define TOMOGRAPHER_TOOLS_EIGENUTIL_H

#include <complex>
#include <vector>

#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/StdVector>

#include <boost/serialization/serialization.hpp>
#include <boost/serialization/split_free.hpp>
#include <boost/serialization/complex.hpp>


// -----------------------------------------------------------------------------
// Eigen and std::vector
// -----------------------------------------------------------------------------

namespace Tomographer {
namespace Tools {

template<typename EigenType>
struct EigenStdVector
{
  typedef std::vector<EigenType, Eigen::aligned_allocator<EigenType> > type;
};

} // namespace Tools
} // namespace Tomographer

// -----------------------------------------------------------------------------
// Random matrices in Eigen
// -----------------------------------------------------------------------------

namespace Tomographer {
namespace Tools {

namespace tomo_internal {

  template<typename Rng, typename RndDist, typename Scalar>
  struct random_generator
  {
    typedef Scalar result_type;

    Rng & rng;
    RndDist & rnddist;

    random_generator(Rng & rng_, RndDist & rnddist_)
      : rng(rng_), rnddist(rnddist_)
    {
    }

    template<typename Index>
    inline const result_type operator() (Index, Index = 0) const {
      return rnddist(rng);
    }
  };

  template<typename Rng, typename RndDist, typename RealScalar>
  struct random_generator<Rng, RndDist, std::complex<RealScalar> >
  {
    typedef std::complex<RealScalar> result_type;

    Rng & rng;
    RndDist & rnddist;

    random_generator(Rng & rng_, RndDist & rnddist_)
      : rng(rng_), rnddist(rnddist_)
    {
    }

    template<typename Index>
    inline const result_type operator() (Index, Index = 0) const {
      return result_type(rnddist(rng), rnddist(rng));
    }
  };
} // end namespace tomo_internal

} // namespace Tools
} // namespace Tomographer

namespace Eigen {
  namespace internal {
    template<typename Rng, typename RndDist, typename Scalar>
    struct functor_traits<Tomographer::Tools::tomo_internal::random_generator<Rng, RndDist, Scalar> >
    { enum { Cost = 50 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; };
  }
} // end namespace Eigen


namespace Tomographer {
namespace Tools {


template<typename Der, typename Rng, typename RndDist, typename... IndexTypes>
inline auto denseRandom(Rng & rng, RndDist &rnddist, IndexTypes... sizes)
  -> const Eigen::CwiseNullaryOp<
    tomo_internal::random_generator<Rng, RndDist, typename Eigen::internal::traits<Der>::Scalar>,
    Der
    >
{
  typedef typename Der::Scalar Scalar;

  return Eigen::DenseBase<Der>::NullaryExpr(
      sizes..., tomo_internal::random_generator<Rng, RndDist, Scalar>(rng, rnddist)
      );
}



// ---------------------------



namespace tomo_internal {
  template<typename Scalar, typename IndexType>
  struct can_basis_vec_generator
  {
    typedef Scalar result_type;

    const IndexType k;
    const IndexType j;

    can_basis_vec_generator(IndexType k_, IndexType j_ = 0)
      : k(k_), j(j_)
    {
    }

    inline const result_type operator() (IndexType a, IndexType b = 0) const {
      return (a == k) && (b == j) ? result_type(1) : result_type(0);
    }
  };
} // namespace tomo_internal
} // namespace Tools
} // namespace Tomographer
namespace Eigen {
  namespace internal {
    template<typename Scalar, typename IndexType>
    struct functor_traits<Tomographer::Tools::tomo_internal::can_basis_vec_generator<Scalar, IndexType> >
    { enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = true }; };
    template<typename Scalar, typename IndexType>
    struct functor_has_linear_access<Tomographer::Tools::tomo_internal::can_basis_vec_generator<Scalar, IndexType> >
    { enum { ret = 0 }; };
  }
} // end namespace Eigen

namespace Tomographer {
namespace Tools {

template<typename Der, typename IndexType>
inline auto canonicalBasisVec(IndexType k, IndexType size)
  -> const Eigen::CwiseNullaryOp<
    tomo_internal::can_basis_vec_generator<typename Eigen::internal::traits<Der>::Scalar, IndexType>,
    Der
    >
{
  typedef typename Der::Scalar Scalar;

  return Eigen::DenseBase<Der>::NullaryExpr(
      size, tomo_internal::can_basis_vec_generator<Scalar, IndexType>(k)
      );
}

template<typename Der, typename IndexType>
inline auto canonicalBasisVec(IndexType k, IndexType j, IndexType rows, IndexType cols)
  -> const Eigen::CwiseNullaryOp<
    tomo_internal::can_basis_vec_generator<typename Eigen::internal::traits<Der>::Scalar, IndexType>,
    Der
    >
{
  typedef typename Der::Scalar Scalar;

  return Eigen::DenseBase<Der>::NullaryExpr(
      rows, cols, tomo_internal::can_basis_vec_generator<Scalar, IndexType>(k, j)
      );
}




// ---------------------------



namespace tomo_internal {
  template<typename Scalar>
  struct powers_of_two_generator
  {
    typedef typename Eigen::NumTraits<Scalar>::Real result_type;

    powers_of_two_generator() { }

    template<typename IndexType>
    inline const result_type operator() (IndexType a) const {
      return std::ldexp(result_type(1), (int)a);
    }

    // Don't expose an operator()(i,j) method, because otherwise Eigen might think that we
    // don't have linear access
    // https://eigen.tuxfamily.org/dox/NullaryFunctors_8h_source.html#l00147
    //
    // template<typename IndexType>
    // inline const result_type operator() (IndexType a, IndexType b) const {
    //   eigen_assert(b == 0 && "powers_of_two_generator may only be used with 1-D objects or with linear access!");
    //   (void)b; // silence unused variable warning if eigen_assert is optimized out
    //   return std::ldexp(result_type(1), a);
    // }

  };
} // namespace tomo_internal
} // namespace Tools
} // namespace Tomographer
namespace Eigen {
  namespace internal {
    template<typename Scalar>
    struct functor_traits<Tomographer::Tools::tomo_internal::powers_of_two_generator<Scalar> >
    { enum { Cost = 8 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = true }; };
  }
} // end namespace Eigen

namespace Tomographer {
namespace Tools {


template<typename Der, typename... IndexTypes>
inline auto powersOfTwo(IndexTypes... sizes)
  -> const Eigen::CwiseNullaryOp<
    tomo_internal::powers_of_two_generator<typename Eigen::internal::traits<Der>::Scalar>,
    Der
    >
{
  typedef typename Der::Scalar Scalar;

  return Eigen::DenseBase<Der>::NullaryExpr(
      sizes..., tomo_internal::powers_of_two_generator<Scalar>()
      );
}





// -----------------------------------------------------------------------------

template<int RowFactorCTime, int ColFactorCTime, typename Derived,
         typename std::enable_if<(RowFactorCTime == Eigen::Dynamic || ColFactorCTime == Eigen::Dynamic),
                                 bool>::type dummy = true>
inline auto replicated(const Eigen::DenseBase<Derived> & x, int row_factor, int col_factor)
  -> const Eigen::Replicate<Derived, Eigen::Dynamic, Eigen::Dynamic>
{
  eigen_assert(RowFactorCTime == Eigen::Dynamic || row_factor == RowFactorCTime);
  eigen_assert(ColFactorCTime == Eigen::Dynamic || col_factor == ColFactorCTime);
  return x.replicate(row_factor, col_factor);
}
template<int RowFactorCTime, int ColFactorCTime, typename Derived,
         typename std::enable_if<(RowFactorCTime != Eigen::Dynamic && ColFactorCTime != Eigen::Dynamic),
                                 bool>::type dummy2 = true>
inline auto replicated(const Eigen::DenseBase<Derived> & x, int row_factor, int col_factor)
  -> const Eigen::Replicate<Derived, RowFactorCTime, ColFactorCTime>
{
  eigen_assert(row_factor == RowFactorCTime); (void)row_factor; // "unused argument" warning
  eigen_assert(col_factor == ColFactorCTime); (void)col_factor; // "unused argument" warning
  return x.template replicate<RowFactorCTime, ColFactorCTime>();
}




} // namespace Tools
} // namespace Tomographer





// -----------------------------------------------------------------------------
// Eigen Matrix/Array and boost::serialization
// -----------------------------------------------------------------------------
namespace boost {
namespace serialization {

template<typename Archive, typename Derived>
void eigen_helper_save(Archive & a, const Derived & m, unsigned int /* version */)
{
  Eigen::Index rows = m.rows();
  Eigen::Index cols = m.cols();
  a << rows;
  a << cols;
  for (Eigen::Index i = 0; i < m.rows(); ++i) {
    for (Eigen::Index j = 0; j < m.cols(); ++j) {
      a << m(i,j);
    }
  }
}

template<typename Archive, typename Derived>
void eigen_helper_load(Archive & a, Derived & m, unsigned int /* version */)
{
  Eigen::Index rows, cols;
  a >> rows;
  a >> cols;
  m.resize(rows, cols); // legal also on fixed-size matrices
  for (Eigen::Index i = 0; i < m.rows(); ++i) {
    for (Eigen::Index j = 0; j < m.cols(); ++j) {
      a >> m(i,j);
    }
  }

}


template<typename Archive,
         typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void save(Archive & a, const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
          const unsigned int version)
{
  eigen_helper_save(a, m, version);
}
template<typename Archive,
         typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void save(Archive & a, const Eigen::Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
          const unsigned int version)
{
  eigen_helper_save(a, m, version);
}
template<typename Archive,
         typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void load(Archive & a, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
          const unsigned int version)
{
  eigen_helper_load(a, m, version);
}
template<typename Archive,
         typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void load(Archive & a, Eigen::Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
          const unsigned int version)
{
  eigen_helper_load(a, m, version);
}

template<typename Archive,
         typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
inline void serialize(Archive & a, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
                      const unsigned int version)
{
  split_free(a, m, version);
}
template<typename Archive,
         typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
inline void serialize(Archive & a, Eigen::Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
                      const unsigned int version)
{
  split_free(a, m, version);
}


} // serialization
} // boost






#endif