add FitResult fit overlays and close leaked HDF5 files after experiment stop

FitResult is a new Drawable in damaris.data with x-values independent of
the measured points, so fits can be sampled on a dense grid. It is stored
as its own data pool entry (e.g. "measurement/fit"), written to and read
from the HDF5 data pool (damaris_type "FitResult", params as attributes),
and drawn by the GUI as overlay on its measurement, including a fit
parameter text box. Fits get no entry in the display tree, so measurement
rows stay directly selectable. Result scripts create them via
FitResult(x, y, title, source, params).

Additionally, HDF5 files that experiment or result scripts left open
(e.g. after a script failure) are closed when the experiment stops.
pytables keeps such handles registered process-wide, which blocked
reopening the dump file and the scripts' own files on the next run.

Generated by Mistral Vibe.
Co-Authored-By: Mistral Vibe <vibe@mistral.ai>

doc/api_reference.rst

@@ -32,6 +32,7 @@ The data module contains result types, data handling, and signal processing clas
    damaris.data.Drawable
    damaris.data.Errorable
    damaris.data.Error_Result
+   damaris.data.FitResult
    damaris.data.MeasurementResult
    damaris.data.Persistance
    damaris.data.Resultable

doc/manual.rst

@@ -285,6 +285,35 @@ A typical method to already extract data like amplitudes etc. is possible by usi
 
 Essentially, a dictionary "measurement" is created with *tau* as the key. For each new key the current magnetization is added. This dictionary is then added to the data pool. The advantage over the previous listing is that the magnetization vs. τ curve can be directly displayed in order to gain an overview of the running experiment, because the "measurement" dictionary is added to the data pool. In this example, the value also stores statistics from the accumulation. The data dictionary is stored automatically as described before, including the result and experiment scripts, information about the used back end and a timeline of the experiment.
 
+**Displaying fits**
+
+A fit to a measurement can be shown as an overlay on the measurement curve using the FitResult class. In contrast to the measurement itself, a fit has its own x-values, independent of the measured data points. Therefore the fit can be sampled on a much denser x-grid, which makes the plotted curve smooth. The fit is stored in the data pool under its own key and is drawn by the GUI on top of the measurement it belongs to. Fit parameters are shown in the plot and are stored together with the fit.
+
+The fit is associated with its measurement in two ways: by the ``source`` argument (the data pool key of the measurement) and, optionally, by storing it under a nested key such as ``"Magnetization(tau)/fit"``, so that it is stored next to the measurement in the HDF5 data pool file. The GUI shows the fit only as an overlay on its measurement; fits do not get their own entry in the display tree.
+
+.. code-block:: python
+   :linenos:
+
+   import numpy
+   from damaris import FitResult
+
+   def result(results):
+       # ... fill "measurement" as in the listing above ...
+       # fit an exponential decay to the measured points
+       from scipy.optimize import curve_fit
+       [k, v, e] = measurement.get_errorplotdata()
+       popt, pcov = curve_fit(lambda t, m0, t1: m0 * numpy.exp(-t / t1),
+                              k, v, p0=(v[0], 1.0), sigma=e)
+       # sample the fit on a dense x-grid for a smooth curve
+       xfit = numpy.linspace(k.min(), k.max(), 200)
+       fit = FitResult(xfit, popt[0] * numpy.exp(-xfit / popt[1]),
+                       title="Magnetization(tau) fit",
+                       source="Magnetization(tau)",
+                       params={"M0": popt[0], "T1": popt[1]})
+       data["Magnetization(tau)/fit"] = fit
+
+The fit is written to the HDF5 data pool file as its own table (columns *x* and *y*, plus *y_err* if given), together with the fit parameters as attributes, and is restored when the data pool is read back.
+
 Note that we can also create a HDF5 file from scratch! In this case experiment and result script are not saved in the resulting file. The advantage is that one can create arbitrary complex structured files at the cost of having a more complicated result script.
 
 .. code-block:: python

src/damaris/data/DataPool.py

@@ -394,6 +394,12 @@ class DataPool(collections.abc.MutableMapping):
                         if obj is not None:
                             self[full_path] = obj
 
+                    elif damaris_type == "FitResult":
+                        from damaris.data.FitResult import read_from_hdf
+                        obj = read_from_hdf(node)
+                        if obj is not None:
+                            self[full_path] = obj
+
                     elif damaris_type == "TemperatureResult":
                         from damaris.data.Temperature import read_from_hdf
                         obj = read_from_hdf(node)

src/damaris/data/FitResult.py

@@ -0,0 +1,127 @@
+# fit curve with x-values independent of the measurement data points
+
+import numpy
+import tables
+from . import Drawable
+
+
+class FitResult(Drawable.Drawable):
+    """
+    plottable fit curve belonging to a MeasurementResult
+
+    In contrast to MeasurementResult, the x-values are a plain array
+    and are independent of the measured data points, so a fit can be
+    sampled on a much denser x-grid for a smooth plot. The fit is
+    stored in the data pool under its own key, e.g. as a child of the
+    measurement ("Magnetization(tau)/fit"), and is drawn by the GUI
+    as an overlay on the measurement it belongs to.
+    """
+
+    def __init__(self, x, y, title=None, source=None, y_err=None, params=None):
+        """
+        x: array-like of x-values (independent of the measurement)
+        y: array-like of y-values, same length as x
+        title: plot title of the fit
+        source: data pool key of the MeasurementResult this fit belongs to
+        y_err: optional array-like of y errors, same length as x
+        params: optional dictionary of fit parameters, shown in the plot
+        """
+        Drawable.Drawable.__init__(self)
+        self.x = numpy.array(x, dtype="float64").reshape(-1)
+        self.y = [numpy.array(y, dtype="float64").reshape(-1)]
+        if y_err is not None:
+            self.y_err = numpy.array(y_err, dtype="float64").reshape(-1)
+        else:
+            self.y_err = None
+        self.title = title
+        self.source = source
+        self.params = dict(params) if params is not None else {}
+
+    def get_lineplotdata(self):
+        """
+        x and y data for plotting the fit curve
+        """
+        return [self.x, self.y[0]]
+
+    def uses_statistics(self):
+        """
+        drawable interface method returns False
+        """
+        return False
+
+    def write_to_hdf(self, hdffile, where, name, title, complib=None, complevel=None):
+        h5_table_format = {
+            "x": tables.Float64Col(),
+            "y": tables.Float64Col()
+            }
+        if self.y_err is not None:
+            h5_table_format["y_err"] = tables.Float64Col()
+        filter = None
+        if complib is not None:
+            if complevel is None:
+                complevel = 9
+            filter = tables.Filters(complevel=complevel, complib=complib, shuffle=1)
+
+        fit_table = hdffile.create_table(where=where, name=name,
+                                         description=h5_table_format,
+                                         title=title,
+                                         filters=filter,
+                                         expectedrows=self.x.shape[0])
+        fit_table.flavor = "numpy"
+        fit_table.attrs.damaris_type = "FitResult"
+        fit_table.attrs.quantity_name = self.get_title()
+        if self.xlabel is not None:
+            fit_table.attrs.xlabel = self.xlabel
+        if self.ylabel is not None:
+            fit_table.attrs.ylabel = self.ylabel
+        if self.source is not None:
+            fit_table.attrs.source = self.source
+        for (key, value) in self.params.items():
+            if isinstance(value, (int, float)) and not isinstance(value, bool):
+                fit_table.attrs.__setattr__("param_" + str(key), float(value))
+            else:
+                fit_table.attrs.__setattr__("param_" + str(key), str(value))
+
+        row = fit_table.row
+        has_y_err = self.y_err is not None
+        for i in range(self.x.shape[0]):
+            row["x"] = self.x[i]
+            row["y"] = self.y[0][i]
+            if has_y_err:
+                row["y_err"] = self.y_err[i]
+            row.append()
+        fit_table.flush()
+
+
+def read_from_hdf(hdf_node):
+    """
+    reads a FitResult 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 != "FitResult":
+        return None
+
+    xdata = hdf_node.col("x")
+    ydata = hdf_node.col("y")
+    if "y_err" in hdf_node.colnames:
+        y_err = hdf_node.col("y_err")
+    else:
+        y_err = None
+
+    fr = FitResult(xdata, ydata,
+                   title=getattr(hdf_node._v_attrs, "quantity_name", None),
+                   source=getattr(hdf_node._v_attrs, "source", None),
+                   y_err=y_err)
+    fr.xlabel = getattr(hdf_node._v_attrs, "xlabel", None)
+    fr.ylabel = getattr(hdf_node._v_attrs, "ylabel", None)
+
+    params = {}
+    for attr_name in hdf_node._v_attrs._f_list("all"):
+        if attr_name.startswith("param_"):
+            params[attr_name[6:]] = getattr(hdf_node._v_attrs, attr_name)
+    fr.params = params
+
+    return fr

src/damaris/data/__init__.py

@@ -1,10 +1,11 @@
 from damaris.data.ADC_Result import ADC_Result
 from damaris.data.Accumulation import Accumulation
 from damaris.data.MeasurementResult import MeasurementResult, AccumulatedValue
+from damaris.data.FitResult import FitResult
 from damaris.data.DataPool import DataPool
 from damaris.data.Error_Result import Error_Result
 from damaris.data.Config_Result import Config_Result
 from damaris.data.Temperature import TemperatureResult
 
-__all__=["ADC_Result", "Accumulation", "MeasurementResult", "AccumulatedValue", "DataPool", "Error_Result", "Config_Result", "TemperatureResult" ]
+__all__=["ADC_Result", "Accumulation", "MeasurementResult", "AccumulatedValue", "FitResult", "DataPool", "Error_Result", "Config_Result", "TemperatureResult" ]
 

src/damaris/gui/DamarisGUI.py

@@ -90,7 +90,7 @@ from damaris.gui import ResultReader, ResultHandling
 from damaris.gui import BackendDriver
 #from damaris.gui.gtkcodebuffer import CodeBuffer, SyntaxLoader
 #from damaris.data import Drawable # this is a base class, it should be used...
-from damaris.data import DataPool, Accumulation, ADC_Result, MeasurementResult, TemperatureResult
+from damaris.data import DataPool, Accumulation, ADC_Result, MeasurementResult, TemperatureResult, FitResult
 # TemperatureMonitor imported lazily to avoid circular import issues
 
 # default, can be set to true from start script
@@ -895,6 +895,11 @@ class DamarisGUI:
                     dump_size = os.stat( self.dump_filename ).st_size / 1e6
                     logger.info("Done (%.1f s, %.1f MB)", time.time( ) - self.dump_start_time, dump_size)
 
+            # close hdf5 files that the scripts left open (e.g. after a
+            # script failure), otherwise pytables keeps them registered
+            # as open and they block the next experiment run
+            self.close_left_open_hdf5_files( )
+
             # keep data to display but throw away everything else
             self.si = None
             # delete lock file so that other experiment can start
@@ -1113,6 +1118,26 @@ class DamarisGUI:
 
         return True
 
+    def close_left_open_hdf5_files( self ):
+        """
+        close hdf5 files that experiment or result scripts left open,
+        e.g. after a script failure; pytables keeps such files registered
+        as open, so they cannot be reopened in the next experiment run
+        """
+        open_files = getattr( tables.file, "_open_files", None )
+        if open_files is None:
+            return
+        if hasattr( open_files, "handlers" ):
+            filehandles = list( open_files.handlers )
+        else:
+            filehandles = list( open_files.values( ) )
+        for filehandle in filehandles:
+            try:
+                logger.warning( "Closing hdf5 file left open by script: %s", filehandle.filename )
+                filehandle.close( )
+            except Exception:
+                pass
+
     def pause_experiment( self, widget, data=None ):
         """
         pause experiment execution (that means delay backend and let others run)
@@ -3101,6 +3126,8 @@ class MonitorWidgets:
         self.update_counter = 0
         self.update_counter_lock = threading.Lock( )
         self.description_text = None
+        self.fitresultlines = [ ]
+        self.fit_source_keys = set( )
         self.clipping_marker = None
         self.source_list_cache = {}
 
@@ -3260,6 +3287,7 @@ class MonitorWidgets:
         self.source_list_reset( )
         self.__rescale = True
         self.displayed_data = [ None, None ]
+        self.fit_source_keys = set( )
         self.display_source_path_label.set_label( "" )
         self.clear_display( )
 
@@ -3389,9 +3417,18 @@ class MonitorWidgets:
                         self.renew_display( )
                 new_data_struct = None
         elif event.what == DataPool.Event.new_key:
-            # update combo-box by inserting and rely on consistent information
-            self.source_list_add( event.subject )
+            if isinstance( self.data_pool.get( event.subject ), FitResult ):
+                # fits are shown as overlay on their measurement only,
+                # they do not get their own entry in the display tree
+                self.fit_source_keys.add( event.subject )
+            else:
+                # update combo-box by inserting and rely on consistent information
+                self.source_list_add( event.subject )
         elif event.what == DataPool.Event.deleted_key:
+            if event.subject in self.fit_source_keys:
+                # fit entries are not in the display tree
+                self.fit_source_keys.discard( event.subject )
+                return
             # update combo-box by removing and rely on consistent information
             if (self.displayed_data[ 0 ] is not None and
                         self.displayed_data[ 0 ] == event.subject):
@@ -3724,6 +3761,15 @@ class MonitorWidgets:
             self.measurementresultgraph = None
             self.matplot_axes.clear( )
             self.matplot_axes.grid( True )
+        if not hasattr( self, "fitresultlines" ):
+            self.fitresultlines = [ ]
+        elif self.fitresultlines:
+            for line in self.fitresultlines:
+                try:
+                    line.remove( )
+                except:
+                    pass
+            self.fitresultlines = [ ]
         if not hasattr( self, "graphen" ):
             self.graphen = [ ]
         elif self.graphen:
@@ -3966,10 +4012,36 @@ class MonitorWidgets:
                     self.matplot_axes.lines.remove( self.measurementresultline[ 0 ] )
                     self.measurementresultline = None
 
+            # remove fit overlay lines
+            for line in self.fitresultlines:
+                try:
+                    line.remove( )
+                except:
+                    pass
+            self.fitresultlines = [ ]
+
+            # collect fits belonging to this measurement:
+            # either by source attribute or by nested data pool key
+            displayed_key = self.displayed_data[ 0 ]
+            fit_results = [ ]
+            for key in self.data_pool.keys( ):
+                if key == displayed_key:
+                    continue
+                fit = self.data_pool.get( key )
+                if isinstance( fit, FitResult ) and \
+                        (fit.source == displayed_key or key.startswith( displayed_key + "/" )):
+                    fit_results.append( (key, fit) )
+            fit_results.sort( key=lambda key_fit: key_fit[ 0 ] )
+            fit_xydata = [ fit.get_lineplotdata( ) for key, fit in fit_results ]
+
             [ k, v, e ] = in_result.get_errorplotdata( )
             if k.shape[ 0 ] != 0:
                 xmin = k.min( )
                 ymin = (v - e).min( )
+                for fx, fy in fit_xydata:
+                    if fx.shape[ 0 ] != 0:
+                        xmin = min( xmin, fx.min( ) )
+                        ymin = min( ymin, fy.min( ) )
 
                 if xmin > 0:
                     self.display_x_scaling_combobox.set_sensitive( True )
@@ -3997,6 +4069,10 @@ class MonitorWidgets:
                 if self.__rescale:  # or self.display_autoscaling_checkbutton.get_active():
                     xmax = k.max( )
                     ymax = (v + e).max( )
+                    for fx, fy in fit_xydata:
+                        if fx.shape[ 0 ] != 0:
+                            xmax = max( xmax, fx.max( ) )
+                            ymax = max( ymax, fy.max( ) )
                     # is there a range problem?
                     if xmin == xmax or ymin == ymax:
                         # fix range and scaling problems
@@ -4038,6 +4114,44 @@ class MonitorWidgets:
             if k.shape[ 0 ] != 0 and v.shape == k.shape:
                 self.measurementresultline = self.matplot_axes.plot( k, v, 'r-' )
 
+            # draw the fits as overlay lines
+            fit_colors = [ "green", "purple", "darkorange" ]
+            for i, (fx, fy) in enumerate( fit_xydata ):
+                if fx.shape[ 0 ] != 0 and fy.shape == fx.shape:
+                    self.fitresultlines.extend(
+                        self.matplot_axes.plot( fx, fy, "-",
+                                                color=fit_colors[ i % len( fit_colors ) ],
+                                                linewidth=1.5 ) )
+
+            # show fit parameters
+            description_lines = [ ]
+            for fit_key, fit in fit_results:
+                if fit.params:
+                    if len( fit_results ) > 1:
+                        description_lines.append( str( fit.get_title( ) or fit_key ) )
+                    for param_key in fit.params:
+                        param_value = fit.params[ param_key ]
+                        if isinstance( param_value, float ):
+                            description_lines.append( "%s = %.4g" % (param_key, param_value) )
+                        else:
+                            description_lines.append( "%s = %s" % (param_key, param_value) )
+            if description_lines:
+                description_string = "\n".join( description_lines )
+                if self.description_text is None:
+                    self.description_text = self.matplot_axes.text(0.98, 0.95, description_string,
+                                            size=10,
+                                            transform=self.matplot_axes.transAxes,
+                                            va='top', ha='right', ma='left',
+                                            bbox=dict(facecolor='lightgray', alpha=0.5, boxstyle='round,pad=0.2'))
+                else:
+                    self.description_text.set_text( description_string )
+            elif self.description_text is not None:
+                try:
+                    self.description_text.remove()
+                except:
+                    pass
+                self.description_text = None
+
 
             # Any title to be set?
             title = in_result.get_title( ) + ""

tests/res_test.py

@@ -1,6 +1,6 @@
 import tables
 import numpy
-
+# from damaris import FitResult
 r = MeasurementResult("test")
 
 def result():
@@ -16,6 +16,19 @@ def result():
         ts.write_to_hdf(h5, "/", f"adc_{num}", f"adc_{num}")
         data["fft"] = (ts+0).fft()
         r[num] = ts.y[0][0]
+        if num > 4:
+            x = r.get_xdata()
+            y = r.get_ydata()
+            coeffs = numpy.polyfit(x,  y, 2)
+            x_line = numpy.linspace(x.min(), x.max(), 100)
+            y_line = numpy.polyval(coeffs, x_line)
+            
+            fit = FitResult(x=x_line,
+                    y=y_line,
+                    title="Magnetization(tau) fit",
+                    source="test",
+                    params={"a": coeffs[0], "b": coeffs[1], "c": coeffs[2]})
+            data["test/fit"] = fit
         data["test"] = r
         accu+=ts
         print(accu)