dmmhrw.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 DMLLHINTEGRATOR_H
#define DMLLHINTEGRATOR_H

#include <cstddef>
#include <cmath>

#include <random>

#include <tomographer/qit/util.h>
#include <tomographer/qit/dist.h>
#include <tomographer/tomoproblem.h>
#include <tomographer/tools/loggers.h>
#include <tomographer/mhrw.h>

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

namespace Tomographer {


template<typename TomoProblem, typename Rng, typename Log>
class DMStateSpaceLLHMHWalker
{
public:
  typedef typename TomoProblem::MatrQ MatrQ;
  typedef typename TomoProblem::LLHValueType LLHValueType;
  typedef typename MatrQ::MatrixType MatrixType;
  typedef typename MatrQ::VectorParamType VectorParamType;
  typedef typename MatrQ::RealScalar RealScalar;
  typedef typename MatrQ::ComplexScalar ComplexScalar;

  typedef MatrixType PointType;
  typedef LLHValueType FnValueType;
  enum {
    UseFnSyntaxType = MHUseFnLogValue
  };

private:

  const TomoProblem & _tomo;
  Rng & _rng;
  std::normal_distribution<RealScalar> _normal_distr_rnd;

  Log & _log;
  
  MatrixType _startpt;

public:

  DMStateSpaceLLHMHWalker(const MatrixType & startpt, const TomoProblem & tomo, Rng & rng,
                          Log & log_)
    : _tomo(tomo),
      _rng(rng),
      _normal_distr_rnd(0.0, 1.0),
      _log(log_),
      _startpt(startpt)
  {
  }

  inline void init()
  {
    _log.debug("DMStateSpaceLLHMHWalker", "Starting random walk");
  }

  inline const MatrixType & startpoint()
  {
    if (_startpt.norm() > 1e-3) {
      // nonzero matrix given: that's the starting point
      return _startpt;
    }

    // zero matrix given: means to choose random starting point
    MatrixType T = _tomo.matq.initMatrixType();
    T = dense_random<MatrixType>(
        _rng, _normal_distr_rnd, _tomo.matq.dim(), _tomo.matq.dim()
        );
    _startpt = T/T.norm(); // normalize to be on surface of the sphere

    _log.debug("DMStateSpaceLLHMHWalker", [&T](std::ostream & str) {
        str << "Chosen random start point T = \n" << T;
      });

    // return start point
    return _startpt;
  }

  inline void thermalizing_done()
  {
  }

  inline void done()
  {
  }

  inline LLHValueType fnlogval(const MatrixType & T)
  {
    MatrixType rho = _tomo.matq.initMatrixType();
    rho.noalias() = T*T.adjoint();

    VectorParamType x = _tomo.matq.initVectorParamType();
    param_herm_to_x(x, rho);

    LLHValueType llhval = -boost::math::constants::half<RealScalar>() * _tomo.calc_llh(x);
    // _log.longdebug("fnlogval(%s) = %g\n", streamcstr(x.transpose()), llhval);
    return llhval;
  }

  inline MatrixType jump_fn(const MatrixType& cur_T, RealScalar step_size)
  {
    MatrixType DeltaT = dense_random<MatrixType>(
        _rng, _normal_distr_rnd, _tomo.matq.dim(), _tomo.matq.dim()
        );

    MatrixType new_T = _tomo.matq.initMatrixType();

    new_T.noalias() = cur_T + step_size * DeltaT;

    // renormalize to "project" onto the large T-space sphere
    new_T /= new_T.norm(); //Matrix<>.norm() is Frobenius norm.

    //    _log.longdebug("jump_fn(): step_size=%g, cur_T =\n%s\nDeltaT = \n%s\nnew_T = \n%s",
    //                   step_size, streamstr(cur_T).c_str(), streamstr(DeltaT).c_str(),
    //                   streamstr(new_T).c_str());

    // hope for Return Value Optimization by compiler
    return new_T;
  }

};





template<typename TomoProblem_, typename FidValueType_ = double>
class FidelityToRefCalculator
{
public:
  typedef TomoProblem_ TomoProblem;
  typedef typename TomoProblem::MatrQ MatrQ;
  typedef typename MatrQ::MatrixType MatrixType;

  typedef FidValueType_ ValueType;

private:
  MatrixType _ref_T;

public:
  FidelityToRefCalculator(const TomoProblem & tomo)
    : _ref_T(tomo.matq.initMatrixType())
  {
    _ref_T = tomo.T_MLE;
  }
  FidelityToRefCalculator(const TomoProblem & tomo, const MatrixType & T_ref)
    : _ref_T(tomo.matq.initMatrixType())
  {
    _ref_T = T_ref;
  }

  inline ValueType getValue(const MatrixType & T) const
  {
    return fidelity_T<ValueType>(T, _ref_T);
  }
};

template<typename TomoProblem_, typename ValueType_ = double>
class PurifDistToRefCalculator
{
public:
  typedef TomoProblem_ TomoProblem;
  typedef typename TomoProblem::MatrQ MatrQ;
  typedef typename MatrQ::MatrixType MatrixType;

  typedef ValueType_ ValueType;

private:
  MatrixType _ref_T;

public:
  PurifDistToRefCalculator(const TomoProblem & tomo)
    : _ref_T(tomo.matq.initMatrixType())
  {
    _ref_T = tomo.T_MLE;
  }
  PurifDistToRefCalculator(const TomoProblem & tomo, const MatrixType & T_ref)
    : _ref_T(tomo.matq.initMatrixType())
  {
    _ref_T = T_ref;
  }

  inline ValueType getValue(const MatrixType & T) const
  {
    auto F = fidelity_T<ValueType>(T, _ref_T);
    return std::sqrt(ValueType(1) - F*F);
  }
};

template<typename TomoProblem_, typename TrDistValueType_ = double>
class TrDistToRefCalculator
{
public:
  typedef TomoProblem_ TomoProblem;
  typedef typename TomoProblem::MatrQ MatrQ;
  typedef typename MatrQ::MatrixType MatrixType;

  typedef TrDistValueType_ ValueType;

private:
  MatrixType _ref_rho;

public:
  TrDistToRefCalculator(const TomoProblem & tomo)
    : _ref_rho(tomo.matq.initMatrixType())
  {
    _ref_rho = tomo.rho_MLE;
  }
  TrDistToRefCalculator(const TomoProblem & tomo, const MatrixType& rho_ref)
    : _ref_rho(tomo.matq.initMatrixType())
  {
    _ref_rho = rho_ref;
  }

  inline ValueType getValue(const MatrixType & T) const
  {
    return boost::math::constants::half<ValueType>() *
      (T*T.adjoint() - _ref_rho).template selfadjointView<Eigen::Lower>()
      .eigenvalues().cwiseAbs().sum();
  }
};

template<typename TomoProblem_>
class ObservableValueCalculator
{
public:
  typedef TomoProblem_ TomoProblem;
  typedef typename TomoProblem::MatrQ MatrQ;
  typedef typename MatrQ::MatrixType MatrixType;
  typedef typename MatrQ::VectorParamType VectorParamType;

  typedef typename MatrQ::RealScalar ValueType;

private:
  const TomoProblem & _tomo;

  VectorParamType _A_x;

  //  MatrixType _A; // not needed

public:
  ObservableValueCalculator(const TomoProblem & tomo, const MatrixType & A)
    : _tomo(tomo),
      _A_x(tomo.matq.initVectorParamType())
  {
    param_herm_to_x(_A_x, A);
  }

  ObservableValueCalculator(const TomoProblem & tomo, const VectorParamType & A_x)
    : _tomo(tomo),
      _A_x(tomo.matq.initVectorParamType())
  {
    _A_x = A_x;
  }

  inline ValueType getValue(const MatrixType & T) const
  {
    MatrixType rho = _tomo.matq.initMatrixType();
    rho = T*T.adjoint();
    VectorParamType x = _tomo.matq.initVectorParamType();
    param_herm_to_x(x, rho);
    return _A_x.transpose() * x;
  }
};



} // namespace Tomographer


#endif