This class allows the representation of the MALA transition kernel with a scaled covariance matrix for the purposes of delayed rejection. More...

#include <MetropolisAdjustedLangevinTK.h>

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

[legend]

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

[legend]

Public Member Functions
Constructor/Destructor methods
	MetropolisAdjustedLangevinTK (const char *prefix, const BayesianJointPdf< V, M > &targetPdf, const std::vector< double > &scales, const M &covMatrix)
	Default constructor. More...

	~MetropolisAdjustedLangevinTK ()
	Destructor. More...

Statistical/Mathematical methods
bool	symmetric () const
	Whether or not the matrix is symmetric. Always 'true'. 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...

void	updateLawCovMatrix (const M &covMatrix)
	Scales the covariance matrix. 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...

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...

virtual unsigned int	set_dr_stage (unsigned int stageId)
	Does nothing. Subclasses may re-implement. Returns the current stage id. More...

Private Member Functions
void	setRVsWithZeroMean ()
	Sets the mean of the RVs to zero. More...

Private Attributes
M	m_originalCovMatrix

const BayesianJointPdf< V, M > &	m_targetPdf

double	m_time_step

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::MetropolisAdjustedLangevinTK< V, M >

This class allows the representation of the MALA transition kernel with a scaled covariance matrix for the purposes of delayed rejection.

Definition at line 43 of file MetropolisAdjustedLangevinTK.h.

Constructor & Destructor Documentation

template<class V , class M >

QUESO::MetropolisAdjustedLangevinTK< V, M >::MetropolisAdjustedLangevinTK	(	const char *	prefix,
		const BayesianJointPdf< V, M > &	targetPdf,
		const std::vector< double > &	scales,
		const M &	covMatrix
	)

Default constructor.

Definition at line 34 of file MetropolisAdjustedLangevinTK.C.

References QUESO::BaseEnvironment::displayVerbosity(), QUESO::BaseTKGroup< V, M >::m_env, QUESO::MetropolisAdjustedLangevinTK< V, M >::m_originalCovMatrix, QUESO::BaseTKGroup< V, M >::m_preComputingPositions, QUESO::BaseTKGroup< V, M >::m_rvs, QUESO::BaseTKGroup< V, M >::m_scales, QUESO::MetropolisAdjustedLangevinTK< V, M >::setRVsWithZeroMean(), and QUESO::BaseEnvironment::subDisplayFile().

   :
   BaseTKGroup<V, M>(prefix, targetPdf.domainSet().vectorSpace(), scales),
   m_originalCovMatrix(covMatrix),
   m_targetPdf(targetPdf),
   m_time_step(1.0)
 {
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
     *m_env.subDisplayFile() << "Entering MetropolisAdjustedLangevinTK<V, M>::constructor()"
                            << std::endl;
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
     *m_env.subDisplayFile() << "In MetropolisAdjustedLangevinTK<V, M>::constructor()"
                            << ": m_scales.size() = "                << m_scales.size()
                            << ", m_preComputingPositions.size() = " << m_preComputingPositions.size()
                            << ", m_rvs.size() = "                   << m_rvs.size()
                            << ", m_originalCovMatrix = "            << m_originalCovMatrix
                            << std::endl;
   }
 
   setRVsWithZeroMean();
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
     *m_env.subDisplayFile() << "Leaving MetropolisAdjustedLangevinTK<V, M>::constructor()"
                            << std::endl;
   }
 }

template<class V , class M >

QUESO::MetropolisAdjustedLangevinTK< V, M >::~MetropolisAdjustedLangevinTK ( )

Destructor.

Definition at line 68 of file MetropolisAdjustedLangevinTK.C.

69 {

70 }

Member Function Documentation

template<class V , class M >

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

virtual

Clears the pre-computing positions m_preComputingPositions[stageId].

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

Definition at line 223 of file MetropolisAdjustedLangevinTK.C.

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

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

template<class V , class M >

void QUESO::MetropolisAdjustedLangevinTK< 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 247 of file MetropolisAdjustedLangevinTK.C.

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

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

template<class V , class M >

const GaussianVectorRV< V, M > & QUESO::MetropolisAdjustedLangevinTK< 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 MetropolisAdjustedLangevinTK.C.

References queso_require, queso_require_greater, queso_require_not_equal_to, and QUESO::GaussianVectorRV< V, M >::updateLawExpVector().

 {
   queso_require_not_equal_to(m_rvs.size(), 0);
   queso_require(m_rvs[0]);
   queso_require_greater(m_preComputingPositions.size(), stageId);
   queso_require(m_preComputingPositions[stageId]);
 
   if ((m_env.subDisplayFile()        ) &&
       (m_env.displayVerbosity() >= 10)) {
     *m_env.subDisplayFile() << "In MetropolisAdjustedLangevinTK<V, M>::rv1()"
                             << ", stageId = " << stageId
                             << ": about to call m_rvs[0]->updateLawExpVector()"
                             << ", vector = " << *m_preComputingPositions[stageId] // FIX ME: might demand parallelism
                             << std::endl;
   }
 
   GaussianVectorRV<V, M> * gaussian_rv = dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[0]);
 
   gaussian_rv->updateLawExpVector(*m_preComputingPositions[stageId]);
 
   return (*gaussian_rv);
 }

template<class V , class M >

const GaussianVectorRV< V, M > & QUESO::MetropolisAdjustedLangevinTK< 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 106 of file MetropolisAdjustedLangevinTK.C.

References queso_require, queso_require_greater, queso_require_greater_equal, and QUESO::GaussianVectorRV< V, M >::updateLawExpVector().

 {
   queso_require_greater_equal(m_rvs.size(), stageIds.size());
   queso_require(m_rvs[stageIds.size()-1]);
   queso_require_greater(m_preComputingPositions.size(), stageIds[0]);
   queso_require(m_preComputingPositions[stageIds[0]]);
 
   if ((m_env.subDisplayFile()        ) &&
       (m_env.displayVerbosity() >= 10)) {
     *m_env.subDisplayFile() << "In MetropolisAdjustedLangevinTK<V, M>::rv2()"
                             << ", stageIds.size() = " << stageIds.size()
                             << ", stageIds[0] = "     << stageIds[0]
                             << ": about to call m_rvs[stageIds.size()-1]->updateLawExpVector()"
                             << ", vector = " << *m_preComputingPositions[stageIds[0]] // FIX ME: might demand parallelism
                             << std::endl;
   }
 
   GaussianVectorRV<V, M> * gaussian_rv = dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[stageIds.size()-1]);
 
   gaussian_rv->updateLawExpVector(*m_preComputingPositions[stageIds[0]]);
 
   return (*gaussian_rv);
 }

template<class V , class M >

const GaussianVectorRV< V, M > & QUESO::MetropolisAdjustedLangevinTK< 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 132 of file MetropolisAdjustedLangevinTK.C.

References queso_require, queso_require_not_equal_to, and QUESO::GaussianVectorRV< V, M >::updateLawExpVector().

 {
   queso_require_not_equal_to(m_rvs.size(), 0);
   queso_require(m_rvs[0]);
 
   GaussianVectorRV<V, M> * gaussian_rv = dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[this->m_stageId]);
 
   // 'position' is a position in the chain?  Hopefully?  Anyway, assume it's a
   // position in the chain and then that means we need to modify it slightly so
   // to get the transition distribution for the next state in the chain.
 
   V grad(this->m_targetPdf.domainSet().vectorSpace().zeroVector());
 
   // Get the gradient of the log-posterior.  This is so inefficient it's
   // painful.  We should be caching the gradient evaluations.
   this->m_targetPdf.lnValue(position, NULL, &grad, NULL, NULL);
 
   // Euler time-step
   grad *= 0.5 * this->m_time_step;
 
   // Add on current position
   grad += position;
 
   // Update the mean of the transition kernel.  The vector gets copied, so
   // we're ok.
   gaussian_rv->updateLawExpVector(grad);
 
   return (*gaussian_rv);
 }

template<class V , class M >

bool QUESO::MetropolisAdjustedLangevinTK< 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 185 of file MetropolisAdjustedLangevinTK.C.

References QUESO::GaussianJointPdf< V, M >::lawCovMatrix(), and QUESO::BaseTKGroup< V, M >::setPreComputingPosition().

 {
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
     *m_env.subDisplayFile() << "Entering MetropolisAdjustedLangevinTK<V, M>::setPreComputingPosition()"
                            << ": position = " << position
                            << ", stageId = "  << stageId
                            << std::endl;
   }
 
   BaseTKGroup<V, M>::setPreComputingPosition(position, stageId);
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
     *m_env.subDisplayFile() << "In MetropolisAdjustedLangevinTK<V, M>::setPreComputingPosition()"
                            << ", position = "        << position
                            << ", stageId = "         << stageId
                            << ": preComputingPos = " << *m_preComputingPositions[stageId];
     if (stageId < m_scales.size()) {
       *m_env.subDisplayFile() << ", factor = " << 1./m_scales[stageId]/m_scales[stageId];
     }
     if (stageId < m_rvs.size()) {
       const GaussianJointPdf<V, M>* pdfPtr = dynamic_cast< const GaussianJointPdf<V, M>* >(&(m_rvs[stageId]->pdf()));
       *m_env.subDisplayFile() << ", rvCov = " << pdfPtr->lawCovMatrix(); // FIX ME: might demand parallelism
     }
     *m_env.subDisplayFile() << std::endl;
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
     *m_env.subDisplayFile() << "Leaving MetropolisAdjustedLangevinTK<V, M>::setPreComputingPosition()"
                            << ": position = " << position
                            << ", stageId = "  << stageId
                            << std::endl;
   }
 
   return true;
 }

template<class V , class M >

void QUESO::MetropolisAdjustedLangevinTK< V, M >::setRVsWithZeroMean ( )

private

Sets the mean of the RVs to zero.

Definition at line 230 of file MetropolisAdjustedLangevinTK.C.

References queso_require, queso_require_equal_to, and queso_require_not_equal_to.

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

 {
   queso_require_not_equal_to(m_rvs.size(), 0);
   queso_require_equal_to(m_rvs.size(), m_scales.size());
 
   for (unsigned int i = 0; i < m_scales.size(); ++i) {
     double factor = 1./m_scales[i]/m_scales[i];
     queso_require(!(m_rvs[i]));
     m_rvs[i] = new GaussianVectorRV<V, M>(m_prefix.c_str(),
                                          *m_vectorSpace,
                                          m_vectorSpace->zeroVector(),
                                          factor*m_time_step*m_originalCovMatrix);
   }
 }

template<class V , class M >

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

virtual

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

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

Implements QUESO::BaseTKGroup< V, M >.

Definition at line 74 of file MetropolisAdjustedLangevinTK.C.

 {
   return false;
 }

template<class V , class M >

void QUESO::MetropolisAdjustedLangevinTK< V, M >::updateLawCovMatrix ( const M & covMatrix )

Scales the covariance matrix.

The covariance matrix is scaled by a factor of $ 1/scales^2 $ .

Definition at line 164 of file MetropolisAdjustedLangevinTK.C.

 {
   for (unsigned int i = 0; i < m_scales.size(); ++i) {
     double factor = 1./m_scales[i]/m_scales[i];
     if ((m_env.subDisplayFile()        ) &&
         (m_env.displayVerbosity() >= 10)) {
       *m_env.subDisplayFile() << "In MetropolisAdjustedLangevinTK<V, M>::updateLawCovMatrix()"
                               << ", m_scales.size() = " << m_scales.size()
                               << ", i = "               << i
                               << ", m_scales[i] = "     << m_scales[i]
                               << ", factor = "          << factor
                               << ": about to call m_rvs[i]->updateLawCovMatrix()"
                               << ", covMatrix = \n" << factor*covMatrix // FIX ME: might demand parallelism
                               << std::endl;
     }
     dynamic_cast<GaussianVectorRV<V, M> * >(m_rvs[i])->updateLawCovMatrix(factor*m_time_step*covMatrix);
   }
 }