import threading
import math
import sys
import tables
import numpy
import collections
from typing import IO
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
# internally it keeps the sum, the sum of squares and the number of data points
class AccumulatedValue:
def __init__(self, mean: float | None = None, mean_err: float | None = None, n: int | None = 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: "AccumulatedValue | float") -> "AccumulatedValue":
new_one=AccumulatedValue()
if 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: "AccumulatedValue | float") -> "AccumulatedValue":
if 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) -> "AccumulatedValue":
a=AccumulatedValue()
a.y=self.y
a.y2=self.y2
a.n=self.n
return a
def mean(self) -> float | None:
"""
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) -> float | None:
"""
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) -> float | None:
"""
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) -> str:
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) -> str:
return str(self)
[docs]
class MeasurementResult(Drawable.Drawable, collections.UserDict):
[docs]
def __init__(self, quantity_name: str) -> None:
"""
convenient accumulation and interface to plot functions
The 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: float) -> "AccumulatedValue":
if key not in self:
a=AccumulatedValue()
self.data[float(key)]=a
return a
else:
return self.data[float(key)]
def __setitem__(self,key: float,item: "AccumulatedValue | float")-> None:
if not isinstance(item, AccumulatedValue):
item=AccumulatedValue(float(item))
return collections.UserDict.__setitem__(self,
float(key),
item)
def __add__(self, right_value: int) -> "MeasurementResult":
if right_value==0:
new_one=MeasurementResult(self.quantity_name)
for key in self.data:
new_one[key]=self.data[key].copy()
return new_one
else:
raise Exception("not implemented")
def get_title(self) -> str:
return self.quantity_name
def get_xdata(self) -> numpy.ndarray:
"""
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, channel: int = 0) -> numpy.ndarray:
"""
y data (mean values); the channel argument is accepted for
compatibility with the Drawable interface, MeasurementResult
is single-channel
"""
return self.get_xydata()[1]
def get_xydata(self) -> list[numpy.ndarray]:
k=self.get_xdata()
v=numpy.array([self.data[key].mean() for key in k], dtype="float64")
return [k,v]
def get_errorplotdata(self) -> list[numpy.ndarray]:
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) -> list[numpy.ndarray]:
return self.get_xydata()
def uses_statistics(self) -> bool:
"""
drawable interface method returns True
"""
return True
def write_to_csv(self,destination: str | IO[str] = sys.stdout, delimiter: str = " ") -> None:
"""
writes the data to a file or to sys.stdout
destination can be a file or a filename
suitable for further processing
"""
if isinstance(destination, str):
with open(destination, "w") as f:
self._write_to_csv_dest(f, delimiter)
else:
self._write_to_csv_dest(destination, delimiter)
def _write_to_csv_dest(self, dest: IO[str], delimiter: str) -> None:
dest.write("# quantity:"+str(self.quantity_name)+"\n")
dest.write("# x y ysigma n\n")
for x in self.get_xdata():
y = self.data[x]
dest.write("%e%s%e%s%e%s%d\n" % (x,
delimiter,
y.mean(),
delimiter,
y.mean_error(),
delimiter,
y.n))
def write_to_hdf(self, hdffile: tables.File, where: str, name: str, title: str, complib: str | None = None, complevel: int | None = None) -> None:
h5_table_format= {
"x" : tables.Float32Col(),
"y" : tables.Float32Col(),
"y_err" : tables.Float32Col(),
"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 xdata:
y=self.data[x]
row["x"]=x
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: object) -> "MeasurementResult | None":
"""
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