25 #include <queso/HessianCovMatricesTKGroup.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
32 template<
class V,
class M>
36 const std::vector<double>& scales,
40 m_targetPdfSynchronizer(targetPdfSynchronizer),
41 m_originalNewtonSteps (scales.size()+1,NULL),
42 m_originalCovMatrices (scales.size()+1,NULL)
49 m_rvs.resize(scales.size()+1,NULL);
53 <<
": m_scales.size() = " <<
m_scales.size()
55 <<
", m_rvs.size() = " <<
m_rvs.size()
67 template<
class V,
class M>
72 template<
class V,
class M>
79 template<
class V,
class M>
85 "HessianCovMatricesTKGroup<V,M>::rv1()",
86 "m_rvs.size() <= stageId");
90 "HessianCovMatricesTKGroup<V,M>::rv1()",
91 "m_rvs[stageId] == NULL");
95 "HessianCovMatricesTKGroup<V,M>::rv1()",
96 "m_preComputingPositions.size() <= stageId");
100 "HessianCovMatricesTKGroup<V,M>::rv1()",
101 "m_preComputingPositions[stageId] == NULL");
106 gaussian_rv->
updateLawExpVector(*m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]);
108 if ((m_env.subDisplayFile() ) &&
109 (m_env.displayVerbosity() >= 10)) {
110 *m_env.subDisplayFile() <<
"In HessianCovMatrixTKGroup<V,M>::rv1()"
111 <<
", stageId = " << stageId
112 <<
": about to call m_rvs[stageId]->updateLawCovMatrix()"
113 <<
", covMatrix = \n" << *m_originalCovMatrices[stageId]
117 gaussian_rv->updateLawCovMatrix(*m_originalCovMatrices[stageId]);
122 template<
class V,
class M>
128 "HessianCovMatricesTKGroup<V,M>::rv2()",
129 "m_rvs.size() <= stageIds[0]");
133 "HessianCovMatricesTKGroup<V,M>::rv2()",
134 "m_rvs[stageIds[0]] == NULL");
138 "HessianCovMatricesTKGroup<V,M>::rv2()",
139 "m_preComputingPositions.size() <= stageIds[0]");
143 "HessianCovMatricesTKGroup<V,M>::rv2()",
144 "m_preComputingPositions[stageIds[0]] == NULL");
146 double factor = 1./m_scales[stageIds.size()-1]/m_scales[stageIds.size()-1];
151 gaussian_rv->
updateLawExpVector(*m_preComputingPositions[stageIds[0]] + factor*(*m_originalNewtonSteps[stageIds[0]]));
153 if ((m_env.subDisplayFile() ) &&
154 (m_env.displayVerbosity() >= 10)) {
155 *m_env.subDisplayFile() <<
"In HessianCovMatrixTKGroup<V,M>::rv2()"
156 <<
", stageIds.size() = " << stageIds.size()
157 <<
", stageIds[0] = " << stageIds[0]
158 <<
", factor = " << factor
159 <<
": about to call m_rvs[stageIds[0]]->updateLawCovVector()"
160 <<
", covMatrix = \n" << factor*(*m_originalCovMatrices[stageIds[0]])
168 template<
class V,
class M>
172 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
173 *m_env.subDisplayFile() <<
"Entering HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
174 <<
": position = " << position
175 <<
", stageId = " << stageId
179 bool validPreComputingPosition =
true;
184 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
185 "m_preComputingPositions.size() <= stageId");
189 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
190 "m_preComputingPositions.size() != m_rvs.size()");
194 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
195 "m_preComputingPositions.size() != m_originalNewtonSteps.size()");
199 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
200 "m_preComputingPositions.size() != m_originalCovMatrices.size()");
205 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
206 "m_preComputingPositions[stageId] != NULL");
210 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
211 "m_rvs[stageId] != NULL");
215 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
216 "m_originalNewtonSteps[stageId] != NULL");
220 "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
221 "m_originalCovMatrices[stageId] != NULL");
225 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
226 *m_env.subDisplayFile() <<
"In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
227 <<
", position = " << position
228 <<
", stageId = " << stageId
229 <<
": m_originalNewtonSteps.size() = " << m_originalNewtonSteps.size()
230 <<
", m_originalCovMatrices.size() = " << m_originalCovMatrices.size()
231 <<
", m_preComputingPositions.size() = " << m_preComputingPositions.size()
232 <<
", m_rvs.size() = " << m_rvs.size()
236 if (m_targetPdfSynchronizer.domainSet().contains(position)) {
237 M* tmpHessian = m_vectorSpace->newMatrix();
238 M* tmpCovMat = m_vectorSpace->newMatrix();
239 V* tmpGrad = m_vectorSpace->newVector();
241 double logPrior = 0.;
242 double logLikelihood = 0.;
243 double logTarget = 0.;
244 logTarget = m_targetPdfSynchronizer.callFunction(&position,
254 V unitVector(m_vectorSpace->zeroVector());
255 V multVector(m_vectorSpace->zeroVector());
256 for (
unsigned int j = 0; j < tmpHessian->numCols(); ++j) {
257 if (j > 0) unitVector[j-1] = 0.;
259 tmpHessian->invertMultiply(unitVector, multVector);
260 for (
unsigned int i = 0; i < tmpHessian->numRowsLocal(); ++i) {
261 (*tmpCovMat)(i,j) = multVector[i];
264 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
265 *m_env.subDisplayFile() <<
"In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
266 <<
", position = " << position
267 <<
", stageId = " << stageId
268 <<
":\n H = " << *tmpHessian
269 <<
"\n H^{-1} = " << *tmpCovMat
270 <<
"\n H*H^{-1} = " << (*tmpHessian)*(*tmpCovMat)
271 <<
"\n H^{-1}*H = " << (*tmpCovMat)*(*tmpHessian)
276 *tmpCovMat = .5*((*tmpCovMat) + tmpCovMat->transpose());
279 M lowerChol(*tmpCovMat);
280 if ((m_env.subDisplayFile() ) &&
281 (m_env.displayVerbosity() >= 10)) {
282 *m_env.subDisplayFile() <<
"In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
283 <<
", position = " << position
284 <<
", stageId = " << stageId
285 <<
": calling lowerChol.chol()"
286 <<
", lowerChol = " << lowerChol
289 int iRC = lowerChol.chol();
291 std::cerr <<
"In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition(): chol failed\n";
293 if ((m_env.subDisplayFile() ) &&
294 (m_env.displayVerbosity() >= 10)) {
295 *m_env.subDisplayFile() <<
"In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
296 <<
", position = " << position
297 <<
", stageId = " << stageId
298 <<
": got lowerChol.chol() with iRC = " << iRC
302 bool covIsPositiveDefinite = !iRC;
304 if (covIsPositiveDefinite) {
312 m_originalNewtonSteps[stageId] =
new V(-1.*(*tmpCovMat)*(*tmpGrad));
313 m_originalCovMatrices[stageId] =
new M(*tmpCovMat);
315 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
316 *m_env.subDisplayFile() <<
"In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
317 <<
", position = " << position
318 <<
", stageId = " << stageId
319 <<
", about to instantiate a Gaussian RV"
320 <<
": tmpHessian = " << *tmpHessian
321 <<
", preComputingPos = " << *m_preComputingPositions[stageId]
322 <<
", tmpCovMat = " << *tmpCovMat
323 <<
", tmpGrad = " << *tmpGrad
324 <<
", preComputedPos = " << *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]
331 *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId],
332 *m_originalCovMatrices[stageId]);
335 validPreComputingPosition =
false;
343 validPreComputingPosition =
false;
346 if (validPreComputingPosition ==
false) {
348 V tmpGrad (m_vectorSpace->zeroVector());
349 M tmpCovMat(tmpGrad,1.);
350 m_originalNewtonSteps[stageId] =
new V(-1.*tmpCovMat*tmpGrad);
351 m_originalCovMatrices[stageId] =
new M(tmpCovMat);
354 *m_preComputingPositions[stageId],
358 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
359 *m_env.subDisplayFile() <<
"Leaving HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
360 <<
": position = " << position
361 <<
", stageId = " << stageId
365 return validPreComputingPosition;
368 template<
class V,
class M>
374 "HessianCovMatricesTKGroup<V,M>::clearPreComputingPositions()",
375 "m_preComputingPositions.size() != m_originalNewtonSteps.size()");
379 "HessianCovMatricesTKGroup<V,M>::clearPreComputingPositions()",
380 "m_preComputingPositions.size() != m_originalCovMatrices.size()");
385 for (
unsigned int i = 0; i < m_rvs.size(); ++i) {
392 for (
unsigned int i = 0; i < m_originalNewtonSteps.size(); ++i) {
393 if (m_originalNewtonSteps[i]) {
394 delete m_originalNewtonSteps[i];
395 m_originalNewtonSteps[i] = NULL;
399 for (
unsigned int i = 0; i < m_originalCovMatrices.size(); ++i) {
400 if (m_originalCovMatrices[i]) {
401 delete m_originalCovMatrices[i];
402 m_originalCovMatrices[i] = NULL;
409 template<
class V,
class M>
std::vector< BaseVectorRV< V, M > * > m_rvs
std::vector< double > m_scales
This class allows the representation of a transition kernel with Hessians.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
std::vector< M * > m_originalCovMatrices
virtual void print(std::ostream &os) const
TODO: Prints the transition kernel.
void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
const GaussianVectorRV< V, M > & rv(unsigned int stageId) const
Gaussian increment property to construct a transition kernel.
virtual void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
const BaseEnvironment & m_env
~HessianCovMatricesTKGroup()
Destructor.
This base class allows the representation of a transition kernel.
void updateLawExpVector(const V &newLawExpVector)
Updates the vector that contains the mean values.
void print(std::ostream &os) const
TODO: Prints the transition kernel.
HessianCovMatricesTKGroup(const char *prefix, const VectorSpace< V, M > &vectorSpace, const std::vector< double > &scales, const ScalarFunctionSynchronizer< V, M > &targetPdfSynchronizer)
Default constructor.
A templated class for synchronizing the calls of scalar functions (BaseScalarFunction and derived cla...
A class representing a vector space.
bool setPreComputingPosition(const V &position, unsigned int stageId)
Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position...
unsigned int displayVerbosity() const
std::vector< V * > m_originalNewtonSteps
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
virtual bool setPreComputingPosition(const V &position, unsigned int stageId)
Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position...
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the covariance matrix.
bool symmetric() const
Whether or not the matrix is symmetric. Always 'false'.
std::vector< const V * > m_preComputingPositions
A class representing a Gaussian vector RV.