queso-0.50.1
HessianCovMatricesTKGroup.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,2009,2010,2011,2012,2013 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/HessianCovMatricesTKGroup.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Default constructor ------------------------------
32 template<class V, class M>
33 HessianCovMatricesTKGroup<V,M>::HessianCovMatricesTKGroup(
34  const char* prefix,
35  const VectorSpace<V,M>& vectorSpace,
36  const std::vector<double>& scales,
37  const ScalarFunctionSynchronizer<V,M>& targetPdfSynchronizer)
38  :
39  BaseTKGroup<V,M>(prefix,vectorSpace,scales),
40  m_targetPdfSynchronizer(targetPdfSynchronizer),
41  m_originalNewtonSteps (scales.size()+1,NULL), // Yes, +1
42  m_originalCovMatrices (scales.size()+1,NULL) // Yes, +1
43 {
44  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
45  *m_env.subDisplayFile() << "Entering HessianCovMatricesTKGroup<V,M>::constructor()"
46  << std::endl;
47  }
48 
49  m_rvs.resize(scales.size()+1,NULL); // Yes, +1 (IMPORTANT: overwrite initialization done by BaseTKGroup<V,M>)
50 
51  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
52  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::constructor()"
53  << ": m_scales.size() = " << m_scales.size()
54  << ", m_preComputingPositions.size() = " << m_preComputingPositions.size()
55  << ", m_rvs.size() = " << m_rvs.size()
56  << ", m_originalNewtonSteps.size() = " << m_originalNewtonSteps.size()
57  << ", m_originalCovMatrices.size() = " << m_originalCovMatrices.size()
58  << std::endl;
59  }
60 
61  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
62  *m_env.subDisplayFile() << "Leaving HessianCovMatricesTKGroup<V,M>::constructor()"
63  << std::endl;
64  }
65 }
66 // Destructor ---------------------------------------
67 template<class V, class M>
68 HessianCovMatricesTKGroup<V,M>::~HessianCovMatricesTKGroup()
69 {
70 }
71 // Math/Stats methods--------------------------------
72 template<class V, class M>
73 bool
74 HessianCovMatricesTKGroup<V,M>::symmetric() const
75 {
76  return false;
77 }
78 //---------------------------------------------------
79 template<class V, class M>
80 const GaussianVectorRV<V,M>&
81 HessianCovMatricesTKGroup<V,M>::rv(unsigned int stageId)
82 {
83  UQ_FATAL_TEST_MACRO(m_rvs.size() <= stageId,
84  m_env.worldRank(),
85  "HessianCovMatricesTKGroup<V,M>::rv1()",
86  "m_rvs.size() <= stageId");
87 
88  UQ_FATAL_TEST_MACRO(m_rvs[stageId] == NULL,
89  m_env.worldRank(),
90  "HessianCovMatricesTKGroup<V,M>::rv1()",
91  "m_rvs[stageId] == NULL");
92 
93  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() <= stageId,
94  m_env.worldRank(),
95  "HessianCovMatricesTKGroup<V,M>::rv1()",
96  "m_preComputingPositions.size() <= stageId");
97 
98  UQ_FATAL_TEST_MACRO(m_preComputingPositions[stageId] == NULL,
99  m_env.worldRank(),
100  "HessianCovMatricesTKGroup<V,M>::rv1()",
101  "m_preComputingPositions[stageId] == NULL");
102 
103  m_rvs[stageId]->updateLawExpVector(*m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]);
104  if ((m_env.subDisplayFile() ) &&
105  (m_env.displayVerbosity() >= 10)) {
106  *m_env.subDisplayFile() << "In HessianCovMatrixTKGroup<V,M>::rv1()"
107  << ", stageId = " << stageId
108  << ": about to call m_rvs[stageId]->updateLawCovMatrix()"
109  << ", covMatrix = \n" << *m_originalCovMatrices[stageId] // FIX ME: might demand parallelism
110  << std::endl;
111  }
112  m_rvs[stageId]->updateLawCovMatrix(*m_originalCovMatrices[stageId]);
113 
114  return *m_rvs[stageId];
115 }
116 //---------------------------------------------------
117 template<class V, class M>
118 const GaussianVectorRV<V,M>&
119 HessianCovMatricesTKGroup<V,M>::rv(const std::vector<unsigned int>& stageIds)
120 {
121  UQ_FATAL_TEST_MACRO(m_rvs.size() <= stageIds[0],
122  m_env.worldRank(),
123  "HessianCovMatricesTKGroup<V,M>::rv2()",
124  "m_rvs.size() <= stageIds[0]");
125 
126  UQ_FATAL_TEST_MACRO(m_rvs[stageIds[0]] == NULL,
127  m_env.worldRank(),
128  "HessianCovMatricesTKGroup<V,M>::rv2()",
129  "m_rvs[stageIds[0]] == NULL");
130 
131  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() <= stageIds[0],
132  m_env.worldRank(),
133  "HessianCovMatricesTKGroup<V,M>::rv2()",
134  "m_preComputingPositions.size() <= stageIds[0]");
135 
136  UQ_FATAL_TEST_MACRO(m_preComputingPositions[stageIds[0]] == NULL,
137  m_env.worldRank(),
138  "HessianCovMatricesTKGroup<V,M>::rv2()",
139  "m_preComputingPositions[stageIds[0]] == NULL");
140 
141  double factor = 1./m_scales[stageIds.size()-1]/m_scales[stageIds.size()-1];
142 
143  m_rvs[stageIds[0]]->updateLawExpVector(*m_preComputingPositions[stageIds[0]] + factor*(*m_originalNewtonSteps[stageIds[0]]));
144  if ((m_env.subDisplayFile() ) &&
145  (m_env.displayVerbosity() >= 10)) {
146  *m_env.subDisplayFile() << "In HessianCovMatrixTKGroup<V,M>::rv2()"
147  << ", stageIds.size() = " << stageIds.size()
148  << ", stageIds[0] = " << stageIds[0]
149  << ", factor = " << factor
150  << ": about to call m_rvs[stageIds[0]]->updateLawCovVector()"
151  << ", covMatrix = \n" << factor*(*m_originalCovMatrices[stageIds[0]]) // FIX ME: might demand parallelism
152  << std::endl;
153  }
154  m_rvs[stageIds[0]]->updateLawCovMatrix(factor*(*m_originalCovMatrices[stageIds[0]]));
155 
156  return *m_rvs[stageIds[0]];
157 }
158 // Misc methods--------------------------------------
159 template<class V, class M>
160 bool
161 HessianCovMatricesTKGroup<V,M>::setPreComputingPosition(const V& position, unsigned int stageId)
162 {
163  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
164  *m_env.subDisplayFile() << "Entering HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
165  << ": position = " << position
166  << ", stageId = " << stageId
167  << std::endl;
168  }
169 
170  bool validPreComputingPosition = true;
171 
172  // Verify consistency of sizes
173  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() <= stageId,
174  m_env.worldRank(),
175  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
176  "m_preComputingPositions.size() <= stageId");
177 
178  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() != m_rvs.size(),
179  m_env.worldRank(),
180  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
181  "m_preComputingPositions.size() != m_rvs.size()");
182 
183  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() != m_originalNewtonSteps.size(),
184  m_env.worldRank(),
185  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
186  "m_preComputingPositions.size() != m_originalNewtonSteps.size()");
187 
188  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() != m_originalCovMatrices.size(),
189  m_env.worldRank(),
190  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
191  "m_preComputingPositions.size() != m_originalCovMatrices.size()");
192 
193  // Verify data is not null
194  UQ_FATAL_TEST_MACRO(m_preComputingPositions[stageId] != NULL,
195  m_env.worldRank(),
196  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
197  "m_preComputingPositions[stageId] != NULL");
198 
199  UQ_FATAL_TEST_MACRO(m_rvs[stageId] != NULL,
200  m_env.worldRank(),
201  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
202  "m_rvs[stageId] != NULL");
203 
204  UQ_FATAL_TEST_MACRO(m_originalNewtonSteps[stageId] != NULL,
205  m_env.worldRank(),
206  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
207  "m_originalNewtonSteps[stageId] != NULL");
208 
209  UQ_FATAL_TEST_MACRO(m_originalCovMatrices[stageId] != NULL,
210  m_env.worldRank(),
211  "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
212  "m_originalCovMatrices[stageId] != NULL");
213 
214  BaseTKGroup<V,M>::setPreComputingPosition(position,stageId);
215 
216  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
217  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
218  << ", position = " << position
219  << ", stageId = " << stageId
220  << ": m_originalNewtonSteps.size() = " << m_originalNewtonSteps.size()
221  << ", m_originalCovMatrices.size() = " << m_originalCovMatrices.size()
222  << ", m_preComputingPositions.size() = " << m_preComputingPositions.size()
223  << ", m_rvs.size() = " << m_rvs.size()
224  << std::endl;
225  }
226 
227  if (m_targetPdfSynchronizer.domainSet().contains(position)) {
228  M* tmpHessian = m_vectorSpace->newMatrix();
229  M* tmpCovMat = m_vectorSpace->newMatrix();
230  V* tmpGrad = m_vectorSpace->newVector();
231 
232  double logPrior = 0.;
233  double logLikelihood = 0.;
234  double logTarget = 0.;
235  logTarget = m_targetPdfSynchronizer.callFunction(&position, // Might demand parallel environment
236  NULL,
237  tmpGrad,
238  tmpHessian,
239  NULL,
240  &logPrior,
241  &logLikelihood);
242  if (logTarget) {}; // just to remove compiler warning
243 
244  // IMPORTANT: covariance matrix = (Hessian)^{-1} !!!
245  V unitVector(m_vectorSpace->zeroVector());
246  V multVector(m_vectorSpace->zeroVector());
247  for (unsigned int j = 0; j < tmpHessian->numCols(); ++j) {
248  if (j > 0) unitVector[j-1] = 0.;
249  unitVector[j] = 1.;
250  tmpHessian->invertMultiply(unitVector, multVector);
251  for (unsigned int i = 0; i < tmpHessian->numRowsLocal(); ++i) {
252  (*tmpCovMat)(i,j) = multVector[i];
253  }
254  }
255  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
256  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
257  << ", position = " << position
258  << ", stageId = " << stageId
259  << ":\n H = " << *tmpHessian
260  << "\n H^{-1} = " << *tmpCovMat
261  << "\n H*H^{-1} = " << (*tmpHessian)*(*tmpCovMat)
262  << "\n H^{-1}*H = " << (*tmpCovMat)*(*tmpHessian)
263  << std::endl;
264  }
265 
266  // Force covariance matrix to be symmetric, as the Hessian (supposedly) is
267  *tmpCovMat = .5*((*tmpCovMat) + tmpCovMat->transpose());
268 
269  // Test if covariance matrix is positive definite
270  M lowerChol(*tmpCovMat);
271  if ((m_env.subDisplayFile() ) &&
272  (m_env.displayVerbosity() >= 10)) {
273  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
274  << ", position = " << position
275  << ", stageId = " << stageId
276  << ": calling lowerChol.chol()"
277  << ", lowerChol = " << lowerChol
278  << std::endl;
279  }
280  int iRC = lowerChol.chol();
281  if (iRC) {
282  std::cerr << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition(): chol failed\n";
283  }
284  if ((m_env.subDisplayFile() ) &&
285  (m_env.displayVerbosity() >= 10)) {
286  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
287  << ", position = " << position
288  << ", stageId = " << stageId
289  << ": got lowerChol.chol() with iRC = " << iRC
290  << std::endl;
291  }
292 
293  bool covIsPositiveDefinite = !iRC;
294 
295  if (covIsPositiveDefinite) {
296  //UQ_FATAL_TEST_MACRO(stageId >= m_scales.size(),
297  // m_env.worldRank(),
298  // "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
299  // "stageId is too large for m_scales");
300  //double factor = 1./m_scales[stageId]/m_scales[stageId];
301  //*tmpCovMat *= factor;
302 
303  m_originalNewtonSteps[stageId] = new V(-1.*(*tmpCovMat)*(*tmpGrad));
304  m_originalCovMatrices[stageId] = new M(*tmpCovMat);
305 
306  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
307  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
308  << ", position = " << position
309  << ", stageId = " << stageId
310  << ", about to instantiate a Gaussian RV"
311  << ": tmpHessian = " << *tmpHessian
312  << ", preComputingPos = " << *m_preComputingPositions[stageId]
313  << ", tmpCovMat = " << *tmpCovMat
314  << ", tmpGrad = " << *tmpGrad
315  << ", preComputedPos = " << *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]
316  //<< ", factor = " << factor
317  //<< ", rvCov = " << factor*(*tmpCovMat)
318  << std::endl;
319  }
320  m_rvs[stageId] = new GaussianVectorRV<V,M>(m_prefix.c_str(),
321  *m_vectorSpace,
322  *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId],
323  *m_originalCovMatrices[stageId]);
324  }
325  else {
326  validPreComputingPosition = false;
327  }
328 
329  delete tmpGrad;
330  delete tmpCovMat;
331  delete tmpHessian;
332  }
333  else {
334  validPreComputingPosition = false;
335  }
336 
337  if (validPreComputingPosition == false) {
338  // Put "default" values on variables
339  V tmpGrad (m_vectorSpace->zeroVector());
340  M tmpCovMat(tmpGrad,1.); // = identity matrix
341  m_originalNewtonSteps[stageId] = new V(-1.*tmpCovMat*tmpGrad);
342  m_originalCovMatrices[stageId] = new M(tmpCovMat);
343  m_rvs[stageId] = new GaussianVectorRV<V,M>(m_prefix.c_str(),
344  *m_vectorSpace,
345  *m_preComputingPositions[stageId],
346  tmpCovMat);
347  }
348 
349  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
350  *m_env.subDisplayFile() << "Leaving HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
351  << ": position = " << position
352  << ", stageId = " << stageId
353  << std::endl;
354  }
355 
356  return validPreComputingPosition;
357 }
358 //---------------------------------------------------
359 template<class V, class M>
360 void
361 HessianCovMatricesTKGroup<V,M>::clearPreComputingPositions()
362 {
363  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() != m_originalNewtonSteps.size(),
364  m_env.worldRank(),
365  "HessianCovMatricesTKGroup<V,M>::clearPreComputingPositions()",
366  "m_preComputingPositions.size() != m_originalNewtonSteps.size()");
367 
368  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() != m_originalCovMatrices.size(),
369  m_env.worldRank(),
370  "HessianCovMatricesTKGroup<V,M>::clearPreComputingPositions()",
371  "m_preComputingPositions.size() != m_originalCovMatrices.size()");
372 
373  BaseTKGroup<V,M>::clearPreComputingPositions();
374 
375  // RVs are not deleted in base class because the cov matrices are constant in the case of scaledTK class
376  for (unsigned int i = 0; i < m_rvs.size(); ++i) {
377  if (m_rvs[i]) {
378  delete m_rvs[i];
379  m_rvs[i] = NULL;
380  }
381  }
382 
383  for (unsigned int i = 0; i < m_originalNewtonSteps.size(); ++i) {
384  if (m_originalNewtonSteps[i]) {
385  delete m_originalNewtonSteps[i];
386  m_originalNewtonSteps[i] = NULL;
387  }
388  }
389 
390  for (unsigned int i = 0; i < m_originalCovMatrices.size(); ++i) {
391  if (m_originalCovMatrices[i]) {
392  delete m_originalCovMatrices[i];
393  m_originalCovMatrices[i] = NULL;
394  }
395  }
396 
397  return;
398 }
399 // I/O methods---------------------------------------
400 template<class V, class M>
401 void
402 HessianCovMatricesTKGroup<V,M>::print(std::ostream& os) const
403 {
404  BaseTKGroup<V,M>::print(os);
405  return;
406 }
407 
408 } // End namespace QUESO
409 
410 template class QUESO::HessianCovMatricesTKGroup<QUESO::GslVector, QUESO::GslMatrix>;
QUESO::HessianCovMatricesTKGroup::clearPreComputingPositions
void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
Definition: HessianCovMatricesTKGroup.C:361
QUESO::BaseTKGroup::clearPreComputingPositions
virtual void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
Definition: TKGroup.C:121
QUESO::VectorSpace
A class representing a vector space.
Definition: VectorSet.h:46
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::HessianCovMatricesTKGroup::print
void print(std::ostream &os) const
TODO: Prints the transition kernel.
Definition: HessianCovMatricesTKGroup.C:402
QUESO::BaseTKGroup::m_preComputingPositions
std::vector< const V * > m_preComputingPositions
Definition: TKGroup.h:102
QUESO::HessianCovMatricesTKGroup::rv
const GaussianVectorRV< V, M > & rv(unsigned int stageId)
Gaussian increment property to construct a transition kernel.
Definition: HessianCovMatricesTKGroup.C:81
QUESO::BaseTKGroup::setPreComputingPosition
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:102
QUESO::BaseTKGroup::print
virtual void print(std::ostream &os) const
TODO: Prints the transition kernel.
Definition: TKGroup.C:135
QUESO::BaseTKGroup::m_rvs
std::vector< GaussianVectorRV< V, M > * > m_rvs
Definition: TKGroup.h:103
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::HessianCovMatricesTKGroup::setPreComputingPosition
bool setPreComputingPosition(const V &position, unsigned int stageId)
Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position...
Definition: HessianCovMatricesTKGroup.C:161
QUESO::BaseTKGroup
This base class allows the representation of a transition kernel.
Definition: TKGroup.h:48
QUESO::HessianCovMatricesTKGroup::symmetric
bool symmetric() const
Whether or not the matrix is symmetric. Always &#39;false&#39;.
Definition: HessianCovMatricesTKGroup.C:74
QUESO::HessianCovMatricesTKGroup
This class allows the representation of a transition kernel with Hessians.
Definition: HessianCovMatricesTKGroup.h:41
QUESO::GaussianVectorRV
A class representing a Gaussian vector RV.
Definition: GaussianVectorRV.h:60
QUESO::HessianCovMatricesTKGroup::HessianCovMatricesTKGroup
HessianCovMatricesTKGroup(const char *prefix, const VectorSpace< V, M > &vectorSpace, const std::vector< double > &scales, const ScalarFunctionSynchronizer< V, M > &targetPdfSynchronizer)
Default constructor.
Definition: HessianCovMatricesTKGroup.C:33
QUESO::ScalarFunctionSynchronizer
A templated class for synchronizing the calls of scalar functions (BaseScalarFunction and derived cla...
Definition: ScalarFunctionSynchronizer.h:44
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:436
QUESO::BaseTKGroup::m_scales
std::vector< double > m_scales
Definition: TKGroup.h:101
QUESO::BaseTKGroup::m_env
const BaseEnvironment & m_env
Definition: TKGroup.h:98
Generated on Thu Apr 23 2015 19:18:34 for queso-0.50.1 by  doxygen 1.8.5