queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
QUESO::Algorithm< V, M > Class Template Reference

#include <Algorithm.h>

Public Member Functions

 Algorithm (const BaseEnvironment &env, const BaseTKGroup< V, M > &tk)
 
 ~Algorithm ()
 
double acceptance_ratio (MarkovChainPositionData< V > x, MarkovChainPositionData< V > y, const V &tk_pos_x, const V &tk_pos_y)
 Calculates the finite dimensional Metropolis-Hastings acceptance ratio. More...
 

Private Attributes

const BaseEnvironmentm_env
 
const BaseTKGroup< V, M > & m_tk
 

Detailed Description

template<class V = GslVector, class M = GslMatrix>
class QUESO::Algorithm< V, M >

Definition at line 35 of file Algorithm.h.

Constructor & Destructor Documentation

template<class V , class M >
QUESO::Algorithm< V, M >::Algorithm ( const BaseEnvironment env,
const BaseTKGroup< V, M > &  tk 
)

Definition at line 36 of file Algorithm.C.

38  :
39  m_env(env),
40  m_tk(tk)
41 {
42 }
const BaseEnvironment & m_env
Definition: Algorithm.h:53
const BaseTKGroup< V, M > & m_tk
Definition: Algorithm.h:54
template<class V , class M >
QUESO::Algorithm< V, M >::~Algorithm ( )

Definition at line 45 of file Algorithm.C.

46 {
47 }

Member Function Documentation

template<class V , class M >
double QUESO::Algorithm< V, M >::acceptance_ratio ( MarkovChainPositionData< V >  x,
MarkovChainPositionData< V >  y,
const V &  tk_pos_x,
const V &  tk_pos_y 
)

Calculates the finite dimensional Metropolis-Hastings acceptance ratio.

tk_pos_x is the position of the tk when evaluating for x tk_pos_y is the position of the tk when evaluating for y

This method is called by the delayed rejection procedure.

Definition at line 51 of file Algorithm.C.

References QUESO::MarkovChainPositionData< V >::logTarget(), QUESO::MarkovChainPositionData< V >::outOfTargetSupport(), QUESO::queso_isnan(), and QUESO::MarkovChainPositionData< V >::vecValues().

56 {
57  double alphaQuotient = 0.;
58  if ((x.outOfTargetSupport() == false) &&
59  (y.outOfTargetSupport() == false)) {
60  if ((x.logTarget() == -INFINITY) ||
61  (x.logTarget() == INFINITY) ||
62  ( queso_isnan(x.logTarget()) )) {
63  std::cerr << "WARNING In Algorithm<V,M>::alpha(x,y)"
64  << ", worldRank " << m_env.worldRank()
65  << ", fullRank " << m_env.fullRank()
66  << ", subEnvironment " << m_env.subId()
67  << ", subRank " << m_env.subRank()
68  << ", inter0Rank " << m_env.inter0Rank()
69  << ": x.logTarget() = " << x.logTarget()
70  << ", x.values() = " << x.vecValues()
71  << ", y.values() = " << y.vecValues()
72  << std::endl;
73  }
74  else if ((y.logTarget() == -INFINITY ) ||
75  (y.logTarget() == INFINITY ) ||
76  ( queso_isnan(y.logTarget()) )) {
77  std::cerr << "WARNING In Algorithm<V,M>::alpha(x,y)"
78  << ", worldRank " << m_env.worldRank()
79  << ", fullRank " << m_env.fullRank()
80  << ", subEnvironment " << m_env.subId()
81  << ", subRank " << m_env.subRank()
82  << ", inter0Rank " << m_env.inter0Rank()
83  << ": y.logTarget() = " << y.logTarget()
84  << ", x.values() = " << x.vecValues()
85  << ", y.values() = " << y.vecValues()
86  << std::endl;
87  }
88  else {
89  double yLogTargetToUse = y.logTarget();
90 
91  if (m_tk.symmetric()) {
92  alphaQuotient = std::exp(yLogTargetToUse - x.logTarget());
93 
94  if ((m_env.subDisplayFile() ) &&
95  (m_env.displayVerbosity() >= 3 )) {
96  *m_env.subDisplayFile() << "In Algorithm<V,M>::alpha(x,y)"
97  << ": symmetric proposal case"
98  << ", x = " << x.vecValues()
99  << ", y = " << y.vecValues()
100  << ", yLogTargetToUse = " << yLogTargetToUse
101  << ", x.logTarget() = " << x.logTarget()
102  << ", alpha = " << alphaQuotient
103  << std::endl;
104  }
105  }
106  else {
107  double qyx = m_tk.rv(tk_pos_x).pdf().lnValue(x.vecValues());
108  if ((m_env.subDisplayFile() ) &&
109  (m_env.displayVerbosity() >= 10 )) {
110  *m_env.subDisplayFile() << m_tk.rv(tk_pos_x).pdf() << std::endl;
111  }
112  double qxy = m_tk.rv(tk_pos_y).pdf().lnValue(y.vecValues());
113  if ((m_env.subDisplayFile() ) &&
114  (m_env.displayVerbosity() >= 10 )) {
115  *m_env.subDisplayFile() << m_tk.rv(tk_pos_y).pdf() << std::endl;
116  }
117  alphaQuotient = std::exp(yLogTargetToUse +
118  qyx -
119  x.logTarget() -
120  qxy);
121  if ((m_env.subDisplayFile() ) &&
122  (m_env.displayVerbosity() >= 3 )) {
123  *m_env.subDisplayFile() << "In Algorithm<V,M>::alpha(x,y)"
124  << ": asymmetric proposal case"
125  << ", x = " << x.vecValues()
126  << ", y = " << y.vecValues()
127  << ", yLogTargetToUse = " << yLogTargetToUse
128  << ", q(y,x) = " << qyx
129  << ", x.logTarget() = " << x.logTarget()
130  << ", q(x,y) = " << qxy
131  << ", alpha = " << alphaQuotient
132  << std::endl;
133  }
134  }
135  } // protection logic against logTarget values
136  }
137  else {
138  if ((m_env.subDisplayFile() ) &&
139  (m_env.displayVerbosity() >= 10 )) {
140  *m_env.subDisplayFile() << "In Algorithm<V,M>::alpha(x,y)"
141  << ": x.outOfTargetSupport = " << x.outOfTargetSupport()
142  << ", y.outOfTargetSupport = " << y.outOfTargetSupport()
143  << std::endl;
144  }
145  }
146 
147  return std::min(1.,alphaQuotient);
148 }
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:342
const BaseEnvironment & m_env
Definition: Algorithm.h:53
const BaseTKGroup< V, M > & m_tk
Definition: Algorithm.h:54
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
bool queso_isnan(T arg)
Definition: math_macros.h:39
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
int fullRank() const
Returns the rank of the MPI process in QUESO&#39;s full communicator.
Definition: Environment.C:268
unsigned int displayVerbosity() const
Definition: Environment.C:450
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

Member Data Documentation

template<class V = GslVector, class M = GslMatrix>
const BaseEnvironment& QUESO::Algorithm< V, M >::m_env
private

Definition at line 53 of file Algorithm.h.

template<class V = GslVector, class M = GslMatrix>
const BaseTKGroup<V, M>& QUESO::Algorithm< V, M >::m_tk
private

Definition at line 54 of file Algorithm.h.


The documentation for this class was generated from the following files:

Generated on Tue Jun 5 2018 19:49:21 for queso-0.57.1 by  doxygen 1.8.5