# HG changeset patch
# User tomwalters
# Date 1266326603 0
# Node ID e91769e84be111e5916012da685173cb050f0e85
# Parent  bc394a9850420339c41fa24ed9a43c4d1e1adeeb
- Fixed the gaussian fitting to use doubles internally (like MATLAB)
- Fixed an irritating bug that was causing the Gaussian fitting to be incorrectly initialized, leading to small differences between the AIM-C and MATLAB code
- Finalised the Gaussian fitting test

diff -r bc394a985042 -r e91769e84be1 src/Modules/Features/ModuleGaussians.cc
--- a/src/Modules/Features/ModuleGaussians.cc	Mon Feb 15 20:37:26 2010 +0000
+++ b/src/Modules/Features/ModuleGaussians.cc	Tue Feb 16 13:23:23 2010 +0000
@@ -97,10 +97,10 @@
          ++iSample) {
       m_pSpectralProfile[iChannel] += input[iChannel][iSample];
     }
-    m_pSpectralProfile[iChannel] /= static_cast<float>(input.buffer_length());
+    m_pSpectralProfile[iChannel] /= static_cast<double>(input.buffer_length());
   }
 
-  float spectral_profile_sum = 0.0f;
+  double spectral_profile_sum = 0.0f;
   for (int i = 0; i < input.channel_count(); ++i) {
     spectral_profile_sum += m_pSpectralProfile[i];
   }
@@ -120,22 +120,33 @@
 
   RubberGMMCore(2, true);
 
-  float fMean1 = m_pMu[0];
-  float fMean2 = m_pMu[1];
+  double fMean1 = m_pMu[0];
+  double fMean2 = m_pMu[1];
+  //LOG_INFO(_T("Orig. mean 0 = %f"), m_pMu[0]);
+  //LOG_INFO(_T("Orig. mean 1 = %f"), m_pMu[1]);
+  //LOG_INFO(_T("Orig. prob 0 = %f"), m_pA[0]);
+  //LOG_INFO(_T("Orig. prob 1 = %f"), m_pA[1]);
 
-  float fA1 = 0.05 * m_pA[0];
-  float fA2 = 1.0 - 0.25 * m_pA[1];
+  double fA1 = 0.05 * m_pA[0];
+  double fA2 = 1.0 - 0.25 * m_pA[1];
 
-  float fGradient = (fMean2 - fMean1) / (fA2 - fA1);
-  float fIntercept = fMean2 - fGradient * fA2;
+  //LOG_INFO(_T("fA1 = %f"), fA1);
+  //LOG_INFO(_T("fA2 = %f"), fA2);
+
+  double fGradient = (fMean2 - fMean1) / (fA2 - fA1);
+  double fIntercept = fMean2 - fGradient * fA2;
+
+  //LOG_INFO(_T("fGradient = %f"), fGradient);
+  //LOG_INFO(_T("fIntercept = %f"), fIntercept);
 
   for (int i = 0; i < m_iParamNComp; ++i) {
-    m_pMu[i] = ((float)i / (float)m_iParamNComp - 1.0f)
-                  * -fGradient + fIntercept;
+    m_pMu[i] = ((double)i / ((double)m_iParamNComp - 1.0f))
+                  * fGradient + fIntercept;
+                //LOG_INFO(_T("mean %d = %f"), i, m_pMu[i]);
   }
 
   for (int i = 0; i < m_iParamNComp; ++i) {
-    m_pA[i] = 1.0f / (float)m_iParamNComp;
+    m_pA[i] = 1.0f / (double)m_iParamNComp;
   }
 
   RubberGMMCore(m_iParamNComp, false);
@@ -155,7 +166,7 @@
   int iSizeX = m_iNumChannels;
 
   // Normalise the spectral profile
-  float fSpectralProfileTotal = 0.0f;
+  double fSpectralProfileTotal = 0.0f;
   for (int iCount = 0; iCount < iSizeX; iCount++) {
     fSpectralProfileTotal += m_pSpectralProfile[iCount];
   }
@@ -165,18 +176,18 @@
 
   if (bDoInit) {
     // Uniformly spaced components
-    float dd = (iSizeX - 1.0f) / iNComponents;
+    double dd = (iSizeX - 1.0f) / iNComponents;
     for (int i = 0; i < iNComponents; i++) {
       m_pMu[i] = dd / 2.0f + (i * dd);
       m_pA[i] = 1.0f / iNComponents;
     }
   }
 
-  vector<float> pA_old;
+  vector<double> pA_old;
   pA_old.resize(iNComponents);
-  vector<float> pP_mod_X;
+  vector<double> pP_mod_X;
   pP_mod_X.resize(iSizeX);
-  vector<float> pP_comp;
+  vector<double> pP_comp;
   pP_comp.resize(iSizeX * iNComponents);
 
   for (int iIteration = 0; iIteration < m_iParamMaxIt; iIteration++) {
@@ -189,7 +200,7 @@
     for (int i = 0; i < iNComponents; i++) {
       for (int iCount = 0; iCount < iSizeX; iCount++) {
         pP_mod_X[iCount] += 1.0f / sqrt(2.0f * M_PI * m_fParamVar)
-                            * exp((-0.5f) * pow(((float)iCount-m_pMu[i]), 2)
+                            * exp((-0.5f) * pow(((double)(iCount + 1)-m_pMu[i]), 2)
                                   / m_fParamVar) * m_pA[i];
       }
     }
@@ -202,14 +213,14 @@
       for (int iCount = 0; iCount < iSizeX; iCount++) {
         pP_comp[iCount + i * iSizeX] =
           1.0f / sqrt(2.0f * M_PI * m_fParamVar)
-          * exp((-0.5f) * pow(((float)iCount - m_pMu[i]), 2) / m_fParamVar);
+          * exp((-0.5f) * pow(((double)(iCount + 1) - m_pMu[i]), 2) / m_fParamVar);
         pP_comp[iCount + i * iSizeX] =
           pP_comp[iCount + i * iSizeX] * m_pA[i] / pP_mod_X[iCount];
       }
     }
 
     for (int iCount = 0; iCount < iSizeX; ++iCount) {
-      float fSum = 0.0f;
+      double fSum = 0.0f;
       for (int i = 0; i < iNComponents; ++i) {
         pP_comp[iCount+i*iSizeX] = pow(pP_comp[iCount + i * iSizeX],
                                        m_fParamPosteriorExp); // expansion
@@ -229,25 +240,27 @@
     }
 
     // finish when already converged
-    float fPrdist = 0.0f;
+    double fPrdist = 0.0f;
     for (int i = 0; i < iNComponents; ++i) {
       fPrdist += pow((m_pA[i] - pA_old[i]), 2);
     }
     fPrdist /= iNComponents;
 
     if (fPrdist < m_fParamPriorsConverged) {
-      LOG_INFO("Converged!");
+      //LOG_INFO("Converged!");
       break;
     }
+    //LOG_INFO("Didn't converge!");
+
 
     // update means (positions)
     for (int i = 0 ; i < iNComponents; ++i) {
-      float mu_old = m_pMu[i];
+      double mu_old = m_pMu[i];
       if (m_pA[i] > 0.0f) {
         m_pMu[i] = 0.0f;
         for (int iCount = 0; iCount < iSizeX; ++iCount) {
           m_pMu[i] += m_pSpectralProfile[iCount]
-                      * pP_comp[iCount + i * iSizeX] * (float)iCount;
+                      * pP_comp[iCount + i * iSizeX] * (double)(iCount + 1);
         }
         m_pMu[i] /= m_pA[i];
         if (isnan(m_pMu[i])) {
@@ -263,7 +276,7 @@
     bSorted = true;
     for (int i = 0; i < iNComponents - 1; ++i) {
       if (m_pMu[i] > m_pMu[i + 1]) {
-        float fTemp = m_pMu[i];
+        double fTemp = m_pMu[i];
         m_pMu[i] = m_pMu[i + 1];
         m_pMu[i + 1] = fTemp;
         fTemp = m_pA[i];
diff -r bc394a985042 -r e91769e84be1 src/Modules/Features/ModuleGaussians.h
--- a/src/Modules/Features/ModuleGaussians.h	Mon Feb 15 20:37:26 2010 +0000
+++ b/src/Modules/Features/ModuleGaussians.h	Tue Feb 16 13:23:23 2010 +0000
@@ -66,11 +66,11 @@
 
   /*! \brief Constant variance of Gaussians
    */
-  float m_fParamVar;
+  double m_fParamVar;
 
   /*! \brief posterior probability expansion exponent
    */
-  float m_fParamPosteriorExp;
+  double m_fParamPosteriorExp;
 
   /*! \brief Maximum Number of iterations
    */
@@ -78,19 +78,19 @@
 
   /*! \brief convergence criterion
    */
-  float m_fParamPriorsConverged;
+  double m_fParamPriorsConverged;
 
   /*! \brief The amplitudes of the components (priors)
    */
-  vector<float> m_pA;
+  vector<double> m_pA;
 
   /*! \brief The means of the components (priors)
    */
-  vector<float> m_pMu;
+  vector<double> m_pMu;
 
   /*! \brief The spectral profile of the incoming buffer
    */
-  vector<float> m_pSpectralProfile;
+  vector<double> m_pSpectralProfile;
 
   int m_iNumChannels;
 };
diff -r bc394a985042 -r e91769e84be1 src/Modules/Features/ModuleGaussians_test.py
--- a/src/Modules/Features/ModuleGaussians_test.py	Mon Feb 15 20:37:26 2010 +0000
+++ b/src/Modules/Features/ModuleGaussians_test.py	Tue Feb 16 13:23:23 2010 +0000
@@ -29,17 +29,55 @@
 import aimc
 import matplotlib
 import pylab
-import scipy
+from scipy import io
 
 def main():
   data_file = "src/Modules/Features/testdata/aa153.0p108.1s100.0t+000itd.mat"
-  data = scipy.io.loadmat(data_file)
+  data = io.loadmat(data_file)
+  
+  # The margin of error allowed between the returned values from AIM-C and
+  # the stored MATLAB values.
+  epsilon = 0.000001;
   
   given_profiles = data["Templates"]
   matlab_features = data["feature"]
   
+  (profile_count, channel_count) = given_profiles.shape
   
+  profile_sig = aimc.SignalBank()
+  profile_sig.Initialize(channel_count, 1, 44100)
+  parameters = aimc.Parameters()
+  mod_gauss = aimc.ModuleGaussians(parameters)
+  mod_gauss.Initialize(profile_sig)
   
+  correct_count = 0;
+  incorrect_count  = 0;
+  for p in range(0, profile_count):
+    profile = given_profiles[p]
+    features = matlab_features[p]
+    for i in range(0, channel_count):
+      profile_sig.set_sample(i, 0, profile[i])
+    mod_gauss.Process(profile_sig)
+    out_sig = mod_gauss.GetOutputBank()
+    error = False;
+    for j in range(0, out_sig.channel_count()):
+      if (abs(out_sig.sample(j, 0) - features[j]) > epsilon):
+        error = True;
+        incorrect_count += 1;
+      else:
+        correct_count += 1;
+    if error:
+      print("Mismatch at profile %d" % (p))
+      print("AIM-C values: %f %f %f %f" % (out_sig.sample(0, 0), out_sig.sample(1, 0), out_sig.sample(2, 0), out_sig.sample(3, 0)))
+      print("MATLAB values: %f %f %f %f" % (features[0], features[1], features[2], features[3]))
+      print("")
+    percent_correct = 100 * correct_count / (correct_count + incorrect_count)
+  print("Total correct: %f percent" % (percent_correct))
+  if percent_correct == 100:
+    print("=== TEST PASSED ===")
+  else:
+    print("=== TEST FAILED! ===")
+
   pass