clusters.py for Collective Intelligence Ch3 Ex2

from math import sqrt
from PIL import Image, ImageDraw
import random as random

def pearson(v1, v2):
    # Simple sums
    sum1 = sum(v1)
    sum2 = sum(v2)

    # Sums of the squares
    sum1Sq = sum([pow(v, 2) for v in v1])
    sum2Sq = sum([pow(v, 2) for v in v2])

    # Sum of the products
    pSum = sum([v1[i]*v2[i] for i in range(len(v1))])

    # Calculatr r (Pearson score)
    num = pSum-(sum1*sum2/len(v1))
    den = sqrt((sum1Sq-pow(sum1, 2)/len(v1))*(sum2Sq-pow(sum2, 2)/len(v1)))
    if den == 0: return 0

    return 1.0-num/den

def tanimoto(v1, v2):
    c1, c2, shr = 0, 0, 0
    for i in range(len(v1)):
        if v1[i] != 0: c1 += 1   # in v1
        if v2[i] != 0: c2 += 1   # in v2
        if v1[i] != 0 and v2[i] != 0: shr += 1   # in both

    return 1.0 -(float(shr)/(c1+c2-shr))

def readfile(filename):
    lines = [line for line in file(filename)]

    # First line is the column titles
    colnames = lines[0].strip().split('\t')[1:]
    rownames = []
    data = []
    for line in lines[1:]:
        p = line.strip().split('\t')
        # First column in each row is the rowname
        rownames.append(p[0])
        # The data for this row is the remainder of the row
        data.append([float(x) for x in p[1:]])
    return rownames, colnames, data

class bicluster:
    def __init__(self,vec,left=None,right=None,distance=0.0,id=None):
        self.left = left
        self.right = right
        self.vec = vec
        self.id = id
        self.distance = distance

def hcluster(rows, distance=pearson):
    distances = {}
    currentclustid = -1

    # Clusters are initially just the rows
    clust = [bicluster(rows[i], id=i) for i in range(len(rows))]

    while len(clust) > 1:
        lowestpair = (0, 1)
        closest = distance(clust[0].vec, clust[1].vec)

        # loop through every pair looking for the smallest distance
        for i in range(len(clust)):
            for j in range(i+1, len(clust)):
                # distances is the cache of distance calculations
                if (clust[i].id, clust[j].id) not in distances:
                    distances[(clust[i].id, clust[j].id)] = distance(clust[i].vec, clust[j].vec)

                d = distances[(clust[i].id, clust[j].id)]

                if d < closest:
                    closest = d
                    lowestpair = (i, j)
        # calculate the average of the two clusters
        mergevec = [(clust[lowestpair[0]].vec[i]+clust[lowestpair[1]].vec[i])/2.0 for i in range(len(clust[0].vec))]

        # create the new cluster
        newcluster = bicluster(mergevec, left=clust[lowestpair[0]], right=clust[lowestpair[1]], distance=closest, id=currentclustid)

        # cluster ids that weren't in the original set are negative
        currentclustid -= 1
        del clust[lowestpair[1]]
        del clust[lowestpair[0]]
        clust.append(newcluster)
    return clust[0]

def printclust(clust, labels=None, n=0):
    # indent to make a hierarchy layout
    for i in range(n): print ' ',
    if clust.id < 0:
        # negative id means that this is branch
        print '-'
    else:
        # positive id means that this is an endpoint
        if labels == None: print clust.id
        else: print labels[clust.id]

    # now print the right and left branches
    if clust.left != None: printclust(clust.left, labels=labels, n=n+1)
    if clust.right != None: printclust(clust.right, labels=labels, n=n+1)

def getheight(clust):
    # Is this an endpoint? Then the height is just 1
    if clust.left == None and clust.right == None: return 1

    # Otherwise the height is the sum of the heights of each branch
    return getheight(clust.left)+getheight(clust.right)

def getdepth(clust):
    # The distance of an endpoint is 0.0
    if clust.left == None and clust.right == None: return 0

    # The distance of a branch is the greater of its two sides plus its own distance
    return max(getdepth(clust.left),getdepth(clust.right))+clust.distance

def drawdendrogram(clust, labels, jpeg='clusters.jpg'):
    # height and width
    h = getheight(clust)*20
    w = 1200
    depth = getdepth(clust)

    # width is fixed, so scale distances accordingly
    scaling = float(w-150)/depth

    # Create a new image with a white background
    img = Image.new('RGB',(w, h),(255,255,255))
    draw = ImageDraw.Draw(img)

    draw.line((0,h/2,10,h/2),fill=(255,0,0))

    # Draw the first node
    drawnode(draw,clust,10,(h/2),scaling,labels)
    img.save(jpeg,'JPEG')

def drawnode(draw, clust, x, y, scaling, labels):
    if clust.id<0:
        h1 = getheight(clust.left)*20
        h2 = getheight(clust.right)*20
        top = y-(h1+h2)/2
        bottom = y+(h1+h2)/2
        # Line length
        ll = clust.distance*scaling
        # Vertical line from this cluster to children
        draw.line((x,top+h1/2, x, bottom-h2/2), fill=(255,0,0))

        # Horizontal line to left item
        draw.line((x,top+h1/2, x+ll, top+h1/2), fill=(255,0,0))

        # Horizontal line to right item
        draw.line((x,bottom-h2/2, x+ll, bottom-h2/2), fill=(255,0,0))

        # Call the function to draw the left and right nodes
        drawnode(draw, clust.left, x+ll, top+h1/2, scaling, labels)
        drawnode(draw, clust.right, x+ll, bottom-h2/2, scaling, labels)
    else:
        # If this is an endpoint, draw the item label
        draw.text((x+5,y-7), labels[clust.id],(0,0,0))

def rotatematrix(data):
    newdata = []
    for i in range(len(data[0])):
        newrow = [data[j][i] for j in range(len(data))]
        newdata.append(newrow)
    return newdata

def kcluster(rows, distance=pearson, k=4):
    # Determine the min and max values for each point
    ranges = [(min([row[i] for row in rows]),max([row[i] for row in rows])) for i in range(len(rows[0]))]

    # Create k randomly placed centroids
    clusters = [[random.random()*(ranges[i][1]-ranges[i][0])+ranges[i][0] for i in range(len(rows[0]))] for j in range(k)]

    lastmatches=None
    for t in range(100):
        print 'Iteration %d' % t
        bestmatches = [[] for i in range(k)]

        # find which centroid is the closest for each row
        for j in range(len(rows)):
            row = rows[j]
            bestmatch = 0
            for i in range(k):
                d = distance(clusters[i],row)
                if d < distance(clusters[bestmatch],row): bestmatch = i
            bestmatches[bestmatch].append(j)

        # If the results are the same as last time, this is complete
        if bestmatches == lastmatches: break
        lastmatches = bestmatches

        # Move the centroids to the average of their members
        for i in range(k):
            avgs = [0.0]*len(rows[0])
            if len(bestmatches[i]) > 0:
                for rowid in bestmatches[i]:
                    for m in range(len(rows[rowid])):
                        avgs[m] += rows[rowid][m]
                for j in range(len(avgs)):
                    avgs[j]/=len(bestmatches[i])
                clusters[i] = avgs

    return bestmatches

def scaledown(data, distance=pearson, rate=0.01):
    n = len(data)

    # The real distances between every pair of items
    realdist = [[distance(data[i], data[j]) for j in range(n)] for i in range(0, n)]

    outersum = 0.0

    # Randomly initialize the starting points of the locations in 2D
    loc = [[random.random(),random.random()] for i in range(n)]
    fakedist = [[0.0 for j in range(n)] for i in range(n)]

    lasterror = None
    for m in range(0, 1000):
        # Find projected distances
        for i in range(n):
            for j in range(n):
                fakedist[i][j]=sqrt(sum([pow(loc[i][x]-loc[j][x],2) for x in range(len(loc[i]))]))

        # Move points
        grad = [[0.0, 0.0] for i in range(n)]

        totalerror = 0
        for k in range(n):
            for j in range(n):
                if j == k: continue
                # The error is percent difference between the distances
                errorterm = (fakedist[j][k]-realdist[j][k])/realdist[j][k]

                # Each point needs to be moved away from or towards the other
                # point in proportion to how much error it has
                grad[k][0] += ((loc[k][0]-loc[j][0])/fakedist[j][k])*errorterm
                grad[k][1] += ((loc[k][1]-loc[j][1])/fakedist[j][k])*errorterm

                # Keep track of the total error
                totalerror += abs(errorterm)
        print totalerror

        # If the answer got worse by moving the points, we are done
        if lasterror and lasterror<totalerror: break
        lasterror = totalerror

        # Move each of the points by the learning rate times the gradient
        for k in range(n):
            loc[k][0] -= rate*grad[k][0]
            loc[k][1] -= rate*grad[k][1]

    return loc

def draw2d(data, labels, jpeg='mds2d.jpg'):
    img=Image.new('RGB',(2000,2000),(255,255,255))
    draw=ImageDraw.Draw(img)
    for i in range(len(data)):
        x=(data[i][0]+0.5)*1000
        y=(data[i][1]+0.5)*1000
        draw.text((x,y),labels[i],(0,0,0))
    img.save(jpeg,'JPEG')

"""
What this program does:
  organises the data in blogdata1.txt into a visual representation of the relationships
  between elents within it. ie Words and the frequency of their usage.
  A potential use could be to parse a text find frequencies of word usage
  and determine from the style of the language used which jargon/slang is being used
  and therefore how to translate it better.
  Or evaluate a persons favorite style and recommend writings in the same/similar styles.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@----Hierarchical clustering----@
@-------------------------------@
@ step   result                 @
@ 0      1,2,5,6,11,14          @
@ 1      [1,2] 5,9,17,19        @
@ 2      [1,2] 5,9 [17,19]      @
@ 3      [[1,2] 5] 9 [17,19]    @
@ 4      [[[1,2] 5] 9] [17,19]  @
@ 5      [[[1,2] 5] 9 [17,19]]  @
@                               @
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@         scaledown() + draw2d()
@------------K Means clustering--------------@         produces a visual representation
@--------------------------------------------@         kclustered words in 2D
@ step   result                              @
@ 0      1,2,5,6,11,14                       @
@ 1      [1 2 (k4)5]         [6 (k7) 11 14]  @
@ 2      [1 2 (k8/3) 5  6]   [(k31/3) 11 14] @
@ 3      [1 2 (k14/4) 5 6]   [11 (k23/2) 14] @
@                                            @
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
where k is a randomly placed centroid

    blogdata1.txt is a dataset
it is in the format rows columns
the first row is a list of commonly used words
the rest of the rows are lists of the number of times each word has been used
the indexes of the times used matches the indexes of the words

readfile()       reads the data in blogdata1.txt and returns it as lists (eg blognames,words)
                                                          and list of lists (eg data)
hcluster()       does sim_pearson on the vectors of the words(ie times they appear)
printclust()     prints out the blogs related by word use onscreen
drawdendrogram() creates a jpeg dendrogram, a view of the word use relationships

rotatematrix()   creates a matrix of word - its usage in all blogs
                       rather than   blog - word(usage in this blog)
                 looks at individual words relationships to other words
                 and shows words used with the same frequency with blogs
                 when the rotated matrix is passed to drawdendrogram()


blogs with similar word usage get linked in the dendrogram



@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Visualizing the workings of hcluster:

###############################################################################
##########clust = [bicluster(rows[i], id=i) for i in range(len(rows))]#########
###############################################################################

>>> frog = [[1,2,3],[4,5,6]]
>>> fruit = [bicluster(frog[i],id=i) for i in range(len(frog))]
>>> fruit[1].id
>>> fruit[1].vec
[4, 5, 6]
>>> print fruit[0].__dict__
{'id': 0, 'distance': 0.0, 'right': None, 'vec': [1, 2, 3], 'left': None}
>>> print fruit[1].__dict__
{'id': 1, 'distance': 0.0, 'right': None, 'vec': [4, 5, 6], 'left': None}
>>> print len(fruit)
2

so cluster is created as a list of instances of the class bicluster

###############################################################################
##############################while len(clust) > 1#############################
###############################################################################
Within this section
len(clust) is reduced to 1
The dictionary distances is created
     distances(id1,id2) = sim_pearson score
     where ids are from clust[instances of bicluster]


###############################################################################
###################for i in range(len(clust)):#################################
#######################for j in range(i+1, len(clust)):########################
###############################################################################
Brute-style
i is an item in clust
j is the rest of the items
if the sim_pearson score is less than closest i,j become the lowestpair

###############################################################################
mergevec = [(clust[lowestpair[0]].vec[i]+clust[lowestpair[1]].vec[i])/2.0 for i in range(len(clust[0].vec))]
###############################################################################
(clust[lowestpair[0]].vec[i]
    lowestpair = (A, B)
    A is an instance of clust
    .vec[i] is a vector from list .vec
    therefore gets the averages of all word instances for A, B and gets average usage

###############################################################################
newcluster = bicluster(mergevec, left=clust[lowestpair[0]], right=clust[lowestpair[1]], distance=closest, id=currentclustid)
###############################################################################

newcluster is an instance of bicluster
   .vec is the list of averages from A(for i...) and B(for j...)
   left is instance of bicluster A
   right is instance of bicluster A
   id starts at -1 and goes down every round by -1

   the two lowestclusters having been joined together in the new cluster
   have their originaL instances deleted

   and we return to the beginning of while len(clust) > 1
   the two deleted instances are represented in cluster by
    a single instance with the originals in .left .right
    the new instance has average vecs

   if on the next round of while the lowest matching pair contains newcluster
   then newcluster will be added(together with .right .left) as .left or .right
   of this new newcluster

   this process will continue until all bicluster instances are matrixed together

ATALHEA

Example usage:

import clusters as clusters
blognames, words, data=clusters.readfile('deliciousdata1.txt')
clust=clusters.hcluster(data)
rdata=clusters.rotatematrix(data)
wordclust = clusters.hcluster(rdata)

clusters.drawdendrogram(wordclust, labels=words, jpeg='deliciousclust.jpg')

clusters.drawdendrogram(clust, labels=blognames, jpeg='deliciousclust.jpg')



kclust = clusters.kcluster(data,k=10)
[blognames[r] for r in kclust[0]]
-------------------------------------------------------------------------------
##########################
# what does scaledown do #
##########################
Return a list containig lists of coordinates
First it calculates the distances between all items
Then using a Kcluster-like method it produces coordinates x and y
    for each item
It does the kclustering process m times less times = less accuracy
                                        more times = more time to calculate
Each item is initially given a random coord and is moved towards similar items
All items are moved m times sorting themselves into a more meaningful cluster


"""