queso-0.53.0
GammaJointPdf.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-2015 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/GammaJointPdf.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  const V& a,
37  const V& b)
38  :
39  BaseJointPdf<V,M>(((std::string)(prefix)+"uni").c_str(),domainSet),
40  m_a(a),
41  m_b(b)
42 {
43  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
44  *m_env.subDisplayFile() << "Entering GammaJointPdf<V,M>::constructor()"
45  << ": prefix = " << m_prefix
46  << std::endl;
47  }
48 
49  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
50  *m_env.subDisplayFile() << "Leaving GammaJointPdf<V,M>::constructor()"
51  << ": prefix = " << m_prefix
52  << std::endl;
53  }
54 }
55 // Destructor --------------------------------------
56 template<class V,class M>
58 {
59 }
60 // Math methods-------------------------------------
61 template<class V, class M>
62 double
64  const V& domainVector,
65  const V* domainDirection,
66  V* gradVector,
67  M* hessianMatrix,
68  V* hessianEffect) const
69 {
70  queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
71 
72  queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
73 
74  // No need to multiply by exp(m_logOfNormalizationFactor) because 'lnValue()' is called [PDF-06]
75  return exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
76 }
77 //--------------------------------------------------
78 template<class V, class M>
79 double
81  const V& domainVector,
82  const V* domainDirection,
83  V* gradVector,
84  M* hessianMatrix,
85  V* hessianEffect) const
86 {
87  queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
88 
89  double aux = 0.;
90  double result = 0.;
91  for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
92  if (m_normalizationStyle == 0) {
93  //aux = log(gsl_ran_gamma_pdf(domainVector[i],m_a[i],m_b[i]));
94  aux = log(m_env.basicPdfs()->gammaPdfActualValue(domainVector[i],m_a[i],m_b[i]));
95  }
96  else {
97  aux = (m_a[i]-1.)*log(domainVector[i]) - domainVector[i]/m_b[i];
98  }
99  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
100  *m_env.subDisplayFile() << "In GammaJointPdf<V,M>::lnValue()"
101  << ", m_normalizationStyle = " << m_normalizationStyle
102  << ": domainVector[" << i << "] = " << domainVector[i]
103  << ", m_a[" << i << "] = " << m_a[i]
104  << ", m_b[" << i << "] = " << m_b[i]
105  << ", log(pdf)= " << aux
106  << std::endl;
107  }
108  result += aux;
109  }
110  result += m_logOfNormalizationFactor; // [PDF-06]
111 
112  return result;
113 }
114 //--------------------------------------------------
115 template<class V, class M>
116 double
117 GammaJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
118 {
119  double value = 0.;
120 
121  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
122  *m_env.subDisplayFile() << "Entering GammaJointPdf<V,M>::computeLogOfNormalizationFactor()"
123  << std::endl;
124  }
125  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
126  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
127  *m_env.subDisplayFile() << "Leaving GammaJointPdf<V,M>::computeLogOfNormalizationFactor()"
128  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
129  << std::endl;
130  }
131 
132  return value;
133 }
134 
135 } // End namespace QUESO
136 
unsigned int displayVerbosity() const
Definition: Environment.C:396
A templated class for handling sets.
Definition: VectorSet.h:52
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the Gamma PDF.
Definition: GammaJointPdf.C:80
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:56
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:77
GammaJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &a, const V &b)
Constructor.
Definition: GammaJointPdf.C:33
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
const BaseEnvironment & m_env
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:274
A class for handling Gamma joint PDFs.
Definition: GammaJointPdf.h:52
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the Gamma PDF.
Definition: GammaJointPdf.C:63
~GammaJointPdf()
Destructor.
Definition: GammaJointPdf.C:57

Generated on Thu Jun 11 2015 13:52:32 for queso-0.53.0 by  doxygen 1.8.5