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>
129 os <<
"Start printing GaussianJointPdf<V, M>" << std::endl;
130 os <<
"m_prefix:" << std::endl;
131 os << this->m_prefix << std::endl;
132 os <<
"m_domainSet:" << std::endl;
133 os << this->m_domainSet << std::endl;
134 os <<
"m_normalizationStyle:" << std::endl;
135 os << this->m_normalizationStyle << std::endl;
136 os <<
"m_logOfNormalizationFactor:" << std::endl;
137 os << this->m_logOfNormalizationFactor << std::endl;
138 os <<
"Mean:" << std::endl;
139 os << this->lawExpVector() << std::endl;
140 os <<
"Variance vector:" << std::endl;
141 os << this->lawVarVector() << std::endl;
142 os <<
"Covariance matrix:" << std::endl;
143 os << this->lawCovMatrix() << std::endl;
144 os <<
"Diagonal covariance?" << std::endl;
145 os << this->m_diagonalCovMatrix << std::endl;
146 os <<
"End printing GaussianJointPdf<V, M>" << std::endl;
150 template<
class V,
class M>
153 const V& domainVector,
154 const V* domainDirection,
157 V* hessianEffect)
const
159 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
160 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::actualValue()"
161 <<
", meanVector = " << *m_lawExpVector
162 <<
", lawCovMatrix = " << *m_lawCovMatrix
163 <<
": domainVector = " << domainVector
169 queso_require_msg(!(hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
171 double returnValue = 0.;
173 if (this->m_domainSet.contains(domainVector) ==
false) {
180 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
183 (*gradVector) *= returnValue;
187 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
188 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::actualValue()"
189 <<
", meanVector = " << *m_lawExpVector
190 <<
", lawCovMatrix = " << *m_lawCovMatrix
191 <<
": domainVector = " << domainVector
192 <<
", returnValue = " << returnValue
199 template<
class V,
class M>
202 const V& domainVector,
203 const V* domainDirection,
206 V* hessianEffect)
const
208 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
209 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::lnValue()"
210 <<
", meanVector = " << *m_lawExpVector
211 <<
", lawCovMatrix = " << *m_lawCovMatrix
212 <<
": domainVector = " << domainVector
216 queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
218 double returnValue = 0.;
220 double lnDeterminant = 0.;
221 if (this->m_domainSet.contains(domainVector) ==
false) {
223 returnValue = -INFINITY;
226 V diffVec(domainVector - this->lawExpVector());
227 if (m_diagonalCovMatrix) {
228 returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents();
242 (*gradVector) = diffVec;
243 (*gradVector) /= this->lawVarVector();
244 (*gradVector) *= -1.0;
247 if (m_normalizationStyle == 0) {
248 unsigned int iMax = this->lawVarVector().sizeLocal();
249 for (
unsigned int i = 0; i < iMax; ++i) {
250 lnDeterminant += std::log(this->lawVarVector()[i]);
255 V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
256 returnValue = (diffVec*tmpVec).sumOfComponents();
267 (*gradVector) = tmpVec;
268 (*gradVector) *= -1.0;
271 if (m_normalizationStyle == 0) {
272 lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
275 if (m_normalizationStyle == 0) {
276 returnValue += ((double) this->lawVarVector().sizeLocal()) * std::log(2*M_PI);
277 returnValue += lnDeterminant;
281 returnValue += m_logOfNormalizationFactor;
283 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
284 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::lnValue()"
285 <<
", m_normalizationStyle = " << m_normalizationStyle
286 <<
", m_diagonalCovMatrix = " << m_diagonalCovMatrix
287 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
288 <<
", lnDeterminant = " << lnDeterminant
289 <<
", meanVector = " << *m_lawExpVector
290 <<
", lawCovMatrix = " << *m_lawCovMatrix
291 <<
": domainVector = " << domainVector
292 <<
", returnValue = " << returnValue
299 template<
class V,
class M>
303 meanVector = this->lawExpVector();
306 template<
class V,
class M>
310 queso_assert_equal_to (covMatrix.numCols(), covMatrix.numRowsGlobal());
312 if (m_diagonalCovMatrix) {
313 covMatrix.zeroLower();
314 covMatrix.zeroUpper();
316 unsigned int n_comp = this->lawVarVector().sizeLocal();
317 queso_assert_equal_to (n_comp, covMatrix.numCols());
319 for (
unsigned int i = 0; i < n_comp; ++i) {
320 covMatrix(i,i) = this->lawVarVector()[i];
323 covMatrix = *this->m_lawCovMatrix;
327 template<
class V,
class M>
333 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
334 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
338 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
339 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
340 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
347 template<
class V,
class M>
352 delete m_lawExpVector;
353 m_lawExpVector =
new V(newLawExpVector);
357 template<
class V,
class M>
362 delete m_lawCovMatrix;
363 m_lawCovMatrix =
new M(newLawCovMatrix);
367 template<
class V,
class M>
371 return *m_lawCovMatrix;
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
GaussianJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
const V & lawExpVector() const
Access to the vector of mean values and private attribute: m_lawExpVector.
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 BaseEnvironment & m_env
A templated (base) class for handling joint PDFs.
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
virtual void print(std::ostream &os) const
Print method for informational and logging purposes.
A templated class for handling sets.
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the Gaussian PDF:
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
void updateLawExpVector(const V &newLawExpVector)
Updates the mean with the new value newLawExpVector.
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"))
~GaussianJointPdf()
Destructor.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
A class for handling Gaussian joint PDFs.
unsigned int displayVerbosity() const