param_herm_x.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) 2015 ETH Zurich, Institute for Theoretical Physics, 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_DENSEDM_PARAM_HERM_X_H
#define TOMOGRAPHER_DENSEDM_PARAM_HERM_X_H



#include <cmath>
#include <complex>

#include <boost/math/constants/constants.hpp>

#include <Eigen/Core>

#include <tomographer2/tools/cxxutil.h> // static_or_dynamic, tomographer_assert()


// // heavily inspired by Eigen::Replicate [Eigen/src/Core/Replicate.h]
// namespace Eigen {
// namespace internal {
// template<typename MatrixType,int RowFactor,int ColFactor>
// struct traits<Tomographer::DenseDM::internal::ParamHermToX<MatrixType,RowFactor,ColFactor> >
//   : traits<MatrixType>
// {
//   typedef Eigen::NumTraits<typename MatrixType::Scalar>::Real Scalar;
//   typedef typename traits<MatrixType>::StorageKind StorageKind;
//   typedef typename traits<MatrixType>::XprKind XprKind;
//   enum {
//     Factor = (RowFactor==Dynamic || ColFactor==Dynamic) ? Dynamic : RowFactor*ColFactor
//   };
//   typedef typename nested<MatrixType,1>::type MatrixTypeNested;
//   typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;
//   enum {
//     RowsAtCompileTime = (int(MatrixType::RowsAtCompileTime)==Eigen::Dynamic
//                          ? Eigen::Dynamic
//                          : MatrixType::RowsAtCompileTime*MatrixType::RowsAtCompileTime ),
//     ColsAtCompileTime = 1,
//     MaxRowsAtCompileTime = RowsAtCompileTime,
//     MaxColsAtCompileTime = 1,
//     IsRowMajor = MaxRowsAtCompileTime==1 && MaxColsAtCompileTime!=1 ? 1
//                : MaxColsAtCompileTime==1 && MaxRowsAtCompileTime!=1 ? 0
//                : (MatrixType::Flags & RowMajorBit) ? 1 : 0,
//     Flags = (_MatrixTypeNested::Flags & HereditaryBits & ~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0),
//     CoeffReadCost = _MatrixTypeNested::CoeffReadCost * Eigen::NumTraits<typename MatrixType::Scalar>::MulCost
//   };
// };
// } // namespace internal
// } // namespace Eigen


// namespace Tomographer {
// namespace DenseDM {
// namespace internal {

// template<typename MatrixType>
// class ParamHermToX : public Eigen::MatrixBase<ParamHermToX<MatrixType> >
// {
//   typedef typename internal::traits<ParamHermToX>::Scalar RealScalar;
//   typedef typename internal::traits<ParamHermToX>::MatrixTypeNested MatrixTypeNested;
//   typedef typename internal::traits<ParamHermToX>::_MatrixTypeNested _MatrixTypeNested;
// public:

//   typedef typename MatrixBase<ParamHermToX> Base;
//   EIGEN_DENSE_PUBLIC_INTERFACE(ParamHermToX)
  
//   template<typename OriginalMatrixType>
//   inline explicit ParamHermToX(const OriginalMatrixType& H)
//     : _H(H), _dimtri(H.rows()*(H.rows()-1)/2)
//   {
//     EIGEN_STATIC_ASSERT((internal::is_same<typename internal::remove_const<MatrixType>::type,OriginalMatrixType>::value),
//                         THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE) ;
//     tomographer_assert(H.rows() == H.cols())
//   }
  
//   inline Index rows() const { return _H.rows()*_H.rows(); }
//   inline Index cols() const { return 1; }

//   inline RealScalar coeff(Index rowId, Index colId) const
//   {
//     tomographer_assert(colId == 0);

//     if (rowId < _H.rows()) {
//       return _H.coeff(rowId, rowId).real();
//     }
//     rowId -= _H.rows();
//     if (rowId < _dimtri.value()) {
//       return _H.coeff(......????)
//     }
//
//     ............ MAYBE NOT THE WAY TO GO .............

//     // try to avoid using modulo; this is a pure optimization strategy
//     const Index actual_row  = internal::traits<MatrixType>::RowsAtCompileTime==1 ? 0
//       : RowFactor==1 ? rowId
//       : rowId%m_matrix.rows();
//     const Index actual_col  = internal::traits<MatrixType>::ColsAtCompileTime==1 ? 0
//       : ColFactor==1 ? colId
//       : colId%m_matrix.cols();
    
//     return m_matrix.coeff(actual_row, actual_col);
//   }

//   const _MatrixTypeNested& nestedExpression() const
//   { 
//     return m_matrix;
//   }

// protected:
//   MatrixTypeNested _H;
//   const static_or_dynamic<Index,(MatrixType::RowsAtCompileTime==Eigen::Dynamic),
//                           Index(MatrixType::RowsAtCompileTime*(MatrixType::RowsAtCompileTime-1)/2> _dimtri;
// };

// } // namespace internal 
// } // namespace DenseDM
// } // namespace Tomographer


// ================================================================================


namespace Tomographer {
namespace DenseDM {

template<typename DMTypes_>
class ParamX {
public:
  typedef DMTypes_ DMTypes;
  typedef typename DMTypes::MatrixType MatrixType;
  typedef typename DMTypes::MatrixTypeConstRef MatrixTypeConstRef;
  typedef typename DMTypes::VectorParamType VectorParamType;
  typedef typename DMTypes::VectorParamTypeConstRef VectorParamTypeConstRef;
  typedef typename DMTypes::RealScalar RealScalar;
  typedef typename DMTypes::ComplexScalar ComplexScalar;
  typedef typename MatrixType::Index IndexType;

  ParamX(DMTypes dmt) : _dmt(dmt) { }

  inline VectorParamType HermToX(MatrixTypeConstRef Herm) const
  {
    // hope RVO kicks in
    VectorParamType x(_dmt.initVectorParamType());
    const IndexType dimtri = (_dmt.dim2() - _dmt.dim())/2;

    tomographer_assert((IndexType)_dmt.dim() == Herm.cols()); // assert Herm is (dim x dim)
    
    x.block(0,0,_dmt.dim(),1) = Herm.real().diagonal();

    IndexType k = _dmt.dim();
    IndexType n, m;
    for (n = 1; n < (IndexType)_dmt.dim(); ++n) {
      for (m = 0; m < n; ++m) {
        x(k)          = Herm(n,m).real() * boost::math::constants::root_two<RealScalar>();
        x(dimtri + k) = Herm(n,m).imag() * boost::math::constants::root_two<RealScalar>();
        ++k;
      }
    }

    return x;
  }
  
  template<bool OnlyLowerTri = false>
  inline MatrixType XToHerm(VectorParamTypeConstRef x) const
  {
    // should be optimized by compiler via RVO
    MatrixType Herm(_dmt.initMatrixType());
    
    const IndexType dimtri = (_dmt.dim2()-_dmt.dim())/2;
    tomographer_assert(x.rows() == (IndexType)_dmt.dim2() && x.cols() == 1); // assert x is (dim*dim x 1)

    Herm.diagonal().real() = x.block(0,0,_dmt.dim(),1);
    Herm.diagonal().imag().setZero();
  
    IndexType k = _dmt.dim();
    IndexType n, m;
    for (n = 1; n < (IndexType)_dmt.dim(); ++n) {
      for (m = 0; m < n; ++m) {
        Herm(n,m) = boost::math::constants::half_root_two<RealScalar>() * ComplexScalar(x(k), x(dimtri + k));
        if (!OnlyLowerTri) {
          // complex conj. on opposite triangular part
          Herm(m,n) = boost::math::constants::half_root_two<RealScalar>() * ComplexScalar(x(k), -x(dimtri + k));
        }
        ++k;
      }
    }
    return Herm;
  }
  
protected:
  const DMTypes _dmt;
};




} // namespace DenseDM
} // namespace Tomographer


#endif