andrew@2
|
1 /*
|
andrew@2
|
2 * bayesianArray.cpp
|
andrew@2
|
3 * bayesianTest5
|
andrew@2
|
4 *
|
andrew@2
|
5 * Created by Andrew Robertson on 08/05/2010.
|
andrew@2
|
6 * Copyright 2010 __MyCompanyName__. All rights reserved.
|
andrew@2
|
7 *
|
andrew@2
|
8 */
|
andrew@2
|
9
|
andrew@2
|
10 #include "bayesianArray.h"
|
andrew@2
|
11 #include "math.h"
|
andrew@2
|
12 #include "ofMain.h"
|
andrew@2
|
13
|
andrew@2
|
14 bayesianArray::bayesianArray(){
|
andrew@2
|
15 likelihoodNoise = 0.5;
|
andrew@6
|
16 likelihoodMean = arraySize/2;
|
andrew@6
|
17 likelihoodStdDev = arraySize / 12;
|
andrew@2
|
18 initialiseArray();
|
andrew@2
|
19 }
|
andrew@2
|
20
|
andrew@2
|
21 void bayesianArray::initialiseArray(){
|
andrew@2
|
22
|
andrew@2
|
23 //maximumIndex = 12;//change this
|
andrew@6
|
24 setGaussianPrior(arraySize/2, arraySize/1);
|
andrew@6
|
25 setGaussianLikelihood(arraySize/2, arraySize/1);//likelihoodMean, likelihoodStdDev);
|
andrew@2
|
26
|
andrew@2
|
27 calculatePosterior();
|
andrew@2
|
28 renormalisePosterior();
|
andrew@2
|
29 posteriorDecayRate = 0.06;
|
andrew@2
|
30
|
andrew@2
|
31 eighthNoteProportion = 0.35;//must be less than 0.5 to discriminate - was 0.4
|
andrew@2
|
32 earlySixteenthNoteProportion = 0;
|
andrew@2
|
33 lateSixteenthNoteProportion = 0;
|
andrew@2
|
34 decayNoiseAmount = 0.1;
|
andrew@6
|
35 decayNoiseStdDev = arraySize/24;
|
andrew@2
|
36 standardDeviation = likelihoodStdDev;
|
andrew@6
|
37 setDecayNoiseGaussian(arraySize/2, decayNoiseStdDev);
|
andrew@2
|
38
|
andrew@2
|
39 setGaussianLikelihood(likelihoodMean, likelihoodStdDev);
|
andrew@2
|
40 }
|
andrew@2
|
41
|
andrew@2
|
42
|
andrew@2
|
43 void bayesianArray::setGaussianPrior(float mean, float StdDev){
|
andrew@2
|
44 int i;
|
andrew@6
|
45 for (i=0;i<arraySize;i++){
|
andrew@2
|
46 prior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
|
andrew@2
|
47 //posterior[i] = prior[i];
|
andrew@2
|
48 }
|
andrew@2
|
49 }
|
andrew@2
|
50
|
andrew@2
|
51 void bayesianArray::setGaussianPosterior(float mean, float StdDev){
|
andrew@2
|
52 int i;
|
andrew@6
|
53 for (i=0;i<arraySize;i++){
|
andrew@2
|
54 posterior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
|
andrew@2
|
55 }
|
andrew@2
|
56 }
|
andrew@2
|
57
|
andrew@2
|
58 void bayesianArray::setGaussianLikelihoodForBeats(float mean, float StdDev){
|
andrew@2
|
59 //this has eighth and sixteenth positions included
|
andrew@2
|
60
|
andrew@6
|
61 if (mean >= 0 && mean <= arraySize){
|
andrew@2
|
62 int i; float eighthDifference;
|
andrew@4
|
63 //main beat position plus these three:
|
andrew@6
|
64 int eighthPosition = ((int)mean + arraySize/2)%arraySize;
|
andrew@6
|
65 int earlySixteenthPosition = ((int)mean + (3*arraySize/4))%arraySize;;
|
andrew@6
|
66 int lateSixteenthPosition = ((int)mean + (arraySize/4))%arraySize;;
|
andrew@2
|
67
|
andrew@2
|
68 float mainDifference, sixteenthDifference;
|
andrew@2
|
69 float gaussianProportion = 1 - likelihoodNoise;
|
andrew@2
|
70 float mainProportion = (1 - eighthNoteProportion - earlySixteenthNoteProportion - lateSixteenthNoteProportion);
|
andrew@2
|
71
|
andrew@6
|
72 for (i=0;i < arraySize;i++){
|
andrew@2
|
73
|
andrew@6
|
74 mainDifference = min( fabs(i-mean) , (double)(i + arraySize - mean));
|
andrew@2
|
75 likelihood[i] = gaussianProportion * mainProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
|
andrew@2
|
76
|
andrew@6
|
77 eighthDifference = min( abs(i - eighthPosition) , i + arraySize - eighthPosition);
|
andrew@6
|
78 eighthDifference = min(eighthDifference , (float)(arraySize + eighthPosition - i ));
|
andrew@2
|
79 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
|
andrew@2
|
80 likelihood[i] += gaussianProportion * eighthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(eighthDifference)*(eighthDifference)/(2*StdDev*StdDev)) ;
|
andrew@2
|
81
|
andrew@6
|
82 sixteenthDifference = min( abs(i - earlySixteenthPosition) , i + arraySize - earlySixteenthPosition);
|
andrew@6
|
83 sixteenthDifference = min(sixteenthDifference , (float)(arraySize + earlySixteenthPosition - i ));
|
andrew@2
|
84 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
|
andrew@2
|
85 likelihood[i] += gaussianProportion * earlySixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
|
andrew@2
|
86
|
andrew@6
|
87 sixteenthDifference = min( abs(i - lateSixteenthPosition) , i + arraySize - lateSixteenthPosition);
|
andrew@6
|
88 sixteenthDifference = min(sixteenthDifference , (float)(arraySize + lateSixteenthPosition - i ));
|
andrew@2
|
89 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
|
andrew@2
|
90 likelihood[i] += gaussianProportion * lateSixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
|
andrew@2
|
91
|
andrew@2
|
92
|
andrew@2
|
93
|
andrew@6
|
94 likelihood[i] += (likelihoodNoise / arraySize);
|
andrew@2
|
95 //likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) ,
|
andrew@6
|
96 //(double) (likelihoodNoise / arraySize) );
|
andrew@2
|
97 }
|
andrew@6
|
98 // renormaliseArray(&likelihood[0], arraySize);
|
andrew@2
|
99 }//end if mean within limits
|
andrew@2
|
100 }
|
andrew@2
|
101
|
andrew@2
|
102
|
andrew@2
|
103 void bayesianArray::setGaussianLikelihood(float mean, float StdDev){
|
andrew@6
|
104 if (mean >= 0 && mean <= arraySize){
|
andrew@2
|
105 int i; float eighthDifference;
|
andrew@6
|
106 int eighthPosition = ((int)mean + arraySize/2)%arraySize;
|
andrew@2
|
107 float mainDifference;
|
andrew@2
|
108 float gaussianProportion = 1 - likelihoodNoise;
|
andrew@2
|
109
|
andrew@6
|
110 for (i=0;i < arraySize;i++){
|
andrew@2
|
111
|
andrew@6
|
112 mainDifference = min( fabs(i-mean) , (double)(i + arraySize - mean));
|
andrew@2
|
113 //without * (1 - eighthNoteProportion)
|
andrew@2
|
114 likelihood[i] = gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
|
andrew@2
|
115
|
andrew@6
|
116 likelihood[i] += (likelihoodNoise / arraySize);
|
andrew@2
|
117 //likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) ,
|
andrew@6
|
118 //(double) (likelihoodNoise / arraySize) );
|
andrew@2
|
119 }
|
andrew@6
|
120 // renormaliseArray(&likelihood[0], arraySize);
|
andrew@2
|
121 }//end if mean within limits
|
andrew@2
|
122 }
|
andrew@2
|
123
|
andrew@2
|
124 void bayesianArray::calculatePosterior(){
|
andrew@2
|
125 int i;
|
andrew@6
|
126 for (i=0;i < arraySize;i++){
|
andrew@2
|
127 posterior[i] = likelihood[i] * prior[i];
|
andrew@2
|
128 }
|
andrew@2
|
129 //renormalisePosterior();
|
andrew@2
|
130 }
|
andrew@2
|
131
|
andrew@2
|
132
|
andrew@2
|
133 float bayesianArray::getMaximum(float *ptr, int length){
|
andrew@2
|
134 int i;
|
andrew@2
|
135 float max = 0;
|
andrew@2
|
136 for (i=0;i < length;i++){
|
andrew@2
|
137 if (*(ptr+i)>max)
|
andrew@2
|
138 max = *(ptr+i);
|
andrew@2
|
139 }
|
andrew@2
|
140 maximumValue = max;
|
andrew@2
|
141 return max;
|
andrew@2
|
142 }
|
andrew@2
|
143
|
andrew@2
|
144 float* bayesianArray::getMaximumEstimate(float *ptr, int length){
|
andrew@2
|
145 float returnArray[2];
|
andrew@2
|
146 int i;
|
andrew@2
|
147 float max = 0;
|
andrew@2
|
148 maximumIndex = 0;
|
andrew@2
|
149 for (i=0;i < length;i++){
|
andrew@2
|
150 if (*(ptr+i)>max){
|
andrew@2
|
151 max = *(ptr+i);
|
andrew@2
|
152 maximumIndex = i;
|
andrew@2
|
153 }
|
andrew@2
|
154 }
|
andrew@2
|
155 returnArray[0] = max;
|
andrew@2
|
156 returnArray[1] = maximumIndex;
|
andrew@2
|
157 maximumValue = max;
|
andrew@2
|
158 return &returnArray[0];
|
andrew@2
|
159 }
|
andrew@2
|
160
|
andrew@2
|
161
|
andrew@2
|
162
|
andrew@2
|
163 double bayesianArray::getIntegratedEstimateIndex(){
|
andrew@2
|
164 int i;
|
andrew@2
|
165 float integratedQuantity = 0;
|
andrew@2
|
166 float integratedTotal = 0;
|
andrew@2
|
167 double integratedIndex = 0;
|
andrew@6
|
168 for (i=0;i < arraySize;i++){
|
andrew@2
|
169 integratedQuantity += posterior[i];//the values of the probability distribution
|
andrew@2
|
170 integratedTotal += i*posterior[i];
|
andrew@2
|
171 }
|
andrew@2
|
172 if (integratedQuantity > 0){
|
andrew@2
|
173 integratedIndex = integratedTotal / integratedQuantity;
|
andrew@2
|
174 }
|
andrew@2
|
175 integratedEstimate = (float) integratedIndex;
|
andrew@2
|
176 return integratedIndex;
|
andrew@2
|
177 }
|
andrew@2
|
178
|
andrew@2
|
179
|
andrew@2
|
180 double bayesianArray::calculateStandardDeviation(){
|
andrew@2
|
181
|
andrew@2
|
182 double total = 0;
|
andrew@2
|
183 double pdfSum;
|
andrew@2
|
184 double variance = 0;
|
andrew@6
|
185 for (int i=0;i < arraySize;i++){
|
andrew@2
|
186 //*posterior[i] *
|
andrew@2
|
187 total += posterior[i] * (i - integratedEstimate) * (i - integratedEstimate);//the values of the probability distribution
|
andrew@2
|
188 pdfSum += posterior[i];
|
andrew@2
|
189 }
|
andrew@2
|
190
|
andrew@2
|
191 if (pdfSum > 0)
|
andrew@2
|
192 variance = total / pdfSum;
|
andrew@2
|
193 else
|
andrew@6
|
194 variance = arraySize;
|
andrew@2
|
195
|
andrew@2
|
196 standardDeviation = sqrt(variance);
|
andrew@2
|
197 return standardDeviation;
|
andrew@2
|
198 }
|
andrew@2
|
199
|
andrew@2
|
200
|
andrew@2
|
201
|
andrew@2
|
202 void bayesianArray::renormaliseArray(float *ptr, int length){
|
andrew@2
|
203 int i;
|
andrew@2
|
204 float totalArea = 0;
|
andrew@2
|
205 for (i=0;i < length;i++){
|
andrew@2
|
206 totalArea += *(ptr+i);
|
andrew@2
|
207 }
|
andrew@2
|
208
|
andrew@2
|
209 for (i=0;i < length;i++){
|
andrew@2
|
210 *(ptr+i) /= totalArea;
|
andrew@2
|
211 }
|
andrew@2
|
212
|
andrew@2
|
213 }
|
andrew@2
|
214
|
andrew@2
|
215 void bayesianArray::resetPrior(){
|
andrew@2
|
216 int i;
|
andrew@6
|
217 for (i=0;i<arraySize;i++){
|
andrew@2
|
218 prior[i] = posterior[i];
|
andrew@2
|
219 }
|
andrew@2
|
220 }
|
andrew@2
|
221
|
andrew@2
|
222 void bayesianArray::renormalisePrior(){
|
andrew@2
|
223 int i;
|
andrew@2
|
224 float totalArea = 0;
|
andrew@6
|
225 for (i=0;i < arraySize;i++){
|
andrew@2
|
226 totalArea += prior[i];
|
andrew@2
|
227 }
|
andrew@2
|
228
|
andrew@6
|
229 for (i=0;i < arraySize;i++){
|
andrew@2
|
230 prior[i] /= totalArea;
|
andrew@2
|
231 }
|
andrew@2
|
232 }
|
andrew@2
|
233
|
andrew@2
|
234 void bayesianArray::renormalisePosterior(){
|
andrew@2
|
235 int i;
|
andrew@2
|
236 float totalArea = 0;
|
andrew@6
|
237 for (i=0;i < arraySize;i++){
|
andrew@2
|
238 totalArea += posterior[i];
|
andrew@2
|
239 }
|
andrew@2
|
240
|
andrew@6
|
241 for (i=0;i < arraySize;i++){
|
andrew@2
|
242 posterior[i] /= totalArea;
|
andrew@2
|
243 }
|
andrew@2
|
244 }
|
andrew@2
|
245
|
andrew@2
|
246 void bayesianArray::decayPosterior(){
|
andrew@2
|
247 float *pointer;
|
andrew@6
|
248 pointer = getMaximumEstimate(&posterior[0], arraySize);
|
andrew@2
|
249 float maximum;
|
andrew@2
|
250 maximum = *pointer;
|
andrew@2
|
251 int i;
|
andrew@6
|
252 for (i=0;i<arraySize;i++){
|
andrew@2
|
253 posterior[i] += (maximum - posterior[i]) * posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
|
andrew@2
|
254 }
|
andrew@2
|
255 maximumIndex = *(pointer+1);
|
andrew@2
|
256 }
|
andrew@2
|
257
|
andrew@2
|
258 void bayesianArray::setDecayNoiseGaussian(float mean, float StdDev){
|
andrew@2
|
259 int i;
|
andrew@6
|
260 for (i=0;i<arraySize;i++){
|
andrew@2
|
261 decayNoiseArray[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
|
andrew@2
|
262 }
|
andrew@2
|
263 }
|
andrew@2
|
264
|
andrew@2
|
265 void bayesianArray::decayPosteriorWithGaussianNoise(){
|
andrew@2
|
266
|
andrew@2
|
267 int i;
|
andrew@6
|
268 float currentMaximum = getMaximum(&posterior[0], arraySize);
|
andrew@6
|
269 for (i=0;i<arraySize;i++){
|
andrew@6
|
270 posterior[i] += decayNoiseArray[(i - (int)maximumIndex + ((3*arraySize)/2)) % arraySize] * currentMaximum * decayNoiseAmount;
|
andrew@2
|
271 //posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
|
andrew@2
|
272 }
|
andrew@2
|
273
|
andrew@2
|
274 }
|
andrew@2
|
275
|
andrew@2
|
276 void bayesianArray::resetMaximumPosterior(){
|
andrew@2
|
277 int i;
|
andrew@2
|
278 float max = 0;
|
andrew@6
|
279 for (i=0;i < arraySize;i++){
|
andrew@2
|
280 if (posterior[i]>max){
|
andrew@2
|
281 maximumIndex = i;
|
andrew@2
|
282 max = posterior[i];
|
andrew@2
|
283 }
|
andrew@2
|
284 }
|
andrew@2
|
285 }
|
andrew@2
|
286
|
andrew@2
|
287 void bayesianArray::translateDistribution(int translationIndex){
|
andrew@2
|
288 int tmpIndex;
|
andrew@2
|
289 //copy array
|
andrew@2
|
290 int i;
|
andrew@6
|
291 for (i=0;i < arraySize;i++){
|
andrew@2
|
292 tempPosteriorArray[i] = posterior[i] ;
|
andrew@2
|
293 }
|
andrew@2
|
294 //translate values
|
andrew@6
|
295 for (i=0;i < arraySize;i++){
|
andrew@6
|
296 tmpIndex = (i + translationIndex + arraySize)%arraySize;
|
andrew@2
|
297 posterior[tmpIndex] = tempPosteriorArray[i];
|
andrew@2
|
298 }
|
andrew@2
|
299 //now delete tmp array
|
andrew@2
|
300 }
|
andrew@2
|
301
|
andrew@2
|
302 double bayesianArray::getKLdivergence(){
|
andrew@2
|
303 double KLsum = 0;
|
andrew@2
|
304 //take no chances - renormalise both prior and posterior
|
andrew@2
|
305 renormalisePosterior();
|
andrew@2
|
306 renormalisePrior();
|
andrew@6
|
307 for (int i = 0;i < arraySize;i++){
|
andrew@2
|
308 if (posterior[i] > 0 && prior[i] > 0){
|
andrew@2
|
309 KLsum += (posterior[i]*log(posterior[i]/prior[i]));
|
andrew@2
|
310 }
|
andrew@2
|
311 }
|
andrew@2
|
312 return KLsum;
|
andrew@2
|
313 }
|
andrew@2
|
314
|
andrew@3
|
315 double bayesianArray::getEntropyOfPosterior(){
|
andrew@3
|
316 double entropy = 0;
|
andrew@4
|
317 //make sure normalised? (it is)
|
andrew@6
|
318 for (int i = 0;i < arraySize;i++){
|
andrew@3
|
319 if (posterior[i] > 0){
|
andrew@3
|
320 entropy += posterior[i]*log(posterior[i]);
|
andrew@3
|
321 }
|
andrew@3
|
322 }
|
andrew@3
|
323 return -1*entropy;
|
andrew@3
|
324 }
|
andrew@3
|
325
|
andrew@11
|
326 double bayesianArray::getEntropyOfPrior(){
|
andrew@11
|
327 double entropy = 0;
|
andrew@11
|
328 //make sure normalised? (it is)
|
andrew@11
|
329 for (int i = 0;i < arraySize;i++){
|
andrew@11
|
330 if (posterior[i] > 0){
|
andrew@11
|
331 entropy += prior[i]*log(prior[i]);
|
andrew@11
|
332 }
|
andrew@11
|
333 }
|
andrew@11
|
334 return -1*entropy;
|
andrew@11
|
335 }
|
andrew@3
|
336
|
andrew@11
|
337
|