indepmeasllh.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_INDEPMEASLLH_H
#define TOMOGRAPHER_DENSEDM_INDEPMEASLLH_H

#include <cstddef>
#include <string>
#include <iomanip> // std::setprecision, std::setw and friends.

#include <Eigen/Eigen>

#include <tomographer2/tools/cxxutil.h> // StaticOrDynamic, TOMOGRAPHER_ENABLED_IF
#include <tomographer2/tools/needownoperatornew.h>
#include <tomographer2/tools/fmt.h> // streamstr
#include <tomographer2/densedm/dmtypes.h>
#include <tomographer2/densedm/param_herm_x.h>
#include <tomographer2/densedm/densellh.h>

namespace Tomographer {
namespace DenseDM {




template<typename DMTypes_, typename LLHValueType_ = typename DMTypes_::RealScalar,
         typename IntFreqType_ = int, int FixedMaxParamList_ = Eigen::Dynamic,
         bool UseNMeasAmplifyFactor_ = false>
class IndepMeasLLH
  : public Tools::NeedEigenAlignedOperatorNew::ProviderType
{
public:
  typedef DMTypes_ DMTypes;
  typedef LLHValueType_ LLHValueType;
  typedef IntFreqType_ IntFreqType;
  static constexpr int FixedMaxParamList = FixedMaxParamList_;
  static constexpr bool IsDynamicMaxParamList = (FixedMaxParamList_ == Eigen::Dynamic);
  static constexpr bool UseNMeasAmplifyFactor = UseNMeasAmplifyFactor_;

  enum {
    LLHCalcType = LLHCalcTypeX
  } ;

  typedef Eigen::Matrix<typename DMTypes::RealScalar, Eigen::Dynamic, DMTypes::FixedDim2,
                        Eigen::RowMajor, FixedMaxParamList, DMTypes::FixedDim2>  VectorParamListType;
  typedef const Eigen::Ref<const VectorParamListType> & VectorParamListTypeConstRef;

  typedef typename VectorParamListType::Index IndexType;

  typedef Eigen::Array<IntFreqType, Eigen::Dynamic, 1,
                       0 /*Options*/, FixedMaxParamList, 1>  FreqListType;
  typedef const Eigen::Ref<const FreqListType> & FreqListTypeConstRef;

  
  inline IndepMeasLLH(DMTypes dmt_)
    : dmt(dmt_), _Exn(VectorParamListType::Zero(0, dmt.dim2())), _Nx(FreqListType::Zero(0)),
      _NMeasAmplifyFactor(1)
  {
  }

  inline IndepMeasLLH(DMTypes dmt_, VectorParamListTypeConstRef Exn_, FreqListTypeConstRef Nx_)
    : dmt(dmt_), _Exn(VectorParamListType::Zero(0, dmt.dim2())), _Nx(FreqListType::Zero(0)), _NMeasAmplifyFactor(1)
  {
    setMeas(Exn_, Nx_);
  }

  const DMTypes dmt;

  inline const IndexType numEffects() const { return _Exn.rows(); }

  inline const VectorParamListType & Exn() const { return _Exn; }

  inline Eigen::Ref<const typename DMTypes::VectorParamType>  Exn(IndexType i) const
  {
    return _Exn.row(i).transpose();
  }


  inline const FreqListType & Nx() const { return _Nx; }

  inline const IntFreqType Nx(IndexType i) const { return _Nx(i); }

  inline void resetMeas()
  {
    _Exn.resize(0, dmt.dim2());
    _Nx.resize(0);
  }

  inline void addMeasEffect(typename DMTypes::VectorParamTypeConstRef E_x, IntFreqType n,
                            bool check_validity = true)
  {
    // no change if frequency count for the added effect == 0
    if (n == 0) {
      return;
    }

    tomographer_assert(n > 0);

    if (check_validity) { // check validity of measurement data
      typename DMTypes::MatrixType E_m(dmt.initMatrixType());
      E_m = ParamX<DMTypes>(dmt).XToHerm(E_x);
      _check_effect(E_m);
    }

    IndexType newi = _Exn.rows();
    tomographer_assert(newi == _Nx.rows());

    _Exn.conservativeResize(newi + 1, Eigen::NoChange);
    _Exn.row(newi) = E_x.transpose();
    _Nx.conservativeResize(newi + 1, Eigen::NoChange);
    _Nx(newi) = n;

    tomographer_assert(_Exn.rows() == _Nx.rows());
  }

  inline void addMeasEffect(typename DMTypes::MatrixTypeConstRef E_m, IntFreqType n,
                            bool check_validity = true)
  {
    tomographer_assert(E_m.rows() == E_m.cols());
    tomographer_assert(E_m.rows() == (IndexType)dmt.dim());

    // skip if no frequency count
    if (n == 0) {
      return;
    }

    if (check_validity) { // check validity of measurement data
      tomographer_assert(n > 0);
      _check_effect(E_m);
    }

    addMeasEffect(Tomographer::DenseDM::ParamX<DMTypes>(dmt).HermToX(E_m), n, false);
  }

  inline void setMeas(VectorParamListTypeConstRef Exn_, FreqListTypeConstRef Nx_, bool check_validity = true)
  {
    tomographer_assert(Exn_.cols() == (IndexType)dmt.dim2());
    tomographer_assert(Exn_.rows() == Nx_.rows());
    tomographer_assert(Nx_.cols() == 1);

    if ((Nx_ > 0).all()) {
      // all measurements are OK, so we can just copy the data.
      _Exn.resize(Exn_.rows(), dmt.dim2());
      _Exn = Exn_;
      _Nx.resize(Nx_.rows(), 1);
      _Nx = Nx_;
    } else {
      // otherwise, we have to add all measurements one by one, and filter out the zero-count ones.
      resetMeas();
      for (IndexType i = 0; i < Exn_.rows(); ++i) {
        addMeasEffect(Exn_.row(i), Nx_(i));
      }
    }
    if (check_validity) {
      checkAllMeas();
    }
  }

  inline void checkAllMeas() const
  {
    tomographer_assert(_Exn.cols() == (IndexType)dmt.dim2());
    tomographer_assert(_Exn.rows() == _Nx.rows());
    tomographer_assert(_Nx.cols() == 1);

    for (IndexType i = 0; i < _Exn.rows(); ++i) {

      tomographer_assert(_Nx(i) > 0);

      typename DMTypes::MatrixType E_m(dmt.initMatrixType());
      E_m = ParamX<DMTypes>(dmt).XToHerm(_Exn.row(i).transpose());
      _check_effect(E_m);
    }
  }
  inline void checkEffect(IndexType i) const
  {
    tomographer_assert(_Exn.cols() == (IndexType)dmt.dim2());
    tomographer_assert(_Exn.rows() == _Nx.rows());
    tomographer_assert(_Nx.cols() == 1);
    tomographer_assert(i >= 0 && i < _Exn.rows());
    tomographer_assert(_Nx(i) > 0);

    typename DMTypes::MatrixType E_m(dmt.initMatrixType());
    E_m = ParamX<DMTypes>(dmt).XToHerm(_Exn.row(i).transpose());
    _check_effect(E_m);
  }

private:
  inline void _check_effect(typename DMTypes::MatrixType E_m) const
  {
    if ( ! (double( (E_m - E_m.adjoint()).norm() ) < 1e-8) ) { // matrix not Hermitian
      throw InvalidMeasData(streamstr("POVM effect is not hermitian : E_m =\n"
                                      << std::setprecision(10) << E_m));
    }
    Eigen::SelfAdjointEigenSolver<typename DMTypes::MatrixType> slv(E_m);
    const double mineigval = slv.eigenvalues().minCoeff();
    if ( ! (mineigval >= -1e-12) ) { // not positive semidef
      throw InvalidMeasData(streamstr("POVM effect is not positive semidefinite (min eigval="<<mineigval<<") : E_m =\n"
                                      << std::setprecision(10) << E_m));
    }
    if ( ! (double(E_m.norm()) > 1e-6) ) { // POVM effect is zero
      throw InvalidMeasData(streamstr("POVM effect is zero : E_m =\n" << E_m));
    }
  }

public:

  TOMOGRAPHER_ENABLED_IF(UseNMeasAmplifyFactor)
  inline LLHValueType NMeasAmplifyFactor() const { return _NMeasAmplifyFactor.value; }

  TOMOGRAPHER_ENABLED_IF(!UseNMeasAmplifyFactor)
  inline LLHValueType NMeasAmplifyFactor() const { return LLHValueType(1); }

  TOMOGRAPHER_ENABLED_IF(UseNMeasAmplifyFactor)
  inline void setNMeasAmplifyFactor(LLHValueType val) { _NMeasAmplifyFactor.value = val; }
  

  inline LLHValueType logLikelihoodX(typename DMTypes::VectorParamTypeConstRef x) const
  {
    return _mult_by_nmeasfactor(
        (_Nx.template cast<LLHValueType>() * (_Exn * x).template cast<LLHValueType>().array().log()).sum()
        );
  }

private:
  template<typename Expr, TOMOGRAPHER_ENABLED_IF_TMPL(UseNMeasAmplifyFactor)>
  inline auto _mult_by_nmeasfactor(Expr&& expr) const -> decltype(LLHValueType(1) * expr)
  {
    return _NMeasAmplifyFactor.value * expr;
  }
  template<typename Expr, TOMOGRAPHER_ENABLED_IF_TMPL(!UseNMeasAmplifyFactor)>
  inline auto _mult_by_nmeasfactor(Expr&& expr) const -> decltype(std::forward<Expr>(expr))
  {
    return std::forward<Expr>(expr);
  }

private:
  VectorParamListType _Exn;
  FreqListType _Nx;

  Tomographer::Tools::StoreIfEnabled<LLHValueType, UseNMeasAmplifyFactor> _NMeasAmplifyFactor;
};
// define static members:
template<typename DMTypes_, typename LLHValueType_,
         typename IntFreqType_, int FixedMaxParamList_,
         bool UseNMeasAmplifyFactor_>
constexpr int
IndepMeasLLH<DMTypes_,LLHValueType_,IntFreqType_,FixedMaxParamList_,UseNMeasAmplifyFactor_>::FixedMaxParamList;
template<typename DMTypes_, typename LLHValueType_,
         typename IntFreqType_, int FixedMaxParamList_,
         bool UseNMeasAmplifyFactor_>
constexpr bool
IndepMeasLLH<DMTypes_,LLHValueType_,IntFreqType_,FixedMaxParamList_,UseNMeasAmplifyFactor_>::IsDynamicMaxParamList;
template<typename DMTypes_, typename LLHValueType_,
         typename IntFreqType_, int FixedMaxParamList_,
         bool UseNMeasAmplifyFactor_>
constexpr bool
IndepMeasLLH<DMTypes_,LLHValueType_,IntFreqType_,FixedMaxParamList_,UseNMeasAmplifyFactor_>::UseNMeasAmplifyFactor;





} // namespace DenseDM
} // namespace Tomographer








#endif