Markus Rosenstihl
28e6c9d1a0
* Changed default number types to int16 for ADC_Result and float32 for Accumulations
971 lines
44 KiB
Python
971 lines
44 KiB
Python
# -*- coding: iso-8859-1 -*-
|
|
|
|
#############################################################################
|
|
# #
|
|
# Name: Class Accumulation #
|
|
# #
|
|
# Purpose: Specialised class of Errorable and Drawable #
|
|
# Contains accumulated ADC-Data #
|
|
# #
|
|
#############################################################################
|
|
|
|
from .Errorable import Errorable
|
|
from .Drawable import Drawable
|
|
from .DamarisFFT import DamarisFFT
|
|
from .Signalpath import Signalpath
|
|
|
|
#from DataPool import DataPool
|
|
|
|
import sys
|
|
import threading
|
|
import types
|
|
import tables
|
|
import numpy
|
|
import datetime # added by Oleg Petrov
|
|
import ctypes # added by Oleg Petrov
|
|
import struct # added by Oleg Petrov
|
|
import os # added by Oleg Petrov
|
|
|
|
class Accumulation(Errorable, Drawable, DamarisFFT, Signalpath):
|
|
def __init__(self, x = None, y = None, y_2 = None, n = None, index = None, sampl_freq = None, error = False):
|
|
Errorable.__init__(self)
|
|
Drawable.__init__(self)
|
|
|
|
# Title of this accumulation (plotted in GUI -> look Drawable)
|
|
self.__title_pattern = "Accumulation: n = %d"
|
|
|
|
# Axis-Labels (inherited from Drawable)
|
|
self.xlabel = "Time (s)"
|
|
self.ylabel = "Avg. Samples [Digits]"
|
|
self.lock=threading.RLock()
|
|
|
|
self.common_descriptions=None
|
|
self.time_period=[]
|
|
self.job_id = None # added by Oleg Petrov
|
|
|
|
self.use_error = error
|
|
|
|
if self.uses_statistics():
|
|
if (y_2 is not None):
|
|
self.y_square = y_2
|
|
elif (y_2 is None) :
|
|
self.y_square = []
|
|
else:
|
|
raise ValueError("Wrong usage of __init__!")
|
|
|
|
|
|
if (x is None) and (y is None) and (index is None) and (sampl_freq is None) and (n is None):
|
|
self.sampling_rate = 0
|
|
|
|
self.n = 0
|
|
self.set_title(self.__title_pattern % self.n)
|
|
|
|
self.cont_data = False
|
|
self.index = []
|
|
self.x = []
|
|
self.y = []
|
|
|
|
elif (x is not None) and (y is not None) and (index is not None) and (sampl_freq is not None) and (n is not None):
|
|
self.x = x
|
|
self.y = y
|
|
self.sampling_rate = sampl_freq
|
|
|
|
self.n = n
|
|
self.set_title(self.__title_pattern % self.n)
|
|
|
|
self.index = index
|
|
self.cont_data = True
|
|
|
|
else:
|
|
raise ValueError("Wrong usage of __init__!")
|
|
|
|
|
|
def get_accu_by_index(self, index):
|
|
self.lock.acquire()
|
|
try:
|
|
start = self.index[index][0]
|
|
end = self.index[index][1]
|
|
except:
|
|
self.lock.release()
|
|
raise
|
|
|
|
tmp_x = self.x[start:end+1]
|
|
tmp_y = []
|
|
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_y.append(self.y[i][start:end+1])
|
|
|
|
r = Accumulation(x = tmp_x, y = tmp_y, n = self.n, index = [(0,len(tmp_y[0])-1)], sampl_freq = self.sampling_rate, error = self.use_error)
|
|
self.lock.release()
|
|
return r
|
|
|
|
def get_ysquare(self, channel):
|
|
if self.uses_statistics():
|
|
try:
|
|
return self.y_square[channel]
|
|
except:
|
|
raise
|
|
else: return None
|
|
|
|
def contains_data(self):
|
|
return self.cont_data
|
|
|
|
|
|
def get_sampling_rate(self):
|
|
"Returns the samplingfrequency"
|
|
return self.sampling_rate + 0
|
|
|
|
|
|
def get_index_bounds(self, index):
|
|
"Returns a tuple with (start, end) of the wanted result"
|
|
return self.index[index]
|
|
|
|
|
|
def uses_statistics(self):
|
|
return self.use_error
|
|
|
|
# Schnittstellen nach Außen --------------------------------------------------------------------
|
|
|
|
def get_yerr(self, channel):
|
|
"""
|
|
return error (std.dev/sqrt(n)) of mean
|
|
"""
|
|
|
|
if not self.uses_statistics(): return numpy.zeros((len(self.y[0]),),dtype="float32")
|
|
if not self.contains_data(): return []
|
|
|
|
self.lock.acquire()
|
|
if self.n < 2:
|
|
retval=numpy.zeros((len(self.y[0]),),dtype="float32")
|
|
self.lock.release()
|
|
return retval
|
|
try:
|
|
variance_over_n = (self.y_square[channel] - (self.y[channel]**2 / float(self.n)))/float((self.n-1)*self.n)
|
|
except IndexError:
|
|
print("Warning Accumulation.get_ydata(channel): Channel index does not exist.")
|
|
variance_over_n = numpy.zeros((len(self.y[0]),), dtype="float32")
|
|
self.lock.release()
|
|
# sample standard deviation / sqrt(n)
|
|
return numpy.nan_to_num(numpy.sqrt(variance_over_n))
|
|
|
|
def get_ydata(self, channel):
|
|
"""
|
|
return mean data
|
|
"""
|
|
|
|
if not self.contains_data(): return []
|
|
self.lock.acquire()
|
|
|
|
try:
|
|
tmp_y = self.y[channel] / self.n
|
|
except IndexError:
|
|
print("Warning Accumulation.get_ydata(channel): Channel index does not exist.")
|
|
tmp_y = numpy.zeros((len(self.y[0]),), dtype="float32")
|
|
|
|
self.lock.release()
|
|
return tmp_y
|
|
|
|
|
|
def get_ymin(self):
|
|
|
|
if not self.contains_data(): return 0
|
|
|
|
tmp_min = []
|
|
self.lock.acquire()
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_min.append(self.get_ydata(i).min())
|
|
|
|
if self.uses_statistics() and self.ready_for_drawing_error():
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_min.append((self.get_ydata(i) - self.get_yerr(i)).min())
|
|
self.lock.release()
|
|
|
|
return min(tmp_min)
|
|
|
|
|
|
def get_ymax(self):
|
|
|
|
if not self.contains_data(): return 0
|
|
|
|
tmp_max = []
|
|
self.lock.acquire()
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_max.append(self.get_ydata(i).max())
|
|
|
|
if self.uses_statistics() and self.ready_for_drawing_error():
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_max.append((self.get_ydata(i) + self.get_yerr(i)).max())
|
|
self.lock.release()
|
|
return max(tmp_max)
|
|
|
|
def get_job_id(self):
|
|
# return None
|
|
return self.job_id # modified by Oleg Petrov
|
|
|
|
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
|
|
"""
|
|
|
|
the_destination=destination
|
|
if type(destination) in (str,):
|
|
the_destination=open(destination, "w")
|
|
|
|
the_destination.write("# accumulation %d\n"%self.n)
|
|
self.lock.acquire()
|
|
try:
|
|
if self.common_descriptions is not None:
|
|
for (key,value) in self.common_descriptions.items():
|
|
the_destination.write("# %s : %s\n"%(key, str(value)))
|
|
the_destination.write("# t")
|
|
ch_no=self.get_number_of_channels()
|
|
if self.use_error:
|
|
for i in range(ch_no): the_destination.write(" ch%d_mean ch%d_err"%(i,i))
|
|
else:
|
|
for i in range(ch_no): the_destination.write(" ch%d_mean"%i)
|
|
the_destination.write("\n")
|
|
xdata=self.get_xdata()
|
|
ydata=list(map(self.get_ydata, range(ch_no)))
|
|
yerr=None
|
|
if self.use_error:
|
|
yerr=list(map(self.get_yerr, range(ch_no)))
|
|
for i in range(len(xdata)):
|
|
the_destination.write("%e"%xdata[i])
|
|
for j in range(ch_no):
|
|
if self.use_error:
|
|
the_destination.write("%s%e%s%e"%(delimiter, ydata[j][i], delimiter, yerr[j][i]))
|
|
else:
|
|
the_destination.write("%s%e"%(delimiter,ydata[j][i]))
|
|
the_destination.write("\n")
|
|
the_destination=None
|
|
xdata=yerr=ydata=None
|
|
finally:
|
|
self.lock.release()
|
|
|
|
# ------------- added by Oleg Petrov, 14 Feb 2012 ----------------------
|
|
def write_to_simpson(self, destination=sys.stdout, delimiter=" ", frequency=100e6):
|
|
"""
|
|
writes the data to a text file or sys.stdout in Simpson format,
|
|
for further processing with the NMRnotebook software;
|
|
destination can be a file or a filename
|
|
"""
|
|
# write sorted
|
|
the_destination=destination
|
|
if type(destination) in (str,):
|
|
the_destination=open(destination, "w")
|
|
|
|
self.lock.acquire()
|
|
try:
|
|
xdata=self.get_xdata()
|
|
the_destination.write("SIMP\n")
|
|
the_destination.write("%s%i%s"%("NP=", len(xdata), "\n"))
|
|
the_destination.write("%s%i%s"%("SW=", self.get_sampling_rate(), "\n"))
|
|
the_destination.write("%s%i%s"%("REF=", frequency, "\n"))
|
|
the_destination.write("TYPE=FID\n")
|
|
the_destination.write("DATA\n")
|
|
ch_no=self.get_number_of_channels()
|
|
ydata=list(map(self.get_ydata, range(ch_no)))
|
|
for i in range(len(xdata)):
|
|
for j in range(ch_no):
|
|
the_destination.write("%g%s"%(ydata[j][i], delimiter))
|
|
the_destination.write("\n")
|
|
the_destination.write("END\n")
|
|
the_destination=None
|
|
xdata=ydata=None
|
|
finally:
|
|
self.lock.release()
|
|
|
|
# ------------- added by Oleg Petrov, 10 Sep 2013 -----------------------
|
|
def write_to_tecmag(self, destination=sys.stdout, nrecords=1,\
|
|
frequency=100.,\
|
|
last_delay = 1.,\
|
|
receiver_phase=0.,\
|
|
nucleus='1H'):
|
|
"""
|
|
writes the data to a binary file in TecMag format;
|
|
destination can be a file object or a filename;
|
|
nrecords determines an indirect dimension in 2D experiments;
|
|
"""
|
|
|
|
#TODO: Function is most likely broken in Python 3 because Strings are now unicode and binary files cannot write unicode
|
|
|
|
if self.job_id == None or self.n == 0:
|
|
raise ValueError("write_to_tecmag: cannot get a record number")
|
|
else:
|
|
record = (self.job_id/self.n)%nrecords + 1
|
|
|
|
the_destination=destination
|
|
if type(destination) in (str,):
|
|
if record == 1 and os.path.exists(destination):
|
|
os.rename(destination, os.path.dirname(destination)+'/~'+os.path.basename(destination))
|
|
|
|
self.lock.acquire()
|
|
try:
|
|
npts = [len(self), nrecords, 1, 1]
|
|
data_offset = 2*4*npts[0]*npts[1]*npts[2]*npts[3] # length of data section
|
|
dwell = 1./self.get_sampling_rate()
|
|
sw = 0.5/dwell
|
|
|
|
base_freq = [frequency, frequency, 0., 0.]
|
|
offset_freq = [0., 0., 0., 0.]
|
|
ob_freq = [sum(x) for x in zip(base_freq, offset_freq)]
|
|
date = self.time_period[0].strftime("%Y/%m/%d %H:%M:%S")
|
|
|
|
# data handling:
|
|
ch_no=self.get_number_of_channels()
|
|
|
|
ydata = list(map(self.get_ydata, range(ch_no)))
|
|
if ch_no == 1:
|
|
ydata = [ydata, numpy.zeros(len(ydata))]
|
|
|
|
# data is arranged in RIRIRIRI blocks in linear order:
|
|
data = numpy.append([ydata[0]], [ydata[1]], axis=0)
|
|
data = data.T
|
|
data = data.flatten()
|
|
|
|
if record == 1:
|
|
the_destination=open(destination, "wb")
|
|
|
|
# allocate space for all records in advance:
|
|
buff = ctypes.create_string_buffer(1056+data_offset+2068)
|
|
|
|
struct.pack_into('8s', buff, 0, 'TNT1.005') # 'TNT1.000' version ID
|
|
struct.pack_into('4s', buff, 8, 'TMAG') # 'TMAG' tag
|
|
struct.pack_into('?', buff, 12, True) # BOOLean value
|
|
struct.pack_into('i', buff, 16, 1024) # length of Tecmag struct
|
|
|
|
#Initialize TECMAG structure:
|
|
struct.pack_into('4i', buff, 20, *npts) # npts[4]
|
|
struct.pack_into('4i', buff, 36, *npts) # actual_npts[4]
|
|
struct.pack_into('i', buff, 52, npts[0]) # acq_points
|
|
struct.pack_into('4i', buff, 56, 1, 1, 1, 1) # npts_start[4]
|
|
struct.pack_into('i', buff, 72, self.n) # scans
|
|
struct.pack_into('i', buff, 76, self.n) # actual_scans
|
|
struct.pack_into('i', buff, 88, 1) # sadimension
|
|
|
|
struct.pack_into('4d', buff, 104, *ob_freq) # ob_freq[4]
|
|
struct.pack_into('4d', buff, 136, *base_freq) # base_freq[4]
|
|
struct.pack_into('4d', buff, 168, *offset_freq) # offset_freq[4]
|
|
struct.pack_into('d', buff, 200, 0.0) # ref_freq
|
|
struct.pack_into('h', buff, 216, 1) # obs_channel
|
|
struct.pack_into('42s', buff, 218, 42*'2') # space2[42]
|
|
|
|
struct.pack_into('4d', buff, 260, sw, sw, 0., 0.) # sw[4], sw = 0.5/dwell
|
|
struct.pack_into('4d', buff, 292, dwell, dwell, 0., 0.) # dwell[4]
|
|
struct.pack_into('d', buff, 324, sw) # filter, = 0.5/dwell
|
|
struct.pack_into('d', buff, 340, (npts[0]*dwell)) # acq_time
|
|
struct.pack_into('d', buff, 348, 1.) # last_delay (5*T1 minus sequence length)
|
|
struct.pack_into('h', buff, 356, 1) # spectrum_direction
|
|
struct.pack_into('16s', buff, 372, 16*'2') # space3[16]
|
|
|
|
struct.pack_into('d', buff, 396, receiver_phase) # receiver_phase
|
|
struct.pack_into('4s', buff, 404, 4*'2') # space4[4]
|
|
struct.pack_into('16s', buff, 444, 16*'2') # space5[16]
|
|
struct.pack_into('264s', buff, 608, 264*'2') # space6[264]
|
|
|
|
struct.pack_into('32s', buff, 884, date) # date[32]
|
|
struct.pack_into('16s', buff, 916, nucleus) # nucleus[16]
|
|
# TECMAG Structure total => 1024
|
|
|
|
struct.pack_into('4s', buff, 1044, 'DATA') # 'DATA' tag
|
|
struct.pack_into('?', buff, 1048, True) # BOOLean
|
|
struct.pack_into('i', buff, 1052, data_offset) # length of data
|
|
|
|
struct.pack_into('%sf' % (2*npts[0]), buff, 1056, *data) # actual data (one record)
|
|
|
|
struct.pack_into('4s', buff, 1056+data_offset, 'TMG2') # 'TMG2' tag
|
|
struct.pack_into('?', buff, 1056+data_offset+4, True) # BOOLean
|
|
struct.pack_into('i', buff, 1056+data_offset+8, 2048) # length of Tecmag2 struct
|
|
|
|
# Leave TECMAG2 structure empty:
|
|
struct.pack_into('52s', buff, 1056+data_offset+372, 52*'2') # space[52]
|
|
struct.pack_into('866s', buff, 1056+data_offset+1194, 866*'2') # space[610]+names+strings
|
|
# TECMAG2 Structure total => 2048
|
|
|
|
struct.pack_into('4s', buff, 1056+data_offset+2060, 'PSEQ') # 'PSEQ' tag 658476
|
|
struct.pack_into('?', buff, 1056+data_offset+2064, False) # BOOLean 658480
|
|
|
|
the_destination.write(buff)
|
|
|
|
else:
|
|
the_destination=open(destination, "rb+")
|
|
buff = ctypes.create_string_buffer(4*2*npts[0])
|
|
struct.pack_into('%sf' % (2*npts[0]), buff, 0, *data)
|
|
|
|
the_destination.seek(1056+4*2*npts[0]*(record-1))
|
|
the_destination.write(buff)
|
|
|
|
the_destination = None
|
|
ydata=None
|
|
|
|
finally:
|
|
self.lock.release()
|
|
|
|
# -----------------------------------------------------------------------
|
|
|
|
def write_to_hdf(self, hdffile, where, name, title, complib=None, complevel=None):
|
|
accu_group=hdffile.create_group(where=where,name=name,title=title)
|
|
accu_group._v_attrs.damaris_type="Accumulation"
|
|
if self.contains_data():
|
|
self.lock.acquire()
|
|
try:
|
|
# save time stamps
|
|
if self.time_period is not None and len(self.time_period)>0:
|
|
accu_group._v_attrs.earliest_time="%04d%02d%02d %02d:%02d:%02d.%03d"%(self.time_period[0].year,
|
|
self.time_period[0].month,
|
|
self.time_period[0].day,
|
|
self.time_period[0].hour,
|
|
self.time_period[0].minute,
|
|
self.time_period[0].second,
|
|
self.time_period[0].microsecond/1000)
|
|
accu_group._v_attrs.oldest_time="%04d%02d%02d %02d:%02d:%02d.%03d"%(self.time_period[1].year,
|
|
self.time_period[1].month,
|
|
self.time_period[1].day,
|
|
self.time_period[1].hour,
|
|
self.time_period[1].minute,
|
|
self.time_period[1].second,
|
|
self.time_period[1].microsecond/1000)
|
|
if self.common_descriptions is not None:
|
|
for (key,value) in self.common_descriptions.items():
|
|
accu_group._v_attrs.__setattr__("description_"+key,str(value))
|
|
accu_group._v_attrs.__setattr__("sampling_rate",self.sampling_rate)
|
|
# save interval information
|
|
filter=None
|
|
if complib is not None:
|
|
if complevel is None:
|
|
complevel=9
|
|
filter=tables.Filters(complevel=complevel,complib=complib,shuffle=1)
|
|
|
|
# tried compression filter, but no effect...
|
|
index_table=hdffile.create_table(where=accu_group,
|
|
name="indices",
|
|
description={"start": tables.UInt64Col(),
|
|
"length": tables.UInt64Col(),
|
|
"start_time": tables.Float64Col(),
|
|
"dwelltime": tables.Float64Col(),
|
|
"number": tables.UInt64Col()},
|
|
title="indices of adc data intervals",
|
|
filters=filter,
|
|
expectedrows=len(self.index))
|
|
index_table.flavor="numpy"
|
|
# save interval data
|
|
new_row=index_table.row
|
|
for i in range(len(self.index)):
|
|
new_row["start"]=self.index[i][0]
|
|
new_row["dwelltime"]=1.0/self.sampling_rate
|
|
new_row["start_time"]=1.0/self.sampling_rate*self.index[i][0]
|
|
new_row["length"]=self.index[i][1]-self.index[i][0]+1
|
|
new_row["number"]=self.n
|
|
new_row.append()
|
|
|
|
index_table.flush()
|
|
new_row=None
|
|
index_table=None
|
|
|
|
# prepare saving data
|
|
channel_no=len(self.y)
|
|
timedata=numpy.empty((len(self.y[0]),channel_no*2), dtype = "float32")
|
|
for ch in range(channel_no):
|
|
timedata[:,ch*2]=self.get_ydata(ch)
|
|
if self.uses_statistics():
|
|
timedata[:,ch*2+1]=self.get_yerr(ch)
|
|
else:
|
|
timedata[:,ch*2+1]=numpy.zeros((len(self.y[0]),),dtype = "float32")
|
|
|
|
# save data
|
|
time_slice_data=None
|
|
if filter is not None:
|
|
chunkshape=timedata.shape
|
|
if len(chunkshape) <= 1:
|
|
chunkshape = (min(chunkshape[0],1024*8),)
|
|
else:
|
|
chunkshape = (min(chunkshape[0],1024*8), chunkshape[1])
|
|
if tables.__version__[0]=="1":
|
|
time_slice_data=hdffile.create_carray(accu_group,
|
|
name="accu_data",
|
|
shape=timedata.shape,
|
|
atom=tables.Float32Atom(shape=chunkshape,
|
|
flavor="numpy"),
|
|
filters=filter,
|
|
title="accu data")
|
|
else:
|
|
time_slice_data=hdffile.create_carray(accu_group,
|
|
name="accu_data",
|
|
shape=timedata.shape,
|
|
chunkshape=chunkshape,
|
|
atom=tables.Float32Atom(),
|
|
filters=filter,
|
|
title="accu data")
|
|
|
|
time_slice_data[:]=timedata
|
|
else:
|
|
time_slice_data=hdffile.create_array(accu_group,
|
|
name="accu_data",
|
|
obj=timedata,
|
|
title="accu data")
|
|
|
|
|
|
|
|
finally:
|
|
time_slice_data=None
|
|
accu_group=None
|
|
self.lock.release()
|
|
|
|
# / Schnittstellen nach Außen ------------------------------------------------------------------
|
|
|
|
# Überladen von Operatoren ---------------------------------------------------------------------
|
|
|
|
def __len__(self):
|
|
"""
|
|
return number of samples per channel, 0 if empty
|
|
"""
|
|
if len(self.y)>0:
|
|
return len(self.y[0])
|
|
return 0
|
|
|
|
def __repr__(self):
|
|
"Redefining repr(Accumulation)"
|
|
|
|
if not self.contains_data(): return "Empty"
|
|
|
|
tmp_string = "X: " + repr(self.x) + "\n"
|
|
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_string += ("Y(%d): " % i) + repr(self.y[i]) + "\n"
|
|
if self.uses_statistics(): tmp_string += "y_square(%d): " % i + str(self.y_square[i]) + "\n"
|
|
|
|
tmp_string += "Indexes: " + str(self.index) + "\n"
|
|
|
|
tmp_string += "Samples per Channel: " + str(len(self.y[0])) + "\n"
|
|
tmp_string += "Samplingfrequency: " + str(self.sampling_rate) + "\n"
|
|
|
|
tmp_string += "n: " + str(self.n)
|
|
|
|
return tmp_string
|
|
|
|
|
|
def __add__(self, other):
|
|
"Redefining self + other"
|
|
# Float or int
|
|
if isinstance(other, int) or isinstance(other, float):
|
|
if not self.contains_data(): raise ValueError("Accumulation: You cant add integers/floats to an empty accumulation")
|
|
else:
|
|
|
|
tmp_y = []
|
|
tmp_ysquare = []
|
|
|
|
|
|
self.lock.acquire()
|
|
for i in range(self.get_number_of_channels()):
|
|
# Dont change errors and mean value
|
|
if self.uses_statistics(): tmp_ysquare.append(self.y_square[i] + ( (2*self.y[i]*other) + ((other**2)*self.n) ))
|
|
tmp_y.append(self.y[i] + (other*self.n))
|
|
|
|
if self.uses_statistics():
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, y_2 = tmp_ysquare, n = self.n, index = self.index, sampl_freq = self.sampling_rate, error = self.use_error)
|
|
else:
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, n = self.n, index = self.index, sampl_freq = self.sampling_rate, error = self.use_error)
|
|
r.job_id = self.job_id # added by Oleg Petrov
|
|
self.lock.release()
|
|
return r
|
|
|
|
# ADC_Result
|
|
elif str(other.__class__) == "damaris.data.ADC_Result.ADC_Result":
|
|
|
|
# Other empty (return)
|
|
# todo: this is seems to be bugy!!!! (Achim)
|
|
if not other.contains_data(): return
|
|
|
|
# Self empty (copy)
|
|
if not self.contains_data():
|
|
|
|
tmp_y = []
|
|
tmp_ysquare = []
|
|
|
|
self.lock.acquire()
|
|
|
|
for i in range(other.get_number_of_channels()):
|
|
tmp_y.append(numpy.array(other.y[i], dtype="float32"))
|
|
if self.uses_statistics(): tmp_ysquare.append(tmp_y[i] ** 2)
|
|
|
|
|
|
if self.uses_statistics():
|
|
r = Accumulation(x = numpy.array(other.x, dtype="float32"), y = tmp_y, y_2 = tmp_ysquare, n = 1, index = other.index, sampl_freq = other.sampling_rate, error = True)
|
|
else:
|
|
r = Accumulation(x = numpy.array(other.x, dtype="float32"), y = tmp_y, index = other.index, sampl_freq = other.sampling_rate, n = 1, error = False)
|
|
r.time_period=[other.job_date,other.job_date]
|
|
r.job_id = other.job_id # added by Oleg Petrov
|
|
r.common_descriptions=other.description.copy()
|
|
self.lock.release()
|
|
return r
|
|
|
|
# Other and self not empty (self + other)
|
|
else:
|
|
self.lock.acquire()
|
|
|
|
if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add ADC-Results with diffrent sampling-rates")
|
|
if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add ADC-Results with diffrent number of samples")
|
|
if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add ADC-Results with diffrent number of channels")
|
|
for i in range(len(self.index)):
|
|
if self.index[i] != other.get_index_bounds(i): raise ValueError("Accumulation: You cant add ADC-Results with diffrent indexing")
|
|
|
|
tmp_y = []
|
|
tmp_ysquare = []
|
|
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_y.append(self.y[i] + other.y[i])
|
|
if self.uses_statistics(): tmp_ysquare.append(self.y_square[i] + (numpy.array(other.y[i], dtype="float32") ** 2))
|
|
|
|
if self.uses_statistics():
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, y_2 = tmp_ysquare, n = self.n + 1, index = self.index, sampl_freq = self.sampling_rate, error = True)
|
|
else:
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, n = self.n + 1, index = self.index, sampl_freq = self.sampling_rate, error = False)
|
|
r.time_period=[min(self.time_period[0],other.job_date),
|
|
max(self.time_period[1],other.job_date)]
|
|
r.job_id = other.job_id # added by Oleg Petrov
|
|
if self.common_descriptions is not None:
|
|
r.common_descriptions={}
|
|
for key in list(self.common_descriptions.keys()):
|
|
if (key in other.description and self.common_descriptions[key]==other.description[key]):
|
|
r.common_descriptions[key]=value
|
|
|
|
self.lock.release()
|
|
return r
|
|
|
|
# Accumulation
|
|
elif str(other.__class__) == "damaris.data.Accumulation.Accumulation":
|
|
|
|
# Other empty (return)
|
|
if not other.contains_data(): return
|
|
|
|
# Self empty (copy)
|
|
if not self.contains_data():
|
|
|
|
tmp_y = []
|
|
tmp_ysquare = []
|
|
|
|
self.lock.acquire()
|
|
|
|
if self.uses_statistics():
|
|
r = Accumulation(x = numpy.array(other.x, dtype="float32"), y = tmp_y, y_2 = tmp_ysquare, n = other.n, index = other.index, sampl_freq = other.sampling_rate, error = True)
|
|
else:
|
|
r = Accumulation(x = numpy.array(other.x, dtype="float32"), y = tmp_y, n = other.n, index = other.index, sampl_freq = other.sampling_rate, error = False)
|
|
for i in range(other.get_number_of_channels()):
|
|
tmp_y.append(other.y[i])
|
|
tmp_ysquare.append(other.y_square[i])
|
|
r.time_period=other.time_period[:]
|
|
r.job_id = other.job_id # added by Oleg Petrov
|
|
if other.common_descriptions is not None:
|
|
r.common_descriptions=othter.common_descriptions.copy()
|
|
else:
|
|
r.common_descriptions=None
|
|
|
|
self.lock.release()
|
|
return r
|
|
|
|
# Other and self not empty (self + other)
|
|
else:
|
|
self.lock.acquire()
|
|
|
|
if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add accumulations with diffrent sampling-rates")
|
|
if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add accumulations with diffrent number of samples")
|
|
if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add accumulations with diffrent number of channels")
|
|
for i in range(len(self.index)):
|
|
if self.index[i] != other.get_index_bounds(i): raise ValueError("Accumulation: You cant add accumulations with diffrent indexing")
|
|
if self.uses_statistics() and not other.uses_statistics(): raise ValueError("Accumulation: You cant add non-error accumulations to accumulations with error")
|
|
|
|
tmp_y = []
|
|
tmp_ysquare = []
|
|
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_y.append(self.y[i] + other.y[i])
|
|
tmp_ysquare.append(self.y_square[i] + other.y_square[i])
|
|
|
|
if self.uses_statistics():
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, y_2 = tmp_ysquare, n = other.n + self.n, index = self.index, sampl_freq = self.sampling_rate, error = True)
|
|
else:
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, n = other.n + self.n, index = self.index, sampl_freq = self.sampling_rate, error = False)
|
|
|
|
r.time_period=[min(self.time_period[0],other.time_period[0]),
|
|
max(self.time_period[1],other.time_period[1])]
|
|
r.job_id = other.job_id # added by Oleg Petrov
|
|
r.common_descriptions={}
|
|
if self.common_descriptions is not None and other.common_descriptions is not None:
|
|
for key in list(self.common_descriptions.keys()):
|
|
if (key in other.common_descriptions and
|
|
self.common_descriptions[key]==other.common_descriptions[key]):
|
|
r.common_descriptions[key]=value
|
|
|
|
self.lock.release()
|
|
return r
|
|
|
|
|
|
def __radd__(self, other):
|
|
"Redefining other + self"
|
|
return self.__add__(other)
|
|
|
|
|
|
def __sub__(self, other):
|
|
"Redefining self - other"
|
|
return self.__add__(-other)
|
|
|
|
|
|
def __rsub__(self, other):
|
|
"Redefining other - self"
|
|
return self.__neg__(self.__add__(-other))
|
|
|
|
|
|
def __iadd__(self, other):
|
|
"Redefining self += other"
|
|
# Float or int
|
|
if isinstance(other, int) or isinstance(other, float):
|
|
if not self.contains_data(): raise ValueError("Accumulation: You cant add integers/floats to an empty accumulation")
|
|
else:
|
|
|
|
self.lock.acquire()
|
|
for i in range(self.get_number_of_channels()):
|
|
#Dont change errors and mean value
|
|
if self.uses_statistics(): self.y_square[i] += (2*self.y[i]*other) + ((other**2)*self.n)
|
|
self.y[i] += other*self.n
|
|
self.lock.release()
|
|
|
|
return self
|
|
|
|
# ADC_Result
|
|
elif type(other).__module__+"."+type(other).__name__ == "damaris.data.ADC_Result.ADC_Result":
|
|
|
|
# Other empty (return)
|
|
if not other.contains_data(): return self
|
|
|
|
# Self empty (copy)
|
|
if not self.contains_data():
|
|
self.lock.acquire()
|
|
self.n += 1
|
|
self.index = other.index[0:]
|
|
self.sampling_rate = other.sampling_rate
|
|
self.x = numpy.array(other.x, dtype="float32")
|
|
self.cont_data = True
|
|
|
|
for i in range(other.get_number_of_channels()):
|
|
self.y.append(numpy.array(other.y[i], dtype="float32"))
|
|
if self.uses_statistics(): self.y_square.append(self.y[i] ** 2)
|
|
|
|
self.set_title(self.__title_pattern % self.n)
|
|
self.lock.release()
|
|
|
|
self.time_period=[other.job_date,other.job_date]
|
|
self.job_id = other.job_id # added by Oleg Petrov
|
|
self.common_descriptions=other.description.copy()
|
|
|
|
return self
|
|
|
|
|
|
# Other and self not empty (self + other)
|
|
else:
|
|
self.lock.acquire()
|
|
|
|
if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add ADC-Results with diffrent sampling-rates")
|
|
if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add ADC-Results with diffrent number of samples")
|
|
if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add ADC-Results with diffrent number of channels")
|
|
for i in range(len(self.index)):
|
|
if self.index[i] != other.get_index_bounds(i): raise ValueError("Accumulation: You cant add ADC-Results with diffrent indexing")
|
|
|
|
for i in range(self.get_number_of_channels()):
|
|
self.y[i] += other.y[i]
|
|
if self.uses_statistics(): self.y_square[i] += numpy.array(other.y[i], dtype="float32") ** 2
|
|
|
|
self.n += 1
|
|
self.time_period=[min(self.time_period[0],other.job_date),
|
|
max(self.time_period[1],other.job_date)]
|
|
self.job_id = other.job_id # added by Oleg Petrov
|
|
if self.common_descriptions is not None:
|
|
for key in list(self.common_descriptions.keys()):
|
|
if not (key in other.description and self.common_descriptions[key]==other.description[key]):
|
|
del self.common_descriptions[key]
|
|
|
|
self.set_title(self.__title_pattern % self.n)
|
|
self.lock.release()
|
|
|
|
return self
|
|
|
|
# Accumulation
|
|
elif type(other).__module__+"."+type(other).__name__ == "damaris.data.Accumulation.Accumulation":
|
|
|
|
# Other empty (return)
|
|
if not other.contains_data(): return
|
|
|
|
# Self empty (copy)
|
|
if not self.contains_data():
|
|
if self.uses_statistics() and not other.uses_statistics(): raise ValueError("Accumulation: You cant add non-error accumulations to accumulations with error")
|
|
|
|
self.lock.acquire()
|
|
self.n += other.n
|
|
self.index = other.index[0:]
|
|
self.sampling_rate = other.sampling_rate
|
|
self.x = numpy.array(other.x, dtype="float32")
|
|
self.cont_data = True
|
|
|
|
for i in range(other.get_number_of_channels()):
|
|
self.y.append(numpy.array(other.y[i], dtype="float32"))
|
|
if self.uses_statistics(): self.y_square.append(self.y[i] ** 2)
|
|
|
|
self.set_title(self.__title_pattern % self.n)
|
|
self.common_descriptions=other.common_desriptions.copy()
|
|
self.time_period=other.time_period[:]
|
|
self.job_id = other.job_id # added by Oleg Petrov
|
|
self.lock.release()
|
|
|
|
return self
|
|
|
|
# Other and self not empty (self + other)
|
|
else:
|
|
self.lock.acquire()
|
|
if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add accumulations with diffrent sampling-rates")
|
|
if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add accumulations with diffrent number of samples")
|
|
if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add accumulations with diffrent number of channels")
|
|
for i in range(len(self.index)):
|
|
if self.index[i] != other.get_index_bounds(i): raise ValueError("Accumulation: You cant add accumulations with diffrent indexing")
|
|
if self.uses_statistics() and not other.uses_statistics(): raise ValueError("Accumulation: You cant add non-error accumulations to accumulations with error")
|
|
|
|
for i in range(self.get_number_of_channels()):
|
|
self.y[i] += other.y[i]
|
|
if self.uses_statistics(): self.y_square[i] += other.y_square[i]
|
|
|
|
self.n += other.n
|
|
self.time_period=[min(self.time_period[0],other.time_period[0]),
|
|
max(self.time_period[1],other.time_period[1])]
|
|
self.job_id = other.job_id # added by Oleg Petrov
|
|
if self.common_descriptions is not None and other.common_descriptions is not None:
|
|
for key in list(self.common_descriptions.keys()):
|
|
if not (key in other.description and
|
|
self.common_descriptions[key]==other.common_descriptions[key]):
|
|
del self.common_descriptions[key]
|
|
|
|
self.set_title(self.__title_pattern % self.n)
|
|
self.lock.release()
|
|
|
|
return self
|
|
|
|
elif other is None:
|
|
# Convenience: ignore add of None
|
|
return self
|
|
else:
|
|
raise ValueError("can not add "+repr(type(other))+" to Accumulation")
|
|
|
|
|
|
def __isub__(self, other):
|
|
"Redefining self -= other"
|
|
return self.__iadd__(-other)
|
|
|
|
|
|
def __neg__(self):
|
|
"Redefining -self"
|
|
|
|
if not self.contains_data(): return
|
|
|
|
tmp_y = []
|
|
|
|
self.lock.acquire()
|
|
for i in range(self.get_number_of_channels()):
|
|
tmp_y.append(numpy.array(-self.y[i], dtype="float32"))
|
|
|
|
if self.uses_statistics():
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, y_2 = numpy.array(self.y_square), n = self.n, index = self.index, sampl_freq = self.sampling_rate, error = True)
|
|
else:
|
|
r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, n = self.n, index = self.index, sampl_freq = self.sampling_rate, error = False)
|
|
self.lock.release()
|
|
return r
|
|
|
|
|
|
def read_from_hdf(hdf_node):
|
|
"""
|
|
read accumulation data from HDF node and return it.
|
|
"""
|
|
|
|
# formal checks first
|
|
if not isinstance(hdf_node, tables.Group):
|
|
return None
|
|
|
|
if hdf_node._v_attrs.damaris_type!="Accumulation":
|
|
return None
|
|
|
|
if not (hdf_node.__contains__("indices") and hdf_node.__contains__("accu_data")):
|
|
print("no accu data")
|
|
return None
|
|
|
|
accu=Accumulation()
|
|
|
|
# populate description dictionary
|
|
accu.common_descriptions={}
|
|
for attrname in hdf_node._v_attrs._v_attrnamesuser:
|
|
if attrname.startswith("description_"):
|
|
accu.common_descriptions[attrname[12:]]=hdf_node._v_attrs.__getattr__(attrname)
|
|
|
|
eariliest_time=None
|
|
if "earliest_time" in dir(hdf_node._v_attrs):
|
|
timestring=hdf_node._v_attrs.__getattr__("earliest_time")
|
|
earliest_time=datetime.datetime(int(timestring[:4]), # year
|
|
int(timestring[4:6]), # month
|
|
int(timestring[6:8]), # day
|
|
int(timestring[9:11]), # hour
|
|
int(timestring[12:14]), # minute
|
|
int(timestring[15:17]), # second
|
|
int(timestring[18:21])*1000 # microsecond
|
|
)
|
|
|
|
oldest_time=None
|
|
if "oldest_time" in dir(hdf_node._v_attrs):
|
|
timestring=hdf_node._v_attrs.__getattr__("oldest_time")
|
|
oldest_time=datetime.datetime(int(timestring[:4]), # year
|
|
int(timestring[4:6]), # month
|
|
int(timestring[6:8]), # day
|
|
int(timestring[9:11]), # hour
|
|
int(timestring[12:14]), # minute
|
|
int(timestring[15:17]), # second
|
|
int(timestring[18:21])*1000 # microsecond
|
|
)
|
|
|
|
if oldest_time is None or earliest_time is None:
|
|
accu.time_period=None
|
|
if len(accu.common_descriptions)==0:
|
|
# no accus inside, so no common description expected
|
|
accu.common_descriptions=None
|
|
accu.cont_data=False
|
|
else:
|
|
accu.time_period=[oldest_time, earliest_time]
|
|
accu.cont_data=True
|
|
|
|
# start with indices
|
|
for r in hdf_node.indices.iterrows():
|
|
accu.index.append((r["start"],r["start"]+r["length"]-1))
|
|
accu.n=r["number"]
|
|
accu.sampling_rate=1.0/r["dwelltime"]
|
|
|
|
# now really belief there are no data
|
|
if len(accu.index)==0 or accu.n==0:
|
|
accu.cont_data=False
|
|
return accu
|
|
|
|
# now do the real data
|
|
accu_data=hdf_node.accu_data.read()
|
|
|
|
accu.x=numpy.arange(accu_data.shape[0], dtype="float32")/accu.sampling_rate
|
|
# assume error information, todo: save this information explicitly
|
|
accu.y_square=[]
|
|
accu.use_error=False
|
|
|
|
for ch in range(accu_data.shape[1]/2):
|
|
accu.y.append(accu_data[:,ch*2]*accu.n)
|
|
if accu.n<2 or numpy.all(accu_data[:,ch*2+1]==0.0):
|
|
accu.y_square.append(numpy.zeros((accu_data.shape[0]) ,dtype="float32"))
|
|
else:
|
|
accu.use_error=True
|
|
accu.y_square.append((accu_data[:,ch*2+1]**2)*float((accu.n-1.0)*accu.n)+(accu_data[:,ch*2]**2)*accu.n)
|
|
|
|
if not accu.use_error:
|
|
del accu.y_square
|
|
|
|
return accu
|