BeatRoot Vamp Plugin

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

/*

  Vamp feature extraction plugin for the BeatRoot beat tracker.

  Centre for Digital Music, Queen Mary, University of London.

  This file copyright 2011 Simon Dixon, Chris Cannam and QMUL.

  This program is free software; you can redistribute it and/or

  modify it under the terms of the GNU General Public License as

  published by the Free Software Foundation; either version 2 of the

  License, or (at your option) any later version.  See the file

  COPYING included with this distribution for more information.

*/

#include "Induction.h"

double Induction::clusterWidth = 0.025;

double Induction::minIOI = 0.070;

double Induction::maxIOI = 2.500;

double Induction::minIBI = 0.3; 

double Induction::maxIBI = 1.0;

int Induction::topN = 10;

AgentList Induction::beatInduction(AgentParameters params, EventList events) {

    int i, j, b, bestCount;

    bool submult;

    int intervals = 0;                        // number of interval clusters

    vector<int> bestn;// count of high-scoring clusters

    bestn.resize(topN);

    double ratio, err;

    int degree;

    int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth);

    vector<double> clusterMean;

    clusterMean.resize(maxClusterCount);

    vector<int> clusterSize;

    clusterSize.resize(maxClusterCount);

    vector<int> clusterScore;

    clusterScore.resize(maxClusterCount);

    EventList::iterator ptr1, ptr2;

    Event e1, e2;

    ptr1 = events.begin();

    while (ptr1 != events.end()) {

        e1 = *ptr1;

        ++ptr1;

        ptr2 = events.begin();

        e2 = *ptr2;

        ++ptr2;

        while (e2 != e1 && ptr2 != events.end()) {

            e2 = *ptr2;

            ++ptr2;

        }

        while (ptr2 != events.end()) {

            e2 = *ptr2;

            ++ptr2;

            double ioi = e2.time - e1.time;

            if (ioi < minIOI)                // skip short intervals

                continue;

            if (ioi > maxIOI)                // ioi too long

                break;

            for (b = 0; b < intervals; b++)                // assign to nearest cluster

                if (fabs(clusterMean[b] - ioi) < clusterWidth) {

                    if ((b < intervals - 1) && (

                            fabs(clusterMean[b+1] - ioi) <

                            fabs(clusterMean[b] - ioi)))

                        b++;                // next cluster is closer

                    clusterMean[b] = (clusterMean[b] * clusterSize[b] +ioi)/

                        (clusterSize[b] + 1);

                    clusterSize[b]++;

                    break;

                }

            if (b == intervals) {        // no suitable cluster; create new one

                if (intervals == maxClusterCount) {

//                        System.err.println("Warning: Too many clusters");

                    continue; // ignore this IOI

                }

                intervals++;

                for ( ; (b>0) && (clusterMean[b-1] > ioi); b--) {

                    clusterMean[b] = clusterMean[b-1];

                    clusterSize[b] = clusterSize[b-1];

                }

                clusterMean[b] = ioi;

                clusterSize[b] = 1;

            }

        }

    }

    for (b = 0; b < intervals; b++)        // merge similar intervals

        // TODO: they are now in order, so don't need the 2nd loop

        // TODO: check BOTH sides before averaging or upper gps don't work

        for (i = b+1; i < intervals; i++)

            if (fabs(clusterMean[b] - clusterMean[i]) < clusterWidth) {

                clusterMean[b] = (clusterMean[b] * clusterSize[b] +

                                  clusterMean[i] * clusterSize[i]) /

                    (clusterSize[b] + clusterSize[i]);

                clusterSize[b] = clusterSize[b] + clusterSize[i];

                --intervals;

                for (j = i+1; j <= intervals; j++) {

                    clusterMean[j-1] = clusterMean[j];

                    clusterSize[j-1] = clusterSize[j];

                }

            }

    if (intervals == 0)

        return AgentList();

    for (b = 0; b < intervals; b++)

        clusterScore[b] = 10 * clusterSize[b];

    bestn[0] = 0;

    bestCount = 1;

    for (b = 0; b < intervals; b++)

        for (i = 0; i <= bestCount; i++)

            if ((i < topN) && ((i == bestCount) ||

                               (clusterScore[b] > clusterScore[bestn[i]]))){

                if (bestCount < topN)

                    bestCount++;

                for (j = bestCount - 1; j > i; j--)

                    bestn[j] = bestn[j-1];

                bestn[i] = b;

                break;

            }

    for (b = 0; b < intervals; b++)        // score intervals

        for (i = b+1; i < intervals; i++) {

            ratio = clusterMean[b] / clusterMean[i];

            submult = ratio < 1;

            if (submult)

                degree = (int) nearbyint(1/ratio);

            else

                degree = (int) nearbyint(ratio);

            if ((degree >= 2) && (degree <= 8)) {

                if (submult)

                    err = fabs(clusterMean[b]*degree - clusterMean[i]);

                else

                    err = fabs(clusterMean[b] - clusterMean[i]*degree);

                if (err < (submult? clusterWidth : clusterWidth * degree)) {

                    if (degree >= 5)

                        degree = 1;

                    else

                        degree = 6 - degree;

                    clusterScore[b] += degree * clusterSize[i];

                    clusterScore[i] += degree * clusterSize[b];

                }

            }

        }

    AgentList a;

    for (int index = 0; index < bestCount; index++) {

        b = bestn[index];

        // Adjust it, using the size of super- and sub-intervals

        double newSum = clusterMean[b] * clusterScore[b];

        int newCount = clusterSize[b];

        int newWeight = clusterScore[b];

        for (i = 0; i < intervals; i++) {

            if (i == b)

                continue;

            ratio = clusterMean[b] / clusterMean[i];

            if (ratio < 1) {

                degree = (int) nearbyint(1 / ratio);

                if ((degree >= 2) && (degree <= 8)) {

                    err = fabs(clusterMean[b]*degree - clusterMean[i]);

                    if (err < clusterWidth) {

                        newSum += clusterMean[i] / degree * clusterScore[i];

                        newCount += clusterSize[i];

                        newWeight += clusterScore[i];

                    }

                }

            } else {

                degree = (int) nearbyint(ratio);

                if ((degree >= 2) && (degree <= 8)) {

                    err = fabs(clusterMean[b] - degree*clusterMean[i]);

                    if (err < clusterWidth * degree) {

                        newSum += clusterMean[i] * degree * clusterScore[i];

                        newCount += clusterSize[i];

                        newWeight += clusterScore[i];

                    }

                }

            }

        }

        double beat = newSum / newWeight;

        // Scale within range ... hope the grouping isn't ternary :(

        while (beat < minIBI)                // Maximum speed

            beat *= 2.0;

        while (beat > maxIBI)                // Minimum speed

            beat /= 2.0;

        if (beat >= minIBI) {

            a.push_back(new Agent(params, beat));

        }

    }

#ifdef DEBUG_BEATROOT

    std::cerr << "Induction complete, returning " << a.size() << " agent(s)" << std::endl;

#endif

    return a;

} // beatInduction()