view src/DynamicBayesianArray.cpp @ 11:df9c838d0b8f

fixed bug where the full midi file was sometimes not drawn (since minNoteIndexToPrint) was at the max and did not decrease properly in draw routine

/*
 *  DynamicDynamicBayesianArray.cpp
 *  midiCannamReader
 *
 *  Created by Andrew on 17/07/2011.
 *  Copyright 2011 QMUL. All rights reserved.
 *
 */

#include "DynamicBayesianArray.h"
#include "math.h"
#include "ofMain.h"

DynamicBayesianArray::DynamicBayesianArray(){
	
//	prior.createVector(240);
//	likelihood.createVector(240);
//	posterior.createVector(240);	
	testVector.createVector(240);
	testVector.addGaussianShape(100,10, 0.1);
	
	
	likelihoodNoise = 0.5;
	likelihoodMean = ARRAY_SIZE/2;
	likelihoodStdDev = ARRAY_SIZE / 12;
	initialiseArray();
}

void DynamicBayesianArray::initialiseArray(){
	
	//maximumIndex = 12;//change this
	setGaussianPrior(ARRAY_SIZE/2, ARRAY_SIZE/1);	
	setGaussianLikelihood(ARRAY_SIZE/2, ARRAY_SIZE/1);//likelihoodMean, likelihoodStdDev);
	
	calculatePosterior();
	renormalisePosterior();
	posteriorDecayRate = 0.06;
	
	eighthNoteProportion = 0.35;//must be less than 0.5 to discriminate - was 0.4
	earlySixteenthNoteProportion = 0;
	lateSixteenthNoteProportion = 0;
	decayNoiseAmount = 0.1;
	decayNoiseStdDev = ARRAY_SIZE/24;
	standardDeviation = likelihoodStdDev;
	setDecayNoiseGaussian(ARRAY_SIZE/2, decayNoiseStdDev);
	
	setGaussianLikelihood(likelihoodMean, likelihoodStdDev);
}


void DynamicBayesianArray::setGaussianPrior(float mean, float StdDev){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		prior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
		//posterior[i] = prior[i];
	}
}

void DynamicBayesianArray::setGaussianPosterior(float mean, float StdDev){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		posterior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	}
}


void DynamicBayesianArray::setGaussianLikelihood(float mean, float StdDev){
	if (mean >= 0 && mean <= ARRAY_SIZE){
		int i;	float eighthDifference;
		int eighthPosition = ((int)mean + ARRAY_SIZE/2)%ARRAY_SIZE;
		float mainDifference;
		float gaussianProportion = 1 - likelihoodNoise;
		
		for (i=0;i < ARRAY_SIZE;i++){
			
			mainDifference = min( double(fabs(i-mean)) , (double)(i + ARRAY_SIZE - mean));
			//without * (1 - eighthNoteProportion)
			likelihood[i] = gaussianProportion  * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
			
			likelihood[i] += (likelihoodNoise / ARRAY_SIZE);
			//likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
			//(double) (likelihoodNoise / ARRAY_SIZE) );
		}
		//	renormaliseArray(&likelihood[0], ARRAY_SIZE);
	}//end if mean within limits
}

void DynamicBayesianArray::calculatePosterior(){
	int i;
	for (i=0;i < ARRAY_SIZE;i++){
		posterior[i] = likelihood[i] * prior[i];
	}
	//renormalisePosterior();
}


float DynamicBayesianArray::getMaximum(float *ptr, int length){
	int i;
	float max = 0;
	for (i=0;i < length;i++){	
		if (*(ptr+i)>max)
			max = *(ptr+i);
	}
	maximumValue = max;
	return max;
}

float* DynamicBayesianArray::getMaximumEstimate(float *ptr, int length){
	float returnArray[2];
	int i;
	float max = 0;
	maximumIndex = 0;
	for (i=0;i < length;i++){	
		if (*(ptr+i)>max){
			max = *(ptr+i);
			maximumIndex = i;
		}
	}
	returnArray[0] = max;
	returnArray[1] = maximumIndex;
	maximumValue = max;
	return &returnArray[0];
}



double DynamicBayesianArray::getIntegratedEstimateIndex(){
	int i;
	float integratedQuantity = 0;
	float integratedTotal = 0;
	double integratedIndex = 0;
	for (i=0;i < ARRAY_SIZE;i++){	
		integratedQuantity += posterior[i];//the values of the probability distribution
		integratedTotal += i*posterior[i];
	}
	if (integratedQuantity > 0){
		integratedIndex = integratedTotal / integratedQuantity;
	}
	integratedEstimate = (float) integratedIndex;
	return integratedIndex;
}


double DynamicBayesianArray::calculateStandardDeviation(){
	
	double total = 0;
	double pdfSum;
	double variance = 0;
	for (int i=0;i < ARRAY_SIZE;i++){	
		//*posterior[i] * 		
		total += posterior[i] * (i - integratedEstimate) * (i - integratedEstimate);//the values of the probability distribution
		pdfSum += posterior[i];
	}
	
	if (pdfSum > 0)
		variance = total / pdfSum;
	else 
		variance = ARRAY_SIZE;
	
	standardDeviation = sqrt(variance);
	return standardDeviation;
}



void DynamicBayesianArray::renormaliseArray(float *ptr, int length){
	int i;
	float totalArea = 0;
	for (i=0;i < length;i++){
		totalArea += *(ptr+i);
	}
	
	for (i=0;i < length;i++){
		*(ptr+i) /= totalArea;
	}
	
}

void DynamicBayesianArray::resetPrior(){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		prior[i] = posterior[i];
	}
}

void DynamicBayesianArray::renormalisePosterior(){
	int i;
	float totalArea = 0;
	for (i=0;i < ARRAY_SIZE;i++){
		totalArea += posterior[i];
	}
	
	for (i=0;i < ARRAY_SIZE;i++){
		posterior[i] /= totalArea;
	}
	
}

void DynamicBayesianArray::decayPosterior(){
	float *pointer;
	pointer = getMaximumEstimate(&posterior[0], ARRAY_SIZE);	
	float maximum;
	maximum = *pointer;
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		posterior[i] += (maximum - posterior[i]) * posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
	}
	maximumIndex = *(pointer+1);
}

void DynamicBayesianArray::setDecayNoiseGaussian(float mean, float StdDev){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		decayNoiseArray[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	}
}

void DynamicBayesianArray::decayPosteriorWithGaussianNoise(){
	
	int i;
	float currentMaximum = getMaximum(&posterior[0], ARRAY_SIZE);
	for (i=0;i<ARRAY_SIZE;i++){
		posterior[i] += decayNoiseArray[(i - (int)maximumIndex + ((3*ARRAY_SIZE)/2)) % ARRAY_SIZE] * currentMaximum * decayNoiseAmount;
		//posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
	}
	
}

void DynamicBayesianArray::resetMaximumPosterior(){
	int i;
	float max = 0;
	for (i=0;i < ARRAY_SIZE;i++){	
		if (posterior[i]>max){
			maximumIndex = i;
			max = posterior[i];
		}
	}
}

void DynamicBayesianArray::translateDistribution(int translationIndex){
	int tmpIndex;
	//copy array
	int i;
	for (i=0;i < ARRAY_SIZE;i++){
		tempPosteriorArray[i] =	posterior[i] ;
	}
	//translate values
	for (i=0;i < ARRAY_SIZE;i++){
		tmpIndex = (i + translationIndex + ARRAY_SIZE)%ARRAY_SIZE;
		posterior[tmpIndex] = tempPosteriorArray[i]; 
	}
	//now delete tmp array
}


void DynamicBayesianArray::drawFloatArray(float* arrayToDraw, const int& minIndex, const int& maxIndex){
	
	if (minIndex >= 0){
	
		double stepSize = ofGetWidth() / (double)(maxIndex - minIndex);
		double screenHeight = ofGetHeight();
		double maxVal = getMaximum(&arrayToDraw[0], maxIndex);

		for (int i = minIndex+1;i < maxIndex;i++){
		
			ofLine (stepSize*(i-1), screenHeight * (1 - arrayToDraw[i-1] / maxVal), stepSize*i, screenHeight * (1 - arrayToDraw[i] / maxVal) );
		}
	
	}
		
		
}

 

/*
void DynamicBayesianArray::drawDoubleArray(double[]& arrayToDraw, const int& minIndex, const int& maxIndex){

	if (minIndex >= 0 && maxIndex <= arrayToDraw.size(0))

}
*/
/*
 void DynamicBayesianArray::setGaussianLikelihoodForBeats(float mean, float StdDev){
 //this has eighth and sixteenth positions included
 
 if (mean >= 0 && mean <= ARRAY_SIZE){
 int i;	float eighthDifference;
 int eighthPosition = ((int)mean + ARRAY_SIZE/2)%ARRAY_SIZE;
 int earlySixteenthPosition = ((int)mean + (3*ARRAY_SIZE/4))%ARRAY_SIZE;;
 int lateSixteenthPosition = ((int)mean + (ARRAY_SIZE/4))%ARRAY_SIZE;;
 
 float mainDifference, sixteenthDifference;
 float gaussianProportion = 1 - likelihoodNoise;
 float mainProportion = (1 - eighthNoteProportion - earlySixteenthNoteProportion - lateSixteenthNoteProportion);
 
 for (i=0;i < ARRAY_SIZE;i++){
 
 mainDifference = min( fabs(i-mean) , (double)(i + ARRAY_SIZE - mean));
 likelihood[i] = gaussianProportion * mainProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
 
 eighthDifference = min( abs(i - eighthPosition) , i + ARRAY_SIZE - eighthPosition);
 eighthDifference = min(eighthDifference , (float)(ARRAY_SIZE + eighthPosition - i ));
 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
 likelihood[i] += gaussianProportion * eighthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(eighthDifference)*(eighthDifference)/(2*StdDev*StdDev)) ;
 
 sixteenthDifference = min( abs(i - earlySixteenthPosition) , i + ARRAY_SIZE - earlySixteenthPosition);
 sixteenthDifference = min(sixteenthDifference , (float)(ARRAY_SIZE + earlySixteenthPosition - i ));
 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
 likelihood[i] += gaussianProportion * earlySixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
 
 sixteenthDifference = min( abs(i - lateSixteenthPosition) , i + ARRAY_SIZE - lateSixteenthPosition);
 sixteenthDifference = min(sixteenthDifference , (float)(ARRAY_SIZE + lateSixteenthPosition - i ));
 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
 likelihood[i] += gaussianProportion * lateSixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
 
 
 
 likelihood[i] += (likelihoodNoise / ARRAY_SIZE);
 //likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
 //(double) (likelihoodNoise / ARRAY_SIZE) );
 }
 //	renormaliseArray(&likelihood[0], ARRAY_SIZE);
 }//end if mean within limits
 }
 */