changeset 92:81eaba98985b bqvec

Complete using bqvec for allocation etc, but with unchanged processing logic

--- a/src/EM.cpp	Tue May 06 13:49:52 2014 +0100
+++ b/src/EM.cpp	Tue May 06 14:19:19 2014 +0100
@@ -68,6 +68,11 @@
 
 EM::~EM()
 {
+    deallocate(m_q);
+    deallocate(m_estimate);
+    deallocate_channels(m_sources, m_sourceCount);
+    deallocate_channels(m_shifts, m_shiftCount);
+    deallocate(m_pitches);
 }
 
 void
@@ -86,41 +91,38 @@
 }
 
 void
-EM::normaliseColumn(V &column)
+EM::normaliseColumn(double *column, int size)
 {
-    double sum = 0.0;
-    for (int i = 0; i < (int)column.size(); ++i) {
-        sum += column[i];
-    }
-    for (int i = 0; i < (int)column.size(); ++i) {
-        column[i] /= sum;
-    }
+    double sum = v_sum(column, size);
+    v_scale(column, 1.0 / sum, size);
 }
 
 void
-EM::normaliseGrid(Grid &grid)
+EM::normaliseGrid(double **grid, int size1, int size2)
 {
-    V denominators(grid[0].size());
+    double *denominators = allocate_and_zero<double>(size2);
 
-    for (int i = 0; i < (int)grid.size(); ++i) {
-        for (int j = 0; j < (int)grid[i].size(); ++j) {
+    for (int i = 0; i < size1; ++i) {
+        for (int j = 0; j < size2; ++j) {
             denominators[j] += grid[i][j];
         }
     }
 
-    for (int i = 0; i < (int)grid.size(); ++i) {
-        for (int j = 0; j < (int)grid[i].size(); ++j) {
-            grid[i][j] /= denominators[j];
-        }
+    for (int i = 0; i < size1; ++i) {
+        v_divide(grid[i], denominators, size2);
     }
+
+    deallocate(denominators);
 }
 
 void
-EM::iterate(V column)
+EM::iterate(const double *column)
 {
-    normaliseColumn(column);
-    expectation(column);
-    maximisation(column);
+    double *norm = allocate<double>(m_binCount);
+    v_copy(norm, column, m_binCount);
+    normaliseColumn(norm, m_binCount);
+    expectation(norm);
+    maximisation(norm);
 }
 
 const double *
@@ -130,7 +132,7 @@
 }
 
 void
-EM::expectation(const V &column)
+EM::expectation(const double *column)
 {
 //    cerr << ".";
 
@@ -159,11 +161,20 @@
 }
 
 void
-EM::maximisation(const V &column)
+EM::maximisation(const double *column)
 {
-    V newPitches(m_noteCount, epsilon);
-    Grid newShifts(m_shiftCount, V(m_noteCount, epsilon));
-    Grid newSources(m_sourceCount, V(m_noteCount, epsilon));
+    double *newPitches = allocate<double>(m_noteCount);
+    v_set(newPitches, epsilon, m_noteCount);
+
+    double **newShifts = allocate_channels<double>(m_shiftCount, m_noteCount);
+    for (int i = 0; i < m_shiftCount; ++i) {
+        v_set(newShifts[i], epsilon, m_noteCount);
+    }
+
+    double **newSources = allocate_channels<double>(m_sourceCount, m_noteCount);
+    for (int i = 0; i < m_sourceCount; ++i) {
+        v_set(newSources[i], epsilon, m_noteCount);
+    }
 
     for (int n = 0; n < m_noteCount; ++n) {
 
@@ -210,9 +221,13 @@
         }
     }
 
-    normaliseColumn(newPitches);
-    normaliseGrid(newShifts);
-    normaliseGrid(newSources);
+    normaliseColumn(newPitches, m_noteCount);
+    normaliseGrid(newShifts, m_shiftCount, m_noteCount);
+    normaliseGrid(newSources, m_sourceCount, m_noteCount);
+
+    deallocate(m_pitches);
+    deallocate_channels(m_shifts, m_shiftCount);
+    deallocate_channels(m_sources, m_sourceCount);
 
     m_pitches = newPitches;
     m_shifts = newShifts;

--- a/src/EM.h	Tue May 06 13:49:52 2014 +0100
+++ b/src/EM.h	Tue May 06 14:19:19 2014 +0100
@@ -28,7 +28,7 @@
     int getNoteCount() const { return m_noteCount; } // size of pitch column
     int getSourceCount() const { return m_sourceCount; }
 
-    void iterate(double *column);
+    void iterate(const double *column);
 
     const double *getEstimate() const { // bin count
 	return m_estimate;
@@ -36,7 +36,7 @@
     const double *getPitchDistribution() const { // note count
 	return m_pitches;
     }
-    const double **getSources() const { // source count * note count
+    const double *const *getSources() const { // source count * note count
 	return m_sources; 
     }
 
@@ -60,10 +60,10 @@
     const int m_highestPitch;
 
     void normaliseColumn(double *column, int size);
-    void normaliseGrid(double **grid, int width, int height);
+    void normaliseGrid(double **grid, int size1, int size2);
 
-    void expectation(double *column); // size is m_binCount
-    void maximisation(double *column); // size is m_binCount
+    void expectation(const double *column); // size is m_binCount
+    void maximisation(const double *column); // size is m_binCount
 
     const double *templateFor(int instrument, int note, int shift);
     void rangeFor(int instrument, int &minPitch, int &maxPitch);

--- a/src/Silvet.cpp	Tue May 06 13:49:52 2014 +0100
+++ b/src/Silvet.cpp	Tue May 06 14:19:19 2014 +0100
@@ -398,10 +398,11 @@
 
         EM em;
         for (int j = 0; j < iterations; ++j) {
-            em.iterate(filtered[i]);
+            em.iterate(filtered[i].data());
         }
 
-        vector<double> pitches = em.getPitchDistribution();
+        const double *pd = em.getPitchDistribution();
+        vector<double> pitches(pd, pd + processingNotes);
         
         for (int j = 0; j < processingNotes; ++j) {
             pitches[j] *= sum;