queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
BayesianJointPdf.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/BayesianJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Default constructor -----------------------------
32 template<class V,class M>
34  const char* prefix,
35  const BaseJointPdf <V,M>& priorDensity,
36  const BaseScalarFunction<V,M>& likelihoodFunction,
37  double likelihoodExponent,
38  const VectorSet <V,M>& intersectionDomain)
39  :
40  BaseJointPdf<V,M>(((std::string)(prefix)+"bay").c_str(),intersectionDomain),
41  m_priorDensity (priorDensity),
42  m_likelihoodFunction (likelihoodFunction),
43  m_likelihoodExponent (likelihoodExponent),
44  m_lastComputedLogPrior (0.),
45  m_lastComputedLogLikelihood(0.),
46  m_tmpVector1 (m_domainSet.vectorSpace().zeroVector()),
47  m_tmpVector2 (m_domainSet.vectorSpace().zeroVector()),
48  m_tmpMatrix (m_domainSet.vectorSpace().newMatrix())
49 {
50 }
51 // Destructor ---------------------------------------
52 template<class V,class M>
54 {
55  delete m_tmpMatrix;
56 }
57 // Math methods -------------------------------------
58 template<class V,class M>
59 void
61 {
62  m_priorDensity.setNormalizationStyle(value);
63  return;
64 }
65 // --------------------------------------------------
66 template<class V,class M>
67 double
69 {
70  return m_lastComputedLogPrior;
71 }
72 // --------------------------------------------------
73 template<class V,class M>
74 double
76 {
77  return m_lastComputedLogLikelihood;
78 }
79 // --------------------------------------------------
80 template<class V, class M>
81 double
83  const V& domainVector,
84  const V* domainDirection,
85  V* gradVector,
86  M* hessianMatrix,
87  V* hessianEffect) const
88 {
89  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
90  *m_env.subDisplayFile() << "Entering BayesianJointPdf<V,M>::actualValue()"
91  << ": domainVector = " << domainVector
92  << std::endl;
93  }
94 
95  queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
96 
97  V* gradVLike = NULL;
98  if (gradVector) gradVLike = &m_tmpVector1;
99 
100  M* hessianMLike = NULL;
101  if (hessianMatrix) hessianMLike = m_tmpMatrix;
102 
103  V* hessianELike = NULL;
104  if (hessianEffect) hessianELike = &m_tmpVector2;
105 
106  double value1 = m_priorDensity.actualValue (domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect);
107  double value2 = 1.;
108  if (m_likelihoodExponent != 0.) {
109  value2 = m_likelihoodFunction.actualValue(domainVector,domainDirection,gradVLike ,hessianMLike ,hessianELike );
110  }
111 
112  queso_require_msg(!(gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
113 
114  double returnValue = value1;
115  if (m_likelihoodExponent == 0.) {
116  // Do nothing
117  }
118  else if (m_likelihoodExponent == 1.) {
119  returnValue *= value2;
120  }
121  else {
122  returnValue *= pow(value2,m_likelihoodExponent);
123  }
124  returnValue *= exp(m_logOfNormalizationFactor); // [PDF-02] ???
125 
126  m_lastComputedLogPrior = log(value1);
127  m_lastComputedLogLikelihood = m_likelihoodExponent*log(value2);
128 
129  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
130  *m_env.subDisplayFile() << "Leaving BayesianJointPdf<V,M>::actualValue()"
131  << ": domainVector = " << domainVector
132  << ", returnValue = " << returnValue
133  << std::endl;
134  }
135 
136  return returnValue;
137 }
138 
139 template<class V, class M>
140 double
141 BayesianJointPdf<V,M>::lnValue(const V & domainVector) const
142 {
143  double value1 = m_priorDensity.lnValue(domainVector);
144 
145  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
146  *m_env.subDisplayFile() << "In BayesianJointPdf<V,M>::lnValue()"
147  << ", domainVector = " << domainVector
148  << ": lnPrior = " << value1
149  << std::endl;
150  }
151 
152  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
153  *m_env.subDisplayFile() << "In BayesianJointPdf<V,M>::lnValue()"
154  << ", domainVector = " << domainVector
155  << ": about to call likelihood()"
156  << std::endl;
157  }
158 
159  double value2 = 0.;
160  if (m_likelihoodExponent != 0.) {
161  value2 = m_likelihoodFunction.lnValue(domainVector);
162  }
163 
164  double returnValue = value1;
165  if (m_likelihoodExponent == 0.) {
166  // Do nothing
167  }
168  else if (m_likelihoodExponent == 1.) {
169  returnValue += value2;
170  }
171  else {
172  returnValue += value2*m_likelihoodExponent;
173  } // prudenci 2010/03/05
174  returnValue += m_logOfNormalizationFactor; // [PDF-02] ???
175 
176  m_lastComputedLogPrior = value1;
177  m_lastComputedLogLikelihood = m_likelihoodExponent*value2;
178 
179  return returnValue;
180 }
181 
182 template<class V, class M>
183 double
184 BayesianJointPdf<V,M>::lnValue(const V & domainVector, V & gradVector) const
185 {
186  double value1 = m_priorDensity.lnValue(domainVector, gradVector);
187 
188  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
189  *m_env.subDisplayFile() << "In BayesianJointPdf<V,M>::lnValue()"
190  << ", domainVector = " << domainVector
191  << ": lnPrior = " << value1
192  << std::endl;
193  }
194 
195  double value2 = 0.;
196  if (m_likelihoodExponent != 0.) {
197  value2 = m_likelihoodFunction.lnValue(domainVector, m_tmpVector1);
198  }
199 
200  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
201  *m_env.subDisplayFile() << "In BayesianJointPdf<V,M>::lnValue()"
202  << ", domainVector = " << domainVector
203  << ": value1 = " << value1
204  << ", value2 = " << value2
205  << std::endl;
206  *m_env.subDisplayFile() << "In BayesianJointPdf<V,M>::lnValue()"
207  << ", domainVector = " << domainVector
208  << ": gradVector = " << gradVector
209  << ", gradVLike = " << m_tmpVector1
210  << std::endl;
211  }
212 
213  gradVector += m_tmpVector1;
214 
215  double returnValue = value1;
216  if (m_likelihoodExponent == 0.) {
217  // Do nothing
218  }
219  else if (m_likelihoodExponent == 1.) {
220  returnValue += value2;
221  }
222  else {
223  returnValue += value2*m_likelihoodExponent;
224  } // prudenci 2010/03/05
225  returnValue += m_logOfNormalizationFactor; // [PDF-02] ???
226 
227  m_lastComputedLogPrior = value1;
228  m_lastComputedLogLikelihood = m_likelihoodExponent*value2;
229 
230  return returnValue;
231 }
232 
233 // --------------------------------------------------
234 template<class V, class M>
235 double
236 BayesianJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
237 {
238  double value = 0.;
239 
240  double volume = m_domainSet.volume();
241  if ((queso_isnan(volume)) ||
242  (volume == -INFINITY ) ||
243  (volume == INFINITY ) ||
244  (volume <= 0. )) {
245  // Do nothing
246  }
247  else {
248  queso_error_msg("incomplete code for computeLogOfNormalizationFactor()");
249  }
250 
251  return value;
252 }
253 
254 } // End namespace QUESO
255 
A templated (base) class for handling scalar functions.
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
TODO: Computes the logarithm of the normalization factor.
double lastComputedLogLikelihood() const
Returns the logarithm of the last computed likelihood value. Access to protected attribute m_lastComp...
virtual double lnValue(const V &domainVector) const
Computes the logarithm of the value of the function.
bool queso_isnan(T arg)
Definition: math_macros.h:39
A templated class for handling sets.
Definition: VectorSet.h:52
double lastComputedLogPrior() const
Returns the logarithm of the last computed Prior value. Access to protected attribute m_lastComputedL...
void setNormalizationStyle(unsigned int value) const
Sets a value to be used in the normalization style of the prior density PDF (ie, protected attribute ...
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the PDF (scalar function).
BayesianJointPdf(const char *prefix, const BaseJointPdf< V, M > &priorDensity, const BaseScalarFunction< V, M > &likelihoodFunction, double likelihoodExponent, const VectorSet< V, M > &intersectionDomain)
Default constructor.
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:55
~BayesianJointPdf()
Destructor.
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))

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