From 4093add9eb873cc17125f16b5387f018b4c6a4d6 Mon Sep 17 00:00:00 2001
From: Emmanuel Benazera <emmanuel.benazera@lri.fr>
Date: Fri, 27 Feb 2015 12:11:56 +0100
Subject: [PATCH] added correlation matrix computation, including for vd-cma,
 ref #129

---
 cmasolutions.cc | 27 +++++++++++++++++++++++++++
 cmasolutions.h  |  6 ++++++
 2 files changed, 33 insertions(+)

diff --git a/cmasolutions.cc b/cmasolutions.cc
index 8dec699c..1fd56fa1 100644
--- a/cmasolutions.cc
+++ b/cmasolutions.cc
@@ -125,6 +125,33 @@ namespace libcmaes
       }
   }
 
+  dMat CMASolutions::corr() const
+  {
+    dMat corr, dinvcov;
+    if (_cov.size() > 0) // full cov
+      {
+	dinvcov = _cov.diagonal().cwiseSqrt().cwiseInverse();
+	corr = dMat(_cov.rows(),_cov.cols());
+	for (int i=0;i<_cov.cols();i++)
+	  corr.col(i) = _cov.col(i).cwiseProduct(dinvcov);
+	for (int i=0;i<_cov.rows();i++)
+	  corr.row(i) = corr.row(i).cwiseProduct(dinvcov.transpose());
+      }
+    else if (_v.size() > 0) //vd
+      {
+	// we need to compute the full covariance matrix, which is counter productive in large-scale settings
+	dMat cov =_sepcov.asDiagonal()*(dMat::Identity(_sepcov.rows(),_sepcov.rows())+_v*_v.transpose())*(_sepcov.asDiagonal());
+	dinvcov = cov.diagonal().cwiseSqrt().cwiseInverse();
+	corr = dMat(cov.rows(),cov.cols());
+	for (int i=0;i<cov.cols();i++)
+	  corr.col(i) = cov.col(i).cwiseProduct(dinvcov);
+	for (int i=0;i<cov.rows();i++)
+	  corr.row(i) = corr.row(i).cwiseProduct(dinvcov.transpose());
+      }
+    else return dMat::Constant(_sepcov.rows(),1,1.0); // sep
+    return corr;
+  }
+
   void CMASolutions::update_eigenv(const dVec &eigenvalues,
 				   const dMat &eigenvectors)
   {
diff --git a/cmasolutions.h b/cmasolutions.h
index 8a00ff44..c349f259 100644
--- a/cmasolutions.h
+++ b/cmasolutions.h
@@ -293,6 +293,12 @@ namespace libcmaes
       return dVec(); // should never reach here.
     }
     
+    /**
+     * \brief returns correlation matrix
+     * @return correlation matrix
+     */
+    dMat corr() const;
+
     /**
      * \brief returns current value of step-size sigma
      * @return current step-size