Refactor logit transform logic into TK

src/stats/inc/HessianCovMatricesTKGroup.h

@@ -82,6 +82,8 @@ class HessianCovMatricesTKGroup : public BaseTKGroup<V,M> {
 
   virtual bool covMatrixIsDirty() { return false; }
   virtual void cleanCovMatrix() { }
+
+  virtual void updateLawCovMatrix(const M & covMatrix);
   //@}
 
   //! @name I/O methods

src/stats/inc/MetropolisAdjustedLangevinTK.h

@@ -75,7 +75,7 @@ class MetropolisAdjustedLangevinTK : public BaseTKGroup<V, M> {
 
   //! Scales the covariance matrix.
   /*! The covariance matrix is scaled by a factor of \f$ 1/scales^2 \f$.*/
-  void updateLawCovMatrix(const M & covMatrix);
+  virtual void updateLawCovMatrix(const M & covMatrix);
   //@}
 
   //! @name Misc methods

src/stats/inc/MetropolisHastingsSG.h

@@ -341,9 +341,6 @@ class MetropolisHastingsSG
 	double m_initialLogPriorValue;
 	double m_initialLogLikelihoodValue;
 
-  void transformInitialCovMatrixToGaussianSpace(const BoxSubset<P_V, P_M> &
-      boxSubset);
-
   bool m_userDidNotProvideOptions;
 
   unsigned int m_latestDirtyCovMatrixIteration;

src/stats/inc/ScaledCovMatrixTKGroup.h

@@ -70,7 +70,7 @@ class ScaledCovMatrixTKGroup : public BaseTKGroup<V,M> {
 
   //! Scales the covariance matrix.
   /*! The covariance matrix is scaled by a factor of \f$ 1/scales^2 \f$.*/
-  void                          updateLawCovMatrix        (const M& covMatrix);
+  virtual void updateLawCovMatrix(const M & covMatrix);
   //@}
 
   //! @name Misc methods

src/stats/inc/TKGroup.h

@@ -115,6 +115,14 @@ class BaseTKGroup {
    * After this method is called, covMatrixIsDirty should return \c true.
    */
   virtual void cleanCovMatrix() = 0;
+
+  //! Sets and proposal covariance matrix (and for all DR stage RVs as well)
+  /*!
+   * If the TK has a proposal covariance matrix, set the covariance matrix.
+   * If the TK support delayed rejection, derived classes should scale by a
+   * factor of \f$ 1/scales^2 \f$.
+   */
+  virtual void updateLawCovMatrix(const M & covMatrix) = 0;
   //@}
 
   //! @name I/O methods

src/stats/inc/TransformedScaledCovMatrixTKGroup.h

@@ -74,7 +74,7 @@ class TransformedScaledCovMatrixTKGroup : public BaseTKGroup<V, M> {
 
   //! Scales the covariance matrix of the underlying Gaussian distribution.
   /*! The covariance matrix is scaled by a factor of \f$ 1/scales^2 \f$.*/
-  void updateLawCovMatrix(const M & covMatrix);
+  virtual void updateLawCovMatrix(const M & covMatrix);
   //@}
 
   //! @name Misc methods
@@ -113,6 +113,8 @@ class TransformedScaledCovMatrixTKGroup : public BaseTKGroup<V, M> {
   //! Sets the mean of the underlying Gaussian RVs to zero.
   void setRVsWithZeroMean();
 
+  void transformCovMatrixToGaussianSpace(M & covMatrix);
+
   using BaseTKGroup<V, M>::m_env;
   using BaseTKGroup<V, M>::m_prefix;
   using BaseTKGroup<V, M>::m_vectorSpace;

src/stats/src/HessianCovMatricesTKGroup.C

@@ -372,6 +372,13 @@ HessianCovMatricesTKGroup<V,M>::clearPreComputingPositions()
   return;
 }
 
+template <class V, class M>
+void
+HessianCovMatricesTKGroup<V, M>::updateLawCovMatrix(const M & covMatrix)
+{
+  queso_error_msg(std::string("Stochastic Newton does not support adaptivity"));
+}
+
 template <class V, class M>
 unsigned int
 HessianCovMatricesTKGroup<V, M>::set_dr_stage(unsigned int stageId)

src/stats/src/MetropolisHastingsSG.C

@@ -525,17 +525,6 @@ MetropolisHastingsSG<P_V,P_M>::commonConstructor()
     queso_require_msg(!(m_nullInputProposalCovMatrix), "proposal cov matrix should have been passed by user, since, according to the input algorithm options, local Hessians will not be used in the proposal");
   }
 
-  // Only transform prop cov matrix if we're doing a logit random walk.
-  // Also note we're transforming *after* we potentially read it from the input
-  // file.
-  if ((m_optionsObj->m_algorithm == "logit_random_walk") &&
-      (m_optionsObj->m_tk        == "logit_random_walk") &&
-      m_optionsObj->m_doLogitTransform) {
-    // Variable transform initial proposal cov matrix
-    transformInitialCovMatrixToGaussianSpace(
-        dynamic_cast<const BoxSubset<P_V, P_M> & >(m_targetPdf.domainSet()));
-  }
-
   // This instantiates all the transition kernels with their associated
   // factories
   TKFactoryInitializer tk_factory_initializer;
@@ -2243,20 +2232,7 @@ MetropolisHastingsSG<P_V, P_M>::adapt(unsigned int positionId,
       }
     }
 
-    // Transform the proposal covariance matrix if we have Logit transforms
-    // turned on.
-    //
-    // This logic should really be done inside the kernel itself
-    if (this->m_optionsObj->m_tk == "logit_random_walk") {
-      (dynamic_cast<TransformedScaledCovMatrixTKGroup<P_V,P_M>* >(m_tk.get()))
-        ->updateLawCovMatrix(tmpMatrix);
-    }
-    else {
-      // Assume that if we're not doing a logit transform, we're doing a random
-      // something sensible
-      (dynamic_cast<ScaledCovMatrixTKGroup<P_V,P_M>* >(m_tk.get()))
-        ->updateLawCovMatrix(tmpMatrix);
-    }
+    m_tk->updateLawCovMatrix(tmpMatrix);
 
 #ifdef UQ_DRAM_MCG_REQUIRES_INVERTED_COV_MATRICES
     queso_require_msg(!(UQ_INCOMPLETE_IMPLEMENTATION_RC), "need to code the update of m_upperCholProposalPrecMatrices");
@@ -2527,51 +2503,6 @@ MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix(
   return;
 }
 
-template <class P_V, class P_M>
-void
-MetropolisHastingsSG<P_V, P_M>::transformInitialCovMatrixToGaussianSpace(
-    const BoxSubset<P_V, P_M> & boxSubset)
-{
-  P_V min_domain_bounds(boxSubset.minValues());
-  P_V max_domain_bounds(boxSubset.maxValues());
-
-  for (unsigned int i = 0; i < min_domain_bounds.sizeLocal(); i++) {
-    double min_val = min_domain_bounds[i];
-    double max_val = max_domain_bounds[i];
-
-    if (queso_isfinite(min_val) && queso_isfinite(max_val)) {
-      if (m_initialProposalCovMatrix(i, i) >= max_val - min_val) {
-        // User is trying to specify a uniform proposal distribution, which
-        // is unsupported.  Throw an error for now.
-        std::cerr << "Proposal variance element "
-                  << i
-                  << " is "
-                  << m_initialProposalCovMatrix(i, i)
-                  << " but domain is of size "
-                  << max_val - min_val
-                  << std::endl;
-        std::cerr << "QUESO does not support uniform-like proposal "
-                  << "distributions.  Try making the proposal variance smaller"
-                  << std::endl;
-      }
-
-      // The jacobian at the midpoint of the domain
-      double transformJacobian = 4.0 / (max_val - min_val);
-
-      // Just do the multiplication by hand for now.  There's no method in
-      // Gsl(Vector|Matrix) to do this for me.
-      for (unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
-        // Multiply column j by element j
-        m_initialProposalCovMatrix(j, i) *= transformJacobian;
-      }
-      for (unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
-        // Multiply row j by element j
-        m_initialProposalCovMatrix(i, j) *= transformJacobian;
-      }
-    }
-  }
-}
-
 template class MetropolisHastingsSG<GslVector, GslMatrix>;
 
 }  // End namespace QUESO

src/stats/src/TransformedScaledCovMatrixTKGroup.C

@@ -53,6 +53,10 @@ TransformedScaledCovMatrixTKGroup<V,M>::TransformedScaledCovMatrixTKGroup(
                            << std::endl;
   }
 
+  // Transform prop cov matrix since we're doing a logit random walk.
+  // Note we're transforming *after* we potentially read it from the input file.
+  transformCovMatrixToGaussianSpace(m_originalCovMatrix);
+
   // Set RVs to have zero mean in the Gaussian space
   setRVsWithZeroMean();
 
@@ -310,6 +314,51 @@ TransformedScaledCovMatrixTKGroup<V, M>::transformToGaussianSpace(
   }
 }
 
+template <class V, class M>
+void
+TransformedScaledCovMatrixTKGroup<V, M>::transformCovMatrixToGaussianSpace(
+    M & covMatrix)
+{
+  V min_domain_bounds(m_boxSubset.minValues());
+  V max_domain_bounds(m_boxSubset.maxValues());
+
+  for (unsigned int i = 0; i < min_domain_bounds.sizeLocal(); i++) {
+    double min_val = min_domain_bounds[i];
+    double max_val = max_domain_bounds[i];
+
+    if (queso_isfinite(min_val) && queso_isfinite(max_val)) {
+      if (covMatrix(i, i) >= max_val - min_val) {
+        // User is trying to specify a uniform proposal distribution, which
+        // is unsupported.  Throw an error for now.
+        std::cerr << "Proposal variance element "
+                  << i
+                  << " is "
+                  << covMatrix(i, i)
+                  << " but domain is of size "
+                  << max_val - min_val
+                  << std::endl;
+        std::cerr << "QUESO does not support uniform-like proposal "
+                  << "distributions.  Try making the proposal variance smaller"
+                  << std::endl;
+      }
+
+      // The jacobian at the midpoint of the domain
+      double transformJacobian = 4.0 / (max_val - min_val);
+
+      // Just do the multiplication by hand for now.  There's no method in
+      // Gsl(Vector|Matrix) to do this for me.
+      for (unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
+        // Multiply column j by element j
+        covMatrix(j, i) *= transformJacobian;
+      }
+      for (unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
+        // Multiply row j by element j
+        covMatrix(i, j) *= transformJacobian;
+      }
+    }
+  }
+}
+
 }  // End namespace QUESO
 
 template class QUESO::TransformedScaledCovMatrixTKGroup<QUESO::GslVector, QUESO::GslMatrix>;