Tomographer  v5.4
Tomographer C++ Framework Documentation
tspacellhwalker.h
Go to the documentation of this file.
1 /* This file is part of the Tomographer project, which is distributed under the
2  * terms of the MIT license.
3  *
4  * The MIT License (MIT)
5  *
6  * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
7  * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
8  *
9  * Permission is hereby granted, free of charge, to any person obtaining a copy
10  * of this software and associated documentation files (the "Software"), to deal
11  * in the Software without restriction, including without limitation the rights
12  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13  * copies of the Software, and to permit persons to whom the Software is
14  * furnished to do so, subject to the following conditions:
15  *
16  * The above copyright notice and this permission notice shall be included in
17  * all copies or substantial portions of the Software.
18  *
19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  * SOFTWARE.
26  */
27 
28 #ifndef TOMOGRAPHER_DENSEDM_TSPACELLHWALKER_H
29 #define TOMOGRAPHER_DENSEDM_TSPACELLHWALKER_H
30 
31 #include <cstddef>
32 #include <cmath>
33 
34 #include <random>
35 
36 #include <boost/math/constants/constants.hpp>
37 
43 #include <tomographer/mhrw.h>
44 
53 namespace Tomographer {
54 namespace DenseDM {
55 namespace TSpace {
56 
57 
58 
59 
60 namespace tomo_internal {
61 
62 template<typename DenseLLHType, typename = void>
63 struct DenseLLHInvoker
64 {
65  typedef typename DenseLLHType::DMTypes DMTypes;
66  typedef typename DMTypes::MatrixType MatrixType;
67  typedef typename DMTypes::VectorParamType VectorParamType;
68  typedef typename DenseLLHType::LLHValueType LLHValueType;
69 
70  const DenseLLHType & llh;
71  const Tools::StoreIfEnabled<const ParamX<DMTypes>, ((int)DenseLLHType::LLHCalcType == (int)LLHCalcTypeX)> param_x;
72 
73  DenseLLHInvoker(const DenseLLHType & llh_)
74  : llh(llh_),
75  param_x(llh.dmt)
76  {
77  }
78 
79  // This implementation deals for \a DenseLLH objects which expose a \a logLikelihoodX()
80  // function (see \ref pageInterfaceDenseLLH)
81  TOMOGRAPHER_ENABLED_IF(DenseLLHType::LLHCalcType == LLHCalcTypeX)
82  inline LLHValueType fnLogVal(const MatrixType & T) const
83  {
84  MatrixType rho(T*T.adjoint());
85  VectorParamType x = param_x.value.HermToX(rho);
86  LLHValueType llhval = llh.logLikelihoodX(x);
87  return llhval;
88  }
89 
90  // This implementation deals for \a DenseLLH objects which expose a \a logLikelihoodRho()
91  // function (see \ref pageInterfaceDenseLLH)
92  TOMOGRAPHER_ENABLED_IF(DenseLLHType::LLHCalcType == LLHCalcTypeRho)
93  inline LLHValueType fnLogVal(const MatrixType & T) const
94  {
95  LLHValueType llhval = llh.logLikelihoodRho(T*T.adjoint());
96  return llhval;
97  }
98 
99 };
100 
101 } // namespace tomo_internal
102 
103 
104 
105 
119 template<typename DenseLLHType_, typename RngType_, typename LoggerType_>
120 class TOMOGRAPHER_EXPORT LLHMHWalker
121  : public Tools::NeedOwnOperatorNew<typename DenseLLHType_::DMTypes::MatrixType>::ProviderType
122 {
123 public:
125  typedef DenseLLHType_ DenseLLHType;
127  typedef RngType_ RngType;
129  typedef LoggerType_ LoggerType;
130 
132  typedef typename DenseLLHType::DMTypes DMTypes;
134  typedef typename DenseLLHType::LLHValueType LLHValueType;
136  typedef typename DMTypes::MatrixType MatrixType;
140  typedef typename DMTypes::RealScalar RealScalar;
143 
148 
150  typedef MatrixType PointType;
152  typedef LLHValueType FnValueType;
154  enum {
158  UseFnSyntaxType = MHUseFnLogValue
159  };
160 
161 private:
162 
163  const DenseLLHType & _llh;
164  const tomo_internal::DenseLLHInvoker<DenseLLHType> _llhinvoker;
165  RngType & _rng;
166  std::normal_distribution<RealScalar> _normal_distr_rnd;
167 
168  LoggerType & _log;
169 
170  MatrixType _startpt;
171 
172 public:
173 
184  LLHMHWalker(const MatrixType & startpt, const DenseLLHType & llh, RngType & rng, LoggerType & log_)
185  : _llh(llh),
186  _llhinvoker(llh),
187  _rng(rng),
188  _normal_distr_rnd(0.0, 1.0),
189  _log(log_),
190  _startpt(startpt)
191  {
192  }
193 
194 
204  inline void init()
205  {
206  _log.debug("TSpace::LLHMHWalker", "Starting random walk");
207  }
208 
210  inline const MatrixType & startPoint()
211  {
212  // It's fine to hard-code "1e-3" because for any type, valid T-matrices have norm == 1
213  if (_startpt.norm() > 1e-3) {
214  // nonzero matrix given: that's the starting point.
215  return _startpt;
216  }
217 
218  // zero matrix given: means to choose random starting point
219  MatrixType T(_llh.dmt.initMatrixType());
220  T = Tools::denseRandom<MatrixType>(
221  _rng, _normal_distr_rnd, (Eigen::Index)_llh.dmt.dim(), (Eigen::Index)_llh.dmt.dim()
222  );
223  _startpt = T/T.norm(); // normalize to be on surface of the sphere
224 
225  _log.debug("TSpace::LLHMHWalker", [&](std::ostream & str) {
226  str << "Chosen random start point T = \n" << _startpt;
227  });
228 
229  // return start point
230  return _startpt;
231  }
232 
234  inline void thermalizingDone()
235  {
236  }
237 
239  inline void done()
240  {
241  }
242 
248  inline LLHValueType fnLogVal(const MatrixType & T) const
249  {
250  return _llhinvoker.fnLogVal(T);
251  }
252 
254  inline MatrixType jumpFn(const MatrixType& cur_T, WalkerParams params)
255  {
256  MatrixType DeltaT(Tools::denseRandom<MatrixType>(
257  _rng, _normal_distr_rnd,
258  (Eigen::Index)_llh.dmt.dim(), (Eigen::Index)_llh.dmt.dim()
259  ));
260 
261  MatrixType new_T(cur_T + params.step_size * DeltaT);
262 
263  // renormalize to "project" onto the large T-space sphere
264  new_T /= new_T.norm(); //Matrix<>.norm() is Frobenius norm.
265 
266  // _log.longdebug("jump_fn(): step_size=%g, cur_T =\n%s\nDeltaT = \n%s\nnew_T = \n%s",
267  // step_size, streamstr(cur_T).c_str(), streamstr(DeltaT).c_str(),
268  // streamstr(new_T).c_str());
269 
270  // hope for Return Value Optimization by compiler
271  return new_T;
272  }
273 
274 };
275 
276 
277 
278 // /** \brief Jump parameters (step size) for \ref LLHMHWalkerLight
279 // *
280 // */
281 // template<typename AngleRealType_, typename JumpsCountIntType_ = std::int_fast16_t>
282 // struct LLHMHWalkerLightParams
283 // {
284 // typedef AngleRealType_ AngleRealType;
285 // typedef JumpsCountIntType_ JumpsCountIntType;
286 
287 // LLHMHWalkerLightParams(AngleRealType angle_step_size_, JumpsCountIntType num_elem_jumps_)
288 // : angle_step_size(angle_step_size_),
289 // num_elem_jumps(num_elem_jumps_)
290 // {
291 // }
292 
293 // // initialize from MHWalkerParamsStepSize
294 // template<typename RealType>
295 // LLHMHWalkerLightParams(MHWalkerParamsStepSize<RealType> angle_step_size_param)
296 // : angle_step_size(step_size_param.step_size),
297 // num_elem_jumps(1)
298 // {
299 // }
300 
301 // AngleRealType angle_step_size;
302 // JumpsCountIntType num_elem_jumps;
303 
304 // };
305 
306 template<typename StepRealType>
308 
309 
310 // template<typename AngleRealType_, typename JumpsCountIntType_>
311 // inline std::ostream & operator<<(std::ostream & stream,
312 // const LLHMHWalkerLightParams<AngleRealType_,JumpsCountIntType_> & p)
313 // {
314 // return stream << "angle_step_size=" << p.angle_step_size << ", num_elem_jumps="<<p.num_elem_jumps;
315 // }
316 
317 
318 
319 
340 template<typename DenseLLHType_, typename RngType_, typename LoggerType_>
341  //typename JumpsCountIntType_ = std::int_fast16_t>
342 class TOMOGRAPHER_EXPORT LLHMHWalkerLight
343  : public Tools::NeedOwnOperatorNew<typename DenseLLHType_::DMTypes::MatrixType>::ProviderType
344 {
345 public:
347  typedef DenseLLHType_ DenseLLHType;
349  typedef RngType_ RngType;
351  typedef LoggerType_ LoggerType;
352 // //! Type for the number of jumps in a "step jump"
353 // typedef JumpsCountIntType_ JumpsCountIntType;
354 
356  typedef typename DenseLLHType::DMTypes DMTypes;
358  typedef typename DenseLLHType::LLHValueType LLHValueType;
360  typedef typename DMTypes::MatrixType MatrixType;
364  typedef typename DMTypes::RealScalar RealScalar;
367 
371  typedef LLHMHWalkerLightParams<RealScalar//,JumpsCountIntType
373 
375  typedef MatrixType PointType;
377  typedef LLHValueType FnValueType;
379  enum {
383  UseFnSyntaxType = MHUseFnLogValue
384  };
385 
386 private:
387 
388  const DenseLLHType & _llh;
389  const tomo_internal::DenseLLHInvoker<DenseLLHType> _llhinvoker;
390  RngType & _rng;
391  std::normal_distribution<RealScalar> _normal_distr_rnd;
394 
396 
397  MatrixType _startpt;
398 
399 public:
400 
411  LLHMHWalkerLight(const MatrixType & startpt, const DenseLLHType & llh, RngType & rng, LoggerType & baselogger)
412  : _llh(llh),
413  _llhinvoker(llh),
414  _rng(rng),
415  _normal_distr_rnd(0.0, 1.0),
416  _jumptype_distr_rnd(0, 2), // choice in {0, 1, 2}
417  _jumpdir_distr_rnd(0, (Eigen::Index)llh.dmt.dim2()-1), // choice in {0,1,...,dim^2-1}
418  _llogger("Tomographer::DenseDM::TSpace::LLHMHWalkerLight", baselogger),
419  _startpt(startpt)
420  {
421  }
422 
423 
433  inline void init()
434  {
435  auto logger = _llogger.subLogger(TOMO_ORIGIN) ;
436  logger.debug("Starting random walk");
437  }
438 
440  inline const MatrixType & startPoint()
441  {
442  auto logger = _llogger.subLogger(TOMO_ORIGIN) ;
443 
444  // It's fine to hard-code "1e-3" because for any type, valid T-matrices have norm == 1
445  if (_startpt.norm() > 1e-3) {
446  // nonzero matrix given: that's the starting point.
447  return _startpt;
448  }
449 
450  // zero matrix given: means to choose random starting point
451  MatrixType T(_llh.dmt.initMatrixType());
452  T = Tools::denseRandom<MatrixType>(
453  _rng, _normal_distr_rnd, (Eigen::Index)_llh.dmt.dim(), (Eigen::Index)_llh.dmt.dim()
454  );
455  _startpt = T/T.norm(); // normalize to be on surface of the sphere
456 
457  logger.debug([&](std::ostream & str) {
458  str << "Chosen random start point T = \n" << _startpt;
459  });
460 
461  // return start point
462  return _startpt;
463  }
464 
466  inline void thermalizingDone()
467  {
468  }
469 
471  inline void done()
472  {
473  }
474 
480  inline LLHValueType fnLogVal(const MatrixType & T) const
481  {
482  return _llhinvoker.fnLogVal(T);
483  }
484 
486  inline MatrixType jumpFn(const MatrixType& cur_T, WalkerParams params)
487  {
488  MatrixType new_T(cur_T);
489 
490  auto logger = _llogger.subLogger(TOMO_ORIGIN) ;
491 
492  // repeat several times, to have some chance of the effect not just rotating the purification ...
493  for (int j = 0; j < (int)_llh.dmt.dim(); ++j) {
494 
495  // select two random indices, and randomly select whether we apply an elementary x, y or z rotation
496  Eigen::Index k1 = _jumpdir_distr_rnd(_rng);
497  Eigen::Index k2;
498  do { k2 = _jumpdir_distr_rnd(_rng); } while (k1 == k2); // choose also randomly, but different than k1.
499  if (k1 > k2) { std::swap(k1, k2); } // ensure that k1 < k2
500  std::int_fast8_t xyz = _jumptype_distr_rnd(_rng);
501 
502  RealScalar sina = params.step_size * _normal_distr_rnd(_rng);
503  if (sina < -1) { sina = -1; }
504  if (sina > 1) { sina = 1; }
505  RealScalar cosa = std::sqrt(1 - sina*sina);
506  // apply transformation
508  // remember: e^{i\phi(\vec{n}\cdot\vec{\sigma})} = \cos\phi \Ident + i\sin\phi (\vec{n}\cdot\vec{\sigma})
509  switch (xyz) {
510  case 0: // X-type rotation
511  tr2d(0,0) = cosa; tr2d(0,1) = ComplexScalar(0,sina);
512  tr2d(1,0) = ComplexScalar(0,sina); tr2d(1,1) = cosa;
513  break;
514  case 1: // Y-type rotation
515  tr2d(0,0) = cosa; tr2d(0,1) = sina;
516  tr2d(1,0) = -sina; tr2d(1,1) = cosa;
517  break;
518  case 2: // Z-type rotation
519  tr2d(0,0) = ComplexScalar(cosa, sina); tr2d(0,1) = 0;
520  tr2d(1,0) = 0; tr2d(1,1) = ComplexScalar(cosa, -sina);
521  break;
522  default:
523  tomographer_assert(false && "Invalid rotation type number sampled!");
524  }
525  // apply the elementary rotation onto new_T, seen as a vector
526  const Eigen::Index i1 = k1 / (Eigen::Index)_llh.dmt.dim();
527  const Eigen::Index j1 = k1 % (Eigen::Index)_llh.dmt.dim();
528  const Eigen::Index i2 = k2 / (Eigen::Index)_llh.dmt.dim();
529  const Eigen::Index j2 = k2 % (Eigen::Index)_llh.dmt.dim();
530 
531  logger.longdebug([&](std::ostream & stream) {
532  stream << "Elementary jump rotation: "
533  << "k1="<<k1<<" -> i1="<<i1<<" j1="<<j1<<" k2="<<k2<<" -> i2="<<i2<<" j2="<<j2<<"\n"
534  << "tr2d=" << tr2d;
535  }) ;
536 
537  const auto x = tr2d(0,0) * new_T(i1,j1) + tr2d(0,1) * new_T(i2,j2) ;
538  const auto y = tr2d(1,0) * new_T(i1,j1) + tr2d(1,1) * new_T(i2,j2) ;
539  new_T(i1,j1) = x;
540  new_T(i2,j2) = y;
541  }
542 
543  // ensure continued normalization
544  new_T /= new_T.norm(); // norm is Frobenius norm
545 
546  // new_T is ready, return it
547  return new_T;
548  }
549 
550 };
551 
552 
553 
554 
555 
556 
557 } // namespace TSpace
558 } // namespace DenseDM
559 } // namespace Tomographer
560 
561 
562 #endif
MatrixType jumpFn(const MatrixType &cur_T, WalkerParams params)
Decides of a new point to jump to for the random walk.
MatrixType PointType
Provided for MHRandomWalk. A point in our random walk = a density matrix.
Tomographer::DenseDM::DMTypes< Eigen::Dynamic, RealScalar > DMTypes
The Tomographer::DenseDM::DMTypes we should use by default, with dynamic sized matrices.
Definition: pydensedm.h:43
const MatrixType & startPoint()
Return the starting point given in the constructor, or a random start point.
Eigen::Matrix< RealScalar, FixedDim2, 1 > VectorParamType
Real vector with dim*dim elements.
Definition: dmtypes.h:139
MHWalkerParamsStepSize< RealScalar > WalkerParams
The real scalar type for the step size — this scalar type is the only parameter we need (see MHWalke...
Base namespace for the Tomographer project.
Definition: densellh.h:45
MatrixType PointType
Provided for MHRandomWalk. A point in our random walk = a density matrix.
LoggerType_ LoggerType
The logger type.
T swap(T... args)
C++ types for describing dense density matrices in various parameterizations.
DMTypes::RealScalar RealScalar
The real scalar corresponding to our data types. Usually a double.
RngType_ RngType
The random number generator type.
DenseLLHType_ DenseLLHType
The DenseLLH interface object type.
DMTypes::VectorParamType VectorParamType
Type of an X-parameterization of a density operator (see X Parameterization)
An MHWalkerParams type which just stores a step size.
Definition: mhrw.h:97
void init()
Provided for MHRandomWalk. Initializes some fields and prepares for a random walk.
MatrixType jumpFn(const MatrixType &cur_T, WalkerParams params)
Decides of a new point to jump to for the random walk.
DMTypes::MatrixType MatrixType
The matrix type for a density operator on our quantum system.
Provide appropriate operator new() definitions for a structure which has a member of the given stored...
STL class.
LLHMHWalkerLightParams< RealScalar > WalkerParams
The real scalar type for the step size — this scalar type is the only parameter we need (see MHWalke...
void done()
Callback for after random walk is finished. No-op.
#define TOMO_ORIGIN
Use this as argument for a Tomographer::Logger::LocalLogger constructor .
Definition: loggers.h:1608
LocalLogger< LocalLogger< BaseLoggerType > > subLogger(std::string new_prefix)
Create a sub-logger.
Definition: loggers.h:1846
DMTypes::MatrixType MatrixType
The matrix type for a density operator on our quantum system.
DMTypes::ComplexScalar ComplexScalar
The complex real scalar corresponding to our data types. Usually a std::complex<double>.
DenseLLHType::LLHValueType LLHValueType
The loglikelihood function value type (see DenseLLH Interface e.g. IndepMeasLLH)
DenseLLHType::DMTypes DMTypes
The data types of our problem.
LLHMHWalker(const MatrixType &startpt, const DenseLLHType &llh, RngType &rng, LoggerType &log_)
Constructor which just initializes the given fields.
Provides the logarithm MH function value at each point (see Role of UseFnSyntaxType) ...
Definition: mhrw.h:72
void init()
Provided for MHRandomWalk. Initializes some fields and prepares for a random walk.
Tools for parameterizing hermitian matrices with the X Parameterization.
Some declarations for the DenseLLH type interface.
LLHValueType fnLogVal(const MatrixType &T) const
Calculate the logarithm of the Metropolis-Hastings function value.
A random walk in the density matrix space of a Hilbert state space of a quantum system (uniform rando...
Managing the need for specific overrides to operator new() for some types (especially Eigen types) ...
RngType_ RngType
The random number generator type.
LLHMHWalkerLight(const MatrixType &startpt, const DenseLLHType &llh, RngType &rng, LoggerType &baselogger)
Constructor which just initializes the given fields.
The DenseLLH-compatible object exposes a method logLikelihoodRho(), taking as argument a (const ref t...
Definition: densellh.h:61
Eigen::Matrix< ComplexScalar, FixedDim, FixedDim, Eigen::Matrix< ComplexScalar, FixedDim, FixedDim >::Options, MaxFixedDim_, MaxFixedDim_ > MatrixType
Matrix type, to store the density operator as a dense matrix.
Definition: dmtypes.h:134
A random walk in the density matrix space of a Hilbert state space of a quantum system (elementary ro...
LLHValueType FnValueType
Provided for MHRandomWalk. The function value type is the loglikelihood value type.
DenseLLHType::LLHValueType LLHValueType
The loglikelihood function value type (see DenseLLH Interface e.g. IndepMeasLLH)
void done()
Callback for after random walk is finished. No-op.
DMTypes::ComplexScalar ComplexScalar
The complex real scalar corresponding to our data types. Usually a std::complex<double>.
DenseLLHType::DMTypes DMTypes
The data types of our problem.
DenseLLHType_ DenseLLHType
The DenseLLH interface object type.
void thermalizingDone()
Callback for after thermalizing is done. No-op.
DMTypes::RealScalar RealScalar
The real scalar corresponding to our data types. Usually a double.
T sqrt(T... args)
Routines for performing a Metropolis-Hastings random walk.
const MatrixType & startPoint()
Return the starting point given in the constructor, or a random start point.
LLHValueType fnLogVal(const MatrixType &T) const
Calculate the logarithm of the Metropolis-Hastings function value.
STL class.
#define tomographer_assert(...)
Assertion test macro.
Definition: cxxdefs.h:84
The DenseLLH-compatible object exposes a method logLikelihoodX(), taking as argument a (const ref to ...
Definition: densellh.h:66
void thermalizingDone()
Callback for after thermalizing is done. No-op.
Utilities for logging messages.
LLHValueType FnValueType
Provided for MHRandomWalk. The function value type is the loglikelihood value type.
RealScalar_ RealScalar
Real scalar type, given in template parameter. Usually double is fine.
Definition: dmtypes.h:121
DMTypes::VectorParamType VectorParamType
Type of an X-parameterization of a density operator (see X Parameterization)