Mercurial > hg > silvet
changeset 112:2169e7a448c5 bqvec-openmp
Merge from default branch
| author | Chris Cannam |
|---|---|
| date | Wed, 07 May 2014 09:01:32 +0100 |
| parents | 9b299b087dd4 (diff) e282930cfca7 (current diff) |
| children | c4eae816bdb3 |
| files | src/EM.cpp src/EM.h src/Silvet.cpp testdata/timing/results.txt |
| diffstat | 8 files changed, 197 insertions(+), 162 deletions(-) [+] |
line wrap: on
line diff
--- a/Makefile.inc Tue May 06 18:55:11 2014 +0100 +++ b/Makefile.inc Wed May 07 09:01:32 2014 +0100 @@ -3,6 +3,7 @@ QMDSP_DIR ?= ../qm-dsp VAMPSDK_DIR ?= ../vamp-plugin-sdk +BQVEC_DIR ?= ../bqvec CQ_DIR ?= constant-q-cpp/cpp-qm-dsp @@ -12,7 +13,7 @@ CC ?= gcc CFLAGS := $(CFLAGS) -CXXFLAGS := -I. -I$(VAMPSDK_DIR) -I$(QMDSP_DIR) $(CXXFLAGS) +CXXFLAGS := -I. -I$(VAMPSDK_DIR) -I$(QMDSP_DIR) -I$(BQVEC_DIR) $(CXXFLAGS) LDFLAGS := $(LDFLAGS) PLUGIN_LDFLAGS := $(LDFLAGS) $(PLUGIN_LDFLAGS)
--- a/Makefile.linux Tue May 06 18:55:11 2014 +0100 +++ b/Makefile.linux Wed May 07 09:01:32 2014 +0100 @@ -1,11 +1,12 @@ -CFLAGS := -Wall -O3 -ffast-math -msse -mfpmath=sse -ftree-vectorize -fPIC -I../vamp-plugin-sdk/ +CFLAGS := -Wall -O3 -fopenmp -ffast-math -msse -mfpmath=sse -ftree-vectorize -fPIC -I../vamp-plugin-sdk/ + #CFLAGS := -g -fPIC -I../vamp-plugin-sdk CXXFLAGS := $(CFLAGS) VAMPSDK_DIR := ../vamp-plugin-sdk -PLUGIN_LDFLAGS := -shared -Wl,--version-script=vamp-plugin.map +PLUGIN_LDFLAGS := -lgomp -shared -Wl,--version-script=vamp-plugin.map PLUGIN_EXT := .so
--- a/data/include/templates.h Tue May 06 18:55:11 2014 +0100 +++ b/data/include/templates.h Wed May 07 09:01:32 2014 +0100 @@ -15,7 +15,7 @@ const char *name; int lowest; int highest; - float data[SILVET_TEMPLATE_NOTE_COUNT][SILVET_TEMPLATE_SIZE]; + double data[SILVET_TEMPLATE_NOTE_COUNT][SILVET_TEMPLATE_SIZE]; } silvet_template_t; static int silvet_templates_lowest_note = 15;
--- a/src/EM.cpp Tue May 06 18:55:11 2014 +0100 +++ b/src/EM.cpp Wed May 07 09:01:32 2014 +0100 @@ -22,12 +22,15 @@ #include <iostream> -#include <vector> +#include "VectorOps.h" +#include "Allocators.h" using std::vector; using std::cerr; using std::endl; +using namespace breakfastquay; + static double epsilon = 1e-16; EM::EM(bool useShifts) : @@ -35,21 +38,21 @@ m_noteCount(SILVET_TEMPLATE_NOTE_COUNT), m_shiftCount(useShifts ? SILVET_TEMPLATE_MAX_SHIFT * 2 + 1 : 1), m_binCount(SILVET_TEMPLATE_HEIGHT), - m_instrumentCount(SILVET_TEMPLATE_COUNT), + m_sourceCount(SILVET_TEMPLATE_COUNT), m_pitchSparsity(1.1), - m_sourceSparsity(1.3) + m_sourceSparsity(1.3), + m_lowestPitch(silvet_templates_lowest_note), + m_highestPitch(silvet_templates_highest_note) { - m_lowestPitch = silvet_templates_lowest_note; - m_highestPitch = silvet_templates_highest_note; - - m_pitches = V(m_noteCount); + 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_shifts = Grid(m_shiftCount); + m_shifts = allocate_channels<double>(m_shiftCount, m_noteCount); + m_updateShifts = allocate_channels<double>(m_shiftCount, m_noteCount); for (int f = 0; f < m_shiftCount; ++f) { - m_shifts[f] = V(m_noteCount); for (int n = 0; n < m_noteCount; ++n) { if (m_useShifts) { m_shifts[f][n] = drand48(); @@ -59,20 +62,28 @@ } } - m_sources = Grid(m_instrumentCount); - for (int i = 0; i < m_instrumentCount; ++i) { - m_sources[i] = V(m_noteCount); + 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 = V(m_binCount); - m_q = V(m_binCount); + m_estimate = allocate<double>(m_binCount); + m_q = allocate<double>(m_binCount); } EM::~EM() { + deallocate(m_q); + deallocate(m_estimate); + deallocate_channels(m_sources, m_sourceCount); + deallocate_channels(m_updateSources, m_sourceCount); + deallocate_channels(m_shifts, m_shiftCount); + deallocate_channels(m_updateShifts, m_shiftCount); + deallocate(m_pitches); + deallocate(m_updatePitches); } void @@ -91,44 +102,42 @@ } 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); + deallocate(norm); } -const float * +const double * EM::templateFor(int instrument, int note, int shift) { if (m_useShifts) { @@ -140,24 +149,21 @@ } void -EM::expectation(const V &column) +EM::expectation(const double *column) { // cerr << "."; - for (int i = 0; i < m_binCount; ++i) { - m_estimate[i] = epsilon; - } + v_set(m_estimate, epsilon, m_binCount); - for (int i = 0; i < m_instrumentCount; ++i) { + for (int i = 0; i < m_sourceCount; ++i) { for (int n = 0; n < m_noteCount; ++n) { + const double pitch = m_pitches[n]; + const double source = m_sources[i][n]; for (int f = 0; f < m_shiftCount; ++f) { - const float *w = templateFor(i, n, f); - double pitch = m_pitches[n]; - double source = m_sources[i][n]; - double shift = m_shifts[f][n]; - for (int j = 0; j < m_binCount; ++j) { - m_estimate[j] += w[j] * pitch * source * shift; - } + const double *w = templateFor(i, n, f); + const double shift = m_shifts[f][n]; + const double factor = pitch * source * shift; + v_add_with_gain(m_estimate, w, factor, m_binCount); } } } @@ -168,77 +174,78 @@ } void -EM::maximisation(const V &column) +EM::maximisation(const double *column) { - V newPitches = m_pitches; + v_set(m_updatePitches, epsilon, m_noteCount); + for (int i = 0; i < m_shiftCount; ++i) { + v_set(m_updateShifts[i], epsilon, m_noteCount); + } + for (int i = 0; i < m_sourceCount; ++i) { + v_set(m_updateSources[i], epsilon, m_noteCount); + } + + double *contributions = allocate<double>(m_binCount); for (int n = 0; n < m_noteCount; ++n) { - newPitches[n] = epsilon; - if (n >= m_lowestPitch && n <= m_highestPitch) { - for (int i = 0; i < m_instrumentCount; ++i) { - for (int f = 0; f < m_shiftCount; ++f) { - const float *w = templateFor(i, n, f); - double pitch = m_pitches[n]; - double source = m_sources[i][n]; - double shift = m_shifts[f][n]; - for (int j = 0; j < m_binCount; ++j) { - newPitches[n] += w[j] * m_q[j] * pitch * source * shift; + + const double pitch = m_pitches[n]; + + for (int f = 0; f < m_shiftCount; ++f) { + + const double shift = m_shifts[f][n]; + + for (int i = 0; i < m_sourceCount; ++i) { + + const double source = m_sources[i][n]; + const double factor = pitch * source * shift; + const double *w = templateFor(i, n, f); + + v_copy(contributions, w, m_binCount); + v_multiply(contributions, m_q, m_binCount); + v_scale(contributions, factor, m_binCount); + + double total = v_sum(contributions, m_binCount); + + if (n >= m_lowestPitch && n <= m_highestPitch) { + + m_updatePitches[n] += total; + + if (inRange(i, n)) { + m_updateSources[i][n] += total; } } - } - } - if (m_pitchSparsity != 1.0) { - newPitches[n] = pow(newPitches[n], m_pitchSparsity); - } - } - normaliseColumn(newPitches); - Grid newShifts = m_shifts; - - if (m_useShifts) { - for (int f = 0; f < m_shiftCount; ++f) { - for (int n = 0; n < m_noteCount; ++n) { - newShifts[f][n] = epsilon; - for (int i = 0; i < m_instrumentCount; ++i) { - const float *w = templateFor(i, n, f); - double pitch = m_pitches[n]; - double source = m_sources[i][n]; - double shift = m_shifts[f][n]; - for (int j = 0; j < m_binCount; ++j) { - newShifts[f][n] += w[j] * m_q[j] * pitch * source * shift; - } - } - } - } - normaliseGrid(newShifts); - } - - Grid newSources = m_sources; - - for (int i = 0; i < m_instrumentCount; ++i) { - for (int n = 0; n < m_noteCount; ++n) { - newSources[i][n] = epsilon; - if (inRange(i, n)) { - for (int f = 0; f < m_shiftCount; ++f) { - const float *w = templateFor(i, n, f); - double pitch = m_pitches[n]; - double source = m_sources[i][n]; - double shift = m_shifts[f][n]; - for (int j = 0; j < m_binCount; ++j) { - newSources[i][n] += w[j] * m_q[j] * pitch * source * shift; - } - } - } - if (m_sourceSparsity != 1.0) { - newSources[i][n] = pow(newSources[i][n], m_sourceSparsity); + m_updateShifts[f][n] += total; } } } - normaliseGrid(newSources); - m_pitches = newPitches; - m_shifts = newShifts; - m_sources = newSources; + if (m_pitchSparsity != 1.0) { + for (int n = 0; n < m_noteCount; ++n) { + m_updatePitches[n] = + pow(m_updatePitches[n], m_pitchSparsity); + } + } + + if (m_sourceSparsity != 1.0) { + for (int n = 0; n < m_noteCount; ++n) { + for (int i = 0; i < m_sourceCount; ++i) { + m_updateSources[i][n] = + pow(m_updateSources[i][n], m_sourceSparsity); + } + } + } + + normaliseColumn(m_updatePitches, m_noteCount); + std::swap(m_pitches, m_updatePitches); + + if (m_useShifts) { + normaliseGrid(m_updateShifts, m_shiftCount, m_noteCount); + std::swap(m_shifts, m_updateShifts); + } + + normaliseGrid(m_updateSources, m_sourceCount, m_noteCount); + std::swap(m_sources, m_updateSources); }
--- a/src/EM.h Tue May 06 18:55:11 2014 +0100 +++ b/src/EM.h Wed May 07 09:01:32 2014 +0100 @@ -24,48 +24,54 @@ EM(bool useShifts); ~EM(); - void iterate(std::vector<double> column); + int getBinCount() const { return m_binCount; } // size of input column + int getNoteCount() const { return m_noteCount; } // size of pitch column + int getSourceCount() const { return m_sourceCount; } - const std::vector<double> &getEstimate() const { + void iterate(const double *column); + + const double *getEstimate() const { // bin count return m_estimate; } - const std::vector<double> &getPitchDistribution() const { + const double *getPitchDistribution() const { // note count return m_pitches; } - const std::vector<std::vector<double> > &getSources() const { + const double *const *getSources() const { // source count * note count return m_sources; } private: - typedef std::vector<double> V; - typedef std::vector<std::vector<double> > Grid; + double *m_pitches; + double **m_shifts; + double **m_sources; + + double *m_updatePitches; + double **m_updateShifts; + double **m_updateSources; bool m_useShifts; - V m_pitches; - Grid m_shifts; - Grid m_sources; + double *m_estimate; + double *m_q; + + const int m_noteCount; + const int m_shiftCount; // 1 + 2 * max template shift + const int m_binCount; + const int m_sourceCount; + + const double m_pitchSparsity; + const double m_sourceSparsity; - V m_estimate; - V m_q; - - int m_noteCount; - int m_shiftCount; // 1 + 2 * max template shift - int m_binCount; - int m_instrumentCount; - - double m_pitchSparsity; - double m_sourceSparsity; + const int m_lowestPitch; + const int m_highestPitch; - int m_lowestPitch; - int m_highestPitch; + void normaliseColumn(double *column, int size); + void normaliseGrid(double **grid, int size1, int size2); - void normaliseColumn(V &column); - void normaliseGrid(Grid &grid); - void expectation(const V &column); - void maximisation(const V &column); + void expectation(const double *column); // size is m_binCount + void maximisation(const double *column); // size is m_binCount - const float *templateFor(int instrument, int note, int shift); + const double *templateFor(int instrument, int note, int shift); void rangeFor(int instrument, int &minPitch, int &maxPitch); bool inRange(int instrument, int pitch); };
--- a/src/Silvet.cpp Tue May 06 18:55:11 2014 +0100 +++ b/src/Silvet.cpp Wed May 07 09:01:32 2014 +0100 @@ -397,6 +397,8 @@ FeatureSet fs; + if (filtered.empty()) return fs; + for (int i = 0; i < (int)filtered.size(); ++i) { Feature f; for (int j = 0; j < processingHeight; ++j) { @@ -409,37 +411,56 @@ int iterations = 12; - for (int i = 0; i < width; ++i) { + int stride = 8; - double sum = 0.0; - for (int j = 0; j < processingHeight; ++j) { - sum += filtered[i][j]; + for (int i = 0; i < width; i += stride) { + + int chunk = stride; + if (i + chunk > width) { + chunk = width - i; } - if (sum < 1e-5) continue; + vector<vector<double> > pitchSubMatrix + (chunk, vector<double>(processingNotes)); - EM em(m_hqMode); - for (int j = 0; j < iterations; ++j) { - em.iterate(filtered[i]); +#pragma omp parallel for + for (int k = 0; k < chunk; ++k) { + + double sum = 0.0; + for (int j = 0; j < processingHeight; ++j) { + sum += filtered[i + k][j]; + } + + if (sum < 1e-5) continue; + + EM em(m_hqMode); + for (int j = 0; j < iterations; ++j) { + em.iterate(filtered[i + k].data()); + } + + const double *pitches = em.getPitchDistribution(); + + for (int j = 0; j < processingNotes; ++j) { + pitchSubMatrix[k][j] = pitches[j] * sum; + } } + + for (int k = 0; k < chunk; ++k) { - vector<double> pitches = em.getPitchDistribution(); - - for (int j = 0; j < processingNotes; ++j) { - pitches[j] *= sum; - } + const vector<double> &pitches = pitchSubMatrix[k]; - Feature f; - for (int j = 0; j < processingNotes; ++j) { - f.values.push_back(float(pitches[j])); - } - fs[m_pitchOutputNo].push_back(f); + Feature f; + for (int j = 0; j < processingNotes; ++j) { + f.values.push_back(float(pitches[j])); + } + fs[m_pitchOutputNo].push_back(f); - FeatureList noteFeatures = postProcess(pitches); + FeatureList noteFeatures = postProcess(pitches); - for (FeatureList::const_iterator fi = noteFeatures.begin(); - fi != noteFeatures.end(); ++fi) { - fs[m_notesOutputNo].push_back(*fi); + for (FeatureList::const_iterator fi = noteFeatures.begin(); + fi != noteFeatures.end(); ++fi) { + fs[m_notesOutputNo].push_back(*fi); + } } }
--- a/testdata/timing/results.txt Tue May 06 18:55:11 2014 +0100 +++ b/testdata/timing/results.txt Wed May 07 09:01:32 2014 +0100 @@ -91,7 +91,6 @@ conclusion: supports the previous test - OPENMP: commit:62b7be1226d5, as commit:ce64d11ef336 but with OpenMP parallel
--- a/yeti/scratch/generateTemplatesC.yeti Tue May 06 18:55:11 2014 +0100 +++ b/yeti/scratch/generateTemplatesC.yeti Wed May 07 09:01:32 2014 +0100 @@ -106,7 +106,7 @@ " const char *name;", " int lowest;", " int highest;", - " float data[SILVET_TEMPLATE_NOTE_COUNT][SILVET_TEMPLATE_SIZE];", + " double data[SILVET_TEMPLATE_NOTE_COUNT][SILVET_TEMPLATE_SIZE];", "} silvet_template_t;", "", "static int silvet_templates_lowest_note = \(overallLowest);",