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(!(gradVector || hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
144 double returnValue = 0.;
146 if (this->m_domainSet.contains(domainVector) ==
false) {
150 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
154 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
155 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::actualValue()"
156 <<
", meanVector = " << *m_lawExpVector
157 <<
", lawCovMatrix = " << *m_lawCovMatrix
158 <<
": domainVector = " << domainVector
159 <<
", returnValue = " << returnValue
166 template<
class V,
class M>
169 const V& domainVector,
170 const V* domainDirection,
173 V* hessianEffect)
const
175 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
176 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::lnValue()"
177 <<
", meanVector = " << *m_lawExpVector
178 <<
", lawCovMatrix = " << *m_lawCovMatrix
179 <<
": domainVector = " << domainVector
183 queso_require_msg(!(gradVector || hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
185 if (domainDirection) {};
187 double returnValue = 0.;
189 double lnDeterminant = 0.;
190 if (this->m_domainSet.contains(domainVector) ==
false) {
191 returnValue = -INFINITY;
194 V diffVec(domainVector - this->lawExpVector());
195 if (m_diagonalCovMatrix) {
196 returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents();
197 if (m_normalizationStyle == 0) {
198 unsigned int iMax = this->lawVarVector().sizeLocal();
199 for (
unsigned int i = 0; i < iMax; ++i) {
200 lnDeterminant += log(this->lawVarVector()[i]);
205 V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
206 returnValue = (diffVec*tmpVec).sumOfComponents();
207 if (m_normalizationStyle == 0) {
208 lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
211 if (m_normalizationStyle == 0) {
212 returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2*M_PI);
213 returnValue += lnDeterminant;
217 returnValue += m_logOfNormalizationFactor;
219 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
220 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::lnValue()"
221 <<
", m_normalizationStyle = " << m_normalizationStyle
222 <<
", m_diagonalCovMatrix = " << m_diagonalCovMatrix
223 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
224 <<
", lnDeterminant = " << lnDeterminant
225 <<
", meanVector = " << *m_lawExpVector
226 <<
", lawCovMatrix = " << *m_lawCovMatrix
227 <<
": domainVector = " << domainVector
228 <<
", returnValue = " << returnValue
235 template<
class V,
class M>
241 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
242 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
246 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
247 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
248 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
255 template<
class V,
class M>
260 delete m_lawExpVector;
261 m_lawExpVector =
new V(newLawExpVector);
265 template<
class V,
class M>
270 delete m_lawCovMatrix;
271 m_lawCovMatrix =
new M(newLawCovMatrix);
275 template<
class V,
class M>
279 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.