25 #include <queso/GaussianVectorRV.h>
26 #include <queso/GaussianVectorRealizer.h>
27 #include <queso/GaussianJointPdf.h>
28 #include <queso/GslVector.h>
29 #include <queso/GslMatrix.h>
34 template<
class V,
class M>
38 const V& lawExpVector,
39 const V& lawVarVector)
41 BaseVectorRV<V,M>(((std::string)(prefix)+
"gau").c_str(),imageSet)
51 "GaussianVectorRV<V,M>::constructor() [1]",
52 "Covariance matrix is not symmetric positive definite.");
59 V cholDiag(lawVarVector);
61 M lowerCholLawCovMatrix(cholDiag);
62 lowerCholLawCovMatrix.zeroUpper(
false);
67 lowerCholLawCovMatrix);
80 template<
class V,
class M>
84 const V& lawExpVector,
85 const M& lawCovMatrix)
87 BaseVectorRV<V,M>(((std::string)(prefix)+
"gau").c_str(),imageSet)
100 M lowerCholLawCovMatrix(lawCovMatrix);
101 int iRC = lowerCholLawCovMatrix.chol();
102 lowerCholLawCovMatrix.zeroUpper(
false);
104 std::cerr <<
"In GaussianVectorRV<V,M>::constructor() [2]: chol failed, will use svd\n";
106 *
m_env.
subDisplayFile() <<
"In GaussianVectorRV<V,M>::constructor() [2]: chol failed; will use svd; lawCovMatrix contents are\n";
110 M matU (lawCovMatrix);
111 M matVt(
m_imageSet.vectorSpace().zeroVector());
112 V vecS (
m_imageSet.vectorSpace().zeroVector());
113 iRC = lawCovMatrix.svd(matU,vecS,matVt);
116 "GaussianVectorRV<V,M>::constructor() [2]",
117 "Cholesky decomposition of covariance matrix failed.");
131 lowerCholLawCovMatrix);
145 template<
class V,
class M>
155 template<
class V,
class M>
165 template<
class V,
class M>
172 M newLowerCholLawCovMatrix(newLawCovMatrix);
173 int iRC = newLowerCholLawCovMatrix.chol();
174 newLowerCholLawCovMatrix.zeroUpper(
false);
176 std::cerr <<
"In GaussianVectorRV<V,M>::updateLawCovMatrix(): chol failed, will use svd\n";
177 if (m_env.subDisplayFile()) {
178 *m_env.subDisplayFile() <<
"In GaussianVectorRV<V,M>::updateLawCovMatrix(): chol failed; will use svd; newLawCovMatrix contents are\n";
179 *m_env.subDisplayFile() << newLawCovMatrix;
180 *m_env.subDisplayFile() << std::endl;
182 M matU (newLawCovMatrix);
183 M matVt(m_imageSet.vectorSpace().zeroVector());
184 V vecS (m_imageSet.vectorSpace().zeroVector());
185 iRC = newLawCovMatrix.svd(matU,vecS,matVt);
188 "GaussianVectorRV<V,M>::updateLawCovMatrix()",
189 "Cholesky decomposition of covariance matrix failed.");
202 template <
class V,
class M>
206 os <<
"GaussianVectorRV<V,M>::print() says, 'Please implement me.'" << std::endl;
214 template<
class V,
class M>
225 M& sigmaMat11_cond_on_2)
228 unsigned int dim1 = muVec1.sizeLocal();
229 unsigned int dim2 = muVec2.sizeLocal();
233 "ComputeConditionalGaussianVectorRV()",
234 "invalid sigmaMat11");
238 "ComputeConditionalGaussianVectorRV()",
239 "invalid sigmaMat12");
243 "ComputeConditionalGaussianVectorRV()",
244 "invalid sigmaMat21");
248 "ComputeConditionalGaussianVectorRV()",
249 "invalid sigmaMat22");
252 M mat_tt(sigmaMat12);
254 mat_tt.fillWithTranspose(0,0,sigmaMat21,
true,
true);
255 double auxNorm = (mat_tt - sigmaMat12).normFrob();
256 if (auxNorm >= 1.e-12) {
258 *env.
subDisplayFile() <<
"In ComputeConditionalGaussianVectorRV()"
259 <<
": WARNING, ||sigmaMat21^T - sigmaMat12||_2 = " << auxNorm
265 "ComputeConditionalGaussianVectorRV()",
266 "sigmaMat12 and sigmaMat21 are not transpose of each other");
270 "ComputeConditionalGaussianVectorRV()",
271 "invalid sampleVec2");
275 "ComputeConditionalGaussianVectorRV()",
276 "invalid muVec1_cond_on_2");
278 UQ_FATAL_TEST_MACRO((sigmaMat11_cond_on_2.numRowsLocal() != dim1) || (sigmaMat11_cond_on_2.numCols() != dim1),
280 "ComputeConditionalGaussianVectorRV()",
281 "invalid sigmaMat11_cond_on_2");
283 muVec1_cond_on_2 = muVec1 + sigmaMat12 * sigmaMat22.invertMultiply(sampleVec2 - muVec2);
284 sigmaMat11_cond_on_2 = sigmaMat11 - sigmaMat12 * sigmaMat22.invertMultiply(sigmaMat21);
Class for matrix operations using GSL library.
int worldRank() const
Returns the process world rank.
void ComputeConditionalGaussianVectorRV(const V &muVec1, const V &muVec2, const M &sigmaMat11, const M &sigmaMat12, const M &sigmaMat21, const M &sigmaMat22, const V &sampleVec2, V &muVec1_cond_on_2, M &sigmaMat11_cond_on_2)
A templated class for handling sets.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
A class for handling sampling from Gaussian probability density distributions.
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
A class for handling Gaussian joint PDFs.
void print(std::ostream &os) const
TODO: Prints the vector RV.
const BaseVectorMdf< V, M > * m_mdf
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
void updateLawExpVector(const V &newLawExpVector)
Updates the vector that contains the mean values.
Class for vector operations using GSL library.
const VectorSet< V, M > & m_imageSet
A templated base class for handling vector RV.
BaseJointPdf< V, M > * m_pdf
const BaseVectorCdf< V, M > * m_subCdf
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the covariance matrix.
const BaseVectorCdf< V, M > * m_unifiedCdf
const BaseEnvironment & m_env
virtual ~GaussianVectorRV()
Virtual destructor.
unsigned int displayVerbosity() const
A class representing a Gaussian vector RV.
GaussianVectorRV(const char *prefix, const VectorSet< V, M > &imageSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
BaseVectorRealizer< V, M > * m_realizer