Chris@34: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@34: 
Chris@34: /*
Chris@34:   Silvet
Chris@34: 
Chris@34:   A Vamp plugin for note transcription.
Chris@34:   Centre for Digital Music, Queen Mary University of London.
Chris@34:     
Chris@34:   This program is free software; you can redistribute it and/or
Chris@34:   modify it under the terms of the GNU General Public License as
Chris@34:   published by the Free Software Foundation; either version 2 of the
Chris@34:   License, or (at your option) any later version.  See the file
Chris@34:   COPYING included with this distribution for more information.
Chris@34: */
Chris@34: 
Chris@34: #include "EM.h"
Chris@34: 
Chris@36: #include <cstdlib>
Chris@42: #include <cmath>
Chris@36: 
Chris@36: #include <iostream>
Chris@36: 
Chris@310: #include "bqvec/VectorOps.h"
Chris@310: #include "bqvec/Allocators.h"
Chris@161: #include "Instruments.h"
Chris@36: 
Chris@36: using std::vector;
Chris@36: using std::cerr;
Chris@36: using std::endl;
Chris@36: 
Chris@91: using namespace breakfastquay;
Chris@91: 
Chris@151: static float epsilon = 1e-10;
Chris@35: 
Chris@161: EM::EM(const InstrumentPack *pack, bool useShifts) :
Chris@161:     m_pack(pack),
Chris@161:     m_noteCount(pack->templateNoteCount),
Chris@161:     m_shiftCount(useShifts ? pack->templateMaxShift * 2 + 1 : 1),
Chris@161:     m_binCount(pack->templateHeight),
Chris@161:     m_sourceCount(pack->templates.size()),
Chris@42:     m_pitchSparsity(1.1),
Chris@174:     m_shiftSparsity(1.1),
Chris@273:     m_sourceSparsity(1.2)
Chris@35: {
Chris@151:     m_pitches = allocate<float>(m_noteCount);
Chris@151:     m_updatePitches = allocate<float>(m_noteCount);
Chris@338:     v_set(m_pitches, 1.f, m_noteCount);
Chris@35: 
Chris@113:     if (useShifts) {
Chris@151:         m_shifts = allocate_channels<float>(m_shiftCount, m_noteCount);
Chris@151:         m_updateShifts = allocate_channels<float>(m_shiftCount, m_noteCount);
Chris@113:         for (int f = 0; f < m_shiftCount; ++f) {
Chris@338:             v_set(m_shifts[f], 1.f, m_noteCount);
Chris@55:         }
Chris@113:     } else {
Chris@113:         m_shifts = 0;
Chris@113:         m_updateShifts = 0;
Chris@35:     }
Chris@35:     
Chris@151:     m_sources = allocate_channels<float>(m_sourceCount, m_noteCount);
Chris@151:     m_updateSources = allocate_channels<float>(m_sourceCount, m_noteCount);
Chris@91:     for (int i = 0; i < m_sourceCount; ++i) {
Chris@55:         for (int n = 0; n < m_noteCount; ++n) {
Chris@338:             m_sources[i][n] = (inRange(i, n) ? 1.f : 0.f);
Chris@35:         }
Chris@35:     }
Chris@35: 
Chris@151:     m_estimate = allocate<float>(m_binCount);
Chris@151:     m_q = allocate<float>(m_binCount);
Chris@35: }
Chris@35: 
Chris@35: EM::~EM()
Chris@35: {
Chris@92:     deallocate(m_q);
Chris@92:     deallocate(m_estimate);
Chris@92:     deallocate_channels(m_sources, m_sourceCount);
Chris@100:     deallocate_channels(m_updateSources, m_sourceCount);
Chris@92:     deallocate_channels(m_shifts, m_shiftCount);
Chris@100:     deallocate_channels(m_updateShifts, m_shiftCount);
Chris@92:     deallocate(m_pitches);
Chris@100:     deallocate(m_updatePitches);
Chris@35: }
Chris@35: 
Chris@45: void
Chris@45: EM::rangeFor(int instrument, int &minPitch, int &maxPitch)
Chris@45: {
Chris@161:     minPitch = m_pack->templates[instrument].lowestNote;
Chris@161:     maxPitch = m_pack->templates[instrument].highestNote;
Chris@45: }
Chris@45: 
Chris@35: bool
Chris@45: EM::inRange(int instrument, int pitch)
Chris@35: {
Chris@45:     int minPitch, maxPitch;
Chris@45:     rangeFor(instrument, minPitch, maxPitch);
Chris@45:     return (pitch >= minPitch && pitch <= maxPitch);
Chris@35: }
Chris@35: 
Chris@36: void
Chris@151: EM::normaliseColumn(float *column, int size)
Chris@36: {
Chris@151:     float sum = v_sum(column, size);
Chris@92:     v_scale(column, 1.0 / sum, size);
Chris@36: }
Chris@36: 
Chris@36: void
Chris@151: EM::normaliseGrid(float **grid, int size1, int size2)
Chris@53: {
Chris@151:     float *denominators = allocate_and_zero<float>(size2);
Chris@122: 
Chris@92:     for (int i = 0; i < size1; ++i) {
Chris@122:         for (int j = 0; j < size2; ++j) {
Chris@122:             denominators[j] += grid[i][j];
Chris@122:         }
Chris@53:     }
Chris@122: 
Chris@122:     for (int i = 0; i < size1; ++i) {
Chris@122:         v_divide(grid[i], denominators, size2);
Chris@122:     }
Chris@122: 
Chris@122:     deallocate(denominators);
Chris@53: }
Chris@53: 
Chris@53: void
Chris@92: EM::iterate(const double *column)
Chris@36: {
Chris@151:     float *norm = allocate<float>(m_binCount);
Chris@151:     v_convert(norm, column, m_binCount);
Chris@92:     normaliseColumn(norm, m_binCount);
Chris@92:     expectation(norm);
Chris@354:     maximisation();
Chris@95:     deallocate(norm);
Chris@36: }
Chris@36: 
Chris@151: const float *
Chris@55: EM::templateFor(int instrument, int note, int shift)
Chris@45: {
Chris@161:     const float *base = m_pack->templates.at(instrument).data.at(note).data();
Chris@113:     if (m_shifts) {
Chris@161:         return base + shift;
Chris@110:     } else {
Chris@161:         return base + m_pack->templateMaxShift;
Chris@110:     }
Chris@45: }
Chris@45: 
Chris@36: void
Chris@151: EM::expectation(const float *column)
Chris@36: {
Chris@62: //    cerr << ".";
Chris@36: 
Chris@99:     v_set(m_estimate, epsilon, m_binCount);
Chris@36: 
Chris@130:     for (int f = 0; f < m_shiftCount; ++f) {
Chris@130: 
Chris@151:         const float *shiftIn = m_shifts ? m_shifts[f] : 0;
Chris@130: 
Chris@130:         for (int i = 0; i < m_sourceCount; ++i) {
Chris@130: 
Chris@151:             const float *sourceIn = m_sources[i];
Chris@130: 
Chris@130:             int lowest, highest;
Chris@130:             rangeFor(i, lowest, highest);
Chris@130: 
Chris@130:             for (int n = lowest; n <= highest; ++n) {
Chris@130: 
Chris@151:                 const float source = sourceIn[n];
Chris@151:                 const float shift = shiftIn ? shiftIn[n] : 1.0;
Chris@151:                 const float pitch = m_pitches[n];
Chris@130: 
Chris@151:                 const float factor = pitch * source * shift;
Chris@151:                 const float *w = templateFor(i, n, f);
Chris@130: 
Chris@111:                 v_add_with_gain(m_estimate, w, factor, m_binCount);
Chris@36:             }
Chris@36:         }
Chris@36:     }
Chris@36: 
Chris@45:     for (int i = 0; i < m_binCount; ++i) {
Chris@36:         m_q[i] = column[i] / m_estimate[i];
Chris@36:     }
Chris@164: 
Chris@164: /*
Chris@164:     double l2norm = 0.0;
Chris@164:     
Chris@164:     for (int i = 0; i < m_binCount; ++i) {
Chris@164:         l2norm += (column[i] - m_estimate[i]) * (column[i] - m_estimate[i]);
Chris@164:     }
Chris@164: 
Chris@164:     l2norm = sqrt(l2norm);
Chris@164:     cerr << "l2norm = " << l2norm << endl;
Chris@164: */
Chris@36: }
Chris@36: 
Chris@36: void
Chris@354: EM::maximisation()
Chris@36: {
Chris@100:     v_set(m_updatePitches, epsilon, m_noteCount);
Chris@113: 
Chris@92:     for (int i = 0; i < m_sourceCount; ++i) {
Chris@100:         v_set(m_updateSources[i], epsilon, m_noteCount);
Chris@92:     }
Chris@62: 
Chris@113:     if (m_shifts) {
Chris@113:         for (int i = 0; i < m_shiftCount; ++i) {
Chris@113:             v_set(m_updateShifts[i], epsilon, m_noteCount);
Chris@113:         }
Chris@113:     }
Chris@113: 
Chris@151:     float *contributions = allocate<float>(m_binCount);
Chris@36: 
Chris@130:     for (int f = 0; f < m_shiftCount; ++f) {
Chris@85: 
Chris@151:         const float *shiftIn = m_shifts ? m_shifts[f] : 0;
Chris@151:         float *shiftOut = m_shifts ? m_updateShifts[f] : 0;
Chris@85: 
Chris@130:         for (int i = 0; i < m_sourceCount; ++i) {
Chris@85: 
Chris@151:             const float *sourceIn = m_sources[i];
Chris@151:             float *sourceOut = m_updateSources[i];
Chris@85: 
Chris@130:             int lowest, highest;
Chris@130:             rangeFor(i, lowest, highest);
Chris@85: 
Chris@130:             for (int n = lowest; n <= highest; ++n) {
Chris@130: 
Chris@151:                 const float shift = shiftIn ? shiftIn[n] : 1.0;
Chris@151:                 const float source = sourceIn[n];
Chris@151:                 const float pitch = m_pitches[n];
Chris@130: 
Chris@151:                 const float factor = pitch * source * shift;
Chris@151:                 const float *w = templateFor(i, n, f);
Chris@85: 
Chris@94:                 v_copy(contributions, w, m_binCount);
Chris@95:                 v_multiply(contributions, m_q, m_binCount);
Chris@94: 
Chris@151:                 float total = factor * v_sum(contributions, m_binCount);
Chris@94: 
Chris@130:                 m_updatePitches[n] += total;
Chris@130:                 sourceOut[n] += total;
Chris@85: 
Chris@130:                 if (shiftOut) {
Chris@130:                     shiftOut[n] += total;
Chris@113:                 }
Chris@42:             }
Chris@36:         }
Chris@36:     }
Chris@36: 
Chris@359:     deallocate(contributions);
Chris@359:     
Chris@103:     if (m_pitchSparsity != 1.0) {
Chris@103:         for (int n = 0; n < m_noteCount; ++n) {
Chris@103:             m_updatePitches[n] = 
Chris@174:                 powf(m_updatePitches[n], m_pitchSparsity);
Chris@174:         }
Chris@174:     }
Chris@174: 
Chris@174:     if (m_shifts && m_shiftSparsity != 1.0) {
Chris@174:         for (int i = 0; i < m_shiftCount; ++i) {
Chris@174:             for (int n = 0; n < m_noteCount; ++n) {
Chris@174:                 m_updateShifts[i][n] =
Chris@174:                     powf(m_updateShifts[i][n], m_shiftSparsity);
Chris@174:             }
Chris@62:         }
Chris@103:     }
Chris@103: 
Chris@103:     if (m_sourceSparsity != 1.0) {
Chris@130:         for (int i = 0; i < m_sourceCount; ++i) {
Chris@130:             for (int n = 0; n < m_noteCount; ++n) {
Chris@103:                 m_updateSources[i][n] =
Chris@174:                     powf(m_updateSources[i][n], m_sourceSparsity);
Chris@62:             }
Chris@62:         }
Chris@62:     }
Chris@85: 
Chris@100:     normaliseColumn(m_updatePitches, m_noteCount);
Chris@112:     std::swap(m_pitches, m_updatePitches);
Chris@112: 
Chris@113:     normaliseGrid(m_updateSources, m_sourceCount, m_noteCount);
Chris@113:     std::swap(m_sources, m_updateSources);
Chris@113: 
Chris@113:     if (m_shifts) {
Chris@112:         normaliseGrid(m_updateShifts, m_shiftCount, m_noteCount);
Chris@112:         std::swap(m_shifts, m_updateShifts);
Chris@112:     }
Chris@36: }
Chris@36: 
Chris@36: