From 4d98ca69be75c58733b1c7be6a30ec0892e8cdc3 Mon Sep 17 00:00:00 2001
From: Emmanuel Benazera <emmanuel.benazera@lri.fr>
Date: Thu, 22 Jan 2015 14:20:58 +0100
Subject: [PATCH] improved normalization scheme for TPA samples, ref #88

---
 cmastrategy.cc | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/cmastrategy.cc b/cmastrategy.cc
index 678579d1..3afbdf76 100644
--- a/cmastrategy.cc
+++ b/cmastrategy.cc
@@ -220,10 +220,17 @@ namespace libcmaes
     if (eostrat<TGenoPheno>::_parameters._tpa && eostrat<TGenoPheno>::_niter > 0)
       {
 	dVec mean_shift = eostrat<TGenoPheno>::_solutions._xmean - eostrat<TGenoPheno>::_solutions._xmean_prev;
-	double mean_shift_norm = mean_shift.norm();
-	dVec z = eostrat<TGenoPheno>::_parameters._chi * (mean_shift / mean_shift_norm);
-	eostrat<TGenoPheno>::_solutions._tpa_x1 = eostrat<TGenoPheno>::_solutions._xmean + eostrat<TGenoPheno>::_solutions._sigma * z;
-	eostrat<TGenoPheno>::_solutions._tpa_x2 = eostrat<TGenoPheno>::_solutions._xmean - eostrat<TGenoPheno>::_solutions._sigma * z;
+	double mean_shift_norm = 1.0;
+	if (!eostrat<TGenoPheno>::_parameters._sep && !eostrat<TGenoPheno>::_parameters._vd)
+	  mean_shift_norm = (_esolver._eigenSolver.eigenvalues().cwiseInverse().cwiseProduct(_esolver._eigenSolver.eigenvectors().transpose()*mean_shift)).norm() / eostrat<TGenoPheno>::_solutions._sigma;
+	else mean_shift_norm = eostrat<TGenoPheno>::_solutions._sepcov.cwiseInverse().cwiseProduct(mean_shift).norm() / eostrat<TGenoPheno>::_solutions._sigma;
+	//std::cout << "mean_shift_norm=" << mean_shift_norm << " / sqrt(N)=" << std::sqrt(std::sqrt(eostrat<TGenoPheno>::_parameters._dim)) << std::endl;
+
+	dMat rz = _esolver.samples_ind(1);
+	double mfactor = rz.norm();
+	dVec z = mfactor * (mean_shift / mean_shift_norm);
+	eostrat<TGenoPheno>::_solutions._tpa_x1 = eostrat<TGenoPheno>::_solutions._xmean + z;
+	eostrat<TGenoPheno>::_solutions._tpa_x2 = eostrat<TGenoPheno>::_solutions._xmean - z;
 	
 	// if gradient is in col 0, move tpa vectors to pos 1 & 2
 	int p1 = 0, p2 = 1;