comparison src/EM.cpp @ 114:b2f0967cb8d1 bqvec-openmp

Pull out shifted templates into separate arrays with proper alignment (trading off memory occupied against getting the alignments right)

 
 using namespace breakfastquay;
 
 static double epsilon = 1e-16;
 
+bool EM::m_initialised = false;
+double ****EM::m_templates = 0;
+
 EM::EM(bool useShifts) :
     m_noteCount(SILVET_TEMPLATE_NOTE_COUNT),
     m_shiftCount(useShifts ? SILVET_TEMPLATE_MAX_SHIFT * 2 + 1 : 1),
     m_binCount(SILVET_TEMPLATE_HEIGHT),
     m_sourceCount(SILVET_TEMPLATE_COUNT),
     m_pitchSparsity(1.1),
     m_sourceSparsity(1.3),
     m_lowestPitch(silvet_templates_lowest_note),
     m_highestPitch(silvet_templates_highest_note)
 {
+    if (!m_initialised) {
+        cerr << "ERROR: You must call EM::initialise() before constructing any EM objects" << endl;
+        abort();
+    }
+
     m_pitches = allocate<double>(m_noteCount);
     m_updatePitches = allocate<double>(m_noteCount);
     for (int n = 0; n < m_noteCount; ++n) {
         m_pitches[n] = drand48();
+    }
+    
+    m_sources = allocate_channels<double>(m_sourceCount, m_noteCount);
+    m_updateSources = allocate_channels<double>(m_sourceCount, m_noteCount);
+    for (int i = 0; i < m_sourceCount; ++i) {
+        for (int n = 0; n < m_noteCount; ++n) {
+            m_sources[i][n] = (inRange(i, n) ? 1.0 : 0.0);
+        }
     }
 
     if (useShifts) {
         m_shifts = allocate_channels<double>(m_shiftCount, m_noteCount);
         m_updateShifts = allocate_channels<double>(m_shiftCount, m_noteCount);

             }
         }
     } else {
         m_shifts = 0;
         m_updateShifts = 0;
-    }
-    
-    m_sources = allocate_channels<double>(m_sourceCount, m_noteCount);
-    m_updateSources = allocate_channels<double>(m_sourceCount, m_noteCount);
-    for (int i = 0; i < m_sourceCount; ++i) {
-        for (int n = 0; n < m_noteCount; ++n) {
-            m_sources[i][n] = (inRange(i, n) ? 1.0 : 0.0);
-        }
     }
 
     m_estimate = allocate<double>(m_binCount);
     m_q = allocate<double>(m_binCount);
 }

     deallocate(m_pitches);
     deallocate(m_updatePitches);
 }
 
 void
+EM::initialise()
+{
+    //!!! need mutex
+
+    if (m_initialised) return;
+    m_templates = new double ***[SILVET_TEMPLATE_COUNT];
+    for (int i = 0; i < SILVET_TEMPLATE_COUNT; ++i) {
+        m_templates[i] = new double **[SILVET_TEMPLATE_NOTE_COUNT];
+        for (int n = 0; n < SILVET_TEMPLATE_NOTE_COUNT; ++n) {
+            m_templates[i][n] = new double *[SILVET_TEMPLATE_MAX_SHIFT * 2 + 1];
+            for (int f = 0; f < SILVET_TEMPLATE_MAX_SHIFT * 2 + 1; ++f) {
+                m_templates[i][n][f] = allocate<double>(SILVET_TEMPLATE_HEIGHT);
+                const float *t = silvet_templates[i].data[n] + f;
+                v_convert(m_templates[i][n][f], t, SILVET_TEMPLATE_HEIGHT);
+            }
+        }
+    }
+    m_initialised = true;
+}
+
+void
 EM::rangeFor(int instrument, int &minPitch, int &maxPitch)
 {
     minPitch = silvet_templates[instrument].lowest;
     maxPitch = silvet_templates[instrument].highest;
 }

 
 const double *
 EM::templateFor(int instrument, int note, int shift)
 {
     if (m_shifts) {
-        return silvet_templates[instrument].data[note] + shift;
+        return m_templates[instrument][note][shift];
     } else {
-        return silvet_templates[instrument].data[note] + 
-            SILVET_TEMPLATE_MAX_SHIFT;
+        return m_templates[instrument][note][SILVET_TEMPLATE_MAX_SHIFT];
     }
 }
 
 void
 EM::expectation(const double *column)