[70] | 1 | #include <cpgplot.h> |
---|
[3] | 2 | #include <iostream> |
---|
| 3 | #include <math.h> |
---|
| 4 | #include <Utils/utils.hh> |
---|
| 5 | |
---|
[70] | 6 | void findMinMax(const float *array, const int size, float &min, float &max) |
---|
| 7 | { |
---|
| 8 | min = max = array[0]; |
---|
| 9 | for(int i=1;i<size;i++) { |
---|
| 10 | if(array[i]<min) min=array[i]; |
---|
| 11 | if(array[i]>max) max=array[i]; |
---|
| 12 | } |
---|
| 13 | } |
---|
| 14 | |
---|
[53] | 15 | float findMean(float *&array, int size) |
---|
| 16 | { |
---|
| 17 | float mean = array[0]; |
---|
| 18 | for(int i=1;i<size;i++) mean += array[i]; |
---|
| 19 | mean /= float(size); |
---|
| 20 | return mean; |
---|
| 21 | } |
---|
| 22 | |
---|
[78] | 23 | float findStddev(float *&array, int size) |
---|
[53] | 24 | { |
---|
| 25 | float mean = findMean(array,size); |
---|
[78] | 26 | float stddev = (array[0]-mean) * (array[0]-mean); |
---|
| 27 | for(int i=1;i<size;i++) stddev += (array[i]-mean)*(array[i]-mean); |
---|
[79] | 28 | stddev = sqrt(stddev/float(size-1)); |
---|
[78] | 29 | return stddev; |
---|
[53] | 30 | } |
---|
| 31 | |
---|
| 32 | float findMedian(float *&array, int size) |
---|
| 33 | { |
---|
| 34 | // NOTE: madfm = median absolute deviation from median |
---|
| 35 | float *newarray = new float[size]; |
---|
| 36 | float median; |
---|
| 37 | for(int i=0;i<size;i++) newarray[i] = array[i]; |
---|
[67] | 38 | sort(newarray,0,size); |
---|
[53] | 39 | if((size%2)==0) median = 0.5*(newarray[size/2-1]+newarray[size/2]); |
---|
| 40 | else median = newarray[size/2]; |
---|
| 41 | delete [] newarray; |
---|
| 42 | return median; |
---|
| 43 | } |
---|
| 44 | |
---|
| 45 | float findMADFM(float *&array, int size) |
---|
| 46 | { |
---|
| 47 | // NOTE: madfm = median absolute deviation from median |
---|
| 48 | float *newarray = new float[size]; |
---|
| 49 | float median = findMedian(array,size); |
---|
| 50 | float madfm; |
---|
| 51 | for(int i=0;i<size;i++) newarray[i] = fabs(array[i]-median); |
---|
[67] | 52 | sort(newarray,0,size); |
---|
[53] | 53 | if((size%2)==0) madfm = 0.5*(newarray[size/2-1]+newarray[size/2]); |
---|
| 54 | else madfm = newarray[size/2]; |
---|
| 55 | delete [] newarray; |
---|
| 56 | return madfm; |
---|
| 57 | } |
---|
| 58 | |
---|
[3] | 59 | void findMedianStats(float *&array, int size, float &median, float &madfm) |
---|
| 60 | { |
---|
| 61 | // NOTE: madfm = median absolute deviation from median |
---|
| 62 | float *newarray = new float[size]; |
---|
| 63 | |
---|
| 64 | for(int i=0;i<size;i++) newarray[i] = array[i]; |
---|
[67] | 65 | sort(newarray,0,size); |
---|
[3] | 66 | if((size%2)==0) median = 0.5*(newarray[size/2-1]+newarray[size/2]); |
---|
| 67 | else median = newarray[size/2]; |
---|
| 68 | |
---|
| 69 | for(int i=0;i<size;i++) newarray[i] = fabs(array[i]-median); |
---|
[67] | 70 | sort(newarray,0,size); |
---|
[3] | 71 | if((size%2)==0) madfm = 0.5*(newarray[size/2-1]+newarray[size/2]); |
---|
| 72 | else madfm = newarray[size/2]; |
---|
| 73 | |
---|
| 74 | delete [] newarray; |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | void findMedianStats(float *&array, long size, float &median, float &madfm) |
---|
| 78 | { |
---|
| 79 | // NOTE: madfm = median absolute deviation from median |
---|
| 80 | float *newarray = new float[size]; |
---|
| 81 | |
---|
| 82 | for(int i=0;i<size;i++) newarray[i] = array[i]; |
---|
[67] | 83 | sort(newarray,0,size); |
---|
[3] | 84 | if((size%2)==0) median = 0.5*(newarray[size/2-1]+newarray[size/2]); |
---|
| 85 | else median = newarray[size/2]; |
---|
| 86 | |
---|
| 87 | for(int i=0;i<size;i++) newarray[i] = fabs(array[i]-median); |
---|
[67] | 88 | sort(newarray,0,size); |
---|
[3] | 89 | if((size%2)==0) madfm = 0.5*(newarray[size/2-1]+newarray[size/2]); |
---|
| 90 | else madfm = newarray[size/2]; |
---|
| 91 | |
---|
| 92 | delete [] newarray; |
---|
| 93 | } |
---|
| 94 | |
---|
| 95 | |
---|
[79] | 96 | void findNormalStats(float *&array, int size, float &mean, float &stddev) |
---|
[3] | 97 | { |
---|
| 98 | mean = array[0]; |
---|
| 99 | for(int i=1;i<size;i++){ |
---|
| 100 | mean += array[i]; |
---|
| 101 | } |
---|
| 102 | mean /= float(size); |
---|
| 103 | |
---|
[79] | 104 | stddev = (array[0]-mean) * (array[0]-mean); |
---|
| 105 | for(int i=1;i<size;i++) stddev += (array[i]-mean)*(array[i]-mean); |
---|
| 106 | stddev = sqrt(stddev/float(size-1)); |
---|
[3] | 107 | |
---|
| 108 | } |
---|
| 109 | |
---|
[70] | 110 | void findTrimmedHistStatsOLD(float *array, const int size, float &tmean, float &tsigma) |
---|
| 111 | { |
---|
| 112 | |
---|
| 113 | const int nbin = 100; |
---|
| 114 | float *num = new float[nbin]; |
---|
| 115 | bool *keep = new bool[nbin]; |
---|
| 116 | |
---|
| 117 | // HOW MANY VALUES IN EACH BIN? |
---|
| 118 | float min,max; |
---|
| 119 | findMinMax(array,size,min,max); |
---|
| 120 | for(int i=0; i<nbin; i++) num[i]=0; |
---|
| 121 | for(int i=0; i<size; i++){ |
---|
| 122 | float fraction = (array[i] - min) / (max - min); |
---|
| 123 | int bin = (int)( floor(fraction*nbin) ); |
---|
| 124 | if((bin>=0)&&(bin<nbin)) num[bin]++; |
---|
| 125 | } |
---|
| 126 | int binmax = 0; |
---|
| 127 | for(int i=1; i<nbin; i++) if(num[i]>num[binmax]) binmax = i; |
---|
| 128 | for(int i=0; i<nbin; i++) keep[i] = (num[i]>=num[binmax]/2); |
---|
| 129 | float *newarray = new float[size]; |
---|
| 130 | int newsize = 0; |
---|
| 131 | for(int i=0; i<size; i++){ |
---|
| 132 | float fraction = (array[i] - min) / (max - min); |
---|
| 133 | int bin = (int)( floor(fraction*nbin) ); |
---|
| 134 | if(keep[bin]) newarray[newsize++] = array[i]; |
---|
| 135 | } |
---|
| 136 | |
---|
| 137 | findNormalStats(newarray,newsize,tmean,tsigma); |
---|
| 138 | delete [] num,keep,newarray; |
---|
| 139 | |
---|
| 140 | } |
---|
| 141 | void findTrimmedHistStats2(float *array, const int size, float &tmean, float &tsigma) |
---|
| 142 | { |
---|
| 143 | |
---|
| 144 | const int nbin = 200; |
---|
| 145 | float *num = new float[nbin]; |
---|
| 146 | bool *keep = new bool[nbin]; |
---|
| 147 | |
---|
| 148 | // HOW MANY VALUES IN EACH BIN? |
---|
| 149 | float min,max; |
---|
| 150 | findMinMax(array,size,min,max); |
---|
| 151 | for(int i=0; i<nbin; i++) num[i]=0; |
---|
| 152 | for(int i=0; i<size; i++){ |
---|
| 153 | float fraction = (array[i] - min) / (max - min); |
---|
| 154 | int bin = (int)( floor(fraction*nbin) ); |
---|
| 155 | if((bin>=0)&&(bin<nbin)) num[bin]++; |
---|
| 156 | } |
---|
| 157 | int binmax = 0; |
---|
| 158 | for(int i=1; i<nbin; i++) if(num[i]>num[binmax]) binmax = i; |
---|
| 159 | for(int i=0; i<nbin; i++) keep[i] = (num[i]>=num[binmax]/2); |
---|
| 160 | float *newarray = new float[size]; |
---|
| 161 | int newsize = 0; |
---|
| 162 | for(int i=0; i<size; i++){ |
---|
| 163 | float fraction = (array[i] - min) / (max - min); |
---|
| 164 | int bin = (int)( floor(fraction*nbin) ); |
---|
| 165 | if(keep[bin]) newarray[newsize++] = array[i]; |
---|
| 166 | } |
---|
| 167 | |
---|
| 168 | tmean = findMean(newarray,newsize); |
---|
| 169 | |
---|
| 170 | tsigma = newsize * (max-min)/float(nbin) / |
---|
| 171 | (num[binmax] * erf(sqrt(M_LN2)) * sqrt(2.*M_PI)); |
---|
| 172 | |
---|
| 173 | } |
---|
| 174 | |
---|
| 175 | void findTrimmedHistStats(float *array, const int size, float &tmean, float &tsigma) |
---|
| 176 | { |
---|
| 177 | float datamin,datamax; |
---|
| 178 | findMinMax(array,size,datamin,datamax); |
---|
| 179 | float *sorted = new float[size]; |
---|
| 180 | float *cumul = new float[size]; |
---|
| 181 | float *angle = new float[size]; |
---|
| 182 | for(int i=0;i<size;i++) sorted[i] = array[i]/(datamax-datamin); |
---|
| 183 | sort(sorted,0,size); |
---|
| 184 | for(int i=0;i<size;i++) cumul[i] = (float)i/(float)size; |
---|
| 185 | int width =(int)( 20. * log10((float)size)); |
---|
| 186 | for(int i=0;i<size;i++){ |
---|
| 187 | int beg,end; |
---|
| 188 | float slope,eSlope,inter,eInter,r; |
---|
| 189 | if(i<width/2) beg=0; else beg=i-width/2; |
---|
| 190 | if(i>=size-width/2) end=size-1; else end=i+width/2; |
---|
| 191 | if(linear_regression(size,sorted,cumul,beg,end,slope,eSlope,inter,eInter,r)==0) |
---|
[107] | 192 | angle[i] = atan( fabs(slope) ) * 180. / M_PI; |
---|
[70] | 193 | else angle[i] = 90.; |
---|
| 194 | } |
---|
| 195 | |
---|
| 196 | // int start = 0; |
---|
| 197 | // while(angle[start] < 45.) start++; |
---|
| 198 | // int finish = size-1; |
---|
| 199 | // while(angle[finish] < 45.) finish--; |
---|
| 200 | |
---|
| 201 | int maxpt = 0; |
---|
| 202 | for(int i = 1; i<size; i++) if(angle[i]>angle[maxpt]) maxpt=i; |
---|
| 203 | int start = maxpt; |
---|
| 204 | while((start>0)&&(angle[start]>45.)) start--; |
---|
| 205 | int finish = maxpt; |
---|
| 206 | while((finish < size-1)&&(angle[finish]>45.)) finish++; |
---|
| 207 | |
---|
[188] | 208 | // std::cerr << "npts = " << size << ", start = " << start << ", finish = " << finish << std::endl; |
---|
[70] | 209 | |
---|
| 210 | int trimSize=0; |
---|
| 211 | float *newarray = new float[finish-start+1]; |
---|
| 212 | for(int i=0;i<finish-start+1;i++) newarray[i] = sorted[i+start]*(datamax-datamin); |
---|
| 213 | |
---|
| 214 | findNormalStats(newarray,finish-start+1,tmean,tsigma); |
---|
| 215 | |
---|
| 216 | } |
---|