queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
InvLogitGaussianVectorRV.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/InvLogitGaussianVectorRV.h>
26 #include <queso/InvLogitGaussianVectorRealizer.h>
27 #include <queso/InvLogitGaussianJointPdf.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 BoxSubset<V, M> & imageBoxSubset,
38  const V & lawExpVector,
39  const V & lawVarVector)
40  : BaseVectorRV<V, M>(((std::string)(prefix)+"invlogit_gau").c_str(),
41  imageBoxSubset)
42 {
43  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
44  *m_env.subDisplayFile() << "Entering InvLogitGaussianVectorRV<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 
52  dynamic_cast<const BoxSubset<V, M> & >(m_imageSet), lawExpVector,
53  lawVarVector);
54 
55  V cholDiag(lawVarVector);
56  cholDiag.cwSqrt();
57  M lowerCholLawCovMatrix(cholDiag);
58  lowerCholLawCovMatrix.zeroUpper(false);
59 
61  dynamic_cast<const BoxSubset<V, M> & >(m_imageSet), lawExpVector,
62  lowerCholLawCovMatrix);
63 
64  m_subCdf = NULL; // FIX ME: complete code
65  m_unifiedCdf = NULL; // FIX ME: complete code
66  m_mdf = NULL; // FIX ME: complete code
67 
68  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
69  *m_env.subDisplayFile() << "Leaving InvLogitGaussianVectorRV<V,M>::constructor() [1]"
70  << ": prefix = " << m_prefix
71  << std::endl;
72  }
73 }
74 
75 template<class V, class M>
77  const char * prefix,
78  const BoxSubset<V, M> & imageBoxSubset,
79  const V & lawExpVector,
80  const M & lawCovMatrix)
81  : BaseVectorRV<V, M>(((std::string)(prefix)+"invlogit_gau").c_str(),
82  imageBoxSubset)
83 {
84  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
85  *m_env.subDisplayFile() << "Entering InvLogitGaussianVectorRV<V,M>::constructor() [2]"
86  << ": prefix = " << m_prefix
87  << std::endl;
88  }
89 
91  dynamic_cast<const BoxSubset<V, M> & >(m_imageSet),
92  lawExpVector, lawCovMatrix);
93 
94  M lowerCholLawCovMatrix(lawCovMatrix);
95  int iRC = lowerCholLawCovMatrix.chol();
96  lowerCholLawCovMatrix.zeroUpper(false);
97  if (iRC) {
98  std::cerr << "In InvLogitGaussianVectorRV<V,M>::constructor() [2]: chol failed, will use svd\n";
99  if (m_env.subDisplayFile()) {
100  *m_env.subDisplayFile() << "In InvLogitGaussianVectorRV<V,M>::constructor() [2]: chol failed; will use svd; lawCovMatrix contents are\n";
101  *m_env.subDisplayFile() << lawCovMatrix; // FIX ME: might demand parallelism
102  *m_env.subDisplayFile() << std::endl;
103  }
104  M matU (lawCovMatrix);
105  M matVt(m_imageSet.vectorSpace().zeroVector());
106  V vecS (m_imageSet.vectorSpace().zeroVector());
107  iRC = lawCovMatrix.svd(matU,vecS,matVt);
108  queso_require_msg(!(iRC), "Cholesky decomposition of covariance matrix failed.");
109 
110  vecS.cwSqrt();
112  dynamic_cast<const BoxSubset<V, M> & >(m_imageSet), lawExpVector, matU,
113  vecS, matVt);
114  }
115  else {
117  dynamic_cast<const BoxSubset<V, M> & >(m_imageSet), lawExpVector, lowerCholLawCovMatrix);
118  }
119 
120  m_subCdf = NULL; // FIX ME: complete code
121  m_unifiedCdf = NULL; // FIX ME: complete code
122  m_mdf = NULL; // FIX ME: complete code
123 
124  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
125  *m_env.subDisplayFile() << "Leaving InvLogitGaussianVectorRV<V,M>::constructor() [2]"
126  << ": prefix = " << m_prefix
127  << std::endl;
128  }
129 }
130 
131 template<class V, class M>
133 {
134  delete m_mdf;
135  delete m_unifiedCdf;
136  delete m_subCdf;
137  delete m_realizer;
138  delete m_pdf;
139 }
140 
141 template<class V, class M>
142 void
144 {
145  // We are sure that m_pdf (and m_realizer, etc) point to associated Gaussian
146  // classes, so all is well
147  (dynamic_cast<InvLogitGaussianJointPdf<V, M> * >(m_pdf))->updateLawExpVector(
148  newLawExpVector);
149  (dynamic_cast<InvLogitGaussianVectorRealizer<V, M> * >(m_realizer))->
150  updateLawExpVector(newLawExpVector);
151  return;
152 }
153 
154 template<class V, class M>
155 void
157 {
158  // We are sure that m_pdf (and m_realizer, etc) point to associated Gaussian
159  // classes, so all is well
160  (dynamic_cast<InvLogitGaussianJointPdf<V,M> * >(m_pdf))->updateLawCovMatrix(
161  newLawCovMatrix);
162 
163  M newLowerCholLawCovMatrix(newLawCovMatrix);
164  int iRC = newLowerCholLawCovMatrix.chol();
165  newLowerCholLawCovMatrix.zeroUpper(false);
166  if (iRC) {
167  std::cerr << "In InvLogitGaussianVectorRV<V,M>::updateLawCovMatrix(): chol failed, will use svd\n";
168  if (m_env.subDisplayFile()) {
169  *m_env.subDisplayFile() << "In InvLogitGaussianVectorRV<V,M>::updateLawCovMatrix(): chol failed; will use svd; newLawCovMatrix contents are\n";
170  *m_env.subDisplayFile() << newLawCovMatrix; // FIX ME: might demand parallelism
171  *m_env.subDisplayFile() << std::endl;
172  }
173  M matU (newLawCovMatrix);
174  M matVt(m_imageSet.vectorSpace().zeroVector());
175  V vecS (m_imageSet.vectorSpace().zeroVector());
176  iRC = newLawCovMatrix.svd(matU,vecS,matVt);
177  queso_require_msg(!(iRC), "Cholesky decomposition of covariance matrix failed.");
178 
179  vecS.cwSqrt();
180  (dynamic_cast<InvLogitGaussianVectorRealizer<V, M> * >(m_realizer))->
181  updateLowerCholLawCovMatrix(matU, vecS, matVt);
182  }
183  else {
184  (dynamic_cast<InvLogitGaussianVectorRealizer<V, M> * >(m_realizer))->
185  updateLowerCholLawCovMatrix(newLowerCholLawCovMatrix);
186  }
187  return;
188 }
189 
190 template <class V, class M>
191 void
193 {
194  os << "InvLogitGaussianVectorRV<V,M>::print() says, 'Please implement me.'" << std::endl;
195  return;
196 }
197 
198 } // End namespace QUESO
199 
std::string m_prefix
Definition: VectorRV.h:118
void print(std::ostream &os) const
TODO: Prints the vector RV.
A class representing a (transformed) Gaussian vector RV with bounds.
InvLogitGaussianVectorRV(const char *prefix, const BoxSubset< V, M > &imageBoxSubset, const V &lawExpVector, const V &lawVarVector)
Constructor.
const BaseVectorCdf< V, M > * m_unifiedCdf
Definition: VectorRV.h:123
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the covariance matrix.
const BaseEnvironment & m_env
Definition: VectorRV.h:117
A class for handling sampling from (transformed) Gaussian probability density distributions with boun...
BaseVectorRealizer< V, M > * m_realizer
Definition: VectorRV.h:121
Class representing a subset of a vector space shaped like a hypercube.
Definition: BoxSubset.h:44
virtual ~InvLogitGaussianVectorRV()
Virtual destructor.
A class for handling hybrid (transformed) Gaussians with bounds.
BaseJointPdf< V, M > * m_pdf
Definition: VectorRV.h:120
A templated base class for handling vector RV.
Definition: VectorRV.h:53
const BaseVectorCdf< V, M > * m_subCdf
Definition: VectorRV.h:122
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
void updateLawExpVector(const V &newLawExpVector)
Updates the vector that contains the mean values for the underlying Gaussian.
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