queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
InvLogitGaussianVectorRealizer.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 <cmath>
26 #include <limits>
27 #include <queso/InvLogitGaussianVectorRealizer.h>
28 #include <queso/GslVector.h>
29 #include <queso/GslMatrix.h>
30 #include <queso/math_macros.h>
31 
32 namespace QUESO {
33 
34 template<class V, class M>
36  const char * prefix,
37  const BoxSubset<V, M> & unifiedImageBoxSubset,
38  const V & lawExpVector,
39  const M & lowerCholLawCovMatrix)
40  : BaseVectorRealizer<V, M>(((std::string)(prefix)+"invlogit_gau").c_str(),
41  unifiedImageBoxSubset, std::numeric_limits<unsigned int>::max()),
42  m_unifiedLawExpVector(new V(lawExpVector)),
43  m_unifiedLawVarVector(
44  unifiedImageBoxSubset.vectorSpace().newVector(INFINITY)), // FIX ME
45  m_lowerCholLawCovMatrix(new M(lowerCholLawCovMatrix)),
46  m_matU(NULL),
47  m_vecSsqrt(NULL),
48  m_matVt(NULL),
49  m_unifiedImageBoxSubset(unifiedImageBoxSubset)
50 {
51  *m_unifiedLawExpVector = lawExpVector;
52 }
53 
54 template<class V, class M>
56  const char * prefix,
57  const BoxSubset<V, M> & unifiedImageBoxSubset,
58  const V & lawExpVector,
59  const M & matU,
60  const V & vecSsqrt,
61  const M & matVt)
62  : BaseVectorRealizer<V, M>(((std::string)(prefix)+"invlogit_gau").c_str(),
63  unifiedImageBoxSubset, std::numeric_limits<unsigned int>::max()),
64  m_unifiedLawExpVector(new V(lawExpVector)),
65  m_unifiedLawVarVector(
66  unifiedImageBoxSubset.vectorSpace().newVector( INFINITY)), // FIX ME
67  m_lowerCholLawCovMatrix(NULL),
68  m_matU(new M(matU)),
69  m_vecSsqrt(new V(vecSsqrt)),
70  m_matVt(new M(matVt)),
71  m_unifiedImageBoxSubset(unifiedImageBoxSubset)
72 {
73  *m_unifiedLawExpVector = lawExpVector; // ????
74 }
75 
76 template<class V, class M>
78 {
79  delete m_matVt;
80  delete m_vecSsqrt;
81  delete m_matU;
82  delete m_lowerCholLawCovMatrix;
83  delete m_unifiedLawVarVector;
84  delete m_unifiedLawExpVector;
85 }
86 
87 template <class V, class M>
88 const V &
90 {
91  return *m_unifiedLawExpVector;
92 }
93 
94 template <class V, class M>
95 const V &
97 {
98  return *m_unifiedLawVarVector;
99 }
100 
101 template<class V, class M>
102 void
104 {
105  V iidGaussianVector(m_unifiedImageSet.vectorSpace().zeroVector());
106 
107  iidGaussianVector.cwSetGaussian(0.0, 1.0);
108 
109  if (m_lowerCholLawCovMatrix) {
110  nextValues = (*m_unifiedLawExpVector) +
111  (*m_lowerCholLawCovMatrix) * iidGaussianVector;
112  }
113  else if (m_matU && m_vecSsqrt && m_matVt) {
114  nextValues = (*m_unifiedLawExpVector) +
115  (*m_matU) * ((*m_vecSsqrt) * ((*m_matVt) * iidGaussianVector));
116  }
117  else {
118  queso_error_msg("GaussianVectorRealizer<V,M>::realization() inconsistent internal state");
119  }
120 
121  V min_domain_bounds(this->m_unifiedImageBoxSubset.minValues());
122  V max_domain_bounds(this->m_unifiedImageBoxSubset.maxValues());
123 
124  for (unsigned int i = 0; i < nextValues.sizeLocal(); i++) {
125  double temp = std::exp(nextValues[i]);
126  double min_val = min_domain_bounds[i];
127  double max_val = max_domain_bounds[i];
128 
129  if (queso_isfinite(min_val) &&
130  queso_isfinite(max_val)) {
131  // Left- and right-hand sides are finite. Do full transform.
132  nextValues[i] = (max_val * temp + min_val) / (1.0 + temp);
133  }
134  else if (queso_isfinite(min_val) &&
135  !queso_isfinite(max_val)) {
136  // Left-hand side finite, but right-hand side is not.
137  // Do only left-hand transform.
138  nextValues[i] = temp + min_val;
139  }
140  else if (!queso_isfinite(min_val) &&
141  queso_isfinite(max_val)) {
142  // Right-hand side is finite, but left-hand side is not.
143  // Do only right-hand transform.
144  nextValues[i] = (max_val * temp - 1.0) / temp;
145  }
146  }
147 }
148 
149 template<class V, class M>
150 void
152  const V & newLawExpVector)
153 {
154  // delete old expected values (allocated at construction or last call to this function)
155  delete m_unifiedLawExpVector;
156 
157  m_unifiedLawExpVector = new V(newLawExpVector);
158 }
159 
160 template<class V, class M>
161 void
163  const M & newLowerCholLawCovMatrix)
164 {
165  // delete old expected values (allocated at construction or last call to this function)
166  delete m_lowerCholLawCovMatrix;
167  delete m_matU;
168  delete m_vecSsqrt;
169  delete m_matVt;
170 
171  m_lowerCholLawCovMatrix = new M(newLowerCholLawCovMatrix);
172  m_matU = NULL;
173  m_vecSsqrt = NULL;
174  m_matVt = NULL;
175 }
176 
177 template<class V, class M>
178 void
180  const M & matU,
181  const V & vecSsqrt,
182  const M & matVt)
183 {
184  // delete old expected values (allocated at construction or last call to this function)
185  delete m_lowerCholLawCovMatrix;
186  delete m_matU;
187  delete m_vecSsqrt;
188  delete m_matVt;
189 
190  m_lowerCholLawCovMatrix = NULL;
191  m_matU = new M(matU);
192  m_vecSsqrt = new V(vecSsqrt);
193  m_matVt = new M(matVt);
194 }
195 
196 } // End namespace QUESO
197 
void updateLawExpVector(const V &newLawExpVector)
Updates the mean of the Gaussian with the new value newLawExpVector.
A templated (base) class for handling sampling from vector RVs.
bool queso_isfinite(T arg)
Definition: math_macros.h:51
void updateLowerCholLawCovMatrix(const M &newLowerCholLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
const V & unifiedLawVarVector() const
Access to the vector of variance values and private attribute: m_unifiedLawVarVector.
A class for handling sampling from (transformed) Gaussian probability density distributions with boun...
Class representing a subset of a vector space shaped like a hypercube.
Definition: BoxSubset.h:44
InvLogitGaussianVectorRealizer(const char *prefix, const BoxSubset< V, M > &unifiedImageBoxSubset, const V &lawExpVector, const M &lowerCholLawCovMatrix)
Constructor.
void realization(V &nextValues) const
Draws a realization.
const V & unifiedLawExpVector() const
Access to the vector of mean values of the Gaussian and private attribute: m_unifiedLawExpVector.

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