Source code for damaris.data.Accumulation

# -*- coding: iso-8859-1 -*-

from __future__ import annotations

#############################################################################
#                                                                           #
# 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 .ADC_Result import ADC_Result
from .MeasurementData import MeasurementData

import sys
import threading
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
from typing import IO, Any, Literal, cast

_XArray = numpy.ndarray | list
_YChannels = list[numpy.ndarray]
_IndexList = list[tuple[int, int]]
_Description = dict[str, str]
_CompressionLibrary = Literal["blosc", "blosc2", "bzip2", "lzo", "zlib"]


[docs] class Accumulation(MeasurementData, Errorable, Drawable, DamarisFFT, Signalpath): """ Represents an accumulation of sampled data with various attributes and methods for data manipulation, analysis, and exporting. This class handles statistical computations and provides mechanisms for thread-safe operations and data interaction. Attributes: xlabel (str): Label for the x-axis. ylabel (str): Label for the y-axis. lock (threading.RLock): A reentrant lock to ensure thread-safe operations. common_descriptions (dict|None): A dictionary containing common descriptions of the data. time_period (list): A list to store the periods of time relevant to the accumulation. job_id (int|None): An identifier for the job associated with this accumulation. use_error (bool): Indicates whether statistical error calculations are enabled. sampling_rate (float): Sampling frequency of the data. n (int): Number of accumulations. cont_data (bool): Indicates whether the object contains valid data. index (list): List of index bounds for data segments. x (list): The x-axis data points. y (list): The y-axis data points. """
[docs] def __init__(self, x: _XArray | None = None, y: _YChannels | None = None, y_2: _YChannels | numpy.ndarray | None = None, n: int | None = None, index: _IndexList | None = None, sampl_freq: float | None = None, error: bool = False) -> None: 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, overridden) self.xlabel = "Time / s" self.ylabel = "Avg. Samples [Digits]" self.lock=threading.RLock() self.is_clipped = False self.common_descriptions=None self.time_period: list[datetime.datetime] | None = [] self.job_id = None # this does not make sense for an object accumulated from multiple job_ids self.job_ids = {} self.use_error = error if self.uses_statistics(): if (y_2 is not None): # y_square is always kept as a list of per-channel arrays self.y_square: _YChannels = [numpy.asarray(ch) for ch in 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 = numpy.array([], dtype="float32") 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 = numpy.asarray(x) self.y = [numpy.asarray(ch) for ch in 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 _copy_provenance(self, r: "Accumulation") -> None: r.is_clipped = self.is_clipped r.job_id = self.job_id r.job_ids = self.job_ids.copy() r.time_period = self.time_period[:] if self.time_period is not None else None if self.common_descriptions is not None: r.common_descriptions = self.common_descriptions.copy() def get_accu_by_index(self, index: int) -> "Accumulation": 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].copy() tmp_y = [] for i in range(self.get_number_of_channels()): tmp_y.append(self.y[i][start:end+1].copy()) 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) if self.uses_statistics(): r.y_square = [] for i in range(self.get_number_of_channels()): r.y_square.append(self.y_square[i][start:end+1].copy()) if self.common_descriptions is not None: r.common_descriptions = self.common_descriptions.copy() r.time_period = self.time_period[:] if self.time_period is not None else None r.job_ids = self.job_ids.copy() r.is_clipped = self.is_clipped self.lock.release() return r def get_result_by_index(self, index: int) -> "Accumulation": return self.get_accu_by_index(index) def get_Result_by_index(self, index: int) -> "Accumulation": return self.get_accu_by_index(index) def get_ysquare(self, channel: int) -> numpy.ndarray | None: if self.uses_statistics(): try: return self.y_square[channel] except: raise else: return None def contains_data(self) -> bool: return self.cont_data def get_sampling_rate(self) -> float: """Returns the sampling frequency""" return self.sampling_rate def get_index_bounds(self, index) -> tuple[int,int]: """Returns a tuple with (start, end) of the wanted result""" return self.index[index] def uses_statistics(self) -> bool: return self.use_error # external interface -------------------------------------------------------------------- def get_yerr(self, channel: int) -> numpy.ndarray: """ return error (std.dev/sqrt(n)) of mean """ if not self.contains_data(): return numpy.array([], dtype="float32") if not self.uses_statistics(): return numpy.zeros((len(self.y[0]),),dtype="float32") 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: int = 0) -> numpy.ndarray: """ return mean data """ if not self.contains_data(): return numpy.array([], dtype="float32") 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) -> float|int: 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) -> float|int: """ Returns the global maximum value of all channels. """ 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) -> int | None: """ Returns the job identifier associated with this accumulation. """ return self.job_id # modified by Oleg Petrov def write_to_csv(self, destination: str | IO[str] = sys.stdout, delimiter: str = " ") -> None: """ Write the accumulation data to a whitespace-delimited text file with comment headers. The file starts with comment lines containing accumulation metadata and a header line with channel names. Data rows contain time and channel means, optionally with errors. Example: accu = Accumulation(...) accu.write_to_csv('data.txt', delimiter=' ') :param destination: file handle or filename to write to. Defaults to sys.stdout. :param delimiter: delimiter to use between columns. Defaults to " ". """ 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("# accumulation %d\n"%self.n) self.lock.acquire() try: if self.common_descriptions is not None: for (key,value) in self.common_descriptions.items(): dest.write("# %s : %s\n"%(key, str(value))) dest.write("# t") ch_no=self.get_number_of_channels() if self.use_error: for i in range(ch_no): dest.write(" ch%d_mean ch%d_err"%(i,i)) else: for i in range(ch_no): dest.write(" ch%d_mean"%i) dest.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)): dest.write("%e"%xdata[i]) for j in range(ch_no): if self.use_error and yerr is not None: dest.write("%s%e%s%e"%(delimiter, ydata[j][i], delimiter, yerr[j][i])) else: dest.write("%s%e"%(delimiter,ydata[j][i])) dest.write("\n") finally: self.lock.release() # ------------- added by Oleg Petrov, 14 Feb 2012 ---------------------- def write_to_simpson(self, destination: str | IO[str] = sys.stdout, delimiter: str = " ", frequency: float = 100e6) -> None: """ Write the accumulation data to a text file in Simpson format for NMRnotebook. The output contains a simple header with NP, SW, REF and TYPE=FID followed by the data block. Example: accu = Accumulation(...) accu.write_to_simpson('data.simp', frequency=100e6) :param destination: file handle or filename. Defaults to sys.stdout. :param delimiter: column delimiter for data rows. Defaults to space. :param frequency: reference frequency. Defaults to 100e6. """ if isinstance(destination, str): with open(destination, "w") as f: self._write_to_simpson_dest(f, delimiter, frequency) else: self._write_to_simpson_dest(destination, delimiter, frequency) def _write_to_simpson_dest(self, dest: IO[str], delimiter: str, frequency: float) -> None: self.lock.acquire() try: xdata=self.get_xdata() dest.write("SIMP\n") dest.write("%s%i%s"%("NP=", len(xdata), "\n")) dest.write("%s%i%s"%("SW=", self.get_sampling_rate(), "\n")) dest.write("%s%i%s"%("REF=", frequency, "\n")) dest.write("TYPE=FID\n") dest.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): dest.write("%g%s"%(ydata[j][i], delimiter)) dest.write("\n") dest.write("END\n") finally: self.lock.release() # ------------- added by Oleg Petrov, 10 Sep 2013 ----------------------- def write_to_tecmag(self, destination: str | IO[bytes] = sys.stdout.buffer, nrecords: int = 1,\ frequency: float = 100.,\ last_delay: float = 1.,\ receiver_phase: float = 0.,\ nucleus: str = '1H') -> None: """ Write the accumulation data to a binary TecMag file. The file is written in TNT1.005 format compatible with TecMag processing. Requires job_id and time_period to be set. Example: accu = Accumulation(...) accu.write_to_tecmag('data.tnt', nrecords=1, frequency=100.) :param destination: file handle or filename. :param nrecords: number of indirect dimension records for 2D experiments. :param frequency: spectrometer frequency. :param last_delay: last delay value. :param receiver_phase: receiver phase. :param nucleus: nucleus label, defaults to '1H'. """ #TODO: Function is most likely broken in Python 3 because Strings are now unicode and binary files cannot write unicode if self.job_id is 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: IO[bytes] | None = None if isinstance(destination, 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)] if self.time_period is None or len(self.time_period)==0: raise ValueError("write_to_tecmag: no time information available") 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: if isinstance(destination, str): the_destination = open(destination, "wb") else: the_destination = destination # 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.raw) else: if isinstance(destination, str): the_destination = open(destination, "rb+") else: the_destination = destination 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.raw) ydata=None finally: if isinstance(destination, str) and the_destination is not None: the_destination.close() self.lock.release() # ----------------------------------------------------------------------- def write_to_hdf(self, hdffile: tables.File, where: str, name: str, title: str, complib: str | None = None, complevel: int | None = None) -> None: """ Write the accumulation data to an HDF5 file using PyTables. Creates a group with indices, job_ids and accu_data arrays. Metadata and sampling rate are stored as attributes. Example: accu = Accumulation(...) accu.write_to_hdf(h5file, '/data', 'accu1', 'Accumulation') :param hdffile: PyTables File object to write to. :param where: location within the HDF5 file. :param name: name of the dataset group. :param title: title/description of the dataset. :param complib: compression library, e.g. 'zlib'. :param complevel: compression level. """ accu_group=hdffile.create_group(where=where,name=name,title=title) attrs = cast(Any, accu_group._v_attrs) 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: attrs.earliest_time = self.time_period[0].isoformat(sep=' ', timespec='milliseconds') attrs.oldest_time = self.time_period[1].isoformat(sep=' ', timespec='milliseconds') if self.common_descriptions is not None: for (key,value) in self.common_descriptions.items(): attrs.__setattr__("description_"+key,str(value)) attrs.__setattr__("sampling_rate",self.sampling_rate) filter=None if complib is not None: if complevel is None: complevel=3 filter=tables.Filters(complevel=complevel, complib=cast(_CompressionLibrary, complib), shuffle=True) # save interval information # start save index_table # 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=tables.Filters(complib="zlib"), expectedrows=len(self.index)) index_table.flavor="numpy" # save interval data new_row = cast(Any, 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() # end save index_table # start save job_ids jobid_table=hdffile.create_table(where=accu_group, name="job_ids", description={"job_id": tables.UInt64Col(), "job_date": tables.StringCol(128), }, title="job_ids used for accumulation", filters=filter, expectedrows=len(self.job_ids)) new_row = cast(Any, jobid_table.row) for i,job_id in enumerate(self.job_ids): new_row["job_date"] = self.job_ids[job_id].isoformat(timespec='milliseconds') new_row["job_id"]=job_id new_row.append() jobid_table.flush() # end save job_ids # 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 = (min(timedata.shape[0], 1024 * 8), timedata.shape[1]) 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() # External interfaces ------------------------------------------------------------------ # Overloaded operators ----------------------------------------------------------------- def __repr__(self) -> str: """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) + "\n" tmp_string += "jobids: " + str(self.job_ids) + "\n" return tmp_string def __add__(self, other: "int | float | ADC_Result | Accumulation | None") -> "Accumulation": """ Redefining self + other We define this for different types of other: 1) Accumulation + number 2) Accumulation + ADC_Result 3) Accumulation + Accumulation This returns a new Accumulation object. """ # Float or int: add to each channel 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) self._copy_provenance(r) self.lock.release() return r # ADC_Result elif isinstance(other, ADC_Result): # Other empty (return copy of self) if not other.contains_data(): return self + 0 # Self empty (copy ADC_Result) 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) if hasattr(other, "is_clipped"): r.is_clipped = other.is_clipped r.time_period=[other.job_date,other.job_date] if other.job_date is not None else None r.job_id = other.job_id r.common_descriptions=other.description.copy() r.job_ids[other.job_id] = other.get_job_date() self.lock.release() return r # Other and self not empty (self + other ADC_Result) else: if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add ADC_Results with different sampling-rates") if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add ADC_Results with different number of samples") if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add ADC_Results with different 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 different indexing") self.lock.acquire() 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.is_clipped = self.is_clipped or (hasattr(other, "is_clipped") and other.is_clipped) if self.time_period and other.job_date is not None: r.time_period=[min(self.time_period[0],other.job_date), max(self.time_period[1],other.job_date)] else: r.time_period=None r.job_id = other.job_id r.job_ids = self.job_ids.copy() r.job_ids[other.job_id] = other.get_job_date() 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]=self.common_descriptions[key] self.lock.release() return r # Accumulation elif isinstance(other, Accumulation): # Other empty (return copy of self) if not other.contains_data(): return self + 0 # 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() tmp_y = [] tmp_ysquare = [] for i in range(other.get_number_of_channels()): tmp_y.append(numpy.array(other.y[i], dtype="float32")) if other.uses_statistics(): tmp_ysquare.append(numpy.array(other.y_square[i], dtype="float32")) 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) other._copy_provenance(r) self.lock.release() return r # Other and self not empty (self + other) else: if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add accumulations with different sampling-rates") if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add accumulations with different number of samples") if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add accumulations with different 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 different indexing") 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() 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] + 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.is_clipped = self.is_clipped or (hasattr(other, "is_clipped") and other.is_clipped) if self.time_period and other.time_period: r.time_period=[min(self.time_period[0],other.time_period[0]), max(self.time_period[1],other.time_period[1])] else: r.time_period=None r.job_id = other.job_id r.job_ids = self.job_ids.copy() for i in other.job_ids: r.job_ids[i] = other.job_ids[i] 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]=self.common_descriptions[key] self.lock.release() return r elif other is None: # Convenience: ignore add of None return self + 0 else: raise ValueError("can not add "+repr(type(other))+" to Accumulation") def __radd__(self, other: "int | float | ADC_Result | Accumulation | None") -> "Accumulation": """Redefining other + self""" return self.__add__(other) def __sub__(self, other: "int | float | ADC_Result | Accumulation") -> "Accumulation": """Redefining self - other""" return self.__add__(-other) def __rsub__(self, other: "int | float | ADC_Result | Accumulation") -> "Accumulation": """Redefining other - self""" return (-self) + other def __mul__(self, other: int | float) -> "Accumulation": """Redefining self * other (scalar)""" if isinstance(other, (int, float)): if not self.contains_data(): raise ValueError("Accumulation: You cant multiply an empty accumulation") self.lock.acquire() tmp_y = [] tmp_ysquare = [] for i in range(self.get_number_of_channels()): tmp_y.append(self.y[i] * other) if self.uses_statistics(): tmp_ysquare.append(self.y_square[i] * (other**2)) 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 r.job_ids = self.job_ids.copy() r.time_period = self.time_period[:] if self.time_period is not None else None if self.common_descriptions is not None: r.common_descriptions = self.common_descriptions.copy() self.lock.release() return r else: raise ValueError(f"ValueError: Cannot multiply \"{other.__class__}\" to Accumulation!") def __rmul__(self, other: int | float) -> "Accumulation": """Redefining other (scalar) * self""" return self.__mul__(other) def __imul__(self, other: int | float) -> "Accumulation": """Redefining self *= other (scalar)""" if isinstance(other, (int, float)): if not self.contains_data(): raise ValueError("Accumulation: You cant multiply an empty accumulation") self.lock.acquire() for i in range(self.get_number_of_channels()): self.y[i] *= other if self.uses_statistics(): self.y_square[i] *= (other**2) self.lock.release() return self else: raise ValueError(f"ValueError: Cannot multiply \"{other.__class__}\" to Accumulation!") def __truediv__(self, other: int | float) -> "Accumulation": """Redefining self / other (scalar)""" if isinstance(other, (int, float)): if not self.contains_data(): raise ValueError("Accumulation: You cant divide an empty accumulation") self.lock.acquire() tmp_y = [] tmp_ysquare = [] for i in range(self.get_number_of_channels()): tmp_y.append(self.y[i] / other) if self.uses_statistics(): tmp_ysquare.append(self.y_square[i] / (other**2)) 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 r.job_ids = self.job_ids.copy() r.time_period = self.time_period[:] if self.time_period is not None else None if self.common_descriptions is not None: r.common_descriptions = self.common_descriptions.copy() self.lock.release() return r else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to Accumulation!") def __itruediv__(self, other: int | float) -> "Accumulation": """Redefining self /= other (scalar)""" if isinstance(other, (int, float)): if not self.contains_data(): raise ValueError("Accumulation: You cant divide an empty accumulation") self.lock.acquire() for i in range(self.get_number_of_channels()): self.y[i] /= other if self.uses_statistics(): self.y_square[i] /= (other**2) self.lock.release() return self else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to Accumulation!") def __floordiv__(self, other: int) -> "Accumulation": """Redefining self // other (scalar)""" if isinstance(other, float): raise ValueError("ValueError: Cannot use floor division (//) on floats! Use \"/\" instead of \"//\"! ") if isinstance(other, int): if not self.contains_data(): raise ValueError("Accumulation: You cant divide an empty accumulation") self.lock.acquire() tmp_y: list[numpy.ndarray] = [] for i in range(self.get_number_of_channels()): tmp_y.append(self.y[i] // other) # Note: Statistics (y_square) cannot be easily maintained with floor division 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) r.job_id = self.job_id r.job_ids = self.job_ids.copy() r.time_period = self.time_period[:] if self.time_period is not None else None if self.common_descriptions is not None: r.common_descriptions = self.common_descriptions.copy() self.lock.release() return r else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to Accumulation!") def __ifloordiv__(self, other: int) -> "Accumulation": """Redefining self //= other (scalar)""" if isinstance(other, float): raise ValueError("ValueError: Cannot use floor division (//) on floats! Use \"/\" instead of \"//\"! ") if isinstance(other, int): if not self.contains_data(): raise ValueError("Accumulation: You cant divide an empty accumulation") self.lock.acquire() for i in range(self.get_number_of_channels()): self.y[i] //= other if self.uses_statistics(): # Invalidate statistics self.y_square = [] self.use_error = False self.lock.release() return self else: raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to Accumulation!") def __iadd__(self, other: "int | float | ADC_Result | Accumulation | None") -> "Accumulation": """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 isinstance(other, 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() if hasattr(other, "is_clipped"): self.is_clipped = other.is_clipped self.time_period=[other.job_date,other.job_date] if other.job_date is not None else None self.job_id = other.job_id # added by Oleg Petrov self.job_ids[other.job_id] = other.get_job_date() self.common_descriptions=other.description.copy() return self # Other and self not empty (self + other/ADC_Result) else: if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You can't add ADC-Results with different sampling-rates") if len(self.y[0]) != len(other): raise ValueError("Accumulation: You can't add ADC-Results with different number of samples") if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You can't add ADC-Results with different number of channels") for i in range(len(self.index)): if self.index[i] != other.get_index_bounds(i): raise ValueError("Accumulation: You can't add ADC-Results with different indexing") self.lock.acquire() 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 if hasattr(other, "is_clipped"): self.is_clipped = self.is_clipped or other.is_clipped if self.time_period and other.job_date is not None: self.time_period=[min(self.time_period[0],other.job_date), max(self.time_period[1],other.job_date)] else: self.time_period=None self.job_id = other.job_id self.job_ids[other.job_id] = other.get_job_date() 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 isinstance(other, Accumulation): # Other empty (return) if not other.contains_data(): return self # 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) if hasattr(other, "is_clipped"): self.is_clipped = other.is_clipped self.common_descriptions = (other.common_descriptions.copy() if other.common_descriptions is not None else None) self.time_period=other.time_period[:] if other.time_period is not None else None self.job_id = other.job_id # added by Oleg Petrov self.job_ids.update(other.job_ids) self.lock.release() return self # Other and self not empty (self + other) else: if self.sampling_rate != other.get_sampling_rate(): raise ValueError("Accumulation: You cant add accumulations with different sampling-rates") if len(self.y[0]) != len(other): raise ValueError("Accumulation: You cant add accumulations with different number of samples") if len(self.y) != other.get_number_of_channels(): raise ValueError("Accumulation: You cant add accumulations with different 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 different indexing") 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() 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 if hasattr(other, "is_clipped"): self.is_clipped = self.is_clipped or other.is_clipped if self.time_period and other.time_period: self.time_period=[min(self.time_period[0],other.time_period[0]), max(self.time_period[1],other.time_period[1])] else: self.time_period=None self.job_id = other.job_id self.job_ids.update(other.job_ids) # Removes mismatched common description keys if self.common_descriptions is not None and other.common_descriptions is not None: # Get all keys that exist in both dictionaries common_keys = set(self.common_descriptions.keys()) & set(other.common_descriptions.keys()) # Find keys where values also match matching_keys = { key for key in common_keys if self.common_descriptions[key] == other.common_descriptions[key] } # Remove keys that don't match for key in set(self.common_descriptions.keys()) - matching_keys: 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: "int | float | ADC_Result | Accumulation") -> "Accumulation": """Redefining self -= other""" return self.__iadd__(-other) def __neg__(self) -> "Accumulation": """Redefining -self""" if not self.contains_data(): raise ValueError("Accumulation: You cant negate an empty accumulation") 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(): tmp_y_square = [numpy.array(sq) for sq in self.y_square] r = Accumulation(x = numpy.array(self.x, dtype="float32"), y = tmp_y, y_2 = tmp_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._copy_provenance(r) self.lock.release() return r
def read_from_hdf(hdf_node: object) -> "Accumulation | None": """ read accumulation data from HDF node and return it. """ # formal checks first if not isinstance(hdf_node, tables.Group): return None attrs = cast(Any, hdf_node._v_attrs) if 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 attrs._v_attrnamesuser: if attrname.startswith("description_"): accu.common_descriptions[attrname[12:]]=attrs.__getattr__(attrname) earliest_time=None if "earliest_time" in dir(attrs): timestring=attrs.__getattr__("earliest_time") # Handle both formats: "20260705 19:13:59.162" and "2026-07-05 19:13:59.162" if '-' in timestring: # New format: "2026-07-05 19:13:59.162" try: earliest_time=datetime.datetime.strptime(str(timestring), '%Y-%m-%d %H:%M:%S.%f') except ValueError: # Fallback to manual parsing parts = timestring.split() if len(parts) >= 2: date_parts = parts[0].split('-') time_parts = parts[1].split(':') if len(date_parts) == 3 and len(time_parts) >= 2: year, month, day = int(date_parts[0]), int(date_parts[1]), int(date_parts[2]) hour, minute = int(time_parts[0]), int(time_parts[1]) second_parts = time_parts[2].split('.') if len(time_parts) > 2 else ['0'] second = int(second_parts[0]) microsecond = int(second_parts[1].ljust(6, '0')[:6]) if len(second_parts) > 1 else 0 earliest_time=datetime.datetime(year, month, day, hour, minute, second, microsecond) else: # Old format: "20260705 19:13:59.162" 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(attrs): timestring=attrs.__getattr__("oldest_time") # Handle both formats: "20260705 19:13:59.162" and "2026-07-05 19:13:59.162" if '-' in timestring: # New format: "2026-07-05 19:13:59.162" try: oldest_time=datetime.datetime.strptime(str(timestring), '%Y-%m-%d %H:%M:%S.%f') except ValueError: # Fallback to manual parsing parts = timestring.split() if len(parts) >= 2: date_parts = parts[0].split('-') time_parts = parts[1].split(':') if len(date_parts) == 3 and len(time_parts) >= 2: year, month, day = int(date_parts[0]), int(date_parts[1]), int(date_parts[2]) hour, minute = int(time_parts[0]), int(time_parts[1]) second_parts = time_parts[2].split('.') if len(time_parts) > 2 else ['0'] second = int(second_parts[0]) microsecond = int(second_parts[1].ljust(6, '0')[:6]) if len(second_parts) > 1 else 0 oldest_time=datetime.datetime(year, month, day, hour, minute, second, microsecond) else: # Old format: "20260705 19:13:59.162" 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"] accu.set_title("Accumulation: n = %d" % accu.n) # 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