Flow. Use new method to interpolate stacking move-up events: match expected block's height

cfatools/processor/flow.py

@@ -270,7 +270,7 @@ def get_datasets_data(dataset: DatasetReader, retrieve_iceblocks_from_processed_
             os.mkdir(processed_dirname)
         iceblock_df = dataset.get_timeseries("ICBKCTRL_instant")
         iceblock_df = iceblock.get_clean_iceblock_df(iceblock_df)
-        iceblock_df.to_csv(iceblock_filename, index=False, date_format="%Y-%m-%dT%H:%M:%SZ")
+        iceblock_df.to_csv(iceblock_filename, index=False, date_format="%Y-%m-%d %H:%M:%S.%fZ")
 
     # Conductivity
     conduct_df = dataset.get_timeseries("CONDUCTI_periodic")
@@ -391,33 +391,21 @@ def get_arrival_df(encoder_df: pd.DataFrame,
     return arrival_df
 
 
-def add_iceblock_info(arrival_df: pd.DataFrame, encoder_df: pd.DataFrame, iceblock_df: pd.DataFrame,
-                      normal_speed_range: typing.Tuple[float, float] = (-3.0, -0.01)) -> pd.DataFrame:
-    """Add iceblock information: melted block height and iceblock code & name.
+def add_iceblock_info(arrival_df: pd.DataFrame, iceblock_df: pd.DataFrame) -> pd.DataFrame:
+    """Add iceblock code & name.
 
     Parameters
     ----------
     arrival_df: pd.DataFrame
         Output of :func:`get_arrival_df` function.
-    encoder_df: pd.DataFrame
-        Encoder's timestamped data, as obtained with :func:`get_datasets_data` function.
     iceblock_df: pd.DataFrame
         Iceblock's timestamped data, as obtained with :func:`get_datasets_data` function.
-    normal_speed_range: typing.Tuple[float, float]
-        Min and max encoder speed for "normal melting". Negative speed means encoder is going down. Speed unit: mm/s
 
     Returns
     -------
     pd.DataFrame
-        Same as input `arrival_df`, with melted block height and iceblock code & name.
+        Same as input `arrival_df`, with iceblock's code & name.
     """
-    moveup_events_df = encoder.get_moveup_events(encoder_df, normal_speed_range)
-    stacked_iceblock_df = iceblock.get_total_stacked_height(iceblock_df)
-    encoder_df = get_absolute_melted_height(encoder_df, stacked_iceblock_df, moveup_events_df,
-                                            starting_on_plexi=True)
-
-    arrival_df = pd.merge(arrival_df, encoder_df[["melted_height"]], left_index=True, right_index=True)
-
     arrival_df = pd.merge_asof(arrival_df,
                                iceblock_df[["id", "name", "datetime_start"]],
                                left_index=True,
@@ -426,11 +414,62 @@ def add_iceblock_info(arrival_df: pd.DataFrame, encoder_df: pd.DataFrame, iceblo
     arrival_df = arrival_df.drop(columns="datetime_start")
     arrival_df = arrival_df.rename(columns={'id': 'icbk_code', 'name': 'icbk_name'})
 
-    def _get_melted_height_by_iceblock_(df):
-        df["melted_height_icbk"] = df["melted_height"] - df.iloc[0]["melted_height"]
+    return arrival_df
+
+def add_melted_height(arrival_df: pd.DataFrame, encoder_df: pd.DataFrame, iceblock_df: pd.DataFrame,
+                      normal_speed_range: typing.Tuple[float, float] = (-3.0, -0.01)) -> pd.DataFrame:
+    """Add melted block height.
+
+    Parameters
+    ----------
+    arrival_df: pd.DataFrame
+        Output of :func:`get_arrival_df` function, assuming that :func:`add_iceblock_info` has been executed beforehand.
+    encoder_df: pd.DataFrame
+        Encoder's timestamped data, as obtained with :func:`get_datasets_data` function.
+    iceblock_df: pd.DataFrame
+        Iceblock's timestamped data, as obtained with :func:`get_datasets_data` function.
+    normal_speed_range: typing.Tuple[float, float]
+        Min and max encoder speed for "normal melting". Negative speed means encoder is going down. Speed unit: mm/s
+
+    Returns
+    -------
+    pd.DataFrame
+        Same as input `arrival_df`, with melted block height.
+    """
+    moveup_events_df = encoder.get_moveup_events(encoder_df, normal_speed_range)
+    stacked_iceblock_df = iceblock.get_total_stacked_height(iceblock_df)
+    arrival_df = get_absolute_melted_height_interp(arrival_df, stacked_iceblock_df, moveup_events_df,
+                                                   starting_on_plexi=True)
+    arrival_df = get_absolute_melted_height_expected(arrival_df, iceblock_df, moveup_events_df,
+                                                 starting_on_plexi=True)
+
+    def _get_melted_height_by_iceblock_(df, suffix: str):
+        df["melted_height_" + suffix + "_icbk"] = df["melted_height_" + suffix] - df.iloc[0]["melted_height_" + suffix]
         return df
 
-    arrival_df = arrival_df.groupby("icbk_code").apply(_get_melted_height_by_iceblock_)
+    arrival_df = arrival_df.groupby("icbk_code").apply(_get_melted_height_by_iceblock_, ("interpstack"))
+    arrival_df = arrival_df.groupby("icbk_code").apply(_get_melted_height_by_iceblock_, ("expected"))
+
+    # DEBUG
+    # comp_df = pd.merge(arrival_df.groupby("icbk_code").max(),
+    #                    iceblock_df,
+    #                    left_index=True,
+    #                    right_on="id")[["id",
+    #                                    "initial_height",
+    #                                    "melted_height_interpstack_icbk",
+    #                                    "melted_height_expected_icbk"]]
+    # comp_df["diff_interpstack"] = comp_df["initial_height"] - comp_df["melted_height_interpstack_icbk"]
+    # comp_df["diff_expected"] = comp_df["initial_height"] - comp_df["melted_height_expected_icbk"]
+    # arrival_df.loc[abs(arrival_df["position_diff"]) < 1, "position_diff"].hist(bins=100)
+    # print(comp_df["diff_expected"].sum())
+    # print(comp_df["diff_expected"].mean())
+    # print(comp_df["diff_expected"].std())
+    # /DEBUG
+
+    # Keep only the columns related to the finally-selected method (matching expected height)
+    arrival_df = arrival_df.drop(columns=["melted_height_interpstack_icbk", "melted_height_interpstack"])
+    arrival_df = arrival_df.rename(columns={"melted_height_expected_icbk": "melted_height_icbk",
+                                            "melted_height_expected": "melted_height"})
 
     return arrival_df
 
@@ -544,16 +583,19 @@ def __compute_mm__(df: pd.DataFrame) -> pd.DataFrame:
     return df
 
 
-def get_absolute_melted_height(encoder_df: pd.DataFrame,
+def get_absolute_melted_height_interp(arrival_df: pd.DataFrame,
                                stacked_iceblocks_df: pd.DataFrame,
                                moveup_event_df: pd.DataFrame,
                                starting_on_plexi: bool) -> pd.DataFrame:
     """Get the total absolute melting height, for each encoder's timestep.
 
+    The "holes" caused by encoder's move-up events are filled by interoplating the underlying melting speed using
+    [height of the ice melted during the move-up event] / [duration of the move-up event].
+
     Parameters
     ----------
-    encoder_df: pd.DataFrame
-        ``ENCODER_periodic``'s dataset.
+    arrival_df: pd.DataFrame
+        Output of :func:`get_arrival_df` function, assuming that :func:`add_iceblock_info` has been executed beforehand.
     stacked_iceblocks_df: pd.DataFrame
         Datetime-indexed DataFrame containing the total height of ice blocks stacked. In other words, the height of all
         the ice blocks stacked at the same time should be summed.
@@ -566,10 +608,10 @@ def get_absolute_melted_height(encoder_df: pd.DataFrame,
     Returns
     -------
     pd.DataFrame
-        Same as ``encoder_df``, with an additional ``melted_height`` column.
+        Same as ``arrival_df``, with an additional ``melted_height`` column.
     """
-    encoder_df["time_diff"] = encoder_df.index.to_series().diff().dt.total_seconds()
-    encoder_df["ice_speed"] = encoder_df["position"].diff() / encoder_df["time_diff"]
+    arrival_df["time_diff"] = arrival_df.index.to_series().diff().dt.total_seconds()
+    arrival_df["ice_speed"] = arrival_df["position"].diff() / arrival_df["time_diff"]
 
     if len(moveup_event_df.index) != 0:
         stack_events_df = moveup_event_df[moveup_event_df["event_type"] == "stacking"]
@@ -580,8 +622,8 @@ def get_absolute_melted_height(encoder_df: pd.DataFrame,
                                             left_on="start_datetime", right_on="datetime_stacked",
                                             direction="forward")
 
-            # Compute the height of already-stacked ice which has been melted during the stacking event, while the
-            # encoder was moving or at its parking position
+            # Compute the height of already-stacked ice which has been melted during the stacking event, while the encoder was
+            # moving or at its parking position
             stack_events_df["melted_while_event"] = \
                 stack_events_df["tot_stacked_height"] - (stack_events_df["end_position"] - stack_events_df["start_position"])
 
@@ -591,27 +633,105 @@ def get_absolute_melted_height(encoder_df: pd.DataFrame,
 
             # Compute speed for each timestep of each event
             for index, stack_event in stack_events_df.iterrows():
-                encoder_df.loc[stack_event["start_datetime"]:stack_event["end_datetime"], "ice_speed"] = \
+                arrival_df.loc[stack_event["start_datetime"]:stack_event["end_datetime"], "ice_speed"] = \
                     stack_event["avg_melting_speed"]
 
         endofmelting_events_df = moveup_event_df[moveup_event_df["event_type"] == "end_of_melting"].reset_index(drop=True)
         if len(endofmelting_events_df.index) == 1:  # Keep only the encoder data which occured before the end-of-melting
-            encoder_df = encoder_df[:endofmelting_events_df.loc[0, "start_datetime"]].copy()
+            arrival_df = arrival_df[:endofmelting_events_df.loc[0, "start_datetime"]].copy()
         elif len(endofmelting_events_df.index) > 1:
             raise ValueError("More than one end-of-melting move-up events!")
 
-    encoder_df["melted_height"] = (-encoder_df["ice_speed"] * encoder_df["time_diff"]).cumsum()
+    arrival_df["melted_height"] = (-arrival_df["ice_speed"] * arrival_df["time_diff"]).cumsum()
 
     if starting_on_plexi:
-        encoder_df["melted_height"] = encoder_df["melted_height"] - 7.5
-        encoder_df.loc[encoder_df["melted_height"] < 0, "melted_height"] = 0
+        arrival_df["melted_height"] = arrival_df["melted_height"] - 7.5
+        arrival_df.loc[arrival_df["melted_height"] < 0, "melted_height"] = 0
 
     # Set all-but-last 0-value melted height to NaN
-    if len(encoder_df[encoder_df["melted_height"] == 0].index) > 2:
-        encoder_df.loc[:encoder_df[encoder_df["melted_height"] == 0].tail(2).index[0], "melted_height"] = np.nan
+    if len(arrival_df[arrival_df["melted_height"] == 0].index) > 2:
+        arrival_df.loc[:arrival_df[arrival_df["melted_height"] == 0].tail(2).index[0], "melted_height"] = np.nan
+
+    arrival_df = arrival_df.drop(columns=["time_diff", "ice_speed"])
+    arrival_df = arrival_df.rename(columns={"melted_height": "melted_height_interpstack"})
+
+    return arrival_df
 
-    return encoder_df
 
+def get_absolute_melted_height_expected(arrival_df: pd.DataFrame,
+                                    iceblock_df: pd.DataFrame,
+                                    moveup_event_df: pd.DataFrame,
+                                    starting_on_plexi: bool) -> pd.DataFrame:
+    """Get the total absolute melting height, for each encoder's timestep.
+
+    The "holes" caused by encoder's move-up events are filled in a way that the final melted height of each block is
+    exactly equal to the `initial_height` the user defined when stacking the ice block.
+
+    Parameters
+    ----------
+    arrival_df: pd.DataFrame
+        Output of :func:`get_arrival_df` function, assuming that :func:`add_iceblock_info` has been executed beforehand.
+    iceblock_df: pd.DataFrame
+        Iceblock's timestamped data, as obtained with :func:`get_datasets_data` function.
+    moveup_event_df: pd.DataFrame
+        Output of :func:`cfatools.processor.encoder.get_moveup_events` function.
+    starting_on_plexi: bool
+        If ``True``, the plexi-plate is in ``PLACED`` position at the beginning of the ``encoder_df`` dataset. In this
+        case, 7.5mm (plexi-plate height) are removed from the ``melted_height``.
+
+    Returns
+    -------
+    pd.DataFrame
+        Same as ``arrival_df``, with an additional ``melted_height`` column.
+    """
+    arrival_df["time_diff"] = arrival_df.index.to_series().diff().dt.total_seconds()
+    arrival_df["position_diff"] = arrival_df["position"].diff()
+
+    arrival_df["ice_speed"] = arrival_df["position"].diff() / arrival_df["time_diff"]
+
+    if len(moveup_event_df.index) != 0:
+        stack_events_df = moveup_event_df[moveup_event_df["event_type"] == "stacking"]
+        if len(stack_events_df.index) != 0:
+            # Compute speed for each timestep of each event
+            for index, stack_event in stack_events_df.iterrows():
+                melting_block_code = arrival_df.loc[stack_event["start_datetime"]:stack_event["end_datetime"], "icbk_code"].unique()
+                if len(melting_block_code) == 0:
+                    raise ValueError("No melting block found")
+                elif len(melting_block_code) > 1:
+                    raise ValueError("More than melting block found")
+                else:
+                    melting_block_code = melting_block_code[0]
+                melting_block = iceblock_df[iceblock_df["id"] == melting_block_code].squeeze()
+                expected_block_height = iceblock_df.loc[iceblock_df["id"] == melting_block_code, "initial_height"].unique()[0]
+                melted_before = arrival_df.iloc[arrival_df.index.get_loc(melting_block["datetime_start"], method='nearest')]["position"]\
+                                - arrival_df.loc[stack_event["start_datetime"], "position"]
+                melted_after = arrival_df.loc[stack_event["end_datetime"], "position"] \
+                               - arrival_df.iloc[arrival_df.index.get_loc(melting_block["datetime_end"], method='nearest')]["position"]
+                melted_while_stak_event = expected_block_height - melted_before - melted_after
+                avg_melting_speed = melted_while_stak_event / (stack_event["end_datetime"]-stack_event["start_datetime"]).total_seconds()
+                arrival_df.loc[stack_event["start_datetime"]:stack_event["end_datetime"], "ice_speed"] = -avg_melting_speed
+
+        endofmelting_events_df = moveup_event_df[moveup_event_df["event_type"] == "end_of_melting"].reset_index(
+            drop=True)
+        if len(endofmelting_events_df.index) == 1:  # Keep only the encoder data which occured before the end-of-melting
+            arrival_df = arrival_df[:endofmelting_events_df.loc[0, "start_datetime"]].copy()
+        elif len(endofmelting_events_df.index) > 1:
+            raise ValueError("More than one end-of-melting move-up events!")
+
+    arrival_df["melted_height"] = (-arrival_df["ice_speed"] * arrival_df["time_diff"]).cumsum()
+
+    if starting_on_plexi:
+        arrival_df["melted_height"] = arrival_df["melted_height"] - 7.5
+        arrival_df.loc[arrival_df["melted_height"] < 0, "melted_height"] = 0
+
+    # Set all-but-last 0-value melted height to NaN
+    if len(arrival_df[arrival_df["melted_height"] == 0].index) > 2:
+        arrival_df.loc[:arrival_df[arrival_df["melted_height"] == 0].tail(2).index[0], "melted_height"] = np.nan
+
+    arrival_df = arrival_df.drop(columns=["time_diff", "ice_speed"])
+    arrival_df = arrival_df.rename(columns={"melted_height": "melted_height_expected"})
+
+    return arrival_df
 
 def get_tubing_volume_dict(filepath: str,
                            max_datetime: datetime.datetime = datetime.datetime.now(datetime.timezone.utc)) -> dict: