<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def otsu(data, min_threshold=None, max_threshold=None,bins=256):
 """Compute a threshold using Otsu's method data - an array of intensity values between zero and one min_threshold - only consider thresholds above this minimum value max_threshold - only consider thresholds below this maximum value bins - we bin the data into this many equally-spaced bins, then pick the bin index that optimizes the metric """

    assert min_threshold is None or max_threshold is None or min_threshold < max_threshold
    def constrain(threshold):
        if not min_threshold is None and threshold < min_threshold:
            threshold = min_threshold
        if not max_threshold is None and threshold > max_threshold:
            threshold = max_threshold
        return threshold
    
    data = np.atleast_1d(data)
    data = data[~ np.isnan(data)]
    if len(data) == 0:
        return (min_threshold if not min_threshold is None
                else max_threshold if not max_threshold is None
                else 0)
    elif len(data) == 1:
        return constrain(data[0])
    if bins > len(data):
        bins = len(data)
    data.sort()
    var = running_variance(data)
    rvar = np.flipud(running_variance(np.flipud(data))) 
    thresholds = data[1:len(data):len(data)//bins]
    score_low = (var[0:len(data)-1:len(data)//bins] * 
                 np.arange(0,len(data)-1,len(data)//bins))
    score_high = (rvar[1:len(data):len(data)//bins] *
                  (len(data) - np.arange(1,len(data),len(data)//bins)))
    scores = score_low + score_high
    if len(scores) == 0:
        return constrain(thresholds[0])
    index = np.argwhere(scores == scores.min()).flatten()
    if len(index)==0:
        return constrain(thresholds[0])
    #
    # Take the average of the thresholds to either side of
    # the chosen value to get an intermediate in cases where there is
    # a steep step between the background and foreground
    index = index[0]
    if index == 0:
        index_low = 0
    else:
        index_low = index-1
    if index == len(thresholds)-1:
        index_high = len(thresholds)-1
    else:
        index_high = index+1 
    return constrain((thresholds[index_low]+thresholds[index_high]) / 2)