<commit_before>import datetime
from dateutil.parser import parse
import math

from openapscontrib.mmhistorytools.models import Unit


class Schedule(object):
    def __init__(self, entries):
        """

        :param entries:
        :type entries: list(dict)
        :return:
        :rtype:
        """
        self.entries = entries

    def at(self, time):
        """

        :param time:
        :type time: datetime.time
        :return:
        :rtype: dict
        """
        result = {}

        for entry in self.entries:
            if parse(entry['start']).time() > time:
                break
            result = entry

        return result


def carb_effect_curve(t, absorption_time, sensor_delay):
    """Returns the fraction of total carbohydrate effect with a given absorption time on blood
    glucose at the specified number of minutes after eating.

    This is the integral of Carbs on Board (COB), defined by a Scheiner GI curve from Think Link a
    Pancreas, fig 7-8. This is based on an algorithm that first appeared in GlucoDyn

    See: https://github.com/kenstack/GlucoDyn

    :param t: The time in t since the carbs were eaten
    :type t: float
    :param absorption_time: The total absorption time of the carbohydrates
    :type absorption_time: int
    :return: A percentage of the initial carb intake, from 0 to 1
    :rtype: float
    """

    # time delay between carb consumption and body/sensor reaction
    # sensor_delay = 10

    if t <= sensor_delay:
        return 0.0
    elif t <= absorption_time / 2.0 + sensor_delay:
        return 2.0 / (absorption_time ** 2) * ((t - sensor_delay) ** 2)
    elif t < absorption_time + sensor_delay:
        return -1.0 + 4.0 / absorption_time * ((t-sensor_delay) - (t - sensor_delay) ** 2 / (2.0 * absorption_time))
    else:
        return 1.0


def walsh_iob_curve(t, insulin_action_duration, sensor_delay):
    """Returns the fraction of a single insulin dosage remaining at the specified number of minutes
    after delivery; also known as Insulin On Board (IOB).

    This is a Walsh IOB curve, and is based on an algorithm that first appeared in GlucoDyn

    See: https://github.com/kenstack/GlucoDyn

    :param t: time in minutes since the dose began
    :type t: Int
    :param insulin_action_duration: The duration of insulin action (DIA) of the patient, in minutes
    :type insulin_action_duration: Int
    :return: The fraction of a insulin dosage remaining at the specified time
    :rtype: float
    """
    assert insulin_action_duration in (3 * 60, 4 * 60, 5 * 60, 6 * 60)
    iob = 0
    
    # time delay between insulin delivery and body/sensor reaction
    # sensor_delay = 10	

    if t >= insulin_action_duration + sensor_delay:
        iob = 0.0
    elif t <= sensor_delay:
        iob = 1.0
    elif insulin_action_duration == 3 * 60:
        t -= sensor_delay
        iob = -3.2030e-9 * (t**4) + 1.354e-6 * (t**3) - 1.759e-4 * (t**2) + 9.255e-4 * t + 0.99951
    elif insulin_action_duration == 4 * 60:
	t -= sensor_delay
        iob = -3.310e-10 * (t**4) + 2.530e-7 * (t**3) - 5.510e-5 * (t**2) - 9.086e-4 * t + 0.99950
    elif insulin_action_duration == 5 * 60:
        t -= sensor_delay
        iob = -2.950e-10 * (t**4) + 2.320e-7 * (t**3) - 5.550e-5 * (t**2) + 4.490e-4 * t + 0.99300
    elif insulin_action_duration == 6 * 60:
        t -= sensor_delay
        iob = -1.493e-10 * (t**4) + 1.413e-7 * (t**3) - 4.095e-5 * (t**2) + 6.365e-4 * t + 0.99700

    return iob


def integrate_iob(t0, t1, insulin_action_duration, t, sensor_delay):
    """Integrates IOB using Simpson's rule for spread-out (basal-like) doses

    TODO: Clean this up and use scipy.integrate.simps

    :param t0: The start time in minutes of the dose
    :type t0: float
    :param t1: The end time in minutes of the dose
    :type t1: float
    :param insulin_action_duration: The duration of insulin action (DIA) of the patient, in minutes
    :type insulin_action_duration: Int
    :param t: The current time in minutes
    :type t: float
    :return:
    :rtype: float
    """
    nn = 50  # nn needs to be even

    # initialize with first and last terms of simpson series
    dx = (t1 - t0) / nn
    integral = walsh_iob_curve(t - t0, insulin_action_duration, sensor_delay) + walsh_iob_curve(t - t1, insulin_action_duration, sensor_delay)

    for i in range(1, nn - 1, 2):
        integral = integral + 4 * walsh_iob_curve(t - (t0 + i * dx), insulin_action_duration, sensor_delay) + 2 * walsh_iob_curve(t - (t0 + (i + 1) * dx), insulin_action_duration, sensor_delay)

    integral = integral * dx / 3.0
    return integral


def bolus_effect_at_datetime(event, t, insulin_sensitivity, insulin_action_duration, sensor_delay):
    return -event['amount'] * insulin_sensitivity * (1 - walsh_iob_curve(t, insulin_action_duration * 60.0, sensor_delay))


def carb_effect_at_datetime(event, t, insulin_sensitivity, carb_ratio, absorption_rate, sensor_delay):
    return insulin_sensitivity / carb_ratio * event['amount'] * carb_effect_curve(t, absorption_rate, sensor_delay)


def temp_basal_effect_at_datetime(event, t, t0, t1, insulin_sensitivity, insulin_action_duration, sensor_delay):
    int_iob = integrate_iob(t0, t1, insulin_action_duration * 60.0, t, sensor_delay)

    return -event['amount'] / 60.0 * insulin_sensitivity * ((t1 - t0) - int_iob)


def future_glucose(
    normalized_history,
    recent_glucose,
    insulin_action_curve,
    insulin_sensitivity_schedule,
    carb_ratio_schedule,
    dt=5,
    sensor_delay=10,
    basal_dosing_end=None
):
    if len(recent_glucose) == 0:
        return []

    last_glucose_entry = recent_glucose[0]
    last_glucose_value = last_glucose_entry.get('sgv') or last_glucose_entry.get('amount') or last_glucose_entry['glucose']
    last_glucose_datetime = parse(last_glucose_entry.get('date') or last_glucose_entry['display_time'])

    # Determine our simulation time.
    simulation_start = last_glucose_datetime
    simulation_end = last_glucose_datetime

    if len(normalized_history) > 0:
        last_history_event = sorted(normalized_history, key=lambda e: e['end_at'])[-1]
        last_history_datetime = parse(last_history_event['end_at'])
        simulation_end = max(simulation_end, last_history_datetime)

    simulation_end += datetime.timedelta(minutes=(insulin_action_curve * 60 + sensor_delay))

    # For each incremental minute from the simulation start time, calculate the effect values
    simulation_minutes = range(0, int(math.ceil((simulation_end - simulation_start).total_seconds() / 60.0)) + dt, dt)
    simulation_timestamps = [simulation_start + datetime.timedelta(minutes=m) for m in simulation_minutes]
    simulation_count = len(simulation_minutes)

    carb_effect = [0.0] * simulation_count
    insulin_effect = [0.0] * simulation_count

    for history_event in normalized_history:
        initial_effect = 0
        start_at = parse(history_event['start_at'])
        end_at = parse(history_event['end_at'])

        insulin_end_datetime = end_at + datetime.timedelta(hours=insulin_action_curve)
        absorption_rate = 180
        absorption_end_datetime = end_at + datetime.timedelta(minutes=absorption_rate)

        for i, timestamp in enumerate(simulation_timestamps):
            t = (timestamp - start_at).total_seconds() / 60.0

            # Cap the time used to determine the sensitivity so it doesn't fluctuate
            # after completion
            sensitivity_time = min(insulin_end_datetime, timestamp)
            insulin_sensitivity = insulin_sensitivity_schedule.at(sensitivity_time.time())['sensitivity']

            if history_event['unit'] == Unit.grams:
                # Cap the time used to determine the carb ratio to absorption end so it doesn't
                # fluctuate after completion
                ratio_time = min(absorption_end_datetime, timestamp)
                carb_ratio = carb_ratio_schedule.at(ratio_time.time())['ratio']

                effect = carb_effect_at_datetime(history_event, t, insulin_sensitivity, carb_ratio, absorption_rate, sensor_delay)
                apply_to = carb_effect
            elif history_event['unit'] == Unit.units:
                effect = bolus_effect_at_datetime(history_event, t, insulin_sensitivity, insulin_action_curve, sensor_delay)
                apply_to = insulin_effect
            elif history_event['unit'] == Unit.units_per_hour:
                end_at = parse(history_event['end_at'])

                if history_event['type'] == 'TempBasal' and basal_dosing_end and end_at > basal_dosing_end:
                    end_at = basal_dosing_end

                t1 = (end_at - start_at).total_seconds() / 60.0

                effect = temp_basal_effect_at_datetime(history_event, t, 0, t1, insulin_sensitivity, insulin_action_curve, sensor_delay)
                apply_to = insulin_effect
            else:
                raise ValueError('Unknown event %s', history_event)

            if i == 0:
                initial_effect = effect

            effect -= initial_effect
            apply_to[i] += effect

    return [{
        'date': timestamp.isoformat(),
        'glucose': last_glucose_value + carb_effect[i] + insulin_effect[i]
    } for i, timestamp in enumerate(simulation_timestamps)]
<commit_msg>Moving sensor_delay out of walsh eq<commit_after>import datetime
from dateutil.parser import parse
import math

from openapscontrib.mmhistorytools.models import Unit


class Schedule(object):
    def __init__(self, entries):
        """

        :param entries:
        :type entries: list(dict)
        :return:
        :rtype:
        """
        self.entries = entries

    def at(self, time):
        """

        :param time:
        :type time: datetime.time
        :return:
        :rtype: dict
        """
        result = {}

        for entry in self.entries:
            if parse(entry['start']).time() > time:
                break
            result = entry

        return result


def carb_effect_curve(t, absorption_time, sensor_delay):
    """Returns the fraction of total carbohydrate effect with a given absorption time on blood
    glucose at the specified number of minutes after eating.

    This is the integral of Carbs on Board (COB), defined by a Scheiner GI curve from Think Link a
    Pancreas, fig 7-8. This is based on an algorithm that first appeared in GlucoDyn

    See: https://github.com/kenstack/GlucoDyn

    :param t: The time in t since the carbs were eaten
    :type t: float
    :param absorption_time: The total absorption time of the carbohydrates
    :type absorption_time: int
    :return: A percentage of the initial carb intake, from 0 to 1
    :rtype: float
    """

    # time delay between carb consumption and body/sensor reaction
    # sensor_delay = 10

    if t <= sensor_delay:
        return 0.0
    elif t <= absorption_time / 2.0 + sensor_delay:
        return 2.0 / (absorption_time ** 2) * ((t - sensor_delay) ** 2)
    elif t < absorption_time + sensor_delay:
        return -1.0 + 4.0 / absorption_time * ((t-sensor_delay) - (t - sensor_delay) ** 2 / (2.0 * absorption_time))
    else:
        return 1.0


def walsh_iob_curve(t, insulin_action_duration):
    """Returns the fraction of a single insulin dosage remaining at the specified number of minutes
    after delivery; also known as Insulin On Board (IOB).

    This is a Walsh IOB curve, and is based on an algorithm that first appeared in GlucoDyn

    See: https://github.com/kenstack/GlucoDyn

    :param t: time in minutes since the dose began
    :type t: Int
    :param insulin_action_duration: The duration of insulin action (DIA) of the patient, in minutes
    :type insulin_action_duration: Int
    :return: The fraction of a insulin dosage remaining at the specified time
    :rtype: float
    """
    assert insulin_action_duration in (3 * 60, 4 * 60, 5 * 60, 6 * 60)
    iob = 0

    if t >= insulin_action_duration:
        iob = 0.0
    elif t <= 0:
        iob = 1.0
    elif insulin_action_duration == 3 * 60:
        iob = -3.2030e-9 * (t**4) + 1.354e-6 * (t**3) - 1.759e-4 * (t**2) + 9.255e-4 * t + 0.99951
    elif insulin_action_duration == 4 * 60:
        iob = -3.310e-10 * (t**4) + 2.530e-7 * (t**3) - 5.510e-5 * (t**2) - 9.086e-4 * t + 0.99950
    elif insulin_action_duration == 5 * 60:
        iob = -2.950e-10 * (t**4) + 2.320e-7 * (t**3) - 5.550e-5 * (t**2) + 4.490e-4 * t + 0.99300
    elif insulin_action_duration == 6 * 60:
        iob = -1.493e-10 * (t**4) + 1.413e-7 * (t**3) - 4.095e-5 * (t**2) + 6.365e-4 * t + 0.99700

    return iob


def integrate_iob(t0, t1, insulin_action_duration, t, sensor_delay):
    """Integrates IOB using Simpson's rule for spread-out (basal-like) doses

    TODO: Clean this up and use scipy.integrate.simps

    :param t0: The start time in minutes of the dose
    :type t0: float
    :param t1: The end time in minutes of the dose
    :type t1: float
    :param insulin_action_duration: The duration of insulin action (DIA) of the patient, in minutes
    :type insulin_action_duration: Int
    :param t: The current time in minutes
    :type t: float
    :return:
    :rtype: float
    """
    nn = 50  # nn needs to be even

    # initialize with first and last terms of simpson series
    dx = (t1 - t0) / nn
    integral = walsh_iob_curve(t - t0 - sensor_delay, insulin_action_duration) + walsh_iob_curve(t - t1 - sensor_delay, insulin_action_duration)

    for i in range(1, nn - 1, 2):
        integral = integral + 4 * walsh_iob_curve(t - (t0 + i * dx) - sensor_delay, insulin_action_duration) + 2 * walsh_iob_curve(t - (t0 + (i + 1) * dx) - sensor_delay, insulin_action_duration)

    integral = integral * dx / 3.0
    return integral


def bolus_effect_at_datetime(event, t, insulin_sensitivity, insulin_action_duration, sensor_delay):
    return -event['amount'] * insulin_sensitivity * (1 - walsh_iob_curve(t - sensor_delay, insulin_action_duration * 60.0))


def carb_effect_at_datetime(event, t, insulin_sensitivity, carb_ratio, absorption_rate, sensor_delay):
    return insulin_sensitivity / carb_ratio * event['amount'] * carb_effect_curve(t, absorption_rate, sensor_delay)


def temp_basal_effect_at_datetime(event, t, t0, t1, insulin_sensitivity, insulin_action_duration, sensor_delay):
    int_iob = integrate_iob(t0, t1, insulin_action_duration * 60.0, t, sensor_delay)

    return -event['amount'] / 60.0 * insulin_sensitivity * ((t1 - t0) - int_iob)


def future_glucose(
    normalized_history,
    recent_glucose,
    insulin_action_curve,
    insulin_sensitivity_schedule,
    carb_ratio_schedule,
    dt=5,
    sensor_delay=10,
    basal_dosing_end=None
):
    if len(recent_glucose) == 0:
        return []

    last_glucose_entry = recent_glucose[0]
    last_glucose_value = last_glucose_entry.get('sgv') or last_glucose_entry.get('amount') or last_glucose_entry['glucose']
    last_glucose_datetime = parse(last_glucose_entry.get('date') or last_glucose_entry['display_time'])

    # Determine our simulation time.
    simulation_start = last_glucose_datetime
    simulation_end = last_glucose_datetime

    if len(normalized_history) > 0:
        last_history_event = sorted(normalized_history, key=lambda e: e['end_at'])[-1]
        last_history_datetime = parse(last_history_event['end_at'])
        simulation_end = max(simulation_end, last_history_datetime)

    simulation_end += datetime.timedelta(minutes=(insulin_action_curve * 60 + sensor_delay))

    # For each incremental minute from the simulation start time, calculate the effect values
    simulation_minutes = range(0, int(math.ceil((simulation_end - simulation_start).total_seconds() / 60.0)) + dt, dt)
    simulation_timestamps = [simulation_start + datetime.timedelta(minutes=m) for m in simulation_minutes]
    simulation_count = len(simulation_minutes)

    carb_effect = [0.0] * simulation_count
    insulin_effect = [0.0] * simulation_count

    for history_event in normalized_history:
        initial_effect = 0
        start_at = parse(history_event['start_at'])
        end_at = parse(history_event['end_at'])

        insulin_end_datetime = end_at + datetime.timedelta(hours=insulin_action_curve)
        absorption_rate = 180
        absorption_end_datetime = end_at + datetime.timedelta(minutes=absorption_rate)

        for i, timestamp in enumerate(simulation_timestamps):
            t = (timestamp - start_at).total_seconds() / 60.0

            # Cap the time used to determine the sensitivity so it doesn't fluctuate
            # after completion
            sensitivity_time = min(insulin_end_datetime, timestamp)
            insulin_sensitivity = insulin_sensitivity_schedule.at(sensitivity_time.time())['sensitivity']

            if history_event['unit'] == Unit.grams:
                # Cap the time used to determine the carb ratio to absorption end so it doesn't
                # fluctuate after completion
                ratio_time = min(absorption_end_datetime, timestamp)
                carb_ratio = carb_ratio_schedule.at(ratio_time.time())['ratio']

                effect = carb_effect_at_datetime(history_event, t, insulin_sensitivity, carb_ratio, absorption_rate, sensor_delay)
                apply_to = carb_effect
            elif history_event['unit'] == Unit.units:
                effect = bolus_effect_at_datetime(history_event, t, insulin_sensitivity, insulin_action_curve, sensor_delay)
                apply_to = insulin_effect
            elif history_event['unit'] == Unit.units_per_hour:
                end_at = parse(history_event['end_at'])

                if history_event['type'] == 'TempBasal' and basal_dosing_end and end_at > basal_dosing_end:
                    end_at = basal_dosing_end

                t1 = (end_at - start_at).total_seconds() / 60.0

                effect = temp_basal_effect_at_datetime(history_event, t, 0, t1, insulin_sensitivity, insulin_action_curve, sensor_delay)
                apply_to = insulin_effect
            else:
                raise ValueError('Unknown event %s', history_event)

            if i == 0:
                initial_effect = effect

            effect -= initial_effect
            apply_to[i] += effect

    return [{
        'date': timestamp.isoformat(),
        'glucose': last_glucose_value + carb_effect[i] + insulin_effect[i]
    } for i, timestamp in enumerate(simulation_timestamps)]

{'fraction_non_alphanumeric': 0.09451302409048888, 'fraction_numerical': 0.026374462994181305, 'mean_word_length': 2.9144316730523627, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7228242', 'n_tokens_mistral': 6252, 'n_tokens_neox': 5583, 'n_words': 1595}