Jenk's Natural Breaks for Classification

This script, with sample data, can be run to demonstrate how a set of values can be subdivided into categories, based upon the natural breaks (user specifies the number) within the data set.

'''GoodVarFit.py
Goodness of variance fit via Jenks Natural Breaks
Source:    http://danieljlewis.org/files/2010/06/Jenks.pdf
Original:  https://stat.ethz.ch/pipermail/r-sig-geo/2006-March/000811.html
Modifications: Dan Patterson July 2013
'''
def getJenksBreaks( dataList, numClass ):
  dataList.sort()
  mat1 = []
  for i in range(0,len(dataList)+1):
    temp = []
    for j in range(0,numClass+1):
      temp.append(0)
    mat1.append(temp)
  mat2 = []
  for i in range(0,len(dataList)+1):
    temp = []
    for j in range(0,numClass+1):
      temp.append(0)
    mat2.append(temp)
  for i in range(1,numClass+1):
    mat1[1][i] = 1
    mat2[1][i] = 0
    for j in range(2,len(dataList)+1):
      mat2[j][i] = float('inf')
  v = 0.0
  for l in range(2,len(dataList)+1):
    s1 = 0.0
    s2 = 0.0
    w = 0.0
    for m in range(1,l+1):
      i3 = l - m + 1
      val = float(dataList[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,numClass+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(dataList)
  kclass = []
  for i in range(0,numClass+1):
    kclass.append(0)
  kclass[numClass] = float(dataList[len(dataList) - 1])
  countNum = numClass
  while countNum >= 2:
    #print "rank = " + str(mat1[k][countNum])
    id = int((mat1[k][countNum]) - 2)
    #print "val = " + str(dataList[id])
    kclass[countNum - 1] = dataList[id]
    k = int((mat1[k][countNum] - 1))
    countNum -= 1
  return kclass

def getGVF( dataList, numClass ):
  '''
  The Goodness of Variance Fit (GVF) is found by taking the difference between the squared deviations
  from the array mean (SDAM) and the squared deviations from the class means (SDCM), and dividing by the SDAM
  Returns:
    the class start, end, average and values from which they were determined
  '''
  breaks = getJenksBreaks(dataList, numClass)
  dataList.sort()
  listMean = sum(dataList)/float(len(dataList))
  #print listMean
  SDAM = 0.0
  for i in range(0,len(dataList)):
    sqDev = (dataList[i] - listMean)**2
    SDAM += sqDev
  SDCM = 0.0
  breaksData = []
  for i in range(0,numClass):
    if breaks[i] == 0:
      classStart = 0
    else:
      classStart = dataList.index(breaks[i])
      classStart += 1
    classEnd = dataList.index(breaks[i+1])
    classList = dataList[classStart:classEnd+1]
    classMean = sum(classList)/float(len(classList))
    breaksData.append( [classStart, classEnd, classMean, classList] )
    preSDCM = 0.0
    for j in range(0,len(classList)):
      sqDev2 = (classList[j] - classMean)**2
      preSDCM += sqDev2
    SDCM += preSDCM
    goodVarFit = (SDAM - SDCM)/SDAM
  return [breaksData, goodVarFit]
if __name__ == "__main__":
  dataList = [0,1,2,3, 10,11,12,13, 20,21,22,23, 30,31,32,33, 40,41,42,43,
              50,51,52,53,60,61,62,63,70,71,72,73,80,81,82,83,90,91,92,93,
              100,101,102,103,110,111,112,113]
  plotPoints = []
  for numClass in range(2,12):
    breaksData, goodVarFit = getGVF( dataList, numClass )
    breaks = len(breaksData)
    for i in range(0,breaks):
      print "\nBreak IDs (start, end): ", breaksData[i][0], breaksData[i][1]
      print "Class mean:   ", breaksData[i][2]
      print "Class values: ", breaksData[i][3]
    print "Goodness of Var Fit: ", goodVarFit
    plotPoints.append([numClass, goodVarFit])
  print "Final plot points: ", plotPoints