queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
WignerJointPdf.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/WignerJointPdf.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 VectorSet<V,M>& domainSet,
36  const V& centerPos,
37  double radius)
38  :
39  BaseJointPdf<V,M>(((std::string)(prefix)+"uni").c_str(),
40  domainSet),
41  m_centerPos(new V(centerPos)),
42  m_radius (radius)
43 {
44  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
45  *m_env.subDisplayFile() << "Entering WignerJointPdf<V,M>::constructor()"
46  << ": prefix = " << m_prefix
47  << std::endl;
48  }
49 
50  queso_require_greater_msg(m_radius, 0., "invalid radius");
51 
52  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
53  *m_env.subDisplayFile() << "Leaving WignerJointPdf<V,M>::constructor()"
54  << ": prefix = " << m_prefix
55  << std::endl;
56  }
57 }
58 // Destructor --------------------------------------
59 template<class V,class M>
61 {
62  delete m_centerPos;
63 }
64 // Math methods-------------------------------------
65 template<class V, class M>
66 double
68  const V& domainVector,
69  const V* domainDirection,
70  V* gradVector,
71  M* hessianMatrix,
72  V* hessianEffect) const
73 {
74  queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
75 
76  if (gradVector ) *gradVector = m_domainSet.vectorSpace().zeroVector();
77  if (hessianMatrix) *hessianMatrix *= 0.;
78  if (hessianEffect) *hessianEffect = m_domainSet.vectorSpace().zeroVector();
79 
80  double returnValue = 0.;
81  double distanceRatio = (domainVector - *m_centerPos).norm2()/m_radius;
82  if (distanceRatio < 1.) {
83  returnValue = 2.*m_radius*m_radius*sqrt(1. - distanceRatio*distanceRatio)/M_PI;
84  }
85  returnValue *= exp(m_logOfNormalizationFactor); // [PDF-09]
86 
87  return returnValue;
88 }
89 //--------------------------------------------------
90 template<class V, class M>
91 double
93  const V& domainVector,
94  const V* domainDirection,
95  V* gradVector,
96  M* hessianMatrix,
97  V* hessianEffect) const
98 {
99  if (gradVector ) *gradVector = m_domainSet.vectorSpace().zeroVector();
100  if (hessianMatrix) *hessianMatrix *= 0.;
101  if (hessianEffect) *hessianEffect = m_domainSet.vectorSpace().zeroVector();
102 
103  // No need to add m_logOfNormalizationFactor because 'actualValue()' is called [PDF-09]
104  return log(this->actualValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
105 }
106 //--------------------------------------------------
107 template<class V, class M>
108 void
110 {
111  meanVector = *m_centerPos;
112 }
113 
114 //---------------------------------------------------
115 template<class V,class M>
116 void
118 {
119  queso_not_implemented();
120 }
121 
122 
123 //--------------------------------------------------
124 template<class V, class M>
125 double
126 WignerJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
127 {
128  double value = 0.;
129 
130  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
131  *m_env.subDisplayFile() << "Entering WignerJointPdf<V,M>::computeLogOfNormalizationFactor()"
132  << std::endl;
133  }
134  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
135  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
136  *m_env.subDisplayFile() << "Leaving WignerJointPdf<V,M>::computeLogOfNormalizationFactor()"
137  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
138  << std::endl;
139  }
140 
141  return value;
142 }
143 
144 } // End namespace QUESO
145 
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the PDF (scalar function).
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
~WignerJointPdf()
Destructor.
A templated class for handling sets.
Definition: VectorSet.h:52
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
const BaseEnvironment & m_env
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
A class for handling Wigner joint PDFs.
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:55
WignerJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &centerPos, double radius)
Constructor.
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Computes the logarithm of the value of the function.
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"))
unsigned int displayVerbosity() const
Definition: Environment.C:450
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:56 for queso-0.57.1 by  doxygen 1.8.5