diff --git a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h index cc0ae5422..7c2e46476 100644 --- a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h +++ b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h @@ -60,6 +60,32 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, const GslBlockMatrix & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodBlockDiagonalCovariance(const char * prefix, + const VectorSet & domainSet, + const V & observations, + const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodBlockDiagonalCovariance(); //@} @@ -70,12 +96,16 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase { //! Get (const) multiplicative coefficient for block \c i const double & getBlockCoefficient(unsigned int i) const; - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: std::vector m_covarianceCoefficients; const GslBlockMatrix & m_covariance; + + void checkDimConsistency(const V & observations, const GslBlockMatrix & covariance) const; }; } // End namespace QUESO diff --git a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h index afee9e955..426863092 100644 --- a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h +++ b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h @@ -62,24 +62,50 @@ class GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients : public Likel const VectorSet & domainSet, const V & observations, const GslBlockMatrix & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(const char * prefix, + const VectorSet & domainSet, + const V & observations, + const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(); //@} - //! Logarithm of the value of the scalar function. +protected: + + //! Logarithm of the value of the likelihood function. /*! * The last \c n elements of \c domainVector, where \c n is the number of * blocks in the block diagonal covariance matrix, are treated as * hyperparameters and will be sample in a statistical inversion. - * - * The user need not concern themselves with handling these in the model - * evaluation routine, since they are handled by the likelihood evaluation - * routine. */ - virtual double lnValue(const V & domainVector) const; + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: const GslBlockMatrix & m_covariance; + + void checkDimConsistency(const V & observations, const GslBlockMatrix & covariance) const; }; } // End namespace QUESO diff --git a/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h b/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h index 4c291a079..da7a6dcbb 100644 --- a/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h +++ b/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h @@ -54,12 +54,40 @@ class GaussianLikelihoodDiagonalCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, const V & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodDiagonalCovariance(const char * prefix, + const VectorSet & domainSet, + const V & observations, + const V & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodDiagonalCovariance(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: const V & m_covariance; diff --git a/src/stats/inc/GaussianLikelihoodFullCovariance.h b/src/stats/inc/GaussianLikelihoodFullCovariance.h index ff6edbf43..00051f312 100644 --- a/src/stats/inc/GaussianLikelihoodFullCovariance.h +++ b/src/stats/inc/GaussianLikelihoodFullCovariance.h @@ -58,12 +58,39 @@ class GaussianLikelihoodFullCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, const M & covariance, double covarianceCoefficient=1.0); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodFullCovariance(const char * prefix, + const VectorSet & domainSet, const V & observations, + const M & covariance, double covarianceCoefficient, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodFullCovariance(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: double m_covarianceCoefficient; diff --git a/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h b/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h index 2ed5c7ef1..0daef6882 100644 --- a/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h +++ b/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h @@ -63,12 +63,39 @@ class GaussianLikelihoodFullCovarianceRandomCoefficient : public LikelihoodBase< const VectorSet & domainSet, const V & observations, const M & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodFullCovarianceRandomCoefficient(const char * prefix, + const VectorSet & domainSet, const V & observations, + const M & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodFullCovarianceRandomCoefficient(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: const M & m_covariance; diff --git a/src/stats/inc/GaussianLikelihoodScalarCovariance.h b/src/stats/inc/GaussianLikelihoodScalarCovariance.h index ab5a459bd..83d134205 100644 --- a/src/stats/inc/GaussianLikelihoodScalarCovariance.h +++ b/src/stats/inc/GaussianLikelihoodScalarCovariance.h @@ -54,12 +54,39 @@ class GaussianLikelihoodScalarCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, double covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodScalarCovariance(const char * prefix, + const VectorSet & domainSet, const V & observations, + double covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodScalarCovariance(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: double m_covariance; diff --git a/src/stats/inc/LikelihoodBase.h b/src/stats/inc/LikelihoodBase.h index fb53754e3..20a39286c 100644 --- a/src/stats/inc/LikelihoodBase.h +++ b/src/stats/inc/LikelihoodBase.h @@ -28,9 +28,16 @@ #include #include #include +#include namespace QUESO { +template +class BaseVectorRV; + +template +class MultiDQuadratureBase; + class GslVector; class GslMatrix; @@ -59,8 +66,34 @@ class LikelihoodBase : public BaseScalarFunction { const VectorSet & domainSet, const V & observations); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + LikelihoodBase(const char * prefix, + const VectorSet & domainSet, + const V & observations, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor, pure to make this class abstract virtual ~LikelihoodBase() =0; + //@} //! Deprecated. Evaluates the user's model at the point \c domainVector @@ -79,7 +112,8 @@ class LikelihoodBase : public BaseScalarFunction { * be compared to actual observations when computing the likelihood functional. * * The first \c n components of \c domainVector are the model parameters. - * The rest of \c domainVector contains the hyperparameters, if any. For + * The rest of \c domainVector contains the hyperparameters, if any. Mathematically, this + * function is \f$ f(m) \f$. For * example, in \c GaussianLikelihoodFullCovarainceRandomCoefficient, the last * component of \c domainVector contains the multiplicative coefficient of * the observational covariance matrix. In this case, the user need not @@ -94,13 +128,51 @@ class LikelihoodBase : public BaseScalarFunction { virtual void evaluateModel(const V & domainVector, V & modelOutput) const { this->evaluateModel(domainVector,NULL,modelOutput,NULL,NULL,NULL); } + //! Evaluates the user's model at the point \c (domainVector,marginalVector) + /*! + * Subclass implementations will fill up the \c modelOutput vector with output + * from the model. This function will be passed the current value of the parameters in domainVector + * and the current values of the marginalization parameters in marginalVector. Mathematically, this + * function is \f$ f(m,q) \f$. This function + * is only called if the likelihood includes marginalization. Otherwise, + * evaluateModel(const V & domainVector, V & modelOutput) will be called by lnValue(). + * + * Not every user will do marginalization, so this is not a pure function, but there is + * no default, so we error if marginalization is attempted without overriding this function. + */ + virtual void evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput) const; + //! Actual value of the scalar function. virtual double actualValue(const V & domainVector, const V * /*domainDirection*/, V * /*gradVector*/, M * /*hessianMatrix*/, V * /*hessianEffect*/) const { return std::exp(this->lnValue(domainVector)); } + //! Logarithm of the value of the scalar function. + /*! + * This method well evaluate the users model (evaluateModel()) and then + * call lnLikelihood() and, if there is marginalization, will handle the + * numerical integration over the marginal parameter space. + */ + virtual double lnValue(const V & domainVector) const; + protected: + const V & m_observations; + + typename SharedPtr >::Type m_marg_param_pdf; + + typename SharedPtr >::Type m_marg_integration; + + bool m_marg_pdf_is_weight_func; + + //! Compute log-likelihood value given the current parameters and model output + /*! + * Subclasses should override this method to compute the likelihood value, given + * the current parameters in domainVector and the output from the model evaluated + * at those parameter values. Note that this function may alter the values of + * modelOutput to reduce copying. + */ + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const =0; }; } // End namespace QUESO diff --git a/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C b/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C index ee3813a30..8e032ebba 100644 --- a/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C +++ b/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C @@ -39,15 +39,21 @@ GaussianLikelihoodBlockDiagonalCovariance::GaussianLikelihoodBlockDiagonal m_covarianceCoefficients(covariance.numBlocks(), 1.0), m_covariance(covariance) { - unsigned int totalDim = 0; - - for (unsigned int i = 0; i < this->m_covariance.numBlocks(); i++) { - totalDim += this->m_covariance.getBlock(i).numRowsLocal(); - } + this->checkDimConsistency(observations,m_covariance); +} - if (totalDim != observations.sizeLocal()) { - queso_error_msg("Covariance matrix not same dimension as observation vector"); - } +template +GaussianLikelihoodBlockDiagonalCovariance::GaussianLikelihoodBlockDiagonalCovariance( + const char * prefix, const VectorSet & domainSet, + const V & observations, const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covarianceCoefficients(covariance.numBlocks(), 1.0), + m_covariance(covariance) +{ + this->checkDimConsistency(observations,m_covariance); } template @@ -73,13 +79,10 @@ GaussianLikelihoodBlockDiagonalCovariance::getBlockCoefficient( template double -GaussianLikelihoodBlockDiagonalCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodBlockDiagonalCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; @@ -108,6 +111,22 @@ GaussianLikelihoodBlockDiagonalCovariance::lnValue(const V & domainVector) return -0.5 * norm2_squared; } +template +void +GaussianLikelihoodBlockDiagonalCovariance::checkDimConsistency( + const V & observations, const GslBlockMatrix & covariance) const +{ + unsigned int totalDim = 0; + + for (unsigned int i = 0; i < covariance.numBlocks(); i++) { + totalDim += covariance.getBlock(i).numRowsLocal(); + } + + if (totalDim != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same dimension as observation vector"); + } +} + } // End namespace QUESO template class QUESO::GaussianLikelihoodBlockDiagonalCovariance; diff --git a/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C b/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C index 3933392fb..8a959ad90 100644 --- a/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C +++ b/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C @@ -36,15 +36,20 @@ GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::GaussianLikel : LikelihoodBase(prefix, domainSet, observations), m_covariance(covariance) { - unsigned int totalDim = 0; - - for (unsigned int i = 0; i < this->m_covariance.numBlocks(); i++) { - totalDim += this->m_covariance.getBlock(i).numRowsLocal(); - } + this->checkDimConsistency(observations,m_covariance); +} - if (totalDim != observations.sizeLocal()) { - queso_error_msg("Covariance matrix not same dimension as observation vector"); - } +template +GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients( + const char * prefix, const VectorSet & domainSet, + const V & observations, const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ + this->checkDimConsistency(observations,m_covariance); } template @@ -54,13 +59,11 @@ GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::~GaussianLike template double -GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::lnValue(const V & domainVector) const +GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::lnLikelihood(const V & domainVector, + V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; @@ -105,6 +108,22 @@ GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::lnValue(const return -0.5 * norm2_squared - cov_norm_factor; } +template +void +GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::checkDimConsistency( + const V & observations, const GslBlockMatrix & covariance) const +{ + unsigned int totalDim = 0; + + for (unsigned int i = 0; i < covariance.numBlocks(); i++) { + totalDim += covariance.getBlock(i).numRowsLocal(); + } + + if (totalDim != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same dimension as observation vector"); + } +} + } // End namespace QUESO template class QUESO::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients; diff --git a/src/stats/src/GaussianLikelihoodDiagonalCovariance.C b/src/stats/src/GaussianLikelihoodDiagonalCovariance.C index 639e45d1d..c93d31023 100644 --- a/src/stats/src/GaussianLikelihoodDiagonalCovariance.C +++ b/src/stats/src/GaussianLikelihoodDiagonalCovariance.C @@ -43,6 +43,21 @@ GaussianLikelihoodDiagonalCovariance::GaussianLikelihoodDiagonalCovariance } } +template +GaussianLikelihoodDiagonalCovariance::GaussianLikelihoodDiagonalCovariance( + const char * prefix, const VectorSet & domainSet, + const V & observations, const V & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ + if (covariance.sizeLocal() != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same size as observation vector"); + } +} + template GaussianLikelihoodDiagonalCovariance::~GaussianLikelihoodDiagonalCovariance() { @@ -50,12 +65,8 @@ GaussianLikelihoodDiagonalCovariance::~GaussianLikelihoodDiagonalCovarianc template double -GaussianLikelihoodDiagonalCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodDiagonalCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy - - this->evaluateModel(domainVector, modelOutput); - modelOutput -= this->m_observations; // Compute misfit modelOutput *= modelOutput; modelOutput /= this->m_covariance; // Multiply by inverse covriance matrix diff --git a/src/stats/src/GaussianLikelihoodFullCovariance.C b/src/stats/src/GaussianLikelihoodFullCovariance.C index 3e8a0bd2f..76cb70126 100644 --- a/src/stats/src/GaussianLikelihoodFullCovariance.C +++ b/src/stats/src/GaussianLikelihoodFullCovariance.C @@ -44,6 +44,22 @@ GaussianLikelihoodFullCovariance::GaussianLikelihoodFullCovariance( } } +template +GaussianLikelihoodFullCovariance::GaussianLikelihoodFullCovariance(const char * prefix, + const VectorSet & domainSet, const V & observations, + const M & covariance, double covarianceCoefficient, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covarianceCoefficient(covarianceCoefficient), + m_covariance(covariance) +{ + if (covariance.numRowsLocal() != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same size as observation vector"); + } +} + template GaussianLikelihoodFullCovariance::~GaussianLikelihoodFullCovariance() { @@ -51,13 +67,10 @@ GaussianLikelihoodFullCovariance::~GaussianLikelihoodFullCovariance() template double -GaussianLikelihoodFullCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodFullCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; diff --git a/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C b/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C index 349c51a5c..7b858a535 100644 --- a/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C +++ b/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C @@ -42,6 +42,21 @@ GaussianLikelihoodFullCovarianceRandomCoefficient::GaussianLikelihoodFullC } } +template +GaussianLikelihoodFullCovarianceRandomCoefficient::GaussianLikelihoodFullCovarianceRandomCoefficient( + const char * prefix, const VectorSet & domainSet, + const V & observations, const M & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ + if (covariance.numRowsLocal() != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same size as observation vector"); + } +} + template GaussianLikelihoodFullCovarianceRandomCoefficient::~GaussianLikelihoodFullCovarianceRandomCoefficient() { @@ -49,13 +64,10 @@ GaussianLikelihoodFullCovarianceRandomCoefficient::~GaussianLikelihoodFull template double -GaussianLikelihoodFullCovarianceRandomCoefficient::lnValue(const V & domainVector) const +GaussianLikelihoodFullCovarianceRandomCoefficient::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; diff --git a/src/stats/src/GaussianLikelihoodScalarCovariance.C b/src/stats/src/GaussianLikelihoodScalarCovariance.C index a6ca88e82..007af6094 100644 --- a/src/stats/src/GaussianLikelihoodScalarCovariance.C +++ b/src/stats/src/GaussianLikelihoodScalarCovariance.C @@ -40,6 +40,18 @@ GaussianLikelihoodScalarCovariance::GaussianLikelihoodScalarCovariance( { } +template +GaussianLikelihoodScalarCovariance::GaussianLikelihoodScalarCovariance( + const char * prefix, const VectorSet & domainSet, + const V & observations, double covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ +} + template GaussianLikelihoodScalarCovariance::~GaussianLikelihoodScalarCovariance() { @@ -47,12 +59,8 @@ GaussianLikelihoodScalarCovariance::~GaussianLikelihoodScalarCovariance() template double -GaussianLikelihoodScalarCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodScalarCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy - - this->evaluateModel(domainVector, modelOutput); - modelOutput -= this->m_observations; // Compute misfit double norm2_squared = modelOutput.norm2Sq(); // Compute square of 2-norm diff --git a/src/stats/src/LikelihoodBase.C b/src/stats/src/LikelihoodBase.C index f81a0271e..d45ff1a0e 100644 --- a/src/stats/src/LikelihoodBase.C +++ b/src/stats/src/LikelihoodBase.C @@ -28,6 +28,8 @@ #include #include #include +#include +#include namespace QUESO { @@ -40,6 +42,25 @@ LikelihoodBase::LikelihoodBase( { } +template +LikelihoodBase::LikelihoodBase(const char * prefix, + const VectorSet & domainSet, + const V & observations, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : BaseScalarFunction(prefix, domainSet), + m_observations(observations), + m_marg_param_pdf(marg_param_pdf), + m_marg_integration(marg_integration), + m_marg_pdf_is_weight_func(marg_pdf_is_weight_func) +{ + // The dimension of the parameter space had better match the dimension in the integration + queso_require_equal_to_msg(this->m_marg_param_pdf->imageSet().vectorSpace().dimGlobal(), + this->m_marg_integration->getDomain().vectorSpace().dimGlobal(), + "Mismatched marginal parameter space dimension and quadrature dimension!"); +} + template LikelihoodBase::~LikelihoodBase() { @@ -59,6 +80,74 @@ void LikelihoodBase::evaluateModel(const V & domainVector, const V * domai queso_error_msg(ss.str()); } +template +void LikelihoodBase::evaluateModel(const V & domainVector, + const V & marginalVector, + V & modelOutput) const +{ + std::stringstream ss; + ss << "ERROR: evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput)" + << std::endl + << " is not implemented! Please override this function in your subclass." + << std::endl; + + queso_error_msg(ss.str()); +} + +template +double +LikelihoodBase::lnValue(const V & domainVector) const +{ + V modelOutput(this->m_observations, 0, 0); + + double lnLikelihood_value = 0.0; + + // If we're not marginalizing, then we just need f(m) + if( !m_marg_param_pdf ) + { + this->evaluateModel(domainVector, modelOutput); + + // Note modelOutput made be modified in lnLikelihood() + lnLikelihood_value = this->lnLikelihood(domainVector, modelOutput); + } + // Otherwise we're integrating over marginal parameter space + else + { + queso_assert(m_marg_integration); + + const std::vector< typename QUESO::SharedPtr< V >::Type > & x = + this->m_marg_integration->positions(); + + const std::vector< double > & w = this->m_marg_integration->weights(); + + unsigned int n_qpoints = x.size(); + + for( unsigned int q = 0; q < n_qpoints; q++ ) + { + this->evaluateModel(domainVector, *(x[q]), modelOutput); + + // Note modelOutput made be modified in lnLikelihood() + double ln_pi_m_q = this->lnLikelihood(domainVector, modelOutput); + + double ln_pi_q = 0.0; + + // If we're not treating the marginal parameters pdf as + // the weighting function in the quadrature evaluation, + // then we need to also evaluate it. + if(!m_marg_pdf_is_weight_func) + this->m_marg_param_pdf->pdf().lnValue( *(x[q]) ); + + /*! \todo [PB]: We might want to play games with a log(\sum) identity + if precision starts to become an issue. */ + lnLikelihood_value += std::exp(ln_pi_m_q + ln_pi_q)*w[q]; + } + + lnLikelihood_value = std::log(lnLikelihood_value); + } + + return lnLikelihood_value; +} + } // End namespace QUESO template class QUESO::LikelihoodBase; diff --git a/test/Makefile.am b/test/Makefile.am index 68ff4c4b6..d4e6e6516 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -132,6 +132,7 @@ unit_driver_SOURCES += unit/rng_cxx11.C unit_driver_SOURCES += unit/rng_boost.C unit_driver_SOURCES += unit/basic_pdfs_cxx11.C unit_driver_SOURCES += unit/basic_pdfs_boost.C +unit_driver_SOURCES += unit/marg_likelihood.C test_boxsubset_centroid_SOURCES = test_centroids/test_boxsubset_centroid.C test_concatenation_centroid_SOURCES = test_centroids/test_concatenation_centroid.C diff --git a/test/unit/marg_likelihood.C b/test/unit/marg_likelihood.C new file mode 100644 index 000000000..d67c701f4 --- /dev/null +++ b/test/unit/marg_likelihood.C @@ -0,0 +1,327 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2017 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include "config_queso.h" + +#ifdef QUESO_HAVE_CPPUNIT + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace QUESOTesting +{ + template + class TestlingLikelihoodBase : public QUESO::LikelihoodBase + { + public: + + TestlingLikelihoodBase(const char * prefix, + const QUESO::VectorSet & domainSet, + const V & observations, + typename QUESO::SharedPtr >::Type & marg_param_pdf, + typename QUESO::SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : QUESO::LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func) + {} + + //! Evaluate the exact ln value of the marginalized likelihood + virtual double lnMargLikelihood( const V & domainVector ) const =0; + + void testLikelihoodValue( const V & domainVector, + const QUESO::LikelihoodBase & likelihood, + double tol ) const + { + double computed_value = likelihood.lnValue(domainVector); + + double exact_value = this->lnMargLikelihood(domainVector); + + double rel_error = std::abs( (computed_value-exact_value)/exact_value ); + + CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, rel_error, tol ); + } + + protected: + + virtual double lnLikelihood(const V & /*domainVector*/, + V & modelOutput) const + { + queso_require_equal_to(1,modelOutput.sizeGlobal()); + return std::log(modelOutput[0]); + } + + }; + + template + class LinearModelLikelihood : public TestlingLikelihoodBase + { + public: + + LinearModelLikelihood(const char * prefix, + const QUESO::VectorSet & domainSet, + const V & observations, + typename QUESO::SharedPtr >::Type & marg_param_pdf, + typename QUESO::SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : TestlingLikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func) + {} + + virtual void evaluateModel(const V & domainVector, + const V & marginalVector, + V & modelOutput) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + queso_require_equal_to(3,marginalVector.sizeGlobal()); + queso_require_equal_to(1,modelOutput.sizeGlobal()); + + modelOutput[0] = 3.14*domainVector[0] + 1.1*marginalVector[0] + + 2.2*marginalVector[1] + 3.3*marginalVector[2]; + } + + virtual double lnMargLikelihood( const V & domainVector ) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + + return std::log(3.14*domainVector[0] + 1.1/2.0 + 2.2/2.0 + 3.3/2.0); + } + + }; + + template + class GaussModelLikelihood : public TestlingLikelihoodBase + { + public: + + GaussModelLikelihood(const char * prefix, + const QUESO::VectorSet & domainSet, + const V & observations, + typename QUESO::SharedPtr >::Type & marg_param_pdf, + typename QUESO::SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : TestlingLikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func) + {} + + virtual void evaluateModel(const V & domainVector, + const V & marginalVector, + V & modelOutput) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + queso_require_equal_to(2,marginalVector.sizeGlobal()); + queso_require_equal_to(1,modelOutput.sizeGlobal()); + + double q1 = marginalVector[0]; + double q2 = marginalVector[1]; + double m = domainVector[0]; + + modelOutput[0] = 3.14*m + 1.1*q1 + 2.2*q2*q2; + } + + virtual double lnMargLikelihood( const V & domainVector ) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + + double m = domainVector[0]; + + return std::log(3.14*m*M_PI + 2.2*M_PI/2.0); + } + + }; + + template + class MarginalLikelihoodTestBase : public CppUnit::TestCase + { + public: + + void setUp() + { + this->init_env(); + } + + void test_linear_func_uniform_marg_space() + { + // Instantiate the parameter space + unsigned int param_dim = 1; + unsigned int marg_dim = 3; + QUESO::VectorSpace param_space( (*this->_env), "param_", param_dim, NULL); + QUESO::VectorSpace marg_space( (*this->_env), "marg_", marg_dim, NULL); + + double param_min_domain_value = 0.0; + double param_max_domain_value = 1.0; + + double marg_min_domain_value = 0.0; + double marg_max_domain_value = 1.0; + + typename QUESO::ScopedPtr::Type param_min_values( param_space.newVector(param_min_domain_value) ); + typename QUESO::ScopedPtr::Type param_max_values( param_space.newVector(param_max_domain_value) ); + + typename QUESO::ScopedPtr::Type marg_min_values( marg_space.newVector(marg_min_domain_value) ); + typename QUESO::ScopedPtr::Type marg_max_values( marg_space.newVector(marg_max_domain_value) ); + + QUESO::BoxSubset param_domain( "param_domain_", param_space, (*param_min_values), (*param_max_values) ); + QUESO::BoxSubset marg_domain( "marg_domain_", marg_space, (*marg_min_values), (*marg_max_values) ); + + typename QUESO::SharedPtr >::Type + marg_param_rv( new QUESO::UniformVectorRV("marg_param_rv_", marg_domain) ); + + const V & data = param_space.zeroVector(); + + unsigned int int_order = 1; + + QUESO::SharedPtr::Type + qrule_1d( new QUESO::UniformLegendre1DQuadrature(marg_min_domain_value, marg_max_domain_value, + int_order, false) ); + + std::vector::Type> all_1d_qrules(marg_dim,qrule_1d); + + typename QUESO::SharedPtr >::Type + marg_integration( new QUESO::TensorProductQuadrature(marg_domain,all_1d_qrules) ); + + LinearModelLikelihood likelihood( "likelihood_test_", + param_domain, + data, + marg_param_rv, + marg_integration, + false ); + + // Test marginalized likelihood over a range of points in the parameter domain + unsigned int n_intervals = 20; + double tol = std::numeric_limits::epsilon()*10; + + this->test_likelihood_values_range( n_intervals, param_min_domain_value, param_max_domain_value, + param_space, tol, likelihood ); + } + + void test_linear_func_gaussian_marg_space() + { + // Instantiate the parameter space + unsigned int param_dim = 1; + unsigned int marg_dim = 2; + QUESO::VectorSpace param_space( (*this->_env), "param_", param_dim, NULL); + QUESO::VectorSpace marg_space( (*this->_env), "marg_", marg_dim, NULL); + + double param_min_domain_value = 0.0; + double param_max_domain_value = 1.0; + + double marg_min_domain_value = -INFINITY; + double marg_max_domain_value = INFINITY; + + typename QUESO::ScopedPtr::Type param_min_values( param_space.newVector(param_min_domain_value) ); + typename QUESO::ScopedPtr::Type param_max_values( param_space.newVector(param_max_domain_value) ); + + typename QUESO::ScopedPtr::Type marg_min_values( marg_space.newVector(marg_min_domain_value) ); + typename QUESO::ScopedPtr::Type marg_max_values( marg_space.newVector(marg_max_domain_value) ); + + QUESO::BoxSubset param_domain( "param_domain_", param_space, (*param_min_values), (*param_max_values) ); + QUESO::BoxSubset marg_domain( "marg_domain_", marg_space, (*marg_min_values), (*marg_max_values) ); + + typename QUESO::ScopedPtr::Type var_vec( param_space.newVector(1.0) ); + + // Zero mean, unit variance Gaussian + typename QUESO::SharedPtr >::Type + marg_param_rv( new QUESO::GaussianVectorRV("marg_param_rv_", marg_domain, + marg_space.zeroVector(), + *var_vec) ); + + const V & data = param_space.zeroVector(); + + unsigned int int_order = 1; + + QUESO::SharedPtr::Type + qrule_1d( new QUESO::GaussianHermite1DQuadrature(0.0,1.0,int_order) ); + + std::vector::Type> all_1d_qrules(marg_dim,qrule_1d); + + typename QUESO::SharedPtr >::Type + marg_integration( new QUESO::TensorProductQuadrature(marg_domain,all_1d_qrules) ); + + GaussModelLikelihood likelihood( "likelihood_test_", + param_domain, + data, + marg_param_rv, + marg_integration, + true ); + + // Test marginalized likelihood over a range of points in the parameter domain + unsigned int n_intervals = 20; + double tol = std::numeric_limits::epsilon()*50; + + this->test_likelihood_values_range( n_intervals, param_min_domain_value, param_max_domain_value, + param_space, tol, likelihood ); + } + + protected: + + QUESO::EnvOptionsValues _options; + + typename QUESO::ScopedPtr::Type _env; + + void init_env() + { + _env.reset( new QUESO::FullEnvironment("","",&_options) ); + } + + void test_likelihood_values_range( unsigned int n_intervals, double param_min_domain_value, + double param_max_domain_value, const QUESO::VectorSpace & param_space, + double tol, const TestlingLikelihoodBase & likelihood ) + { + for( unsigned int i = 0; i < n_intervals; i++ ) + { + double x = param_min_domain_value + i*(param_max_domain_value-param_min_domain_value)/(double)n_intervals; + + typename QUESO::ScopedPtr::Type param_vec( param_space.newVector(x) ); + + likelihood.testLikelihoodValue(*param_vec,likelihood,tol); + } + } + }; + + class MarginalLikelihoodGslTest : public MarginalLikelihoodTestBase + { + public: + CPPUNIT_TEST_SUITE( MarginalLikelihoodGslTest ); + + CPPUNIT_TEST( test_linear_func_uniform_marg_space ); + CPPUNIT_TEST( test_linear_func_gaussian_marg_space ); + + CPPUNIT_TEST_SUITE_END(); + }; + + CPPUNIT_TEST_SUITE_REGISTRATION( MarginalLikelihoodGslTest ); + +} // end namespace QUESOTesting + +#endif // QUESO_HAVE_CPPUNIT diff --git a/test/unit/quadrature_1d.C b/test/unit/quadrature_1d.C index 75265ebe6..82826d7b5 100644 --- a/test/unit/quadrature_1d.C +++ b/test/unit/quadrature_1d.C @@ -166,7 +166,6 @@ namespace QUESOTesting void test_1d_hermite_quadrature_x4erf() { - std::cout << std::endl << "Beginning Hermite X4Erf test!" << std::endl; X4Erf func; std::vector testing_orders; this->testing_orders(testing_orders);