queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
InvLogitGaussianJointPdf.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/InvLogitGaussianJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Constructor -------------------------------------
32 template<class V,class M>
34  const char* prefix,
35  const BoxSubset<V,M>& domainBoxSubset,
36  const V& lawExpVector,
37  const V& lawVarVector)
38  :
39  BaseJointPdf<V,M>(((std::string)(prefix)+"invlogit_gau").c_str(),
40  domainBoxSubset),
41  m_lawExpVector(new V(lawExpVector)),
42  m_lawVarVector(new V(lawVarVector)),
43  m_diagonalCovMatrix(true),
44  m_lawCovMatrix(m_domainSet.vectorSpace().newDiagMatrix(lawVarVector)),
45  m_domainBoxSubset(domainBoxSubset)
46 {
47 }
48 
49 
50 template<class V,class M>
52  const char* prefix,
53  const BoxSubset<V,M>& domainBoxSubset,
54  const V& lawExpVector,
55  const M& lawCovMatrix)
56  :
57  BaseJointPdf<V,M>(((std::string)(prefix)+"invlogit_gau").c_str(),
58  domainBoxSubset),
59  m_lawExpVector(new V(lawExpVector)),
60  m_lawVarVector(domainBoxSubset.vectorSpace().newVector(INFINITY)), // FIX ME
61  m_diagonalCovMatrix(false),
62  m_lawCovMatrix(new M(lawCovMatrix)),
63  m_domainBoxSubset(domainBoxSubset)
64 {
65 }
66 
67 template<class V,class M>
69 {
70  delete m_lawCovMatrix;
71  delete m_lawVarVector;
72  delete m_lawExpVector;
73 }
74 
75 template <class V, class M>
76 const V&
78 {
79  return *m_lawExpVector;
80 }
81 
82 template <class V, class M>
83 const V&
85 {
86  return *m_lawVarVector;
87 }
88 
89 template <class V, class M>
90 void
91 InvLogitGaussianJointPdf<V, M>::print(std::ostream & os) const
92 {
93  // Print m_env?
94 
95  os << "Start printing InvLogitGaussianJointPdf<V, M>" << std::endl;
96  os << "m_prefix:" << std::endl;
97  os << this->m_prefix << std::endl;
98  os << "m_domainSet:" << std::endl;
99  os << this->m_domainBoxSubset << std::endl;
100  os << "m_normalizationStyle:" << std::endl;
101  os << this->m_normalizationStyle << std::endl;
102  os << "m_logOfNormalizationFactor:" << std::endl;
103  os << this->m_logOfNormalizationFactor << std::endl;
104  os << "Mean:" << std::endl;
105  os << this->lawExpVector() << std::endl;
106  os << "Variance vector:" << std::endl;
107  os << this->lawVarVector() << std::endl;
108  os << "Covariance matrix:" << std::endl;
109  os << this->lawCovMatrix() << std::endl;
110  os << "Diagonal covariance?" << std::endl;
111  os << this->m_diagonalCovMatrix << std::endl;
112  os << "End printing InvLogitGaussianJointPdf<V, M>" << std::endl;
113 }
114 
115 template<class V, class M>
116 double
118  const V& domainVector,
119  const V* domainDirection,
120  V* gradVector,
121  M* hessianMatrix,
122  V* hessianEffect) const
123 {
124  double returnValue;
125 
126  returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
127 
128  return returnValue;
129 }
130 
131 template<class V, class M>
132 double
134  const V& domainVector,
135  const V* domainDirection,
136  V* gradVector,
137  M* hessianMatrix,
138  V* hessianEffect) const
139 {
140  double returnValue;
141  double lnDeterminant = 0.0;
142  V transformedDomainVector(domainVector);
143 
144  V min_domain_bounds(this->m_domainBoxSubset.minValues());
145  V max_domain_bounds(this->m_domainBoxSubset.maxValues());
146 
147  double lnjacobian = 0.0;
148  for (unsigned int i = 0; i < domainVector.sizeLocal(); i++) {
149  double min_val = min_domain_bounds[i];
150  double max_val = max_domain_bounds[i];
151 
152  if (queso_isfinite(min_val) &&
153  queso_isfinite(max_val)) {
154 
155  if (domainVector[i] == min_val || domainVector[i] == max_val) {
156  // Exit early if we can
157  return -INFINITY;
158  }
159 
160  // Left- and right-hand sides are finite. Do full transform.
161  transformedDomainVector[i] = std::log(domainVector[i] - min_val) -
162  std::log(max_val - domainVector[i]);
163 
164  lnjacobian += std::log(max_val - min_val) -
165  std::log(domainVector[i] - min_val) -
166  std::log(max_val - domainVector[i]);
167  }
168  else if (queso_isfinite(min_val) &&
169  !queso_isfinite(max_val)) {
170 
171  if (domainVector[i] == min_val) {
172  // Exit early if we can
173  return -INFINITY;
174  }
175 
176  // Left-hand side finite, but right-hand side is not.
177  // Do only left-hand transform.
178  transformedDomainVector[i] = std::log(domainVector[i] - min_val);
179 
180  lnjacobian += -std::log(domainVector[i] - min_val);
181  }
182  else if (!queso_isfinite(min_val) &&
183  queso_isfinite(max_val)) {
184 
185  if (domainVector[i] == max_val) {
186  // Exit early if we can
187  return -INFINITY;
188  }
189 
190  // Right-hand side is finite, but left-hand side is not.
191  // Do only right-hand transform.
192  transformedDomainVector[i] = -std::log(max_val - domainVector[i]);
193 
194  lnjacobian += -std::log(max_val - domainVector[i]);
195  }
196  else {
197  // No transform.
198  transformedDomainVector[i] = domainVector[i];
199  }
200  }
201 
202  V diffVec(transformedDomainVector - this->lawExpVector());
203  if (m_diagonalCovMatrix) {
204  returnValue = ((diffVec * diffVec) /
205  this->lawVarVector()).sumOfComponents();
206  if (m_normalizationStyle == 0) {
207  unsigned int iMax = this->lawVarVector().sizeLocal();
208  for (unsigned int i = 0; i < iMax; ++i) {
209  lnDeterminant += log(this->lawVarVector()[i]);
210  }
211  }
212  }
213  else {
214  V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
215  returnValue = (diffVec * tmpVec).sumOfComponents();
216  if (m_normalizationStyle == 0) {
217  lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
218  }
219  }
220  if (m_normalizationStyle == 0) {
221  returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2 * M_PI);
222  returnValue += lnDeterminant;
223  }
224  returnValue *= -0.5;
225  returnValue += m_logOfNormalizationFactor;
226  returnValue += lnjacobian;
227 
228  return returnValue;
229 }
230 
231 template<class V, class M>
232 void
234 {
235  // AFAIK there's no simple closed form here, and taking the inverse
236  // transformation of the mean in the transformed space probably
237  // isn't accurate enough in cases where the mean is too near the
238  // bounds.
239  queso_not_implemented();
240 }
241 
242 //---------------------------------------------------
243 template<class V,class M>
244 void
246 {
247  // AFAIK there's no simple closed form here, and taking the inverse
248  // transformation of the variance in the transformed space probably
249  // isn't accurate enough in cases where the mean is too near the
250  // bounds.
251  queso_not_implemented();
252 }
253 
254 
255 template<class V, class M>
256 double
257 InvLogitGaussianJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
258 {
259  double value = 0.;
260 
261  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
262  *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
263  << std::endl;
264  }
265  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
266  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
267  *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
268  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
269  << std::endl;
270  }
271 
272  return value;
273 }
274 
275 template<class V, class M>
276 void
278 {
279  // delete old expected values (allocated at construction or last call to this function)
280  delete m_lawExpVector;
281  m_lawExpVector = new V(newLawExpVector);
282 }
283 
284 template<class V, class M>
285 void
287 {
288  // delete old expected values (allocated at construction or last call to this function)
289  delete m_lawCovMatrix;
290  m_lawCovMatrix = new M(newLawCovMatrix);
291 }
292 
293 template<class V, class M>
294 const M&
296 {
297  return *m_lawCovMatrix;
298 }
299 
300 } // End namespace QUESO
301 
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
bool queso_isfinite(T arg)
Definition: math_macros.h:51
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the (transformed) Gaussian PDF (scalar function).
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the (transformed) Gaussian PDF.
const V & lawExpVector() const
Access to the vector of mean values of the Gaussian (not transformed) and private attribute: m_lawExp...
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
Definition: JointPdf.C:115
Class representing a subset of a vector space shaped like a hypercube.
Definition: BoxSubset.h:44
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:55
void updateLawExpVector(const V &newLawExpVector)
Updates the mean of the Gaussian (not transformed) with the new value newLawExpVector.
A class for handling hybrid (transformed) Gaussians with bounds.
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
virtual void print(std::ostream &os) const
Prints the distribution.
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
InvLogitGaussianJointPdf(const char *prefix, const BoxSubset< V, M > &domainBoxSubset, const V &lawExpVector, const V &lawVarVector)
Constructor.
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.

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