def compute_median_survival_time(times, surv_function):
    """
    Computes a median survival time estimate by looking for where the survival
    function crosses 1/2.

    Parameters
    ----------
    times : 1D numpy array
        Sorted list of unique times (in ascending order).

    surv_function : 1D numpy array
        A survival function evaluated at each of time in `times`, in the same
        order.

    Returns
    -------
    output : float
        Median survival time estimate.
    """
    t_left = times[0]
    t_right = times[-1]

    if surv_function[-1] > 1/2:
        # survival function never crosses 1/2; just output this last time point
        return t_right

    for t, s in zip(times, surv_function):
        if s >= 0.5:
            t_left = t

    for t, s in zip(reversed(times), reversed(surv_function)):
        if s <= 0.5:
            t_right = t
    return (t_left + t_right) / 2.