25 #include <queso/GaussianJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
32 template<
class V,
class M>
36 const V& lawExpVector,
37 const V& lawVarVector)
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))
67 template<
class V,
class M>
71 const V& lawExpVector,
72 const M& lawCovMatrix)
74 BaseJointPdf<V,M>(((std::string)(prefix)+
"gau").c_str(),domainSet),
75 m_lawExpVector (new V(lawExpVector)),
76 m_lawVarVector (domainSet.vectorSpace().newVector(INFINITY)),
77 m_diagonalCovMatrix(false),
78 m_lawCovMatrix (new M(lawCovMatrix))
90 <<
", Covariance Matrix = " << lawCovMatrix
101 template<
class V,
class M>
104 delete m_lawCovMatrix;
105 delete m_lawVarVector;
106 delete m_lawExpVector;
109 template <
class V,
class M>
113 return *m_lawExpVector;
116 template <
class V,
class M>
120 return *m_lawVarVector;
123 template<
class V,
class M>
126 const V& domainVector,
127 const V* domainDirection,
130 V* hessianEffect)
const
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
140 UQ_FATAL_TEST_MACRO(domainVector.sizeLocal() != this->m_domainSet.vectorSpace().dimLocal(),
142 "GaussianJointPdf<V,M>::actualValue()",
147 "GaussianJointPdf<V,M>::actualValue()",
148 "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
150 double returnValue = 0.;
152 if (this->m_domainSet.contains(domainVector) ==
false) {
156 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
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
172 template<
class V,
class M>
175 const V& domainVector,
176 const V* domainDirection,
179 V* hessianEffect)
const
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
191 "GaussianJointPdf<V,M>::lnValue()",
192 "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
194 if (domainDirection) {};
196 double returnValue = 0.;
198 double lnDeterminant = 0.;
199 if (this->m_domainSet.contains(domainVector) ==
false) {
200 returnValue = -INFINITY;
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]);
214 V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
215 returnValue = (diffVec*tmpVec).sumOfComponents();
216 if (m_normalizationStyle == 0) {
217 lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
220 if (m_normalizationStyle == 0) {
221 returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2*M_PI);
222 returnValue += lnDeterminant;
226 returnValue += m_logOfNormalizationFactor;
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
244 template<
class V,
class M>
250 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
251 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
255 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
256 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
257 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
264 template<
class V,
class M>
269 delete m_lawExpVector;
270 m_lawExpVector =
new V(newLawExpVector);
274 template<
class V,
class M>
279 delete m_lawCovMatrix;
280 m_lawCovMatrix =
new M(newLawCovMatrix);
284 template<
class V,
class M>
288 return *m_lawCovMatrix;
GaussianJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
~GaussianJointPdf()
Destructor.
A templated class for handling sets.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
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...
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
unsigned int displayVerbosity() const
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
A templated (base) class for handling joint PDFs.
const V & lawExpVector() const
Access to the vector of mean values and private attribute: m_lawExpVector.