Edit according to the Ryan's review

CMakeLists.txt

@@ -413,9 +413,7 @@ add_definitions(-DBOOST_TEST_DYN_LINK)
 
 # Detect OpenMP support in a compiler. If the compiler supports OpenMP, flags
 # to compile with OpenMP are returned and added and the HAS_OPENMP definition
-# is added for compilation.  Visual Studio doesn't support unsigned type index
-# in for loop, while OpenMP uses uword in for loop.  So to prevent build crash
-# on Windows, disable OpenMP, even though its support was detected.
+# is added for compilation.
 #
 # This way we can skip calls to functions defined in omp.h with code like:
 # #ifdef HAS_OPENMP
@@ -423,7 +421,7 @@ add_definitions(-DBOOST_TEST_DYN_LINK)
 #   ... openMP code here ...
 # }
 # #endif
-if (NOT WIN32 AND USE_OPENMP)
+if (USE_OPENMP)
   find_package(OpenMP)
 endif ()
 

src/mlpack/core/dists/diagonal_gaussian_distribution.cpp

@@ -116,8 +116,7 @@ void DiagonalGaussianDistribution::Train(const arma::mat& observations,
   mean = observations * normalizedProbs;
 
   // Now calculate the covariance.
-  const arma::mat diffs = observations - mean *
-      arma::ones<arma::rowvec>(observations.n_cols);
+  const arma::mat diffs = observations.each_col() - mean;
   covariance += (diffs % diffs) * normalizedProbs;
 
   // Calculate the sum of each weight squared.

src/mlpack/core/dists/diagonal_gaussian_distribution.hpp

@@ -165,8 +165,7 @@ inline void DiagonalGaussianDistribution::LogProbability(
 
   // Column i of 'diffs' is the difference between observations.col(i) and
   // the mean.
-  arma::mat diffs = observations - mean *
-      arma::ones<arma::rowvec>(observations.n_cols);
+  arma::mat diffs = observations.each_col() - mean;
 
   // Calculates log of exponent equation in multivariate gaussian distribution.
   // We use only diagonal part for faster computation.

src/mlpack/methods/gmm/em_fit_impl.hpp

@@ -46,8 +46,14 @@ Estimate(const arma::mat& observations,
 {
   if (std::is_same<CovarianceConstraintPolicy, DiagonalConstraint>::value)
   {
-    ArmadilloGMMWrapper(observations, dists, weights, useInitialModel);
-    return;
+    #ifdef _WIN32
+      Log::Warn << "Cannot use arma::gmm_diag on Visual Studio due to OpenMP"
+          << " compilation issues! Using slower EMFit::Estimate() instead..."
+          << std::endl;
+    #else
+      ArmadilloGMMWrapper(observations, dists, weights, useInitialModel);
+      return;
+    #endif
   }
 
   // Only perform initial clustering if the user wanted it.
@@ -100,18 +106,29 @@ Estimate(const arma::mat& observations,
       // Don't update if there's no probability of the Gaussian having points.
       if (probRowSums[i] != 0)
         dists[i].Mean() = (observations * condProb.col(i)) / probRowSums[i];
+      else
+        continue;
 
       // Calculate the new value of the covariances using the updated
       // conditional probabilities and the updated means.
-      arma::mat tmp = observations - (dists[i].Mean() *
-          arma::ones<arma::rowvec>(observations.n_cols));
-      arma::mat tmpB = tmp % (arma::ones<arma::vec>(observations.n_rows) *
-          trans(condProb.col(i)));
+      arma::mat tmp = observations.each_col() - dists[i].Mean();
 
-      // Don't update if there's no probability of the Gaussian having points.
-      if (probRowSums[i] != 0.0)
+      // If the distribution is DiagonalGaussianDistribution, calculate the
+      // covariance only with diagonal components.
+      if (std::is_same<Distribution,
+          distribution::DiagonalGaussianDistribution>::value)
       {
+        arma::vec covariance = arma::sum((tmp % tmp) %
+            (arma::ones<arma::vec>(observations.n_rows) *
+            trans(condProb.col(i))), 1) / probRowSums[i];
+        covariance = arma::clamp(covariance, 1e-10, DBL_MAX);
+        dists[i].Covariance(std::move(covariance));
+      }
+      else
+      {
+        arma::mat tmpB = tmp.each_row() % trans(condProb.col(i));
         arma::mat covariance = (tmp * trans(tmpB)) / probRowSums[i];
+
         // Apply covariance constraint.
         constraint.ApplyConstraint(covariance);
         dists[i].Covariance(std::move(covariance));
@@ -190,14 +207,21 @@ Estimate(const arma::mat& observations,
       // model.
       probRowSums[i] = accu(condProb.col(i) % probabilities);
 
-      dists[i].Mean() = (observations * (condProb.col(i) % probabilities)) /
-        probRowSums[i];
+      // Don't update if there's no probability of the Gaussian having points.
+      if (probRowSums[i] != 0)
+      {
+        dists[i].Mean() = (observations * (condProb.col(i) % probabilities)) /
+            probRowSums[i];
+      }
+      else
+        continue;
 
       // Calculate the new value of the covariances using the updated
       // conditional probabilities and the updated means.
-      arma::mat tmp = observations - (dists[i].Mean() *
-          arma::ones<arma::rowvec>(observations.n_cols));
+      arma::mat tmp = observations.each_col() - dists[i].Mean();
 
+      // If the distribution is DiagonalGaussianDistribution, calculate the
+      // covariance only with diagonal components.
       if (std::is_same<Distribution,
           distribution::DiagonalGaussianDistribution>::value)
       {
@@ -209,9 +233,8 @@ Estimate(const arma::mat& observations,
       }
       else
       {
-        arma::mat tmpB = tmp % (arma::ones<arma::vec>(observations.n_rows) *
-            trans(condProb.col(i) % probabilities));
-
+        arma::mat tmpB = tmp.each_row() % trans(condProb.col(i) %
+            probabilities);
         arma::mat cov = (tmp * trans(tmpB)) / probRowSums[i];
 
         // Apply covariance constraint.