queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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-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/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>
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>
69 {
70 }
71 // Math/Stats methods--------------------------------
72 template<class V, class M>
73 bool
75 {
76  return false;
77 }
78 //---------------------------------------------------
79 template<class V, class M>
81 HessianCovMatricesTKGroup<V,M>::rv(unsigned int stageId) const
82 {
83  queso_require_greater_msg(m_rvs.size(), stageId, "m_rvs.size() <= stageId");
84 
85  queso_require_msg(m_rvs[stageId], "m_rvs[stageId] == NULL");
86 
87  queso_require_greater_msg(m_preComputingPositions.size(), stageId, "m_preComputingPositions.size() <= stageId");
88 
89  queso_require_msg(m_preComputingPositions[stageId], "m_preComputingPositions[stageId] == NULL");
90 
91  GaussianVectorRV<V, M> * gaussian_rv =
92  dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[stageId]);
93 
94  gaussian_rv->updateLawExpVector(*m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]);
95 
96  if ((m_env.subDisplayFile() ) &&
97  (m_env.displayVerbosity() >= 10)) {
98  *m_env.subDisplayFile() << "In HessianCovMatrixTKGroup<V,M>::rv1()"
99  << ", stageId = " << stageId
100  << ": about to call m_rvs[stageId]->updateLawCovMatrix()"
101  << ", covMatrix = \n" << *m_originalCovMatrices[stageId] // FIX ME: might demand parallelism
102  << std::endl;
103  }
104 
105  gaussian_rv->updateLawCovMatrix(*m_originalCovMatrices[stageId]);
106 
107  return *gaussian_rv;
108 }
109 //---------------------------------------------------
110 template<class V, class M>
112 HessianCovMatricesTKGroup<V,M>::rv(const std::vector<unsigned int>& stageIds)
113 {
114  queso_require_greater_msg(m_rvs.size(), stageIds[0], "m_rvs.size() <= stageIds[0]");
115 
116  queso_require_msg(m_rvs[stageIds[0]], "m_rvs[stageIds[0]] == 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  double factor = 1./m_scales[stageIds.size()-1]/m_scales[stageIds.size()-1];
123 
124  GaussianVectorRV<V, M> * gaussian_rv =
125  dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[stageIds[0]]);
126 
127  gaussian_rv->updateLawExpVector(*m_preComputingPositions[stageIds[0]] + factor*(*m_originalNewtonSteps[stageIds[0]]));
128 
129  if ((m_env.subDisplayFile() ) &&
130  (m_env.displayVerbosity() >= 10)) {
131  *m_env.subDisplayFile() << "In HessianCovMatrixTKGroup<V,M>::rv2()"
132  << ", stageIds.size() = " << stageIds.size()
133  << ", stageIds[0] = " << stageIds[0]
134  << ", factor = " << factor
135  << ": about to call m_rvs[stageIds[0]]->updateLawCovVector()"
136  << ", covMatrix = \n" << factor*(*m_originalCovMatrices[stageIds[0]]) // FIX ME: might demand parallelism
137  << std::endl;
138  }
139  gaussian_rv->updateLawCovMatrix(factor*(*m_originalCovMatrices[stageIds[0]]));
140 
141  return *gaussian_rv;
142 }
143 
144 template <class V, class M>
146 HessianCovMatricesTKGroup<V, M>::rv(const V & position) const
147 {
148  queso_require_greater_msg(m_rvs.size(), this->m_stageId, "m_rvs.size() <= stageId");
149  queso_require_msg(m_rvs[this->m_stageId], "m_rvs[stageId] == 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  GaussianVectorRV<V, M> * gaussian_rv =
154  dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[this->m_stageId]);
155 
156  gaussian_rv->updateLawExpVector(position + *m_originalNewtonSteps[this->m_stageId]);
157 
158  gaussian_rv->updateLawCovMatrix(*m_originalCovMatrices[this->m_stageId]);
159 
160  return *gaussian_rv;
161 }
162 
163 // Misc methods--------------------------------------
164 template<class V, class M>
165 bool
166 HessianCovMatricesTKGroup<V,M>::setPreComputingPosition(const V& position, unsigned int stageId)
167 {
168  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
169  *m_env.subDisplayFile() << "Entering HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
170  << ": position = " << position
171  << ", stageId = " << stageId
172  << std::endl;
173  }
174 
175  bool validPreComputingPosition = true;
176 
177  // Verify consistency of sizes
178  queso_require_greater_msg(m_preComputingPositions.size(), stageId, "m_preComputingPositions.size() <= stageId");
179 
180  queso_require_equal_to_msg(m_preComputingPositions.size(), m_rvs.size(), "m_preComputingPositions.size() != m_rvs.size()");
181 
182  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalNewtonSteps.size(), "m_preComputingPositions.size() != m_originalNewtonSteps.size()");
183 
184  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalCovMatrices.size(), "m_preComputingPositions.size() != m_originalCovMatrices.size()");
185 
186  // Verify data is not null
187  queso_require_msg(!(m_preComputingPositions[stageId]), "m_preComputingPositions[stageId] != NULL");
188 
189  queso_require_msg(!(m_rvs[stageId]), "m_rvs[stageId] != NULL");
190 
191  queso_require_msg(!(m_originalNewtonSteps[stageId]), "m_originalNewtonSteps[stageId] != NULL");
192 
193  queso_require_msg(!(m_originalCovMatrices[stageId]), "m_originalCovMatrices[stageId] != NULL");
194 
196 
197  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
198  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
199  << ", position = " << position
200  << ", stageId = " << stageId
201  << ": m_originalNewtonSteps.size() = " << m_originalNewtonSteps.size()
202  << ", m_originalCovMatrices.size() = " << m_originalCovMatrices.size()
203  << ", m_preComputingPositions.size() = " << m_preComputingPositions.size()
204  << ", m_rvs.size() = " << m_rvs.size()
205  << std::endl;
206  }
207 
208  if (m_targetPdfSynchronizer.domainSet().contains(position)) {
209  M* tmpHessian = m_vectorSpace->newMatrix();
210  M* tmpCovMat = m_vectorSpace->newMatrix();
211  V* tmpGrad = m_vectorSpace->newVector();
212 
213  double logPrior = 0.;
214  double logLikelihood = 0.;
215  double logTarget = 0.;
216  logTarget = m_targetPdfSynchronizer.callFunction(&position, // Might demand parallel environment
217  NULL,
218  tmpGrad,
219  tmpHessian,
220  NULL,
221  &logPrior,
222  &logLikelihood);
223  if (logTarget) {}; // just to remove compiler warning
224 
225  // IMPORTANT: covariance matrix = (Hessian)^{-1} !!!
226  V unitVector(m_vectorSpace->zeroVector());
227  V multVector(m_vectorSpace->zeroVector());
228  for (unsigned int j = 0; j < tmpHessian->numCols(); ++j) {
229  if (j > 0) unitVector[j-1] = 0.;
230  unitVector[j] = 1.;
231  tmpHessian->invertMultiply(unitVector, multVector);
232  for (unsigned int i = 0; i < tmpHessian->numRowsLocal(); ++i) {
233  (*tmpCovMat)(i,j) = multVector[i];
234  }
235  }
236  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
237  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
238  << ", position = " << position
239  << ", stageId = " << stageId
240  << ":\n H = " << *tmpHessian
241  << "\n H^{-1} = " << *tmpCovMat
242  << "\n H*H^{-1} = " << (*tmpHessian)*(*tmpCovMat)
243  << "\n H^{-1}*H = " << (*tmpCovMat)*(*tmpHessian)
244  << std::endl;
245  }
246 
247  // Force covariance matrix to be symmetric, as the Hessian (supposedly) is
248  *tmpCovMat = .5*((*tmpCovMat) + tmpCovMat->transpose());
249 
250  // Test if covariance matrix is positive definite
251  M lowerChol(*tmpCovMat);
252  if ((m_env.subDisplayFile() ) &&
253  (m_env.displayVerbosity() >= 10)) {
254  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
255  << ", position = " << position
256  << ", stageId = " << stageId
257  << ": calling lowerChol.chol()"
258  << ", lowerChol = " << lowerChol
259  << std::endl;
260  }
261  int iRC = lowerChol.chol();
262  if (iRC) {
263  std::cerr << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition(): chol failed\n";
264  }
265  if ((m_env.subDisplayFile() ) &&
266  (m_env.displayVerbosity() >= 10)) {
267  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
268  << ", position = " << position
269  << ", stageId = " << stageId
270  << ": got lowerChol.chol() with iRC = " << iRC
271  << std::endl;
272  }
273 
274  bool covIsPositiveDefinite = !iRC;
275 
276  if (covIsPositiveDefinite) {
277 
278  // m_env.worldRank(),
279  // "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
280  // "stageId is too large for m_scales");
281  //double factor = 1./m_scales[stageId]/m_scales[stageId];
282  //*tmpCovMat *= factor;
283 
284  m_originalNewtonSteps[stageId] = new V(-1.*(*tmpCovMat)*(*tmpGrad));
285  m_originalCovMatrices[stageId] = new M(*tmpCovMat);
286 
287  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
288  *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
289  << ", position = " << position
290  << ", stageId = " << stageId
291  << ", about to instantiate a Gaussian RV"
292  << ": tmpHessian = " << *tmpHessian
293  << ", preComputingPos = " << *m_preComputingPositions[stageId]
294  << ", tmpCovMat = " << *tmpCovMat
295  << ", tmpGrad = " << *tmpGrad
296  << ", preComputedPos = " << *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]
297  //<< ", factor = " << factor
298  //<< ", rvCov = " << factor*(*tmpCovMat)
299  << std::endl;
300  }
301  m_rvs[stageId] = new GaussianVectorRV<V,M>(m_prefix.c_str(),
302  *m_vectorSpace,
303  *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId],
304  *m_originalCovMatrices[stageId]);
305  }
306  else {
307  validPreComputingPosition = false;
308  }
309 
310  delete tmpGrad;
311  delete tmpCovMat;
312  delete tmpHessian;
313  }
314  else {
315  validPreComputingPosition = false;
316  }
317 
318  if (validPreComputingPosition == false) {
319  // Put "default" values on variables
320  V tmpGrad (m_vectorSpace->zeroVector());
321  M tmpCovMat(tmpGrad,1.); // = identity matrix
322  m_originalNewtonSteps[stageId] = new V(-1.*tmpCovMat*tmpGrad);
323  m_originalCovMatrices[stageId] = new M(tmpCovMat);
324  m_rvs[stageId] = new GaussianVectorRV<V,M>(m_prefix.c_str(),
325  *m_vectorSpace,
326  *m_preComputingPositions[stageId],
327  tmpCovMat);
328  }
329 
330  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
331  *m_env.subDisplayFile() << "Leaving HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
332  << ": position = " << position
333  << ", stageId = " << stageId
334  << std::endl;
335  }
336 
337  return validPreComputingPosition;
338 }
339 //---------------------------------------------------
340 template<class V, class M>
341 void
343 {
344  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalNewtonSteps.size(), "m_preComputingPositions.size() != m_originalNewtonSteps.size()");
345 
346  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalCovMatrices.size(), "m_preComputingPositions.size() != m_originalCovMatrices.size()");
347 
349 
350  // RVs are not deleted in base class because the cov matrices are constant in the case of scaledTK class
351  for (unsigned int i = 0; i < m_rvs.size(); ++i) {
352  if (m_rvs[i]) {
353  delete m_rvs[i];
354  m_rvs[i] = NULL;
355  }
356  }
357 
358  for (unsigned int i = 0; i < m_originalNewtonSteps.size(); ++i) {
359  if (m_originalNewtonSteps[i]) {
360  delete m_originalNewtonSteps[i];
361  m_originalNewtonSteps[i] = NULL;
362  }
363  }
364 
365  for (unsigned int i = 0; i < m_originalCovMatrices.size(); ++i) {
366  if (m_originalCovMatrices[i]) {
367  delete m_originalCovMatrices[i];
368  m_originalCovMatrices[i] = NULL;
369  }
370  }
371 
372  return;
373 }
374 
375 template <class V, class M>
376 unsigned int
378 {
379  unsigned int old_stageId = this->m_stageId;
380  this->m_stageId = stageId;
381  return old_stageId;
382 }
383 
384 
385 // I/O methods---------------------------------------
386 template<class V, class M>
387 void
389 {
391  return;
392 }
393 
394 } // End namespace QUESO
395 
A templated class for synchronizing the calls of scalar functions (BaseScalarFunction and derived cla...
std::vector< const V * > m_preComputingPositions
Definition: TKGroup.h:132
A class representing a vector space.
Definition: VectorSet.h:49
A class representing a Gaussian vector RV.
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
This class allows the representation of a transition kernel with Hessians.
virtual void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
Definition: TKGroup.C:111
HessianCovMatricesTKGroup(const char *prefix, const VectorSpace< V, M > &vectorSpace, const std::vector< double > &scales, const ScalarFunctionSynchronizer< V, M > &targetPdfSynchronizer)
Default constructor.
void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
void updateLawExpVector(const V &newLawExpVector)
Updates the vector that contains the mean values.
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...
std::vector< double > m_scales
Definition: TKGroup.h:131
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"))
void print(std::ostream &os) const
TODO: Prints the transition kernel.
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the covariance matrix.
const GaussianVectorRV< V, M > & rv(unsigned int stageId) const
Gaussian increment property to construct a transition kernel.
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
virtual unsigned int set_dr_stage(unsigned int stageId)
Does nothing. Subclasses may re-implement. Returns the current stage id.
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:55 for queso-0.57.1 by  doxygen 1.8.5