import random


def natural(values, classes=5, maxsize=1000, samples=3):
    """
    Jenks Optimal (Natural Breaks) algorithm implemented in Python.
    The original Python code comes from here:
    http://danieljlewis.org/2010/06/07/jenks-natural-breaks-algorithm-in-python/
    and is based on a JAVA and Fortran code available here:
    https://stat.ethz.ch/pipermail/r-sig-geo/2006-March/000811.html
    
    Returns class breaks such that classes are internally homogeneous while 
    assuring heterogeneity among classes.

    For very large datasets (larger than maxsize), will calculate only on
    subsample to avoid exponential runtimes. Calculated multiple times (samples)
    and takes the average break values for better consistency. Lower and higher
    bounds are kept intact. 
    """

    #values = sorted(values) # maybe not needed as is already done main.py

    # if too few values, just return breakpoints for each unique value, ignoring classes
    if len(values) <= classes:
        return list(values) + [values[-1]]

    def getbreaks(values, classes):
        # the original algorithm by Carson Farmer
        mat1 = []
        for i in range(0,len(values)+1):
            temp = []
            for j in range(0,classes+1):
                temp.append(0)
            mat1.append(temp)
        mat2 = []
        for i in range(0,len(values)+1):
            temp = []
            for j in range(0,classes+1):
                temp.append(0)
            mat2.append(temp)
        for i in range(1,classes+1):
            mat1[1][i] = 1
            mat2[1][i] = 0
            for j in range(2,len(values)+1):
                mat2[j][i] = float('inf')
        v = 0.0
        for l in range(2,len(values)+1):
            s1 = 0.0
            s2 = 0.0
            w = 0.0
            for m in range(1,l+1):
                i3 = l - m + 1
                val = float(values[i3-1])
                s2 += val * val
                s1 += val
                w += 1
                v = s2 - (s1 * s1) / w
                i4 = i3 - 1
                if i4 != 0:
                    for j in range(2,classes+1):
                        if mat2[l][j] >= (v + mat2[i4][j - 1]):
                            mat1[l][j] = i3
                            mat2[l][j] = v + mat2[i4][j - 1]
            mat1[l][1] = 1
            mat2[l][1] = v
        k = len(values)
        kclass = []
        for i in range(0,classes+1):
            kclass.append(0)
        kclass[classes] = float(values[len(values) - 1])
        kclass[0] = float(values[0])
        countNum = classes
        while countNum >= 2:
            id = int((mat1[k][countNum]) - 2)
            kclass[countNum - 1] = values[id]
            k = int((mat1[k][countNum] - 1))
            countNum -= 1
        return kclass

    # Automatic sub sampling for large datasets
    # The idea of using random sampling for large datasets was in the original code. 
    # However, since these samples tend to produce different results,
    # ...to produce more stable results we might as well calculate the
    # ...breaks several times and using the sample means for the final break values.
    
    if len(values) > maxsize:
        allrandomsamples = []
        for _ in range(samples):
            randomsample = sorted(random.sample(values, maxsize))
            
            # include lower and higher bounds to ensure the whole range is considered
            randomsample[0] = values[0] 
            randomsample[-1] = values[-1]
            
            # get sample break
            tempbreaks = getbreaks(randomsample, classes)
            allrandomsamples.append(tempbreaks)
            
        # get average of all sampled break values
        jenksbreaks = [sum(allbreakvalues)/float(len(allbreakvalues))
                       for allbreakvalues in zip(*allrandomsamples)]
        
    else:
        jenksbreaks = getbreaks(values, classes)

    return jenksbreaks