Add new LikelihoodBase constructor for marginal params

This is laying the groundwork for having the likelihood do
marginalization for the user. This add the neccessary data and
constructor to LikelihoodBase. In particular, we'll need the joint
pdf of the marginal parameter space and the user-constructed
integration rule for integrating over that space.

src/stats/inc/LikelihoodBase.h

@@ -28,9 +28,16 @@
 #include <vector>
 #include <cmath>
 #include <queso/ScalarFunction.h>
+#include <queso/SharedPtr.h>
 
 namespace QUESO {
 
+template<class V, class M>
+class BaseVectorRV;
+
+template<class V, class M>
+class MultiDQuadratureBase;
+
 class GslVector;
 class GslMatrix;
 
@@ -59,8 +66,34 @@ class LikelihoodBase : public BaseScalarFunction<V, M> {
       const VectorSet<V, M> & 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<V, M> & domainSet,
+                 const V & observations,
+                 typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+                 typename SharedPtr<MultiDQuadratureBase<V,M> >::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
@@ -99,12 +132,24 @@ class LikelihoodBase : public BaseScalarFunction<V, M> {
                              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<const BaseVectorRV<V,M> >::Type m_marg_param_pdf;
+
+  typename SharedPtr<MultiDQuadratureBase<V,M> >::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

src/stats/src/LikelihoodBase.C

@@ -28,6 +28,8 @@
 #include <queso/GslMatrix.h>
 #include <queso/VectorSet.h>
 #include <queso/LikelihoodBase.h>
+#include <queso/VectorRV.h>
+#include <queso/MultiDQuadratureBase.h>
 
 namespace QUESO {
 
@@ -40,6 +42,25 @@ LikelihoodBase<V, M>::LikelihoodBase(
 {
 }
 
+template<class V, class M>
+LikelihoodBase<V, M>::LikelihoodBase(const char * prefix,
+                                     const VectorSet<V, M> & domainSet,
+                                     const V & observations,
+                                     typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+                                     typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+                                     bool marg_pdf_is_weight_func)
+  : BaseScalarFunction<V, M>(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<class V, class M>
 LikelihoodBase<V, M>::~LikelihoodBase()
 {