queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
GaussianVectorRV.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
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>
30 
31 namespace QUESO {
32 
33 // Constructor---------------------------------------
34 template<class V, class M>
36  const char* prefix,
37  const VectorSet<V,M>& imageSet,
38  const V& lawExpVector,
39  const V& lawVarVector)
40  :
41  BaseVectorRV<V,M>(((std::string)(prefix)+"gau").c_str(),imageSet)
42 {
43  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
44  *m_env.subDisplayFile() << "Entering GaussianVectorRV<V,M>::constructor() [1]"
45  << ": prefix = " << m_prefix
46  << std::endl;
47  }
48 
49  queso_require_greater_msg(lawVarVector.getMinValue(), 0.0, "Covariance matrix is not symmetric positive definite.");
50 
51  m_pdf = new GaussianJointPdf<V,M>(m_prefix.c_str(),
52  m_imageSet,
53  lawExpVector,
54  lawVarVector);
55 
56  V cholDiag(lawVarVector);
57  cholDiag.cwSqrt();
58  M lowerCholLawCovMatrix(cholDiag);
59  lowerCholLawCovMatrix.zeroUpper(false);
60 
62  m_imageSet,
63  lawExpVector,
64  lowerCholLawCovMatrix);
65 
66  m_subCdf = NULL; // FIX ME: complete code
67  m_unifiedCdf = NULL; // FIX ME: complete code
68  m_mdf = NULL; // FIX ME: complete code
69 
70  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
71  *m_env.subDisplayFile() << "Leaving GaussianVectorRV<V,M>::constructor() [1]"
72  << ": prefix = " << m_prefix
73  << std::endl;
74  }
75 }
76 // Constructor---------------------------------------
77 template<class V, class M>
79  const char* prefix,
80  const VectorSet<V,M>& imageSet,
81  const V& lawExpVector,
82  const M& lawCovMatrix)
83  :
84  BaseVectorRV<V,M>(((std::string)(prefix)+"gau").c_str(),imageSet)
85 {
86  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
87  *m_env.subDisplayFile() << "Entering GaussianVectorRV<V,M>::constructor() [2]"
88  << ": prefix = " << m_prefix
89  << std::endl;
90  }
91 
92  m_pdf = new GaussianJointPdf<V,M>(m_prefix.c_str(),
93  m_imageSet,
94  lawExpVector,
95  lawCovMatrix);
96 
97  M lowerCholLawCovMatrix(lawCovMatrix);
98  int iRC = lowerCholLawCovMatrix.chol();
99  lowerCholLawCovMatrix.zeroUpper(false);
100  if (iRC) {
101  std::cerr << "In GaussianVectorRV<V,M>::constructor() [2]: chol failed, will use svd\n";
102  if (m_env.subDisplayFile()) {
103  *m_env.subDisplayFile() << "In GaussianVectorRV<V,M>::constructor() [2]: chol failed; will use svd; lawCovMatrix contents are\n";
104  *m_env.subDisplayFile() << lawCovMatrix; // FIX ME: might demand parallelism
105  *m_env.subDisplayFile() << std::endl;
106  }
107  M matU (lawCovMatrix);
108  M matVt(m_imageSet.vectorSpace().zeroVector());
109  V vecS (m_imageSet.vectorSpace().zeroVector());
110  iRC = lawCovMatrix.svd(matU,vecS,matVt);
111  queso_require_msg(!(iRC), "Cholesky decomposition of covariance matrix failed.");
112 
113  vecS.cwSqrt();
115  m_imageSet,
116  lawExpVector,
117  matU,
118  vecS, // already square rooted
119  matVt);
120  }
121  else {
123  m_imageSet,
124  lawExpVector,
125  lowerCholLawCovMatrix);
126  }
127 
128  m_subCdf = NULL; // FIX ME: complete code
129  m_unifiedCdf = NULL; // FIX ME: complete code
130  m_mdf = NULL; // FIX ME: complete code
131 
132  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
133  *m_env.subDisplayFile() << "Leaving GaussianVectorRV<V,M>::constructor() [2]"
134  << ": prefix = " << m_prefix
135  << std::endl;
136  }
137 }
138 // Destructor ---------------------------------------
139 template<class V, class M>
141 {
142  delete m_mdf;
143  delete m_unifiedCdf;
144  delete m_subCdf;
145  delete m_realizer;
146  delete m_pdf;
147 }
148 // Statistical methods-------------------------------
149 template<class V, class M>
150 void
152 {
153  // We are sure that m_pdf (and m_realizer, etc) point to associated Gaussian classes, so all is well
154  ( dynamic_cast< GaussianJointPdf <V,M>* >(m_pdf ) )->updateLawExpVector(newLawExpVector);
155  ( dynamic_cast< GaussianVectorRealizer<V,M>* >(m_realizer) )->updateLawExpVector(newLawExpVector);
156  return;
157 }
158 //---------------------------------------------------
159 template<class V, class M>
160 void
162 {
163  // We are sure that m_pdf (and m_realizer, etc) point to associated Gaussian classes, so all is well
164  ( dynamic_cast< GaussianJointPdf<V,M>* >(m_pdf) )->updateLawCovMatrix(newLawCovMatrix);
165 
166  M newLowerCholLawCovMatrix(newLawCovMatrix);
167  int iRC = newLowerCholLawCovMatrix.chol();
168  newLowerCholLawCovMatrix.zeroUpper(false);
169  if (iRC) {
170  std::cerr << "In GaussianVectorRV<V,M>::updateLawCovMatrix(): chol failed, will use svd\n";
171  if (m_env.subDisplayFile()) {
172  *m_env.subDisplayFile() << "In GaussianVectorRV<V,M>::updateLawCovMatrix(): chol failed; will use svd; newLawCovMatrix contents are\n";
173  *m_env.subDisplayFile() << newLawCovMatrix; // FIX ME: might demand parallelism
174  *m_env.subDisplayFile() << std::endl;
175  }
176  M matU (newLawCovMatrix);
177  M matVt(m_imageSet.vectorSpace().zeroVector());
178  V vecS (m_imageSet.vectorSpace().zeroVector());
179  iRC = newLawCovMatrix.svd(matU,vecS,matVt);
180  queso_require_msg(!(iRC), "Cholesky decomposition of covariance matrix failed.");
181 
182  vecS.cwSqrt();
183  ( dynamic_cast< GaussianVectorRealizer<V,M>* >(m_realizer) )->updateLowerCholLawCovMatrix(matU,
184  vecS, // already square rooted
185  matVt);
186  }
187  else {
188  ( dynamic_cast< GaussianVectorRealizer<V,M>* >(m_realizer) )->updateLowerCholLawCovMatrix(newLowerCholLawCovMatrix);
189  }
190  return;
191 }
192 // I/O methods---------------------------------------
193 template <class V, class M>
194 void
195 GaussianVectorRV<V,M>::print(std::ostream& os) const
196 {
197  os << "GaussianVectorRV<V,M>::print() says, 'Please implement me.'" << std::endl;
198  return;
199 }
200 
201 
202 //---------------------------------------------------
203 // Method declared outside class definition ---------
204 //---------------------------------------------------
205 template<class V, class M>
206 void
208  const V& muVec1,
209  const V& muVec2,
210  const M& sigmaMat11,
211  const M& sigmaMat12,
212  const M& sigmaMat21,
213  const M& sigmaMat22,
214  const V& sampleVec2,
215  V& muVec1_cond_on_2,
216  M& sigmaMat11_cond_on_2)
217 {
218  const BaseEnvironment& env = muVec1.env();
219  unsigned int dim1 = muVec1.sizeLocal();
220  unsigned int dim2 = muVec2.sizeLocal();
221 
222  queso_require_msg(!((sigmaMat11.numRowsLocal() != dim1) || (sigmaMat11.numCols() != dim1)), "invalid sigmaMat11");
223 
224  queso_require_msg(!((sigmaMat12.numRowsLocal() != dim1) || (sigmaMat12.numCols() != dim2)), "invalid sigmaMat12");
225 
226  queso_require_msg(!((sigmaMat21.numRowsLocal() != dim2) || (sigmaMat21.numCols() != dim1)), "invalid sigmaMat21");
227 
228  queso_require_msg(!((sigmaMat22.numRowsLocal() != dim2) || (sigmaMat22.numCols() != dim2)), "invalid sigmaMat22");
229 
230  // Check transpose operation
231  M mat_tt(sigmaMat12);
232  mat_tt.cwSet(0.);
233  mat_tt.fillWithTranspose(0,0,sigmaMat21,true,true);
234  double auxNorm = (mat_tt - sigmaMat12).normFrob();
235  if (auxNorm >= 1.e-12) {
236  if (env.subDisplayFile()) {
237  *env.subDisplayFile() << "In ComputeConditionalGaussianVectorRV()"
238  << ": WARNING, ||sigmaMat21^T - sigmaMat12||_2 = " << auxNorm
239  << std::endl;
240  }
241  }
242  queso_require_less_msg(auxNorm, 1.e-12, "sigmaMat12 and sigmaMat21 are not transpose of each other");
243 
244  queso_require_equal_to_msg(sampleVec2.sizeLocal(), dim2, "invalid sampleVec2");
245 
246  queso_require_equal_to_msg(muVec1_cond_on_2.sizeLocal(), dim1, "invalid muVec1_cond_on_2");
247 
248  queso_require_msg(!((sigmaMat11_cond_on_2.numRowsLocal() != dim1) || (sigmaMat11_cond_on_2.numCols() != dim1)), "invalid sigmaMat11_cond_on_2");
249 
250  muVec1_cond_on_2 = muVec1 + sigmaMat12 * sigmaMat22.invertMultiply(sampleVec2 - muVec2);
251  sigmaMat11_cond_on_2 = sigmaMat11 - sigmaMat12 * sigmaMat22.invertMultiply(sigmaMat21);
252 
253  return;
254 }
255 
256 } // End namespace QUESO
257 
259 
260 template void QUESO::ComputeConditionalGaussianVectorRV<QUESO::GslVector, QUESO::GslMatrix>(QUESO::GslVector const&, QUESO::GslVector const&, QUESO::GslMatrix const&, QUESO::GslMatrix const&, QUESO::GslMatrix const&, QUESO::GslMatrix const&, QUESO::GslVector const&, QUESO::GslVector&, QUESO::GslMatrix&);
void print(std::ostream &os) const
TODO: Prints the vector RV.
std::string m_prefix
Definition: VectorRV.h:118
A class for handling Gaussian joint PDFs.
const BaseVectorCdf< V, M > * m_unifiedCdf
Definition: VectorRV.h:123
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
A class representing a Gaussian vector RV.
virtual ~GaussianVectorRV()
Virtual destructor.
const BaseEnvironment & m_env
Definition: VectorRV.h:117
Class for matrix operations using GSL library.
Definition: GslMatrix.h:49
A templated class for handling sets.
Definition: VectorSet.h:52
BaseVectorRealizer< V, M > * m_realizer
Definition: VectorRV.h:121
A class for handling sampling from Gaussian probability density distributions.
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)
void updateLawExpVector(const V &newLawExpVector)
Updates the vector that contains the mean values.
Class for vector operations using GSL library.
Definition: GslVector.h:49
BaseJointPdf< V, M > * m_pdf
Definition: VectorRV.h:120
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"))
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the covariance matrix.
A templated base class for handling vector RV.
Definition: VectorRV.h:53
const BaseVectorCdf< V, M > * m_subCdf
Definition: VectorRV.h:122
GaussianVectorRV(const char *prefix, const VectorSet< V, M > &imageSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
unsigned int displayVerbosity() const
Definition: Environment.C:450
const VectorSet< V, M > & m_imageSet
Definition: VectorRV.h:119
const BaseVectorMdf< V, M > * m_mdf
Definition: VectorRV.h:124
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5