queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
TransformedScaledCovMatrixTKGroup.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/TransformedScaledCovMatrixTKGroup.h>
26 #include <queso/InvLogitGaussianJointPdf.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 
30 namespace QUESO {
31 
32 template<class V, class M>
34  const char * prefix,
35  const BoxSubset<V,M> & boxSubset,
36  const std::vector<double> & scales,
37  const M & covMatrix)
38  : BaseTKGroup<V, M>(prefix, boxSubset.vectorSpace(), scales),
39  m_boxSubset(boxSubset),
40  m_originalCovMatrix(covMatrix)
41 {
42  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
43  *m_env.subDisplayFile() << "Entering TransformedScaledCovMatrixTKGroup<V,M>::constructor()"
44  << std::endl;
45  }
46 
47  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
48  *m_env.subDisplayFile() << "In TransformedScaledCovMatrixTKGroup<V,M>::constructor()"
49  << ": m_scales.size() = " << m_scales.size()
50  << ", m_preComputingPositions.size() = " << m_preComputingPositions.size()
51  << ", m_rvs.size() = " << m_rvs.size()
52  << ", m_originalCovMatrix = " << m_originalCovMatrix
53  << std::endl;
54  }
55 
56  // Set RVs to have zero mean in the Gaussian space
58 
59  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
60  *m_env.subDisplayFile() << "Leaving TransformedScaledCovMatrixTKGroup<V,M>::constructor()"
61  << std::endl;
62  }
63 }
64 // Destructor ---------------------------------------
65 template<class V, class M>
67 {
68 }
69 // Math/Stats methods--------------------------------
70 template<class V, class M>
71 bool
73 {
74  return false;
75 }
76 //---------------------------------------------------
77 template<class V, class M>
79 TransformedScaledCovMatrixTKGroup<V,M>::rv(unsigned int stageId) const
80 {
81  queso_require_not_equal_to_msg(m_rvs.size(), 0, "m_rvs.size() = 0");
82 
83  queso_require_msg(m_rvs[0], "m_rvs[0] == NULL");
84 
85  queso_require_greater_msg(m_preComputingPositions.size(), stageId, "m_preComputingPositions.size() <= stageId");
86 
87  queso_require_msg(m_preComputingPositions[stageId], "m_preComputingPositions[stageId] == NULL");
88 
89  if ((m_env.subDisplayFile() ) &&
90  (m_env.displayVerbosity() >= 10)) {
91  *m_env.subDisplayFile() << "In TransformedScaledCovMatrixTKGroup<V,M>::rv1()"
92  << ", stageId = " << stageId
93  << ": about to call m_rvs[0]->updateLawExpVector()"
94  << ", vector = " << *m_preComputingPositions[stageId] // FIX ME: might demand parallelism
95  << std::endl;
96  }
97 
98  InvLogitGaussianVectorRV<V, M> * invlogit_gaussian =
99  dynamic_cast<InvLogitGaussianVectorRV<V, M> * >(m_rvs[0]);
100 
101  V transformedPreComputingPositions(*m_preComputingPositions[stageId]);
102  transformToGaussianSpace(*m_preComputingPositions[stageId],
103  transformedPreComputingPositions);
104 
105  invlogit_gaussian->updateLawExpVector(transformedPreComputingPositions);
106 
107  return (*invlogit_gaussian);
108 }
109 //---------------------------------------------------
110 template<class V, class M>
112 TransformedScaledCovMatrixTKGroup<V,M>::rv(const std::vector<unsigned int>& stageIds)
113 {
114  queso_require_greater_equal_msg(m_rvs.size(), stageIds.size(), "m_rvs.size() < stageIds.size()");
115 
116  queso_require_msg(m_rvs[stageIds.size()-1], "m_rvs[stageIds.size()-1] == NULL");
117 
118  queso_require_greater_msg(m_preComputingPositions.size(), stageIds[0], "m_preComputingPositions.size() <= stageIds[0]");
119 
120  queso_require_msg(m_preComputingPositions[stageIds[0]], "m_preComputingPositions[stageIds[0]] == NULL");
121 
122  if ((m_env.subDisplayFile() ) &&
123  (m_env.displayVerbosity() >= 10)) {
124  *m_env.subDisplayFile() << "In TransformedScaledCovMatrixTKGroup<V,M>::rv2()"
125  << ", stageIds.size() = " << stageIds.size()
126  << ", stageIds[0] = " << stageIds[0]
127  << ": about to call m_rvs[stageIds.size()-1]->updateLawExpVector()"
128  << ", vector = " << *m_preComputingPositions[stageIds[0]] // FIX ME: might demand parallelism
129  << std::endl;
130  }
131 
132  InvLogitGaussianVectorRV<V, M> * invlogit_gaussian =
133  dynamic_cast<InvLogitGaussianVectorRV<V, M> * >(m_rvs[stageIds.size()-1]);
134 
135  V transformedPreComputingPositions(*m_preComputingPositions[stageIds[0]]);
136  transformToGaussianSpace(*m_preComputingPositions[stageIds[0]],
137  transformedPreComputingPositions);
138 
139  invlogit_gaussian->updateLawExpVector(transformedPreComputingPositions);
140 
141  return (*invlogit_gaussian);
142 }
143 
144 template <class V, class M>
147 {
148  queso_require_not_equal_to_msg(m_rvs.size(), 0, "m_rvs.size() = 0");
149  queso_require_msg(m_rvs[0], "m_rvs[0] == NULL");
150  // queso_require_greater_msg(m_preComputingPositions.size(), this->m_stageId, "m_preComputingPositions.size() <= stageId");
151  // queso_require_msg(m_preComputingPositions[this->m_stageId], "m_preComputingPositions[stageId] == NULL");
152 
153  InvLogitGaussianVectorRV<V, M> * invlogit_gaussian =
154  dynamic_cast<InvLogitGaussianVectorRV<V, M> * >(m_rvs[this->m_stageId]);
155 
156  V transformedPreComputingPositions(position);
157  transformToGaussianSpace(position, transformedPreComputingPositions);
158  invlogit_gaussian->updateLawExpVector(transformedPreComputingPositions);
159 
160  return (*invlogit_gaussian);
161 }
162 
163 //---------------------------------------------------
164 template<class V, class M>
165 void
167 {
168  for (unsigned int i = 0; i < m_scales.size(); ++i) {
169  double factor = 1./m_scales[i]/m_scales[i];
170  if ((m_env.subDisplayFile() ) &&
171  (m_env.displayVerbosity() >= 10)) {
172  *m_env.subDisplayFile() << "In TransformedScaledCovMatrixTKGroup<V,M>::updateLawCovMatrix()"
173  << ", m_scales.size() = " << m_scales.size()
174  << ", i = " << i
175  << ", m_scales[i] = " << m_scales[i]
176  << ", factor = " << factor
177  << ": about to call m_rvs[i]->updateLawCovMatrix()"
178  << ", covMatrix = \n" << factor*covMatrix // FIX ME: might demand parallelism
179  << std::endl;
180  }
181 
182  InvLogitGaussianVectorRV<V, M> * invlogit_gaussian =
183  dynamic_cast<InvLogitGaussianVectorRV<V, M> * >(m_rvs[i]);
184 
185  invlogit_gaussian->updateLawCovMatrix(factor*covMatrix);
186  }
187 
188  return;
189 }
190 
191 // Misc methods -------------------------------------
192 template<class V, class M>
193 bool
195 {
196  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
197  *m_env.subDisplayFile() << "Entering TransformedScaledCovMatrixTKGroup<V,M>::setPreComputingPosition()"
198  << ": position = " << position
199  << ", stageId = " << stageId
200  << std::endl;
201  }
202 
204  //setRVsWithZeroMean();
205 
206  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
207  *m_env.subDisplayFile() << "In TransformedScaledCovMatrixTKGroup<V,M>::setPreComputingPosition()"
208  << ", position = " << position
209  << ", stageId = " << stageId
210  << ": preComputingPos = " << *m_preComputingPositions[stageId];
211  if (stageId < m_scales.size()) {
212  *m_env.subDisplayFile() << ", factor = " << 1./m_scales[stageId]/m_scales[stageId];
213  }
214  if (stageId < m_rvs.size()) {
215  const InvLogitGaussianJointPdf<V,M>* pdfPtr = dynamic_cast< const InvLogitGaussianJointPdf<V,M>* >(&(m_rvs[stageId]->pdf()));
216  *m_env.subDisplayFile() << ", rvCov = " << pdfPtr->lawCovMatrix(); // FIX ME: might demand parallelism
217  }
218  *m_env.subDisplayFile() << std::endl;
219  }
220 
221  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
222  *m_env.subDisplayFile() << "Leaving TransformedScaledCovMatrixTKGroup<V,M>::setPreComputingPosition()"
223  << ": position = " << position
224  << ", stageId = " << stageId
225  << std::endl;
226  }
227 
228  return true;
229 }
230 //---------------------------------------------------
231 template<class V, class M>
232 void
234 {
236  return;
237 }
238 
239 template <class V, class M>
240 unsigned int
242 {
243  unsigned int old_stageId = this->m_stageId;
244  this->m_stageId = stageId;
245  return old_stageId;
246 }
247 
248 // Private methods------------------------------------
249 template<class V, class M>
250 void
252 {
253  queso_require_not_equal_to_msg(m_rvs.size(), 0, "m_rvs.size() = 0");
254 
255  queso_require_equal_to_msg(m_rvs.size(), m_scales.size(), "m_rvs.size() != m_scales.size()");
256 
257  for (unsigned int i = 0; i < m_scales.size(); ++i) {
258  double factor = 1./m_scales[i]/m_scales[i];
259  queso_require_msg(!(m_rvs[i]), "m_rvs[i] != NULL");
260  m_rvs[i] = new InvLogitGaussianVectorRV<V,M>(m_prefix.c_str(),
261  m_boxSubset, m_vectorSpace->zeroVector(),
262  factor*m_originalCovMatrix);
263  }
264 
265  return;
266 }
267 
268 template<class V, class M>
269 void
271 {
273  return;
274 }
275 
276 template<class V, class M>
277 void
279  const V & physicalPoint, V & transformedPoint) const
280 {
281  V min_domain_bounds(this->m_boxSubset.minValues());
282  V max_domain_bounds(this->m_boxSubset.maxValues());
283 
284  for (unsigned int i = 0; i < transformedPoint.sizeLocal(); i++) {
285  double min_val = min_domain_bounds[i];
286  double max_val = max_domain_bounds[i];
287 
288  if (queso_isfinite(min_val) &&
289  queso_isfinite(max_val)) {
290  // Left- and right-hand sides are finite. Do full transform.
291  transformedPoint[i] = std::log(physicalPoint[i] - min_val) -
292  std::log(max_val - physicalPoint[i]);
293  }
294  else if (queso_isfinite(min_val) &&
295  !queso_isfinite(max_val)) {
296  // Left-hand side finite, but right-hand side is not.
297  // Do only left-hand transform.
298  transformedPoint[i] = std::log(physicalPoint[i] - min_val);
299  }
300  else if (!queso_isfinite(min_val) &&
301  queso_isfinite(max_val)) {
302  // Right-hand side is finite, but left-hand side is not.
303  // Do only right-hand transform.
304  transformedPoint[i] = -std::log(max_val - physicalPoint[i]);
305  }
306  else {
307  // No transform.
308  transformedPoint[i] = physicalPoint[i];
309  }
310  }
311 }
312 
313 } // End namespace QUESO
314 
A class representing a (transformed) Gaussian vector RV with bounds.
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
std::vector< const V * > m_preComputingPositions
Definition: TKGroup.h:132
virtual unsigned int set_dr_stage(unsigned int stageId)
Does nothing. Subclasses may re-implement. Returns the current stage id.
bool queso_isfinite(T arg)
Definition: math_macros.h:51
std::vector< BaseVectorRV< V, M > * > m_rvs
Definition: TKGroup.h:133
const BaseEnvironment & m_env
Definition: TKGroup.h:128
virtual void print(std::ostream &os) const
TODO: Prints the transition kernel.
Definition: TKGroup.C:133
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the covariance matrix.
void updateLawCovMatrix(const M &covMatrix)
Scales the covariance matrix of the underlying Gaussian distribution.
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
virtual void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
Definition: TKGroup.C:111
void print(std::ostream &os) const
TODO: Prints the transition kernel.
This class represents a transition kernel with a scaled covariance matrix on hybrid bounded/unbounded...
bool symmetric() const
Whether or not the matrix is symmetric. Always &#39;false&#39;.
bool setPreComputingPosition(const V &position, unsigned int stageId)
Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position...
Class representing a subset of a vector space shaped like a hypercube.
Definition: BoxSubset.h:44
const InvLogitGaussianVectorRV< V, M > & rv(unsigned int stageId) const
InvLogitGaussian increment property to construct a transition kernel.
void setRVsWithZeroMean()
Sets the mean of the underlying Gaussian RVs to zero.
std::vector< double > m_scales
Definition: TKGroup.h:131
A class for handling hybrid (transformed) Gaussians with bounds.
This base class allows the representation of a transition kernel.
Definition: Algorithm.h:32
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"))
TransformedScaledCovMatrixTKGroup(const char *prefix, const BoxSubset< V, M > &boxSubset, const std::vector< double > &scales, const M &covMatrix)
Default constructor.
unsigned int displayVerbosity() const
Definition: Environment.C:450
virtual bool setPreComputingPosition(const V &position, unsigned int stageId)
Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position...
Definition: TKGroup.C:98
void transformToGaussianSpace(const V &physicalPoint, V &transformedPoint) const
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:56 for queso-0.57.1 by  doxygen 1.8.5