From 93caea9ff842ba6fcaa706533df787340455bd76 Mon Sep 17 00:00:00 2001 From: Bruno Quint Date: Thu, 14 Sep 2023 06:09:00 -0700 Subject: [PATCH] Implement intertia_compensation_system.py analysis --- .../m1m3/inertia_compensation_system.py | 275 ++++++++++++++++-- 1 file changed, 250 insertions(+), 25 deletions(-) diff --git a/python/lsst/summit/testing/analysis/m1m3/inertia_compensation_system.py b/python/lsst/summit/testing/analysis/m1m3/inertia_compensation_system.py index d8b8f9e..2ae54cc 100644 --- a/python/lsst/summit/testing/analysis/m1m3/inertia_compensation_system.py +++ b/python/lsst/summit/testing/analysis/m1m3/inertia_compensation_system.py @@ -1,10 +1,16 @@ #!/usr/bin/python import logging +import numpy as np import pandas as pd +from astropy import units as u +from astropy.time import Time +from datetime import timedelta + +from lsst.summit.utils.efdUtils import getEfdData, makeEfdClient from lsst.summit.utils.tmaUtils import TMAEventMaker, TMAEvent -# TODO @bquint: fix logger - not working now. + formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s") formatter.datefmt = "%Y-%m-%d %H:%M:%S" @@ -30,12 +36,163 @@ class InertiaCompensationSystemAnalysis: Abtract representation of a slew event. """ - # ToDo @bquint - def __init__(event: TMAEvent): - raise NotImplementedError + def __init__( + self, + event: TMAEvent, + inner_pad: float = 0.0, + outer_pad: float = 0.0, + n_sigma: float = 1.0, + ): + self.event = event + self.inner_pad = inner_pad * u.second + self.outer_pad = outer_pad * u.second + self.n_sigma = n_sigma + self.client = makeEfdClient() + + # TODO: Find a better way to implement the logger inside the class + self.logger = logger + + self.number_of_hardpoints = 6 + self.measured_forces_topics = [ + f"measuredForce{i}" for i in range(self.number_of_hardpoints) + ] + + self.df = None + self.stats = None + + async def run(self): + """Performs all the measurements""" + self.df = self.query_dataset() + self.stats = self.get_stats() + + def query_dataset(self): + """ + Queries all the relevant data, resample them to have the same requency and merge them in a single dataframe. + + Returns + ------- + pd.DataFrame + """ + self.logger.info("Querying dataset: m1m3 hp meadures forces") + hp_measured_forces = getEfdData( + self.client, + "lsst.sal.MTM1M3.hardpointActuatorData", + columns=self.measured_forces_topics, + event=self.event, + prePadding=self.outer_pad, + postPadding=self.outer_pad, + warn=False, + ) + + self.logger.info("Querying dataset: mtmount azimuth torque and velocity") + tma_az = getEfdData( + self.client, + "lsst.sal.MTMount.azimuth", + columns=["timestamp", "actualTorque", "actualVelocity"], + event=self.event, + prePadding=self.outer_pad, + postPadding=self.outer_pad, + warn=False, + ) + + tma_az = tma_az.rename( + columns={ + "actualTorque": "az_actual_torque", + "actualVelocity": "az_actual_velocity", + } + ) + tma_az["timestamp"] = Time( + tma_az["timestamp"], format="unix_tai", scale="utc" + ).datetime + tma_az.set_index("timestamp", inplace=True) + tma_az.index = tma_az.index.tz_localize("UTC") + + self.logger.info("Querying dataset: mtmount elevation torque and velocity") + tma_el = getEfdData( + self.client, + "lsst.sal.MTMount.elevation", + columns=["timestamp", "actualTorque", "actualVelocity"], + event=self.event, + prePadding=self.outer_pad, + postPadding=self.outer_pad, + warn=False, + ) + + tma_el = tma_el.rename( + columns={ + "actualTorque": "el_actual_torque", + "actualVelocity": "el_actual_velocity", + } + ) + tma_el["timestamp"] = Time( + tma_el["timestamp"], format="unix_tai", scale="utc" + ).datetime + tma_el.set_index("timestamp", inplace=True) + tma_el.index = tma_el.index.tz_localize("UTC") + + merge_cfg = { + "left_index": True, + "right_index": True, + "tolerance": timedelta(seconds=1), + "direction": "nearest", + } + + merged_df = pd.merge_asof(hp_measured_forces, tma_az, **merge_cfg) + merged_df = pd.merge_asof(merged_df, tma_el, **merge_cfg) + merged_df[["az_actual_torque", "el_actual_torque"]] = ( + 1e-3 * merged_df[["az_actual_torque", "el_actual_torque"]] + ) + + return merged_df + + def get_stats(self): + """ + Calculate statistics for each column in a given dataset. + + Returns + ------- + pandas.DataFrame + A DataFrame containing calculated statistics for each column in the dataset. + For each column, the statistics include minimum, maximum, and peak-to-peak values. + + Notes + ----- + This function computes statistics for each column in the provided dataset. It utilizes the `get_minmax` function + to calculate minimum, maximum, and peak-to-peak values for each column's data. + """ + cols = self.measured_forces_topics + full_slew_stats = pd.DataFrame( + data=[self.get_slew_minmax(self.df[col]) for col in cols], index=cols + ) + self.logger.info("Find stable time window") + begin, end = self.find_stable_region() + + self.logger.debug("Update begin and end times") + begin = begin + self.inner_pad + end = end - self.inner_pad + + self.logger.debug("Calculate statistics in stable time window") + stable_slew_stats = pd.DataFrame( + data=[ + self.get_stats_in_torqueless_interval( + self.df[col].loc[begin.isot : end.isot] + ) + for col in cols + ], + index=cols, + ) + + self.logger.debug("Concatenate statistics") + stats = pd.concat((full_slew_stats, stable_slew_stats), axis=1) + + return stats + + def get_midppoint(self): + """Returns the halfway point between begin and end""" + return self.df.index[len(self.df.index) // 2] @staticmethod - def get_minmax_from_series(s): + def get_slew_minmax(s): """ Calculate minimum, maximum, and peak-to-peak values for a data-series. @@ -47,8 +204,7 @@ def get_minmax_from_series(s): Returns ------- pandas.Series - A Series containing the following calculated values for the two - halves of the input Series: + A Series containing the following calculated values for the two halves of the input Series: - min: Minimum value of the Series. - max: Maximum value of the Series. - ptp: Peak-to-peak (ptp) value of the Series (abs(max - min)). @@ -61,23 +217,26 @@ def get_minmax_from_series(s): return result @staticmethod - def get_avgmedstd_from_series(s): + def get_stats_in_torqueless_interval(s): """ - Calculate average, median, and standard deviation for a data-series. + Calculate statistical measures within a torqueless interval. - Parameters - ---------- + This static method computes descriptive statistics for a given pandas Series within + a torqueless interval. The torqueless interval represents a period of the data + analysis when no external torque is applied. + + Parameters: + ----------- s : pandas.Series - The input pandas Series containing data. + A pandas Series containing data values for analysis. - Returns - ------- + Returns: + -------- pandas.Series - A Series containing the following calculated values for the two - halves of the input Series: - - mean: The mean value of the Series. - - median: The median value of the Series. - - std: The standard deviation of the Series. + A pandas Series containing the following statistical measures: + - Mean: The arithmetic mean of the data. + - Median: The median value of the data. + - Standard Deviation (Std): The standard deviation of the data. """ result = pd.Series( data=[s.mean(), s.median(), s.std()], @@ -86,6 +245,56 @@ def get_avgmedstd_from_series(s): ) return result + def find_stable_region(self): + """ + ToDo @b1quint: add docstring + """ + az_torque = self.df["az_actual_torque"] + az_torque_regions = find_adjacent_true_regions( + np.abs(az_torque - az_torque.mean()) < self.n_sigma * az_torque.std() + ) + + el_torque = self.df["el_actual_torque"] + el_torque_regions = find_adjacent_true_regions( + np.abs(el_torque - el_torque.mean()) < self.n_sigma * el_torque.std() + ) + + stable_begin = max([reg[0] for reg in az_torque_regions + el_torque_regions]) + stable_begin = Time(stable_begin, scale="utc") + + stable_end = min([reg[-1] for reg in az_torque_regions + el_torque_regions]) + stable_end = Time(stable_end, scale="utc") + + return stable_begin, stable_end + + +def find_adjacent_true_regions(series, min_adjacent=None): + """ + Find regions in a boolean Series containing adjacent True values. + + Parameters + ---------- + series : pd.Series + The boolean Series to search for regions. + + min_adjacent : int, optional + Minimum number of adjacent True values in a region. + Defaults to half size of the series. + + Returns + ------- + list + A list of tuples representing the start and end indices of regions + containing more than or equal to min_adjacent adjacent True values. + """ + min_adjacent = min_adjacent if min_adjacent else 0.5 * series.size + regions = [] + for key, group in series.groupby((series != series.shift()).cumsum()): + if key and len(group) >= min_adjacent: + region_indices = group.index + regions.append((region_indices.min(), region_indices.max())) + return regions + def get_tma_slew_event(day_obs, seq_number): """ @@ -111,16 +320,16 @@ def get_tma_slew_event(day_obs, seq_number): the specified day of observation (dayObs). The events are filtered to include only those that start after 1 second before the specified start time and end before 1 second after the specified end time. - - Returs + + Returs ------ lsst.summit.utils.tmaUtils.TMAEvent - A TMA slew event that occurred within the specified time range. + A TMA slew event that occurred within the specified time range. Raises ------ ValueError - If no events are found for the specified time range. + If no events are found for the specified time range. ValueError If more than one event is found for the specified time range. """ @@ -142,7 +351,7 @@ def get_tma_slew_event(day_obs, seq_number): if len(single_event) == 0: raise ValueError(f"Could not find any events for {day_obs}. ") - return single_event[0] + return single_event[0] def evaluate_single_slew(day_obs, seq_number): @@ -159,7 +368,23 @@ def evaluate_single_slew(day_obs, seq_number): """ logger.info("Retriving TMA slew event.") event = get_tma_slew_event(day_obs, seq_number) - + + logger.info("Start inertia compensation system analysis.") + performance_analysis = InertiaCompensationSystemAnalysis( + event, + inner_pad=1.0, + outer_pad=1.0, + n_sigma=1.0, + ) + + logger.info("Query datasets") + performance_analysis.df = performance_analysis.query_dataset() + + logger.info("Calculate statistics") + performance_analysis.stats = performance_analysis.get_stats() + + print(performance_analysis.stats) + if __name__ == "__main__": logger.info("Start")