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
142 queso_require_msg(!(hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
144 double returnValue = 0.;
146 if (this->m_domainSet.contains(domainVector) ==
false) {
153 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
156 (*gradVector) *= returnValue;
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
189 queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
191 double returnValue = 0.;
193 double lnDeterminant = 0.;
194 if (this->m_domainSet.contains(domainVector) ==
false) {
196 returnValue = -INFINITY;
199 V diffVec(domainVector - this->lawExpVector());
200 if (m_diagonalCovMatrix) {
201 returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents();
215 (*gradVector) = diffVec;
216 (*gradVector) /= this->lawVarVector();
217 (*gradVector) *= -1.0;
220 if (m_normalizationStyle == 0) {
221 unsigned int iMax = this->lawVarVector().sizeLocal();
222 for (
unsigned int i = 0; i < iMax; ++i) {
223 lnDeterminant += std::log(this->lawVarVector()[i]);
228 V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
229 returnValue = (diffVec*tmpVec).sumOfComponents();
240 (*gradVector) = tmpVec;
241 (*gradVector) *= -1.0;
244 if (m_normalizationStyle == 0) {
245 lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
248 if (m_normalizationStyle == 0) {
249 returnValue += ((double) this->lawVarVector().sizeLocal()) * std::log(2*M_PI);
250 returnValue += lnDeterminant;
254 returnValue += m_logOfNormalizationFactor;
256 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
257 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::lnValue()"
258 <<
", m_normalizationStyle = " << m_normalizationStyle
259 <<
", m_diagonalCovMatrix = " << m_diagonalCovMatrix
260 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
261 <<
", lnDeterminant = " << lnDeterminant
262 <<
", meanVector = " << *m_lawExpVector
263 <<
", lawCovMatrix = " << *m_lawCovMatrix
264 <<
": domainVector = " << domainVector
265 <<
", returnValue = " << returnValue
272 template<
class V,
class M>
278 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
279 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
283 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
284 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
285 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
292 template<
class V,
class M>
297 delete m_lawExpVector;
298 m_lawExpVector =
new V(newLawExpVector);
302 template<
class V,
class M>
307 delete m_lawCovMatrix;
308 m_lawCovMatrix =
new M(newLawCovMatrix);
312 template<
class V,
class M>
316 return *m_lawCovMatrix;
unsigned int displayVerbosity() const
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
~GaussianJointPdf()
Destructor.
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).
A templated class for handling sets.
A templated (base) class for handling joint PDFs.
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
GaussianJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
#define queso_require_msg(asserted, msg)
#define queso_require_equal_to_msg(expr1, expr2, msg)
const V & lawExpVector() const
Access to the vector of mean values and private attribute: m_lawExpVector.
const BaseEnvironment & m_env
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
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:
A class for handling Gaussian joint PDFs.
void updateLawExpVector(const V &newLawExpVector)
Updates the mean with the new value newLawExpVector.