Source code for FEPdiaGrabber.target._FEPmatrixMethods
# -*- coding: utf-8 *-*
import numpy
from copy import deepcopy
#import bottleneck as bn
#from sys import exit
#from copy import deepcopy
#from scipy.interpolate import griddata
#own modules
#from diaGrabber import source
from diaGrabber import _utils
[docs]class FEPmatrixMethods(object):
def __init__(self):
pass
[docs] def eval2dGaussianDistribution(self, dEB_types=[], merge = None):
print "eval2dGaussianDistribution"
if self.nBasis != 2:
exit("ERROR: eval2dGaussianDistribution works only for 2d-matrixes")
min_dEB_types = [0.25,0.5,0.75]#are at least necassary to get deviation to normaldistr.
for i in min_dEB_types:
if i not in dEB_types:
dEB_types.append(i)#important to the the standard deviation
#from smalles to biggest diameter
dEB_types.sort()
dEB_types.reverse()#e.q.d_EB_75, d_EB_50, d_EB_25
dEB_types.insert(0,0.95)
if merge == None:
merge = self._merge_dim[0]
else:
if merge not in self._merge_dim:
raise KeyError("given merge is not a member of merge-dims im target")
w = merge._mergeIndex
peak = numpy.nanmax(self.mergeMatrix[w])
peak_location = list(numpy.where(self.mergeMatrix[w]==peak))
peak_location[0] = peak_location[0][0]
peak_location[1] = peak_location[1][0]
##removing all nan-values in array:
clean = self.mergeMatrix[w][~numpy.isnan(self.mergeMatrix[w]).any(1)]
##building histogramm (with 5 equal sized bins) of the distribution of the values
histo = numpy.histogram(clean,5)
##the environmental value is defines as the value of the first (lowest) bin
environment = histo[1][0]
dEB_limit = []
dEB = []
border_locations = []
deviationTh2RealList = []
for i in range(len(dEB_types)):
dEB_limit.append((peak-environment)*dEB_types[i] + environment )
dEB.append(numpy.nan)
border_locations.append([])
deviationTh2RealList.append(0)
n_circle_points = 40#has to be int(multiple of 2)
n_part_diameters = n_circle_points/2
delta_radian = 2*numpy.pi / n_circle_points
for d in range(len(dEB_limit)):
radian = 0.0
for i in range(n_circle_points):
n = 0#at peak_location
at_border = False
while not at_border:
x = peak_location[0] + int(n*numpy.cos(radian))
y = peak_location[1] + int(n*numpy.sin(radian))
try:
if not numpy.isnan(self.mergeMatrix[w][x][y]) and self.mergeMatrix[w][x][y] < dEB_limit[d]:#found border
at_border = True
#if [x,y] not in border_locations[d]:
border_locations[d].append([x,y])
except IndexError:
exit("ERROR: peak of merge %s in matrix is to close to border")
n += 1
radian += delta_radian
border_coords = deepcopy(border_locations)
distances = deepcopy(border_locations)
#reset the middle of the gaussiandistr. as the middle of the first (95%) polygon
for i in range(len(border_locations[0])):
border_coords[0][i] = [self.basisMatrix[0][border_locations[0][i][0]],
self.basisMatrix[1][border_locations[0][i][1]]]
polygon = zip(*border_coords[0])
##build mean ... doest work good
##peak_coord= [sum(polygon[0])/len(polygon[0]),sum(polygon[1])/len(polygon[1])]
##build half of max-min:
minx=min(polygon[0])
maxx=max(polygon[0])
miny=min(polygon[1])
maxy=max(polygon[1])
peak_coord= [minx+(maxx-minx)/2, miny+(maxy-miny)/2]
##calc the x-y-coord of all borders in d_EB_types
##and the distance from eBeammiddle to border
print "Gausian-Diameters are..."
for d in range(1,len(dEB_limit)):
for i in range(len(border_locations[d])):
border_coords[d][i] = [self.basisMatrix[0][border_locations[d][i][0]],
self.basisMatrix[1][border_locations[d][i][1]]]
distances[d][i] = _utils.pktAbs2(border_coords[d][i],peak_coord)
#values for dEB_types[i]
dEB[d] = sum(distances[d]) / len(distances[d])#mean
if dEB_types[d] == 0.50:
st_deviation = dEB[d]
print " for %s -> %s" %(dEB_types[d],dEB[d])
##get max/min diameter-factor and the radian(max_diameter) for each d_EB-area
minMax_dEB = []
for d in range(len(dEB_limit)-1):#for each diameter
minMax_dEB.append([None,None,None,None])#min, max, max/min, radian_of_max
for i in range(n_part_diameters):
# == radius + radius in opposite dir.
part_diameter = distances[d+1][i] + distances[d+1][i+n_part_diameters]
#print minMax_dEB, d
if part_diameter < minMax_dEB[d][0] or minMax_dEB[d][0] == None:#min
minMax_dEB[d][0] = part_diameter
if part_diameter > minMax_dEB[d][1] or minMax_dEB[d][1] == None:#max
minMax_dEB[d][1] = part_diameter
minMax_dEB[d][3] = delta_radian * i
minMax_dEB[d][2] = minMax_dEB[d][1] / minMax_dEB[d][0] #max/min
#get mean max/min-factor and mean radian_of_max
mean_maxMin_dEB = 0
mean_max_radian = 0
for d in range(len(dEB_limit)-1):#for each diameter
mean_maxMin_dEB += minMax_dEB[d][2]
mean_max_radian += minMax_dEB[d][3]
mean_maxMin_dEB /= (len(dEB_limit)-1)
mean_max_radian /= (len(dEB_limit)-1)
print "The mean factor max/min-diameter is %s" % mean_maxMin_dEB
print "The mean radian of the max diameter is %s" % mean_max_radian
#calc deviation to a normaldisribution defined as:
#normal_distribution_fn = \
corrF = 1/ (
1/(st_deviation* ((2*numpy.pi)**0.5) ) ) * numpy.e**(-0.5 *
( ( 0 / st_deviation)**2 ) ) * (peak-environment)
for d in range(1,len(dEB_limit)):
for x in distances[d]:
dEB_th = corrF * (
1/(st_deviation* ((2*numpy.pi)**0.5) ) ) * numpy.e**(-0.5 *
( ( x / st_deviation)**2 ) ) * (peak-environment)
deviationTh2RealList[d] += abs( (dEB_th - dEB_limit[d]) / dEB_th)#dEB_limit[d]/dEB_th#
deviationTh2RealList[d] /= (len(distances[d])-1) #middle fo difference between theory and reality
##middle of all deviantion-rings
relDevianceTh2Real = sum(deviationTh2RealList) / ( len(deviationTh2RealList))
print "The relative deviation/variance between theoretical an real distribution is %s" %relDevianceTh2Real
#add point middle of gaussian:
peak_coord.append("Mittelpunkt")
self.plot_overlay[w][0].append(peak_coord)
#self.plot_overlay[w][0].append("Mittelpunkt")#x-,y-list
###redundant!!!->
self.plot_overlay[w][6].append("Mittelpunkt")#x-,y-list
self.plot_overlay[w][5].append("deviation= %s" %relDevianceTh2Real)
self.plot_overlay[w][5].append("maxRad= %s" %mean_max_radian)
self.plot_overlay[w][5].append("max/min-d-factor= %s" %mean_maxMin_dEB)
self.plot_overlay[w][5].append("d_EB50= %s" %dEB[dEB_types.index(0.5)])
#add rings
for d in range(1,len(dEB_limit)):
self.plot_overlay[w][1].append(zip(*border_coords[d]))#x-,y-list
self.plot_overlay[w][1][-1].append("d_EB= %s" %(dEB_types[d]) )#names for legend
###redundant!!!->
self.plot_overlay[w][6].append("d_EB= %s" %(dEB_types[d]) )#x-,y-list
print "done"
