queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
Algorithm.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/Environment.h>
26 #include <queso/math_macros.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 #include <queso/Algorithm.h>
30 #include <queso/TKGroup.h>
31 #include <queso/InvLogitGaussianJointPdf.h>
32 
33 namespace QUESO {
34 
35 template <class V, class M>
37  const BaseTKGroup<V, M> & tk)
38  :
39  m_env(env),
40  m_tk(tk)
41 {
42 }
43 
44 template <class V, class M>
46 {
47 }
48 
49 template <class V, class M>
50 double
54  const V & tk_pos_x,
55  const V & tk_pos_y)
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 }
149 
150 template class Algorithm<GslVector, GslMatrix>;
151 
152 } // End namespace QUESO
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
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.
Definition: Algorithm.C:51
bool queso_isnan(T arg)
Definition: math_macros.h:39
const V & vecValues() const
Values of the chain (vector); access to private attribute m_vecValues.
Algorithm(const BaseEnvironment &env, const BaseTKGroup< V, M > &tk)
Definition: Algorithm.C:36
double logTarget() const
Logarithm of the value of the target; access to private attribute m_logTarget.
This base class allows the representation of a transition kernel.
Definition: Algorithm.h:32
A templated class that represents a Markov Chain.
bool outOfTargetSupport() const
Whether or not a position is out of target support; access to private attribute m_outOfTargetSupport...

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5