Merge pull request #886 from BCDA-APS/884-signal-stats

new lineup2() plan can be used with queueserver, console, and notebooks

Author: prjemian

CHANGES.rst

@@ -34,6 +34,7 @@ New Features
 ------------
 
 * Add ophyd device support for DG-645 digital delay/pulse generator.
+* New lineup2() plan can be used in console, notebooks, and queueserver.
 
 1.6.17
 ******

apstools/callbacks/__init__.py

@@ -5,6 +5,8 @@
 from .nexus_writer import NEXUS_RELEASE
 from .nexus_writer import NXWriter
 from .nexus_writer import NXWriterAPS
+from .scan_signal_statistics import factor_fwhm
+from .scan_signal_statistics import SignalStatsCallback
 from .spec_file_writer import SCAN_ID_RESET_VALUE
 from .spec_file_writer import SPEC_TIME_FORMAT
 from .spec_file_writer import SpecWriterCallback

apstools/callbacks/scan_signal_statistics.py

@@ -0,0 +1,185 @@
+"""
+Collect statistics on the signals used in 1-D scans.
+====================================================
+
+.. autosummary::
+
+    ~factor_fwhm
+    ~SignalStatsCallback
+"""
+
+__all__ = """
+factor_fwhm
+SignalStatsCallback
+""".split()
+
+import math
+import logging
+
+import pyRestTable
+import pysumreg
+
+logger = logging.getLogger(__name__)
+logger.info(__file__)
+
+factor_fwhm = 2 * math.sqrt(2 * math.log(2))
+r"""
+FWHM :math:`=2\sqrt{2\ln{2}}\cdot\sigma_c`
+
+see: https://statproofbook.github.io/P/norm-fwhm.html
+"""
+
+
+class SignalStatsCallback:
+    """
+    Callback: Collect peak (& other) statistics during a scan.
+
+    .. caution:: This is an early draft and is subject to change!
+
+    Subscribe the ``receiver()`` method. Use with step scan plans such as
+    ``bp.scan()`` and ``bp.rel_scan()``.
+
+    .. caution:: It is recommended to subscribe this callback to specific plans.  It
+        should not be run with just any plan (it could easily raise exceptions).
+
+    .. rubric:: Basic example
+    .. code-block::
+        :linenos:
+
+        from bluesky import plans as bp
+        from bluesky import preprocessors as bpp
+
+        signal_stats = SignalStats()
+
+        def my_plan(detectors, mover, rel_start, rel_stop, points, md={}):
+
+            @bpp.subs_decorator(signal_stats.receiver)  # collect the data
+            def _inner():
+                yield from bp.rel_scan(detectors, mover, rel_start, rel_end, points, md)
+
+            yield from _inner()  # run the scan
+            signal_stats.report()  # print the statistics
+
+    .. rubric:: Public API
+    .. autosummary::
+
+        ~receiver
+        ~report
+        ~data_stream
+        ~stop_report
+
+    .. rubric:: Internal API
+    .. autosummary::
+
+        ~clear
+        ~descriptor
+        ~event
+        ~start
+        ~stop
+        ~_scanning
+        ~_registers
+    """
+
+    data_stream: str = "primary"
+    """RunEngine document with signals to to watch."""
+
+    stop_report: bool = True
+    """If ``True`` (default), call the ``report()`` method when a ``stop`` document is received."""
+
+    _scanning: bool = False
+    """Is a run *in progress*?"""
+
+    _registers: dict = {}
+    """Dictionary (keyed on Signal name) of ``SummationRegister()`` objects."""
+
+    # TODO: What happens when the run is paused?
+
+    def __repr__(self):
+        if "_motor" not in dir(self):  # not initialized
+            self.clear()
+        args = f"motor={self._motor!r},  detectors={self._detectors!r}"
+        return f"{self.__class__.__name__}({args})"
+
+    def clear(self):
+        """Clear the internal memory for the next run."""
+        self._scanning = False
+        self._detectors = []
+        self._motor = ""
+        self._registers = {}
+        self._descriptor_uid = None
+        self._x_name = None
+        self._y_names = []
+
+    def descriptor(self, doc):
+        """Receives 'descriptor' documents from the RunEngine."""
+        if not self._scanning:
+            return
+        if doc["name"] != self.data_stream:
+            return
+
+        # Remember now, to match with later events.
+        self._descriptor_uid = doc["uid"]
+
+        # Pick the first motor signal.
+        self._x_name = doc["hints"][self._motor]["fields"][0]
+        # Get the signals for each detector object.s
+        for d in self._detectors:
+            self._y_names += doc["hints"][d]["fields"]
+
+        # Keep statistics for each of the Y signals (vs. the one X signal).
+        self._registers = {y: pysumreg.SummationRegisters() for y in self._y_names}
+
+    def event(self, doc):
+        """Receives 'event' documents from the RunEngine."""
+        if not self._scanning:
+            return
+        if doc["descriptor"] != self._descriptor_uid:
+            return
+
+        # Collect the data for the signals.
+        x = doc["data"][self._x_name]
+        for yname in self._y_names:
+            self._registers[yname].add(x, doc["data"][yname])
+
+    def receiver(self, key, document):
+        """Client method used to subscribe to the RunEngine."""
+        handlers = "start stop descriptor event".split()
+        if key in handlers:
+            getattr(self, key)(document)
+        else:
+            logger.debug("%s: unhandled document type: %s", self.__class__.__name__, key)
+
+    def report(self):
+        """Print a table with the collected statistics for each signal."""
+        if len(self._registers) == 0:
+            return
+        keys = "n centroid sigma x_at_max_y max_y min_y mean_y stddev_y".split()
+        table = pyRestTable.Table()
+        table.labels = ["detector"] + keys
+        for yname, stats in self._registers.items():
+            row = [yname]
+            for k in keys:
+                try:
+                    v = getattr(stats, k)
+                except ZeroDivisionError:
+                    v = 0
+                row.append(v)
+            table.addRow(row)
+        print(f"Motor: {self._x_name}")
+        print(table)
+
+    def start(self, doc):
+        """Receives 'start' documents from the RunEngine."""
+        self.clear()
+        self._scanning = True
+        # These command arguments might each have many signals.
+        self._detectors = doc["detectors"]
+        self._motor = doc["motors"][0]  # just keep the first one
+
+    def stop(self, doc):
+        """Receives 'stop' documents from the RunEngine."""
+        if not self._scanning:
+            return
+        self._scanning = False
+        if self.stop_report:
+            self.report()

apstools/plans/alignment.py

@@ -5,6 +5,7 @@
 .. autosummary::
 
    ~lineup
+   ~lineup2
    ~tune_axes
    ~TuneAxis
    ~TuneResults
@@ -42,68 +43,59 @@ def lineup(
     # fmt: on
 ):
     """
-    Lineup and center a given axis, relative to current position.
+    (use ``lineup2()`` now) Lineup and center a given axis, relative to current position.
 
     If first run identifies a peak, makes a second run to fine tune the result.
 
+    ..caution:: ``lineup()`` does not work in the queueserver.  Use
+        :func:`~apstools.plans.alignment.lineup2()` instead.
+
     .. index:: Bluesky Plan; lineup
 
     PARAMETERS
 
     detectors
-        *[object]* or *object* :
-        Instance(s) of ophyd.Signal (or subclass such as ophyd.scaler.ScalerChannel)
-        dependent measurement to be maximized.  If a list, the first signal in the
-        list will be used.
+        *[object]* or *object* : Instance(s) of ophyd.Signal (or subclass such
+        as ophyd.scaler.ScalerChannel) dependent measurement to be maximized.
+        If a list, the first signal in the list will be used.
 
     axis
-        movable *object* :
-        instance of ophyd.Signal (or subclass such as EpicsMotor)
-        independent axis to use for alignment
+        movable *object* : instance of ophyd.Signal (or subclass such as
+        EpicsMotor) independent axis to use for alignment
 
     minus
-        *float* :
-        first point of scan at this offset from starting position
+        *float* : first point of scan at this offset from starting position
 
     plus
-        *float* :
-        last point of scan at this offset from starting position
+        *float* : last point of scan at this offset from starting position
 
     npts
-        *int* :
-        number of data points in the scan
+        *int* : number of data points in the scan
 
     time_s
-        *float* :
-        Count time per step (if detectors[0] is ScalerChannel object),
-        other object types not yet supported.
-        (default: 0.1)
+        *float* : Count time per step (if detectors[0] is ScalerChannel object),
+        other object types not yet supported. (default: 0.1)
 
     peak_factor
-        *float* :
-        peak maximum must be greater than ``peak_factor*minimum``
+        *float* : peak maximum must be greater than ``peak_factor*minimum``
         (default: 4)
 
     width_factor
-        *float* :
-        fwhm must be less than ``width_factor*plot_range``
-        (default: 0.8)
+        *float* : fwhm must be less than ``width_factor*plot_range`` (default:
+        0.8)
 
     feature
-        *str* :
-        One of the parameters returned by BestEffortCallback peak stats.
-        (``cen``, ``com``, ``max``, ``min``)
-        (default: ``cen``)
+        *str* : One of the parameters returned by BestEffortCallback peak stats.
+        (``cen``, ``com``, ``max``, ``min``) (default: ``cen``)
 
     rescan
-        *bool* :
-        If first scan indicates a peak, should a second scan refine the result?
-        (default: ``True``)
+        *bool* : If first scan indicates a peak, should a second scan refine the
+        result? (default: ``True``)
 
     bec
-        *object* :
-        Instance of ``bluesky.callbacks.best_effort.BestEffortCallback``.
-        (default: ``None``, meaning look for it from global namespace)
+        *object* : Instance of
+        ``bluesky.callbacks.best_effort.BestEffortCallback``. (default:
+        ``None``, meaning look for it from global namespace)
 
     EXAMPLE::
 
@@ -221,6 +213,165 @@ def peak_analysis():
         scaler.stage_sigs = old_sigs
 
 
+def lineup2(
+    # fmt: off
+    detectors, mover, rel_start, rel_end, points,
+    peak_factor=2.5, width_factor=0.8,
+    feature="centroid",
+    nscans=2,
+    signal_stats=None,
+    md={},
+    # fmt: on
+):
+    """
+    Lineup and center a given mover, relative to current position.
+
+    This plan can be used in the queueserver, Jupyter notebooks, and IPython
+    consoles.  It does not require the bluesky BestEffortCallback.  Instead, it
+    uses *PySumReg*  [#pysumreg]_ to compute statistics for each signal in a 1-D
+    scan.
+
+    New in release 1.6.18
+
+    .. caution:: This is an early draft and is subject to change!
+
+    .. index:: Bluesky Plan; lineup2; lineup
+
+    PARAMETERS
+
+    detectors *Readable* or [*Readable*]:
+        Detector object or list of detector objects (each is a Device or
+        Signal). If a list, the first Signal will be used for alignment.
+
+    mover *Movable*:
+        Mover object, such as motor or other positioner.
+
+    rel_start *float*:
+        Starting point for the scan, relative to the current mover position.
+
+    rel_end *float*:
+        Ending point for the scan, relative to the current mover position.
+
+    points *int*:
+        Number of points in the scan.
+
+    peak_factor *float* :
+        peak maximum must be greater than ``peak_factor*minimum`` (default: 2.5)
+
+    width_factor *float* :
+        fwhm must be less than ``width_factor*plot_range`` (default: 0.8)
+
+    feature *str*:
+        Use this statistical measure (default: centroid) to set the mover
+        position after a peak has been found.  Must be one of these values:
+
+        ==========  ====================
+        feature     description
+        ==========  ====================
+        centroid    center of mass
+        x_at_max_y  x location of y maximum
+        x_at_min_y  x location of y minimum
+        ==========  ====================
+
+        Statistical analysis provided by *PySumReg*.  [#pysumreg]_
+
+        .. [#pysumreg] https://prjemian.github.io/pysumreg/latest/
+
+    nscans *int*:
+        Number of scans.  (default: 2)  Scanning will stop if any scan cannot
+        find a peak.
+
+    signal_stats *object*:
+        Caller could provide an object of
+        :class:`~apstools.callbacks.scan_signal_statistics.SignalStatsCallback`.
+        If ``None``, this function will search for the ``signal_stats`` signal
+        in the global namespace. If not still found, a new one will be created
+        for the brief lifetime of this function.
+
+    md *dict*:
+        User-supplied metadata for this scan.
+    """
+    from ..callbacks import SignalStatsCallback
+    from ..callbacks import factor_fwhm
+
+    if not isinstance(detectors, (tuple, list)):
+        detectors = [detectors]
+
+    _md = dict(purpose="alignment")
+    _md.update(md or {})
+
+    if signal_stats is None:
+        signal_stats = utils.ipython_shell_namespace().get("signal_stats")
+        if signal_stats is None:
+            signal_stats = SignalStatsCallback()
+            signal_stats.stop_report = True
+    else:
+        signal_stats.stop_report = False  # Turn this automation off.
+
+    # Allow for feature to be defined using a name from PeakStats.
+    xref_PeakStats = {
+        "com": "centroid",
+        "cen": "x_at_max_y",
+        "max": "x_at_max_y",
+        "min": "x_at_min_y",
+    }
+    # translate from PeakStats feature to SignalStats
+    feature = xref_PeakStats.get(feature, feature)
+
+    def get_x_by_feature():
+        """Return the X value of the specified ``feature``."""
+        stats = principal_signal_stats()
+        if strong_peak(stats) and not too_wide(stats):
+            return getattr(stats, feature)
+
+    def principal_signal_stats() -> str:
+        """Return the name of the first detector Signal."""
+        return signal_stats._registers[signal_stats._y_names[0]]
+
+    def strong_peak(stats) -> bool:
+        """Determine if the peak is strong."""
+        try:
+            value = (stats.max_y - stats.min_y) / stats.sigma
+            return value > peak_factor
+        except ZeroDivisionError:  # not enough samples
+            try:
+                value = abs(stats.max_y / stats.min_y)
+                return value > peak_factor
+            except ZeroDivisionError:
+                return False
+
+    def too_wide(stats):
+        """Does the measured peak width fill the full range of X?"""
+        try:
+            x_range = stats.max_x - stats.min_x
+            fwhm = stats.sigma * factor_fwhm
+            return fwhm > width_factor * x_range
+        except ZeroDivisionError:  # not enough samples
+            return True
+
+    @bpp.subs_decorator(signal_stats.receiver)
+    def _inner():
+        """Run the scan, collecting statistics at each step."""
+        # TODO: save signal stats into separate stream
+        yield from bp.rel_scan(detectors, mover, rel_start, rel_end, points, md=_md)
+
+    while nscans > 0:  # allow for repeated scans
+        yield from _inner()  # Run the scan.
+        nscans -= 1
+
+        target = get_x_by_feature()
+        if target is None:
+            nscans = 0  # Nothing found, no point scanning again.
+        else:
+            yield from bps.mv(mover, target)  # Move to the feature position.
+            logger.info("Moved %s to %s: %f", mover.name, feature, target)
+
+            if nscans > 0:
+                # move the end points for the next scan
+                rel_end = principal_signal_stats().sigma * factor_fwhm
+                rel_start = -rel_end
+
+
 class TuneAxis(object):
     """
     tune an axis with a signal

apstools/plans/tests/test_alignment.py

@@ -169,6 +169,33 @@ def test_lineup(signals, mover, start, finish, npts, feature, rescan):
     #     # assert lo <= position <= hi, f"{bec=} {bec.peaks=} {position=} {center=} {width=}"
 
 
+@pytest.mark.parametrize(
+    "signals, mover, start, finish, npts, feature, nscans",
+    [
+        [noisy, m1, -1.2, 1.2, 11, "max", 2],  # slower, is motor
+        [pvoigt, axis, -1.2, 1.2, 11, "cen", 2],  # faster, is ao (float)
+        [pvoigt, axis, -1.2, 1.2, 11, "max", 2],
+        [[pvoigt], axis, -1.2, 1.2, 11, "max", 2],  # list
+        [[pvoigt, noisy, scaler1], axis, -1.2, 1.2, 11, "max", 2],  # more than one detector
+        [(pvoigt), axis, -1.2, 1.2, 11, "max", 2],  # tuple
+        [pvoigt, axis, -1.2, 1.2, 11, "cen", 1],
+        [pvoigt, axis, -1.2, 1.2, 11, "com", 1],
+        [pvoigt, axis, -1.2, 1.2, 11, "max", 1],
+        [pvoigt, axis, -1.2, 1.2, 11, "min", 1],  # pathological
+    ],
+)
+def test_lineup2(signals, mover, start, finish, npts, feature, nscans):
+    if isinstance(signals, SynPseudoVoigt):
+        signals.randomize_parameters(scale=250_000, bkg=0.000_000_000_1)
+    else:
+        change_noisy_parameters()
+
+    RE(bps.mv(mover, 0))
+    assert get_position(mover) == 0.0
+
+    RE(alignment.lineup2(signals, mover, start, finish, npts, feature=feature, nscans=nscans))
+
+
 def test_TuneAxis():
     signal = SynPseudoVoigt(name="signal", motor=m1, motor_field=m1.name)
     signal.kind = "hinted"

docs/source/api/_callbacks.rst

@@ -8,6 +8,10 @@ Callbacks
 .. automodule:: apstools.callbacks.doc_collector
     :members:
 
+.. automodule:: apstools.callbacks.scan_signal_statistics
+    :members:
+    :private-members:
+
 
 File Writers
 ------------

docs/source/api/_plans.rst

@@ -32,6 +32,7 @@ Custom Scans
    ~apstools.plans.doc_run.documentation_run
    ~apstools.plans.labels_to_streams.label_stream_decorator
    ~apstools.plans.alignment.lineup
+   ~apstools.plans.alignment.lineup2
    ~apstools.plans.nscan_support.nscan
    ~apstools.plans.sscan_support.sscan_1D
    ~apstools.plans.alignment.TuneAxis
@@ -50,6 +51,7 @@ Overall
    ~apstools.plans.command_list.execute_command_list
    ~apstools.plans.command_list.get_command_list
    ~apstools.plans.alignment.lineup
+   ~apstools.plans.alignment.lineup2
    ~apstools.plans.nscan_support.nscan
    ~apstools.plans.command_list.parse_Excel_command_file
    ~apstools.plans.command_list.parse_text_command_file