QUESO::HessianCovMatricesTKGroup< V, M > Class Template Reference

This class allows the representation of a transition kernel with Hessians. More...

#include <HessianCovMatricesTKGroup.h>

Inheritance diagram for QUESO::HessianCovMatricesTKGroup< V, M >:

Collaboration diagram for QUESO::HessianCovMatricesTKGroup< V, M >:

Public Member Functions
Constructor/Destructor methods
	HessianCovMatricesTKGroup (const char *prefix, const VectorSpace< V, M > &vectorSpace, const std::vector< double > &scales, const ScalarFunctionSynchronizer< V, M > &targetPdfSynchronizer)
	Default constructor. More...
	~HessianCovMatricesTKGroup ()
	Destructor. More...
Statistical/Mathematical methods
bool	symmetric () const
	Whether or not the matrix is symmetric. Always 'false'. More...
const GaussianVectorRV< V, M > &	rv (unsigned int stageId) const
	Gaussian increment property to construct a transition kernel. More...
const GaussianVectorRV< V, M > &	rv (const std::vector< unsigned int > &stageIds)
	Gaussian increment property to construct a transition kernel. More...
virtual const GaussianVectorRV < V, M > &	rv (const V &position) const
	Constructs transition kernel pdf based on internal m_stageId variable. More...
Misc methods
bool	setPreComputingPosition (const V &position, unsigned int stageId)
	Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position. More...
void	clearPreComputingPositions ()
	Clears the pre-computing positions m_preComputingPositions[stageId]. More...
virtual unsigned int	set_dr_stage (unsigned int stageId)
	Does nothing. Subclasses may re-implement. Returns the current stage id. More...
I/O methods
void	print (std::ostream &os) const
	TODO: Prints the transition kernel. More...
Public Member Functions inherited from QUESO::BaseTKGroup< V, M >
	BaseTKGroup ()
	Default constructor. More...
	BaseTKGroup (const char *prefix, const VectorSpace< V, M > &vectorSpace, const std::vector< double > &scales)
	Constructor. More...
virtual	~BaseTKGroup ()
	Destructor. More...
const BaseEnvironment &	env () const
	QUESO's environment. More...
const V &	preComputingPosition (unsigned int stageId) const
	Pre-computing position; access to protected attribute *m_preComputingPositions[stageId]. More...

Private Attributes
const ScalarFunctionSynchronizer< V, M > &	m_targetPdfSynchronizer
std::vector< V * >	m_originalNewtonSteps
std::vector< M * >	m_originalCovMatrices

Additional Inherited Members
Protected Attributes inherited from QUESO::BaseTKGroup< V, M >
const EmptyEnvironment *	m_emptyEnv
const BaseEnvironment &	m_env
std::string	m_prefix
const VectorSpace< V, M > *	m_vectorSpace
std::vector< double >	m_scales
std::vector< const V * >	m_preComputingPositions
std::vector< BaseVectorRV< V, M > * >	m_rvs
unsigned int	m_stageId

Detailed Description

template<class V = GslVector, class M = GslMatrix>
class QUESO::HessianCovMatricesTKGroup< V, M >

This class allows the representation of a transition kernel with Hessians.

Definition at line 44 of file HessianCovMatricesTKGroup.h.

Constructor & Destructor Documentation

template<class V , class M >

QUESO::HessianCovMatricesTKGroup< V, M >::HessianCovMatricesTKGroup	(	const char *	prefix,
		const VectorSpace< V, M > &	vectorSpace,
		const std::vector< double > &	scales,
		const ScalarFunctionSynchronizer< V, M > &	targetPdfSynchronizer
	)

Default constructor.

Definition at line 33 of file HessianCovMatricesTKGroup.C.

References QUESO::BaseEnvironment::displayVerbosity(), QUESO::BaseTKGroup< V, M >::m_env, QUESO::HessianCovMatricesTKGroup< V, M >::m_originalCovMatrices, QUESO::HessianCovMatricesTKGroup< V, M >::m_originalNewtonSteps, QUESO::BaseTKGroup< V, M >::m_preComputingPositions, QUESO::BaseTKGroup< V, M >::m_rvs, QUESO::BaseTKGroup< V, M >::m_scales, and QUESO::BaseEnvironment::subDisplayFile().

  :
  BaseTKGroup<V,M>(prefix,vectorSpace,scales),
  m_targetPdfSynchronizer(targetPdfSynchronizer),
  m_originalNewtonSteps  (scales.size()+1,NULL), // Yes, +1
  m_originalCovMatrices  (scales.size()+1,NULL)  // Yes, +1
{
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "Entering HessianCovMatricesTKGroup<V,M>::constructor()"
                           << std::endl;
  }

  m_rvs.resize(scales.size()+1,NULL); // Yes, +1 (IMPORTANT: overwrite initialization done by BaseTKGroup<V,M>)

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::constructor()"
                           << ": m_scales.size() = "                   << m_scales.size()
                           << ", m_preComputingPositions.size() = "    << m_preComputingPositions.size()
                           << ", m_rvs.size() = "                      << m_rvs.size()
                           << ", m_originalNewtonSteps.size() = "      << m_originalNewtonSteps.size()
                           << ", m_originalCovMatrices.size() = "      << m_originalCovMatrices.size()
                           << std::endl;
  }

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "Leaving HessianCovMatricesTKGroup<V,M>::constructor()"
                           << std::endl;
  }
}

template<class V , class M >

QUESO::HessianCovMatricesTKGroup< V, M >::~HessianCovMatricesTKGroup

(

)

Destructor.

Definition at line 68 of file HessianCovMatricesTKGroup.C.

{
}

Member Function Documentation

template<class V , class M >

void QUESO::HessianCovMatricesTKGroup< V, M >::clearPreComputingPositions ( )

virtual

Clears the pre-computing positions m_preComputingPositions[stageId].

Reimplemented from QUESO::BaseTKGroup< V, M >.

Definition at line 342 of file HessianCovMatricesTKGroup.C.

References QUESO::BaseTKGroup< V, M >::clearPreComputingPositions(), and queso_require_equal_to_msg.

{
  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalNewtonSteps.size(), "m_preComputingPositions.size() != m_originalNewtonSteps.size()");

  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalCovMatrices.size(), "m_preComputingPositions.size() != m_originalCovMatrices.size()");

  BaseTKGroup<V,M>::clearPreComputingPositions();

  // RVs are not deleted in base class because the cov matrices are constant in the case of scaledTK class
  for (unsigned int i = 0; i < m_rvs.size(); ++i) {
    if (m_rvs[i]) {
      delete m_rvs[i];
      m_rvs[i] = NULL;
    }
  }

  for (unsigned int i = 0; i < m_originalNewtonSteps.size(); ++i) {
    if (m_originalNewtonSteps[i]) {
      delete m_originalNewtonSteps[i];
      m_originalNewtonSteps[i] = NULL;
    }
  }

  for (unsigned int i = 0; i < m_originalCovMatrices.size(); ++i) {
    if (m_originalCovMatrices[i]) {
      delete m_originalCovMatrices[i];
      m_originalCovMatrices[i] = NULL;
    }
  }

  return;
}

template<class V , class M >

void QUESO::HessianCovMatricesTKGroup< V, M >::print ( std::ostream &  os ) const

virtual

TODO: Prints the transition kernel.

Todo:

: implement me!

Reimplemented from QUESO::BaseTKGroup< V, M >.

Definition at line 388 of file HessianCovMatricesTKGroup.C.

References QUESO::BaseTKGroup< V, M >::print().

{
  BaseTKGroup<V,M>::print(os);
  return;
}

template<class V , class M >

const GaussianVectorRV< V, M > & QUESO::HessianCovMatricesTKGroup< V, M >::rv ( unsigned int  stageId ) const

virtual

Gaussian increment property to construct a transition kernel.

Implements QUESO::BaseTKGroup< V, M >.

Definition at line 81 of file HessianCovMatricesTKGroup.C.

References queso_require_greater_msg, queso_require_msg, QUESO::GaussianVectorRV< V, M >::updateLawCovMatrix(), and QUESO::GaussianVectorRV< V, M >::updateLawExpVector().

{
  queso_require_greater_msg(m_rvs.size(), stageId, "m_rvs.size() <= stageId");

  queso_require_msg(m_rvs[stageId], "m_rvs[stageId] == NULL");

  queso_require_greater_msg(m_preComputingPositions.size(), stageId, "m_preComputingPositions.size() <= stageId");

  queso_require_msg(m_preComputingPositions[stageId], "m_preComputingPositions[stageId] == NULL");

  GaussianVectorRV<V, M> * gaussian_rv =
    dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[stageId]);

  gaussian_rv->updateLawExpVector(*m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]);

  if ((m_env.subDisplayFile()        ) &&
      (m_env.displayVerbosity() >= 10)) {
    *m_env.subDisplayFile() << "In HessianCovMatrixTKGroup<V,M>::rv1()"
                            << ", stageId = " << stageId
                            << ": about to call m_rvs[stageId]->updateLawCovMatrix()"
                            << ", covMatrix = \n" << *m_originalCovMatrices[stageId] // FIX ME: might demand parallelism
                            << std::endl;
  }

  gaussian_rv->updateLawCovMatrix(*m_originalCovMatrices[stageId]);

  return *gaussian_rv;
}

template<class V , class M >

const GaussianVectorRV< V, M > & QUESO::HessianCovMatricesTKGroup< V, M >::rv ( const std::vector< unsigned int > &  stageIds )

virtual

Gaussian increment property to construct a transition kernel.

Implements QUESO::BaseTKGroup< V, M >.

Definition at line 112 of file HessianCovMatricesTKGroup.C.

References queso_require_greater_msg, queso_require_msg, QUESO::GaussianVectorRV< V, M >::updateLawCovMatrix(), and QUESO::GaussianVectorRV< V, M >::updateLawExpVector().

{
  queso_require_greater_msg(m_rvs.size(), stageIds[0], "m_rvs.size() <= stageIds[0]");

  queso_require_msg(m_rvs[stageIds[0]], "m_rvs[stageIds[0]] == NULL");

  queso_require_greater_msg(m_preComputingPositions.size(), stageIds[0], "m_preComputingPositions.size() <= stageIds[0]");

  queso_require_msg(m_preComputingPositions[stageIds[0]], "m_preComputingPositions[stageIds[0]] == NULL");

  double factor = 1./m_scales[stageIds.size()-1]/m_scales[stageIds.size()-1];

  GaussianVectorRV<V, M> * gaussian_rv =
    dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[stageIds[0]]);

  gaussian_rv->updateLawExpVector(*m_preComputingPositions[stageIds[0]] + factor*(*m_originalNewtonSteps[stageIds[0]]));

  if ((m_env.subDisplayFile()        ) &&
      (m_env.displayVerbosity() >= 10)) {
    *m_env.subDisplayFile() << "In HessianCovMatrixTKGroup<V,M>::rv2()"
                            << ", stageIds.size() = " << stageIds.size()
                            << ", stageIds[0] = "     << stageIds[0]
                            << ", factor = "          << factor
                            << ": about to call m_rvs[stageIds[0]]->updateLawCovVector()"
                            << ", covMatrix = \n" << factor*(*m_originalCovMatrices[stageIds[0]]) // FIX ME: might demand parallelism
                            << std::endl;
  }
  gaussian_rv->updateLawCovMatrix(factor*(*m_originalCovMatrices[stageIds[0]]));

  return *gaussian_rv;
}

template<class V , class M >

const GaussianVectorRV< V, M > & QUESO::HessianCovMatricesTKGroup< V, M >::rv ( const V &  position ) const

virtual

Constructs transition kernel pdf based on internal m_stageId variable.

Implements QUESO::BaseTKGroup< V, M >.

Definition at line 146 of file HessianCovMatricesTKGroup.C.

References queso_require_greater_msg, queso_require_msg, QUESO::GaussianVectorRV< V, M >::updateLawCovMatrix(), and QUESO::GaussianVectorRV< V, M >::updateLawExpVector().

{
  queso_require_greater_msg(m_rvs.size(), this->m_stageId, "m_rvs.size() <= stageId");
  queso_require_msg(m_rvs[this->m_stageId], "m_rvs[stageId] == NULL");
  // queso_require_greater_msg(m_preComputingPositions.size(), this->m_stageId, "m_preComputingPositions.size() <= stageId");
  // queso_require_msg(m_preComputingPositions[this->m_stageId], "m_preComputingPositions[stageId] == NULL");

  GaussianVectorRV<V, M> * gaussian_rv =
    dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[this->m_stageId]);

  gaussian_rv->updateLawExpVector(position + *m_originalNewtonSteps[this->m_stageId]);

  gaussian_rv->updateLawCovMatrix(*m_originalCovMatrices[this->m_stageId]);

  return *gaussian_rv;
}

template<class V , class M >

unsigned int QUESO::HessianCovMatricesTKGroup< V, M >::set_dr_stage ( unsigned int  stageId )

virtual

Does nothing. Subclasses may re-implement. Returns the current stage id.

Reimplemented from QUESO::BaseTKGroup< V, M >.

Definition at line 377 of file HessianCovMatricesTKGroup.C.

{
  unsigned int old_stageId = this->m_stageId;
  this->m_stageId = stageId;
  return old_stageId;
}

template<class V , class M >

bool QUESO::HessianCovMatricesTKGroup< V, M >::setPreComputingPosition ( const V &  position, unsigned int  stageId  )

virtual

Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position.

Reimplemented from QUESO::BaseTKGroup< V, M >.

Definition at line 166 of file HessianCovMatricesTKGroup.C.

References queso_require_equal_to_msg, queso_require_greater_msg, queso_require_msg, and QUESO::BaseTKGroup< V, M >::setPreComputingPosition().

{
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "Entering HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                           << ": position = " << position
                           << ", stageId = "  << stageId
                           << std::endl;
  }

  bool validPreComputingPosition = true;

  // Verify consistency of sizes
  queso_require_greater_msg(m_preComputingPositions.size(), stageId, "m_preComputingPositions.size() <= stageId");

  queso_require_equal_to_msg(m_preComputingPositions.size(), m_rvs.size(), "m_preComputingPositions.size() != m_rvs.size()");

  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalNewtonSteps.size(), "m_preComputingPositions.size() != m_originalNewtonSteps.size()");

  queso_require_equal_to_msg(m_preComputingPositions.size(), m_originalCovMatrices.size(), "m_preComputingPositions.size() != m_originalCovMatrices.size()");

  // Verify data is not null
  queso_require_msg(!(m_preComputingPositions[stageId]), "m_preComputingPositions[stageId] != NULL");

  queso_require_msg(!(m_rvs[stageId]), "m_rvs[stageId] != NULL");

  queso_require_msg(!(m_originalNewtonSteps[stageId]), "m_originalNewtonSteps[stageId] != NULL");

  queso_require_msg(!(m_originalCovMatrices[stageId]), "m_originalCovMatrices[stageId] != NULL");

  BaseTKGroup<V,M>::setPreComputingPosition(position,stageId);

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                           << ", position = "                          << position
                           << ", stageId = "                           << stageId
                           << ": m_originalNewtonSteps.size() = "      << m_originalNewtonSteps.size()
                           << ", m_originalCovMatrices.size() = "      << m_originalCovMatrices.size()
                           << ", m_preComputingPositions.size() = "    << m_preComputingPositions.size()
                           << ", m_rvs.size() = "                      << m_rvs.size()
                           << std::endl;
  }

  if (m_targetPdfSynchronizer.domainSet().contains(position)) {
    M* tmpHessian = m_vectorSpace->newMatrix();
    M* tmpCovMat  = m_vectorSpace->newMatrix();
    V* tmpGrad    = m_vectorSpace->newVector();

    double logPrior = 0.;
    double logLikelihood = 0.;
    double logTarget = 0.;
    logTarget = m_targetPdfSynchronizer.callFunction(&position, // Might demand parallel environment
                                                     NULL,
                                                     tmpGrad,
                                                     tmpHessian,
                                                     NULL,
                                                     &logPrior,
                                                     &logLikelihood);
    if (logTarget) {}; // just to remove compiler warning

    // IMPORTANT: covariance matrix = (Hessian)^{-1} !!!
    V unitVector(m_vectorSpace->zeroVector());
    V multVector(m_vectorSpace->zeroVector());
    for (unsigned int j = 0; j < tmpHessian->numCols(); ++j) {
      if (j > 0) unitVector[j-1] = 0.;
      unitVector[j] = 1.;
      tmpHessian->invertMultiply(unitVector, multVector);
      for (unsigned int i = 0; i < tmpHessian->numRowsLocal(); ++i) {
        (*tmpCovMat)(i,j) = multVector[i];
      }
    }
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
      *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                             << ", position = "  << position
                             << ", stageId = "   << stageId
                             << ":\n H = "       << *tmpHessian
                             << "\n H^{-1} = "   << *tmpCovMat
                             << "\n H*H^{-1} = " << (*tmpHessian)*(*tmpCovMat)
                             << "\n H^{-1}*H = " << (*tmpCovMat)*(*tmpHessian)
                             << std::endl;
    }

    // Force covariance matrix to be symmetric, as the Hessian (supposedly) is
    *tmpCovMat = .5*((*tmpCovMat) + tmpCovMat->transpose());

    // Test if covariance matrix is positive definite
    M lowerChol(*tmpCovMat);
    if ((m_env.subDisplayFile()        ) &&
        (m_env.displayVerbosity() >= 10)) {
      *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                              << ", position = "  << position
                              << ", stageId = "   << stageId
                              << ": calling lowerChol.chol()"
                              << ", lowerChol = " << lowerChol
                              << std::endl;
    }
    int iRC = lowerChol.chol();
    if (iRC) {
      std::cerr << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition(): chol failed\n";
    }
    if ((m_env.subDisplayFile()        ) &&
        (m_env.displayVerbosity() >= 10)) {
      *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                              << ", position = "  << position
                              << ", stageId = "   << stageId
                              << ": got lowerChol.chol() with iRC = " << iRC
                              << std::endl;
    }

    bool covIsPositiveDefinite = !iRC;

    if (covIsPositiveDefinite) {

      //                    m_env.worldRank(),
      //                    "HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()",
      //                    "stageId is too large for m_scales");
      //double factor = 1./m_scales[stageId]/m_scales[stageId];
      //*tmpCovMat *= factor;

      m_originalNewtonSteps[stageId] = new V(-1.*(*tmpCovMat)*(*tmpGrad));
      m_originalCovMatrices[stageId] = new M(*tmpCovMat);

      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
        *m_env.subDisplayFile() << "In HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                               << ", position = "        << position
                               << ", stageId = "         << stageId
                               << ", about to instantiate a Gaussian RV"
                               << ": tmpHessian = "      << *tmpHessian
                               << ", preComputingPos = " << *m_preComputingPositions[stageId]
                               << ", tmpCovMat = "       << *tmpCovMat
                               << ", tmpGrad = "         << *tmpGrad
                               << ", preComputedPos = "  << *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId]
                             //<< ", factor = "          << factor
                             //<< ", rvCov = "           << factor*(*tmpCovMat)
                               << std::endl;
      }
      m_rvs[stageId] = new GaussianVectorRV<V,M>(m_prefix.c_str(),
                                                        *m_vectorSpace,
                                                        *m_preComputingPositions[stageId] + *m_originalNewtonSteps[stageId],
                                                        *m_originalCovMatrices[stageId]);
    }
    else {
      validPreComputingPosition = false;
    }

    delete tmpGrad;
    delete tmpCovMat;
    delete tmpHessian;
  }
  else {
    validPreComputingPosition = false;
  }

  if (validPreComputingPosition == false) {
    // Put "default" values on variables
    V tmpGrad  (m_vectorSpace->zeroVector());
    M tmpCovMat(tmpGrad,1.); // = identity matrix
    m_originalNewtonSteps[stageId] = new V(-1.*tmpCovMat*tmpGrad);
    m_originalCovMatrices[stageId] = new M(tmpCovMat);
    m_rvs[stageId] = new GaussianVectorRV<V,M>(m_prefix.c_str(),
                                                      *m_vectorSpace,
                                                      *m_preComputingPositions[stageId],
                                                      tmpCovMat);
  }

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "Leaving HessianCovMatricesTKGroup<V,M>::setPreComputingPosition()"
                           << ": position = " << position
                           << ", stageId = "  << stageId
                           << std::endl;
  }

  return validPreComputingPosition;
}

template<class V , class M >

bool QUESO::HessianCovMatricesTKGroup< V, M >::symmetric ( ) const

virtual

Whether or not the matrix is symmetric. Always 'false'.

Todo:

: It only returns 'false', thus a test for its symmetricity must be included.

Implements QUESO::BaseTKGroup< V, M >.

Definition at line 74 of file HessianCovMatricesTKGroup.C.

{
  return false;
}

Member Data Documentation

template<class V = GslVector, class M = GslMatrix>

std::vector<M*> QUESO::HessianCovMatricesTKGroup< V, M >::m_originalCovMatrices

private

Definition at line 100 of file HessianCovMatricesTKGroup.h.

Referenced by QUESO::HessianCovMatricesTKGroup< V, M >::HessianCovMatricesTKGroup().

template<class V = GslVector, class M = GslMatrix>

std::vector<V*> QUESO::HessianCovMatricesTKGroup< V, M >::m_originalNewtonSteps

private

Definition at line 99 of file HessianCovMatricesTKGroup.h.

Referenced by QUESO::HessianCovMatricesTKGroup< V, M >::HessianCovMatricesTKGroup().

template<class V = GslVector, class M = GslMatrix>

const ScalarFunctionSynchronizer<V,M>& QUESO::HessianCovMatricesTKGroup< V, M >::m_targetPdfSynchronizer

private

Definition at line 98 of file HessianCovMatricesTKGroup.h.

---

The documentation for this class was generated from the following files:

• src/stats/inc/HessianCovMatricesTKGroup.h
• src/stats/src/HessianCovMatricesTKGroup.C