tspacellhwalker.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
 *
 * 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_TSPACELLHWALKER_H
#define TOMOGRAPHER_DENSEDM_TSPACELLHWALKER_H

#include <cstddef>
#include <cmath>

#include <random>

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

#include <tomographer2/tools/loggers.h>
#include <tomographer2/tools/needownoperatornew.h>
#include <tomographer2/densedm/densellh.h>
#include <tomographer2/densedm/dmtypes.h>
#include <tomographer2/densedm/param_herm_x.h>
#include <tomographer2/mhrw.h>

namespace Tomographer {
namespace DenseDM {
namespace TSpace {

template<typename DenseLLH_, typename Rng, typename LoggerType>
class LLHMHWalker
  : public Tools::NeedOwnOperatorNew<typename DenseLLH_::DMTypes::MatrixType>::ProviderType
{
public:
  typedef DenseLLH_ DenseLLH;
  typedef typename DenseLLH::DMTypes DMTypes;
  typedef typename DenseLLH::LLHValueType LLHValueType;
  typedef typename DMTypes::MatrixType MatrixType;
  typedef typename DMTypes::VectorParamType VectorParamType;
  typedef typename DMTypes::RealScalar RealScalar;
  typedef typename DMTypes::ComplexScalar ComplexScalar;

  typedef RealScalar StepRealType;

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

private:

  const DenseLLH & _llh;
  Rng & _rng;
  std::normal_distribution<RealScalar> _normal_distr_rnd;
  const ParamX<DMTypes> _px;

  LoggerType & _log;
  
  MatrixType _startpt;

public:

  LLHMHWalker(const MatrixType & startpt, const DenseLLH & llh, Rng & rng,
              LoggerType & log_)
    : _llh(llh),
      _rng(rng),
      _normal_distr_rnd(0.0, 1.0),
      _px(llh.dmt),
      _log(log_),
      _startpt(startpt)
  {
  }

  LLHMHWalker(LLHMHWalker && moveother)
    : _llh(std::move(moveother._llh)),
      _rng(moveother._rng),
      _normal_distr_rnd(moveother._normal_distr_rnd),
      _px(std::move(moveother._px)),
      _log(moveother._log),
      _startpt(moveother._startpt)
  {
  }


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

  inline const MatrixType & startPoint()
  {
    // It's fine to hard-code "1e-3" because for any type, valid T-matrices have norm == 1
    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(_llh.dmt.initMatrixType());
    T = Tools::denseRandom<MatrixType>(
        _rng, _normal_distr_rnd, _llh.dmt.dim(), _llh.dmt.dim()
        );
    _startpt = T/T.norm(); // normalize to be on surface of the sphere

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

    // return start point
    return _startpt;
  }

  inline void thermalizingDone()
  {
  }

  inline void done()
  {
  }

  TOMOGRAPHER_ENABLED_IF(DenseLLH::LLHCalcType == LLHCalcTypeX)
  inline LLHValueType fnLogVal(const MatrixType & T) const
  {
    MatrixType rho(T*T.adjoint());
    VectorParamType x = _px.HermToX(rho);
    LLHValueType llhval =  _llh.logLikelihoodX(x);
    // _log.longdebug("fnlogval(%s) = %g\n", streamcstr(x.transpose()), llhval);
    return llhval;
  }

  TOMOGRAPHER_ENABLED_IF(DenseLLH::LLHCalcType == LLHCalcTypeRho)
  inline LLHValueType fnLogVal(const MatrixType & T) const
  {
    LLHValueType llhval =  _llh.logLikelihoodRho(T*T.adjoint());
    // _log.longdebug("fnlogval( %s ) = %g\n", streamcstr(T*T.adjoint()), llhval);
    return llhval;
  }


  inline MatrixType jumpFn(const MatrixType& cur_T, RealScalar step_size)
  {
    MatrixType DeltaT(Tools::denseRandom<MatrixType>(
                          _rng, _normal_distr_rnd, _llh.dmt.dim(), _llh.dmt.dim()
                          ));

    MatrixType new_T(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;
  }

};



} // namespace TSpace
} // namespace DenseDM
} // namespace Tomographer


#endif