Coverage for python/lsst/summit/utils/tmaUtils.py: 14%
750 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-19 02:51 -0800
1# This file is part of summit_utils.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21from __future__ import annotations
23import datetime
24import enum
25import itertools
26import logging
27import re
28from collections.abc import Iterable
29from dataclasses import dataclass, field
30from typing import TYPE_CHECKING, Any, Callable
32import astropy
33import humanize
34import matplotlib
35import matplotlib.dates as mdates
36import matplotlib.pyplot as plt
37import numpy as np
38import pandas as pd
39from astropy.time import Time
40from matplotlib.ticker import FuncFormatter
41from scipy.optimize import minimize
43from lsst.utils.iteration import ensure_iterable
44from lsst.utils.plotting.figures import make_figure
46from .blockUtils import BlockParser
47from .efdUtils import (
48 COMMAND_ALIASES,
49 clipDataToEvent,
50 efdTimestampToAstropy,
51 getCommands,
52 getDayObsEndTime,
53 getDayObsForTime,
54 getDayObsStartTime,
55 getEfdData,
56 makeEfdClient,
57)
58from .enums import AxisMotionState, PowerState
59from .utils import dayObsIntToString, getCurrentDayObs_int
61if TYPE_CHECKING:
62 try:
63 from lsst_efd_client import EfdClient
64 except ImportError:
65 EfdClient = None
67__all__ = (
68 "TMAStateMachine",
69 "TMAEvent",
70 "TMAEventMaker",
71 "TMAState",
72 "AxisMotionState",
73 "PowerState",
74 "getSlewsFromEventList",
75 "getTracksFromEventList",
76 "getTorqueMaxima",
77 "filterBadValues",
78)
80# we don't want to use `None` for a no data sentinel because dict.get('key')
81# returns None if the key isn't present, and also we need to mark that the data
82# was queried for and no data was found, whereas the key not being present
83# means that we've not yet looked for the data.
84NO_DATA_SENTINEL = "NODATA"
86# The known time difference between the TMA demand position and the TMA
87# position when tracking. 20Hz data times three points = 150ms.
88TRACKING_RESIDUAL_TAIL_CLIP = -0.15 # seconds
90MOUNT_IMAGE_WARNING_LEVEL = 0.01 # this determines the colouring of the cells in the table, yellow for this
91MOUNT_IMAGE_BAD_LEVEL = 0.05 # and red for this
92HARDPOINT_FORCE_WARNING_LEVEL = 500.0 # for colouring of the cells in the table, yellow for this
93HARDPOINT_FORCE_BAD_LEVEL = 1000.0 # and red for this
96def getSlewsFromEventList(
97 events: list[TMAEvent],
98) -> list[TMAEvent]:
99 """Get the slew events from a list of TMAEvents.
101 Parameters
102 ----------
103 events : `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
104 The list of events to filter.
106 Returns
107 -------
108 events : `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
109 The filtered list of events.
110 """
111 return [e for e in events if e.type == TMAState.SLEWING]
114def getTracksFromEventList(
115 events: list[TMAEvent],
116) -> list[TMAEvent]:
117 """Get the tracking events from a list of TMAEvents.
119 Parameters
120 ----------
121 events : `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
122 The list of events to filter.
124 Returns
125 -------
126 events : `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
127 The filtered list of events.
128 """
129 return [e for e in events if e.type == TMAState.TRACKING]
132def getTorqueMaxima(table: pd.DataFrame) -> None:
133 """Print the maximum positive and negative azimuth and elevation torques.
135 Designed to be used with the table as downloaded from RubinTV.
137 Parameters
138 ----------
139 table : `pd.DataFrame`
140 The table of data to use, as generated by Rapid Analysis.
141 """
142 for axis in ["elevation", "azimuth"]:
143 col = f"Largest {axis} torque"
144 maxPos = np.argmax(table[col])
145 maxVal = table[col].iloc[maxPos]
146 print(f"Max positive {axis:9} torque during seqNum {maxPos:>4}: {maxVal / 1000:>7.1f}kNm")
147 minPos = np.argmin(table[col])
148 minVal = table[col].iloc[minPos]
149 print(f"Max negative {axis:9} torque during seqNum {minPos:>4}: {minVal / 1000:>7.1f}kNm")
150 return None
153def getAzimuthElevationDataForEvent(
154 client: EfdClient,
155 event: TMAEvent,
156 prePadding: float = 0,
157 postPadding: float = 0,
158 maxDeltaT: float = 1.0e-3,
159) -> tuple[pd.DataFrame, pd.DataFrame]:
160 """Get the data for the az/el telemetry topics for a given TMAEvent.
162 The error between the actual and demanded positions is calculated and added
163 to the dataframes in the az/elError columns. For TRACKING type events, this
164 error should be extremely close to zero, whereas for SLEWING type events,
165 this error represents the how far the TMA is from the demanded position,
166 and is therefore arbitrarily large, and tends to zero as the TMA get closer
167 to tracking the sky.
169 Parameters
170 ----------
171 client : `lsst_efd_client.efd_helper.EfdClient`
172 The EFD client to use.
173 event : `lsst.summit.utils.tmaUtils.TMAEvent`
174 The event to get the data for.
175 prePadding : `float`, optional
176 The amount of time to pad the event with before the start time, in
177 seconds.
178 postPadding : `float`, optional
179 The amount of time to pad the event with after the end time, in
180 seconds.
181 maxDeltaT : `float`, optional
182 The maximum allowed time adjustment to bring the median error to
183 zero, in seconds.
185 Returns
186 -------
187 azimuthData : `pd.DataFrame`
188 The azimuth data for the specified event.
189 elevationData : `pd.DataFrame`
190 The elevation data for the specified event.
191 """
193 def calcDeltaT(params, args):
194 # This calculates the deltaT needed
195 # to make the median(error) = 0
196 [values, valTimes, demand, demTimes] = args
197 [deltaT] = params
198 demandInterp = np.interp(valTimes, demTimes + deltaT, demand)
199 error = (values - demandInterp) * 3600
200 value = abs(np.median(error))
201 return value
203 azimuthData = getEfdData(
204 client,
205 "lsst.sal.MTMount.azimuth",
206 event=event,
207 prePadding=prePadding,
208 postPadding=postPadding,
209 raiseIfTopicNotInSchema=False,
210 )
211 elevationData = getEfdData(
212 client,
213 "lsst.sal.MTMount.elevation",
214 event=event,
215 prePadding=prePadding,
216 postPadding=postPadding,
217 raiseIfTopicNotInSchema=False,
218 )
220 azValues = np.asarray(azimuthData["actualPosition"])
221 azValTimes = np.asarray(azimuthData["actualPositionTimestamp"])
222 azDemand = np.asarray(azimuthData["demandPosition"])
223 azDemTimes = np.asarray(azimuthData["demandPositionTimestamp"])
224 elValues = np.asarray(elevationData["actualPosition"])
225 elValTimes = np.asarray(elevationData["actualPositionTimestamp"])
226 elDemand = np.asarray(elevationData["demandPosition"])
227 elDemTimes = np.asarray(elevationData["demandPositionTimestamp"])
229 # Calculate the deltaT needed to drive the median(error) to zero
230 args = [azValues, azValTimes, azDemand, azDemTimes]
231 x0 = [0.0]
232 result = minimize(calcDeltaT, x0, args=args, method="Powell", bounds=[(-maxDeltaT, maxDeltaT)])
233 deltaTAz = result.x[0]
235 args = [elValues, elValTimes, elDemand, elDemTimes]
236 x0 = [0.0]
237 result = minimize(calcDeltaT, x0, args=args, method="Powell", bounds=[(-maxDeltaT, maxDeltaT)])
238 deltaTEl = result.x[0]
240 azDemandInterp = np.interp(azValTimes, azDemTimes + deltaTAz, azDemand)
241 elDemandInterp = np.interp(elValTimes, elDemTimes + deltaTEl, elDemand)
243 azError = (azValues - azDemandInterp) * 3600
244 elError = (elValues - elDemandInterp) * 3600
246 azimuthData["azError"] = azError
247 elevationData["elError"] = elError
249 return azimuthData, elevationData
252def getM1M3HardpointDataForEvent(
253 client: EfdClient,
254 event: TMAEvent,
255 prePadding: float = 0,
256 postPadding: float = 0,
257) -> pd.DataFrame:
258 """Get the data for the M1M3 hardpoint telemetry for a given TMAEvent.
259 During slews, the M1M3 mirror moves in response to the accelerations. The
260 hardpoint forces need to stay within accepted limits to prevent long term
261 cracking on the mirror.
263 Parameters
264 ----------
265 client : `lsst_efd_client.efd_helper.EfdClient`
266 The EFD client to use.
267 event : `lsst.summit.utils.tmaUtils.TMAEvent`
268 The event to get the data for.
269 prePadding : `float`, optional
270 The amount of time to pad the event with before the start time, in
271 seconds.
272 postPadding : `float`, optional
273 The amount of time to pad the event with after the end time, in
274 seconds.
276 Returns
277 -------
278 hardpointData : `pd.DataFrame`
279 The hardpoint forces for the specified event.
280 """
281 hardpointData = getEfdData(
282 client,
283 "lsst.sal.MTM1M3.hardpointActuatorData",
284 event=event,
285 prePadding=prePadding,
286 postPadding=postPadding,
287 raiseIfTopicNotInSchema=False,
288 )
290 # Add a column with the absolute max force
291 maxForce = [
292 np.nanmax([abs(hardpointData.iloc[i][f"measuredForce{j}"]) for j in range(6)])
293 for i in range(len(hardpointData))
294 ]
296 hardpointData["maxForce"] = maxForce
298 return hardpointData
301def filterBadValues(values: list | np.ndarray, maxDelta: float = 0.1, maxConsecutiveValues: int = 3) -> int:
302 """Filter out bad values from a dataset, replacing them in-place.
304 This function replaces non-physical points in the dataset with an
305 extrapolation of the preceding two values. No more than 3 successive data
306 points are allowed to be replaced. Minimum length of the input is 3 points.
308 Parameters
309 ----------
310 values : `list` or `np.ndarray`
311 The dataset containing the values to be filtered.
312 maxDelta : `float`, optional
313 The maximum allowed difference between consecutive values. Values with
314 a difference greater than `maxDelta` will be considered as bad values
315 and replaced with an extrapolation.
316 maxConsecutiveValues : `int`, optional
317 The maximum number of consecutive values to replace. Defaults to 3.
319 Returns
320 -------
321 nBadPoints : `int`
322 The number of bad values that were replaced out.
323 """
324 # Find non-physical points and replace with extrapolation. No more than
325 # maxConsecutiveValues successive data points can be replaced.
326 badCounter = 0
327 consecutiveCounter = 0
329 log = logging.getLogger(__name__)
331 median = np.nanmedian(values)
332 # if either of the the first two points are more than maxDelta away from
333 # the median, replace them with the median
334 for i in range(2):
335 if abs(values[i] - median) > maxDelta:
336 log.warning(f"Replacing bad value of {values[i]} at index {i} with {median=}")
337 values[i] = median
338 badCounter += 1
340 # from the second element of the array, walk through and calculate the
341 # difference between each element and the previous one. If the difference
342 # is greater than maxDelta, replace the element with the average of the
343 # previous two known good values, i.e. ones which have not been replaced.
344 # if the first two points differ from the median by more than maxDelta,
345 # replace them with the median
346 lastGoodValue1 = values[1] # the most recent good value
347 lastGoodValue2 = values[0] # the second most recent good value
348 replacementValue = (lastGoodValue1 + lastGoodValue2) / 2.0 # in case we have to replace the first value
349 for i in range(2, len(values)):
350 if abs(values[i] - lastGoodValue1) >= maxDelta:
351 if consecutiveCounter < maxConsecutiveValues:
352 consecutiveCounter += 1
353 badCounter += 1
354 log.warning(f"Replacing value at index {i} with {replacementValue}")
355 values[i] = replacementValue
356 else:
357 log.warning(
358 f"More than 3 consecutive replacements at index {i}. Stopping replacements"
359 " until the next good value."
360 )
361 else:
362 lastGoodValue2 = lastGoodValue1
363 lastGoodValue1 = values[i]
364 replacementValue = (lastGoodValue1 + lastGoodValue2) / 2.0
365 consecutiveCounter = 0
366 return badCounter
369def plotEvent(
370 client: EfdClient,
371 event: TMAEvent,
372 fig: matplotlib.figure.Figure | None = None,
373 prePadding: float = 0,
374 postPadding: float = 0,
375 commands: dict[pd.Timestamp | datetime.datetime, str] = {},
376 azimuthData: pd.DataFrame | None = None,
377 elevationData: pd.DataFrame | None = None,
378 hardpointData: pd.DataFrame | None = None,
379 doFilterResiduals: bool = False,
380 maxDelta: float = 0.1,
381 metadataWriter: Callable | None = None,
382) -> matplotlib.figure.Figure:
383 """Plot the TMA axis positions over the course of a given TMAEvent.
385 Plots the axis motion profiles for the given event, with optional padding
386 at the start and end of the event. If the data is provided via the
387 azimuthData and elevationData parameters, it will be used, otherwise it
388 will be queried from the EFD. For slews, the M1M3 hardpoint forces are
389 also plotted. Similarly, it will query the EFD if the hardpoint forces are
390 not provided.
392 Optionally plots any commands issued during or around the event, if these
393 are supplied. Commands are supplied as a dictionary of tje times at which
394 the commands were issued, with values as command strings.
396 Due to a problem with the way the data is uploaded to the EFD, there are
397 occasional points in the tracking error plots that are very much larger
398 than the typical mount jitter. These points are unphysical, since it is not
399 possible for the mount to move that fast. We don't want these points, which
400 are not true mount problems, to distract from any real mount problems, and
401 these can be filtered out via the ``doFilterResiduals`` kwarg, which
402 replaces these non-physical points with an extrapolation of the average of
403 the preceding two known-good points. If the first two points are bad these
404 are replaced with the median of the dataset. The maximum difference between
405 the model and the actual data, in arcseconds, to allow before filtering a
406 data point can be set with the ``maxDelta`` kwarg.
408 Parameters
409 ----------
410 client : `lsst_efd_client.efd_helper.EfdClient`
411 The EFD client to use.
412 event : `lsst.summit.utils.tmaUtils.TMAEvent`
413 The event to plot.
414 fig : `matplotlib.figure.Figure`, optional
415 The figure to plot on. If not specified, a new figure will be created.
416 prePadding : `float`, optional
417 The amount of time to pad the event with before the start time, in
418 seconds.
419 postPadding : `float`, optional
420 The amount of time to pad the event with after the end time, in
421 seconds.
422 commands : `dict` [`pd.Timestamp`, `str`], or
423 `dict` [`datetime.datetime`, `str`], or optional
424 A dictionary of commands to plot on the figure. The keys are the times
425 at which a command was issued, and the value is the command string, as
426 returned by efdUtils.getCommands().
427 azimuthData : `pd.DataFrame`, optional
428 The azimuth data to plot. If not specified, it will be queried from the
429 EFD.
430 elevationData : `pd.DataFrame`, optional
431 The elevation data to plot. If not specified, it will be queried from
432 the EFD.
433 hardpointData : `pd.DataFrame`, optional
434 The hardpoint forces to plot. If not specified, it will be queried from
435 the EFD.
436 doFilterResiduals : 'bool', optional
437 Enables filtering of unphysical data points in the tracking residuals.
438 maxDelta : `float`, optional
439 The maximum difference between the model and the actual data, in
440 arcseconds, to allow before filtering the data point. Ignored if
441 ``doFilterResiduals`` is `False`.
442 metadataWriter : `callable`, optional
443 Should be a callable
444 ``lsst.rubintv.production.utils.writeMetadataShard`` function that has
445 had the path filled in with ``functools.patrial`` so that it will just
446 write out the data when called with the event's dayObs and a dictionary
447 containing the row data that should be written.
449 Returns
450 -------
451 fig : `matplotlib.figure.Figure`
452 The figure on which the plot was made.
453 """
455 def tickFormatter(value: float, tick_number: float) -> str:
456 # Convert the value to a string without subtracting large numbers
457 # tick_number is unused.
458 return f"{value:.2f}"
460 def getPlotTime(time: pd.Timestamp | Time | datetime.datetime) -> datetime.datetime:
461 """Get the right time to plot a point from the various time formats."""
462 match time:
463 case pd.Timestamp():
464 return time.to_pydatetime()
465 case astropy.time.Time():
466 return time.utc.datetime
467 case datetime.datetime():
468 return time
469 case _:
470 raise ValueError(f"Unknown type for commandTime: {type(time)}")
472 # plot any commands we might have
473 if not isinstance(commands, dict):
474 raise TypeError("commands must be a dict of command names with values as" " astropy.time.Time values")
476 if fig is None:
477 fig = make_figure(figsize=(10, 8))
478 log = logging.getLogger(__name__)
479 log.warning(
480 "Making new matplotlib figure - if this is in a loop you're going to have a bad time."
481 " Pass in a figure with fig = plt.figure(figsize=(10, 8)) to avoid this warning."
482 )
484 fig.clear()
485 ax1, ax1p5, ax2 = fig.subplots(
486 3, sharex=True, gridspec_kw={"wspace": 0, "hspace": 0, "height_ratios": [2.5, 1, 1]}
487 ) # type: ignore
489 if azimuthData is None or elevationData is None:
490 azimuthData, elevationData = getAzimuthElevationDataForEvent(
491 client, event, prePadding=prePadding, postPadding=postPadding
492 )
494 assert azimuthData is not None
495 assert elevationData is not None
497 # Use the native color cycle for the lines. Because they're on different
498 # axes they don't cycle by themselves
499 lineColors = [p["color"] for p in plt.rcParams["axes.prop_cycle"]]
500 nColors = len(lineColors)
501 colorCounter = 0
503 ax1.plot(azimuthData["actualPosition"], label="Azimuth position", c=lineColors[colorCounter % nColors])
504 colorCounter += 1
505 ax1.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
506 ax1.set_ylabel("Azimuth (degrees)")
508 ax1_twin = ax1.twinx()
509 ax1_twin.plot(
510 elevationData["actualPosition"], label="Elevation position", c=lineColors[colorCounter % nColors]
511 )
512 colorCounter += 1
513 ax1_twin.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
514 ax1_twin.set_ylabel("Elevation (degrees)")
515 ax1.set_xticks([]) # remove x tick labels on the hidden upper x-axis
517 ax2_twin = ax2.twinx()
518 ax2.plot(azimuthData["actualTorque"], label="Azimuth torque", c=lineColors[colorCounter % nColors])
519 colorCounter += 1
520 ax2_twin.plot(
521 elevationData["actualTorque"], label="Elevation torque", c=lineColors[colorCounter % nColors]
522 )
523 colorCounter += 1
524 ax2.set_ylabel("Azimuth torque (Nm)")
525 ax2_twin.set_ylabel("Elevation torque (Nm)")
526 ax2.set_xlabel("Time (UTC)") # yes, it really is UTC, matplotlib converts this automatically!
528 if event.type == TMAState.TRACKING:
529 # returns a copy
530 clippedAzimuthData = clipDataToEvent(azimuthData, event, postPadding=TRACKING_RESIDUAL_TAIL_CLIP)
531 clippedElevationData = clipDataToEvent(elevationData, event, postPadding=TRACKING_RESIDUAL_TAIL_CLIP)
533 azError = clippedAzimuthData["azError"].to_numpy()
534 elError = clippedElevationData["elError"].to_numpy()
535 elVals = clippedElevationData["actualPosition"].to_numpy()
537 if doFilterResiduals:
538 # Filtering out bad values
539 nReplacedAz = filterBadValues(azError, maxDelta)
540 nReplacedEl = filterBadValues(elError, maxDelta)
541 clippedAzimuthData["azError"] = azError
542 clippedElevationData["elError"] = elError
543 replacementMsg = (
544 f"{nReplacedAz} bad azimuth values and {nReplacedEl} bad elevation values were replaced"
545 )
546 else:
547 # currently unused, but needed for mypy, and better than an
548 # alternative in case someone redesigns that part later
549 replacementMsg = "No values were replaced because doFilterResiduals is False"
551 # Calculate RMS
552 az_rms = np.sqrt(np.mean(azError * azError))
553 el_rms = np.sqrt(np.mean(elError * elError))
555 # Calculate Image impact RMS
556 # We are less sensitive to Az errors near the zenith
557 image_az_rms = az_rms * np.cos(elVals[0] * np.pi / 180.0)
558 image_el_rms = el_rms
559 image_impact_rms = np.sqrt(image_az_rms**2 + image_el_rms**2)
561 ax1p5.plot(
562 clippedAzimuthData["azError"],
563 label="Azimuth tracking error",
564 c=lineColors[colorCounter % nColors],
565 )
566 colorCounter += 1
567 ax1p5.plot(
568 clippedElevationData["elError"],
569 label="Elevation tracking error",
570 c=lineColors[colorCounter % nColors],
571 )
572 colorCounter += 1
573 ax1p5.axhline(0.01, ls="-.", color="black")
574 ax1p5.axhline(-0.01, ls="-.", color="black")
575 ax1p5.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
576 ax1p5.set_ylabel("Tracking error (arcsec)")
577 ax1p5.set_xticks([]) # remove x tick labels on the hidden upper x-axis
578 ax1p5.set_ylim(-0.05, 0.05)
579 ax1p5.set_yticks([-0.04, -0.02, 0.0, 0.02, 0.04])
580 ax1p5.legend()
581 ax1p5.text(0.1, 0.9, f"Image impact RMS = {image_impact_rms:.3f} arcsec", transform=ax1p5.transAxes)
582 if doFilterResiduals:
583 ax1p5.text(
584 0.1,
585 0.8,
586 replacementMsg,
587 transform=ax1p5.transAxes,
588 )
589 if metadataWriter is not None:
590 md = {"Tracking image impact": f"{image_impact_rms:.3f}"}
591 flagKey = "_Tracking image impact"
592 if image_impact_rms > MOUNT_IMAGE_BAD_LEVEL:
593 md.update({flagKey: "bad"})
594 elif image_impact_rms > MOUNT_IMAGE_WARNING_LEVEL:
595 md.update({flagKey: "warning"})
597 rowData = {event.seqNum: md}
598 metadataWriter(dayObs=event.dayObs, mdDict=rowData)
600 elif event.type == TMAState.SLEWING:
601 if hardpointData is None:
602 hardpointData = getM1M3HardpointDataForEvent(
603 client, event, prePadding=prePadding, postPadding=postPadding
604 )
606 if not hardpointData.empty:
607 # Calculate Max hardpoint force
608 maxForce = np.nanmax(hardpointData["maxForce"].to_numpy())
610 for hp_index in range(6):
611 ax1p5.plot(
612 hardpointData[f"measuredForce{hp_index}"],
613 label=f"HP_{hp_index}",
614 c=lineColors[colorCounter % nColors],
615 )
616 colorCounter += 1
617 ax1p5.axhline(500.0, ls="-.", color="goldenrod")
618 ax1p5.axhline(-500.0, ls="-.", color="goldenrod")
619 ax1p5.axhline(1000.0, ls="-.", color="red")
620 ax1p5.axhline(-1000.0, ls="-.", color="red")
621 ax1p5.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
622 ax1p5.set_ylabel("HP force (N)")
623 ax1p5.set_xticks([]) # remove x tick labels on the hidden upper x-axis
624 if maxForce > 2000.0:
625 ax1p5.set_ylim(-4000.0, 4000.0)
626 ax1p5.set_yticks([-2000, 0.0, 2000])
627 elif maxForce > 1000.0:
628 ax1p5.set_ylim(-2000.0, 2000.0)
629 ax1p5.set_yticks([-1000, 0.0, 1000])
630 else:
631 ax1p5.set_ylim(-1000.0, 1000.0)
632 ax1p5.set_yticks([-500, 0.0, 500])
633 ax1p5.legend()
634 ax1p5.text(0.1, 0.9, f"Max HP force = {maxForce:.1f} N", transform=ax1p5.transAxes)
636 if metadataWriter is not None:
637 md = {"HP Max force": f"{maxForce:.2f}"}
638 flagKey = "_HP Max force"
639 if maxForce > HARDPOINT_FORCE_BAD_LEVEL:
640 md.update({flagKey: "bad"})
641 elif maxForce > HARDPOINT_FORCE_WARNING_LEVEL:
642 md.update({flagKey: "warning"})
644 rowData = {event.seqNum: md}
645 metadataWriter(dayObs=event.dayObs, mdDict=rowData)
647 if prePadding or postPadding:
648 # note the conversion to utc because the x-axis from the dataframe
649 # already got automagically converted when plotting before, so this is
650 # necessary for things to line up
651 ax1_twin.axvline(event.begin.utc.datetime, c="k", ls="--", alpha=0.5, label="Event begin/end")
652 ax1_twin.axvline(event.end.utc.datetime, c="k", ls="--", alpha=0.5)
653 # extend lines down across lower plot, but do not re-add label
654 ax2_twin.axvline(event.begin.utc.datetime, c="k", ls="--", alpha=0.5)
655 ax2_twin.axvline(event.end.utc.datetime, c="k", ls="--", alpha=0.5)
656 ax1p5.axvline(event.begin.utc.datetime, c="k", ls="--", alpha=0.5)
657 ax1p5.axvline(event.end.utc.datetime, c="k", ls="--", alpha=0.5)
659 for commandTime, command in commands.items():
660 plotTime = getPlotTime(commandTime)
661 ax1_twin.axvline(
662 plotTime, c=lineColors[colorCounter % nColors], ls="--", alpha=0.75, label=f"{command}"
663 )
664 # extend lines down across lower plot, but do not re-add label
665 ax2_twin.axvline(plotTime, c=lineColors[colorCounter % nColors], ls="--", alpha=0.75)
666 ax1p5.axvline(plotTime, c=lineColors[colorCounter % nColors], ls="--", alpha=0.75)
667 colorCounter += 1
669 # combine the legends and put inside the plot
670 handles1a, labels1a = ax1.get_legend_handles_labels()
671 handles1b, labels1b = ax1_twin.get_legend_handles_labels()
672 handles2a, labels2a = ax2.get_legend_handles_labels()
673 handles2b, labels2b = ax2_twin.get_legend_handles_labels()
675 handles = handles1a + handles1b + handles2a + handles2b
676 labels = labels1a + labels1b + labels2a + labels2b
677 # ax2 is "in front" of ax1 because it has the vlines plotted on it, and
678 # vlines are on ax2 so that they appear at the bottom of the legend, so
679 # make sure to plot the legend on ax2, otherwise the vlines will go on top
680 # of the otherwise-opaque legend.
681 ax1_twin.legend(handles, labels, facecolor="white", framealpha=1)
683 # Add title with the event name, type etc
684 dayObsStr = dayObsIntToString(event.dayObs)
685 title = (
686 # top line is the event title, the details go on the line below
687 f"{dayObsStr} - seqNum {event.seqNum} (version {event.version})"
688 f"\nDuration = {event.duration:.2f}s"
689 f" Event type: {event.type.name}"
690 f" End reason: {event.endReason.name}"
691 )
692 ax1_twin.set_title(title)
694 # Moving this block to the end recovers the time ticks and ticklabels.
695 # put the ticks at an angle, and right align with the tick marks
696 ax2.set_xticks(ax2.get_xticks()) # needed to supress a user warning
697 xlabels = ax2.get_xticks()
698 ax2.set_xticklabels(xlabels, rotation=40, ha="right")
699 ax2.xaxis.set_major_locator(mdates.AutoDateLocator())
700 ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
702 return fig
705def getCommandsDuringEvent(
706 client: EfdClient,
707 event: TMAEvent,
708 commands: Iterable[str] = ("raDecTarget",),
709 prePadding: float = 0,
710 postPadding: float = 0,
711 timeFormat: str = "python",
712 log: logging.Logger | None = None,
713 doLog: bool = True,
714) -> dict[datetime.datetime, str]:
715 """Get the commands issued during an event.
717 Get the times at which the specified commands were issued during the event.
719 Parameters
720 ----------
721 client : `lsst_efd_client.efd_helper.EfdClient`
722 The EFD client to use.
723 event : `lsst.summit.utils.tmaUtils.TMAEvent`
724 The event to plot.
725 commands : `list` of `str`, optional
726 The commands or command aliases to look for. Defaults to
727 ['raDecTarget'].
728 prePadding : `float`, optional
729 The amount of time to pad the event with before the start time, in
730 seconds.
731 postPadding : `float`, optional
732 The amount of time to pad the event with after the end time, in
733 seconds.
734 timeFormat : `str`, optional
735 One of 'pandas' or 'astropy' or 'python'. If 'pandas', the dictionary
736 keys will be pandas timestamps, if 'astropy' they will be astropy times
737 and if 'python' they will be python datetimes.
738 log : `logging.Logger`, optional
739 The logger to use. If not specified, a new logger will be created if
740 needed.
741 doLog : `bool`, optional
742 Whether to log messages. Defaults to True.
744 Returns
745 -------
746 commandTimes : `dict` [`time`, `str`]
747 A dictionary of the times at which the commands where issued. The type
748 that `time` takes is determined by the format key, and defaults to
749 python datetime.
750 """
751 commands = list(ensure_iterable(commands))
752 fullCommands = [c if c not in COMMAND_ALIASES else COMMAND_ALIASES[c] for c in commands]
753 del commands # make sure we always use their full names
755 commandTimes = getCommands(
756 client,
757 fullCommands,
758 begin=event.begin,
759 end=event.end,
760 prePadding=prePadding,
761 postPadding=postPadding,
762 timeFormat=timeFormat,
763 )
765 if not commandTimes and doLog:
766 log = logging.getLogger(__name__)
767 log.info(f"Found no commands in {fullCommands} issued during event {event.seqNum}")
769 return commandTimes
772def _initializeTma(tma: TMAStateMachine) -> None:
773 """Helper function to turn a TMA into a valid state for testing.
775 Do not call directly in normal usage or code, as this just arbitrarily
776 sets values to make the TMA valid.
778 Parameters
779 ----------
780 tma : `lsst.summit.utils.tmaUtils.TMAStateMachine`
781 The TMA state machine model to initialize.
782 """
783 tma._parts["azimuthInPosition"] = False
784 tma._parts["azimuthMotionState"] = AxisMotionState.STOPPED
785 tma._parts["azimuthSystemState"] = PowerState.ON
786 tma._parts["elevationInPosition"] = False
787 tma._parts["elevationMotionState"] = AxisMotionState.STOPPED
788 tma._parts["elevationSystemState"] = PowerState.ON
791@dataclass(kw_only=True, frozen=True)
792class TMAEvent:
793 """A movement event for the TMA.
795 Contains the dayObs on which the event occured, using the standard
796 observatory definition of the dayObs, and the sequence number of the event,
797 which is unique for each event on a given dayObs.
799 The event type can be either 'SLEWING' or 'TRACKING', defined as:
800 - SLEWING: some part of the TMA is in motion
801 - TRACKING: both axes are in position and tracking the sky
803 The end reason can be 'STOPPED', 'TRACKING', 'FAULT', 'SLEWING', or 'OFF'.
804 - SLEWING: The previous event was a TRACKING event, and one or more of
805 the TMA components either stopped being in position, or stopped
806 moving, or went into fault, or was turned off, and hence we are now
807 only slewing and no longer tracking the sky.
808 - TRACKING: the TMA started tracking the sky when it wasn't previously.
809 Usualy this would always be preceded by directly by a SLEWING
810 event, but this is not strictly true, as the EUI seems to be able
811 to make the TMA start tracking the sky without slewing first.
812 - STOPPED: the components of the TMA transitioned to the STOPPED state.
813 - FAULT: the TMA went into fault.
814 - OFF: the TMA components were turned off.
816 Note that this class is not intended to be instantiated directly, but
817 rather to be returned by the ``TMAEventMaker.getEvents()`` function.
819 Parameters
820 ----------
821 dayObs : `int`
822 The dayObs on which the event occured.
823 seqNum : `int`
824 The sequence number of the event,
825 type : `lsst.summit.utils.tmaUtils.TMAState`
826 The type of the event, either 'SLEWING' or 'TRACKING'.
827 endReason : `lsst.summit.utils.tmaUtils.TMAState`
828 The reason the event ended, either 'STOPPED', 'TRACKING', 'FAULT',
829 'SLEWING', or 'OFF'.
830 duration : `float`
831 The duration of the event, in seconds.
832 begin : `astropy.time.Time`
833 The time the event began.
834 end : `astropy.time.Time`
835 The time the event ended.
836 blockInfos : `list` of `lsst.summit.utils.tmaUtils.BlockInfo`, or `None`
837 The block infomation, if any, relating to the event. Could be `None`,
838 or one or more block informations.
839 version : `int`
840 The version of the TMAEvent class. Equality between events is only
841 valid for a given version of the class. If the class definition
842 changes, the time ranges can change, and hence the equality between
843 events is ``False``.
844 _startRow : `int`
845 The first row in the merged EFD data which is part of the event.
846 _endRow : `int`
847 The last row in the merged EFD data which is part of the event.
848 """
850 dayObs: int
851 seqNum: int
852 type: TMAState # can be 'SLEWING', 'TRACKING'
853 endReason: TMAState # can be 'STOPPED', 'TRACKING', 'FAULT', 'SLEWING', 'OFF'
854 duration: float # seconds
855 begin: Time
856 end: Time
857 blockInfos: list = field(default_factory=list)
858 version: int = 0 # update this number any time a code change which could change event definitions is made
859 _startRow: int
860 _endRow: int
862 def __lt__(self, other: TMAEvent) -> bool:
863 if self.version != other.version:
864 raise ValueError(
865 f"Cannot compare TMAEvents with different versions: {self.version} != {other.version}"
866 )
867 if self.dayObs < other.dayObs:
868 return True
869 elif self.dayObs == other.dayObs:
870 return self.seqNum < other.seqNum
871 return False
873 def __eq__(self, other: object) -> bool:
874 if not isinstance(other, TMAEvent):
875 return NotImplemented
876 if self.version != other.version:
877 return False
878 return (
879 self.dayObs == other.dayObs
880 and self.seqNum == other.seqNum
881 and self.type == other.type
882 and self.endReason == other.endReason
883 and self.duration == other.duration
884 and self.begin == other.begin
885 and self.end == other.end
886 # deliberately do not compare on blockInfos, so skip explicitly
887 # and self.blockInfos == other.blockInfos
888 and self._startRow == other._startRow
889 and self._endRow == other._endRow
890 )
892 def __repr__(self) -> str:
893 return (
894 f"TMAEvent(dayObs={self.dayObs}, seqNum={self.seqNum}, type={self.type!r},"
895 f" endReason={self.endReason!r}, duration={self.duration}, begin={self.begin!r},"
896 f" end={self.end!r}"
897 )
899 def __hash__(self) -> int:
900 # deliberately don't hash the blockInfos here, as they are not
901 # a core part of the event itself, and are listy and cause problems
902 return hash(
903 (
904 self.dayObs,
905 self.seqNum,
906 self.type,
907 self.endReason,
908 self.duration,
909 self.begin,
910 self.end,
911 self.version,
912 self._startRow,
913 self._endRow,
914 )
915 )
917 def _ipython_display_(self) -> None:
918 print(self.__str__())
920 def __str__(self) -> str:
921 def indent(string: str) -> str:
922 return "\n" + "\n".join([" " + s for s in string.splitlines()])
924 blockInfoStr = "None"
925 if self.blockInfos is not None:
926 blockInfoStr = "".join(indent(str(i)) for i in self.blockInfos)
928 return (
929 f"dayObs: {self.dayObs}\n"
930 f"seqNum: {self.seqNum}\n"
931 f"type: {self.type.name}\n"
932 f"endReason: {self.endReason.name}\n"
933 f"duration: {self.duration}\n"
934 f"begin: {self.begin!r}\n"
935 f"end: {self.end!r}\n"
936 f"blockInfos: {blockInfoStr}"
937 )
939 def associatedWith(
940 self,
941 block: int | str | None = None,
942 blockSeqNum: int | None = None,
943 ticket: str | None = None,
944 salIndex: int | None = None,
945 ) -> bool:
946 """Check whether an event is associated with a set of parameters.
948 Check if an event is associated with a specific block and/or ticket
949 and/or salIndex. All specified parameters must match for the function
950 to return True. If checking if an event is in a block, the blockSeqNum
951 can also be specified to identify events which related to a given
952 running the specified block.
954 Parameters
955 ----------
956 block : `int`, optional
957 The block number to check for.
958 blockSeqNum : `int`, optional
959 The block sequence number to check for, if the block is specified.
960 ticket : `str`, optional
961 The ticket number to check for.
962 salIndex : `int`, optional
963 The salIndex to check for.
965 Returns
966 -------
967 relates : `bool`
968 Whether the event is associated with the specified block, ticket,
969 and salIndex.
970 """
971 if isinstance(block, int):
972 block = str(block)
974 if all([block is None, ticket is None, salIndex is None]):
975 raise ValueError("Must specify at least one of block, ticket, or salIndex")
977 if blockSeqNum is not None and block is None:
978 raise ValueError("block must be specified if blockSeqNum is specified")
980 for blockInfo in self.blockInfos:
981 # "X is None or" is used for each parameter to allow it to be None
982 # in the kwargs
983 blockMatches = False
984 if block is not None:
985 if blockSeqNum is None and blockInfo.blockNumber == block:
986 blockMatches = True
987 elif (
988 blockSeqNum is not None
989 and blockInfo.blockNumber == block
990 and blockInfo.seqNum == blockSeqNum
991 ):
992 blockMatches = True
993 else:
994 blockMatches = True # no block specified at all, so it matches
996 salIndexMatches = salIndex is None or salIndex in blockInfo.salIndices
997 ticketMatches = ticket is None or ticket in blockInfo.tickets
999 if blockMatches and salIndexMatches and ticketMatches:
1000 return True
1002 return False
1005class TMAState(enum.IntEnum):
1006 """Overall state of the TMA.
1008 States are defined as follows:
1010 UNINITIALIZED
1011 We have not yet got data for all relevant components, so the overall
1012 state is undefined.
1013 STOPPED
1014 All components are on, and none are moving.
1015 TRACKING
1016 We are tracking the sky.
1017 SLEWING
1018 One or more components are moving, and one or more are not tracking the
1019 sky. This should probably be called MOVING, as it includes: slewing,
1020 MOVING_POINT_TO_POINT, and JOGGING.
1021 FAULT
1022 All (if engineeringMode) or any (if not engineeringMode) components are
1023 in fault.
1024 OFF
1025 All components are off.
1026 """
1028 UNINITIALIZED = -1
1029 STOPPED = 0
1030 TRACKING = 1
1031 SLEWING = 2
1032 FAULT = 3
1033 OFF = 4
1035 def __repr__(self) -> str:
1036 return f"TMAState.{self.name}"
1039def getAxisAndType(rowFor: str) -> tuple[str, str]:
1040 """Get the axis the data relates to, and the type of data it contains.
1042 Parameters
1043 ----------
1044 rowFor : `str`
1045 The column in the dataframe denoting what this row is for, e.g.
1046 "elevationMotionState" or "azimuthInPosition", etc.
1048 Returns
1049 -------
1050 axis : `str`
1051 The axis the row is for, e.g. "azimuth", "elevation".
1052 rowType : `str`
1053 The type of the row, e.g. "MotionState", "SystemState", "InPosition".
1054 """
1055 regex = r"(azimuth|elevation)(InPosition|MotionState|SystemState)$" # matches the end of the line
1056 matches = re.search(regex, rowFor)
1057 if matches is None:
1058 raise ValueError(f"Could not parse axis and rowType from {rowFor=}")
1059 axis = matches.group(1)
1060 rowType = matches.group(2)
1062 assert rowFor.endswith(f"{axis}{rowType}")
1063 return axis, rowType
1066class ListViewOfDict:
1067 """A class to allow making lists which contain references to an underlying
1068 dictionary.
1070 Normally, making a list of items from a dictionary would make a copy of the
1071 items, but this class allows making a list which contains references to the
1072 underlying dictionary items themselves. This is useful for making a list of
1073 components, such that they can be manipulated in their logical sets.
1074 """
1076 def __init__(self, underlyingDictionary: dict, keysToLink: list):
1077 self.dictionary = underlyingDictionary
1078 self.keys = keysToLink
1080 def __getitem__(self, index: Any) -> Any:
1081 return self.dictionary[self.keys[index]]
1083 def __setitem__(self, index: int, value: Any) -> None:
1084 self.dictionary[self.keys[index]] = value
1086 def __len__(self) -> int:
1087 return len(self.keys)
1089 def __iter__(self) -> Any:
1090 for key in self.keys:
1091 yield self.dictionary[key]
1094class TMAStateMachine:
1095 """A state machine model of the TMA.
1097 Note that this is currently only implemented for the azimuth and elevation
1098 axes, but will be extended to include the rotator in the future.
1100 Note that when used for event generation, changing ``engineeringMode`` to
1101 False might change the resulting list of events, and that if the TMA moves
1102 with some axis in fault, then these events will be missed. It is therefore
1103 thought that ``engineeringMode=True`` should always be used when generating
1104 events. The option, however, is there for completeness, as this will be
1105 useful for knowing is the CSC would consider the TMA to be in fault in the
1106 general case.
1108 Parameters
1109 ----------
1110 engineeringMode : `bool`, optional
1111 Whether the TMA is in engineering mode. Defaults to True. If False,
1112 then the TMA will be in fault if any component is in fault. If True,
1113 then the TMA will be in fault only if all components are in fault.
1114 debug : `bool`, optional
1115 Whether to log debug messages. Defaults to False.
1116 """
1118 _UNINITIALIZED_VALUE: int = -999
1120 def __init__(self, engineeringMode: bool = True, debug: bool = False):
1121 self.engineeringMode = engineeringMode
1122 self.log = logging.getLogger("lsst.summit.utils.tmaUtils.TMA")
1123 if debug:
1124 self.log.level = logging.DEBUG
1125 self._mostRecentRowTime = -1
1127 # the actual components of the TMA
1128 self._parts = {
1129 "azimuthInPosition": self._UNINITIALIZED_VALUE,
1130 "azimuthMotionState": self._UNINITIALIZED_VALUE,
1131 "azimuthSystemState": self._UNINITIALIZED_VALUE,
1132 "elevationInPosition": self._UNINITIALIZED_VALUE,
1133 "elevationMotionState": self._UNINITIALIZED_VALUE,
1134 "elevationSystemState": self._UNINITIALIZED_VALUE,
1135 }
1136 systemKeys = ["azimuthSystemState", "elevationSystemState"]
1137 positionKeys = ["azimuthInPosition", "elevationInPosition"]
1138 motionKeys = ["azimuthMotionState", "elevationMotionState"]
1140 # references to the _parts as conceptual groupings
1141 self.system = ListViewOfDict(self._parts, systemKeys)
1142 self.motion = ListViewOfDict(self._parts, motionKeys)
1143 self.inPosition = ListViewOfDict(self._parts, positionKeys)
1145 # tuples of states for state collapsing. Note that STOP_LIKE +
1146 # MOVING_LIKE must cover the full set of AxisMotionState enums
1147 self.STOP_LIKE = (AxisMotionState.STOPPING, AxisMotionState.STOPPED, AxisMotionState.TRACKING_PAUSED)
1148 self.MOVING_LIKE = (
1149 AxisMotionState.MOVING_POINT_TO_POINT,
1150 AxisMotionState.JOGGING,
1151 AxisMotionState.TRACKING,
1152 )
1153 # Likewise, ON_LIKE + OFF_LIKE must cover the full set of PowerState
1154 # enums
1155 self.OFF_LIKE = (PowerState.OFF, PowerState.TURNING_OFF)
1156 self.ON_LIKE = (PowerState.ON, PowerState.TURNING_ON)
1157 self.FAULT_LIKE = (PowerState.FAULT,) # note the trailing comma - this must be an iterable
1159 def apply(self, row: pd.Series) -> None:
1160 """Apply a row of data to the TMA state.
1162 Checks that the row contains data for a later time than any data
1163 previously applied, and applies the relevant column entry to the
1164 relevant component.
1166 Parameters
1167 ----------
1168 row : `pd.Series`
1169 The row of data to apply to the state machine.
1170 """
1171 timestamp = row["private_efdStamp"]
1172 if timestamp < self._mostRecentRowTime: # NB equals is OK, technically, though it never happens
1173 raise ValueError(
1174 "TMA evolution must be monotonic increasing in time, tried to apply a row which"
1175 " predates the most previous one"
1176 )
1177 self._mostRecentRowTime = timestamp
1179 rowFor = row["rowFor"] # e.g. elevationMotionState
1180 axis, rowType = getAxisAndType(rowFor) # e.g. elevation, MotionState
1181 value = self._getRowPayload(row, rowType, rowFor)
1182 self.log.debug(f"Setting {rowFor} to {repr(value)}")
1183 self._parts[rowFor] = value
1184 try:
1185 # touch the state property as this executes the sieving, to make
1186 # sure we don't fall through the sieve at any point in time
1187 _ = self.state
1188 except RuntimeError as e:
1189 # improve error reporting, but always reraise this, as this is a
1190 # full-blown failure
1191 raise RuntimeError(f"Failed to apply {value} to {axis}{rowType} with state {self._parts}") from e
1193 def _getRowPayload(
1194 self, row: pd.Series, rowType: str, rowFor: str
1195 ) -> bool | AxisMotionState | PowerState:
1196 """Get the relevant value from the row.
1198 Given the row, and which component it relates to, get the relevant
1199 value, as a bool or cast to the appropriate enum class.
1201 Parameters
1202 ----------
1203 row : `pd.Series`
1204 The row of data from the dataframe.
1205 rowType : `str`
1206 The type of the row, e.g. "MotionState", "SystemState",
1207 "InPosition".
1208 rowFor : `str`
1209 The component the row is for, e.g. "azimuth", "elevation".
1211 Returns
1212 -------
1213 value : `bool` or `enum`
1214 The value of the row, as a bool or enum, depending on the
1215 component, cast to the appropriate enum class or bool.
1216 """
1217 match rowType:
1218 case "MotionState":
1219 value = row[f"state_{rowFor}"]
1220 return AxisMotionState(value)
1221 case "SystemState":
1222 value = row[f"powerState_{rowFor}"]
1223 return PowerState(value)
1224 case "InPosition":
1225 value = row[f"inPosition_{rowFor}"]
1226 return bool(value)
1227 case _:
1228 raise ValueError(f"Failed to get row payload with {rowType=} and {row=}")
1230 @property
1231 def _isValid(self) -> bool:
1232 """Has the TMA had a value applied to all its components?
1234 If any component has not yet had a value applied, the TMA is not valid,
1235 as those components will be in an unknown state.
1237 Returns
1238 -------
1239 isValid : `bool`
1240 Whether the TMA is fully initialized.
1241 """
1242 return not any([v == self._UNINITIALIZED_VALUE for v in self._parts.values()])
1244 # state inspection properties - a high level way of inspecting the state as
1245 # an API
1246 @property
1247 def isMoving(self) -> bool:
1248 return self.state in [TMAState.TRACKING, TMAState.SLEWING]
1250 @property
1251 def isNotMoving(self) -> bool:
1252 return not self.isMoving
1254 @property
1255 def isTracking(self) -> bool:
1256 return self.state == TMAState.TRACKING
1258 @property
1259 def isSlewing(self) -> bool:
1260 return self.state == TMAState.SLEWING
1262 @property
1263 def canMove(self) -> bool:
1264 badStates = [PowerState.OFF, PowerState.TURNING_OFF, PowerState.FAULT, PowerState.UNKNOWN]
1265 return bool(
1266 self._isValid
1267 and self._parts["azimuthSystemState"] not in badStates
1268 and self._parts["elevationSystemState"] not in badStates
1269 )
1271 # Axis inspection properties, designed for internal use. These return
1272 # iterables so that they can be used in any() and all() calls, which make
1273 # the logic much easier to read, e.g. to see if anything is moving, we can
1274 # write `if not any(_axisInMotion):`
1275 @property
1276 def _axesInFault(self) -> list[bool]:
1277 return [x in self.FAULT_LIKE for x in self.system]
1279 @property
1280 def _axesOff(self) -> list[bool]:
1281 return [x in self.OFF_LIKE for x in self.system]
1283 @property
1284 def _axesOn(self) -> list[bool]:
1285 return [not x for x in self._axesOn]
1287 @property
1288 def _axesInMotion(self) -> list[bool]:
1289 return [x in self.MOVING_LIKE for x in self.motion]
1291 @property
1292 def _axesTRACKING(self) -> list[bool]:
1293 """Note this is deliberately named _axesTRACKING and not _axesTracking
1294 to make it clear that this is the AxisMotionState type of TRACKING and
1295 not the normal conceptual notion of tracking (the sky, i.e. as opposed
1296 to slewing).
1297 """
1298 return [x == AxisMotionState.TRACKING for x in self.motion]
1300 @property
1301 def _axesInPosition(self) -> list[bool]:
1302 return [x is True for x in self.inPosition]
1304 @property
1305 def state(self) -> TMAState:
1306 """The overall state of the TMA.
1308 Note that this is both a property, and also the method which applies
1309 the logic sieve to determine the state at a given point in time.
1311 Returns
1312 -------
1313 state : `lsst.summit.utils.tmaUtils.TMAState`
1314 The overall state of the TMA.
1315 """
1316 # first, check we're valid, and if not, return UNINITIALIZED state, as
1317 # things are unknown
1318 if not self._isValid:
1319 return TMAState.UNINITIALIZED
1321 # if we're not in engineering mode, i.e. we're under normal CSC
1322 # control, then if anything is in fault, we're in fault. If we're
1323 # engineering then some axes will move when others are in fault
1324 if not self.engineeringMode:
1325 if any(self._axesInFault):
1326 return TMAState.FAULT
1327 else:
1328 # we're in engineering mode, so return fault state if ALL are in
1329 # fault
1330 if all(self._axesInFault):
1331 return TMAState.FAULT
1333 # if all axes are off, the TMA is OFF
1334 if all(self._axesOff):
1335 return TMAState.OFF
1337 # we know we're valid and at least some axes are not off, so see if
1338 # we're in motion if no axes are moving, we're stopped
1339 if not any(self._axesInMotion):
1340 return TMAState.STOPPED
1342 # now we know we're initialized, and that at least one axis is moving
1343 # so check axes for motion and in position. If all axes are tracking
1344 # and all are in position, we're tracking the sky
1345 if all(self._axesTRACKING) and all(self._axesInPosition):
1346 return TMAState.TRACKING
1348 # we now know explicitly that not everything is in position, so we no
1349 # longer need to check that. We do actually know that something is in
1350 # motion, but confirm that's the case and return SLEWING
1351 if any(self._axesInMotion):
1352 return TMAState.SLEWING
1354 # if we want to differentiate between MOVING_POINT_TO_POINT moves,
1355 # JOGGING moves and regular slews, the logic in the step above needs to
1356 # be changed and the new steps added here.
1358 raise RuntimeError("State error: fell through the state sieve - rewrite your logic!")
1361class TMAEventMaker:
1362 """A class to create per-dayObs TMAEvents for the TMA's movements.
1364 If this class is being used in tests, make sure to pass the EFD client in,
1365 and create it with `makeEfdClient(testing=True)`. This ensures that the
1366 USDF EFD is "used" as this is the EFD which has the recorded data available
1367 in the test suite via `vcr`.
1369 Example usage:
1370 >>> dayObs = 20230630
1371 >>> eventMaker = TMAEventMaker()
1372 >>> events = eventMaker.getEvents(dayObs)
1373 >>> print(f'Found {len(events)} for {dayObs=}')
1375 Parameters
1376 ----------
1377 client : `lsst_efd_client.efd_helper.EfdClient`, optional
1378 The EFD client to use, created if not provided.
1379 """
1381 # the topics which need logical combination to determine the overall mount
1382 # state. Will need updating as new components are added to the system.
1384 # relevant column: 'state'
1385 _movingComponents = [
1386 "lsst.sal.MTMount.logevent_azimuthMotionState",
1387 "lsst.sal.MTMount.logevent_elevationMotionState",
1388 ]
1390 # relevant column: 'inPosition'
1391 _inPositionComponents = [
1392 "lsst.sal.MTMount.logevent_azimuthInPosition",
1393 "lsst.sal.MTMount.logevent_elevationInPosition",
1394 ]
1396 # the components which, if in fault, put the TMA into fault
1397 # relevant column: 'powerState'
1398 _stateComponents = [
1399 "lsst.sal.MTMount.logevent_azimuthSystemState",
1400 "lsst.sal.MTMount.logevent_elevationSystemState",
1401 ]
1403 def __init__(self, client: EfdClient | None = None):
1404 if client is not None:
1405 self.client = client
1406 else:
1407 self.client = makeEfdClient()
1408 self.log = logging.getLogger(__name__)
1409 self._data: dict = {}
1411 @dataclass(frozen=True)
1412 class ParsedState:
1413 eventStart: Time
1414 eventEnd: int
1415 previousState: TMAState
1416 state: TMAState
1418 @staticmethod
1419 def isToday(dayObs: int) -> bool:
1420 """Find out if the specified dayObs is today, or in the past.
1422 If the day is today, the function returns ``True``, if it is in the
1423 past it returns ``False``. If the day is in the future, a
1424 ``ValueError`` is raised, as this indicates there is likely an
1425 off-by-one type error somewhere in the logic.
1427 Parameters
1428 ----------
1429 dayObs : `int`
1430 The dayObs to check, in the format YYYYMMDD.
1432 Returns
1433 -------
1434 isToday : `bool`
1435 ``True`` if the dayObs is today, ``False`` if it is in the past.
1437 Raises
1438 ValueError: if the dayObs is in the future.
1439 """
1440 todayDayObs = getCurrentDayObs_int()
1441 if dayObs == todayDayObs:
1442 return True
1443 if dayObs > todayDayObs:
1444 raise ValueError("dayObs is in the future")
1445 return False
1447 @staticmethod
1448 def _shortName(topic: str) -> str:
1449 """Get the short name of a topic.
1451 Parameters
1452 ----------
1453 topic : `str`
1454 The topic to get the short name of.
1456 Returns
1457 -------
1458 shortName : `str`
1459 The short name of the topic, e.g. 'azimuthInPosition'
1460 """
1461 # get, for example 'azimuthInPosition' from
1462 # lsst.sal.MTMount.logevent_azimuthInPosition
1463 return topic.split("_")[-1]
1465 def _mergeData(self, data: dict[str, pd.DataFrame]) -> pd.DataFrame:
1466 """Merge a dict of dataframes based on private_efdStamp, recording
1467 where each row came from.
1469 Given a dict or dataframes, keyed by topic, merge them into a single
1470 dataframe, adding a column to record which topic each row came from.
1472 Parameters
1473 ----------
1474 data : `dict` of `str` : `pd.DataFrame`
1475 The dataframes to merge.
1477 Returns
1478 -------
1479 merged : `pd.DataFrame`
1480 The merged dataframe.
1481 """
1482 excludeColumns = ["private_efdStamp", "rowFor"]
1484 mergeArgs = {
1485 "how": "outer",
1486 "sort": True,
1487 }
1489 merged = None
1490 originalRowCounter = 0
1492 # Iterate over the keys and merge the corresponding DataFrames
1493 for key, df in data.items():
1494 if df.empty:
1495 # Must skip the df if it's empty, otherwise the merge will fail
1496 # due to lack of private_efdStamp. Because other axes might
1497 # still be in motion, so we still want to merge what we have
1498 continue
1500 originalRowCounter += len(df)
1501 component = self._shortName(key) # Add suffix to column names to identify the source
1502 suffix = "_" + component
1504 df["rowFor"] = component
1506 columnsToSuffix = [col for col in df.columns if col not in excludeColumns]
1507 df_to_suffix = df[columnsToSuffix].add_suffix(suffix)
1508 df = pd.concat([df[excludeColumns], df_to_suffix], axis=1)
1510 if merged is None:
1511 merged = df.copy()
1512 else:
1513 merged = pd.merge(merged, df, **mergeArgs)
1514 assert merged is not None, "merged data is None"
1515 merged = merged.loc[:, ~merged.columns.duplicated()] # Remove duplicate columns after merge
1517 if len(merged) != originalRowCounter:
1518 self.log.warning(
1519 "Merged data has a different number of rows to the original data, some"
1520 " timestamps (rows) will contain more than one piece of actual information."
1521 )
1523 # if the index is still a DatetimeIndex here then we didn't actually
1524 # merge any data, so there is only data from a single component.
1525 # This is likely to result in no events, but not necessarily, and for
1526 # generality, instead we convert to a range index to ensure consistency
1527 # in the returned data, and allow processing to continue.
1528 if isinstance(merged.index, pd.DatetimeIndex):
1529 self.log.warning("Data was only found for a single component in the EFD.")
1530 merged.reset_index(drop=True, inplace=True)
1532 return merged
1534 def getEvent(self, dayObs: int, seqNum: int) -> TMAEvent | None:
1535 """Get a specific event for a given dayObs and seqNum.
1537 Repeated calls for the same ``dayObs`` will use the cached data if the
1538 day is in the past, and so will be much quicker. If the ``dayObs`` is
1539 the current day then the EFD will be queried for new data for each
1540 call, so a call which returns ``None`` on the first try might return an
1541 event on the next, if the TMA is still moving and thus generating
1542 events.
1544 Parameters
1545 ----------
1546 dayObs : `int`
1547 The dayObs to get the event for.
1548 seqNum : `int`
1549 The sequence number of the event to get.
1551 Returns
1552 -------
1553 event : `lsst.summit.utils.tmaUtils.TMAEvent`
1554 The event for the specified dayObs and seqNum, or `None` if the
1555 event was not found.
1556 """
1557 events = self.getEvents(dayObs)
1558 if seqNum <= len(events):
1559 event = events[seqNum]
1560 if event.seqNum != seqNum:
1561 # it's zero-indexed and contiguous so this must be true but
1562 # a sanity check doesn't hurt.
1563 raise AssertionError(f"Event sequence number mismatch: {event.seqNum} != {seqNum}")
1564 return event
1565 else:
1566 self.log.warning(f"Event {seqNum} not found for {dayObs}")
1567 return None
1569 def getEvents(self, dayObs: int, addBlockInfo: bool = True) -> list[TMAEvent]:
1570 """Get the TMA events for the specified dayObs.
1572 Gets the required mount data from the cache or the EFD as required,
1573 handling whether we're working with live vs historical data. The
1574 dataframes from the EFD is merged and applied to the TMAStateMachine,
1575 and that series of state changes is used to generate a list of
1576 TmaEvents for the day's data.
1578 If the data is for the current day, i.e. if new events can potentially
1579 land, then if the last event is "open" (meaning that the TMA appears to
1580 be in motion and thus the event is growing with time), then that event
1581 is excluded from the event list as it is expected to be changing with
1582 time, and will likely close eventually. However, if that situation
1583 occurs on a day in the past, then that event can never close, and the
1584 event is therefore included, but a warning about the open event is
1585 logged.
1587 Parameters
1588 ----------
1589 dayObs : `int`
1590 The dayObs for which to get the events.
1591 addBlockInfo : `bool`, optional
1592 Whether to add block information to the events. This allows
1593 skipping this step for speed when generating events for purposes
1594 which don't need block information.
1596 Returns
1597 -------
1598 events : `list` of `lsst.summit.utils.tmaUtils.TMAState`
1599 The events for the specified dayObs.
1600 """
1601 workingLive = self.isToday(dayObs)
1602 data = None
1604 if workingLive:
1605 # it's potentially updating data, so we must update the date
1606 # regarless of whether we have it already or not
1607 self.log.info(f"Updating mount data for {dayObs} from the EFD")
1608 self._getEfdDataForDayObs(dayObs)
1609 data = self._data[dayObs]
1610 elif dayObs in self._data:
1611 # data is in the cache and it's not being updated, so use it
1612 data = self._data[dayObs]
1613 elif dayObs not in self._data:
1614 # we don't have the data yet, but it's not growing, so put it in
1615 # the cache and use it from there
1616 self.log.info(f"Retrieving mount data for {dayObs} from the EFD")
1617 self._getEfdDataForDayObs(dayObs)
1618 data = self._data[dayObs]
1619 else:
1620 raise RuntimeError("This should never happen")
1622 # if we don't have something to work with, log a warning and return
1623 if not self.dataFound(data):
1624 self.log.warning(f"No EFD data found for {dayObs=}")
1625 return []
1627 # applies the data to the state machine, and generates events from the
1628 # series of states which results
1629 events = self._calculateEventsFromMergedData(
1630 data, dayObs, dataIsForCurrentDay=workingLive, addBlockInfo=addBlockInfo
1631 )
1632 if not events:
1633 self.log.warning(f"Failed to calculate any events for {dayObs=} despite EFD data existing!")
1634 return events
1636 @staticmethod
1637 def dataFound(data: pd.DataFrame) -> bool:
1638 """Check if any data was found.
1640 Parameters
1641 ----------
1642 data : `pd.DataFrame`
1643 The merged dataframe to check.
1645 Returns
1646 -------
1647 dataFound : `bool`
1648 Whether data was found.
1649 """
1650 # You can't just compare to with data == NO_DATA_SENTINEL because
1651 # `data` is usually a dataframe, and you can't compare a dataframe to a
1652 # string directly.
1653 return not (isinstance(data, str) and data == NO_DATA_SENTINEL)
1655 def _getEfdDataForDayObs(self, dayObs: int) -> None:
1656 """Get the EFD data for the specified dayObs and store it in the cache.
1658 Gets the EFD data for all components, as a dict of dataframes keyed by
1659 component name. These are then merged into a single dataframe in time
1660 order, based on each row's `private_efdStamp`. This is then stored in
1661 self._data[dayObs].
1663 If no data is found, the value is set to ``NO_DATA_SENTINEL`` to
1664 differentiate this from ``None``, as this is what you'd get if you
1665 queried the cache with `self._data.get(dayObs)`. It also marks that we
1666 have already queried this day.
1668 Parameters
1669 ----------
1670 dayObs : `int`
1671 The dayObs to query.
1672 """
1673 data = {}
1674 for component in itertools.chain(
1675 self._movingComponents, self._inPositionComponents, self._stateComponents
1676 ):
1677 data[component] = getEfdData(
1678 self.client,
1679 component,
1680 dayObs=dayObs,
1681 warn=False,
1682 raiseIfTopicNotInSchema=False,
1683 )
1684 self.log.debug(f"Found {len(data[component])} for {component}")
1686 if all(dataframe.empty for dataframe in data.values()):
1687 # if every single dataframe is empty, set the sentinel and don't
1688 # try to merge anything, otherwise merge all the data we found
1689 self.log.debug(f"No data found for {dayObs=}")
1690 # a sentinel value that's not None
1691 self._data[dayObs] = NO_DATA_SENTINEL
1692 else:
1693 merged = self._mergeData(data)
1694 self._data[dayObs] = merged
1696 def _calculateEventsFromMergedData(
1697 self, data: pd.DataFrame, dayObs: int, dataIsForCurrentDay: bool, addBlockInfo: bool
1698 ) -> list[TMAEvent]:
1699 """Calculate the list of events from the merged data.
1701 Runs the merged data, row by row, through the TMA state machine (with
1702 ``tma.apply``) to get the overall TMA state at each row, building a
1703 dict of these states, keyed by row number.
1705 This time-series of TMA states are then looped over (in
1706 `_statesToEventTuples`), building a list of tuples representing the
1707 start and end of each event, the type of the event, and the reason for
1708 the event ending.
1710 This list of tuples is then passed to ``_makeEventsFromStateTuples``,
1711 which actually creates the ``TMAEvent`` objects.
1713 Parameters
1714 ----------
1715 data : `pd.DataFrame`
1716 The merged dataframe to use.
1717 dayObs : `int`
1718 The dayObs for the data.
1719 dataIsForCurrentDay : `bool`
1720 Whether the data is for the current day. Determines whether to
1721 allow an open last event or not.
1722 addBlockInfo : `bool`
1723 Whether to add block information to the events. This allows
1724 skipping this step for speed when generating events for purposes
1725 which don't need block information.
1727 Returns
1728 -------
1729 events : `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
1730 The events for the specified dayObs.
1731 """
1732 engineeringMode = True
1733 tma = TMAStateMachine(engineeringMode=engineeringMode)
1735 # For now, we assume that the TMA starts each day able to move, but
1736 # stationary. If this turns out to cause problems, we will need to
1737 # change to loading data from the previous day(s), and looking back
1738 # through it in time until a state change has been found for every
1739 # axis. For now though, Bruno et. al think this is acceptable and
1740 # preferable.
1741 _initializeTma(tma)
1743 tmaStates = {}
1744 for rowNum, row in data.iterrows():
1745 assert isinstance(rowNum, int)
1746 tma.apply(row)
1747 tmaStates[rowNum] = tma.state
1749 stateTuples = self._statesToEventTuples(tmaStates, dataIsForCurrentDay)
1750 events = self._makeEventsFromStateTuples(stateTuples, dayObs, data)
1751 if addBlockInfo:
1752 self.addBlockDataToEvents(dayObs, events)
1753 return events
1755 def _statesToEventTuples(
1756 self, states: dict[int, TMAState], dataIsForCurrentDay: bool
1757 ) -> list[ParsedState]:
1758 """Get the event-tuples from the dictionary of TMAStates.
1760 Chunks the states into blocks of the same state, so that we can create
1761 an event for each block in `_makeEventsFromStateTuples`. Off-type
1762 states are skipped over, with each event starting when the telescope
1763 next resumes motion or changes to a different type of motion state,
1764 i.e. from non-tracking type movement (MOVE_POINT_TO_POINT, JOGGING,
1765 TRACKING-but-not-in-position, i.e. slewing) to a tracking type
1766 movement, or vice versa.
1768 Parameters
1769 ----------
1770 states : `dict` of `int` : `lsst.summit.utils.tmaUtils.TMAState`
1771 The states of the TMA, keyed by row number.
1772 dataIsForCurrentDay : `bool`
1773 Whether the data is for the current day. Determines whether to
1774 allow and open last event or not.
1776 Returns
1777 -------
1778 parsedStates : `list` of `tuple`
1779 The parsed states, as a list of tuples of the form:
1780 ``(eventStart, eventEnd, eventType, endReason)``
1781 """
1782 # Consider rewriting this with states as a list and using pop(0)?
1783 skipStates = (TMAState.STOPPED, TMAState.OFF, TMAState.FAULT)
1785 parsedStates = []
1786 eventStart = None
1787 rowNum = 0
1788 nRows = len(states)
1789 while rowNum < nRows:
1790 previousState = None
1791 state = states[rowNum]
1792 # if we're not in an event, fast forward through off-like rows
1793 # until a new event starts
1794 if eventStart is None and state in skipStates:
1795 rowNum += 1
1796 continue
1798 # we've started a new event, so walk through it and find the end
1799 eventStart = rowNum
1800 previousState = state
1801 rowNum += 1 # move to the next row before starting the while loop
1802 if rowNum == nRows:
1803 # we've reached the end of the data, and we're still in an
1804 # event, so don't return this presumably in-progress event
1805 self.log.warning("Reached the end of the data while starting a new event")
1806 break
1807 state = states[rowNum]
1808 while state == previousState:
1809 rowNum += 1
1810 if rowNum == nRows:
1811 break
1812 state = states[rowNum]
1813 parsedStates.append(
1814 self.ParsedState(
1815 eventStart=eventStart, eventEnd=rowNum, previousState=previousState, state=state
1816 )
1817 )
1818 if state in skipStates:
1819 eventStart = None
1821 # done parsing, just check the last event is valid
1822 if parsedStates: # ensure we have at least one event
1823 lastEvent = parsedStates[-1]
1824 if lastEvent.eventEnd == nRows:
1825 # Generally, you *want* the timespan for an event to be the
1826 # first row of the next event, because you were in that state
1827 # right up until that state change. However, if that event is
1828 # a) the last one of the day and b) runs right up until the end
1829 # of the dataframe, then there isn't another row, so this will
1830 # overrun the array.
1831 #
1832 # If the data is for the current day then this isn't a worry,
1833 # as we're likely still taking data, and this event will likely
1834 # close yet, so we don't issue a warning, and simply drop the
1835 # event from the list.
1837 # However, if the data is for a past day then no new data will
1838 # come to close the event, so allow the event to be "open", and
1839 # issue a warning
1840 if dataIsForCurrentDay:
1841 self.log.info("Discarding open (likely in-progess) final event from current day's events")
1842 parsedStates = parsedStates[:-1]
1843 else:
1844 self.log.warning("Last event ends open, forcing it to end at end of the day's data")
1845 # it's a tuple, so (deliberately) awkward to modify
1846 parsedStates[-1] = self.ParsedState(
1847 eventStart=lastEvent.eventStart,
1848 eventEnd=lastEvent.eventEnd - 1,
1849 previousState=lastEvent.previousState,
1850 state=lastEvent.state,
1851 )
1853 return parsedStates
1855 def addBlockDataToEvents(
1856 self,
1857 dayObs: int,
1858 events: list[TMAEvent] | TMAEvent,
1859 ) -> None:
1860 """Find all the block data in the EFD for the specified events.
1862 Finds all the block data in the EFD relating to the events, parses it,
1863 from the rows of the dataframe, and adds it to the events in place.
1865 Parameters
1866 ----------
1867 events : `lsst.summit.utils.tmaUtils.TMAEvent` or
1868 `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
1869 One or more events to get the block data for.
1870 """
1871 try:
1872 blockParser = BlockParser(dayObs, client=self.client)
1873 except Exception as e:
1874 # adding the block data should never cause a failure so if we can't
1875 # get the block data, log a warning and return. It is, however,
1876 # never expected, so use log.exception to get the full traceback
1877 # and scare users so it gets reported
1878 self.log.exception(f"Failed to parse block data for {dayObs=}, {e}")
1879 return
1880 if not isinstance(events, list):
1881 events = [events]
1882 blocks = blockParser.getBlockNums()
1883 blockDict = {}
1884 for block in blocks:
1885 blockDict[block] = blockParser.getSeqNums(block)
1887 for block, seqNums in blockDict.items():
1888 for seqNum in seqNums:
1889 blockInfo = blockParser.getBlockInfo(block=block, seqNum=seqNum)
1891 relatedEvents = blockParser.getEventsForBlock(events, block=block, seqNum=seqNum)
1892 for event in relatedEvents:
1893 toSet = [blockInfo]
1894 if event.blockInfos is not None:
1895 existingInfo = event.blockInfos
1896 existingInfo.append(blockInfo)
1897 toSet = existingInfo
1899 # Add the blockInfo to the TMAEvent. Because this is a
1900 # frozen dataclass, use object.__setattr__ to set the
1901 # attribute. This is the correct way to set a frozen
1902 # dataclass attribute after creation.
1903 object.__setattr__(event, "blockInfos", toSet)
1905 def _makeEventsFromStateTuples(
1906 self, states: list[ParsedState], dayObs: int, data: pd.DataFrame
1907 ) -> list[TMAEvent]:
1908 """For the list of state-tuples, create a list of ``TMAEvent`` objects.
1910 Given the underlying data, and the start/stop points for each event,
1911 create the TMAEvent objects for the dayObs.
1913 Parameters
1914 ----------
1915 states : `list` of `tuple`
1916 The parsed states, as a list of tuples of the form:
1917 ``(eventStart, eventEnd, eventType, endReason)``
1918 dayObs : `int`
1919 The dayObs for the data.
1920 data : `pd.DataFrame`
1921 The merged dataframe.
1923 Returns
1924 -------
1925 events : `list` of `lsst.summit.utils.tmaUtils.TMAEvent`
1926 The events for the specified dayObs.
1927 """
1928 seqNum = 0
1929 events = []
1930 for parsedState in states:
1931 begin = data.iloc[parsedState.eventStart]["private_efdStamp"]
1932 end = data.iloc[parsedState.eventEnd]["private_efdStamp"]
1933 beginAstropy = efdTimestampToAstropy(begin)
1934 endAstropy = efdTimestampToAstropy(end)
1935 duration = end - begin
1936 event = TMAEvent(
1937 dayObs=dayObs,
1938 seqNum=seqNum,
1939 type=parsedState.previousState,
1940 endReason=parsedState.state,
1941 duration=duration,
1942 begin=beginAstropy,
1943 end=endAstropy,
1944 blockInfos=[], # this is added later
1945 _startRow=parsedState.eventStart,
1946 _endRow=parsedState.eventEnd,
1947 )
1948 events.append(event)
1949 seqNum += 1
1950 return events
1952 @staticmethod
1953 def printTmaDetailedState(tma: TMAStateMachine) -> None:
1954 """Print the full state of all the components of the TMA.
1956 Currently this is the azimuth and elevation axes' power and motion
1957 states, and their respective inPosition statuses.
1959 Parameters
1960 ----------
1961 tma : `lsst.summit.utils.tmaUtils.TMAStateMachine`
1962 The TMA state machine in the state we want to print.
1963 """
1964 axes = ["azimuth", "elevation"]
1965 p = tma._parts
1966 axisPad = len(max(axes, key=len)) # length of the longest axis string == 9 here, but this is general
1967 motionPad = max(len(s.name) for s in AxisMotionState)
1968 powerPad = max(len(s.name) for s in PowerState)
1970 # example output to show what's being done with the padding:
1971 # azimuth - Power: ON Motion: STOPPED InPosition: True # noqa: W505
1972 # elevation - Power: ON Motion: MOVING_POINT_TO_POINT InPosition: False # noqa: W505
1973 for axis in axes:
1974 powerState = p[f"{axis}SystemState"]
1975 motionState = p[f"{axis}MotionState"]
1976 assert isinstance(powerState, PowerState)
1977 assert isinstance(motionState, AxisMotionState)
1978 print(
1979 f"{axis:>{axisPad}} - "
1980 f"Power: {powerState.name:>{powerPad}} "
1981 f"Motion: {motionState.name:>{motionPad}} "
1982 f"InPosition: {p[f'{axis}InPosition']}"
1983 )
1984 print(f"Overall system state: {tma.state.name}")
1986 def printFullDayStateEvolution(self, dayObs: int, taiOrUtc: str = "utc") -> None:
1987 """Print the full TMA state evolution for the specified dayObs.
1989 Replays all the data from the EFD for the specified dayObs through
1990 the TMA state machine, and prints both the overall and detailed state
1991 of the TMA for each row.
1993 Parameters
1994 ----------
1995 dayObs : `int`
1996 The dayObs for which to print the state evolution.
1997 taiOrUtc : `str`, optional
1998 Whether to print the timestamps in TAI or UTC. Default is UTC.
1999 """
2000 # create a fake event which spans the whole day, and then use
2001 # printEventDetails code while skipping the header to print the
2002 # evolution.
2003 _ = self.getEvents(dayObs) # ensure the data has been retrieved from the EFD
2004 data = self._data[dayObs]
2005 lastRowNum = len(data) - 1
2007 fakeEvent = TMAEvent(
2008 dayObs=dayObs,
2009 seqNum=-1, # anything will do
2010 type=TMAState.OFF, # anything will do
2011 endReason=TMAState.OFF, # anything will do
2012 duration=-1, # anything will do
2013 begin=efdTimestampToAstropy(data.iloc[0]["private_efdStamp"]),
2014 end=efdTimestampToAstropy(data.iloc[-1]["private_efdStamp"]),
2015 _startRow=0,
2016 _endRow=lastRowNum,
2017 )
2018 self.printEventDetails(fakeEvent, taiOrUtc=taiOrUtc, printHeader=False)
2020 def printEventDetails(
2021 self,
2022 event: TMAEvent,
2023 taiOrUtc: str = "tai",
2024 printHeader: bool = False,
2025 ) -> None:
2026 """Print a detailed breakdown of all state transitions during an event.
2028 Note: this is not the most efficient way to do this, but it is much the
2029 cleanest with respect to the actual state machine application and event
2030 generation code, and is easily fast enough for the cases it will be
2031 used for. It is not worth complicating the normal state machine logic
2032 to try to use this code.
2034 Parameters
2035 ----------
2036 event : `lsst.summit.utils.tmaUtils.TMAEvent`
2037 The event to display the details of.
2038 taiOrUtc : `str`, optional
2039 Whether to display time strings in TAI or UTC. Defaults to TAI.
2040 Case insensitive.
2041 printHeader : `bool`, optional
2042 Whether to print the event summary. Defaults to True. The primary
2043 reason for the existence of this option is so that this same
2044 printing function can be used to show the evolution of a whole day
2045 by supplying a fake event which spans the whole day, but this event
2046 necessarily has a meaningless summary, and so needs suppressing.
2047 """
2048 taiOrUtc = taiOrUtc.lower()
2049 if taiOrUtc not in ["tai", "utc"]:
2050 raise ValueError(f"Got unsuppoted value for {taiOrUtc=}")
2051 useUtc = taiOrUtc == "utc"
2053 if printHeader:
2054 print(
2055 f"Details for {event.duration:.2f}s {event.type.name} event dayObs={event.dayObs}"
2056 f" seqNum={event.seqNum}:"
2057 )
2058 print(f"- Event began at: {event.begin.utc.isot if useUtc else event.begin.isot}")
2059 print(f"- Event ended at: {event.end.utc.isot if useUtc else event.end.isot}")
2061 dayObs = event.dayObs
2062 data = self._data[dayObs]
2063 startRow = event._startRow
2064 endRow = event._endRow
2065 nRowsToApply = endRow - startRow + 1
2066 print(f"\nTotal number of rows in the merged dataframe: {len(data)}")
2067 if printHeader:
2068 print(f"of which rows {startRow} to {endRow} (inclusive) relate to this event.")
2070 # reconstruct all the states
2071 tma = TMAStateMachine(engineeringMode=True)
2072 _initializeTma(tma)
2074 tmaStates = {}
2075 firstAppliedRow = True # flag to print a header on the first row that's applied
2076 for rowNum, row in data.iterrows(): # must replay rows right from start to get full correct state
2077 if rowNum == startRow:
2078 # we've not yet applied this row, so this is the state just
2079 # before event
2080 print(f"\nBefore the event the TMA was in state {tma.state.name}:")
2081 self.printTmaDetailedState(tma)
2083 if rowNum >= startRow and rowNum <= endRow:
2084 if firstAppliedRow: # only print this intro on the first row we're applying
2085 print(
2086 f"\nThen, applying the {nRowsToApply} rows of data for this event, the state"
2087 " evolved as follows:\n"
2088 )
2089 firstAppliedRow = False
2091 # break the row down and print its details
2092 rowFor = row["rowFor"]
2093 axis, rowType = getAxisAndType(rowFor) # e.g. elevation, MotionState
2094 value = tma._getRowPayload(row, rowType, rowFor)
2095 valueStr = f"{str(value) if isinstance(value, bool) else value.name}"
2096 rowTime = efdTimestampToAstropy(row["private_efdStamp"])
2097 print(
2098 f"On row {rowNum} the {axis} axis had the {rowType} set to {valueStr} at"
2099 f" {rowTime.utc.isot if useUtc else rowTime.isot}"
2100 )
2102 # then apply it as usual, printing the state right afterwards
2103 tma.apply(row)
2104 tmaStates[rowNum] = tma.state
2105 self.printTmaDetailedState(tma)
2106 print()
2108 else:
2109 # if it's not in the range of interest then just apply it
2110 # silently as usual
2111 tma.apply(row)
2112 tmaStates[rowNum] = tma.state
2114 def findEvent(self, time: astropy.time.Time) -> TMAEvent | None:
2115 """Find the event which contains the specified time.
2117 If the specified time lies within an event, that event is returned. If
2118 it is at the exact start, that is logged, and if that start point is
2119 shared by the end of the previous event, that is logged too. If the
2120 event lies between events, the events either side are logged, but
2121 ``None`` is returned. If the time lies before the first event of the
2122 day a warning is logged, as for times after the last event of the day.
2124 Parameters
2125 ----------
2126 time : `astropy.time.Time`
2127 The time.
2129 Returns
2130 -------
2131 event : `lsst.summit.utils.tmaUtils.TMAEvent` or `None`
2132 The event which contains the specified time, or ``None`` if the
2133 time doesn't fall during an event.
2134 """
2135 # there are five possible cases:
2136 # 1) the time lies before the first event of the day
2137 # 2) the time lies after the last event of the day
2138 # 3) the time lies within an event
2139 # 3a) the time is exactly at the start of an event
2140 # 3b) if so, time can be shared by the end of the previous event if
2141 # they are contiguous
2142 # 4) the time lies between two events
2143 # 5) the time is exactly at end of the last event of the day. This is
2144 # an issue because event end times are exclusive, so this time is
2145 # not technically in that event, it's the moment it closes (and if
2146 # there *was* an event which followed contiguously, it would be in
2147 # that event instead, which is what motivates this definition of
2148 # lies within what event)
2150 dayObs = getDayObsForTime(time)
2151 # we know this is on the right day, and definitely before the specified
2152 # time, but sanity check this before continuing as this needs to be
2153 # true for this to give the correct answer
2154 assert getDayObsStartTime(dayObs) <= time
2155 assert getDayObsEndTime(dayObs) > time
2157 # command start to many log messages so define once here
2158 logStart = f"Specified time {time.isot} falls on {dayObs=}"
2160 events = self.getEvents(dayObs, addBlockInfo=False)
2161 if len(events) == 0:
2162 self.log.warning(f"There are no events found for {dayObs}")
2163 return None
2165 # check case 1)
2166 if time < events[0].begin:
2167 self.log.warning(f"{logStart} and is before the first event of the day")
2168 return None
2170 # check case 2)
2171 if time > events[-1].end:
2172 self.log.warning(f"{logStart} and is after the last event of the day")
2173 return None
2175 # check case 5)
2176 if time == events[-1].end:
2177 self.log.warning(
2178 f"{logStart} and is exactly at the end of the last event of the day"
2179 f" (seqnum={events[-1].seqNum}). Because event intervals are half-open, this"
2180 " time does not technically lie in any event"
2181 )
2182 return None
2184 # we are now either in an event, or between events. Walk through the
2185 # events, and if the end of the event is after the specified time, then
2186 # we're either in it or past it, so check if we're in.
2187 for eventNum, event in enumerate(events):
2188 if event.end > time: # case 3) we are now into or past the right event
2189 # the event end encloses the time, so note the > and not >=,
2190 # this must be strictly greater, we check the overlap case
2191 # later
2192 if time >= event.begin: # we're fully inside the event, so return it.
2193 # 3a) before returning, check if we're exactly at the start
2194 # of the event, and if so, log it. Then 3b) also check if
2195 # we're at the exact end of the previous event, and if so,
2196 # log that too.
2197 if time == event.begin:
2198 self.log.info(f"{logStart} and is exactly at the start of event" f" {eventNum}")
2199 if eventNum == 0: # I think this is actually impossible, but check anyway
2200 return event # can't check the previous event so return here
2201 previousEvent = events[eventNum - 1]
2202 if previousEvent.end == time:
2203 self.log.info(
2204 "Previous event is contiguous, so this time is also at the exact"
2205 f" end of {eventNum - 1}"
2206 )
2207 return event
2208 else: # case 4)
2209 # the event end is past the time, but it's not inside the
2210 # event, so we're between events. Log which we're between
2211 # and return None
2212 previousEvent = events[eventNum - 1]
2213 timeAfterPrev = (time - previousEvent.end).to_datetime()
2214 naturalTimeAfterPrev = humanize.naturaldelta(timeAfterPrev, minimum_unit="MICROSECONDS")
2215 timeBeforeCurrent = (event.begin - time).to_datetime()
2216 naturalTimeBeforeCurrent = humanize.naturaldelta(
2217 timeBeforeCurrent, minimum_unit="MICROSECONDS"
2218 )
2219 self.log.info(
2220 f"{logStart} and lies"
2221 f" {naturalTimeAfterPrev} after the end of event {previousEvent.seqNum}"
2222 f" and {naturalTimeBeforeCurrent} before the start of event {event.seqNum}."
2223 )
2224 return None
2226 raise RuntimeError(
2227 "Event finding logic fundamentally failed, which should never happen - the code" " needs fixing"
2228 )