Home » Posts tagged 'python'

# Tag Archives: python

## Understanding scan() in Theano

This is just short post since the meat is in an ipython notebook. I discuss the scan function in Theano, used for looping, iteration over tensors and recurrences. If you are a complete beginner to theano, check out my previous post. It seemed a bit obscure when I first started trying to understand it, so I thought I would share my learning process and hopefully speed up yours!

My sources for learning were these two great resources: this and that.

OK check it out here and let me know what you think.

## Theano, Autoencoders and MNIST

Recently I have finally upgraded my ancient laptop. Other than the ability for me to play the occasional video game, this means that I now have a dedicated Nvidia graphics card, in particular one that supports something called CUDA.

Now that I have the hardware, I thought that it would be fun to check out this python library theano. As I understand it, it has two very interesting aspects.

One is that it does symbolic computation, that is it represents variables as in mathematics, and their manipulations.

A variable in computer science is something like a integer or a string or a list etc…, a piece of information stored somewhere in memory with a name made up by the programmer, and the compiler takes care of associating the human understandable name to the address of the information in memory.

A variable in mathematics is like in , an indeterminate quantity that can be manipulated using arithmetic operations, or using calculus and so on. So what in mathematics we call a variable corresponds roughly to a function in computer science.

The difference between symbolic computation and its variables, and regular functions, is that the programming language is aware of the mathematical nature of a variable, such that when you apply operations to variables to create new variables, it does this in an efficient way. The best way to understand this is a simple example. Let’s suppose we have two functions . Now we can combine these functions by adding them to make a new function .

Now in python say, when you call with an argument, it will pass this argument to both and and then add the result and return it.

If we are doing symbolic computation, and were variables with , then the variable constantly. And so the computer does no work at all when using , because it knows that it must always be zero.

So one important aspect is that symbolic computation can simplify and hence often optimize functions automatically.

Another is that theano can automagically differentiate variables with respect to variables! This might not sound like a big deal, but it is. In machine learning, we often create a complicated function, like in logistic regression or neural networks, and then derive a helpful form for the derivative to do gradient descent to optimize it. In theano, *all you have to do is specify the objective function*, and it will let you do gradient descent automatically.

The second amazing thing about theano is that it compiles its theano functions into C, and so it often runs much much faster than native python or numpy. It can also make it easier to run computations on a gpu which can mean a 100x speedup for the right cases!

Hopefully I’ve made you curious, so click here for an installation guide. On this linux machine it was very simple. I am dual booting windows 8.1 and Ubuntu at the moment, it is very tricky to get the gpu aspect to work on windows, but here is a guide for the brave. To check that the GPU aspect is working, run the test code here.

**Introduction to Theano**

I recommend that you read the tutorials here, and also this excellent ipython worksheet here.

Let’s little play with the basic functionality.

import theano as th from theano import tensor as T #To make a variable, specify the type of variable from the tensor #library. #We can have scalars. x=T.dscalar(name='scalar_x') #Note that the name is just a string that theano will use #to communicate to us about the variable. print x print x.type #We can also have vectors. v = T.dvector(name='vector_v') print v print v.type #And matrices. A = T.dmatrix(name='matrix_A') print A print A.type #We can also make new variables using standard mathematical operations. x_2 = x*x print x_2 print x_2.type #We can also make python functions which return new variables. def power(variable, n): return variable**n x_10 = power(x,10) print x_10 print x_10.type #We can of course do standard linear algebra operations also. Av = T.dot(A,v) print Av print Av.type #Of course when a variable is actually evaluated, you must #ensure that the dimensions are correct, or you will #get an error. #To see the value of a variable for a particular value of the variables #comprising it, we make a theano function. f = th.function([A,v], [Av]) #The syntax is a list of input variables, and then a list of output variables. #The python code takes a little longer to run initially, because thenao #compiles the function into C, but thereafter it will run extremely fast. #Let's try using the function. import numpy as np m=np.ones((2,2)) print m vec = np.asarray([1,2]) print vec print f(m,vec) #Now we can try computing gradients. #First let's make a scalar variable by taking an inner product. w=T.dvector(name='w_vector') vTw = T.dot(v,w) #Now we take the gradient with respect to w. vTw_grad = T.grad(vTw,w) print vTw_grad print vTw_grad.type #Now let's test it. vec1 = np.asarray([1,2]) vec2 = np.asarray([0,0]) #To evaulate a variable given inputs, there is another syntax #in additon to creating a thenano function. print vTw_grad.eval({w:vec1,v:vec2}) #Next we discuss theano's shared variables, #these differ from regular theano variables in that #theano variables only have a value within a theano function #whereas shared theano variables have a value independent of #being called in a function. w=T.shared(name='shared_matrix', value=np.ones((2,3))) print w print w.type print w.get_value() #You can also set a shared variable's value. w.set_value(np.zeros((2,3))) print w.get_value() #You can also have theano functions update shared variables, #we illustrate this with a silly example that updates a matrix #like in gradient descent. x=T.dvector('x') wx=T.dot(w,x) cost = wx.mean() cost_grad = T.grad(cost, w) f=th.function(inputs=[x], outputs=[wx, cost], updates=[(w, w-0.1*cost_grad)]) #Notice the syntax of updates argument, should be a list #of two tuples of the form: (variable_to_be_updated, updated_variable). print f([1,1,1]) print w.get_value()

If you want to tinker with this, get the sourcecode here.

**Autoencoders**

I don’t want to go into too much detail here, but an autoencoder is a feedforward neural network with a single hidden layer with the twist being that its target output is exactly the same as its input. If the number of hidden neurons is less than the number of input/output neurons, then the network will be forced to learn some sort of compression.

Even when the number of hidden neurons is large, there are various tricks in cost functions you can use to ensure that the network does something more interesting than learning the identity function.

The hope of an autoencoder is that it learns some structure in the data, so that the hidden layer can be thought of as a feature space that one can use for a supervised learning algorithm.

Basically then we have the hidden layer, then , where are bias vectors, are matrices, and are element-wise non-linearities like the logistic function.

A simplification that is often made is to enforce . This gives the model less free parameters to learn, and so is a form of regularization. It also uses less memory. This is referred to as an autoencoder with tied weights.

**Autoencoders with Theano**

Thanks to automatic differentiation, backprogation is effortless with theano, and thanks to gpu computation, much faster also! I recommend having a look at my python implementation of backpropagation just so you can see the effort saved. This implementation benefited greatly from this tutorial.

import numpy as np import theano as th from theano import tensor as T from numpy import random as rng class AutoEncoder(object): def __init__(self, X, hidden_size, activation_function, output_function): #X is the data, an m x n numpy matrix #where rows correspond to datapoints #and columns correspond to features. assert type(X) is np.ndarray assert len(X.shape)==2 self.X=X self.X=th.shared(name='X', value=np.asarray(self.X, dtype=th.config.floatX),borrow=True) #The config.floatX and borrow=True stuff is to get this to run #fast on the gpu. I recommend just doing this without thinking about #it until you understand the code as a whole, then learning more #about gpus and theano. self.n = X.shape[1] self.m = X.shape[0] #Hidden_size is the number of neurons in the hidden layer, an int. assert type(hidden_size) is int assert hidden_size > 0 self.hidden_size=hidden_size initial_W = np.asarray(rng.uniform( low=-4 * np.sqrt(6. / (self.hidden_size + self.n)), high=4 * np.sqrt(6. / (self.hidden_size + self.n)), size=(self.n, self.hidden_size)), dtype=th.config.floatX) self.W = th.shared(value=initial_W, name='W', borrow=True) self.b1 = th.shared(name='b1', value=np.zeros(shape=(self.hidden_size,), dtype=th.config.floatX),borrow=True) self.b2 = th.shared(name='b2', value=np.zeros(shape=(self.n,), dtype=th.config.floatX),borrow=True) self.activation_function=activation_function self.output_function=output_function def train(self, n_epochs=100, mini_batch_size=1, learning_rate=0.1): index = T.lscalar() x=T.matrix('x') params = [self.W, self.b1, self.b2] hidden = self.activation_function(T.dot(x, self.W)+self.b1) output = T.dot(hidden,T.transpose(self.W))+self.b2 output = self.output_function(output) #Use cross-entropy loss. L = -T.sum(x*T.log(output) + (1-x)*T.log(1-output), axis=1) cost=L.mean() updates=[] #Return gradient with respect to W, b1, b2. gparams = T.grad(cost,params) #Create a list of 2 tuples for updates. for param, gparam in zip(params, gparams): updates.append((param, param-learning_rate*gparam)) #Train given a mini-batch of the data. train = th.function(inputs=[index], outputs=[cost], updates=updates, givens={x:self.X[index:index+mini_batch_size,:]}) import time start_time = time.clock() for epoch in xrange(n_epochs): print "Epoch:",epoch for row in xrange(0,self.m, mini_batch_size): train(row) end_time = time.clock() print "Average time per epoch=", (end_time-start_time)/n_epochs def get_hidden(self,data): x=T.dmatrix('x') hidden = self.activation_function(T.dot(x,self.W)+self.b1) transformed_data = th.function(inputs=[x], outputs=[hidden]) return transformed_data(data) def get_weights(self): return [self.W.get_value(), self.b1.get_value(), self.b2.get_value()]

And that is it! Now we should test this puppy out, so let’s find some data.

**MNIST**

The MNIST data is famous in machine learning circles, it consists of single handwritten digits. Nearly every paper written on neural networks and so on tests their contribution on the data. The task is to classify the digits, but we will just test our autoencoder. First we will steal a very nice helper function to load the data set, this will download it for you if you don’t have it already. You can find it here.

import cPickle import gzip import os

def load_data(dataset): ''' Loads the dataset :type dataset: string :param dataset: the path to the dataset (here MNIST) ''' ############# # LOAD DATA # ############# # Download the MNIST dataset if it is not present data_dir, data_file = os.path.split(dataset) if data_dir == "" and not os.path.isfile(dataset): # Check if dataset is in the data directory. new_path = os.path.join(os.path.split(__file__)[0], "..", "data", dataset) if os.path.isfile(new_path) or data_file == 'mnist.pkl.gz': dataset = new_path if (not os.path.isfile(dataset)) and data_file == 'mnist.pkl.gz': import urllib origin = 'http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz' print 'Downloading data from %s' % origin urllib.urlretrieve(origin, dataset) print '... loading data' # Load the dataset f = gzip.open(dataset, 'rb') train_set, valid_set, test_set = cPickle.load(f) f.close() #train_set, valid_set, test_set format: tuple(input, target) #input is an numpy.ndarray of 2 dimensions (a matrix) #which row's correspond to an example. target is a #numpy.ndarray of 1 dimensions (vector)) that have the same length as #the number of rows in the input. It should give the target #target to the example with the same index in the input. return (train_set, valid_set, test_set) path="/home/harri/Dropbox/Work/Blogs/triangleinequality/Theano/data/mnist.pkl.gz" data=load_data(path)

We also want a plotting function to visualize the features learned, I won’t explain this here, but if you google “visualizing features of neural network” you can learn plenty about it.

def plot_first_k_numbers(X,k): from matplotlib import mpl,pyplot m=X.shape[0] k=min(m,k) j = int(round(k / 10.0)) fig, ax = pyplot.subplots(j,10) for i in range(k): w=X[i,:] w=w.reshape(28,28) ax[i/10, i%10].imshow(w,cmap=pyplot.cm.gist_yarg, interpolation='nearest', aspect='equal') ax[i/10, i%10].axis('off') pyplot.tick_params(\ axis='x', # changes apply to the x-axis which='both', # both major and minor ticks are affected bottom='off', # ticks along the bottom edge are off top='off', # ticks along the top edge are off labelbottom='off') pyplot.tick_params(\ axis='y', # changes apply to the x-axis which='both', # both major and minor ticks are affected left='off', right='off', # ticks along the top edge are off labelleft='off') fig.show()

Now we put it together.

def m_test(data): X=data[0][0] activation_function = T.nnet.sigmoid output_function=activation_function A = AutoEncoder(X, 500, activation_function, output_function) A.train(20,20) W=np.transpose(A.get_weights()[0]) plot_first_k_numbers(W, 100, False) plot_first_k_numbers(W,100, True) m_test(data)

Running this, with 500 hidden units for 20 epochs with minibatch size of 20, the following features were learned:

Nothing too interesting yet, in order to get more useful features we could change the cost function to encourage sparsity, or introduce noise into the inputs, or use a different activation function. Hopefully this has been helpful for those just starting out with theano. Here is the source code.

## How to Balance your Binary Search Trees – AVL Trees

Last time we introduced the binary search tree (BST) and saw that they could do inserts and deletions in time where is the height of the tree. We also saw that it could find the successor or predecessor to a node in the same time, and hence that it could sort a list in time, where is the length of the list.

Ideally we would have , and this is often the case ‘on average’, but the worst case is , which will occur if we start with a sorted list! So today we are going to discuss how to have our tree try to balance itself automagically so as to keep its height .

**AVL-Trees**

Recall that the defining characteristic of a BST is the so-called ‘search property’: for any subtree rooted at a node, the stuff on the left is smaller, and the stuff on the right is greater.

We now introduce another property or invariant that we will insist upon and call it the ‘plus or minus property’ (pm) or ‘balance property’, and it is satisfied by a node if the height of its left subtree differs from the height of its right subtree by no more than one. We call a BST that satisfies pm an ‘AVL tree’ (AVL).

So we will give each node an extra field node.height, and node.balance, where node.balance is equal to the height on the left minus the height on the right. Hence it will be 0 when they are equal, +1 when the left is larger, and -1 when the right is larger.

Before we go any further, you might be wondering whether or not this property forces the tree to have height, and indeed it does, so let’s prove that to ourselves before coercing our trees to have this property.

**Proof That AVL Trees are Balanced**

Define to be the smallest number of nodes you can build a AVL of height with. That is, any AVL of height must have at least nodes.

Now suppose that we had shown that , or in particular that for all . Now take any AVL tree of height and having nodes. Then as desired.

So we just need to show that . The way to do this is by setting up a recurrence relation. Indeed we will see that .

Now using this notice that:

.

And so if then applying this repeatedly yields:

.

On the other hand if then we apply till we get to giving:

.

Combining these gives . Now taking logs base gives

.

And we are done.

** Keeping Track of the AVL Property**

Now that we agree that the AVL or pm property is a thing worth having, we should alter our code for a BST so that the nodes have height and balance fields, and so that the insertion and deletion methods update this accordingly. Once we have done that we will look at how to enforce the balance property.

Now we could use class inheritance and modify the BST class we built last time, but some of the code is a bit kludgy and so I will implement this from scratch, although I won’t go into detail about the stuff that is the same from last time.

We begin with the AVLNode and gift it with a sense of balance.

class AVLNode(object): def __init__(self, key): self.key=key self.right_child=None self.left_child=None self.parent=None self.height=0 self.balance=0 def update_height(self, upwards=True): #If upwards we go up the tree correcting heights and balances, #if not we just correct the given node. if self.left_child is None: #Empty left tree. left_height = 0 else: left_height = self.left_child.height+1 if self.right_child is None: #Empty right tree. right_height = 0 else: right_height = self.right_child.height+1 #Note that the balance can change even when the height does not, #so change it before checking to see if height needs updating. self.balance = left_height-right_height height = max(left_height, right_height) if self.height != height: self.height = height if self.parent is not None: #We only need to go up a level if the height changes. if upwards: self.parent.update_height() def is_left(self): #Handy to find out whether a node is a left or right child or neither. if self.parent is None: return self.parent else: return self is self.parent.left_child

I think the code, and comments, should explain everything if you read my last post. Let’s also pretend that we had a magic function balance(node), which will balance a node and then move up ancestors and do the same, we won’t worry about what exactly it will do, but add it into our insert routine.

class AVLTree(object): def __init__(self): self.root =None def insert(self, key, node=None): #The first call is slightly different. if node is None: #First call, start node at root. node = self.root if node is None: #Empty tree, create root. node = AVLNode(key=key) self.root=node return node else: ret= self.insert(key=key, node=node) self.balance(ret) return ret #Not a first call. if node.key ==key: #No need to insert, key already present. return node elif node.key >key: child = node.left_child if child is None: #Reached the bottom, insert node and update heights. child = AVLNode(key=key) child.parent=node node.left_child = child node.update_height() return child else: return self.insert(key=key, node=child) elif node.key < key: child = node.right_child if child is None: #Reached the bottom, insert node and update heights. child = AVLNode(key=key) child.parent=node node.right_child = child return child else: return self.insert(key=key, node=child) else: print "This shouldn't happen." def find(self, key, node=None): if node is None: #First call. node=self.root if self.root is None: return None else: return self.find(key, self.root) #Now we handle nonfirst calls. elif node.key == key: #Found the node. return node elif key < node.key: if node.left_child is None: #If key not in tree, we return a node that would be its parent. return node else: return self.find(key,node.left_child) else: if node.right_child is None: return node else: return self.find(key, node.right_child) def delete(self, key, node=None): #Delete key from tree. if node is None: #Initial call. node = self.find(key) if (node is None) or (node.key != key): #Empty tree or key not in tree. return if (node.left_child is None) and (node.right_child is not None): #Has one right child. right_child = node.right_child left = node.is_left() if left is not None: parent=node.parent if not left: parent.right_child=right_child else: parent.left_child=right_child right_child.parent =parent self.balance(parent) else: right_child.parent=None self.root = right_child #No need to update heights or rebalance. elif (node.left_child is not None) and (node.right_child is None): #Has one left child. left_child = node.left_child left= node.is_left() if left is not None: parent=node.parent if left: parent.left_child=left_child else: parent.right_child=right_child left_child.parent =parent self.balance(parent) else: left_child.parent=None self.root=left_child elif node.left_child is None: #Has no children. parent = node.parent if parent is None: #Deleting a lone root, set tree to empty. self.root = None else: if parent.left_child is node: parent.left_child =None else: parent.right_child=None self.balance(parent) else: #Node has two childen, swap keys with successor node #and delete successor node. right_most_child = self.find_leftmost(node.right_child) node.key = right_most_child.key self.delete(key=node.key,node=right_most_child) #Note that updating the heights will be handled in the next #call of delete. def find_rightmost(self, node): if node.right_child is None: return node else: return self.find_rightmost(node.right_child) def find_leftmost(self, node): if node.left_child is None: return node else: return self.find_leftmost(node.left_child) def find_next(self, key): node = self.find(key) if (node is None) or (node.key != key): #Key not in tree. return None else: right_child = node.right_child if right_child is not None: node= self.find_leftmost(right_child) else: parent = node.parent while(parent is not None): if node is parent.left_child: break node = parent parent = node.parent node=parent if node is None: #Key is largest in tree. return node else: return node.key def find_prev(self, key): node = self.find(key) if (node is None) or (node.key != key): #Key not in tree. return None else: left_child = node.left_child if left_child is not None: node= self.find_rightmost(left_child) else: parent = node.parent while(parent is not None): if node is parent.right_child: break node = parent parent = node.parent node=parent if node is None: #Key is largest in tree. return node else: return node.key #I also include a new plotting routine to show the balances or keys of the node. def plot(self, balance=False): #Builds a copy of the BST in igraphs for plotting. #Since exporting the adjacency lists loses information about #left and right children, we build it using a queue. import igraph as igraphs G = igraphs.Graph() if self.root is not None: G.add_vertices(1) queue = [[self.root,0]] #Queue has a pointer to the node in our BST, and its index #in the igraphs copy. index=0 not_break=True while(not_break): #At each iteration, we label the head of the queue with its key, #then add any children into the igraphs graph, #and into the queue. node=queue[0][0] #Select front of queue. node_index = queue[0][1] if not balance: G.vs[node_index]['label']=node.key else: G.vs[node_index]['label']=node.balance if index ==0: #Label root green. G.vs[node_index]['color']='green' if node.left_child is not None: G.add_vertices(1) G.add_edges([(node_index, index+1)]) queue.append([node.left_child,index+1]) G.vs[index+1]['color']='red' #Left children are red. index+=1 if node.right_child is not None: G.add_vertices(1) G.add_edges([(node_index, index+1)]) G.vs[index+1]['color']='blue' queue.append([node.right_child, index+1]) index += 1 queue.pop(0) if len(queue)==0: not_break=False layout = G.layout_reingold_tilford(root=0) igraphs.plot(G, layout=layout)

If you want to make sure it works, I have written a small test function.

def test(): lst= [1,4,2,5,1,3,7,11,4.5] B = AVLTree() for item in lst: print "inserting", item B.insert(item) B.plot(True) print "End of inserts" print "Deleting 5" B.plot(True) B.delete(5) print "Deleting 1" B.plot(True) B.delete(1) B.plot(False) print B.root.key ==4 print B.find_next(3) ==4 print B.find_prev(7)==4.5 print B.find_prev(1) is None print B.find_prev(7)==4.5 print B.find_prev(2) is None print B.find_prev(11) == 7

Great now we know the balance of each node, but how to balance the tree?

**Tree Rotations**

Before solving our balancing problem I want to introduce a new tool, an operation called a tree rotation that respects the search property. Here’s what a tree rotation looks like:

So there are two moves that are inverses of one another, a right rotation and a left rotation. Staring at the diagram for a little bit you can also see that these moves preserve the search property. Let’s implement this, and then consider when to use it.

We’ll first create a right rotation, then swap the lefts for rights and vice versa to make a left rotation.

The argument to the function will always be the root, that is the parent of the two nodes, so in our diagram above, the label root is appropriate for the right-rotation.

We then move around the the parent-child pointers. Notice that the relation between the root and , and the relation between pivot and , is unchanged. The changes we have to make are that changes parents, and that root and pivot swap places from parent/child to child/parent.

So the left-child of root, pivot, is changed to the right-child of pivot, , the right-child of pivot is changed to root, the parent of root is changed to pivot and finally the previous parent of root, if any, becomes the new parent of pivot.

def right_rotation(self, root): left=root.is_left() pivot = root.left_child if pivot is None: return root.left_child = pivot.right_child if pivot.right_child is not None: root.left_child.parent = root pivot.right_child = root pivot.parent = root.parent root.parent=pivot if left is None: self.root = pivot elif left: pivot.parent.left_child=pivot else: pivot.parent.right_child=pivot root.update_height(False) pivot.update_height(False)

And analogously:

def left_rotation(self, root): left=root.is_left() pivot = root.right_child if pivot is None: return root.right_child = pivot.left_child if pivot.left_child is not None: root.right_child.parent = root pivot.left_child = root pivot.parent = root.parent root.parent=pivot if left is None: self.root = pivot elif left: pivot.parent.left_child=pivot else: pivot.parent.right_child=pivot root.update_height(False) pivot.update_height(False)

Now we ought to test it so using the same list as last time:

def test_rotation(): lst= [1,4,2,5,1,3,7,11,4.5] print "List is ",lst B = AVLTree() for item in lst: print "inserting", item B.insert(item) node=B.find(4) B.plot() B.left_rotation(node) B.plot() B.right_rotation(node.parent) B.plot() test_rotation()

Producing the following two plots:

Finally we show how to use tree rotations to keep a AVL tree balanced after insertions and deletions.

**Bringing Balance to the Force
**

After an insertion or deletion, heights can change by at most 1. So assuming that the tree was balanced prior to the operation, if a node becomes unbalanced, it will have a balance of .

Let’s look at the case where we have a node with balance with some pictures, the other case is the same modulo lefts and rights.

This is called the left-left case. Alternatively we could have the left-right case:

This takes to a similar case to before:

To implement this is simple now that we have the right tools. We want to function to start with the changed node, either the inserted node in insertion, or the parent of the deleted node in deletion, and move up through the ancestors, updating heights and balancing as it goes.

def balance(self, node): node.update_height(False) if node.balance == 2: if node.left_child.balance != -1: #Left-left case. self.right_rotation(node) if node.parent.parent is not None: #Move up a level. self.balance(node.parent.parent) else: #Left-right case. self.left_rotation(node.left_child) self.balance(node) elif node.balance ==-2: if node.right_child.balance != 1: #Right-right case. self.left_rotation(node) if node.parent.parent is not None: self.balance(node.parent.parent) else: #Right-left case. self.right_rotation(node.right_child) self.balance(node) else: if node.parent is not None: self.balance(node.parent)

We might as well make a sorting routine to wrap the AVL Tree now.

def sort(lst, ascending=True): A = AVLTree() for item in lst: A.insert(item) ret=[] if ascending: node=A.find_leftmost(A.root) if node is not None: key = node.key else: key=node while (key is not None): ret.append(key) key=A.find_next(key) else: node=A.find_rightmost(A.root) if node is not None: key = node.key else: key=node while (key is not None): ret.append(key) key=A.find_prev(key) return ret

And that’s it! Please let me know if you spot any mistakes. You can download the source here.

## Binary Search Trees

So in preparation for starting a PhD in Machine Learning/Data Science at the University of Edinburgh .

I have been watching some of the lectures on algorithms and data structures at MIT here.

One of the data structures discussed is the *binary search tree* (BST), and so in this post I will explain what they are and give a python implementation.

**Basics**

As you might expect, a BST is a way of organizing data in a binary tree. It consists of a collection of *nodes*, each of which has a *key*, corresponding to a piece of data, say an integer. Importantly the keys must be comparable, so that for any two nodes either or .

Each node also has (at most) three pointers: parent, left child and right child.

Most importantly we have the *invariant* of the BST. An invariant is, of course, something that stays the same, a property that is not altered under any of the data structures allowable operations. The invariant or* search property* is that for any given node all of the nodes in the left subtree rooted at have keys less than , and all of the nodes in the right subtree have keys greater than .

Any useful data structure has operations. Some obvious things to want to do are to *insert* and *delete. *Of course in order to do these things we have to take care not to mess up the search property. To insert we begin at the root and at each stage we ask whether the key to be inserted is lesser, greater or equal to the parent, in which case we move right, left or nowhere respectively. It’s easiest just to see an example:

Now I think you get the idea, so let’s start coding this up. First we will make the BST node.

class TreeNode(object): #A tree node has a key, which can be compared with other keys, #and possibly 'pointers' to a parent, left child or right child. def __init__(self, key, parent=None, left_child=None, right_child=None): self.parent=parent self.left_child=left_child self.right_child=right_child self.key=key

Next we will start making the BST. The BST remembers the root, and the pointers do the rest of the work. We begin by implementing a helpful method called ‘find’, which returns a node with a given key, or that could (but doesn’t) have a child with that key.

class BinarySearchTree(object): def __init__(self, root=None): self.root=root def find(self, key): #Returns a node corresponding to a key. #If key is the key of some node, returns that node. #If tree is empty, returns None. #Otherwise returns a leaf that could accept #a node with that key as a child. #This function wraps a call to a recursive function self._find. if self.root is None: return None else: return self._find(key, self.root) def _find(self, key, node): #This is a recursive function that does all the work for #self.find. if key == node.key: return node elif key < node.key: if node.left_child is None: return node else: return self._find(key, node.left_child) else: if node.right_child is None: return node else: return self._find(key, node.right_child)

Using ‘find’ we can then implement ‘insert’.

def insert(self, key): #Inserts a node with given key. #If key already in tree, then returns the node with that key. #Otherwise creates a new node and returns it. #This takes time of order the height of the tree. node = self.find(key) if node is None: self.root=TreeNode(key) return self.root elif key == node.key: return node elif key < node.key: left_child = TreeNode(key, parent=node) node.left_child=left_child return left_child else: right_child = TreeNode(key, parent=node) node.right_child=right_child return right_child

So far so good!

The next thing we want to do is be able to delete keys.

If a key belongs to a leaf node then this is simple we just delete the node.

If the node has one child, then if it has a parent then when we delete the node we must connect the parent to the child of the deleted node.

If the node has two children things are a bit more tricky. If we delete the node then we break the tree. What we need to do instead is replace the node’s key with something else in the tree. What we replace it with is the rightmost node of the subtree rooted at the node’s left child, as we shall see later this is the successor of the node in the BST. The following diagram should explain:

So we swap the node’s key for this new rightmode node’s, then delete the rightmost node. Since it is rightmost it cannot have a right child, and so it is one of our base cases for the deletion operation. If this is still not clear to you, do some examples on paper and you will soon see the idea. Let’s implement this.

def delete(self, key): #Delete key from tree. #If key is not in BST does nothing. #Otherwise it calls a semi-recursive function _delete. node = self.find(key) if (node is None) or (node.key != key): return else: self._delete(node) def _delete(self, node): #If the node has less than two children we can delete the node #directly by removing it and then gluing the tree back together. #If node has two children, it swaps a key lower down in the tree to #replace the deleted node, and deletes the lower down node. if (node.left_child is None) and (node.right_child is not None): #Has one right child. right_child = node.right_child parent = node.parent if parent is not None: parent.right_child=right_child right_child.parent =parent else: right_child.parent=None if node is self.root: self.root=right_child elif (node.left_child is not None) and (node.right_child is None): #Has one left child. left_child = node.left_child parent = node.parent if parent is not None: parent.left_child=left_child left_child.parent =parent else: left_child.parent=None if node is self.root: self.root=left_child elif node.left_child is None: #Has no children. parent = node.parent if parent is None: self.root = None else: if parent.left_child is not None: left = parent.left_child is node else: left = False if left: parent.left_child =None else: parent.right_child=None else: right_most_child = self.find_leftmost(node.right_child) node.key = right_most_child.key self._delete(right_most_child) def find_rightmost(self, node): if node.right_child is None: return node else: return self.find_rightmost(node.right_child) def find_leftmost(self, node): if node.left_child is None: return node else: return self.find_leftmost(node.left_child)

**Sorting**

Ok great, we have now covered the basics of how to modify a BST. Now we want to know how to extract information from it. In particular given a key, we want to be able to find the next largest or smallest key in the BST, if any. Calling the successor (predecessor) operation repeatedly, starting from the left(right)-most node, will give us the sorted keys in ascending(descending) order.

I will explain how to find the successor, the predecessor case is of course very similar. Let’s suppose we want to find the successor of node with right-child . In this case we choose the left-most child of which we will call .

Note that in a subtree rooted at some node, the left-most node is the smallest element in the subtree, and the right-most node is the largest.

So in the subtree rooted at , the right-subtree rooted at contains all the elements in the subtree greater than . The left-most element of the right-subtree is then the smallest element in the subtree greater than and is a good candidate for a successor. Let’s draw a picture and try to see why this must in fact be the correct node.

OK, but what if has no right child? In this case we move up through the ancestors of , searching for the first ancestor that is a left-child of its parent. If there is no such node, ie we get to the root, then is the right-most, and hence largest, element in the tree, and has no successor. To see that this is correct employ a similar argument to the previous case.

Let’s code that up.

def find_next(self, key): node = self.find(key) if (node is None) or (node.key != key): #Key not in tree. return else: ret= self._find_next(node) if ret is None: #Key is largest in tree. return ret else: return ret.key def _find_next(self, node): right_child = node.right_child if right_child is not None: return self.find_leftmost(right_child) else: parent = node.parent while(parent is not None): if node is parent.left_child: break node = parent parent = node.parent return parent def find_prev(self, key): node = self.find(key) if (node is None) or (node.key != key): #Key not in tree. return else: ret= self._find_prev(node) if ret is None: #Key is smallest in tree. return ret else: return ret.key def _find_prev(self, node): left_child = node.left_child if left_child is not None: return self.find_rightmost(left_child) else: parent = node.parent while(parent is not None): if node is parent.right_child: break node = parent parent = node.parent return parent

By now I think we deserve to see the fruits of our hard work, so let’s whip up a plotting routine using igraphs. Don’t worry about the details: since I can’t find a binary tree plotting style in igraphs, I have colored the root green, left children red and right children blue.

In order to preserve information about left and right children, the plotting routine builds the graph in igraphs using a queue.

def plot(self): #Builds a copy of the BST in igraphs for plotting. #Since exporting the adjacency lists loses information about #left and right children, we build it using a queue. import igraph as igraphs G = igraphs.Graph() if self.root is not None: G.add_vertices(1) queue = [[self.root,0]] #Queue has a pointer to the node in our BST, and its index #in the igraphs copy. index=0 not_break=True while(not_break): #At each iteration, we label the head of the queue with its key, #then add any children into the igraphs graph, #and into the queue. node=queue[0][0] #Select front of queue. node_index = queue[0][1] G.vs[node_index]['label']=node.key if index ==0: #Label root green. G.vs[node_index]['color']='green' if node.left_child is not None: G.add_vertices(1) G.add_edges([(node_index, index+1)]) queue.append([node.left_child,index+1]) G.vs[index+1]['color']='red' #Left children are red. index+=1 if node.right_child is not None: G.add_vertices(1) G.add_edges([(node_index, index+1)]) G.vs[index+1]['color']='blue' queue.append([node.right_child, index+1]) index += 1 queue.pop(0) if len(queue)==0: not_break=False layout = G.layout_reingold_tilford(root=0) igraphs.plot(G, layout=layout)

Now let’s test it!

def test(): lst= [1,4,2,5,1,3,7,11,4.5] B = BinarySearchTree() for item in lst: B.insert(item) B.plot() B.delete(5) B.delete(1) B.plot() print B.root.key print B.find_next(3) print B.find_prev(7) print B.find_prev(1) test()

This produced the following output:

and

“4

4

4.5

None” to the console, as desired.

You can probably see that the operations take time where is the height of the tree. This means that we can sort a length list in time. On average , but if given a monotonic list we will get a linear tree with . So next time we will see how to ‘balance’ the tree using tree rotations to achieve fast operations.

Please feel free to leave any comments, suggestions or corrections! Here is the source code.

## Enter the Perceptron

**Introduction**

Some time ago in my post on logistic regression, I talked about the problem where you have a matrix representing data points with features, and a vector of classifications.

In linear regression you would produce a vector so that , but the problem here is that our output could be any real number, not . The solution we investigated was using the logistic link function which turns any real number into something we interpret as a probability.

That was one solution, another even simpler one is the function defined . It is convenient in this case to think of the response as taking values in not .

Previously we added a column of ones to in order to create an intercept term, here we will write it separately so so that where is called the *bias*.

This post will take the book ‘The Elements of Statistical Learning’ as its guide (the section starting on p. 130), which is freely available online, but there is also a great blog post about it here.

**The Algorithm**

So define , we classify a point as being in the positive class if and negative if . Geometrically we have an affine hyperplane , and points are classified depending on whether they are ‘above’ or ‘below’ the hyperplane.

The unit normal to is given by , points above then mean points with positive projection onto this vector.

The signed distance of a point to is found by projecting the vector onto , where is the ‘origin’ of the hyperplane (since an affine hyperplane is a translation of a linear subspace, there is a point corresponding to the translation of the origin). So we get

,

notice this means that is proportional to the signed distance of from .

If is the set of misclassified points, then we can measure the error as the sum of the distances (up to a constant of proportionality), of the misclassified points to , that is

,

where is the response of , notice that this makes each summand negative, which is reversed by the minus sign outside.

If we fix we can differentiate with respect to and , giving

,

and

.

Just like in logistic regression, we can use these partial derivatives to search for a minimum of . The algorithm purposed by Frank Rosenblatt, the inventor of the perceptron, is to modify the hyperplane after looking at each individual misclassified item in as it comes, that is , where is the learning rate. This has the advantage of being much faster than computing the global gradient informed by all of the data, and is an example of stochastic gradient descent.

Of course by changing the hyperplane as we go we may change , and so what we actually do is iterate over each datapoint in a fixed order, and check if it is misclassified under the latest hyperplane, then update and move onto the next point. This is then very easy to code.

**The Code**

First let’s recycle some code we used for logistic regression to get some test data.

import numpy as np from scipy import linalg import random as rd import matplotlib.pyplot as plt def test_set(a,b, no_samples, intercept, slope, variance, offset): #First we grab some random x-values in the range [a,b]. X1 = [rd.uniform(a,b) for i in range(no_samples)] #Sorting them seems to make the plots work better later on. X1.sort() #Now we define our linear function that will underly our synthetic dataset: the true decision boundary. def g(x): return intercept + slope*x #We create two classes of points R,B. R points are given by g(x) plus a Gaussian noise term with positive mean, #B points are given by g(x) minus a Gaussian norm term with positive mean. R_turn=True X2=[] #This will hold the y-coordinates of our data. Y=[] #This will hold the classification of our data. for x in X1: x2=g(x) if R_turn: x2 += rd.gauss(offset, variance) Y.append(1) else: x2 -= rd.gauss(offset, variance) Y.append(-1) R_turn=not R_turn X2.append(x2) #Now we combine the input data into a single matrix. X1 = np.array([X1]).T X2 = np.array([X2]).T X = np.hstack((X1, X2)) #Now we return the input, output, and true decision boundary. return [ X,np.array(Y).T, map(g, X1)]

Fairly straightforward stuff I hope! If not then it should make sense when we plot it.

X,Y,g = test_set(a=0, b=5, no_samples=200, intercept=10, slope=0, variance=0.5, offset=1) fig, ax = plt.subplots() R = X[::2] #Even terms B = X[1::2] #Odd terms ax.scatter(R[:,0],R[:,1],marker='o', facecolor='red', label="Outcome 1") ax.scatter(B[:,0],B[:,1],marker='o', facecolor='blue', label="Outcome 0") ax.plot(X[:,0], g, color='green', label="True Decision Boundary") ax.legend(loc="upper left", prop={'size':8}) plt.show()

This will produce something like:

Now to implement the algorithm.

def perceptron(X,y, learning_rate=0.1, max_iter=1000): #Expecting X = n x m inputs, y = n x 1 outputs taking values in -1,1 m=X.shape[1] n=X.shape[0] weight = np.zeros((m)) bias = 0 n_iter=0 index = list(range(n)) while n_iter <= max_iter: rd.shuffle(index) for row in index: if (X[row,:].dot(weight)+bias)*y[row] <=0: #Misclassified point. weight += learning_rate*X[row,:]*y[row] bias += learning_rate*y[row] n_iter+=1 return {'weight':weight, 'bias':bias}

we shuffle the order we take the data points in at each step to get rid of any peculiarities in the order we update. Now let’s try it on our test data.

result= perceptron(X,Y,0.1) def decision_boundary(x, weight, bias): return (-bias - x*weight[0])/weight[1] fig, ax = plt.subplots() R = X[::2] #Even terms B = X[1::2] #Odd terms ax.scatter(R[:,0],R[:,1],marker='o', facecolor='red', label="Outcome 1") ax.scatter(B[:,0],B[:,1],marker='o', facecolor='blue', label="Outcome -1") ax.plot(X[:,0], g, color='green', label="True Decision Boundary") ax.plot(X[:,0], [decision_boundary(x, result['weight'], result['bias']) for x in X[:,0]], label='Predicted Decision Boundary') ax.legend(loc="upper left", prop={'size':8}) plt.show()

Producing:

Not bad! Get the code here.

## Snakes on a Manifold: Part 2

Welcome back. Last time we implemented the main components of the game. This time we will implement the main game loop. The job of the loop is to: get keyboard input, update the game state, update the display.

We will now be using pygame, so you might want to check out a tutorial before proceeding, but I think the best way to learn is just to dive right in. You will however need to install pygame before running the code. If you haven’t already, you should download the code from last time too. Finally if you want to skip to playing with the code yourself, here it is.

**The Code**

Ok here it is.

import pygame import snake as sk from pygame.locals import * class Game(object): def __init__(self, size, topology, speed=6): pygame.init() #Resolution in pixels. self.screen_resolution =600 self.topology=topology self.speed = speed #Controls how fast the snake moves. self.screen = pygame.display.set_mode((self.screen_resolution, self.screen_resolution)) pygame.display.set_caption('SNAKES ON A MANIFOLD!') self.board = sk.Board(size, topology) #Put the snake in the center. self.snake = sk.Snake([size/2, size/2], self.board.opposite_direction, self.board.get_next_coord) #Number of grid cells. self.size=size #Width in pixels of a grid cell. self.square_size=int(self.screen_resolution/self.size) #Now we name some colors in RGB for ease of use. self.black = (0,0,0) self.red = (255,0,0) self.green = (0,255,0) self.blue=(0,0,255) self.darkgrey=(40,40,40) self.background_colour=self.black self.boundary_colour=self.blue self.snake_colour = self.green self.apple_colour=self.red #Initialises screen to background color. self.screen.fill(self.background_colour) #Then we draw the boundary cells. self.draw_boundary() #This command updates the changes we have made to the screen. pygame.display.update() #Main game loop. self.main() def main(self): score=0 Quit=False first=True while True: # main game loop #Displays score in window caption. pygame.display.set_caption('Topology: '+self.board.topology+' | Score: '+ str(score)) #Gets keyboard input direction=None for event in pygame.event.get(): if event.type == QUIT: #User has closed window. Quit=True break if event.type==pygame.KEYDOWN: #User has pressed key. if event.key==K_q: Quit=True break if event.key==K_UP: direction='UP' if event.key==K_DOWN: direction='DOWN' if event.key==K_LEFT: direction='LEFT' if event.key==K_RIGHT: direction='RIGHT' if Quit: break #To see if anything has happened, we check where the head is. head_loc = self.snake.head_position #Now we find out what, apart from the head, is in the grid cell. head_on = self.board.grid[head_loc[0], head_loc[1]] if (head_on == 1) or (head_on==2): #The snake has hit the boundary or part of itself. Quit=True break if head_on==3: #The snake has found an apple. self.snake.grow() x,y=self.board.generate_apple() self.draw_apple(x,y) score+=1 self.update_board() if first: #When the game starts, the snake stands still waiting #for user input to start the game. x,y=self.board.generate_apple() first=False self.draw_apple(x,y) self.draw_snake() pygame.display.update() while True: got_direction=False for event in pygame.event.get(): if event.type == QUIT: Quit=True break if event.type==pygame.KEYDOWN: if event.key==K_q: Quit=True break if event.key==K_UP: direction='UP' got_direction=True break if event.key==K_DOWN: direction='DOWN' got_direction=True break if event.key==K_LEFT: direction='LEFT' got_direction=True break if event.key==K_RIGHT: direction='RIGHT' got_direction=True break if Quit: break elif got_direction: break if Quit: break self.snake.move(direction) self.draw_snake() pygame.display.update() #This pauses the game to slow it down, adjust to make faster/slower. pygame.time.Clock().tick(self.speed) #Waiting to quit after died. Quit=False restart=False while True: # main game loop #Gets keyboard input for event in pygame.event.get(): if event.type == QUIT: Quit=True break if event.type==pygame.KEYDOWN: if event.key==K_q: Quit=True break if event.key==K_r: restart=True break if Quit: break if restart: break if restart: G=Game(self.size, self.topology) pygame.quit() def draw_boundary(self): #Draws any boundary grid cells. for row in range(self.size): for col in range(self.size): if self.board.grid[row, col]==1: y=row*self.square_size x=col*self.square_size rect=pygame.Rect(x,y, self.square_size, self.square_size) pygame.draw.rect(self.screen,self.boundary_colour, rect) def draw_snake(self): #We add a snake colored square where the head has moved, #and add a background colored square where the tail was. old_tail_position=self.snake.segments[0].previous_position new_head_position=self.snake.head_position top=old_tail_position[0]*self.square_size left=old_tail_position[1]*self.square_size background_rect=pygame.Rect(left,top, self.square_size, self.square_size) pygame.draw.rect(self.screen, self.background_colour,background_rect) top=new_head_position[0]*self.square_size left=new_head_position[1]*self.square_size snake_rect=pygame.Rect(left,top, self.square_size, self.square_size) pygame.draw.rect(self.screen, self.snake_colour,snake_rect) def draw_apple(self, x, y): apple_rect=pygame.Rect(y*self.square_size,x*self.square_size, self.square_size, self.square_size) pygame.draw.rect(self.screen, self.apple_colour, apple_rect) def update_board(self): #Updates the numpy array of the board. old_tail_position = self.snake.segments[0].previous_position new_head_position=self.snake.head_position self.board.grid[old_tail_position[0], old_tail_position[1]]=0 self.board.grid[new_head_position[0], new_head_position[1]]=2

And we’re done! Now you can play a game by typing something like this.

G=Game(20, 'Mobius Strip')

Please let me know what you think. The code is here. If anyone is feeling ambitious, you could add a menu screen where you could select a speed and topology.

## Computing Homology

Last time we introduced abstract simplicial complexes (abs) and how to represent them in Python. Today we will look at the homology groups of an abs.

For an excellent primer on homology I recommend checking out this post , I will just quickly recap the basics.

**Homology: A quick reminder**

Let be an abs.

Write for the dimensional skeleton of , formed by restricting to those elements of with dimension .

If is a -simplex with vertices labelled we may write to denote this. To denote the -simplex consisting of the vertices of except the vertex we write .

For each dimension we write for the set of -dimensional simplices, and write for its cardinality.

For each dimension we write for the –*dimensional chain group* given by the group presentation . What this says is that the generators of the group are the -simplices, the group is abelian and each element is its own inverse.

An alternative way to think of this is as a vector space, where has basis and is over the field (mod 2 arithmetic). So every element can be represented as a binary vector of length , and addition is bitwise xor.

Now we define boundary maps . These are both group homomorphisms and linear maps. For a -simplex define:

, the sum of its -dimensional faces.

Since is a a basis for this determines a linear mapping.

Define called the –*dimensional cycle group* (subspace).

Define called the –*dimensional boundary group* (subspace).

Now we can define the –*dimensional homology group* .

This is both a quotient group and a quotient vector space.

Finally define the *Betti number* .

**Computation**

We are going to be interested in the Betti numbers as these can be thought of as counting the -dimensional holes in our space. For this just means the number of connected components.

In order to calculate the ranks of interest, we must construct the matrices representing the boundary mappings. In order to do this we must agree upon an ordering of the simplices. So we simply use the one generated by the data structure, which is why each simplex has an index attribute.

To find the rank of , which we have called , we simply count the number of -simplices. So we can add that feature to our SimplicialComplex class easily.

def chain_rank(self,p): return len(self.n_simplex_dict[p])

Now the matrix of will be have rows, one for each -simplex, and columns, one for each -simplex. Then the column of the matrix will contain the coefficients of the image of the simplex expressed in terms of our basis for .

So to find the column vector, we ask each of a simplex’s boundary faces for their index. These will be the positions of the 1s in the column.

def boundary_column(self,simplex): column_size = self.chain_rank(simplex.dimension-1) faces = simplex.faces column = [0 for i in range(column_size)] for face in faces: column[face.index]=1 return column

Now we simply combine these columns to form the boundary matrix. Note that for greater than the dimension of the simplicial complex, the trivial group, and similarly for negative . We don’t have to calculate anything in this case.

This also means that the kernel of is . This is handled by us having the empty set as a kind of -dimensional simplex.

def get_boundary_matrix(self,p, skeleton=None): if skeleton==None: skeleton=self.dimension if ((p>skeleton)or(p<0)): return np.asarray([]) mat=[self.boundary_column(simplex) for simplex in self.n_simplex_dict[p]] return np.asarray(mat, dtype=int).transpose()

The skeleton argument is used when you want to only consider simplices up to a certain dimension. For instance if is given by then is a representation of a triangle. Then the zeroth Betti number as there is one connected component. However the zeorth Betti number of is 3, one for each vertex.

Let’s test what we’ve got so far. In this post Jeremy Kun computes the homology of the tetrahedron with an ‘arm’, he provides a nice picture of it here. If we label the vertices of the tetrahedron , then this structure can be specified with the maximal simplices then restricting to the -skeleton, since we want there to be a hole in the tetrahedron.

Let’s see what the boundary matrices are:

def jeremy_example(): maximal_simplices = [[0,1,2,3],[0,4],[2,4]] S = SimplicialComplex(maximal_simplices) S.plot(direction='up', label_color='blue') for i in range(3): print "Boundary matrix of dimension ", i print S.get_boundary_matrix(i,skeleton=2) jeremy_example()

Which produced the following plot (of the 3-skeleton):

And outputted to the console:

**Boundary matrix of dimension 0**

** [[1 1 1 1 1]]**

** Boundary matrix of dimension 1**

** [[1 1 1 0 0 0 1 0]**

** [1 0 0 1 1 0 0 0]**

** [0 1 0 1 0 1 0 1]**

** [0 0 1 0 1 1 0 0]**

** [0 0 0 0 0 0 1 1]]**

** Boundary matrix of dimension 2**

** [[1 1 0 0]**

** [1 0 1 0]**

** [0 1 1 0]**

** [1 0 0 1]**

** [0 1 0 1]**

** [0 0 1 1]**

** [0 0 0 0]**

** [0 0 0 0]]**

If you check his post, you will see that the matrices are the same, up to a reordering of basis, and noticing that

So far so good. Now all that remains is to calculate the rank and nullity of each matrix. But this is easily done, we simply diagonalize the matrix and count the pivots. For the number of pivots is equal to and by the rank-nullity theorem.

The form we will put the matrix in is called Smith-Normal form I won’t go through the details of how the algorithm works, they are on the wiki page. If you have seen Gaussian elimination before this is basically that.

def reduce_matrix(self, matrix): #Returns [reduced_matrix, rank, nullity] if np.size(matrix)==0: return [matrix,0,0] m=matrix.shape[0] n=matrix.shape[1] def _reduce(x): #We recurse through the diagonal entries. #We move a 1 to the diagonal entry, then #knock out any other 1s in the same col/row. #The rank is the number of nonzero pivots, #so when we run out of nonzero diagonal entries, we will #know the rank. nonzero=False #Searching for a nonzero entry then moving it to the diagonal. for i in range(x,m): for j in range(x,n): if matrix[i,j]==1: matrix[[x,i],:]=matrix[[i,x],:] matrix[:,[x,j]]=matrix[:,[j,x]] nonzero=True break if nonzero: break #Knocking out other nonzero entries. if nonzero: for i in range(x+1,m): if matrix[i,x]==1: matrix[i,:] = np.logical_xor(matrix[x,:], matrix[i,:]) for i in range(x+1,n): if matrix[x,i]==1: matrix[:,i] = np.logical_xor(matrix[:,x], matrix[:,i]) #Proceeding to next diagonal entry. return _reduce(x+1) else: #Run out of nonzero entries so done. return x rank=_reduce(0) return [matrix, rank, n-rank] def boundary_rank(self,p, skeleton=None): if skeleton==None: skeleton=self.dimension if (p >= skeleton) or (p<0): #Trivial case. return 0 else: return self.reduce_matrix(self.get_boundary_matrix(p+1, skeleton))[1] def cycle_rank(self, p, skeleton=None): if skeleton==None: skeleton=self.dimension if p==0: #Mapping to trivial group so kernel is everything in this case. return self.chain_rank(0) elif (p>skeleton) or (p<0): #Trivial case. return 0 else: return self.reduce_matrix(self.get_boundary_matrix(p, skeleton))[2] def betti_number(self, p, skeleton=None): if skeleton==None: skeleton=self.dimension return self.cycle_rank(p,skeleton) - self.boundary_rank(p, skeleton)

Now let’s just test it out on the same example. We expect since there is one connected component, from the -dimensional hole made by the arm and for the ball inside the tetrahedron’s shell.

def jeremy_example(skeleton=2): maximal_simplices = [[0,1,2,3],[0,4],[2,4]] S = SimplicialComplex(maximal_simplices) S.plot(direction='up', label_color='blue') for i in range(4): print "Boundary matrix of dimension", i print S.get_boundary_matrix(i,skeleton) print "Reduced boundary matrix of dimension", i print S.reduce_matrix(S.get_boundary_matrix(i,skeleton))[0] print "Betti number",i,"=", S.betti_number(i,skeleton) jeremy_example()

The console outputs:

**Boundary matrix of dimension 0**

** [[1 1 1 1 1]]**

** Reduced boundary matrix of dimension 0**

** [[1 0 0 0 0]]**

** Betti number 0 = 1**

** Boundary matrix of dimension 1**

** [[1 1 1 0 0 0 1 0]**

** [1 0 0 1 1 0 0 0]**

** [0 1 0 1 0 1 0 1]**

** [0 0 1 0 1 1 0 0]**

** [0 0 0 0 0 0 1 1]]**

** Reduced boundary matrix of dimension 1**

** [[1 0 0 0 0 0 0 0]**

** [0 1 0 0 0 0 0 0]**

** [0 0 1 0 0 0 0 0]**

** [0 0 0 1 0 0 0 0]**

** [0 0 0 0 0 0 0 0]]**

** Betti number 1 = 1**

** Boundary matrix of dimension 2**

** [[1 1 0 0]**

** [1 0 1 0]**

** [0 1 1 0]**

** [1 0 0 1]**

** [0 1 0 1]**

** [0 0 1 1]**

** [0 0 0 0]**

** [0 0 0 0]]**

** Reduced boundary matrix of dimension 2**

** [[1 0 0 0]**

** [0 1 0 0]**

** [0 0 1 0]**

** [0 0 0 0]**

** [0 0 0 0]**

** [0 0 0 0]**

** [0 0 0 0]**

** [0 0 0 0]]**

** Betti number 2 = 1**

** Boundary matrix of dimension 3**

** []**

** Reduced boundary matrix of dimension 3**

** []**

** Betti number 3 = 0**

Just as expected. So we’re done. Get the code here.