Shuffling a List

Here are two algorithms that randomly shuffle a list. There are also two implementations of the swap function. One is standard, the other is bitwise. There is a demo() function that you should run on short lists to convince yourself that shuffle is indeed doing what it claims to do.

#! /usr/bin/env python

from random import randint
import time
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

def swap(lst, i, j):
 temp = lst[i]
 lst[i] = lst[j]
 lst[j] = temp

def swap_bitwise(lst, i, j):
 lst[i] = lst[i] ^ lst[j]
 lst[j] = lst[i] ^ lst[j]
 lst[i] = lst[i] ^ lst[j]

def shuffle(lst, swap_func, should_print=False):
 for i in xrange(len(lst)):
  j = randint(0, len(lst)-1)
  swap_func(lst, i, j)
  if should_print:
   print i, j, lst

# more or less due to Fisher, Yates, Durstenfeld, and Knuth
def shuffle2(lst, swap_func, should_print=False):
 for i in range(len(lst))[::-1]:
  j = randint(0, i)
  swap_func(lst, i, j)
  if should_print:
   print i, j, lst

def demo(shuffle_func):
 lst = range(10)
 print 'list before:', lst
 print 'shuffling ...'
 shuffle_func(lst, swap, True)
 print '...done'
 print 'list after:', lst

def function_times(fs, maxlen):
 res = []
 for f in fs:
  r = []
  for i in xrange(maxlen):
   data = range(i)
   start = time.time()
   shuffle(data, f)
   end = time.time()
   r.append((i, end-start))
  res.append(r)
 return res

def function_times2(fs, maxlen):
 res = []
 for f in fs:
  r = []
  for i in [1000*i for i in xrange(maxlen)]:
   data = range(i)
   start = time.time()
   shuffle(data, f)
   end = time.time()
   r.append((i, end-start))
  res.append(r)
 return res

def plot_prob():
 tests = 10000
 n_elements = 20
 m = np.array([[0.]*n_elements]*n_elements)
 for k in xrange(tests):
  initial = range(n_elements)
  shuffle2(initial, swap)
  for i in range(len(initial)):
   m[i][initial[i]] += 1
 m = m / tests
 plt.imshow(m)
 plt.colorbar()
 plt.title('Color Map')
 plt.show()
 plt.clf()
 #print m.reshape(-1,).tolist()
 plt.title('Histogram')
 plt.hist(m.reshape(-1,).tolist())
 plt.show()
 plt.clf()

def plot_times():
 functions = [swap, swap_bitwise]
 mlen = 1000
 res = function_times(functions, mlen)
 fig = plt.figure(dpi=90)
 ax = fig.add_subplot(111)
 for r in res:
  ax.plot([rr[0] for rr in r], [rr[1] for rr in r])
 ax.legend(['swap', 'bitwise swap'], 
  loc='upper center', 
  #bbox_to_anchor=(0.5, 1.05),
  ncol=2, fancybox=True, shadow=True)
 #ax.set_yscale('log')
 plt.title('Shuffle function runtimes versus list length')
 plt.xlabel('list length')
 plt.ylabel('time (seconds)')
 plt.show()

if __name__ == '__main__':
 #demo(shuffle)
 #plot_times()
 plot_prob()

Below is the printout of the demo() function. 

 Below is a plot of the time to shuffle, the output of the plot_times() function. As you can see, the bitwise swap is a little slower (it takes more time) than the regular swap. 

Now, we change the function function_times() to shuffle much larger lists. Since it takes quite a while to iterate through every size 1, 2, ..., n, we will change the list sizes to 1000, 2000, ..., 200000. This is accomplished by changing one line in the function_times function, from

for i in xrange(maxlen):

to 

for i in [1000*i for i in xrange(maxlen)]:

This functionality is incorporated into the function_times2() function. Here is the plot:

Notice the fluctuations near x = 150,000. This is most likely due to the fact that the list is so large that it cannot fit into cache, i.e., elements of the list are being read from RAM. If we continue to increase the list size, there will be another spike where the list is being read from disk, which is even slower. If this doesn't make sense to you, look up 'paging' and the 'principle of locality' to see how data is moved between cache, RAM, and disk. 
Finally, let's compare shuffle() and shuffle2(). There are subtle differences that you can see, but let's look at the output in two ways: a colormap and a histogram. First, the colormap for shuffle():

Since we are using 20 values, the probability that i goes to j should be 1/20 = 0.05. There is, however, a non-uniform distribution of values. The histogram below confirms this: if there was no bias, the values would assume a normal (gaussian) distribution around 0.05, but the data is skewed right. 
The histogram for shuffle():

Here's the colormap for shuffle2():

This looks much more uniform. Indeed, as the histogram for shuffle2() function below shows, the data is almost normally distributed. Notice the nice round shape and lack of (or less) skew. 

For this reason, I recommend using shuffle2().  

Permutations

Below are five functions to calculate all permutations of a list. Some are mine, some are borrowed, all are valid. 

#! /usr/bin/env python

import time
import matplotlib
import matplotlib.pyplot as plt

def perm1(lst):
 if len(lst) == 0:
  return []
 elif len(lst) == 1:
  return [lst]
 else:
  l = []
  for i in range(len(lst)):
   x = lst[i]
   xs = lst[:i] + lst[i+1:]
   for p in perm1(xs):
    l.append([x] + p)
  return l

def perm2(lst):
 if len(lst) == 0:
  yield []
 elif len(lst) == 1:
  yield lst
 else:
  for i in range(len(lst)):
   x = lst[i]
   xs = lst[:i] + lst[i+1:]
   for p in perm2(xs):
    yield [x] + p

# http://stackoverflow.com/questions/2710713/algorithm-to-generate-all-possible-permutations-of-a-list
# from zhywu at stackoverflow.com
def perm3(lst):
 if len(lst) <= 2:
  yield lst
  lst = lst[::-1]
  yield lst
 else:
  l = []
  for i in range(len(lst)):
   x = lst[i]
   xs = lst[:i] + lst[i+1:]
   for p in perm3(xs):
    yield [x]+p

# from Peter Breuer at stackoverflow (ported from his Haskell implementation)
def perm4(lst):
 if len(lst) <= 0:
  yield []
 elif len(lst) == 1:
  yield lst
 else:
  for p in perm4(lst[1:]):
   for a, b in splits(p):
    yield a + [lst[0]] + b

# list of ways of splitting a list into two parts
def splits(lst):
 for i in xrange(len(lst)+1):
  yield [lst[:i], lst[i:]]

# found this on internet, not sure where
def perm5(xs, low=0):
    if low + 1 >= len(xs):
        yield xs[:]
    else:
        for p in perm5(xs, low + 1):
            yield p
        for i in range(low + 1, len(xs)):
            xs[low], xs[i] = xs[i], xs[low]
            for p in perm5(xs, low + 1):
                yield p        
            xs[low], xs[i] = xs[i], xs[low]

def function_times(fs, maxlen=5):
 res = []
 for f in fs:
  r = []
  for i in range(maxlen):
   data = range(i)
   start = time.time()
   for p in f(data):
    pass
   end = time.time()
   r.append(end-start)
  res.append(r)
 return res

def all_results_equal(fs, data):
 p1 = sorted(perm1(data))
 p2 = sorted([p for p in perm2(data)])
 p3 = sorted([p for p in perm3(data)])
 p4 = sorted([p for p in perm4(data)])
 p5 = sorted([p for p in perm5(data)])
 return p1 == p2 == p3 == p4 == p5

if __name__ == '__main__':
 functions = [perm1, perm2, perm3, perm4, perm5]
 assert all_results_equal(functions, range(4))

 res = function_times(functions, 9)
 fig = plt.figure(dpi=90)
 ax = fig.add_subplot(111)
 for r in res:
  ax.plot(r)
 ax.legend(['perm1', 'perm2', 'perm3', 'perm4', 'perm5'], 
  loc='upper center', 
  #bbox_to_anchor=(0.5, 1.05),
  ncol=2, fancybox=True, shadow=True)
 ax.set_yscale('log')
 plt.title('Permutation function runtimes versus list length')
 plt.xlabel('list length')
 plt.ylabel('time (seconds)')
 plt.show()

Below are two plots of the same data--one using linear y axis, the other using logarithmic y axis. As you can see, the fourth function (teal line) performs slightly better than the rest, for large lists, although none of the functions performs that well on large lists. 

Combinations

A combination is an instance of one possible way to choose k things from n things. For example, all combinations of size 2 from the list ['a', 'b', 'c'] are ['a', 'b'], ['b', 'c'], and ['a', 'c']. The number of combinations is just famed n-choose-k function: $ \binom{n}{k} = \frac{n!}{k!(n-k!)} $. In the example above, n = 3 and k = 2, so  the number of combinations is $ \frac{3!}{2!(3-2!)} = \frac{3!}{2! 1!} = \frac{3!}{2!} = \frac{3*2*1}{2*1} = 3 $. The following functions compute all combinations of a list. Each has advantages and disadvantages, depending on your needs. The first function, combs1(), makes copies of the list in the recursive call. This will effect performance for large lists. However, the second function, combs2(), makes no copies, which is why it requires a little more manipulation.

#! /usr/bin/env python

def combs1(k, available, used=[]):
 if len(used) == k:
  yield tuple(used)
 elif len(available) <= 0:
  pass
 else:
  for c in combs1(k, available[1:], used[:]+[available[0]]):
   yield c
  for c in combs1(k, available[1:], used[:]):
   yield c

def combs2(k, available, used=[]):
 if len(used) == k:
  yield tuple(used)
 elif len(available) <= 0:
  pass
 else:
  head = available.pop(0)
  used.append(head)
  for c in combs2(k, available, used):
   yield c
  used.pop(-1)
  for c in combs2(k, available, used):
   yield c
  available.insert(0, head)

def testfunctions(testlist, k=2):
 mycombs = []
 testfunctions = [combs1, combs2]
 for tf in testfunctions:
  mycombs.append([s for s in tf(k, testlist)])
  
 # this is for validation
 import itertools
 combs_from_lib = [i for i in itertools.combinations(testlist, k)]
 
 results = []
 for mycomb in mycombs:
  if mycomb == combs_from_lib:
   results.append('correct')
  else:
   results.append('wrong')
 return results

if __name__ == '__main__':
 print testfunctions(list('abc'))

For validation, I import the itertools.combinations and compare the results from my functions. They match!

Majority Algorithm

This algorithm is for finding the majority item in a list or stream of items. Here, for simplicity, all items are integers. This algorithm is known as the Majority Algorithm, and was first published by Fisher and Salzberg in 1982. The point of this algorithm is to read a stream of n items, and return the majority item if it exists. A majority item is defined as an item that occurs more than n/2 times in a stream with n items. The algorithm does this with one counter. 
The algorithm first initializes the counter to 0 and the majority element to null. Then it looks at each item in the stream. For each item, if the counter is 0, the item is set to the majority and the counter is set to 1. If the item is the current majority item, the counter is incremented by 1. If the item is not the current majority item, the counter is decremented by 1. 
For testing, I have created a little stream whose majority item is the integer 5. 

#! /usr/bin/env python

from random import choice

def find_majority(data_stream):
 majority = None
 counter = 0
 for item in data_stream:
  if item == majority:
   counter += 1
  else:
   if counter == 0:
    majority = item
    counter = 1
   else:
    counter -= 1
  print 'item:', item, 'majority:', majority, 'counter:', counter
 return majority

def stream(items, maxi=2000):
 i = 0
 while True:
  yield choice(items)
  i += 1
  if i > maxi:
   break

if __name__ == '__main__':
 l = range(10) + [5]*10
 print 'stream:', l
 s = stream(l, len(l))
 maj = find_majority(s)
 print 'majority:', maj

Sample:

Since 5 is barely the majority, the counter hovers around 0 and 1, but the final result is correct. If 5 was the standout majority, the counter would increase, as so:

Note that if the counter is 0 at the end of the algorithm, the last item seen by the algorithm could have occurred up to n/2 times, but not more than n/2. Also note that the stream() function above is stochastic, i.e., it takes a list and emits random values from that list. If you do not want a stochastic stream, just pass a regular old list or array to the algorithm. It can also be easily extended to read streams from disk. Just open a file and pass the reference to the majority algorithm. In python, the code

f = open('myfile')
for item in f:
    # do something

is different than

f = open('myfile')
for item in f.readlines():
    # do something

The first code snippet will read from disk one line at a time, therefore memory requirements will be low, but at the cost of repeated disk accesses. This is slower, but if the data stream is so large that it does not fit into RAM, then this is a good solution. The second code snippet reads the entire file into memory first, then starts iterating, which is only good if the entire file can fit into RAM. 

Anagrams using MapReduce and mrjob

This code discovers anagrams. The input is a text file dictionary with one word per line. The output is a text file.  

#! /usr/bin/env python

from mrjob.job import MRJob

class Anagram(MRJob):
 def mapper(self, _, word):
  letters = list(word)
  letters.sort()
  if word != '':
   yield letters, word

 def reducer(self, _, words):
  anagrams = [w for w in words]
  if len(anagrams) >= 2:
   yield len(anagrams), anagrams

if __name__ == '__main__':
 Anagram.run()

Call it like this:

Sample output:

The output file consists of two columns: the first column is the number of word in the set, and the second column is the set itself. 

Perceptron

The perceptron algorithm attempts to find a line of separation between two classes of data--hence it is a classification algorithm. 
Below is the python code for the perceptron algorithm. We first generate a dataset with two dimensional input and one-dimensional output. The input is normally distributed along both axes. The output, y, is either +1 or -1. The training set is a list of tuples of the form ([x1, x2], y) where the first element of the tuple is a list (actually a numpy array in the code, but for simplicity, I write it as a list), and the second element is the output. 
Then we feed the training set to the learning algorithm. The perceptron function initializes the weight vector, w, randomly. It then iterates over the training set data. At each iteration, we multiply w by x, resulting in the dot product over the two vectors. If this dot product is positive, then the point x is on the correct side of the dividing line. If the dot product is negative, then the point x is on the wrong side of the dividing line, and therefore the classifier is wrong. To correct this wrong, we nudge the line slightly.  

#! /usr/bin/env python

import math
import numpy as np
from numpy.random import normal
from pylab import plot, cla, show, norm, rand, xlabel, ylabel, title, legend, get_cmap, savefig, xlim, ylim
from random import shuffle, random
import matplotlib.pyplot as plt

def generateDataNormal(n):
 data = []
 mu, sigma = 0.15, 0.15 # mean and standard deviation
 x1 = np.random.normal(mu, sigma, n)
 y1 = np.random.normal(mu, sigma, n)
 data.extend([(np.array([x1[i], y1[i]]), 1) for i in range(n)])
 mu, sigma = -0.15, 0.15 # mean and standard deviation
 x2 = np.random.normal(mu, sigma, n)
 y2 = np.random.normal(mu, sigma, n)
 data.extend([(np.array([x2[i], y2[i]]), -1) for i in range(n)])
 return data

def test_error(w, test_data):
 correct = 0
 for x_i, y_i in test_data:
  y_i_prime = np.dot(w, x_i)
  if y_i * y_i_prime >= 0:
   correct += 1
 e = (len(test_data)-float(correct))/len(test_data)
 return e

def perceptron_weights(training_data, lrate=0.2):
 w = []
 xlen = len(training_data[0][0])
 for i in xrange(xlen):
  w.append(random())
 w = np.array(w)
 for x_i, y_i in training_data:
  y_i_prime = np.dot(w, x_i)
  if y_i * y_i_prime <= 0.0:
   w = w + lrate*y_i*x_i
   yield w.copy()

def perceptron_weight(training_data, lrate=0.2):
 w = []
 xlen = len(training_data[0][0])
 for i in xrange(xlen):
  w.append(random())
 w = np.array(w)
 for x_i, y_i in training_data:
  y_i_prime = np.dot(w, x_i)
  if y_i * y_i_prime <= 0.0:
   w = w + lrate*y_i*x_i
 return w

def plot_data(training_data, ws, lrate):
 for x in training_data:
  if x[-1] == 1:
   plot(x[0][0], x[0][1], 'ob')
  else:
   plot(x[0][0], x[0][1], 'or')
 colors = np.linspace(0.1, 1, len(ws))
 mymap = get_cmap("Greys")
 mycolors = mymap(colors)
 i = 0
 for w in ws:
  n = norm(w)
  ww = w/n
  ww1 = [ww[1], -ww[0]]
  ww2 = [-ww[1], ww[0]]
  plot([ww1[0], ww2[0]],[ww1[1], ww2[1]],color=mycolors[i])
  i += 1
 xlim(-1, 1)
 ylim(-1, 1)
 savefig('perceptron_'+str(lrate)+'.png')
 cla()
 #show()

if __name__ == '__main__':
 training_data = generateDataNormal(100)
 shuffle(training_data)
 for lrate in [0.1, 0.2, 0.4, 0.8]:
  weights = [w for w in perceptron_weights(training_data, lrate)]
  plot_data(training_data, weights, lrate)

Below are three plots of how w changes on each iteration (where it is wrong and w must be changed) for different learning rates, represented by the Greek letter lambda. A note about the color mapping for the lines (which are actually the lines orthogonal to w): there is a line for each time w is updated, and the first line is the lightest, and the last line is the darkest. As you can see, small learning rates like 0.1 result in small changes in w and larger learning rates like 1.0 result in larger change. Caution should be taken when choosing lambda to somewhat guarantee convergence and avoid thrashing or cycling of w. 
The blue circles are one class and the red circles are the other. 

Counting Words using MapReduce

This class counts words in a document.

#! /usr/bin/env python

import re
from mrjob.job import MRJob

class WordCount(MRJob):
 def mapper(self, key, value):
  for word in value.split():
   yield word.lower(), 1
 def reducer(self, key, value):
  yield key, sum(value)

if __name__ == '__main__':
 WordCount.run()

Assuming you have Jane Austin's Sense and Sensibility in a .txt file like I do, call it like so: 

Each output line contains two tab-delimited elements: the word, and the count of that word. 

Oh, dear, that's a lot of 'dear's. As you can see, the mapper function simply defines words as anything between two spaces, or space characters. Additional work can clean that up a bit. The following fix uses regular expressions--regex. The special character '\w' means 'a word character', which means a-z or A-Z or 0-9. The '+' means 'one or more'. Putting these together, '\w+' means one or more word characters, which comprises a word. The 'ws' variable contains all matching words in the 'word' string. So for example, if 'word' is "dear--sure!", then 'ws' is all occurrences of one or more word characters, which for this example is "dear" and "sure."

def mapper(self, key, value):
  for word in value.split():
   ws = re.findall(r'\w+', word)
   for w in ws:
    yield w.lower(), 1

The occurrences of the word 'dear' and its kin have been stripped of the non-word characters, giving better knowledge of the words themselves, but loosing some of the peripheral info contained in the punctuation; that is, 'dear' does not equal 'dear!'