def METHOD_NAME(self, time, mol_count, output_times):
    """
    Takes raw times and inputs and mol counts and outputs interpolated
    values at output_times
    """
    if len(time) == 0:
        raise ValueError('At least one time must be given.')
    if len(time) != len(mol_count):
        raise ValueError(
            'The number of entries in time must match mol_count')
    # Check output times
    output_times = np.asarray(output_times)
    if not np.all(output_times[1:] >= output_times[:-1]):
        raise ValueError('The output_times must be non-decreasing.')
    # Interpolate as step function, decreasing mol_count by 1 at each
    # reaction time point.
    if len(time) == 1:
        # Need at least 2 values to interpolate
        return np.ones(len(output_times)) * mol_count[0]
    else:
        # Note: Can't use fill_value='extrapolate' here as:
        #  1. This require scipy >= 0.17
        #  2. There seems to be a bug in some scipy versions
        interp_func = scipy.interpolate.interp1d(
            time, mol_count, kind='previous', axis=0, fill_value=np.nan,
            bounds_error=False)
        values = interp_func(output_times)
    # At any point past the final time, repeat the last value
    values[output_times >= time[-1]] = mol_count[-1]
    values[output_times < time[0]] = mol_count[0]
    return values