Replace ad-hoc DiagonalWrapper with built-in Eigen approach

As far as I can tell, this is mathematically and operationally
equivalent, but using Eigen's built-in functionality for producing a
diagonal matrix expression from a 1-d array solves issue
#419.

include/albatross/src/core/transformed_distribution.hpp

@@ -36,17 +36,12 @@ namespace albatross {
 
 namespace details {
 
-template <typename X> inline auto diagonal_wrapper(X &x) {
-  return Eigen::DiagonalWrapper<X>(x);
-}
-
 template <typename MatrixType, typename DiagType>
 inline Eigen::MatrixXd
 product_sqrt(const Eigen::MatrixBase<MatrixType> &matrix,
              const Eigen::DiagonalBase<DiagType> &diag_matrix) {
-
   return matrix.derived() *
-         diagonal_wrapper(diag_matrix.diagonal().array().sqrt().eval());
+         diag_matrix.diagonal().array().sqrt().matrix().asDiagonal();
 }
 
 template <typename Lhs, typename Rhs>
@@ -61,7 +56,7 @@ inline Eigen::MatrixXd
 product_sqrt(const Eigen::SparseMatrixBase<MatrixType> &matrix,
              const Eigen::DiagonalBase<DiagType> &diag_matrix) {
   return matrix.derived() *
-         diagonal_wrapper(diag_matrix.diagonal().array().sqrt().eval());
+         diag_matrix.diagonal().array().sqrt().matrix().asDiagonal();
 }
 
 template <typename SparseType, typename MatrixType>