diff --git a/nlib.py b/nlib.py
index 93b1a45..4bf2e19 100644
--- a/nlib.py
+++ b/nlib.py
@@ -1,6 +1,39 @@
 # Created by Massimo Di Pierro - BSD License
 
-__all__ = ['CONSTANT', 'CUBIC', 'Canvas', 'Cholesky', 'Cluster', 'D', 'DD', 'Dijkstra', 'DisjointSets', 'E', 'Ellipse', 'Figure', 'FigureCanvasAgg', 'HAVE_MATPLOTLIB', 'Jacobi_eigenvalues', 'Kruskal', 'LINEAR', 'MCEngine', 'MCG', 'Markowitz', 'MarsenneTwister', 'Matrix', 'NeuralNetwork', 'POLYNOMIAL', 'PersistentDictionary', 'Prim', 'PrimVertex', 'QUADRATIC', 'QUARTIC', 'QuadratureIntegrator', 'RandomSource', 'Trader', 'YStock', 'bootstrap', 'breadth_first_search', 'compute_correlation', 'condition_number', 'confidence_intervals', 'continuum_knapsack', 'correlation', 'covariance', 'decode_huffman', 'depth_first_search', 'encode_huffman', 'fib', 'fit', 'fit_least_squares', 'gradient', 'hessian', 'integrate', 'integrate_naive', 'integrate_quadrature_naive', 'invert_bicgstab', 'invert_minimum_residual', 'is_almost_symmetric', 'is_almost_zero', 'is_positive_definite', 'jacobian', 'lcs', 'leapfrog', 'make_maze', 'mean', 'memoize', 'memoize_persistent', 'needleman_wunsch', 'norm', 'optimize_bisection', 'optimize_golden_search', 'optimize_newton', 'optimize_newton_multi', 'optimize_newton_multi_imporved', 'optimize_secant', 'partial', 'resample', 'sd', 'solve_bisection', 'solve_fixed_point', 'solve_newton', 'solve_newton_multi', 'solve_secant', 'variance']
+__all__ = ['CONSTANT', 'CUBIC', 'Canvas', 'Cholesky', 'Cluster', 'D', 'DD',
+           'Dijkstra', 'DisjointSets', 'E', 'Ellipse',
+           'Figure', 'FigureCanvasAgg',
+           'HAVE_MATPLOTLIB', 'Jacobi_eigenvalues',
+           'Kruskal', 'LINEAR', 'MCEngine',
+           'MCG', 'Markowitz', 'MarsenneTwister',
+           'Matrix', 'NeuralNetwork',
+           'POLYNOMIAL', 'PersistentDictionary',
+           'Prim', 'PrimVertex', 'QUADRATIC',
+           'QUARTIC', 'QuadratureIntegrator',
+           'RandomSource', 'Trader', 'YStock',
+           'bootstrap', 'breadth_first_search',
+           'compute_correlation', 'condition_number',
+           'confidence_intervals', 'continuum_knapsack',
+           'correlation', 'covariance',
+           'decode_huffman', 'depth_first_search',
+           'encode_huffman', 'fib', 'fit',
+           'fit_least_squares', 'gradient', 'hessian',
+           'integrate', 'integrate_naive',
+           'integrate_quadrature_naive',
+           'invert_bicgstab', 'invert_minimum_residual',
+           'is_almost_symmetric', 'is_almost_zero',
+           'is_positive_definite', 'jacobian',
+           'lcs', 'leapfrog', 'make_maze', 'mean',
+           'memoize', 'memoize_persistent',
+           'needleman_wunsch', 'norm',
+           'optimize_bisection', 'optimize_golden_search',
+           'optimize_newton', 'optimize_newton_multi',
+           'optimize_newton_multi_imporved',
+           'optimize_secant', 'partial',
+           'resample', 'sd', 'solve_bisection',
+           'solve_fixed_point', 'solve_newton',
+           'solve_newton_multi', 'solve_secant', 'variance']
+
 
 class YStock:
     """
@@ -16,7 +49,8 @@ class YStock:
     """
     URL_CURRENT = 'http://finance.yahoo.com/d/quotes.csv?s=%(symbol)s&f=%(columns)s'
     URL_HISTORICAL = 'http://ichart.yahoo.com/table.csv?s=%(s)s&a=%(a)s&b=%(b)s&c=%(c)s&d=%(d)s&e=%(e)s&f=%(f)s'
-    def __init__(self,symbol):
+
+    def __init__(self, symbol):
         self.symbol = symbol.upper()
 
     def current(self):
@@ -45,24 +79,26 @@ def current(self):
         url = self.URL_CURRENT % dict(symbol=self.symbol, columns=columns)
         raw_data = urllib.urlopen(url).read().strip().strip('"').split(',')
         current = dict()
-        for i,row in enumerate(FIELDS):
+        for i, row in enumerate(FIELDS):
             try:
                 current[row[0]] = float(raw_data[i])
             except:
                 current[row[0]] = raw_data[i]
         return current
 
-    def historical(self,start=None, stop=None):
-        import datetime, time, urllib, math
-        start =  start or datetime.date(1900,1,1)
+    def historical(self, start=None, stop=None):
+        import datetime
+        import urllib
+        import math
+        start = start or datetime.date(1900, 1, 1)
         stop = stop or datetime.date.today()
         url = self.URL_HISTORICAL % dict(
             s=self.symbol,
-            a=start.month-1,b=start.day,c=start.year,
-            d=stop.month-1,e=stop.day,f=stop.year)
+            a=start.month - 1, b=start.day, c=start.year,
+            d=stop.month - 1, e=stop.day, f=stop.year)
         # Date,Open,High,Low,Close,Volume,Adj Close
         lines = urllib.urlopen(url).readlines()
-        raw_data = [row.split(',') for row in lines[1:] if row.count(',')==6]
+        raw_data = [row.split(',') for row in lines[1:] if row.count(',') == 6]
         previous_adjusted_close = 0
         series = []
         raw_data.reverse()
@@ -72,37 +108,38 @@ def historical(self,start=None, stop=None):
             adjusted_close = float(row[6])
             adjustment = adjusted_close/close
             if previous_adjusted_close:
-                arithmetic_return = adjusted_close/previous_adjusted_close-1.0
+                arithmetic_return = adjusted_close / previous_adjusted_close - 1.0
 
                 log_return = math.log(adjusted_close/previous_adjusted_close)
             else:
                 arithmetic_return = log_return = None
             previous_adjusted_close = adjusted_close
             series.append(dict(
-               date = datetime.datetime.strptime(row[0],'%Y-%m-%d'),
-               open = open,
-               high = high,
-               low = low,
-               close = close,
-               volume = vol,
-               adjusted_close = adjusted_close,
-               adjusted_open = open*adjustment,
-               adjusted_high = high*adjustment,
-               adjusted_low = low*adjustment,
-               adjusted_vol = vol/adjustment,
-               arithmetic_return = arithmetic_return,
-               log_return = log_return))
+                date=datetime.datetime.strptime(row[0], '%Y-%m-%d'),
+                open=open,
+                high=high,
+                low=low,
+                close=close,
+                volume=vol,
+                adjusted_close=adjusted_close,
+                adjusted_open=open*adjustment,
+                adjusted_high=high*adjustment,
+                adjusted_low=low*adjustment,
+                adjusted_vol=vol/adjustment,
+                arithmetic_return=arithmetic_return,
+                log_return=log_return))
         return series
 
     @staticmethod
-    def download(symbol='goog',what='adjusted_close',start=None,stop=None):
-        return [d[what] for d in YStock(symbol).historical(start,stop)]
+    def download(symbol='goog', what='adjusted_close', start=None, stop=None):
+        return [d[what] for d in YStock(symbol).historical(start, stop)]
 
 import os
 import uuid
 import sqlite3
 import cPickle as pickle
 
+
 class PersistentDictionary(object):
     """
     A sqlite based key,value storage.
@@ -134,8 +171,8 @@ def __init__(self,
         self.path = path
         self.autocommit = autocommit
         create_table = not os.path.exists(path)
-        self.connection  = sqlite3.connect(path)
-        self.connection.text_factory = str # do not use unicode
+        self.connection = sqlite3.connect(path)
+        self.connection.text_factory = str  # do not use unicode
         self.cursor = self.connection.cursor()
         if create_table:
             self.cursor.execute(self.CREATE_TABLE)
@@ -144,19 +181,19 @@ def __init__(self,
     def uuid(self):
         return str(uuid.uuid4())
 
-    def keys(self,pattern='*'):
-        "returns a list of keys filtered by a pattern, * is the wildcard"
-        self.cursor.execute(self.SELECT_KEYS,(pattern,))
+    def keys(self, pattern='*'):
+        """returns a list of keys filtered by a pattern, * is the wildcard"""
+        self.cursor.execute(self.SELECT_KEYS, (pattern, ))
         return [row[0] for row in self.cursor.fetchall()]
 
-    def __contains__(self,key):
+    def __contains__(self, key):
         return True if self[key] else False
 
     def __iter__(self):
         for key in self:
             yield key
 
-    def __setitem__(self,key,value):
+    def __setitem__(self, key, value):
         if value is None:
             del self[key]
             return
@@ -164,19 +201,19 @@ def __setitem__(self,key,value):
                             (key, pickle.dumps(value)))
         if self.autocommit: self.connection.commit()
 
-    def __getitem__(self,key):
+    def __getitem__(self, key):
         self.cursor.execute(self.SELECT_VALUE, (key,))
         row = self.cursor.fetchone()
         return pickle.loads(row[0]) if row else None
 
-    def __delitem__(self,pattern):
+    def __delitem__(self, pattern):
         self.cursor.execute(self.DELETE_KEY_VALUE, (pattern,))
         if self.autocommit: self.connection.commit()
 
-    def items(self,pattern='*'):
+    def items(self, pattern='*'):
         self.cursor.execute(self.SELECT_KEY_VALUE, (pattern,))
-        return [(row[0],pickle.loads(row[1])) \
-                    for row in self.cursor.fetchall()]
+        return [(row[0], pickle.loads(row[1]))
+                for row in self.cursor.fetchall()]
 
 import math
 import cmath
@@ -194,6 +231,7 @@ def items(self,pattern='*'):
 except ImportError:
     HAVE_MATPLOTLIB = False
 
+
 class Canvas(object):
 
     def __init__(self, title='', xlab='x', ylab='y', xrange=None, yrange=None):
@@ -226,36 +264,36 @@ def binary(self):
 
     def hist(self, data, bins=20, color='blue', legend=None):
         q = self.ax.hist(data, bins)
-        #if legend:
+        # if legend:
         #    self.legend.append((q[0], legend))
         return self
 
     def plot(self, data, color='blue', style='-', width=2,
              legend=None, xrange=None):
         if callable(data) and xrange:
-            x = [xrange[0]+0.01*i*(xrange[1]-xrange[0]) for i in xrange(0,101)]
+            x = [xrange[0] + 0.01*i*(xrange[1] - xrange[0]) for i in xrange(0, 101)]
             y = [data(p) for p in x]
-        elif data and isinstance(data[0],(int,float)):
+        elif data and isinstance(data[0], (int, float)):
             x, y = xrange(len(data)), data
         else:
             x, y = [p[0] for p in data], [p[1] for p in data]
         q = self.ax.plot(x, y, linestyle=style, linewidth=width, color=color)
         if legend:
-            self.legend.append((q[0],legend))
+            self.legend.append((q[0], legend))
         return self
 
     def errorbar(self, data, color='black', marker='o', width=2, legend=None):
-        x,y,dy = [p[0] for p in data], [p[1] for p in data], [p[2] for p in data]
+        x, y, dy = [p[0] for p in data], [p[1] for p in data], [p[2] for p in data]
         q = self.ax.errorbar(x, y, yerr=dy, fmt=marker, linewidth=width, color=color)
         if legend:
-            self.legend.append((q[0],legend))
+            self.legend.append((q[0], legend))
         return self
 
     def ellipses(self, data, color='blue', width=0.01, height=0.01, legend=None):
         for point in data:
             x, y = point[:2]
-            dx = point[2] if len(point)>2 else width
-            dy = point[3] if len(point)>3 else height
+            dx = point[2] if len(point) > 2 else width
+            dy = point[3] if len(point) > 3 else height
             ellipse = Ellipse(xy=(x, y), width=dx, height=dy)
             self.ax.add_artist(ellipse)
             ellipse.set_clip_box(self.ax.bbox)
@@ -269,11 +307,14 @@ def imshow(self, data, interpolation='bilinear'):
         self.ax.imshow(data).set_interpolation(interpolation)
         return self
 
+
 class memoize(object):
-    def __init__ (self, f):
+
+    def __init__(self, f):
         self.f = f
         self.storage = {}
-    def __call__ (self, *args, **kwargs):
+
+    def __call__(self, *args, **kwargs):
         key = str((self.f.__name__, args, kwargs))
         try:
             value = self.storage[key]
@@ -282,16 +323,20 @@ def __call__ (self, *args, **kwargs):
             self.storage[key] = value
         return value
 
+
 @memoize
 def fib(n):
-    return n if n<2 else fib(n-1)+fib(n-2)
+    return n if n < 2 else fib(n-1)+fib(n-2)
+
 
 class memoize_persistent(object):
     STORAGE = 'memoize.sqlite'
-    def __init__ (self, f):
+
+    def __init__(self, f):
         self.f = f
         self.storage = PersistentDictionary(memoize_persistent.STORAGE)
-    def __call__ (self, *args, **kwargs):
+
+    def __call__(self, *args, **kwargs):
         key = str((self.f.__name__, args, kwargs))
         if key in self.storage:
             value = self.storage[key]
@@ -300,16 +345,18 @@ def __call__ (self, *args, **kwargs):
             self.storage[key] = value
         return value
 
-def timef(f, ns=1000, dt = 60):
+
+def timef(f, ns=1000, dt=60):
     import time
     t = t0 = time.time()
-    for k in xrange(1,ns):
+    for k in xrange(1, ns):
         f()
         t = time.time()
-        if t-t0>dt: break
-    return (t-t0)/k
+        if t-t0 > dt: break
+    return (t-t0) / k
+
 
-def breadth_first_search(graph,start):
+def breadth_first_search(graph, start):
     vertices, link = graph
     blacknodes = []
     graynodes = [start]
@@ -319,12 +366,13 @@ def breadth_first_search(graph,start):
     while graynodes:
         current = graynodes.pop()
         for neighbor in neighbors[current]:
-            if not neighbor in blacknodes+graynodes:
-                graynodes.insert(0,neighbor)
+            if neighbor not in blacknodes + graynodes:
+                graynodes.insert(0, neighbor)
         blacknodes.append(current)
     return blacknodes
 
-def depth_first_search(graph,start):
+
+def depth_first_search(graph, start):
     vertices, link = graph
     blacknodes = []
     graynodes = [start]
@@ -334,232 +382,255 @@ def depth_first_search(graph,start):
     while graynodes:
         current = graynodes.pop()
         for neighbor in neighbors[current]:
-            if not neighbor in blacknodes+graynodes:
+            if neighbor not in blacknodes+graynodes:
                 graynodes.append(neighbor)
         blacknodes.append(current)
     return blacknodes
 
+
 class DisjointSets(object):
-    def __init__(self,n):
+
+    def __init__(self, n):
         self.sets = [-1]*n
         self.counter = n
-    def parent(self,i):
+
+    def parent(self, i):
         while True:
             j = self.sets[i]
-            if j<0:
+            if j < 0:
                 return i
             i = j
-    def join(self,i,j):
-        i,j = self.parent(i),self.parent(j)
-        if i!=j:
+
+    def join(self, i, j):
+        i, j = self.parent(i), self.parent(j)
+        if i != j:
             self.sets[i] += self.sets[j]
             self.sets[j] = i
-            self.counter-=1
-            return True # they have been joined
+            self.counter -= 1
+            return True  # they have been joined
         return False    # they were already joined
-    def joined(self,i,j):
-       return self.parent(i) == self.parent(j)
+
+    def joined(self, i, j):
+        return self.parent(i) == self.parent(j)
+
     def __len__(self):
         return self.counter
 
-def make_maze(n,d):
-    walls = [(i,i+n**j) for i in xrange(n**2) for j in xrange(d) if (i/n**j)%n+1<n]
+
+def make_maze(n, d):
+    walls = [(i, i+n**j) for i in xrange(n**2) for j in xrange(d) if (i/n**j) % n+1 < n]
     torn_down_walls = []
     ds = DisjointSets(n**d)
     random.shuffle(walls)
-    for i,wall in enumerate(walls):
-        if ds.join(wall[0],wall[1]):
+    for i, wall in enumerate(walls):
+        if ds.join(wall[0], wall[1]):
             torn_down_walls.append(wall)
-        if len(ds)==1:
+        if len(ds) == 1:
             break
-    walls = [wall for wall in walls if not wall in torn_down_walls]
+    walls = [wall for wall in walls if wall not in torn_down_walls]
     return walls, torn_down_walls
 
+
 def Kruskal(graph):
     vertices, links = graph
     A = []
     S = DisjointSets(len(vertices))
-    links.sort(cmp=lambda a,b: cmp(a[2],b[2]))
-    for source,dest,length in links:
-        if S.join(source,dest):
-            A.append((source,dest,length))
+    links.sort(cmp=lambda a, b: cmp(a[2], b[2]))
+    for source, dest, length in links:
+        if S.join(source, dest):
+            A.append((source, dest, length))
     return A
 
+
 class PrimVertex(object):
     INFINITY = 1e100
-    def __init__(self,id,links):
+
+    def __init__(self, id, links):
         self.id = id
         self.closest = None
         self.closest_dist = PrimVertex.INFINITY
-        self.neighbors = [link[1:] for link in links if link[0]==id]
-    def __cmp__(self,other):
+        self.neighbors = [link[1:] for link in links if link[0] == id]
+
+    def __cmp__(self, other):
         return cmp(self.closest_dist, other.closest_dist)
 
+
 def Prim(graph, start):
     from heapq import heappush, heappop, heapify
     vertices, links = graph
-    P = [PrimVertex(i,links) for i in vertices]
-    Q = [P[i] for i in vertices if not i==start]
+    P = [PrimVertex(i, links) for i in vertices]
+    Q = [P[i] for i in vertices if not i == start]
     vertex = P[start]
     while Q:
-        for neighbor_id,length in vertex.neighbors:
+        for neighbor_id, length in vertex.neighbors:
             neighbor = P[neighbor_id]
-            if neighbor in Q and length<neighbor.closest_dist:
-                 neighbor.closest = vertex
-                 neighbor.closest_dist = length
+            if neighbor in Q and length < neighbor.closest_dist:
+                neighbor.closest = vertex
+                neighbor.closest_dist = length
         heapify(Q)
         vertex = heappop(Q)
-    return [(v.id,v.closest.id,v.closest_dist) for v in P if not v.id==start]
+    return [(v.id, v.closest.id, v.closest_dist) for v in P if not v.id == start]
+
 
 def Dijkstra(graph, start):
     from heapq import heappush, heappop, heapify
     vertices, links = graph
-    P = [PrimVertex(i,links) for i in vertices]
-    Q = [P[i] for i in vertices if not i==start]
+    P = [PrimVertex(i, links) for i in vertices]
+    Q = [P[i] for i in vertices if not i == start]
     vertex = P[start]
     vertex.closest_dist = 0
     while Q:
-        for neighbor_id,length in vertex.neighbors:
+        for neighbor_id, length in vertex.neighbors:
             neighbor = P[neighbor_id]
-            dist = length+vertex.closest_dist
-            if neighbor in Q and dist<neighbor.closest_dist:
-                 neighbor.closest = vertex
-                 neighbor.closest_dist = dist
+            dist = length + vertex.closest_dist
+            if neighbor in Q and dist < neighbor.closest_dist:
+                neighbor.closest = vertex
+                neighbor.closest_dist = dist
         heapify(Q)
         vertex = heappop(Q)
-    return [(v.id,v.closest.id,v.closest_dist) for v in P if not v.id==start]
+    return [(v.id, v.closest.id, v.closest_dist) for v in P if not v.id == start]
+
 
 def encode_huffman(input):
     from heapq import heappush, heappop
 
     def inorder_tree_walk(t, key, keys):
-        (f,ab) = t
-        if isinstance(ab,tuple):
-            inorder_tree_walk(ab[0],key+'0',keys)
-            inorder_tree_walk(ab[1],key+'1',keys)
+        (f, ab) = t
+        if isinstance(ab, tuple):
+            inorder_tree_walk(ab[0], key+'0', keys)
+            inorder_tree_walk(ab[1], key+'1', keys)
         else:
             keys[ab] = key
 
     symbols = {}
     for symbol in input:
-        symbols[symbol] = symbols.get(symbol,0)+1
+        symbols[symbol] = symbols.get(symbol, 0)+1
     heap = []
-    for (k,f) in symbols.items():
-        heappush(heap,(f,k))
-    while len(heap)>1:
-        (f1,k1) = heappop(heap)
-        (f2,k2) = heappop(heap)
-        heappush(heap,(f1+f2,((f1,k1),(f2,k2))))
+    for (k, f) in symbols.items():
+        heappush(heap, (f, k))
+    while len(heap) > 1:
+        (f1, k1) = heappop(heap)
+        (f2, k2) = heappop(heap)
+        heappush(heap, (f1 + f2, ((f1, k1), (f2, k2))))
     symbol_map = {}
-    inorder_tree_walk(heap[0],'',symbol_map)
+    inorder_tree_walk(heap[0], '', symbol_map)
     encoded = ''.join(symbol_map[symbol] for symbol in input)
     return symbol_map, encoded
 
+
 def decode_huffman(keys, encoded):
-    reversed_map = dict((v,k) for (k,v) in keys.items())
+    reversed_map = dict((v, k) for (k, v) in keys.items())
     i, output = 0, []
-    for j in xrange(1,len(encoded)+1):
+    for j in xrange(1, len(encoded)+1):
         if encoded[i:j] in reversed_map:
-           output.append(reversed_map[encoded[i:j]])
-           i=j
+            output.append(reversed_map[encoded[i:j]])
+            i = j
     return ''.join(output)
 
+
 def lcs(a, b):
     previous = [0]*len(a)
-    for i,r in enumerate(a):
+    for i, r in enumerate(a):
         current = []
-        for j,c in enumerate(b):
-            if r==c:
-                e = previous[j-1]+1 if i*j>0 else 1
+        for j, c in enumerate(b):
+            if r == c:
+                e = previous[j-1]+1 if i*j > 0 else 1
             else:
-                e = max(previous[j] if i>0 else 0,
-                        current[-1] if j>0 else 0)
+                e = max(previous[j] if i > 0 else 0,
+                        current[-1] if j > 0 else 0)
             current.append(e)
-        previous=current
+        previous = current
     return current[-1]
 
-def needleman_wunsch(a,b,p=0.97):
-    z=[]
-    for i,r in enumerate(a):
+
+def needleman_wunsch(a, b, p=0.97):
+    z = []
+    for i, r in enumerate(a):
         z.append([])
-        for j,c in enumerate(b):
-            if r==c:
-                e = z[i-1][j-1]+1 if i*j>0 else 1
+        for j, c in enumerate(b):
+            if r == c:
+                e = z[i-1][j-1]+1 if i*j > 0 else 1
             else:
-                e = p*max(z[i-1][j] if i>0 else 0,
-                          z[i][j-1] if j>0 else 0)
+                e = p*max(z[i-1][j] if i > 0 else 0,
+                          z[i][j-1] if j > 0 else 0)
             z[-1].append(e)
     return z
 
-def continuum_knapsack(a,b,c):
-    table = [(a[i]/b[i],i) for i in xrange(len(a))]
+
+def continuum_knapsack(a, b, c):
+    table = [(a[i]/b[i], i) for i in xrange(len(a))]
     table.sort()
     table.reverse()
-    f=0.0
-    for (y,i) in table:
-        quantity = min(c/b[i],1)
-        x.append((i,quantity))
+    f = 0.0
+    for (y, i) in table:
+        quantity = min(c / b[i], 1)
+        x.append((i, quantity))
         c = c-b[i]*quantity
         f = f+a[i]*quantity
-    return (f,x)
+    return (f, x)
+
 
 class Cluster(object):
-    def __init__(self,points,metric,weights=None):
+
+    def __init__(self, points, metric, weights=None):
         self.points, self.metric = points, metric
         self.k = len(points)
         self.w = weights or [1.0]*self.k
-        self.q = dict((i,[i]) for i,e in enumerate(points))
+        self.q = dict((i, [i]) for i, e in enumerate(points))
         self.d = []
         for i in xrange(self.k):
-            for j in xrange(i+1,self.k):
-                m = metric(points[i],points[j])
-                if not m is None:
-                    self.d.append((m,i,j))
+            for j in xrange(i+1, self.k):
+                m = metric(points[i], points[j])
+                if m is not None:
+                    self.d.append((m, i, j))
         self.d.sort()
         self.dd = []
-    def parent(self,i):
-        while isinstance(i,int): (parent, i) = (i, self.q[i])
+
+    def parent(self, i):
+        while isinstance(i, int): (parent, i) = (i, self.q[i])
         return parent, i
+
     def step(self):
-        if self.k>1:
+        if self.k > 1:
             # find new clusters to join
-            (self.r,i,j),self.d = self.d[0],self.d[1:]
+            (self.r, i, j), self.d = self.d[0], self.d[1:]
             # join them
-            i,x = self.parent(i) # find members of cluster i
-            j,y = self.parent(j) # find members if cluster j
+            i, x = self.parent(i)  # find members of cluster i
+            j, y = self.parent(j)  # find members if cluster j
             x += y               # join members
             self.q[j] = i        # make j cluster point to i
             self.k -= 1          # decrease cluster count
             # update all distances to new joined cluster
-            new_d = [] # links not related to joined clusters
-            old_d = {} # old links related to joined clusters
-            for (r,h,k) in self.d:
-                if h in (i,j):
-                    a,b = old_d.get(k,(0.0,0.0))
-                    old_d[k] = a+self.w[k]*r,b+self.w[k]
-                elif k in (i,j):
-                    a,b = old_d.get(h,(0.0,0.0))
-                    old_d[h] = a+self.w[h]*r,b+self.w[h]
+            new_d = []  # links not related to joined clusters
+            old_d = {}  # old links related to joined clusters
+            for (r, h, k) in self.d:
+                if h in (i, j):
+                    a, b = old_d.get(k, (0.0, 0.0))
+                    old_d[k] = a+self.w[k]*r, b+self.w[k]
+                elif k in (i, j):
+                    a, b = old_d.get(h, (0.0, 0.0))
+                    old_d[h] = a + self.w[h]*r, b+self.w[h]
                 else:
-                    new_d.append((r,h,k))
-            new_d += [(a/b,i,k) for k,(a,b) in old_d.items()]
+                    new_d.append((r, h, k))
+            new_d += [(a/b, i, k) for k, (a, b) in old_d.items()]
             new_d.sort()
             self.d = new_d
             # update weight of new cluster
             self.w[i] = self.w[i]+self.w[j]
             # get new list of cluster members
-            self.v = [s for s in self.q.values() if isinstance(s,list)]
-            self.dd.append((self.r,len(self.v)))
+            self.v = [s for s in self.q.values() if isinstance(s, list)]
+            self.dd.append((self.r, len(self.v)))
         return self.r, self.v
 
-    def find(self,k):
+    def find(self, k):
         # if necessary start again
-        if self.k<k: self.__init__(self.points,self.metric)
+        if self.k < k: self.__init__(self.points, self.metric)
         # step until we get k clusters
-        while self.k>k: self.step()
+        while self.k > k: self.step()
         # return list of cluster members
         return self.r, self.v
 
+
 class NeuralNetwork:
     """
     Back-Propagation Neural Networks
@@ -572,7 +643,7 @@ class NeuralNetwork:
     @staticmethod
     def rand(a, b):
         """ calculate a random number where:  a <= rand < b """
-        return (b-a)*random.random() + a
+        return (b - a) * random.random() + a
 
     @staticmethod
     def sigmoid(x):
@@ -613,7 +684,7 @@ def update(self, inputs):
 
         # hidden activations
         for j in xrange(self.nh):
-            s = sum(self.ai[i] * self.wi[i,j] for i in xrange(self.ni))
+            s = sum(self.ai[i] * self.wi[i, j] for i in xrange(self.ni))
             self.ah[j] = self.sigmoid(s)
 
         # output activations
@@ -650,8 +721,8 @@ def back_propagate(self, targets, N, M):
         for i in xrange(self.ni):
             for j in xrange(self.nh):
                 change = hidden_deltas[j]*self.ai[i]
-                self.wi[i,j] = self.wi[i,j] + N*change + M*self.ci[i,j]
-                self.ci[i,j] = change
+                self.wi[i, j] = self.wi[i, j] + N*change + M*self.ci[i, j]
+                self.ci[i, j] = change
 
         # calculate error
         error = sum(0.5*(targets[k]-self.ao[k])**2 for k in xrange(len(targets)))
@@ -683,24 +754,28 @@ def train(self, patterns, iterations=1000, N=0.5, M=0.1, check=False):
             if check and i % 100 == 0:
                 print 'error %-14f' % error
 
+
 def D(f,h=1e-6): # first derivative of f
     return lambda x,f=f,h=h: (f(x+h)-f(x-h))/2/h
 
+
 def DD(f,h=1e-6): # second derivative of f
     return lambda x,f=f,h=h: (f(x+h)-2.0*f(x)+f(x-h))/(h*h)
 
+
 def myexp(x,precision=1e-6,max_steps=40):
-    if x==0:
-       return 1.0
-    elif x>0:
-       return 1.0/myexp(-x,precision,max_steps)
+    if x == 0:
+        return 1.0
+    elif x > 0:
+        return 1.0/myexp(-x,precision,max_steps)
     else:
-       t = s = 1.0 # first term
-       for k in xrange(1,max_steps):
+        t = s = 1.0 # first term
+        for k in xrange(1,max_steps):
            t = t*x/k   # next term
            s = s + t   # add next term
            if abs(t)<precision: return s
-       raise ArithmeticError('no convergence')
+        raise ArithmeticError('no convergence')
+
 
 def mysin(x,precision=1e-6,max_steps=40):
     pi = math.pi
@@ -725,6 +800,7 @@ def mysin(x,precision=1e-6,max_steps=40):
            if r<precision: return s  # stopping condition
        raise ArithmeticError('no convergence')
 
+
 def mycos(x,precision=1e-6,max_steps=40):
     pi = math.pi
     if x==0:
@@ -748,6 +824,7 @@ def mycos(x,precision=1e-6,max_steps=40):
            if r<precision: return s  # stopping condition
        raise ArithmeticError('no convergence')
 
+
 class Matrix(object):
     def __init__(self,rows,cols=1,fill=0.0):
         """
@@ -776,17 +853,17 @@ def __init__(self,rows,cols=1,fill=0.0):
         self.ncols = len(self.rows[0])
 
     def __getitem__(A,coords):
-        " x = A[0,1]"
-        i,j = coords
+        """ x = A[0,1]"""
+        i, j = coords
         return A.rows[i][j]
 
     def __setitem__(A,coords,value):
-        " A[0,1] = 3.0 "
-        i,j = coords
+        """ A[0,1] = 3.0 """
+        i, j = coords
         A.rows[i][j] = value
 
     def tolist(A):
-        " assert(Matrix([[1,2],[3,4]]).tolist() == [[1,2],[3,4]]) "
+        """ assert(Matrix([[1,2],[3,4]]).tolist() == [[1,2],[3,4]]) """
         return A.rows
 
     def __str__(A):
@@ -801,9 +878,9 @@ def reshape(A,n,m):
         if n*m != A.nrows*A.ncols:
              raise RuntimeError("Impossible reshape")
         flat = A.flatten()
-        return Matrix(n,m,fill=lambda r,c,m=m,flat=flat: flat[r*m+c])
+        return Matrix(n,m, fill=lambda r, c, m=m, flat=flat: flat[r*m+c])
 
-    def swap_rows(A,i,j):
+    def swap_rows(A,i, j):
         " assert(Matrix([[1,2],[3,4]]).swap_rows(1,0).tolist() == [[3,4],[1,2]]) "
         A.rows[i], A.rows[j] = A.rows[j], A.rows[i]
 
@@ -862,10 +939,13 @@ def __sub__(A,B):
             for c in xrange(m):
                 C[r,c] = A[r,c]-B[r,c]
         return C
+
     def __radd__(A,B): #B+A
         return A+B
+
     def __rsub__(A,B): #B-A
         return (-A)+B
+
     def __neg__(A):
         return Matrix(A.nrows,A.ncols,fill=lambda r,c:-A[r,c])
 
@@ -928,7 +1008,6 @@ def __div__(A,B):
         else:
             return (1.0/B)*A # matrix/scalar
 
-
     @property
     def T(A):
         """Transposed of A"""
@@ -939,6 +1018,7 @@ def H(A):
         """Hermitian of A"""
         return Matrix(A.ncols,A.nrows, fill=lambda r,c: A[c,r].conj())
 
+
 def is_almost_symmetric(A, ap=1e-6, rp=1e-4):
     if A.nrows != A.ncols: return False
     for r in xrange(A.nrows):
@@ -948,6 +1028,7 @@ def is_almost_symmetric(A, ap=1e-6, rp=1e-4):
                 return False
     return True
 
+
 def is_almost_zero(A, ap=1e-6, rp=1e-4):
     for r in xrange(A.nrows):
         for c in xrange(A.ncols):
@@ -956,6 +1037,7 @@ def is_almost_zero(A, ap=1e-6, rp=1e-4):
                 return False
     return True
 
+
 def norm(A,p=1):
     if isinstance(A,(list,tuple)):
         return sum(abs(x)**p for x in A)**(1.0/p)
@@ -973,6 +1055,7 @@ def norm(A,p=1):
     else:
         return abs(A)
 
+
 def condition_number(f,x=None,h=1e-6):
     if callable(f) and not x is None:
         return D(f,h)(x)*x/f(x)
@@ -981,6 +1064,7 @@ def condition_number(f,x=None,h=1e-6):
     else:
         raise NotImplementedError
 
+
 def exp(x,ap=1e-6,rp=1e-4,ns=40):
     if isinstance(x,Matrix):
        t = s = Matrix.identity(x.ncols)
@@ -994,8 +1078,10 @@ def exp(x,ap=1e-6,rp=1e-4,ns=40):
     else:
        return math.exp(x)
 
+
 def Cholesky(A):
-    import copy, math
+    import copy
+    import math
     if not is_almost_symmetric(A):
         raise ArithmeticError('not symmetric')
     L = copy.deepcopy(A)
@@ -1008,10 +1094,10 @@ def Cholesky(A):
         for j in xrange(k+1,L.nrows):
             p=float(L[j,k])
             for i in xrange(k+1,L.nrows):
-                L[i,j] -= p*L[i,k]
+                L[i, j] -= p*L[i,k]
     for  i in xrange(L.nrows):
         for j in xrange(i+1,L.ncols):
-            L[i,j]=0
+            L[i, j]=0
     return L
 
 def is_positive_definite(A):
@@ -1068,7 +1154,7 @@ def eval_fitting_function(f,c,x):
         weight = 1.0/points[i][2] if len(points[i])>2 else 1.0
         b[i,0] = weight*float(points[i][1])
         for j in xrange(A.ncols):
-            A[i,j] = weight*f[j](float(points[i][0]))
+            A[i, j] = weight*f[j](float(points[i][0]))
     c = (1.0/(A.T*A))*(A.T*b)
     chi = A*c-b
     chi2 = norm(chi,2)**2
@@ -1200,8 +1286,8 @@ def maxind(M,k):
     # normalize vectors
     U = Matrix(n,n)
     for i in indexes:
-        norm = sqrt(sum(E[i,j]**2 for j in indexes))
-        for j in indexes: U[j,i] = E[i,j]/norm
+        norm = sqrt(sum(E[i, j]**2 for j in indexes))
+        for j in indexes: U[j,i] = E[i, j]/norm
     return U,e
 
 
@@ -1223,7 +1309,7 @@ def compute_correlation(stocks, key='arithmetic_return'):
     mus = [sum(v[i][k] for k in xrange(1,n))/n for i in iter_stocks]
     # fill in the covariance matrix
     var = [sum(v[i][k]**2 for k in xrange(1,n))/n - mus[i]**2 for i in iter_stocks]
-    corr = Matrix(nstocks,nstocks,fill=lambda i,j: \
+    corr = Matrix(nstocks,nstocks,fill=lambda i, j: \
              (sum(v[i][k]*v[j][k] for k in xrange(1,n))/n - mus[i]*mus[j])/ \
              math.sqrt(var[i]*var[j]))
     return corr
@@ -1313,7 +1399,7 @@ def solve_secant(f, x, ap=1e-6, rp=1e-4, ns=20):
     raise ArithmeticError('no convergence')
 
 def optimize_bisection(f, a, b, ap=1e-6, rp=1e-4, ns=100):
-    return solve_bisection(D(f), a, b , ap, rp, ns)
+    return solve_bisection(D(f), a, b, ap, rp, ns)
 
 def optimize_newton(f, x, ap=1e-6, rp=1e-4, ns=20):
     x = float(x) # make sure it is not int
@@ -1465,6 +1551,7 @@ def g(b, data=data, f=fs, constraint=constraint):
         return a, chi2
     else:
         na = len(fs)
+
         def core(b,data=data,fs=fs):
             A = Matrix([[fs[k](b,x)/dy for k in xrange(na)] \
                                   for (x,y,dy) in data])
@@ -1472,6 +1559,7 @@ def core(b,data=data,fs=fs):
             a = (1/(A.T*A))*(A.T*z)
             chi2 = norm(A*a-z)**2
             return a.flatten(), chi2
+
         def g(b,data=data,fs=fs,constraint=constraint):
             a, chi2 = core(b, data, fs)
             if constraint:
@@ -1508,12 +1596,14 @@ class QuadratureIntegrator:
     Calculates the integral of the function f from points a to b
     using n Vandermonde weights and numerical quadrature.
     """
+
     def __init__(self,delta,order=4):
         h = float(delta)/(order-1)
         A = Matrix(order, order, fill = lambda r,c: (c*h)**r)
         s = Matrix(order, 1, fill = lambda r,c: (delta**(r+1))/(r+1))
         w = (1/A)*s
         self.packed = (h, order, w)
+
     def integrate(self,f,a):
         (h, order, w) = self.packed
         return sum(w[i,0]*f(a+i*h) for i in xrange(order))
@@ -1535,12 +1625,15 @@ def correlation(X,Y):
     return covariance(X,Y)/sd(X)/sd(Y)
 
 class MCG(object):
+
     def __init__(self,seed,a=66539,m=2**31):
         self.x = seed
         self.a, self.m = a, m
+
     def next(self):
         self.x = (self.a*self.x) % self.m
         return self.x
+
     def random(self):
         return float(self.next())/self.m
 
@@ -1579,54 +1672,58 @@ def random(self):
         y ^= (y << 7) & 0x9d2c5680
         y ^= (y << 15) & 0xefc60000
         y ^= (y >> 18)
-        return (float(y)/0xffffffff )
+        return (float(y)/0xffffffff)
+
 
-def leapfrog(mcg,k):
+def leapfrog(mcg, k):
     a = mcg.a**k % mcg.m
-    return [MCG(mcg.next(),a,mcg.m) for i in range(k)]
+    return [MCG(mcg.next(), a, mcg.m) for i in range(k)]
+
 
 class RandomSource(object):
-    def __init__(self,generator=None):
+    def __init__(self, generator=None):
         if not generator:
             import random as generator
         self.generator = generator
+
     def random(self):
         return self.generator.random()
-    def randint(self,a,b):
-        return int(a+(b-a+1)*self.random())
 
-    def choice(self,S):
-        return S[self.randint(0,len(S)-1)]
+    def randint(self, a, b):
+        return int(a + (b - a + 1) * self.random())
 
-    def bernoulli(self,p):
-        return 1 if self.random()<p else 0
+    def choice(self, S):
+        return S[self.randint(0, len(S) - 1)]
 
-    def lookup(self,table, epsilon=1e-6):
-        if isinstance(table,dict): table = table.items()
+    def bernoulli(self, p):
+        return 1 if self.random() < p else 0
+
+    def lookup(self, table, epsilon=1e-6):
+        if isinstance(table, dict): table = table.items()
         u = self.random()
-        for key,p in table:
-            if u<p+epsilon:
+        for key, p in table:
+            if u < p+epsilon:
                 return key
             u = u - p
         raise ArithmeticError('invalid probability')
 
-    def binomial(self,n,p,epsilon=1e-6):
+    def binomial(self, n, p, epsilon=1e-6):
         u = self.random()
         q = (1.0-p)**n
         for k in xrange(n+1):
-            if u<q+epsilon:
+            if u < q+epsilon:
                 return k
             else:
                 u = u - q
             q = q*(n-k)/(k+1)*p/(1.0-p)
         raise ArithmeticError('invalid probability')
 
-    def poisson(self,lamb,epsilon=1e-6):
+    def poisson(self, lamb, epsilon=1e-6):
         u = self.random()
         q = exp(-lamb)
-        k=0
+        k = 0
         while True:
-            if u<q+epsilon:
+            if u < q+epsilon:
                 return k
             else:
                 u = u - q
@@ -1634,73 +1731,76 @@ def poisson(self,lamb,epsilon=1e-6):
             k = k+1
         raise ArithmeticError('invalid probability')
 
-    def uniform(self,a,b):
+    def uniform(self, a, b):
         return a+(b-a)*self.random()
 
-    def exponential(self,lamb):
+    def exponential(self, lamb):
         return -log(self.random())/lamb
 
-    def gauss(self,mu=0.0,sigma=1.0):
-        if hasattr(self,'other') and self.other:
+    def gauss(self, mu=0.0, sigma=1.0):
+        if hasattr(self, 'other') and self.other:
             this, other = self.other, None
         else:
             while True:
-                v1 = self.random(-1,1)
-                v2 = self.random(-1,1)
+                v1 = self.random(-1, 1)
+                v2 = self.random(-1, 1)
                 r = v1*v1+v2*v2
-                if r<1: break
+                if r < 1: break
             this = sqrt(-2.0*log(r)/r)*v1
             self.other = sqrt(-2.0*log(r)/r)*v1
         return mu+sigma*this
 
-def confidence_intervals(mu,sigma):
+
+def confidence_intervals(mu, sigma):
     """Computes the normal confidence intervals"""
-    CONFIDENCE=[
-        (0.68,1.0),
-        (0.80,1.281551565545),
-        (0.90,1.644853626951),
-        (0.95,1.959963984540),
-        (0.98,2.326347874041),
-        (0.99,2.575829303549),
-        (0.995,2.807033768344),
-        (0.998,3.090232306168),
-        (0.999,3.290526731492),
-        (0.9999,3.890591886413),
-        (0.99999,4.417173413469)
+    CONFIDENCE = [
+        (0.68, 1.0),
+        (0.80, 1.281551565545),
+        (0.90, 1.644853626951),
+        (0.95, 1.959963984540),
+        (0.98, 2.326347874041),
+        (0.99, 2.575829303549),
+        (0.995, 2.807033768344),
+        (0.998, 3.090232306168),
+        (0.999, 3.290526731492),
+        (0.9999, 3.890591886413),
+        (0.99999, 4.417173413469)
         ]
-    return [(a,mu-b*sigma,mu+b*sigma) for (a,b) in CONFIDENCE]
+    return [(a, mu-b*sigma, mu+b*sigma) for (a, b) in CONFIDENCE]
 
-    def pareto(self,alpha,xm):
+    def pareto(self, alpha, xm):
         u = self.random()
         return xm*(1.0-u)**(-1.0/alpha)
 
     def point_on_circle(self, radius=1.0):
-        angle = 2.0*pi*self.random()
+        angle = 2.0 * pi * self.random()
         return radius*math.cos(angle), radius*math.sin(angle)
 
-    def point_in_circle(self,radius=1.0):
+    def point_in_circle(self, radius=1.0):
         while True:
-            x = self.uniform(-radius,radius)
-            y = self.uniform(-radius,radius)
-            if x*x+y*y < radius*radius:
-                return x,y
+            x = self.uniform(-radius, radius)
+            y = self.uniform(-radius, radius)
+            if x * x + y * y < radius*radius:
+                return x, y
 
-    def point_in_sphere(self,radius=1.0):
+    def point_in_sphere(self, radius=1.0):
         while True:
-            x = self.uniform(-radius,radius)
-            y = self.uniform(-radius,radius)
-            z = self.uniform(-radius,radius)
+            x = self.uniform(-radius, radius)
+            y = self.uniform(-radius, radius)
+            z = self.uniform(-radius, radius)
             if x*x+y*y*z*z < radius*radius:
-                return x,y,z
+                return x, y, z
 
     def point_on_sphere(self, radius=1.0):
-        x,y,z = self.point_in_sphere(radius)
+        x, y, z = self.point_in_sphere(radius)
         norm = math.sqrt(x*x+y*y+z*z)
-        return x/norm,y/norm,z/norm
+        return x / norm, y / norm, z / norm
 
-def resample(S,size=None):
+
+def resample(S, size=None):
     return [random.choice(S) for i in xrange(size or len(S))]
 
+
 def bootstrap(x, confidence=0.68, nsamples=100):
     """Computes the bootstrap errors of the input list."""
     def mean(S): return float(sum(x for x in S))/len(S)
@@ -1710,6 +1810,7 @@ def mean(S): return float(sum(x for x in S))/len(S)
     right_tail = nsamples-1-left_tail
     return means[left_tail], mean(x), means[right_tail]
 
+
 class MCEngine:
     """
     Monte Carlo Engine parent class.
@@ -1723,16 +1824,16 @@ def simulate_many(self, ap=0.1, rp=0.1, ns=1000):
         self.results = []
         s1=s2=0.0
         self.convergence=False
-        for k in xrange(1,ns):
+        for k in xrange(1, ns):
             x = self.simulate_once()
             self.results.append(x)
             s1 += x
             s2 += x*x
-            mu = float(s1)/k
-            variance = float(s2)/k-mu*mu
-            dmu = sqrt(variance/k)
-            if k>10:
-                if abs(dmu)<max(ap,abs(mu)*rp):
+            mu = float(s1) / k
+            variance = float(s2) / k-mu*mu
+            dmu = sqrt(variance / k)
+            if k > 10:
+                if abs(dmu) < max(ap, abs(mu) * rp):
                     self.converence = True
                     break
         self.results.sort()
@@ -1764,7 +1865,7 @@ def test071():
     ...     if not key in storage:
     ...         storage[key] = YStock(symbol).historical(start=date(2011,1,1),
     ...                                                  stop=date(2011,12,31))
-    
+
     """
     pass
 
@@ -1775,7 +1876,7 @@ def test072():
     >>> appl = storage['AAPL/2011']
     >>> points = [(x,y['adjusted_close']) for (x,y) in enumerate(appl)]
     >>> Canvas(title='Apple Stock (2011)',xlab='trading day',ylab='adjusted close').plot(points,legend='AAPL').save('images/aapl2011.png')
-    
+
     """
     pass
 
@@ -1786,7 +1887,7 @@ def test073():
     >>> appl = storage['AAPL/2011'][1:] # skip 1st day
     >>> points = [day['arithmetic_return'] for day in appl]
     >>> Canvas(title='Apple Stock (2011)',xlab='arithmetic return', ylab='frequency').hist(points).save('images/aapl2011hist.png')
-    
+
     """
     pass
 
@@ -1796,7 +1897,7 @@ def test074():
     >>> from random import gauss
     >>> points = [(gauss(0,1),gauss(0,1),gauss(0,0.2),gauss(0,0.2)) for i in xrange(30)]
     >>> Canvas(title='example scatter plot', xrange=(-2,2), yrange=(-2,2)).ellipses(points).save('images/scatter.png')
-    
+
     """
     pass
 
@@ -1814,7 +1915,7 @@ def test075():
     ...      xrange = (min(p[0] for p in points),max(p[0] for p in points)),
     ...      yrange = (min(p[1] for p in points),max(p[1] for p in points))
     ...      ).ellipses(points).save('images/sp100rr.png')
-    
+
     """
     pass
 
@@ -1824,7 +1925,7 @@ def test076():
     >>> def f(x,y): return (x-1)**2+(y-2)**2
     >>> points = [[f(0.1*i-3,0.1*j-3) for i in range(61)] for j in range(61)]
     >>> Canvas(title='example 2d function').imshow(points).save('images/color2d.png')
-    
+
     """
     pass
 
@@ -1833,7 +1934,7 @@ def test077():
     """
     >>> print fib(11)
     89
-    
+
     """
     pass
 
@@ -1841,7 +1942,7 @@ def test077():
 def test078():
     """
     >>> walls, torn_down_walls = make_maze(n=20,d=2)
-    
+
     """
     pass
 
@@ -1849,7 +1950,7 @@ def test078():
 def test079():
     """
     >>> vertices = xrange(10)
-    >>> links = [(i,j,abs(math.sin(i+j+1))) for i in vertices for j in vertices]
+    >>> links = [(i, j,abs(math.sin(i+j+1))) for i in vertices for j in vertices]
     >>> graph = [vertices,links]
     >>> links = Dijkstra(graph,0)
     >>> for link in links: print link
@@ -1862,7 +1963,7 @@ def test079():
     (7, 2, 0.685...)
     (8, 0, 0.412...)
     (9, 0, 0.544...)
-    
+
     """
     pass
 
@@ -1871,11 +1972,11 @@ def test080():
     """
     >>> n,d = 4, 2
     >>> walls, links = make_maze(n,d)
-    >>> symmetrized_links = [(i,j,1) for (i,j) in links]+[(j,i,1) for (i,j) in links]
+    >>> symmetrized_links = [(i, j,1) for (i, j) in links]+[(j,i,1) for (i, j) in links]
     >>> graph = [xrange(n*n),symmetrized_links]
     >>> links = Dijkstra(graph,0)
-    >>> paths = dict((i,(j,d)) for (i,j,d) in links)
-    
+    >>> paths = dict((i,(j,d)) for (i, j,d) in links)
+
     """
     pass
 
@@ -1891,7 +1992,7 @@ def test081():
     True
     >>> print 1.0*len(input)/(len(encoded)/8)
     2.57...
-    
+
     """
     pass
 
@@ -1904,7 +2005,7 @@ def test082():
     >>> E = -sum(wi*log(wi,2) for wi in w)
     >>> print E
     3.23...
-    
+
     """
     pass
 
@@ -1915,7 +2016,7 @@ def test083():
     >>> dna2 = 'GATAGGTACCACAATAATAAGGATAGCTCGCAAATCCTCGA'
     >>> print lcs(dna1,dna2)
     26
-    
+
     """
     pass
 
@@ -1929,7 +2030,7 @@ def test084():
     >>> chromosome2 = ''.join(choice(genes) for i in xrange(10))
     >>> z = needleman_wunsch(chromosome1, chromosome2)
     >>> Canvas(title='Needleman-Wunsch').imshow(z).save('images/needleman.png')
-    
+
     """
     pass
 
@@ -1945,7 +2046,7 @@ def test085():
     ...       ).plot(c.dd[150:]).save('clustering1.png')
     >>> Canvas(title='clustering example (2d projection)',xlab='p[0]',ylab='p[1]'
     ...       ).ellipses([p[:2] for p in points]).save('clustering2.png')
-    
+
     """
     pass
 
@@ -1960,7 +2061,7 @@ def test086():
     [0, 1] -> [0.98...]
     [1, 0] -> [0.98...]
     [1, 1] -> [-0.00...]
-    
+
     """
     pass
 
@@ -1979,7 +2080,7 @@ def test088():
     >>> f2 = D(f1) # second derivative
     >>> print f2(0)
     1.99999...
-    
+
     """
     pass
 
@@ -1997,7 +2098,7 @@ def test089():
     >>> c.plot([(x,x-x**3/6+x**5/120) for x in X[:100]],legend='Taylor 5th')
     <...>
     >>> c.save('images/sin.png')
-    
+
     """
     pass
 
@@ -2016,7 +2117,7 @@ def test090():
     >>> c.plot([(x,1-(x-a)**2/2+(x-a)**4/24-(x-a)**6/720) for x in X[:150]],legend='Taylor 6th')
     <...>
     >>> c.save('images/sin2.png')
-    
+
     """
     pass
 
@@ -2026,7 +2127,7 @@ def test091():
     >>> for i in xrange(10):
     ...     x= 0.1*i
     ...     assert abs(myexp(x) - math.exp(x)) < 1e-4
-    
+
     """
     pass
 
@@ -2036,7 +2137,7 @@ def test092():
     >>> for i in xrange(10):
     ...     x= 0.1*i
     ...     assert abs(mysin(x) - math.sin(x)) < 1e-4
-    
+
     """
     pass
 
@@ -2046,7 +2147,7 @@ def test093():
     >>> for i in xrange(10):
     ...     x = 0.1*i
     ...     assert abs(mycos(x) - math.cos(x)) < 1e-4
-    
+
     """
     pass
 
@@ -2067,7 +2168,7 @@ def test094():
     >>> b = Matrix([[1.0],[2.0],[3.0]])
     >>> print b + 2      # calls b.__add__(2)
     [[3.0], [4.0], [5.0]]
-    
+
     """
     pass
 
@@ -2077,7 +2178,7 @@ def test095():
     >>> A = Matrix([[1,2],[3,4]])
     >>> print A + 1j
     [[(1+1j), (2+0j)], [(3+0j), (4+1j)]]
-    
+
     """
     pass
 
@@ -2092,7 +2193,7 @@ def test096():
     >>> b = Matrix([[1],[2],[3]])
     >>> print(b*b)       # scalar product
     14
-    
+
     """
     pass
 
@@ -2120,7 +2221,7 @@ def test097():
     >>> f(B1, points, 'la5.png')
     >>> B2 = Matrix([[0.2,0.4],[0.5,0.3]])
     >>> f(B2, points, 'la6.png')
-    
+
     """
     pass
 
@@ -2134,7 +2235,7 @@ def test098():
     [[1.0, 0.0], [0.0, 1.0]]
     >>> print A/2
     [[0.5, 1.0], [2.0, 4.5]]
-    
+
     """
     pass
 
@@ -2144,7 +2245,7 @@ def test099():
     >>> A = Matrix([[1,2],[3,4]])
     >>> print A.T
     [[1, 3], [2, 4]]
-    
+
     """
     pass
 
@@ -2156,7 +2257,7 @@ def test100():
     >>> x = (1/A)*b
     >>> print x
     [[-1.0], [3.0], [-1.0]]
-    
+
     """
     pass
 
@@ -2169,7 +2270,7 @@ def test101():
     >>> A = Matrix([[1,2],[3,4]])
     >>> print condition_number(A)
     21.0
-    
+
     """
     pass
 
@@ -2179,7 +2280,7 @@ def test102():
     >>> A = Matrix([[1,2],[3,4]])
     >>> print exp(A)
     [[51.96..., 74.73...], [112.10..., 164.07...]]
-    
+
     """
     pass
 
@@ -2190,7 +2291,7 @@ def test103():
     >>> L = Cholesky(A)
     >>> print is_almost_zero(A - L*L.T)
     True
-    
+
     """
     pass
 
@@ -2215,7 +2316,7 @@ def test104():
     ...      ).errorbar(points[:10],legend='o(t)'
     ...      ).plot([(p[0],fitting_f(p[0])) for p in points[:10]],legend='e(t)'
     ...      ).save('images/polynomialfit.png')
-    
+
     """
     pass
 
@@ -2231,7 +2332,7 @@ def test105():
     1030...
     >>> print 1000.0*data[-1]['adjusted_close']/data[0]['adjusted_close']
     1228...
-    
+
     """
     pass
 
@@ -2246,7 +2347,7 @@ def test106():
     >>> U,e = Jacobi_eigenvalues(A)
     >>> print is_almost_zero(U*Matrix.diagonal(e)*U.T-A)
     True
-    
+
     """
     pass
 
@@ -2261,7 +2362,7 @@ def test107():
     >>> Canvas(title='SP100 eigenvalues',xlab='i',ylab='e[i]'
     ...       ).plot([(i,ei) for i,ei, in enumerate(e)]
     ...       ).save('images/sp100eigen.png')
-    
+
     """
     pass
 
@@ -2288,7 +2389,7 @@ def test108():
     >>> y = y.reshape(m,m)
     >>> Canvas(title="Defocused image").imshow(y.tolist()).save('images/defocused.png')
     >>> Canvas(title="refocus image").imshow(z.tolist()).save('images/refocused.png')
-    
+
     """
     pass
 
@@ -2298,7 +2399,7 @@ def test109():
     >>> def f(x): return (x-2)*(x-5)/10
     >>> print round(solve_fixed_point(f,1.0,rp=0),4)
     2.0
-    
+
     """
     pass
 
@@ -2308,7 +2409,7 @@ def test110():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(solve_bisection(f,1.0,3.0),4)
     2.0
-    
+
     """
     pass
 
@@ -2318,7 +2419,7 @@ def test111():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(solve_newton(f,1.0),4)
     2.0
-    
+
     """
     pass
 
@@ -2328,7 +2429,7 @@ def test112():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(solve_secant(f,1.0),4)
     2.0
-    
+
     """
     pass
 
@@ -2338,7 +2439,7 @@ def test113():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(optimize_bisection(f,2.0,5.0),4)
     3.5
-    
+
     """
     pass
 
@@ -2348,7 +2449,7 @@ def test114():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(optimize_newton(f,3.0),3)
     3.5
-    
+
     """
     pass
 
@@ -2358,7 +2459,7 @@ def test115():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(optimize_secant(f,3.0),3)
     3.5
-    
+
     """
     pass
 
@@ -2368,7 +2469,7 @@ def test116():
     >>> def f(x): return (x-2)*(x-5)
     >>> print round(optimize_golden_search(f,2.0,5.0),3)
     3.5
-    
+
     """
     pass
 
@@ -2382,7 +2483,7 @@ def test117():
     >>> x = (1,1,1)
     >>> print round(df0(x),4), round(df1(x),4), round(df2(x),4)
     2.0 8.0 5.0
-    
+
     """
     pass
 
@@ -2394,7 +2495,7 @@ def test118():
     [[1.999999...], [7.999999...], [4.999999...]]
     >>> print hessian(f, x=(1,1,1))
     [[0.0, 0.0, 0.0], [0.0, 0.0, 5.000000...], [0.0, 5.000000..., 0.0]]
-    
+
     """
     pass
 
@@ -2404,7 +2505,7 @@ def test119():
     >>> def f(x): return (2.0*x[0]+3.0*x[1]+5.0*x[1]*x[2], 2.0*x[0])
     >>> print jacobian(f, x=(1,1,1))
     [[1.9999999..., 7.999999..., 4.9999999...], [1.9999999..., 0.0, 0.0]]
-    
+
     """
     pass
 
@@ -2414,7 +2515,7 @@ def test120():
     >>> def f(x): return [x[0]+x[1], x[0]+x[1]**2-2]
     >>> print solve_newton_multi(f, x=(0,0))
     [1.0..., -1.0...]
-    
+
     """
     pass
 
@@ -2424,7 +2525,7 @@ def test121():
     >>> def f(x): return (x[0]-2)**2+(x[1]-3)**2
     >>> print optimize_newton_multi(f, x=(0,0))
     [2.0, 3.0]
-    
+
     """
     pass
 
@@ -2436,7 +2537,7 @@ def test122():
     >>> ab, chi2 = fit(data,fs,[5])
     >>> print ab, chi2
     [0.999..., 2.000..., 300.000..., 10.000...] ...
-    
+
     """
     pass
 
@@ -2454,7 +2555,7 @@ def test123():
     1.9899...
     >>> print 1.0-cos(3)
     1.9899...
-    
+
     """
     pass
 
@@ -2468,7 +2569,7 @@ def test124():
     1.99373945223
     >>> print integrate_quadrature_naive(sin,0,3,n=2,order=4)
     1.99164529955
-    
+
     """
     pass
 
@@ -2480,7 +2581,7 @@ def test128():
     1.000...
     >>> print sd(S)
     0.4...
-    
+
     """
     pass
 
@@ -2491,7 +2592,7 @@ def test129():
     >>> def payoff(x): return 20.0 if x==6 else -5.0
     >>> print E(payoff,S)
     -0.83333...
-    
+
     """
     pass
 
@@ -2503,7 +2604,7 @@ def test130():
     ...     return sum(f(random.uniform(a,b)) for i in xrange(N))/N*(b-a)
     >>> print integrate_mc(sin,0,pi,N=10000)
     2.000....
-    
+
     """
     pass
 
@@ -2528,7 +2629,7 @@ def test131():
     0.0802804358268
     >>> print correlation(X,Y)
     0.479067813484
-    
+
     """
     pass
 
@@ -2545,11 +2646,12 @@ def test132():
     ...       ).hist(make_set(4),legend='N=4').save('images/central4.png')
     >>> Canvas(title='Central Limit Theorem',xlab='y',ylab='p(y)'
     ...       ).hist(make_set(8),legend='N=8').save('images/central8.png')
-    
+
     """
     pass
 
-if __name__=='__main__':
-    import os,doctest
+if __name__ == '__main__':
+    import os
+    import doctest
     if not os.path.exists('images'): os.mkdir('images')
     doctest.testmod(optionflags=doctest.ELLIPSIS)