birch-algo.py

from datetime import datetime
from math import sqrt
from operator import add
import csv
import sys
import random

# DONE: Added read in of data set
#       Tree building successfully

# NOT DONE: Focus on adding Phase 3 of BIRCH, e.g. agglomerative hierarchical clustering. This implementation currently just builds the tree.

# FEATURES
# - This is a quick test implementation of the BIRCH algorithm (http://www.cs.sfu.ca/cc/459/han/papers/zhang96.pdf)
# - There's no guarantee that this implements the algorithm correctly.
# - It was never meant to be used in production environments.
# - We used it to study how the algorithm reacts to different (smallish) datasets and how to tweak its parameters.
# - It does not rebuild the tree if it hits a memory boundary
# - Because of that it works with a fixed T, there is no heuristic for T_(i+1)
# - Hence, it does not consider outliers because there's no tree rebuilding
# - Overall, it actually is not the whole algorithm, skips phase 2 to 4!

# Tracking statistics
splitcount = 0
nodecount = 0
entrycount = 0
leafcount = 0


class BaseNode(object):
    def __init__(self):
        global nodecount
        nodecount += 1

        self.n = 0
        self.ls = Vector()
        self.squared = 0

        self.parent = None

    @property
    def is_root(self):
        return not self.parent

    @property
    def level(self):
        count = 0
        r = self.parent
        while r:
            r = r.parent
            count += 1
        return count

    def indent(self):
        indent = '\n'
        r = self.parent
        while r:
            r = r.parent
            indent += '       '
        return indent

    @classmethod
    def closest(cls, node, list, force=False):
        "Returns the closest match of node to list"
        min_dist = 0
        min_item = None

        for item in list:

            dist = item.distance(node)
            # if not force:
            #    print "d:",dist,len(node),len(list)

            if not min_item:
                min_item = item
                min_dist = dist
            elif dist < min_dist:
                min_dist = dist
                min_item = item

        return min_item

    def d0(self, other):
        res = 0.0
        for key, val in self.ls.items():
            if key in other.ls:
                res += (val / self.n - other.ls[key] / other.n) ** 2
            else:
                res += (val / self.n) ** 2
        for key, val in other.ls.items():
            if key not in self.ls:
                res += (val / other.n) ** 2
        return res

    def d2(self, other):
        # print "\n\nself%s %s\nn: %i other n: %i" % (hash(self), self, self.n,other.n)
        return (other.n * self.squared + self.n * other.squared - 2 * (self.ls % other.ls)) / (self.n * other.n)

    def d4(self, other):
        dot1, dot2, dot3 = 0.0, 0.0, 0.0
        for val in self.ls.values():
            dot1 += (val / self.n) ** 2
        for val in other.ls.values():
            dot2 += (val / other.n) ** 2

        for key, val in self.ls.items():
            if key in other.ls:
                dot3 += ((val + other.ls[key]) / (self.n + other.n)) ** 2
            else:
                dot3 += (val / (self.n + other.n)) ** 2
        for key, val in other.ls.items():
            if key not in self.ls:
                dot3 += (val / (self.n + other.n)) ** 2

        return self.n * dot1 + other.n * dot2 - (self.n + other.n) * dot3

    distance = d2

    @classmethod
    def farthest_pair(cls, list):

        max_dist = None
        max_pair = None
        for e1 in list:
            for e2 in list:
                if e1 == e2: continue

                dist = e1.distance(e2)
                if not max_pair:
                    max_pair = (e1, e2,)
                    max_dist = dist
                elif dist > max_dist:
                    max_pair = (e1, e2,)
                    max_dist = dist
        return max_pair

    @classmethod
    def closest_pair(cls, list):
        min_dist = None
        min_pair = None
        for e1 in list:
            for e2 in list:
                if e1 == e2: continue

                dist = e1.distance(e2)
                if not min_pair:
                    min_pair = (e1, e2,)
                    min_dist = dist
                elif dist < min_dist:
                    min_pair = (e1, e2,)
                    min_dist = dist
        return min_dist

    def reset_cf(self):
        self.n = 0
        self.ls = Vector()
        self.squared = 0

    def update_cf(self, data):
        self.n += data.n
        self.ls += data.ls
        self.squared += data.squared

    @classmethod
    def calculate_height(self, list):
        if not list:
            return 1
        else:
            cum = 0
            for x in list:
                cum += x.height
            return cum

    @classmethod
    def calculate_depth(self, list):
        if not list:
            return 0
        else:
            return max([(lambda x: x.depth)(x) for x in list])


class Node(BaseNode):
    # has children which are nodes or leafs

    def __init__(self, *args, **kwargs):
        self.children = []
        super(Node, self).__init__(*args, **kwargs)
        global nodecount
        nodecount += 1

    @property
    def childnodes(self):
        return self.children

    @property
    def height(self):
        return self.calculate_height(self.children)

    @property
    def depth(self):
        return self.calculate_depth(self.children)

    def __str__(self):
        return '%sNODE %i (%i)->' % (self.indent(), hash(self), len(self.children)) + ' '.join(
            [str(c) for c in self.children])

    def trickle(self, vector):
        "Gets a vector and hands it down to the closest child, checks for split afterwards"

        # refresh CF vector
        self.update_cf(vector)

        closest = self.closest(vector, self.children)

        if closest:
            closest.trickle(vector)
        else:
            l = Leaf()
            l.trickle(vector)
            self.add_node(l)

    def add_node(self, node, update=False):
        self.children.append(node)
        node.parent = self

        if update:
            self.update_cf(node)

        if len(self.children) > B:
            self.split_node()

    def split_node(self):
        global splitcount
        splitcount += 1

        c1, c2 = self.farthest_pair(self.children)

        # save the old list
        self.reset_cf()
        old_children = self.children
        self.children = []

        old_children.remove(c1)
        old_children.remove(c2)

        # two new leafs
        if self.is_root:
            node1 = Node()
        else:
            node1 = self
        node2 = Node()

        # add the farthest children to the new nodes
        node1.add_node(c1, True)
        node2.add_node(c2, True)

        while old_children:
            c = old_children.pop()
            if node1.distance(c) > node2.distance(c):
                node2.add_node(c, True)
            else:
                node1.add_node(c, True)

        # try to push down nodes if it only has one child...
        if len(node1.children) == 1:
            node1 = node1.children[0]
        if len(node2.children) == 1:
            node2 = node2.children[0]

        # create a new leaf and append it to our parent
        if self.is_root:
            self.add_node(node1, True)
            self.add_node(node2, True)
            # try to re-merge node 1 or 2
        else:
            self.parent.add_node(node2)


class Leaf(BaseNode):
    # has entries

    def __init__(self, *args, **kwargs):
        self.entries = []
        super(Leaf, self).__init__(*args, **kwargs)
        global leafcount
        leafcount += 1

    @property
    def childnodes(self):
        return self.entries

    @property
    def height(self):
        return self.calculate_height(self.entries)

    @property
    def depth(self):
        return self.calculate_depth(self.entries)

    def __str__(self):
        return '%sLEAF %i (%i)->' % (self.indent(), hash(self), len(self.entries)) + ' '.join(
            [str(c) for c in self.entries])

    def trickle(self, vector):
        "Gets a vector and stores it in the closest entry, checks for split afterwards"

        closest = self.closest(vector, self.entries)

        if closest:
            closest.store_vector(vector)
            self.update_cf(vector)
        else:
            e = Entry()
            e.store_vector(vector)
            self.add_entry(e)

    @classmethod
    def closest(cls, node, list, force=False):
        "Returns the closest match of node to list"
        min_dist = 0
        min_item = None

        if not list:
            return

        for item in list:
            dist = item.distance(node)
            # if not force:
            #    print "d:",dist,len(node),len(list)

            if not min_item:
                min_item = item
                min_dist = dist
            elif dist < min_dist:
                min_dist = dist
                min_item = item

        # try to insert into min_item and check T afterwards
        if min_item.test_radius(node) > T:
            return

        return min_item

    def add_entry(self, entry):
        self.entries.append(entry)
        entry.parent = self

        self.update_cf(entry)
        if len(self.entries) > B:
            self.split_leaf()

    def split_leaf(self):
        global splitcount
        splitcount += 1

        e1, e2 = self.farthest_pair(self.entries)

        # save the old list
        self.reset_cf()
        old_entries = self.entries
        self.entries = []

        old_entries.remove(e1)
        old_entries.remove(e2)

        # two new leafs
        leaf1 = self
        leaf2 = Leaf()
        leaf1.add_entry(e1)
        leaf2.add_entry(e2)

        while old_entries:
            e = old_entries.pop()
            if leaf1.distance(e) > leaf2.distance(e):
                leaf2.add_entry(e)
            else:
                leaf1.add_entry(e)

        # create a new leaf and append it to our parent
        self.parent.add_node(leaf2)


class Entry(BaseNode):
    # has vectors
    is_entry = True

    def __init__(self, *args, **kwargs):
        self.vectors = []
        super(Entry, self).__init__(*args, **kwargs)
        global entrycount
        entrycount += 1
        self.radius = 0.0

    def testvolume(self, vector=None):
        if not self.vectors:
            return 0
        vecs = self.vectors[:]

        if vector:
            vecs.append(vector)

        n = len(vecs)

        dist = 0
        for v1 in vecs:
            for v2 in vecs:
                dist += v1.distance(v2)

        vol = (dist / (n * (n - 1))) ** 0.5

        # print "vol:",vol
        return vol

    volume = property(testvolume)

    def test_radius(self, vector):
        if self.n == 0:
            return 0

        new_n = self.n + 1
        new_ls = self.ls + vector.ls
        new_squared = self.squared + vector.squared

        testrad = self.radius + ((new_ls / new_n).distance(vector))
        # print "testrad:",testrad
        return testrad

    @property
    def volume(self):
        return self.radius

    @property
    def refdist(self):
        dist = 0
        for v1 in self.vectors:
            for v2 in self.vectors:
                dist += v1.distance(v2)
        return dist

    @property
    def height(self):
        return len(self.vectors)

    @property
    def depth(self):
        return 1

    def store_vector(self, vector):
        "Stores a point in its list, this is the end!"
        self.vectors.append(vector)
        self.update_cf(vector)
        self.radius += (self.ls / self.n).distance(vector)


class Vector(dict):
    n = 1

    def __init__(self, *args, **kwargs):
        self._squared = None
        self.ls = self
        super(Vector, self).__init__(*args, **kwargs)

    def __hash__(self):
        assert self.item, "Cannot hash a value without its item"
        return hash(self.item)

    def __add__(self, other):
        output = {}
        for key, val in self.items():
            if key in other:
                output[key] = val + other[key]
            else:
                output[key] = val

        for key, val in other.items():
            if key not in self:
                output[key] = val

        return Vector(output)

    def __sub__(self, other):
        output = {}
        for key, val in self.items():
            if key in other:
                output[key] = val - other[key]
            else:
                output[key] = val

        for key, val in other.items():
            if key not in self:
                output[key] = -val

        return Vector(output)

    def __div__(self, scalar):
        # print "div:",scalar
        output = {}
        # python divisions..
        scalar = float(scalar)
        for key, val in self.items():
            output[key] = val / scalar
        return Vector(output)

    def __mod__(self, other):
        "Dot Product"
        result = 0.0
        for key, val in self.items():
            if key in other:
                result += val * other[key]
        return result

    def __pow__(self, p):
        output = {}
        for key, val in self.items():
            output[key] = val ** p
        return Vector(output)

    def sqrt(self):
        output = {}
        for key, val in self.items():
            output[key] = sqrt(val)
        return Vector(output)

    def length(self):
        length = 0
        for val in self.values():
            length += val ** 2
        return sqrt(length)

    @property
    def squared(self):
        if not self._squared:
            self._squared = self % self
        return self._squared

    def distance(self, other):
        distance = 0
        num = 0
        matching = []
        for key, val in self.items():
            num += 1
            if key in other:
                matching.append(key)
                distance += (val - other[key]) ** 2
            else:
                distance += val ** 2

        for key, val in other.items():
            if key not in self:
                num += 1
                distance += val ** 2

        return sqrt(distance)


if __name__ == '__main__':

    ### SETTINGS
    B = 10  # limit the amount of children a node can have
    # (the paper has another one especially for leaves. But we'll leave it at that!)

    T = 5000  # maximum size (threshold) of a cluster before it has to be split

    vectors = []
    dimensions = 3
    fillpercentage = 100

    input_file = "../dataset/noaa-hail-cleaned-index.csv"

    with open(input_file, 'rb') as f:
        reader = csv.reader(f)
        totalpoints = 0
        try:
            for row in reader:
                totalpoints += 1

                if totalpoints < 150000:

                    v = Vector()
                    v[0] = int(row[0])  # id of data point
                    v[1] = float(row[16])  # latitude of data point
                    v[2] = float(row[17])  # longitude of data point
                    vectors.append(v)  # add each record

                elif totalpoints > 150000:
                    totalpoints -= 1

        except csv.Error as e:
            sys.exit('file %s, line %d: %s' % (filename, reader.line_num, e))

        print "creating vectors.. (%i-dimensional, %i%% filled, %i vectors, branch %i, threshold %i)" % (
        dimensions, fillpercentage, totalpoints, B, T)

    start = datetime.now()  # Start measuring the computing time
    print "starting clustering..."

    # Create root of CF tree
    root = Node()

    # Insert all vectors
    for v in vectors:
        root.trickle(v)  # insert this vector to the tree's root

    # print root.children

    time = datetime.now() - start  # Stop measuring the computing time
    print "took %s (%.2fms per point)" % (time, (time.seconds * 1000 + time.microseconds / 1000.0) / totalpoints)
    print "splitcount: %s, nodecount: %s, entrycount: %s, leafcount: %s" % (
    splitcount, nodecount, entrycount, leafcount)