annotate bayesianArraySrc/DynamicVector.cpp @ 4:45b5cf9be377

checking through Bayesian update procedure - working without cross update method.

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


#include "DynamicVector.h"

DynamicVector::DynamicVector(){
	length = 0;
	arraySize = 0;
	maximumValue = 0;
	MAPestimate = 0;
	offset = 0;
	scalar = 1;
	integratedEstimate = length/2;
	
	gaussianLookupMean = (double) GAUSSIAN_LOOKUP_LENGTH/2;
	gaussianLookupStdDev = (double)(GAUSSIAN_LOOKUP_LENGTH/16);
	double factor = 1.0;//(1.0 / (gaussianLookupStdDev*sqrt(2*PI)) );//1.0;//-1.0/(2*PI*sqrt(gaussianLookupStdDev));
	for (int i = 0;i < GAUSSIAN_LOOKUP_LENGTH;i++){
	gaussianLookupTable[i] = factor*exp(-1.0*(i-gaussianLookupMean)*(i-gaussianLookupMean)/(2.0*gaussianLookupStdDev*gaussianLookupStdDev));
	}
	
}

void DynamicVector::copyFromDynamicVector(const DynamicVector& dynamicVec){

	if (dynamicVec.length == length){
		offset = dynamicVec.offset;
		for (int i = 0;i < length;i++)
			array[i] = dynamicVec.array[i];
	}
	else{
		printf("CANNOT COPY VECTORS OF NON SAME LENGTH!!\n");
	}
}

void DynamicVector::createVector(int len){
	array.clear();
	for (int i = 0; i < len;i++){
		array.push_back(0);
	}
	length = len;
	arraySize = array.size();
	integratedEstimate = length/2;
}


double DynamicVector::getMaximum(){
	int i;
	double max = 0;
	for (i=0;i < length;i++){	
		if (array[i] > max){
			max = array[i];
			MAPestimate = i;
		}
	}
	maximumValue = max;
	return max;
}

int DynamicVector::getMAPestimate(){
	int i;
	double max = 0;
	for (i=0;i < length;i++){	
		if (array[i] > max){
			max = array[i];
			MAPestimate = i;
		}
	}
	maximumValue = max;
	return MAPestimate;
}

double DynamicVector::getIntegratedEstimate(){
	//returns the index of the integrated average - where the probability distribution is centred
	integratedEstimate = 0;
	double integratedTotal = 0;
	for (int i = 0;i < length;i++){
		integratedEstimate += array[i]*i;
		integratedTotal += array[i];
	}
	if (integratedTotal > 0){
		integratedEstimate /= integratedTotal;
	}
	return integratedEstimate;
}

void DynamicVector::updateIntegratedEstimate(){
	//returns the index of the integrated average - where the probability distribution is centred
	integratedEstimate = 0;
	double integratedTotal = 0;
	for (int i = 0;i < length;i++){
		integratedEstimate += array[i]*i;
		integratedTotal += array[i];
	}
	if (integratedTotal > 0){
		integratedEstimate /= integratedTotal;
	}

}

void DynamicVector::updateLimitedIntegratedEstimate(){
	//returns the index of the integrated average - where the probability distribution is centred
	//but limited round the MAP estimate
	double tmp = getMaximum();
	int limit = min(MAPestimate, length - MAPestimate);
	int start = max(0, MAPestimate - limit);
	int end = min(MAPestimate + limit, length-1);
	
	integratedEstimate = 0;
	double integratedTotal = 0;
	for (int i = start;i <= end;i++){
		integratedEstimate += array[i]*i;
		integratedTotal += array[i];
	}
	if (integratedTotal > 0){
		integratedEstimate /= integratedTotal;
	}
	
}


void DynamicVector::zero(){
	for (int i = 0;i < array.size();i++)
		array[i] = 0;
}

void DynamicVector::renormalise(){
	double tmpMax = getMaximum();
	if (tmpMax > 0){
//	printf("renormalise : max is %f and size is %i\n", tmpMax, arraySize);		
		for (int i = 0;i < array.size();i++)
			array[i] /= tmpMax;

	}
	//printArray();
}

void DynamicVector::doProduct(DynamicVector& arrayOne, DynamicVector& arrayTwo){

	for (int i = 0;i < arrayOne.length;i++)
		array[i] = arrayOne.array[i] * arrayTwo.array[i];
}


void DynamicVector::printArray(){
	for (int i = 0;i < arraySize;i++){
		printf("[%i] = %f\n", i, array[i]);
	}
}

void DynamicVector::translateDistribution(int translationIndex){
	int tmpIndex;
	DoubleVector tmpArray;
	int i;

	for (i=0;i < arraySize;i++){
		tmpArray.push_back(array[i]);
	}
	//translate values
	for (i=0;i < arraySize;i++){
		tmpIndex = (i + translationIndex + arraySize)%arraySize;
		array[tmpIndex] = tmpArray[i]; 
	}
	tmpArray.clear();
	//now delete tmp array
}



void DynamicVector::addGaussianShapeFromRealTime(const double& actualTime, const double& StdDev, double factor){
	
	double mean = getRealTermsAsIndex(actualTime);
//	printf("Gaussian realtime %f at index %f\n", actualTime, mean);
	int i;
	double std_dev_factor = (2*StdDev*StdDev);
	factor *= (1/(StdDev*sqrt(2*PI)));
	int maxVal = min((int) array.size(), (int)(mean + 4.8*StdDev));
	int minVal = max(0, (int)(mean - 4.8*StdDev));
	
	for (i=minVal;i < maxVal;i++){
		array[i] += factor*exp(-1*(i-mean)*(i-mean)/(std_dev_factor));
	}
	
	//	addGaussianShapeByLookupTable(mean, StdDev, factor);
}


void DynamicVector::addGaussianShape(const double& mean, const double& StdDev, double factor){
	
	int i;
	double std_dev_factor = (2*StdDev*StdDev);
	factor *= (1/(StdDev*sqrt(2*PI)));
	int maxVal = min((int) array.size(), (int)(mean + 4.8*StdDev));
	int minVal = max(0, (int)(mean - 4.8*StdDev));
	
	for (i=minVal;i < maxVal;i++){
		array[i] += factor*exp(-1*(i-mean)*(i-mean)/(std_dev_factor));
	}
	
//	addGaussianShapeByLookupTable(mean, StdDev, factor);
}

void DynamicVector::addGaussianShapeByLookupTable(double& mean, double& StdDev, double factor){
	int i;
	int lookupIndex ;
	factor *= (1/(StdDev*sqrt(2*PI)));
	for (i=0;i<array.size()-1;i++){
		lookupIndex = round(getLookupIndex(i, mean, StdDev));
		array[i] += factor*gaussianLookupTable[lookupIndex];
	}
	//printf("ADDED GAUSSIAN SHAPE %i\n", (int)array.size());
}

double DynamicVector::getLookupIndex(const int& i, const double& mean, const double& StdDev){
	
	double Z = ((double)i - mean)/StdDev;
	double lookupIndex = Z*gaussianLookupStdDev + gaussianLookupMean;
	
	if (lookupIndex < 0)
		lookupIndex = 0;
	
	if (lookupIndex >= GAUSSIAN_LOOKUP_LENGTH)
		lookupIndex = GAUSSIAN_LOOKUP_LENGTH-1;
	
//	(i - mean)*(i-mean)*(GAUSSIAN_LOOKUP_LENGTH*GAUSSIAN_LOOKUP_LENGTH/16.0)/(StdDev*StdDev);
	return lookupIndex;
}

void DynamicVector::addTriangularShape(double mean, double width, double factor){
	int i;

	for (i= max(0., (double)(mean - width));i < min((mean+width), (double)array.size());i++){
		array[i] += factor * abs(i - mean) / mean;		
	}

}

void DynamicVector::addConstant(const double& value){
	for (int i=0;i<array.size();i++){
		array[i] += value;
	}
}


void DynamicVector::addToIndex(const int& index, const double& constant){
	array[index] += constant;
}


double DynamicVector::getIndexInRealTerms(const int& index){
	if (index < arraySize)
		return (offset + scalar*index);
	else
		return 0;
}

double DynamicVector::getRealTermsAsIndex(double value){
	value -= offset;
	value /= scalar;
	
	return value;
	
}

double DynamicVector::getValueAtMillis(const double& millis){
	
	int index = round(getRealTermsAsIndex(millis));
		if (index >= 0 && index < length)
			return array[index];
		else
			return 0;
}

void DynamicVector::drawVector(const int& minIndex, const int& maxIndex){

		
		double stepSize = ofGetWidth() / (double)(maxIndex - minIndex);
		double screenHeight = (double) ofGetHeight();
		double maxVal = getMaximum();
	
		int startInt = max(1,minIndex+1);
		int endInt = min(maxIndex, (int)array.size());
		double heightConstant = screenHeight / maxVal;
		int lastHeightPixel = heightConstant * (maxVal - array[startInt-1]);
		int newHeightPixel;
		for (int i = startInt;i < endInt;i++){
			newHeightPixel = (int) heightConstant * (maxVal - array[i]);
			int xPos = i - startInt;
			ofLine (stepSize*(xPos-1), lastHeightPixel, stepSize*xPos, newHeightPixel);
			lastHeightPixel = newHeightPixel;
		}
	
}


void DynamicVector::drawVector(const int& minIndex, const int& maxIndex, ofxWindowRegion window){
	
	
	double stepSize = window.width / (double)(maxIndex - minIndex);
	double screenHeight = (double) window.height;
	double maxVal = getMaximum();
	
	int startInt = max(1,minIndex+1);
	int endInt = min(maxIndex, (int)array.size());
	double heightConstant = screenHeight / maxVal;
	int lastHeightPixel = heightConstant * (maxVal - array[startInt-1]);
	int newHeightPixel;
	for (int i = startInt;i < endInt;i++){
		newHeightPixel = (int) heightConstant * (maxVal - array[i]);
		int xPos = i - startInt;
		ofLine (window.x+stepSize*(xPos-1), window.y+lastHeightPixel, window.x+stepSize*xPos, window.y+newHeightPixel);
		lastHeightPixel = newHeightPixel;
	}
	ofDrawBitmapString("start index "+ofToString(startInt)+" end "+ofToString(endInt), window.x+20, window.y+20);
}


void DynamicVector::drawConstrainedVector(const int& minIndex, const int& maxIndex, const int& minScreenIndex, const int& maxScreenIndex){
	//constrain the height and width
	
	double stepSize = (maxScreenIndex - minScreenIndex) / (double)(maxIndex - minIndex);//step size in pixels per array bin 
	double screenHeight = ofGetHeight();
	double maxVal = getMaximum();

	//OPTIMIZE!! XXX could just add stepsize each time
	//not add minindex each time
	int i = max(1,minIndex+1);
//	ofDrawBitmapString("i = "+ofToString(i)+"  :: screen min: "+ofToString(minScreenIndex + stepSize*(i-minIndex-1)), 20, 640);
	
	while ((minScreenIndex + stepSize*(i-minIndex)) < 0)
		i++;//only draw what is on the screen
	
	for ( ; i < min(maxIndex+1, (int)array.size());i++){
		ofLine (minScreenIndex + (stepSize*(i-minIndex-1)), screenHeight * (1 - array[i-1] / maxVal), 
				minScreenIndex + (stepSize*(i-minIndex)),	 screenHeight * (1 - array[i] / maxVal) );
		
	}
	
	ofLine(minScreenIndex, screenHeight, minScreenIndex, screenHeight/2);
	ofLine(maxScreenIndex, screenHeight, maxScreenIndex, screenHeight/2);
	
//	ofDrawBitmapString(ofToString(stepSize, 2)+"  "+ofToString(maxScreenIndex - minScreenIndex, 0), 20, 600);

}