Markus Rosenstihl
289 lines
9.1 KiB
Python
289 lines
9.1 KiB
Python
import threading
|
|
import math
|
|
import types
|
|
import sys
|
|
import tables
|
|
import numpy
|
|
import collections
|
|
from . import Drawable
|
|
|
|
## provide gaussian statistics for a series of measured data points
|
|
#
|
|
# AccumulatedValue provides mean and error of mean after being fed with measured data
|
|
# internaly it keeps the sum, the sum of squares and the number of data points
|
|
class AccumulatedValue:
|
|
|
|
def __init__(self, mean=None, mean_err=None, n=None):
|
|
"""
|
|
one value with std. deviation
|
|
can be initialized by:
|
|
No argument: no entries
|
|
one argument: first entry
|
|
two arguments: mean and its error, n is set 2
|
|
three arguments: already existing statistics defined by mean, mean's error, n
|
|
"""
|
|
if mean is None:
|
|
self.y=0.0
|
|
self.y2=0.0
|
|
self.n=0
|
|
elif mean_err is None and n is None:
|
|
self.y=float(mean)
|
|
self.y2=self.y**2
|
|
self.n=1
|
|
elif mean_err is None:
|
|
self.n=max(1, int(n))
|
|
self.y=float(mean)*self.n
|
|
self.y2=(float(mean)**2)*self.n
|
|
elif n is None:
|
|
self.n=2
|
|
self.y=float(mean)*2
|
|
self.y2=(float(mean_err)**2+float(mean)**2)*2
|
|
else:
|
|
self.n=int(n)
|
|
self.y=float(mean)*self.n
|
|
self.y2=float(mean_err)**2*n*(n-1.0)+float(mean)**2*n
|
|
|
|
def __add__(self,y):
|
|
new_one=AccumulatedValue()
|
|
if (type(y) is types.InstanceType and isinstance(y, AccumulatedValue)):
|
|
new_one.y=self.y+y.y
|
|
new_one.y2=self.y2+y.y2
|
|
new_one.n=self.n+y.n
|
|
else:
|
|
new_one.y=self.y+float(y)
|
|
new_one.y2=self.y2+float(y)**2
|
|
new_one.n=self.n+1
|
|
return new_one
|
|
|
|
def __iadd__(self,y):
|
|
if (type(y) is types.InstanceType and isinstance(y, AccumulatedValue)):
|
|
self.y+=y.y
|
|
self.y2+=y.y2
|
|
self.n+=y.n
|
|
else:
|
|
self.y+=float(y)
|
|
self.y2+=float(y)**2
|
|
self.n+=1
|
|
return self
|
|
|
|
def copy(self):
|
|
a=AccumulatedValue()
|
|
a.y=self.y
|
|
a.y2=self.y2
|
|
a.n=self.n
|
|
return a
|
|
|
|
def mean(self):
|
|
"""
|
|
returns the mean of all added/accumulated values
|
|
"""
|
|
if self.n is None or self.n==0:
|
|
return None
|
|
else:
|
|
return self.y/self.n
|
|
|
|
def sigma(self):
|
|
"""
|
|
returns the standard deviation added/accumulated values
|
|
"""
|
|
if self.n>1:
|
|
variance=(self.y2-(self.y**2)/float(self.n))/(self.n-1.0)
|
|
if variance<0:
|
|
if variance<-1e-20:
|
|
print("variance=%g<0! assuming 0"%variance)
|
|
return 0.0
|
|
return math.sqrt(variance)
|
|
elif self.n==1:
|
|
return 0.0
|
|
else:
|
|
return None
|
|
|
|
def mean_error(self):
|
|
"""
|
|
returns the mean's error (=std.dev/sqrt(n)) of all added/accumulated values
|
|
"""
|
|
if self.n>1:
|
|
variance=(self.y2-(self.y**2)/float(self.n))/(self.n-1.0)
|
|
if variance<0:
|
|
if variance<-1e-20:
|
|
print("variance=%g<0! assuming 0"%variance)
|
|
return 0.0
|
|
return math.sqrt(variance/self.n)
|
|
elif self.n==1:
|
|
return 0.0
|
|
else:
|
|
return None
|
|
|
|
def __str__(self):
|
|
if self.n==0:
|
|
return "no value"
|
|
elif self.n==1:
|
|
return str(self.y)
|
|
else:
|
|
return "%g +/- %g (%d accumulations)"%(self.mean(),self.mean_error(),self.n)
|
|
|
|
def __repr__(self):
|
|
return str(self)
|
|
|
|
class MeasurementResult(Drawable.Drawable, collections.UserDict):
|
|
|
|
def __init__(self, quantity_name):
|
|
"""
|
|
convenient accumulation and interface to plot functions
|
|
|
|
dictionary must not contain anything but AccumulatedValue instances
|
|
"""
|
|
Drawable.Drawable.__init__(self)
|
|
collections.UserDict.__init__(self)
|
|
self.quantity_name=quantity_name
|
|
self.lock=threading.RLock()
|
|
|
|
# get the selected item, if it does not exist, create an empty one
|
|
def __getitem__(self, key):
|
|
if key not in self:
|
|
a=AccumulatedValue()
|
|
self.data[float(key)]=a
|
|
return a
|
|
else:
|
|
return self.data[float(key)]
|
|
|
|
def __setitem__(self,key,value):
|
|
if not isinstance(value, AccumulatedValue):
|
|
value=AccumulatedValue(float(value))
|
|
return collections.UserDict.__setitem__(self,
|
|
float(key),
|
|
value)
|
|
|
|
def __add__(self, right_value):
|
|
if right_value==0:
|
|
return self.copy()
|
|
else:
|
|
raise Exception("not implemented")
|
|
|
|
def get_title(self):
|
|
return self.quantity_name
|
|
|
|
def get_xdata(self):
|
|
"""
|
|
sorted array of all dictionary entries without Accumulated Value objects with n==0
|
|
"""
|
|
keys=numpy.array([k for k in list(self.data.keys()) if not (isinstance(self.data[k], AccumulatedValue) and self.data[k].n==0)],
|
|
dtype="float64")
|
|
keys.sort()
|
|
return keys
|
|
|
|
def get_ydata(self):
|
|
return self.get_xydata()[1]
|
|
|
|
def get_xydata(self):
|
|
k=self.get_xdata()
|
|
v=numpy.array([self.data[key].mean() for key in k], dtype="float64")
|
|
return [k,v]
|
|
|
|
def get_errorplotdata(self):
|
|
k=self.get_xdata()
|
|
v=numpy.array([self.data[key].mean() for key in k], dtype="float64")
|
|
e=numpy.array([self.data[key].mean_error() for key in k], dtype="float64")
|
|
return [k,v,e]
|
|
|
|
def get_lineplotdata(self):
|
|
k=self.get_xdata()
|
|
v=numpy.array(self.y, dtype="float64")
|
|
return [k, v]
|
|
|
|
def uses_statistics(self):
|
|
"""
|
|
drawable interface method, returns True
|
|
"""
|
|
return True
|
|
|
|
def write_to_csv(self,destination=sys.stdout, delimiter=" "):
|
|
"""
|
|
writes the data to a file or to sys.stdout
|
|
destination can be a file or a filename
|
|
suitable for further processing
|
|
"""
|
|
# write sorted
|
|
the_destination=destination
|
|
if type(destination) in (str,):
|
|
the_destination=open(destination, "w")
|
|
|
|
the_destination.write("# quantity:"+str(self.quantity_name)+"\n")
|
|
the_destination.write("# x y ysigma n\n")
|
|
for x in self.get_xdata():
|
|
y=self.data[x]
|
|
if type(y) in [float, int, int]:
|
|
the_destination.write("%e%s%e%s0%s1\n"%(x, delimiter, y, delimiter, delimiter))
|
|
else:
|
|
the_destination.write("%e%s%e%s%e%s%d\n"%(x,
|
|
delimiter,
|
|
y.mean(),
|
|
delimiter,
|
|
y.mean_error(),
|
|
delimiter,
|
|
y.n))
|
|
the_destination=None
|
|
|
|
|
|
def write_to_hdf(self, hdffile, where, name, title, complib=None, complevel=None):
|
|
|
|
h5_table_format= {
|
|
"x" : tables.Float64Col(),
|
|
"y" : tables.Float64Col(),
|
|
"y_err" : tables.Float64Col(),
|
|
"n" : tables.Int64Col()
|
|
}
|
|
filter=None
|
|
if complib is not None:
|
|
if complevel is None:
|
|
complevel=9
|
|
filter=tables.Filters(complevel=complevel,complib=complib,shuffle=1)
|
|
|
|
mr_table=hdffile.create_table(where=where,name=name,
|
|
description=h5_table_format,
|
|
title=title,
|
|
filters=filter,
|
|
expectedrows=len(self))
|
|
mr_table.flavor="numpy"
|
|
mr_table.attrs.damaris_type="MeasurementResult"
|
|
self.lock.acquire()
|
|
try:
|
|
mr_table.attrs.quantity_name=self.quantity_name
|
|
row=mr_table.row
|
|
xdata=self.get_xdata()
|
|
if xdata.shape[0]!=0:
|
|
for x in self.get_xdata():
|
|
y=self.data[x]
|
|
row["x"]=x
|
|
if type(y) in [float, int, int]:
|
|
row["y"]=y
|
|
row["y_err"]=0.0
|
|
row["n"]=1
|
|
else:
|
|
row["y"]=y.mean()
|
|
row["y_err"]=y.mean_error()
|
|
row["n"]=y.n
|
|
row.append()
|
|
|
|
finally:
|
|
mr_table.flush()
|
|
self.lock.release()
|
|
|
|
def read_from_hdf(hdf_node):
|
|
"""
|
|
reads a MeasurementResult object from the hdf_node
|
|
or None if the node is not suitable
|
|
"""
|
|
|
|
if not isinstance(hdf_node, tables.Table):
|
|
return None
|
|
|
|
if hdf_node._v_attrs.damaris_type!="MeasurementResult":
|
|
return None
|
|
|
|
mr=MeasurementResult(hdf_node._v_attrs.quantity_name)
|
|
for r in hdf_node.iterrows():
|
|
mr[r["x"]]=AccumulatedValue(r["y"],r["y_err"],r["n"])
|
|
|
|
return mr
|