queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
JeffreysJointPdf.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/JeffreysJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Constructor -------------------------------------
32 template<class V,class M>
34  const char* prefix,
35  const VectorSet<V,M>& domainSet)
36  :
37  BaseJointPdf<V,M>(((std::string)(prefix)+"jef").c_str(),
38  domainSet)
39 {
40  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
41  *m_env.subDisplayFile() << "Entering JeffreysJointPdf<V,M>::constructor()"
42  << ": prefix = " << m_prefix
43  << std::endl;
44  }
45 
46  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
47  *m_env.subDisplayFile() << "Leaving JeffreysJointPdf<V,M>::constructor()"
48  << ": prefix = " << m_prefix
49  << std::endl;
50  }
51 }
52 // Destructor --------------------------------------
53 template<class V,class M>
55 {
56 }
57 // Math methods-------------------------------------
58 template<class V, class M>
59 double
61  const V& domainVector,
62  const V* domainDirection,
63  V* gradVector,
64  M* hessianMatrix,
65  V* hessianEffect) const
66 {
67  if (domainVector.sizeLocal() != this->m_domainSet.vectorSpace().dimLocal()) {
68  queso_error_msg("There is an invalid input while computing JeffreysJointPdf<V,M>::actualValue()");
69  }
70 
71  if (gradVector ) *gradVector = m_domainSet.vectorSpace().zeroVector();
72  if (hessianMatrix) *hessianMatrix *= 0.;
73  if (hessianEffect) *hessianEffect = m_domainSet.vectorSpace().zeroVector();
74 
75  if (domainDirection) {}; // just to remove compiler warning
76 
77  double pdf = 1.0;
78  for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i){
79  if (domainVector[i] < 0.0 ) {
80  queso_error_msg("The domain for Jeffreys prior should be greater than zero.");
81  }
82  else if ((domainVector[i] == -INFINITY ) ||
83  (domainVector[i] == INFINITY ) ||
84  (m_normalizationStyle != 0 )) {//TODO: R: not sure what this is doing?
85  pdf = 0.0;
86  }
87  else {
88  pdf = pdf * (1.0 / domainVector[i]);
89  }
90  }
91  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
92  *m_env.subDisplayFile() << " return pdf " << std::endl;
93  }
94  return pdf; // No need to multiply by exp(m_logOfNormalizationFactor) [PDF-04]
95 }
96 
97 //--------------------------------------------------
98 
99 template<class V, class M>
100 double
102  const V& domainVector,
103  const V* domainDirection,
104  V* gradVector,
105  M* hessianMatrix,
106  V* hessianEffect) const
107 {
108  if (gradVector ) *gradVector = m_domainSet.vectorSpace().zeroVector();
109  if (hessianMatrix) *hessianMatrix *= 0.0;
110  if (hessianEffect) *hessianEffect = m_domainSet.vectorSpace().zeroVector();
111 
112  if (domainVector[0]) {}; // just to remove compiler warning
113  if (domainDirection) {}; // just to remove compiler warning
114 
115  double pdf = 1.0;
116  double result = 0.;
117  for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i){
118  if (domainVector[i] < 0.0) {
119  queso_error_msg("The domain for Jeffreys prior should be greater than zero.");
120  }
121  else if ((domainVector[i] == -INFINITY ) ||
122  (domainVector[i] == INFINITY ) ||
123  (m_normalizationStyle != 0 )) {//TODO: R: not sure what this is doing?
124  pdf = 0.0;
125  result = -INFINITY; //TODO: what do we do here?
126  }
127  else {
128  pdf = pdf * (1.0 / domainVector[i]);
129  result = std::log(pdf);
130  }
131  }
132  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
133  *m_env.subDisplayFile() << " return log(pdf) " << std::endl;
134  }
135  return result; // No need to add m_logOfNormalizationFactor [PDF-04]
136 }
137 //--------------------------------------------------
138 template<class V, class M>
139 void
141 {
142  // FIXME: This is an improper prior; the "mean" makes little sense.
143  // At least we can return something within the domain set.
144  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
145  *m_env.subDisplayFile() << "Warning: JeffreysJointPdf<V,M>::distributionMean() makes little sense"
146  << std::endl;
147  }
148  m_domainSet.centroid(meanVector);
149 }
150 //--------------------------------------------------
151 template<class V, class M>
152 void
154 {
155  // There's no way this is anything like well-defined
156  queso_not_implemented();
157 }
158 //--------------------------------------------------
159 template<class V, class M>
160 double
161 JeffreysJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
162 {
163  double value = 0.;
164 
165  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
166  *m_env.subDisplayFile() << "Entering JeffreysJointPdf<V,M>::computeLogOfNormalizationFactor()"
167  << std::endl;
168  }
169  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
170  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
171  *m_env.subDisplayFile() << "Leaving JeffreysJointPdf<V,M>::computeLogOfNormalizationFactor()"
172  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
173  << std::endl;
174  }
175 
176  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
177  *m_env.subDisplayFile() << " return normalization factor " << std::endl;
178  }
179  return value;
180 }
181 
182 } // End namespace QUESO
183 
virtual void distributionMean(V &meanVector) const
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the jeffreys PDF.
JeffreysJointPdf(const char *prefix, const VectorSet< V, M > &domainSet)
Constructor.
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the jeffreys PDF.
A templated class for handling sets.
Definition: VectorSet.h:52
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
Definition: JointPdf.C:115
const BaseEnvironment & m_env
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:55
A class for handling jeffreys joint PDFs.
~JeffreysJointPdf()
Destructor.
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
unsigned int displayVerbosity() const
Definition: Environment.C:450
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

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