queso-0.51.1
GaussianJointPdf.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,2009,2010,2011,2012,2013 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/GaussianJointPdf.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& lawExpVector,
37  const V& lawVarVector)
38  :
39  BaseJointPdf<V,M>(((std::string)(prefix)+"gau").c_str(),domainSet),
40  m_lawExpVector (new V(lawExpVector)),
41  m_lawVarVector (new V(lawVarVector)),
42  m_diagonalCovMatrix(true),
43  m_lawCovMatrix (m_domainSet.vectorSpace().newDiagMatrix(lawVarVector))
44 {
45 
46  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
47  *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::constructor() [1]"
48  << ": prefix = " << m_prefix
49  << std::endl;
50  }
51 
52  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
53  *m_env.subDisplayFile() << "In GaussianJointPdf<V,M>::constructor()"
54  //<< ", prefix = " << m_prefix
55  << ": meanVector = " << this->lawExpVector()
56  << ", Variances = " << this->lawVarVector()
57  << std::endl;
58  }
59 
60  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
61  *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::constructor() [1]"
62  << ": prefix = " << m_prefix
63  << std::endl;
64  }
65 }
66 // Constructor -------------------------------------
67 template<class V,class M>
69  const char* prefix,
70  const VectorSet<V,M>& domainSet,
71  const V& lawExpVector,
72  const M& lawCovMatrix)
73  :
74  BaseJointPdf<V,M>(((std::string)(prefix)+"gau").c_str(),domainSet),
75  m_lawExpVector (new V(lawExpVector)),
76  m_lawVarVector (domainSet.vectorSpace().newVector(INFINITY)), // FIX ME
77  m_diagonalCovMatrix(false),
78  m_lawCovMatrix (new M(lawCovMatrix))
79 {
80  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
81  *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::constructor() [2]"
82  << ": prefix = " << m_prefix
83  << std::endl;
84  }
85 
86  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
87  *m_env.subDisplayFile() << "In GaussianJointPdf<V,M>::constructor()"
88  //<< ", prefix = " << m_prefix
89  << ": meanVector = " << this->lawExpVector()
90  << ", Covariance Matrix = " << lawCovMatrix
91  << std::endl;
92  }
93 
94  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
95  *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::constructor() [2]"
96  << ": prefix = " << m_prefix
97  << std::endl;
98  }
99 }
100 // Destructor --------------------------------------
101 template<class V,class M>
103 {
104  delete m_lawCovMatrix;
105  delete m_lawVarVector;
106  delete m_lawExpVector;
107 }
108 // Math methods-------------------------------------
109 template <class V, class M>
110 const V&
112 {
113  return *m_lawExpVector;
114 }
115 //--------------------------------------------------
116 template <class V, class M>
117 const V&
119 {
120  return *m_lawVarVector;
121 }
122 //--------------------------------------------------
123 template<class V, class M>
124 double
126  const V& domainVector,
127  const V* domainDirection,
128  V* gradVector,
129  M* hessianMatrix,
130  V* hessianEffect) const
131 {
132  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
133  *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::actualValue()"
134  << ", meanVector = " << *m_lawExpVector
135  << ", lawCovMatrix = " << *m_lawCovMatrix
136  << ": domainVector = " << domainVector
137  << std::endl;
138  }
139 
140  UQ_FATAL_TEST_MACRO(domainVector.sizeLocal() != this->m_domainSet.vectorSpace().dimLocal(),
141  m_env.worldRank(),
142  "GaussianJointPdf<V,M>::actualValue()",
143  "invalid input");
144 
145  UQ_FATAL_TEST_MACRO((gradVector || hessianMatrix || hessianEffect),
146  m_env.worldRank(),
147  "GaussianJointPdf<V,M>::actualValue()",
148  "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
149 
150  double returnValue = 0.;
151 
152  if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04
153  returnValue = 0.;
154  }
155  else {
156  returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
157  }
158  //returnValue *= exp(m_logOfNormalizationFactor); // No need, because 'lnValue()' is called right above [PDF-03]
159 
160  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
161  *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::actualValue()"
162  << ", meanVector = " << *m_lawExpVector
163  << ", lawCovMatrix = " << *m_lawCovMatrix
164  << ": domainVector = " << domainVector
165  << ", returnValue = " << returnValue
166  << std::endl;
167  }
168 
169  return returnValue;
170 }
171 //--------------------------------------------------
172 template<class V, class M>
173 double
175  const V& domainVector,
176  const V* domainDirection,
177  V* gradVector,
178  M* hessianMatrix,
179  V* hessianEffect) const
180 {
181  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
182  *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::lnValue()"
183  << ", meanVector = " << *m_lawExpVector
184  << ", lawCovMatrix = " << *m_lawCovMatrix
185  << ": domainVector = " << domainVector
186  << std::endl;
187  }
188 
189  UQ_FATAL_TEST_MACRO((gradVector || hessianMatrix || hessianEffect),
190  m_env.worldRank(),
191  "GaussianJointPdf<V,M>::lnValue()",
192  "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
193 
194  if (domainDirection) {}; // just to remove compiler warning
195 
196  double returnValue = 0.;
197 
198  double lnDeterminant = 0.;
199  if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04
200  returnValue = -INFINITY;
201  }
202  else {
203  V diffVec(domainVector - this->lawExpVector());
204  if (m_diagonalCovMatrix) {
205  returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents();
206  if (m_normalizationStyle == 0) {
207  unsigned int iMax = this->lawVarVector().sizeLocal();
208  for (unsigned int i = 0; i < iMax; ++i) {
209  lnDeterminant += log(this->lawVarVector()[i]);
210  }
211  }
212  }
213  else {
214  V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
215  returnValue = (diffVec*tmpVec).sumOfComponents();
216  if (m_normalizationStyle == 0) {
217  lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
218  }
219  }
220  if (m_normalizationStyle == 0) {
221  returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2*M_PI); // normalization of pdf
222  returnValue += lnDeterminant; // normalization of pdf
223  }
224  returnValue *= -0.5;
225  }
226  returnValue += m_logOfNormalizationFactor; // [PDF-03]
227 
228  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
229  *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::lnValue()"
230  << ", m_normalizationStyle = " << m_normalizationStyle
231  << ", m_diagonalCovMatrix = " << m_diagonalCovMatrix
232  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
233  << ", lnDeterminant = " << lnDeterminant
234  << ", meanVector = " << *m_lawExpVector
235  << ", lawCovMatrix = " << *m_lawCovMatrix
236  << ": domainVector = " << domainVector
237  << ", returnValue = " << returnValue
238  << std::endl;
239  }
240 
241  return returnValue;
242 }
243 //--------------------------------------------------
244 template<class V, class M>
245 double
246 GaussianJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
247 {
248  double value = 0.;
249 
250  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
251  *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
252  << std::endl;
253  }
254  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
255  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
256  *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
257  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
258  << std::endl;
259  }
260 
261  return value;
262 }
263 //--------------------------------------------------
264 template<class V, class M>
265 void
267 {
268  // delete old expected values (allocated at construction or last call to this function)
269  delete m_lawExpVector;
270  m_lawExpVector = new V(newLawExpVector);
271  return;
272 }
273 
274 template<class V, class M>
275 void
277 {
278  // delete old expected values (allocated at construction or last call to this function)
279  delete m_lawCovMatrix;
280  m_lawCovMatrix = new M(newLawCovMatrix);
281  return;
282 }
283 
284 template<class V, class M>
285 const M&
287 {
288  return *m_lawCovMatrix;
289 }
290 
291 } // End namespace QUESO
292 
GaussianJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
~GaussianJointPdf()
Destructor.
A templated class for handling sets.
Definition: VectorSet.h:49
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the Gaussian PDF (scalar function).
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
void updateLawExpVector(const V &newLawExpVector)
Updates the mean with the new value newLawExpVector.
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the Gaussian PDF:
const BaseEnvironment & m_env
A class for handling Gaussian joint PDFs.
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
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
unsigned int displayVerbosity() const
Definition: Environment.C:436
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:223
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:60
const V & lawExpVector() const
Access to the vector of mean values and private attribute: m_lawExpVector.

Generated on Thu Apr 23 2015 19:26:15 for queso-0.51.1 by  doxygen 1.8.5