annotate DrumTimingLoader_OF/PerformanceAnalyserSrc/TimingAnalyser.cpp @ 3:303edbbcf1bd tip

updated ofxAubioOnsetDetection file

/*
 *  TimingAnalyser.cpp
 *  performanceTimingAnalyser
 *
 *  Created by Andrew on 17/12/2011.
 *  Copyright 2011 QMUL. All rights reserved.
 *
 */

#include "TimingAnalyser.h"

//FROM TestApp.cpp
//change these to your preferred filepaths for exporting the data
//timingDataFilepath = "/Users/andrew/outputFile.txt";//here as decoded tempo and phase
//change this to your filepath for storing the final path as output
//can then load it back in easily without reanalysing using this...
//	importFileFromStoredData(fileToLoad);
//processBeatTimesFilepath = "//Users/andrew/processedBeats";//here as a list of beat times
//this latter will be a smoothed version of the original

//RESTRICT PHASE HOP POSSIBLE


//Correct or DELETE this BeatPosition variable - what is it? Why used?


TimingAnalyser::TimingAnalyser(){

	moveMatrixToOptimalCostPosition = true;
	
	//resolution in ms
	tempoScalar = 1.0;
	phaseScalar = 1.0;
	
	phaseCost = 1000.4;//a phase jump is 1.5 (y-x) where y now becomes the predicted position
	tempoCost = 2.8;//a tempo jump of z msec per period is charged at 2z units
	
	drawIOIintervals = false;
	
	drawExpressiveTimingData = true;
	
	printHistory = false;
	mozartTriplets = false;
	
	movementFactor = 0.5;
	drawBPM = true;
	
	//COST VARIABLES - Assuming that a simple error is costing linear (t-x) where x is predicted time

	
	blackAndWhite = true;
	
	setTempoLimits(400);

	
	//index between which we look at the variation
	tempoVariationStartIndex = 4;
	
	
	DoubleVector v;
	v.assign(phaseRange, 0.0);
	logProbabilityMatrix.assign(tempoRange, v);	
	previousLogMatrix = logProbabilityMatrix;
	
	//[tempo][phase]
	phaseMinimumOffset = -1 * phaseScalar*(phaseRange - 1)/2;//offset in millis
	globalTimeOffset = 0;//global offset in millis
	
	lastBeatPosition = 0;


	numberOfPointsPerPage = 80;
	startPoint = 0;
	
	//	initialiseRoutes();
		
	recentPhaseShift = 0;
	recentTempoShift = 0;

/*	meanGlobalTempo = 440;
	tempoMinimumOffset = meanGlobalTempo - (tempoRange-1)/2;
	meanTempoIndex = (tempoRange - 1)/2;
	meanTempo = tempoMinimumOffset + meanTempoIndex;
*/	
	meanPhaseIndex = (phaseRange - 1)/2;
	
	currentBestTempo = meanTempoIndex;
	currentBestPhase = meanPhaseIndex;
	
	printf("Mean tempo initialised at %i and mean index %i\n", meanTempo, meanTempoIndex);
	int	minDrawTempo = minimumTempo + minimumPhaseHop - 4;
	int	maxDrawTempo = maximumTempo + maximumPhaseHop + 4;
	
	ofTrueTypeFont::setGlobalDpi(72);
	timesFont.loadFont("TimesNewRoman.ttf", 18, true, true);
	timesFont.setLineHeight(18.0f);
	timesFont.setLetterSpacing(1.037);
	
}

void TimingAnalyser::initialiseRoutes(){
}

void TimingAnalyser::clearData(){
	timingData.clear();
	basicInterOnsetInterval.clear();
	beatPosition.clear();
}

void TimingAnalyser::setTempoLimits(const int& newGlobalTempo){
	
	meanGlobalTempo = newGlobalTempo;
	tempoMinimumOffset = meanGlobalTempo - (tempoScalar*(tempoRange-1)/2);
	meanTempoIndex = (tempoRange - 1)/2;
	meanTempo = tempoMinimumOffset + meanTempoIndex*tempoScalar;
	
	printf("Setting tempo to %i units i.e. %f ms\n", newGlobalTempo, getTempo(meanTempoIndex));
}


double TimingAnalyser::getPhaseIndexFromEventTime(const double& eventTime){
	return eventTime / phaseScalar;
}

double TimingAnalyser::getTempoInUnits(const double& tempoInMs){
	return tempoInMs / tempoScalar;
}

//NEW UPDATE METHOD
double TimingAnalyser::getBestMinimumCost(const int& newTempoInUnits, const int& newPhaseInUnits, Route& r){
//new method to check all points in last matrix
	//new phase is the phase point from zero in units
	//including any offsets etc
	
	double cost = 0;
	double bestCost = INFINITY;
	Path* lastPath;
	int minPhase = 0;
	int minTempo = 0;
	r.tempoHop = 0;
	r.phaseHop = 0;

	//int tempoHop, phaseHop; no need
	if (pathHistory.size() > 0)
	lastPath = &pathHistory[pathHistory.size()-1]; 
	else{
	//	printf("Initiating first path...\n");
	}
	
	for (int tempo = 0;tempo < tempoRange;tempo++){
	
		for (int phase = 0;phase < phaseRange;phase++){
		//change tempo
			cost = 0;
			int phasePoint = 0;
			
			if (pathHistory.size() > 0){
				cost = (*lastPath).currentRoute[tempo][phase].totalCost;
				cost += tempoCost * abs(newTempoInUnits - (tempo + (*lastPath).tempoOffset));
				
			//our point then projects forwards at the new tempo
			phasePoint = phase + (*lastPath).phaseOffset + newTempoInUnits;
			//phase cost is
				cost += phaseHopCost(newPhaseInUnits - phasePoint);//i.e. phaseCost * abs(newPhaseInUnits - phasePoint);
			}//else it's the first point and all previous cost is zero
			
			if (cost < bestCost){
				bestCost = cost;
				minPhase = phase;
				minTempo = tempo;

				if (pathHistory.size() > 0){
					r.tempoHop = newTempoInUnits - (tempo + (*lastPath).tempoOffset);
					r.phaseHop = newPhaseInUnits - phasePoint;
				} else{
					r.tempoHop = 0;
					r.phaseHop = 0;
				}
				
				r.previousTempoIndex = tempo;
				r.previousPhaseIndex = phase;
				
			}//end if mew minimum
		}
	}
	
/*	printf("best cost for tempo-phase pair %i, %i from tempo %i and phase %i, hop %i, %i is %.2f\n", 
		   newTempoInUnits, newPhaseInUnits, 
		   minTempo, minPhase,
		   r.tempoHop, r.phaseHop, bestCost);
*/	return bestCost;
}


void TimingAnalyser::printCostMatrix(const RouteMatrix& m){
	printf("\n");
	for (int i = 0; i< m.size();i++){
		printf("Tempo %.1f (index %i)\n", getTempo(i), i);
		for (int k = 0;k < m[i].size();k++){
			printf("%.1f (new %.1f + prev %.1f (%i,%i)), ", m[i][k].totalCost, 
				   m[i][k].addedCost,m[i][k].previousCost, 
				   m[i][k].previousTempoIndex, m[i][k].previousPhaseIndex );
		}
		printf("\n");
	}
}

void TimingAnalyser::updatePlayIndex(const double& seconds){
	double millis = seconds * 1000.0;
	
	while (playingIndex < timingData.size() && getNoteTime(playingIndex) < millis)
		playingIndex++;
	
	while (playingIndex > 0 && getNoteTime(playingIndex) > millis)
		playingIndex--;
	
	if (playingIndex > startPoint + numberOfPointsPerPage)
		startPoint += numberOfPointsPerPage;//scrolling with play
	
}

double TimingAnalyser::getNoteTime(const int& index){
		return offsetToFirstPoint + timingData[index][2];//
}

void TimingAnalyser::updateCostToPoint(const int& eventTime, const int& eventBeatLocation){
	
	//	previousLogMatrix = logProbabilityMatrix;//store previous matrix - settings are in pathHistory vector
	printf("\nTiming Analyser: update cost to %i, \n", eventTime);
	//event beat location %i is 1 so don't bother printing
//	printf("old phase range is %.1f to %.1f\n", getPhase(0), getPhase(phaseRange-1));
//	printf("old tempo range is %.1f to %.1f\n", getTempo(0), getTempo(tempoRange-1));
	
	int interval = 1;
	
	//setting the new path offset points
	RouteMatrix r_mat;
	Path* lastPath;
	
	Path p;//the new path point
	p.tempoOffset = tempoMinimumOffset;//for now keep it the same
	
	if (pathHistory.size() > 0){
		lastPath = &pathHistory[pathHistory.size()-1];
		int ioi = (int) getTempoInUnits(eventTime - (*lastPath).eventTimeObserved);
		p.tempoOffset = ioi - ((tempoRange-1)/2);
		tempoMinimumOffset = p.tempoOffset;
		printf("New time %i, last time %i diff %i: IOI is %i\n", 
			   eventTime, (*lastPath).eventTimeObserved,
			   (eventTime - (*lastPath).eventTimeObserved), ioi);
	}
	
	//and in tempo units!
	p.phaseOffset = round(getPhaseIndexFromEventTime(eventTime)) + phaseMinimumOffset;//in units
	//abandoned - now using the offset too - dont need - (phaseRange/2) as in the get phase fn
	printf("NEW PHASE and tempo offsets %i and %i, \n", p.phaseOffset, p.tempoOffset);
	//end set offset pts
	
	for (int tempoIndex = 0;tempoIndex < logProbabilityMatrix.size();tempoIndex++){
		RouteVector r_vec;
		
		for (int phase = 0;phase < logProbabilityMatrix[tempoIndex].size();phase++){
			

			Route r;
			
			double bestCost = getBestMinimumCost(p.tempoOffset + tempoIndex, phase + p.phaseOffset, r);
			double newCost = fabs(getPhaseIndexFromEventTime(eventTime) - (p.phaseOffset + phase));

			//this is the phase position that would lead to this current tempo and phase iwthout any hopping
			//the bestPreviousCost calculates the best oposition including hopping
			
			//then we compare the past cost
			
			r.previousCost = bestCost;	
			
			//this is the best past cost that can get us to current location
			//then we need new cost
			
			r.addedCost = newCost;

			//r.globalPhaseOffset = globalTimeOffset;
			
			logProbabilityMatrix[tempoIndex][phase] = bestCost + newCost;	
			r.totalCost = bestCost + newCost;					
			
			//just to check we get this right!
			r.tempoIndex = tempoIndex;
			r.phaseIndex = phase;
			
			r_vec.push_back(r);
		}
		r_mat.push_back(r_vec);
	}
	//history.push_back(r_mat);
	
	//printHistory(); - defunct as using pathHistory instead
	//printBestPath();
	//printf("Global offset %i ph min off %i\n", globalTimeOffset, phaseMinimumOffset);
	
	//new method using Path class 
	p.currentRoute = r_mat;
	
	globalTimeOffset = p.phaseOffset;
	
	//store the settings for the matrix at this time step
	//p.globalPhaseOffset = globalTimeOffset;
//	p.tempoOffset = tempoMinimumOffset;

	setBestTempoAndPhase(p);
	lastBeatPosition = eventBeatLocation;
	printf(" error %i ", (int)(eventTime - getPhase(currentBestPhase, globalTimeOffset)));
	p.phaseShiftApplied = recentPhaseShift;//this is what we did in the recent step
	p.tempoShiftApplied = recentTempoShift;
	//so is the shift by which current stored matrix and indices differ from the last path
	
	p.eventTimeObserved = eventTime;
	p.costMatrix = logProbabilityMatrix;
	
	if (pathHistory.size() > 0){
		//then can calculate inter-onset interval
		int ioi = eventTime - pathHistory[pathHistory.size()-1].eventTimeObserved;
		basicInterOnsetInterval.push_back(ioi);
		printf("ioi for time %i is %i\n", eventTime, ioi);
	}
	printf("New PATH offset is %i\n", p.phaseOffset);
	pathHistory.push_back(p);
	
	//prints the complete cost matrix at every step
	//printCostMatrix(p.currentRoute);
	
	
	if (printHistory)
		printPathHistory();
	
	//	printf(", SHIFT (%i, %i)\n", p.phaseShiftApplied,  p.tempoShiftApplied);
	
	//loaded the new prob matrix to the main log one
	//now looking to shift
	
	//change global offset to match new position
	//	printLogMatrix();
	
	
	previousLogMatrix = logProbabilityMatrix;//store previous matrix - settings are in pathHistory vector
	
	//	printLogMatrix();	
	
}



void TimingAnalyser::printIOIdata(){
	for (int i = 0;i < basicInterOnsetInterval.size();i++){
		printf("IOI[%i] : %i\n", i, basicInterOnsetInterval[i]);
	}
}


void TimingAnalyser::checkShiftMatrix(){

	if (moveMatrixToOptimalCostPosition){
		int phaseShiftIndex = (int)(meanPhaseIndex + movementFactor*(currentBestPhase - meanPhaseIndex));
		int tempoShiftIndex = (int)(meanTempoIndex + movementFactor*(currentBestTempo - meanTempoIndex));
		printf("current best tempo %i, best phase %i, mean tempo %i, shift of %i\n", currentBestTempo, currentBestPhase, meanTempoIndex, tempoShiftIndex);
		shiftMatrixToMatchBestPhase(phaseShiftIndex);
		shiftMatrixToMatchBestTempo(tempoShiftIndex);
	}
	
}



void TimingAnalyser::shiftMatrixToMatchBestPhase(const int& bestPhaseIndex){

//	int bestPhaseIndex = pathHistory[pathHistory.size()-1].bestPhase;
	printf("SHIFT:: BEST PHASE %i, ", bestPhaseIndex);
	
	int shift = bestPhaseIndex - meanPhaseIndex;
	DoubleMatrix tmpMat = logProbabilityMatrix;
	
	for (int t = 0;t < tempoRange;t++){
		for (int p = 0;p < phaseRange;p++){
			int switchIndex = p + shift;
			if (switchIndex < phaseRange && switchIndex >= 0){
				logProbabilityMatrix[t][p] = tmpMat[t][switchIndex];
			}else{
				logProbabilityMatrix[t][p] = INFINITY;
			}
		}
	}
	globalTimeOffset += shift * phaseScalar;
	printf("new phase shift of %i and offset is %i \n", shift, globalTimeOffset);
	recentPhaseShift = shift;
}

void TimingAnalyser::shiftMatrixToMatchBestTempo(const int& bestTempoIndex){
	
	//	int bestPhaseIndex = pathHistory[pathHistory.size()-1].bestPhase;
	printf("SHIFT:: BEST TEMPO %i, ", bestTempoIndex);
	
	int shift = bestTempoIndex - meanTempoIndex;
	DoubleMatrix tmpMat = logProbabilityMatrix;
	
	for (int t = 0;t < tempoRange;t++){
		for (int p = 0;p < phaseRange;p++){
			int switchIndex = t + shift;
			if (switchIndex < tempoRange && switchIndex >= 0){
				logProbabilityMatrix[t][p] = tmpMat[switchIndex][p];
			}else{
				logProbabilityMatrix[t][p] = INFINITY;
			}
		}
	}
	
	tempoMinimumOffset += shift * tempoScalar;
	printf("new tempo shift of %i and tempo offset is %i \n", shift, tempoMinimumOffset);
	recentTempoShift = shift;
}


bool TimingAnalyser::checkPhaseRange(const int& phaseToCheck){

	if (phaseToCheck < phaseRange &&  phaseToCheck >= 0)
		return true;
	else 
		return false;
}

bool TimingAnalyser::checkTempoRange(const int& tempoToCheck){
	
	if (tempoToCheck < tempoRange &&  tempoToCheck >= 0)
		return true;
	else 
		return false;
}

double TimingAnalyser::phaseHopCost(const int& phaseHop){
	//printf("phase hop %i returns %f", phaseHop, phaseCost*abs(phaseHop));
	return phaseCost*abs(phaseHop);
}

double TimingAnalyser::tempoHopCost(const int& tempoHop){
	return tempoCost*abs(tempoHop);
}

double TimingAnalyser::getTempo(const int& tempoIndex){
	return (tempoMinimumOffset + tempoIndex*tempoScalar);
}

double TimingAnalyser::getPhase(const int& phaseIndex){
	return globalTimeOffset + (phaseIndex*phaseScalar);
}

double TimingAnalyser::getTempo(const int& tempoIndex, const int& tempoOffset){
	return ((tempoOffset + tempoIndex)*tempoScalar);
}

double TimingAnalyser::getPhase(const int& phaseIndex, const int& phaseOffset){
	//phase offset for middle of matrix
	return (phaseOffset + phaseIndex)*phaseScalar;
}

int TimingAnalyser::getTempoIndex(const double& tempo){
	double index = tempo;
	index -= tempoMinimumOffset;
	index /= tempoScalar;
	return (int)round(index);
}


int TimingAnalyser::getPhaseIndex(const double& location){
	double index = location;
	index -= globalTimeOffset;
	index -= phaseMinimumOffset;
	index /= phaseScalar;
	return (int)round(index);
}


double TimingAnalyser::getMinimumTransitionCost(const int& interval, const int& tempoIndex, const int& phaseIndex){
	
	double minimumCost = 100000000.0;//very big
	
	//need to get the location of where we are now
	int endLocation = 0;//getBeatLocation()
	
//	int zeroLocation = getBeatLocation();	
//	getLocation(0, 0, beatPosition + interval);
	int tempo = 0;
	
//	for (int tempo = 0;tempo < logProbabilityMatrix.size();tempo++){
		for (int phase = 0;phase < logProbabilityMatrix[tempo].size();phase++){
			
			//double tempoTransition = getTempoTransitionCost(tempoIndex - tempo);
			//double phaseTransition = getPhaseTransitionCost(phaseindex - phase);
			
			
			
		}
//	}
			
		return minimumCost;
}


double TimingAnalyser::getMaximum(){
	double maxVal = 1.0;
	for (int tempo = 0;tempo < logProbabilityMatrix.size();tempo++){
		for (int phase = 0;phase < logProbabilityMatrix[tempo].size();phase++){
		if (logProbabilityMatrix[tempo][phase] > maxVal && logProbabilityMatrix[tempo][phase] < INFINITY)
			maxVal = logProbabilityMatrix[tempo][phase];
		}
	}
	return maxVal;
}

void TimingAnalyser::drawCostMatrix(){
	double maximum = getMaximum();
	double width = ofGetWidth() / phaseRange;
	double height = ofGetHeight() / tempoRange;
	
	for (int tempo = 0;tempo < logProbabilityMatrix.size();tempo++){
		for (int phase = 0;phase < logProbabilityMatrix[tempo].size();phase++){
			
			ofSetColor(0,0, 255.0*(maximum - logProbabilityMatrix[tempo][phase])/maximum );
			ofRect(phase*width, tempo*height, width, height);
			ofSetColor(255,255,255);
			
			string infoStr = "T "+ofToString(getTempo(tempo, tempoMinimumOffset));
			infoStr += "\nP "+ofToString(phase)+" ("+ofToString(getPhase(phase, globalTimeOffset))+")";
			infoStr += "\n"+ofToString(logProbabilityMatrix[tempo][phase]);
			
		//	ofDrawBitmapString(infoStr, phase*width + 10, tempo*height + 10);
		}
	}
	
}


void TimingAnalyser::setBestTempoAndPhase(Path& newPath){

	double bestCost = newPath.currentRoute[currentBestTempo][currentBestPhase].totalCost;// INFINITY;
	for (int t = 0;t < tempoRange;t++){
		for (int p = 0;p < phaseRange;p++){
			if (newPath.currentRoute[t][p].totalCost < bestCost){
				bestCost = newPath.currentRoute[t][p].totalCost;
				currentBestPhase = p;
				currentBestTempo = t;
			}
		}
	}
	
	newPath.bestTempo = currentBestTempo;// bestTempo;
	newPath.bestPhase = currentBestPhase;// bestPhase;
	
	printf("BEST TEMPO %i, PHASE %i :: ", currentBestTempo, currentBestPhase);
	string tString = "Tempo "+ofToString(getTempo(currentBestTempo, tempoMinimumOffset));
	tString += ", Phase "+ofToString(getPhase(currentBestPhase, globalTimeOffset));
	tString += ", BEST COST "+ofToString(bestCost);
	//tstring += ", error "+ofToString()
	printf("i.e. %s", tString.c_str());
	
}	

void TimingAnalyser::printPathHistory(){
/*	printf("tempo %i (%3.0f), phase %i (%3.0f)\n", pathHistory[i].bestTempo, 
		   getTempo(pathHistory[i].bestTempo, pathHistory[i].tempoOffset),
		   pathHistory[i].bestPhase, 
		   getPhase(pathHistory[i].bestPhase, pathHistory[i].globalPhaseOffset)
		   );
	*/
	int tempo, phase;
	
	RouteMatrix* tmpRoute;
	int pathIndex = pathHistory.size()-1;
	if ( pathHistory.size() > 0){
		pathHistory[0].phaseShiftApplied = 0;//hack to fix non zero
		tempo = pathHistory[pathHistory.size()-1].bestTempo;
		phase = pathHistory[pathHistory.size()-1].bestPhase;

		printf("PATH BEST[%i] tempo %i (%3.0f), phase %i (%3.0f) at COST %3.3f preCOST %3.3f addCOST %3.3f\n", 
		pathIndex,
		tempo, getTempo(tempo, pathHistory[pathHistory.size()-1].tempoOffset),
		phase, getPhase(phase, pathHistory[pathHistory.size()-1].phaseOffset),
		pathHistory[pathHistory.size()-1].currentRoute[tempo][phase].totalCost,
			   pathHistory[pathHistory.size()-1].currentRoute[tempo][phase].previousCost,
			   pathHistory[pathHistory.size()-1].currentRoute[tempo][phase].addedCost
			 );
	}//end if 
	
	
	while (pathIndex >= 0 && pathIndex < pathHistory.size()){
	//	printf("index %i, history size %i\n", pathIndex, pathHistory.size());
		tmpRoute = &(pathHistory[pathIndex].currentRoute);
		int previousTempo = (*tmpRoute)[tempo][phase].previousTempoIndex;
		int previousPhase = (*tmpRoute)[tempo][phase].previousPhaseIndex;// + pathHistory[pathIndex].phaseShiftApplied;
		//above was error on 24/12 that created wrong results - it is already included!
		
		printf("PTH[%i](%i,%i)\total{%3.2f}\t(prevCost %3.2f (%3.2f, %3.2f, %3.2f)\t:: adddedCost %3.2f\t: hop (%i, %i)\tprevious (%i, %i) applied (%i, %i){%i, %i}, ", 
			   pathIndex,  
			   tempo, phase, 
			   (*tmpRoute)[tempo][phase].totalCost, 
			   (*tmpRoute)[tempo][phase].previousCost,
			   (*tmpRoute)[tempo][phase].preHopCost,
			   (*tmpRoute)[tempo][phase].tempoHopCost,
			   (*tmpRoute)[tempo][phase].phaseHopCost,			   
			    (*tmpRoute)[tempo][phase].addedCost,
			    (*tmpRoute)[tempo][phase].tempoHop, 
			   (*tmpRoute)[tempo][phase].phaseHop,
			   previousTempo, previousPhase,  
			   pathHistory[pathIndex].tempoShiftApplied, pathHistory[pathIndex].phaseShiftApplied,
			  pathHistory[pathIndex].tempoOffset,  pathHistory[pathIndex].phaseOffset	   
			   
			   );
		
		if (pathIndex > 0){
			printf("Tempo %3.0f Phase %3.0f (PRED %3.0f) @event %i (%3.0f)\n", 
				   getTempo(previousTempo, pathHistory[pathIndex-1].tempoOffset),
				   getPhase(previousPhase, pathHistory[pathIndex-1].phaseOffset),
				   getPhase(previousPhase, pathHistory[pathIndex-1].phaseOffset) + getTempo(previousTempo, pathHistory[pathIndex-1].tempoOffset),
				   pathHistory[pathIndex].eventTimeObserved,
				   getPhase(phase, pathHistory[pathIndex].phaseOffset)
				   );
		}

		tempo =	 previousTempo;
		phase =  previousPhase;
		
		pathIndex--;
	}
	

}


void TimingAnalyser::processPathHistory(){

	printf("Process path history...\n");
	
	timingData.clear();
	
	maximumTempo = 0;
	minimumTempo = INFINITY;
	
	minimumPhaseHop = 0;
	maximumPhaseHop = 0;
	
	printf("\n");
	int tempo, phase;
	
	RouteMatrix* tmpRoute;
	int pathIndex = pathHistory.size()-1;
	int observedPhasePoint;
	int observedTempo, observedEventTime;
	int storedTempoHop = 0;
	int storedTempo = 0;
	
	if ( pathHistory.size() > 0){
		pathHistory[0].phaseShiftApplied = 0;//hack to fix non zero
		tempo = pathHistory[pathIndex].bestTempo;
		phase = pathHistory[pathIndex].bestPhase;
		observedPhasePoint = getPhase(phase, pathHistory[pathIndex].phaseOffset);
		
	}
		
	while (pathIndex >= 0){
		tmpRoute = &(pathHistory[pathIndex].currentRoute);
		int previousTempo = (*tmpRoute)[tempo][phase].previousTempoIndex;
		int previousPhase = (*tmpRoute)[tempo][phase].previousPhaseIndex;		
		IntVector v;
	
		observedTempo = storedTempo;
		observedPhasePoint = (int)round(getPhase(phase, pathHistory[pathIndex].phaseOffset));
		observedEventTime =  pathHistory[pathIndex].eventTimeObserved;
		
		if (observedTempo > maximumTempo)
			maximumTempo = observedTempo;
		
		if (observedTempo < minimumTempo && observedTempo > 0)
			minimumTempo = observedTempo;
		
		int phaseHop = (*tmpRoute)[tempo][phase].phaseHop;
		
		if (phaseHop < minimumPhaseHop)
			minimumPhaseHop = phaseHop;
		
		if (phaseHop > maximumPhaseHop)
			maximumPhaseHop = phaseHop;
		
		v.push_back(observedTempo);//timingdata[][0]
		v.push_back(observedPhasePoint);//timingdata[][1]
		v.push_back(observedEventTime);//timingdata[][2]
		v.push_back(storedTempoHop);
		v.push_back( (*tmpRoute)[tempo][phase].phaseHop);//timingData[][4]
		v.push_back(observedEventTime - observedPhasePoint);//timingdata[][5]
		
		timingData.push_back(v);
		
		
		printf("Index %i TEMPO %i, PHASE %i, EVENT %i HOP (%i, %i) error %i \n", 
			   (int) timingData.size(),
			   observedTempo, observedPhasePoint, 
			   observedEventTime, storedTempoHop,
			   (*tmpRoute)[tempo][phase].phaseHop, 
			   (observedEventTime - observedPhasePoint));
		//checking this ned to align tempo, observed phase and observed events
		
		//N.B. MUST STORE THIS FROM FUTURE POINT
		storedTempo = (int)round(getTempo(tempo, pathHistory[pathIndex].tempoOffset));
		storedTempoHop = (*tmpRoute)[tempo][phase].tempoHop;
		
		tempo =	 previousTempo;
		phase =  previousPhase;
		
		pathIndex--;
	}
	
	reverse(timingData.begin(), timingData.end());

	clearHistogram();
	
	for (int i = 0;i < timingData.size();i++){
		printf("Rev index %i tempo %i, phase hop %i, error %i, beat location %i\n", i, 
			   timingData[i][0], timingData[i][4], timingData[i][5], (int) beatPosition[i]);
		
		int error = (int) timingData[i][4] + timingData[i][5];
		
		if (abs(error) < 30 && (int)beatPosition[i] - 1 < timingHistogram.size())
			timingHistogram[(int)beatPosition[i] - 1].push_back(error);
	}
	
	printf("Maximum is %f and min %f hopMin %i hopMax %i\n", maximumTempo, minimumTempo, minimumPhaseHop, maximumPhaseHop);
	minDrawTempo = minimumTempo + minimumPhaseHop - 4;
	maxDrawTempo = maximumTempo + maximumPhaseHop + 4;
	
	getHistogramResults();
	
}

void TimingAnalyser::clearHistogram(){
	
	timingHistogram.clear();
	IntVector v;
	for (int i = 0;i < 4;i++){
		v.clear();
		timingHistogram.push_back(v);
	}
}


void TimingAnalyser::getHistogramResults(){
	DoubleVector results;
	for (int i = 0;i < 4;i++){
		double total;
		for (int k = 0; k < timingHistogram[i].size();k++){
			total += timingHistogram[i][k];
		}
		total /= timingHistogram[i].size();
		results.push_back(total);
		printf("Mean at beat %i is %.2f\n", i, results[i]);
	}
}

#pragma mark -drawTempoCurve
void TimingAnalyser::drawTempoCurve(){
	if (timingData.size() > 0){
		screenWidth = ofGetWidth();
		screenHeight = ofGetHeight();
		double height = screenHeight;
//		int	minDrawTempo = minimumTempo + minimumPhaseHop - 4;
//		int	maxDrawTempo = maximumTempo + maximumPhaseHop + 4;
		
		if (maximumTempo > 0)
			height /= (maxDrawTempo - minDrawTempo);
	
		double bpmMaxTempo = msToBpm(minDrawTempo);
		double bpmMinTempo = msToBpm(maxDrawTempo);
		double bpmHeight = screenHeight / (bpmMaxTempo - bpmMinTempo);
		
	//int minHeightPixels = minHeight * height;
		
		ofBackground(255);
		double width = ofGetWidth()/numberOfPointsPerPage;
	
		int lastHeightPoint = getHeightPoint(height * (timingData[0][0] - minDrawTempo) );
		int newHeightPoint;
		
		double error;
		int Xpoint;
		
		//horizonal lines for tempo
		ofSetColor(150);
		if(!drawBPM){
		int horizontalTempo = minDrawTempo + 20 - minDrawTempo%20;
	//	while (horizontalTempo > minDrawTempo)
	//		horizontalTempo -= 20;
			while (horizontalTempo < maxDrawTempo){
				ofSetColor(150);
				ofLine(0, getHeightPoint(height*(horizontalTempo - minDrawTempo)), screenWidth, getHeightPoint(height*(horizontalTempo - minDrawTempo)));
				ofSetColor(50);
			ofDrawBitmapString(ofToString(horizontalTempo), 20, getHeightPoint(height*(horizontalTempo - minDrawTempo)));
			//	timesFont.drawString(ofToString(horizontalTempo), 20, getHeightPoint(height*(horizontalTempo - minDrawTempo)));//FONT
				horizontalTempo += 20;
			}//end while
		
		} else{
			int resolution = 2;
			
				if (maxDrawTempo - minDrawTempo > 150)
					resolution = 4;
			if (maxDrawTempo - minDrawTempo > 400)
				resolution = 8;
			
			int horizontalTempo = (int)bpmMinTempo + resolution - ((int)bpmMinTempo % resolution);//mod ten above min
			while (horizontalTempo < bpmMaxTempo) {
				int tmp = getHeightPoint(bpmHeight * (horizontalTempo - bpmMinTempo));
				ofSetColor(150);
				ofLine(0, tmp , screenWidth, tmp);
				ofSetColor(50);
				ofDrawBitmapString(ofToString(horizontalTempo), 20, tmp+2);
				std::string number = ofToString(horizontalTempo);
//				timesFont.drawString(number, 20, tmp+2);//FONT
				horizontalTempo += resolution;
			}
			
		}

		//vertical bars 
		int solidBarEvery = 4;//normally every four beats
		
		int counter = 0;
		for (int i = 0;i < timingData.size();i++){
			ofSetColor(200);
			bool drawLine = false;

			
			//this useful for drawing on ISMIR paper
			if (mozartTriplets){
				//stronger lines on the bars
				solidBarEvery = 12;//- pollini triplets
				
				if ((i % solidBarEvery == 0)){// || beatPosition[i] == 0 ){
					ofSetColor(80);
					drawLine = true;
					counter = 0;
				}
				
				//normal every three (eighth notes)
				if (i % 3 == 0 && mozartTriplets){
					drawLine = true;
				}
			} 
			
			
			int beatsPerBar = 4;//set 3 for some beatles tracks
			int mainBeat = 0; 
			
			if ((int)i % beatsPerBar == mainBeat && !mozartTriplets){
				//was BeatPosoiton[i] - dodgy
				drawLine = true;
			}
			/*
			if (i == tempoVariationEndIndex){
				ofSetColor(200, 0,0);
				drawLine = true;
			}
			 */
			
			if (drawLine){
				ofLine((i-startPoint) * width, 0, (i-startPoint) * width, ofGetHeight());
			}
		}
		
		
		for (int i = startPoint+1;i < min((int)timingData.size(), startPoint + numberOfPointsPerPage);i++){
			if (!drawIOIintervals){
			ofSetColor(0,0,0);
			Xpoint = (i-startPoint) * width ;
			if (!drawBPM)
				newHeightPoint = getHeightPoint(height * (timingData[i][0] - minDrawTempo));
			else 
				newHeightPoint = getHeightPoint(bpmHeight * (msToBpm(timingData[i][0]) - bpmMinTempo));

			ofLine((i-startPoint-1) * width, lastHeightPoint,Xpoint, newHeightPoint);
			lastHeightPoint = newHeightPoint;
			
			error = timingData[i][2] - timingData[i][1];
			int phaseHop = timingData[i][4];
			
			if (!blackAndWhite)
				ofSetColor(255,0,0);
			else 
				ofSetColor(155);
			
			if (drawExpressiveTimingData)
			ofLine(Xpoint, newHeightPoint, Xpoint, newHeightPoint - (phaseHop * height));
			//line up from tempo to where we observe point
			
			if (!blackAndWhite)
				ofSetColor(0,0, 155);
			else 
				ofSetColor(60);
				
			if (drawExpressiveTimingData)	
			ofCircle(Xpoint, newHeightPoint - ((phaseHop + error) * height), 2);
//			ofLine((i-startPoint) * width , screenHeight/2, (i-startPoint) * width , (screenHeight/2)*(1+(timingData[i][3]/60.0)));
			
			//if (!blackAndWhite)	
			ofSetColor(0,0,255);
			ofLine(((playingIndex)-startPoint) * width, 0, ((playingIndex)-startPoint) * width, screenHeight);
			
			}//if not draw IOI
			
			if (drawIOIintervals){	
				ofSetColor(160,160,160);
				ofLine((i-startPoint-1) * width, getHeightPoint(height * (basicInterOnsetInterval[i-1]- minDrawTempo)), 
					   (i-startPoint) * width, getHeightPoint(height * (basicInterOnsetInterval[i]- minDrawTempo)) );
			}
	//		ofLine((i-startPoint) * width , screenHeight/2, (i-startPoint) * width , (screenHeight/2)*(1+(timingData[i][4]/60.0)));
			
		}//end for
	
	}
	
}

static const int moveAmount = 2;

void TimingAnalyser::widenDrawWindow(){
	maxDrawTempo += moveAmount;
	minDrawTempo -= moveAmount;
}

void TimingAnalyser::narrowDrawWindow(){
	maxDrawTempo -= moveAmount;
	minDrawTempo += moveAmount;
}

void TimingAnalyser::moveDrawWindowUp(){
	maxDrawTempo += moveAmount;
	minDrawTempo += moveAmount;
}

void TimingAnalyser::moveDrawWindowDown(){
	maxDrawTempo -= moveAmount;
	minDrawTempo -= moveAmount;
}

int TimingAnalyser::getIndexAtMouseXposition(const int& Xpos){
	int numberPointsIn = (int)round(Xpos * numberOfPointsPerPage / ofGetWidth());
	
	printf("number pts %i size %i , start %i no. in from left %i\n", 
		   numberOfPointsPerPage, (int)timingData.size(), 
		   startPoint, numberPointsIn);
	
	numberPointsIn += startPoint;
	
	printf("MOUSE X IS %i == %i\n", Xpos, numberPointsIn);
	
	return min((int)timingData.size(), numberPointsIn);
}

int TimingAnalyser::getHeightPoint(float f){
	return (int)round(screenHeight - f);
}



void TimingAnalyser::printMatrix(DoubleMatrix& logMatrix){
	for (int t = 0;t < logMatrix.size();t++){
	for (int p = 0;p < logMatrix[t].size();p++){
		printf("(%i,%i)= %2.1f ", t, p, logMatrix[t][p]);
	}
		printf("\n");
	}

}

void TimingAnalyser::printLogMatrix(){
	printf("\n");
	for (int p = 0;p < phaseRange;p++){
	printf("LOG M [%i] %2.2f, ", p, logProbabilityMatrix[0][p]);
	}
	printf("offset %i\n", globalTimeOffset);
}

void TimingAnalyser::exportTimingData(){
	ofstream ofs(timingDataFilepath.c_str());
//	ofs << "output timng data size " << (int) timingData.size() << "IOI size " << (int)basicInterOnsetInterval.size() << endl;
	
	for (int i = 0;i < min((int)timingData.size(), (int) basicInterOnsetInterval.size());i++){
		for (int k = 0; k < 5;k++){
			ofs << timingData[i][k] << "\t";
		}
		ofs << basicInterOnsetInterval[i] << "\t";
		ofs << endl;
	}

}

void TimingAnalyser::importTimingData(std::string importFileName){
	ifstream ifs(importFileName.c_str());
	printf("IMPORT FILE %s\n", importFileName.c_str());
	timingData.clear();
	string value;
	int count = 0;
	//Notation n;
	IntVector v;
	while ( ifs.good() )
	{
		getline ( ifs, value, '\n' ); // read a string until next enter command
		//		cout << "disp " <<  string( value, 0, value.length()-1 ) << endl; // display value removing the first and the last character from it
		
		//printf("size of line is %i\n", (int)value.size()); 
		v.clear();
		
		if (value.size() > 0){
			
			string::size_type start = value.find_first_not_of(" \t\v");
			string part = value.substr(start, string::npos);
			string otherPart = value;
			printf("input line%i : '%s'\n", count, otherPart.c_str());
			size_t found;
			found=part.find_first_of("\t\n");
			size_t startC = 0;
			int r = 0;
			while (found!=string::npos){
				string section = otherPart.substr(startC, found - startC);
				int  my_int = atoi(section.c_str());
				v.push_back(my_int);
				printf("timing[%i][%i] ", count, r, my_int);
				printf("section '%s'\n", section.c_str());
				
				startC= otherPart.find_first_not_of("\t\n",found);
				found = otherPart.find_first_of("\t\n",found+1);
				r++;
			}
			timingData.push_back(v);
		}//end if > 0
 
		count++;
	
	 }
	
	
	basicInterOnsetInterval.clear();
	
	maximumTempo = 0;
	minimumTempo = INFINITY;
	
	for (int i = 0;i < max(0, (int)timingData.size()-1);i++){
		printf("TEMPO %i %i IOI %i\n", timingData[i][0], timingData[i][1], timingData[i][4]);
		
		basicInterOnsetInterval.push_back(timingData[i][5]);
		
		int observedTempo = timingData[i][0];
	//	printf("observed %i\n", observedTempo);
		
		if (observedTempo > maximumTempo)
			maximumTempo = observedTempo;
		
		if (observedTempo < minimumTempo && observedTempo > 0)
			minimumTempo = observedTempo;
		
	}
	
	
	printf("min t %f, Max t %f\n", minimumTempo, maximumTempo);
}


void TimingAnalyser::exportProcessedBeatTimes(const double& firstBeatTime){
	ofstream ofs(processedBeatTimesFilepath.c_str());
	
	for (int i = 0;i < min((int)timingData.size(), (int) basicInterOnsetInterval.size());i++){
		double beatTime = (timingData[i][1] + firstBeatTime)/1000.0;
		printf("Beat %i : %f\n", i, beatTime);
		ofs << beatTime;
		ofs << endl;
	}
}


double TimingAnalyser::msToBpm(const double& ms){
	return 60000./ms;
}


void TimingAnalyser::calculateTempoLimits(){
	
	maximumTempo = 0;
	minimumTempo = INFINITY;
	
	for (int i = 0;i < max(0, (int)timingData.size()-1);i++){
		
		int observedTempo = timingData[i][0];
		
		if (observedTempo > maximumTempo)
			maximumTempo = observedTempo;
		
		if (observedTempo < minimumTempo && observedTempo > 0)
			minimumTempo = observedTempo;
		
	}
	
	tempoVariationStartIndex = min(4, (int)timingData.size());
	tempoVariationEndIndex = max(0, (int)timingData.size() - 4);
	printf("VARIATION TO BE CALCULATED BETWEEN %i and %i\n", tempoVariationStartIndex, tempoVariationEndIndex);
}

double TimingAnalyser::calculateTempoVariation(){
	if (timingData.size() > tempoVariationStartIndex && tempoVariationEndIndex + 1< timingData.size()){
		double mean = 0;
		
		for (int i = tempoVariationStartIndex;i <= tempoVariationEndIndex;i++){
			mean += timingData[i][0];
		}
		mean /= (tempoVariationEndIndex - tempoVariationStartIndex);
		double variation = 0;
		double fluctuation = 0;
		int lastValue = timingData[tempoVariationStartIndex][0];
		
		for (int i = tempoVariationStartIndex;i <= tempoVariationEndIndex;i++){
			variation += (timingData[i][0] - mean)*(timingData[i][0] - mean);
			fluctuation += abs(timingData[i][0] - lastValue);
			lastValue = timingData[i][0];
		}
		variation /= (tempoVariationEndIndex - tempoVariationStartIndex);
		variation = sqrt(variation);
		fluctuation /= (tempoVariationEndIndex - tempoVariationStartIndex);
		
		double slope = (double) timingData[tempoVariationStartIndex][0] / timingData[tempoVariationEndIndex][0];
		printf("TEMPO STATS:\nMean: %.2f\nVariance %.2f \nSlope %.1f per cent\nAverage delta Tempo %.2f", 
			   mean, variation, (slope*100.0), fluctuation);
	}

}


void TimingAnalyser::zoomIn(){
	numberOfPointsPerPage /= 2;
}


void TimingAnalyser::zoomOut(){
	numberOfPointsPerPage *= 2;
}


/*
 double TimingAnalyser::getCost(const int& eventTime, const int& eventLocation, const int& tempoIndex, const int& phaseIndex){
 
 double interval = eventLocation - lastBeatPosition;
 double newLocation = getPhase(phaseIndex) + interval*getTempo(tempoIndex);	
 double phaseCost = fabs(newLocation - eventTime);
 
 return phaseCost;
 }
 */