temperature-based-fertility-prediction-thesis/code/utils/feature_functions.py

from datetime import timedelta
from typing import Callable

import numpy as np
from vsm_datascience_common import constants
from vsm_datascience_common.cycle_database_connection.db_utils import get_cycles_collection
from vsm_datascience_common.cycles.sequences import get_timestamps, get_values

from utils.smoothing import get_curve_composition


def get_day_of_week_encoded(cycle: dict, shift: int = 0) -> dict:
    timestamps = get_timestamps(cycle)
    weekdays = np.array([x.weekday() for x in timestamps])
    sine_encoded = np.sin(weekdays * (2 * np.pi / 7))
    cosine_encoded = np.cos(weekdays * (2 * np.pi / 7))
    return {
        "sine": sine_encoded[shift if shift > 0 else 0:shift if shift < 0 else None],
        "cosine": cosine_encoded[shift if shift > 0 else 0:shift if shift < 0 else None]
    }


def get_hour_of_day_encoded(cycle: dict, shift: int = 0) -> dict:
    timestamps = get_timestamps(cycle)
    hours = np.array([x.hour for x in timestamps])
    sine_encoded = np.sin(hours * (2 * np.pi / 24))
    cosine_encoded = np.cos(hours * (2 * np.pi / 24))
    return {
        "sine": sine_encoded[shift if shift > 0 else 0:shift if shift < 0 else None],
        "cosine": cosine_encoded[shift if shift > 0 else 0:shift if shift < 0 else None]
    }


def get_month_of_year_encoded(cycle: dict, shift: int = 0) -> dict:
    timestamps = get_timestamps(cycle)
    months = np.array([x.month for x in timestamps])
    sine_encoded = np.sin(months * (2 * np.pi / 12))
    cosine_encoded = np.cos(months * (2 * np.pi / 12))
    return {
        "sine": sine_encoded[shift if shift > 0 else 0:shift if shift < 0 else None],
        "cosine": cosine_encoded[shift if shift > 0 else 0:shift if shift < 0 else None]
    }


def get_hours_from_start(cycle: dict, shift: int = 0) -> np.ndarray:
    timestamps = get_timestamps(cycle)
    return np.arange(len(timestamps))[shift if shift > 0 else 0:shift if shift < 0 else None]


def get_temperature(cycle: dict, shift: int = 0) -> np.ndarray:
    return get_values(cycle)[shift if shift > 0 else 0:shift if shift < 0 else None]


def get_fertility_curve(cycle_length: int,
                        ov_index: int,
                        fertility_base_curve: np.ndarray | list,
                        base_curve_offset: int) -> np.ndarray:
    fert_prob = np.zeros((cycle_length,))
    if ov_index + base_curve_offset < 0:
        length_in_cycle = len(fertility_base_curve) + ov_index + base_curve_offset
        fert_prob[:length_in_cycle] = fertility_base_curve[-length_in_cycle:]
    elif ov_index + base_curve_offset + len(fertility_base_curve) > cycle_length:
        length_in_cycle = cycle_length - ov_index - base_curve_offset
        fert_prob[ov_index + base_curve_offset:] = fertility_base_curve[:length_in_cycle]
    else:
        fert_prob[ov_index + base_curve_offset:
                  ov_index + base_curve_offset + len(fertility_base_curve)] = fertility_base_curve
    return fert_prob


def get_fertility_probability_base_curve() -> (np.ndarray, int):
    """
    This function returns the base curve for the fertility probability as well as the offset to the ovulation day in measurements
    :return: numpy array with the base curve resampled to the number of measurements per day
    """
    fertilization_chance_curve = [0, 0.07, 0.12, 0.25, 0.3, 0.18, 0]
    offset_to_ov = -6

    new_indices = np.arange(len(fertilization_chance_curve) * constants.MEASUREMENTS_PER_DAY)
    fertility_chance_curve_resampled = np.interp(new_indices,
                                                 np.linspace(0, len(new_indices), num=len(fertilization_chance_curve)),
                                                 fertilization_chance_curve)
    ov_offset_resampled = offset_to_ov * constants.MEASUREMENTS_PER_DAY
    return fertility_chance_curve_resampled, ov_offset_resampled


def get_fertility_probability(cycle: dict, shift: int = 0) -> np.ndarray:
    # check, if biphasic
    if "classification_results" in cycle and "results" in cycle["classification_results"][0]:
        if cycle["classification_results"][0]["results"]["predicted_class"] != "biphasic":
            return np.zeros((len(get_values(cycle)) - abs(shift, )))
    else:
        return np.zeros((len(get_values(cycle)) - abs(shift, )))

    if "ov_detection_results" in cycle and "results" in cycle["ov_detection_results"][0]:
        try:
            ov_timestamp = cycle["ov_detection_results"][0]["results"]["ovulation_timestamp"]
            timestamps = get_timestamps(cycle)
            ov_index = np.where(timestamps >= ov_timestamp)[0][0]
            fertility_chance_curve_resampled, ov_offset_resampled = get_fertility_probability_base_curve()
            return get_fertility_curve(len(get_values(cycle)), ov_index, fertility_chance_curve_resampled,
                                       ov_offset_resampled)[shift if shift > 0 else 0:shift if shift < 0 else None]
        except:
            return np.zeros((len(get_values(cycle)) - abs(shift, )))
    else:
        return np.zeros((len(get_values(cycle)) - abs(shift, )))


def get_ov_over_probability(cycle: dict, shift: int = 0) -> np.ndarray:
    """
    Get the ovulation over probability for a cycle
    :param cycle: cycle data
    :param shift: shift to apply
    :return: ovulation over probability
    """
    if "ov_detection_results" in cycle and "results" in cycle["ov_detection_results"][0]:
        try:
            ov_timestamp = cycle["ov_detection_results"][0]["results"]["ovulation_timestamp"]
            timestamps = get_timestamps(cycle)
            ov_index = np.where(timestamps > ov_timestamp)[0][0]
            ov_over_probability = np.zeros(len(get_values(cycle)))
            ov_over_probability[ov_index:] = 1
            return ov_over_probability[shift if shift > 0 else 0:shift if shift < 0 else None]
        except:
            return np.zeros(len(get_values(cycle)))[shift if shift > 0 else 0:shift if shift < 0 else None]
    else:
        return np.zeros(len(get_values(cycle)))[shift if shift > 0 else 0:shift if shift < 0 else None]


def get_days_relative_to_ov(cycle: dict) -> np.ndarray:
    """
    Get the time in days until the next OV event.
    """
    timestamps = get_timestamps(cycle)
    if ("ov_detection_results" in cycle and "results" in cycle["ov_detection_results"][0]) and (
            "classification_results" in cycle and "results" in cycle["classification_results"][0] and
            cycle["classification_results"][0]["results"]["predicted_class"] == "biphasic"):
        try:
            ov_timestamp = cycle["ov_detection_results"][0]["results"]["ovulation_timestamp"]
            ov_index = np.where(timestamps > ov_timestamp)[0][0]
            ov_timestamp = timestamps[ov_index]
            days_until_ov = np.empty(len(timestamps))
            for i in range(len(timestamps)):
                current_day = timestamps[i]
                day_diff = ov_timestamp - current_day
                if day_diff.days < -1:
                    pass
                # days_until_ov[i] = max(day_diff.days -1)
                # use - days, as days relative to event
                days_until_ov[i] = - day_diff.days

            return np.array(days_until_ov)
        except:
            pass

    return np.full((len(timestamps),), np.nan)


def get_ov_day(cycle: dict) -> np.ndarray:
    """
    Get the ovulation day for a cycle
    :param cycle: cycle data
    :return: ovulation day
    """
    values = get_values(cycle)
    seq_len = len(values)
    is_biphasic = "classification_results" in cycle and "results" in cycle["classification_results"][0] and \
                  cycle["classification_results"][0]["results"]["predicted_class"] == "biphasic"
    has_ov_results = "ov_detection_results" in cycle and "results" in cycle["ov_detection_results"][0]
    if is_biphasic and has_ov_results:
        try:
            ov_day = cycle["ov_detection_results"][0]["results"]["ovulation_day"]
            return np.full((seq_len,), ov_day)
        except:
            pass

    return np.full((seq_len,), np.nan)


def get_is_biphasic(cycle: dict) -> np.ndarray:
    """
    Get the is_biphasic flag for a cycle
    :param cycle: cycle data
    :return: is_biphasic flag
    """
    values = get_values(cycle)
    seq_len = len(values)
    is_biphasic = "classification_results" in cycle and "results" in cycle["classification_results"][0] and \
                  cycle["classification_results"][0]["results"]["predicted_class"] == "biphasic"
    if is_biphasic:
        return np.full((seq_len,), 1)

    return np.full((seq_len,), 0)


def get_curve_composition_as_features(cycle: dict, shift: int = 0) -> dict:
    composition = get_curve_composition(get_values(cycle))
    return {
        "trend": composition[0][shift if shift > 0 else 0:shift if shift < 0 else None],
        "seasonal": composition[1][shift if shift > 0 else 0:shift if shift < 0 else None],
        "residual": composition[2][shift if shift > 0 else 0:shift if shift < 0 else None],
        "smoothed": composition[3][shift if shift > 0 else 0:shift if shift < 0 else None]
    }


def get_rolling_average_with_padding(cycle: dict,
                                     shift: int = 0,
                                     rolling_average_window_length: int = constants.MEASUREMENTS_PER_DAY) -> np.ndarray:
    """
    Get the rolling average of the temperature values with padding, so that the rolling average has the same length as the
    original temperature values
    Args:
        cycle (dict): cycle data
        shift (int): shift to apply
        rolling_average_window_length (int): window length for the rolling average, default is one day (288 measurements)

    Returns:
        np.ndarray: rolling average of the temperature values with padding

    """

    values = get_values(cycle)
    values_padded = np.concatenate(
        [np.full((constants.MEASUREMENTS_PER_DAY,), np.mean(values[:constants.MEASUREMENTS_PER_DAY])), values])
    rolling_average = np.convolve(values_padded, np.ones(rolling_average_window_length) / rolling_average_window_length,
                                  mode="valid")[-len(values):]
    return rolling_average[shift if shift > 0 else 0:shift if shift < 0 else None]


def get_window_fn(cycle: dict,
                  window_size: int,
                  fn: Callable) -> np.ndarray:
    """
    Apply a window function to a time series.
    :param cycle: cycle data
    :param window_size: window size
    :return: array with minimum value in window for each value in input array
    """
    values = get_values(cycle)
    values_padded = np.concatenate([np.full(window_size - 1, values[0]), values])

    windows = np.lib.stride_tricks.sliding_window_view(values_padded, window_shape=window_size)
    return fn(windows)


def get_cycle_length_stats(cycle: dict) -> dict:
    user_id = cycle["user_id"]
    cycle_length = len(get_values(cycle))
    cycles_dates = list(get_cycles_collection().find({
        "user_id": user_id, "ends_at": {"$exists": True},
        "starts_at": {"$lt": cycle["starts_at"]},
    },
        {"starts_at": 1, "ends_at": 1}))
    if len(cycles_dates) == 0:
        # use the average cycle length as fallback
        average_cycle_length = 37.13
        cycle_length_std = 12.42
    else:
        cycle_lengths = [x["ends_at"] - x["starts_at"] for x in cycles_dates]
        average_cycle_length = np.mean(cycle_lengths).total_seconds() / (60 * 60 * 24)
        cycle_length_std = np.std([x.total_seconds() for x in cycle_lengths]) / (60 * 60 * 24)

    return {
        "average_cycle_length": np.full((cycle_length,), average_cycle_length),
        "cycle_length_std": np.full((cycle_length,), cycle_length_std),
    }


def get_average_ovulation_day(cycle: dict) -> float:
    user_id = cycle["user_id"]
    cycle_length = len(get_values(cycle))
    cycles_dates = list(
        get_cycles_collection().find({
            "user_id": user_id,
            "ov_detection_results.0.results": {"$exists": True},
            "starts_at": {"$lt": cycle["starts_at"]},
        },
            {"ov_detection_results": 1}))
    if len(cycles_dates) == 0:
        # use the average ovulation day as fallback
        average_ovulation_day = 18.9
    else:
        ovulation_days = [x["ov_detection_results"][0]["results"]["ovulation_day"] for x in cycles_dates]
        average_ovulation_day = np.mean(ovulation_days)
    return np.full((cycle_length,), average_ovulation_day)


def get_ovulation_std(cycle: dict) -> float:
    user_id = cycle["user_id"]
    cycle_length = len(get_values(cycle))
    cycles_dates = list(
        get_cycles_collection().find({
            "user_id": user_id,
            "ov_detection_results.0.results": {"$exists": True},
            "starts_at": {"$lt": cycle["starts_at"]},
        },
            {"ov_detection_results": 1}))
    if len(cycles_dates) == 0:
        # use the average ovulation std as fallback
        ovulation_std = 4.07
    else:
        ovulation_days = [x["ov_detection_results"][0]["results"]["ovulation_day"] for x in cycles_dates]
        ovulation_std = np.std(ovulation_days)
    return np.full((cycle_length,), ovulation_std)


def get_biphasic_fraction(cycle: dict) -> float:
    user_id = cycle["user_id"]
    cycle_length = len(get_values(cycle))
    cycle_dates = list(
        get_cycles_collection().find({
            "user_id": user_id,
            "classification_results.0.results": {"$exists": True},
            "starts_at": {"$lt": cycle["starts_at"]},
        },
            {"classification_results": 1}))
    if len(cycle_dates) == 0:
        # use the average cycle length as fallback
        biphasic_fraction = 0.9449
    else:
        biphasic = [x for x in cycle_dates if
                    x["classification_results"][0]["results"]["predicted_class"] == "biphasic"]
        biphasic_fraction = len(biphasic) / len(cycle_dates)
    return np.full((cycle_length,), biphasic_fraction)


def get_num_cycles(cycle: dict) -> int:
    user_id = cycle["user_id"]
    cycle_length = len(get_values(cycle))
    cycles_dates = list(get_cycles_collection().find({"user_id": user_id,
                                                      "ends_at": {"$exists": True},
                                                      "starts_at": {"$lt": cycle["starts_at"]}
                                                      },
                                                     {"starts_at": 1, "ends_at": 1}))
    num_cycles = len(cycles_dates)
    return np.full((cycle_length,), num_cycles)


def get_average_temperatures(cycle: dict) -> dict:
    user_id = cycle["user_id"]
    cycle_length = len(get_values(cycle))
    cycles_dates = list(
        get_cycles_collection().find({"user_id": user_id, "ov_detection_results.0.results": {"$exists": True}},
                                     {"starts_at": 1, "ov_detection_results": 1, "measurements": 1}))

    if len(cycles_dates) == 0:
        return {
            "pre_ov_temperatures": np.full((cycle_length,), 37.07),
            "post_ov_temperatures": np.full((cycle_length,), 37.37),
        }

    pre_ov_temperatures = []
    post_ov_temperatures = []
    for x in cycles_dates:
        ov_detection = x["ov_detection_results"][0]["results"]
        values = get_values(x)
        timestamps = get_timestamps(x)
        cutoff_day = x["starts_at"] + timedelta(days=ov_detection["ovulation_day"])
        cutoff_index = np.where(np.array(timestamps) >= cutoff_day)[0][0]
        pre_ov_temperatures.append(values[:cutoff_index])
        post_ov_temperatures.append(values[cutoff_index:])

    pre_ov_temperatures = np.concatenate(pre_ov_temperatures)
    post_ov_temperatures = np.concatenate(post_ov_temperatures)

    pre_ov_temperatures = np.mean(pre_ov_temperatures)
    post_ov_temperatures = np.mean(post_ov_temperatures)

    return {
        "pre_ov_temperatures": np.full((cycle_length,), pre_ov_temperatures),
        "post_ov_temperatures": np.full((cycle_length,), post_ov_temperatures),
    }