comparison DrumTimingLoader_OF/PerformanceAnalyserSrc/TimingAnalyser.cpp @ 1:106bc2d4f702

added timing analyser 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;
+	
+	tempoScalar = 2.0;
+	phaseScalar = 2.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::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
+	
+	}
+	
+}
+
+
+void TimingAnalyser::widenDrawWindow(){
+	maxDrawTempo += 2;
+	minDrawTempo -= 2;
+}
+
+
+
+
+void TimingAnalyser::narrowDrawWindow(){
+	maxDrawTempo -= 2;
+	minDrawTempo += 2;
+}
+
+void TimingAnalyser::moveDrawWindowUp(){
+	maxDrawTempo += 2;
+	minDrawTempo += 2;
+}
+
+void TimingAnalyser::moveDrawWindowDown(){
+	maxDrawTempo -= 2;
+	minDrawTempo -= 2;
+}
+
+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;
+ }
+ */