queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ConcatenatedJointPdf.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/ConcatenatedJointPdf.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 BaseJointPdf<V,M>& density1,
36  const BaseJointPdf<V,M>& density2,
37  const VectorSet <V,M>& concatenatedDomain)
38  :
39  BaseJointPdf<V,M>(((std::string)(prefix)+"concat").c_str(),concatenatedDomain),
40  m_densities (2,(const BaseJointPdf<V,M>*) NULL)
41 {
42  m_densities[0] = &density1;
43  m_densities[1] = &density2;
44 
45  unsigned int size1 = m_densities[0]->domainSet().vectorSpace().dimLocal();
46  unsigned int size2 = m_densities[1]->domainSet().vectorSpace().dimLocal();
47  unsigned int size = concatenatedDomain.vectorSpace().dimLocal();
48 
49  queso_require_equal_to_msg((size1+size2), size, "incompatible dimensions");
50 }
51 // Constructor -------------------------------------
52 template<class V,class M>
54  const char* prefix,
55  const std::vector<const BaseJointPdf<V,M>* >& densities,
56  const VectorSet<V,M>& concatenatedDomain)
57  :
58  BaseJointPdf<V,M>(((std::string)(prefix)+"concat").c_str(),concatenatedDomain),
59  m_densities (densities.size(),(const BaseJointPdf<V,M>*) NULL)
60 {
61  unsigned int sumSizes = 0;
62  for (unsigned i = 0; i < m_densities.size(); ++i) {
63  m_densities[i] = densities[i];
64  sumSizes += m_densities[i]->domainSet().vectorSpace().dimLocal();
65  }
66 
67  unsigned int size = concatenatedDomain.vectorSpace().dimLocal();
68 
69  queso_require_equal_to_msg(sumSizes, size, "incompatible dimensions");
70 }
71 // Destructor --------------------------------------
72 template<class V,class M>
74 {
75 }
76 // Math methods-------------------------------------
77 template<class V,class M>
78 void
80 {
81  for (unsigned i = 0; i < m_densities.size(); ++i) {
82  m_densities[i]->setNormalizationStyle(value);
83  }
84  return;
85 }
86 //--------------------------------------------------
87 template<class V, class M>
88 double
90  const V& domainVector,
91  const V* domainDirection,
92  V* gradVector,
93  M* hessianMatrix,
94  V* hessianEffect) const
95 {
96  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
97  *m_env.subDisplayFile() << "Entering ConcatenatedJointPdf<V,M>::actualValue()"
98  << ": domainVector = " << domainVector
99  << std::endl;
100  }
101 
102  queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
103 
104  queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
105 
106  double returnValue = 1.;
107  unsigned int cumulativeSize = 0;
108  for (unsigned int i = 0; i < m_densities.size(); ++i) {
109  V vec_i(m_densities[i]->domainSet().vectorSpace().zeroVector());
110  domainVector.cwExtract(cumulativeSize,vec_i);
111  double value_i = m_densities[i]->actualValue(vec_i,NULL,NULL,NULL,NULL);
112  returnValue *= value_i;
113  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
114  *m_env.subDisplayFile() << "In ConcatenatedJointPdf<V,M>::actualValue()"
115  << ", vec_" << i << ") = " << vec_i
116  << ": value_" << i << " = " << value_i
117  << ", temporary cumulative value = " << returnValue
118  << std::endl;
119  }
120  cumulativeSize += vec_i.sizeLocal();
121  }
122  //returnValue *= exp(m_logOfNormalizationFactor); // No need, because each PDF should be already normalized [PDF-11]
123 
124  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
125  *m_env.subDisplayFile() << "Leaving ConcatenatedJointPdf<V,M>::actualValue()"
126  << ": domainVector = " << domainVector
127  << ", returnValue = " << returnValue
128  << std::endl;
129  }
130 
131  return returnValue;
132 }
133 //--------------------------------------------------
134 template<class V, class M>
135 double
137  const V& domainVector,
138  const V* domainDirection,
139  V* gradVector,
140  M* hessianMatrix,
141  V* hessianEffect) const
142 {
143  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
144  *m_env.subDisplayFile() << "Entering ConcatenatedJointPdf<V,M>::lnValue()"
145  << ": domainVector = " << domainVector
146  << std::endl;
147  }
148 
149  queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
150 
151  double returnValue = 0.;
152  unsigned int cumulativeSize = 0;
153  for (unsigned int i = 0; i < m_densities.size(); ++i) {
154  V vec_i(m_densities[i]->domainSet().vectorSpace().zeroVector());
155  domainVector.cwExtract(cumulativeSize,vec_i);
156  double value_i = m_densities[i]->lnValue(vec_i,NULL,NULL,NULL,NULL);
157  returnValue += value_i;
158  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) { // gpmsa
159  *m_env.subDisplayFile() << "In ConcatenatedJointPdf<V,M>::lnValue()"
160  << ", vec_" << i << " = " << vec_i
161  << ": value_" << i << " = " << value_i
162  << ", temporary cumulative value = " << returnValue
163  << std::endl;
164  }
165  cumulativeSize += vec_i.sizeLocal();
166  }
167  //returnValue += m_logOfNormalizationFactor; // No need, because each PDF should be already normalized [PDF-11]
168 
169  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
170  *m_env.subDisplayFile() << "Leaving ConcatenatedJointPdf<V,M>::lnValue()"
171  << ": domainVector = " << domainVector
172  << ", returnValue = " << returnValue
173  << std::endl;
174  }
175 
176  return returnValue;
177 }
178 //--------------------------------------------------
179 template<class V, class M>
180 void
182 {
183  covMatrix.zeroLower();
184  covMatrix.zeroUpper();
185 
186  unsigned int cumulativeSize = 0;
187  for (unsigned int i = 0; i < m_densities.size(); ++i) {
188  const Map & map = m_densities[i]->domainSet().vectorSpace().map();
189  const unsigned int n_columns = map.NumGlobalElements();
190  M mat_i(m_densities[i]->domainSet().env(),
191  map, n_columns);
192  m_densities[i]->distributionVariance(mat_i);
193  covMatrix.cwSet(cumulativeSize,cumulativeSize,mat_i);
194 
195  cumulativeSize += n_columns;
196  }
197 }
198 
199 //--------------------------------------------------
200 template<class V, class M>
201 void
203 {
204  unsigned int cumulativeSize = 0;
205  for (unsigned int i = 0; i < m_densities.size(); ++i) {
206  V vec_i(m_densities[i]->domainSet().vectorSpace().zeroVector());
207  m_densities[i]->distributionMean(vec_i);
208  meanVector.cwSet(cumulativeSize,vec_i);
209  cumulativeSize += vec_i.sizeLocal();
210  }
211 }
212 
213 //--------------------------------------------------
214 template<class V, class M>
215 double
216 ConcatenatedJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
217 {
218  double value = 0.;
219 
220  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
221  *m_env.subDisplayFile() << "Entering ConcatenatedJointPdf<V,M>::computeLogOfNormalizationFactor()"
222  << std::endl;
223  }
224  double volume = m_domainSet.volume();
225  if ((queso_isnan(volume)) ||
226  (volume == -INFINITY ) ||
227  (volume == INFINITY ) ||
228  (volume <= 0. )) {
229  // Do nothing
230  }
231  else {
232  for (unsigned int i = 0; i < m_densities.size(); ++i) {
233  m_densities[i]->computeLogOfNormalizationFactor(numSamples, updateFactorInternally);
234  }
235  }
236  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
237  *m_env.subDisplayFile() << "Leaving ConcatenatedJointPdf<V,M>::computeLogOfNormalizationFactor()"
238  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
239  << std::endl;
240  }
241 
242  return value;
243 }
244 
245 } // End namespace QUESO
246 
A class for partitioning vectors and matrices.
Definition: Map.h:49
A class for handling concatenated PDFs.
bool queso_isnan(T arg)
Definition: math_macros.h:39
ConcatenatedJointPdf(const char *prefix, const BaseJointPdf< V, M > &density1, const BaseJointPdf< V, M > &density2, const VectorSet< V, M > &concatenatedDomain)
Constructor.
A templated class for handling sets.
Definition: VectorSet.h:52
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:55
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
int NumGlobalElements() const
Returns the total number of elements across all processors.
Definition: Map.C:85
std::vector< const BaseJointPdf< V, M > * > m_densities
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Calculates the actual values of each density.
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Calculates the logarithm of the values of each density.
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 setNormalizationStyle(unsigned int value) const
Sets the normalization style of all densities to value.
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.

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