Chris@7: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@7: 
Chris@7: /*
Chris@7:   Vamp feature extraction plugin for the BeatRoot beat tracker.
Chris@7: 
Chris@7:   Centre for Digital Music, Queen Mary, University of London.
Chris@7:   This file copyright 2011 Simon Dixon, Chris Cannam and QMUL.
Chris@7:     
Chris@7:   This program is free software; you can redistribute it and/or
Chris@7:   modify it under the terms of the GNU General Public License as
Chris@7:   published by the Free Software Foundation; either version 2 of the
Chris@7:   License, or (at your option) any later version.  See the file
Chris@7:   COPYING included with this distribution for more information.
Chris@7: */
Chris@7: 
Chris@9: #include "Induction.h"
Chris@9: 
Chris@7: double Induction::clusterWidth = 0.025;
Chris@7: double Induction::minIOI = 0.070;
Chris@7: double Induction::maxIOI = 2.500;
Chris@7: double Induction::minIBI = 0.3; 
Chris@7: double Induction::maxIBI = 1.0;
Chris@7: int Induction::topN = 10;
Chris@7: 
Chris@15: 
Chris@23: AgentList Induction::beatInduction(AgentParameters params, EventList events) {
Chris@15:     int i, j, b, bestCount;
Chris@15:     bool submult;
Chris@15:     int intervals = 0;			// number of interval clusters
Chris@15:     vector<int> bestn;// count of high-scoring clusters
Chris@15:     bestn.resize(topN);
Chris@15: 
Chris@15:     double ratio, err;
Chris@15:     int degree;
Chris@15:     int maxClusterCount = (int) ceil((maxIOI - minIOI) / clusterWidth);
Chris@15:     vector<double> clusterMean;
Chris@15:     clusterMean.resize(maxClusterCount);
Chris@15:     vector<int> clusterSize;
Chris@15:     clusterSize.resize(maxClusterCount);
Chris@15:     vector<int> clusterScore;
Chris@15:     clusterScore.resize(maxClusterCount);
Chris@15: 		
Chris@15:     EventList::iterator ptr1, ptr2;
Chris@15:     Event e1, e2;
Chris@15:     ptr1 = events.begin();
Chris@15:     while (ptr1 != events.end()) {
Chris@15:         e1 = *ptr1;
Chris@15:         ++ptr1;
Chris@15:         ptr2 = events.begin();
Chris@15:         e2 = *ptr2;
Chris@15:         ++ptr2;
Chris@15:         while (e2 != e1 && ptr2 != events.end()) {
Chris@15:             e2 = *ptr2;
Chris@15:             ++ptr2;
Chris@15:         }
Chris@15:         while (ptr2 != events.end()) {
Chris@15:             e2 = *ptr2;
Chris@15:             ++ptr2;
Chris@15:             double ioi = e2.time - e1.time;
Chris@15:             if (ioi < minIOI)		// skip short intervals
Chris@15:                 continue;
Chris@15:             if (ioi > maxIOI)		// ioi too long
Chris@15:                 break;
Chris@15:             for (b = 0; b < intervals; b++)		// assign to nearest cluster
Chris@15:                 if (fabs(clusterMean[b] - ioi) < clusterWidth) {
Chris@15:                     if ((b < intervals - 1) && (
Chris@15:                             fabs(clusterMean[b+1] - ioi) <
Chris@15:                             fabs(clusterMean[b] - ioi)))
Chris@15:                         b++;		// next cluster is closer
Chris@15:                     clusterMean[b] = (clusterMean[b] * clusterSize[b] +ioi)/
Chris@15:                         (clusterSize[b] + 1);
Chris@15:                     clusterSize[b]++;
Chris@15:                     break;
Chris@15:                 }
Chris@15:             if (b == intervals) {	// no suitable cluster; create new one
Chris@15:                 if (intervals == maxClusterCount) {
Chris@15: //			System.err.println("Warning: Too many clusters");
Chris@15:                     continue; // ignore this IOI
Chris@15:                 }
Chris@15:                 intervals++;
Chris@15:                 for ( ; (b>0) && (clusterMean[b-1] > ioi); b--) {
Chris@15:                     clusterMean[b] = clusterMean[b-1];
Chris@15:                     clusterSize[b] = clusterSize[b-1];
Chris@15:                 }
Chris@15:                 clusterMean[b] = ioi;
Chris@15:                 clusterSize[b] = 1;
Chris@15:             }
Chris@15:         }
Chris@15:     }
Chris@15:     for (b = 0; b < intervals; b++)	// merge similar intervals
Chris@15:         // TODO: they are now in order, so don't need the 2nd loop
Chris@15:         // TODO: check BOTH sides before averaging or upper gps don't work
Chris@15:         for (i = b+1; i < intervals; i++)
Chris@15:             if (fabs(clusterMean[b] - clusterMean[i]) < clusterWidth) {
Chris@15:                 clusterMean[b] = (clusterMean[b] * clusterSize[b] +
Chris@15:                                   clusterMean[i] * clusterSize[i]) /
Chris@15:                     (clusterSize[b] + clusterSize[i]);
Chris@15:                 clusterSize[b] = clusterSize[b] + clusterSize[i];
Chris@15:                 --intervals;
Chris@15:                 for (j = i+1; j <= intervals; j++) {
Chris@15:                     clusterMean[j-1] = clusterMean[j];
Chris@15:                     clusterSize[j-1] = clusterSize[j];
Chris@15:                 }
Chris@15:             }
Chris@15:     if (intervals == 0)
Chris@15:         return AgentList();
Chris@15:     for (b = 0; b < intervals; b++)
Chris@15:         clusterScore[b] = 10 * clusterSize[b];
Chris@15:     bestn[0] = 0;
Chris@15:     bestCount = 1;
Chris@15:     for (b = 0; b < intervals; b++)
Chris@15:         for (i = 0; i <= bestCount; i++)
Chris@15:             if ((i < topN) && ((i == bestCount) ||
Chris@15:                                (clusterScore[b] > clusterScore[bestn[i]]))){
Chris@15:                 if (bestCount < topN)
Chris@15:                     bestCount++;
Chris@15:                 for (j = bestCount - 1; j > i; j--)
Chris@15:                     bestn[j] = bestn[j-1];
Chris@15:                 bestn[i] = b;
Chris@15:                 break;
Chris@15:             }
Chris@15:     for (b = 0; b < intervals; b++)	// score intervals
Chris@15:         for (i = b+1; i < intervals; i++) {
Chris@15:             ratio = clusterMean[b] / clusterMean[i];
Chris@15:             submult = ratio < 1;
Chris@15:             if (submult)
Chris@15:                 degree = (int) nearbyint(1/ratio);
Chris@15:             else
Chris@15:                 degree = (int) nearbyint(ratio);
Chris@15:             if ((degree >= 2) && (degree <= 8)) {
Chris@15:                 if (submult)
Chris@15:                     err = fabs(clusterMean[b]*degree - clusterMean[i]);
Chris@15:                 else
Chris@15:                     err = fabs(clusterMean[b] - clusterMean[i]*degree);
Chris@15:                 if (err < (submult? clusterWidth : clusterWidth * degree)) {
Chris@15:                     if (degree >= 5)
Chris@15:                         degree = 1;
Chris@15:                     else
Chris@15:                         degree = 6 - degree;
Chris@15:                     clusterScore[b] += degree * clusterSize[i];
Chris@15:                     clusterScore[i] += degree * clusterSize[b];
Chris@15:                 }
Chris@15:             }
Chris@15:         }
Chris@15: 
Chris@15:     AgentList a;
Chris@15:     for (int index = 0; index < bestCount; index++) {
Chris@15:         b = bestn[index];
Chris@15:         // Adjust it, using the size of super- and sub-intervals
Chris@15:         double newSum = clusterMean[b] * clusterScore[b];
Chris@15:         int newCount = clusterSize[b];
Chris@15:         int newWeight = clusterScore[b];
Chris@15:         for (i = 0; i < intervals; i++) {
Chris@15:             if (i == b)
Chris@15:                 continue;
Chris@15:             ratio = clusterMean[b] / clusterMean[i];
Chris@15:             if (ratio < 1) {
Chris@15:                 degree = (int) nearbyint(1 / ratio);
Chris@15:                 if ((degree >= 2) && (degree <= 8)) {
Chris@15:                     err = fabs(clusterMean[b]*degree - clusterMean[i]);
Chris@15:                     if (err < clusterWidth) {
Chris@15:                         newSum += clusterMean[i] / degree * clusterScore[i];
Chris@15:                         newCount += clusterSize[i];
Chris@15:                         newWeight += clusterScore[i];
Chris@15:                     }
Chris@15:                 }
Chris@15:             } else {
Chris@15:                 degree = (int) nearbyint(ratio);
Chris@15:                 if ((degree >= 2) && (degree <= 8)) {
Chris@15:                     err = fabs(clusterMean[b] - degree*clusterMean[i]);
Chris@15:                     if (err < clusterWidth * degree) {
Chris@15:                         newSum += clusterMean[i] * degree * clusterScore[i];
Chris@15:                         newCount += clusterSize[i];
Chris@15:                         newWeight += clusterScore[i];
Chris@15:                     }
Chris@15:                 }
Chris@15:             }
Chris@15:         }
Chris@15:         double beat = newSum / newWeight;
Chris@15:         // Scale within range ... hope the grouping isn't ternary :(
Chris@15:         while (beat < minIBI)		// Maximum speed
Chris@15:             beat *= 2.0;
Chris@15:         while (beat > maxIBI)		// Minimum speed
Chris@15:             beat /= 2.0;
Chris@15:         if (beat >= minIBI) {
Chris@23:             a.push_back(new Agent(params, beat));
Chris@15:         }
Chris@15:     }
Chris@15: #ifdef DEBUG_BEATROOT
Chris@15:     std::cerr << "Induction complete, returning " << a.size() << " agent(s)" << std::endl;
Chris@15: #endif
Chris@15:     return a;
Chris@15: } // beatInduction()
Chris@15: