queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
GaussianVectorRealizer.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 <limits>
26 #include <queso/GaussianVectorRealizer.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 
30 namespace QUESO {
31 
32 // Constructor -------------------------------------
33 template<class V, class M>
35  const VectorSet<V,M>& unifiedImageSet,
36  const V& lawExpVector,
37  const M& lowerCholLawCovMatrix)
38  :
39  BaseVectorRealizer<V,M>( ((std::string)(prefix)+"gau").c_str(), unifiedImageSet, std::numeric_limits<unsigned int>::max()), // 2011/Oct/02 - Correction thanks to Corey
40  m_unifiedLawExpVector (new V(lawExpVector)),
41  m_unifiedLawVarVector (unifiedImageSet.vectorSpace().newVector( INFINITY)), // FIX ME
42  m_lowerCholLawCovMatrix(new M(lowerCholLawCovMatrix)),
43  m_matU (NULL),
44  m_vecSsqrt (NULL),
45  m_matVt (NULL)
46 {
47  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
48  *m_env.subDisplayFile() << "Entering GaussianVectorRealizer<V,M>::constructor() [1]"
49  << ": prefix = " << m_prefix
50  << std::endl;
51  }
52 
53  *m_unifiedLawExpVector = lawExpVector; // ????
54 
55  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
56  *m_env.subDisplayFile() << "Leaving GaussianVectorRealizer<V,M>::constructor() [1]"
57  << ": prefix = " << m_prefix
58  << std::endl;
59  }
60 }
61 // Constructor -------------------------------------
62 template<class V, class M>
64  const VectorSet<V,M>& unifiedImageSet,
65  const V& lawExpVector,
66  const M& matU,
67  const V& vecSsqrt,
68  const M& matVt)
69  :
70  BaseVectorRealizer<V,M>( ((std::string)(prefix)+"gau").c_str(), unifiedImageSet, std::numeric_limits<unsigned int>::max()), // 2011/Oct/02 - Correction thanks to Corey
71  m_unifiedLawExpVector (new V(lawExpVector)),
72  m_unifiedLawVarVector (unifiedImageSet.vectorSpace().newVector( INFINITY)), // FIX ME
73  m_lowerCholLawCovMatrix(NULL),
74  m_matU (new M(matU)),
75  m_vecSsqrt (new V(vecSsqrt)),
76  m_matVt (new M(matVt))
77 {
78  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
79  *m_env.subDisplayFile() << "Entering GaussianVectorRealizer<V,M>::constructor() [2]"
80  << ": prefix = " << m_prefix
81  << std::endl;
82  }
83 
84  *m_unifiedLawExpVector = lawExpVector; // ????
85 
86  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
87  *m_env.subDisplayFile() << "Leaving GaussianVectorRealizer<V,M>::constructor() [2]"
88  << ": prefix = " << m_prefix
89  << std::endl;
90  }
91 }
92 // Destructor --------------------------------------
93 template<class V, class M>
95 {
96  delete m_matVt;
97  delete m_vecSsqrt;
98  delete m_matU;
99  delete m_lowerCholLawCovMatrix;
100  delete m_unifiedLawVarVector;
101  delete m_unifiedLawExpVector;
102 }
103 // Realization-related methods----------------------
104 template <class V, class M>
105 const V&
107 {
108  return *m_unifiedLawExpVector;
109 }
110 //--------------------------------------------------
111 template <class V, class M>
112 const V&
114 {
115  return *m_unifiedLawVarVector;
116 }
117 //--------------------------------------------------
118 template<class V, class M>
119 void
121 {
122  V iidGaussianVector(m_unifiedImageSet.vectorSpace().zeroVector());
123 
124  bool outOfSupport = true;
125  do {
126  iidGaussianVector.cwSetGaussian(0.0, 1.0);
127 
128  if (m_lowerCholLawCovMatrix) {
129  nextValues = (*m_unifiedLawExpVector) + (*m_lowerCholLawCovMatrix)*iidGaussianVector;
130  }
131  else if (m_matU && m_vecSsqrt && m_matVt) {
132  nextValues = (*m_unifiedLawExpVector) + (*m_matU)*( (*m_vecSsqrt) * ((*m_matVt)*iidGaussianVector) );
133  }
134  else {
135  queso_error_msg("inconsistent internal state");
136  }
137 
138  outOfSupport = !(this->m_unifiedImageSet.contains(nextValues));
139  } while (outOfSupport); // prudenci 2011-Oct-04
140 
141  return;
142 }
143 //--------------------------------------------------
144 template<class V, class M>
145 void
147 {
148  // delete old expected values (allocated at construction or last call to this function)
149  delete m_unifiedLawExpVector;
150 
151  m_unifiedLawExpVector = new V(newLawExpVector);
152 
153  return;
154 }
155 //--------------------------------------------------
156 template<class V, class M>
157 void
159 {
160  // delete old expected values (allocated at construction or last call to this function)
161  delete m_lowerCholLawCovMatrix;
162  delete m_matU;
163  delete m_vecSsqrt;
164  delete m_matVt;
165 
166  m_lowerCholLawCovMatrix = new M(newLowerCholLawCovMatrix);
167  m_matU = NULL;
168  m_vecSsqrt = NULL;
169  m_matVt = NULL;
170 
171  return;
172 }
173 //--------------------------------------------------
174 template<class V, class M>
175 void
177  const M& matU,
178  const V& vecSsqrt,
179  const M& matVt)
180 {
181  // delete old expected values (allocated at construction or last call to this function)
182  delete m_lowerCholLawCovMatrix;
183  delete m_matU;
184  delete m_vecSsqrt;
185  delete m_matVt;
186 
187  m_lowerCholLawCovMatrix = NULL;
188  m_matU = new M(matU);
189  m_vecSsqrt = new V(vecSsqrt);
190  m_matVt = new M(matVt);
191 
192  return;
193 }
194 
195 } // End namespace QUESO
196 
A templated (base) class for handling sampling from vector RVs.
const BaseEnvironment & m_env
GaussianVectorRealizer(const char *prefix, const VectorSet< V, M > &unifiedImageSet, const V &lawExpVector, const M &lowerCholLawCovMatrix)
Constructor.
A templated class for handling sets.
Definition: VectorSet.h:52
const V & unifiedLawExpVector() const
Access to the vector of mean values and private attribute: m_unifiedLawExpVector. ...
void updateLowerCholLawCovMatrix(const M &newLowerCholLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
A class for handling sampling from Gaussian probability density distributions.
void realization(V &nextValues) const
Draws a realization.
unsigned int displayVerbosity() const
Definition: Environment.C:450
const V & unifiedLawVarVector() const
Access to the vector of variance values and private attribute: m_unifiedLawVarVector.
void updateLawExpVector(const V &newLawExpVector)
Updates the mean with the new value newLawExpVector.
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