Using numpy¶

# download the data locally
if not os.path.exists('populations.txt'):
    ! wget http://scipy-lectures.github.io/_downloads/populations.txt

# peek at the file to see its structure
! head -n 6 populations.txt

# load data into a numpy array
data = np.loadtxt('populations.txt').astype('int')
data[:5, :]

# provide convenient named variables
populations = data[:, 1:]
year, hare, lynx, carrot = data.T

# The mean and std of the populations of each species for the years in the period
print "Mean (hare, lynx, carrot):", populations.mean(axis=0)
print "Std (hare, lynx, carrot):", populations.std(axis=0)

# Which year each species had the largest population.
print "Year with largest population (hare, lynx, carrot)",
print year[np.argmax(populations, axis=0)]

# Which species has the largest population for each year.
species = ['hare', 'lynx', 'carrot']
zip(year, np.take(species, np.argmax(populations, axis=1)))

# Which years any of the populations is above 50000
print year[np.any(populations > 50000, axis=1)]

# The top 2 years for each species when they had the lowest populations.
print year[np.argsort(populations, axis=0)[:2]]

plt.plot(year, lynx, 'r-', year, np.gradient(hare), 'b--')
plt.legend(['lynx', 'grad(hare)'], loc='best')
print np.corrcoef(lynx, np.gradient(hare))

x = np.array([1,2,3,4,5,6])
print x
print 'dytpe', x.dtype
print 'shape', x.shape
print 'strides', x.strides

x.shape = (2,3)
print x
print 'dytpe', x.dtype
print 'shape', x.shape
print 'strides', x.strides

x = x.astype('complex')
print x
print 'dytpe', x.dtype
print 'shape', x.shape
print 'strides', x.strides

# operations across rows, cols or entire matrix
print(xs.max())
print(xs.max(axis=0)) # max of each col
print(xs.max(axis=1)) # max of each row

# A funcitonal rather than object-oriented approacha also wokrs
print(np.max(xs, axis=0))
print(np.max(xs, axis=1))

# broadcasting
xs = np.arange(12).reshape(2,6)
print(xs, '\n')
print(xs * 10, '\n')

# broadcasting just works when doing column-wise operations
col_means = xs.mean(axis=0)
print(col_means, '\n')
print(xs + col_means, '\n')

# but needs a little more work for row-wise operations
row_means = xs.mean(axis=1)[:, np.newaxis]
print(row_means)
print(xs + row_means)

# convert matrix to have zero mean and unit standard deviation using col summary statistics
print((xs - xs.mean(axis=0))/xs.std(axis=0))

# convert matrix to have zero mean and unit standard deviation using row summary statistics
print((xs - xs.mean(axis=1)[:, np.newaxis])/xs.std(axis=1)[:, np.newaxis])

# broadcasting for outer product
# e.g. create the 12x12 multiplication toable
u = np.arange(1, 13)
u[:,None] * u[None,:]

def distance_matrix_py(pts):
    """Returns matrix of pairwise Euclidean distances. Pure Python version."""
    n = len(pts)
    p = len(pts[0])
    m = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            s = 0
            for k in range(p):
                s += (pts[i,k] - pts[j,k])**2
            m[i, j] = s**0.5
    return m

def distance_matrix_np(pts):
    """Returns matrix of pairwise Euclidean distances. Vectorized numpy version."""
    return np.sum((pts[None,:] - pts[:, None])**2, -1)**0.5

pts = np.array([(0,0), (4,0), (4,3), (0,3)])

distance_matrix_py(pts)

distance_matrix_np(pts)

# Broaccasting and vectorization is faster than looping
%timeit distance_matrix_py(pts)
%timeit distance_matrix_np(pts)

xs = np.linspace(0, 2*np.pi, 100)
ys = np.sin(xs) # np.sin is a universal function
plt.plot(xs, ys);

# operators also perform elementwise operations by default

xs = np.arange(10)
print xs
print -xs
print xs+xs
print xs*xs
print xs**3
print xs < 5

from numpy.core.umath_tests import matrix_multiply

print matrix_multiply.signature

us = np.random.random((5, 2, 3)) # 5 2x3 matrics
vs = np.random.random((5, 3, 4)) # 5 3x4 matrices
# perform matrix multiplication for each of the 5 sets of matrices
ws = matrix_multiply(us, vs)
print ws.shape
print ws

import numpy.random as npr
npr.seed(123) # fix seed for reproducible results

# 10 trials of rolling a fair 6-sided 100 times
roll = 1.0/6
x = npr.multinomial(100, [roll]*6, 10)
x

# uniformly distributed numbers in 2D
x = npr.uniform(-1, 1, (100, 2))
plt.scatter(x[:,0], x[:,1], s=50)
plt.axis([-1.05, 1.05, -1.05, 1.05]);

# ranodmly shuffling a vector
x = np.arange(10)
npr.shuffle(x)
x

# radnom permutations
npr.permutation(10)

# radnom selection without replacement
x = np.arange(10,20)
npr.choice(x, 10, replace=False)

# radnom selection with replacement
npr.choice(x, (5, 10), replace=True) # this is default

# toy example - estimating pi inefficiently
n = 1e6
x = npr.uniform(-1,1,(n,2))
4.0*np.sum(x[:,0]**2 + x[:,1]**2 < 1)/n

import scipy.stats as stats

# Create a "frozen" distribution - i.e. a partially applied function
dist = stats.norm(10, 2)

#  same a rnorm
dist.rvs(10)

# same as pnorm
dist.pdf(np.linspace(5, 15, 10))

# same as dnorm
dist.cdf(np.linspace(5, 15, 11))

# same as qnorm
dist.ppf(dist.cdf(np.linspace(5, 15, 11)))

import scipy.linalg as la

A = np.array([[1,2],[3,4]])
b = np.array([1,4])
print(A)
print(b)

# Matrix operations
import numpy as np
import scipy.linalg as la
from functools import reduce

A = np.array([[1,2],[3,4]])
print(np.dot(A, A))
print(A)
print(la.inv(A))
print(A.T)

x = la.solve(A, b) # do not use x = dot(inv(A), b) as it is inefficient and numerically unstable
print(x)
print(np.dot(A, x) - b)

m = np.arange(12).reshape((3,4))

# YOUR CODE HERE

m = np.arange(12).reshape((3,4))
print m
print

print "OVerall", m.mean()
print "Row", m.mean(1)
print "Columne", m.mean(0)

u = np.array([1,3,5,7])
v = np.array([2,4,6,8])

# YOUR CODE HERE

u = np.array([1,3,5,7])
v = np.array([2,4,6,8])

print np.outer(u, v)
print
print np.array([[u_ * v_ for v_ in v] for u_ in u])
print
print u[:,None] * v[None,:]

array([[ 0.49722235,  0.88833973,  0.07289358,  0.12375223,  0.39659254,
         0.70267114],
       [ 0.3954172 ,  0.889077  ,  0.71286225,  0.06353112,  0.68107965,
         0.17186995],
       [ 0.74821206,  0.92692111,  0.24871227,  0.26904958,  0.80410194,
         0.22304055],
       [ 0.22582605,  0.37671244,  0.96510957,  0.88819053,  0.14654176,
         0.33987323]])

array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.74821206,  0.92692111,  0.24871227,  0.26904958,  0.80410194,
         0.22304055],
       [ 0.22582605,  0.37671244,  0.96510957,  0.88819053,  0.14654176,
         0.33987323]])

# YOUR CODE HERE

xs = np.random.random((10,6))
print xs
print
xs[(xs < 0.1).any(axis=1), :] = 0
print xs

# YOUR CODE HERE

xs = np.linspace(0, 2*np.pi, 100)
print xs[::10]
print
print xs[::-1]
print
idx = np.abs(np.sin(xs)-np.cos(xs)) < 0.1
print xs[idx]
print
plt.scatter(xs[idx], np.sin(xs[idx]))
plt.plot(xs, np.sin(xs), xs, np.cos(xs));

# YOUR CODE HERE

ns = np.arange(1, 11)
m = ns[:, None] * ns[None, :]
print m
print
print m.trace()
print
print np.flipud(m).diagonal()
print
print m.diagonal(offset=1)

A = np.array([
    [1,  2, 1],
    [6, -1, 0],
    [-1,-2,-1]
])

# YOUR CODE HERE

A = np.array([
    [1,  2, 1],
    [6, -1, 0],
    [-1,-2,-1]
])

dotm = lambda *args: reduce(np.dot, args)

u, V = la.eig(A)
P = V
D = np.diag(u)
print P
print
print np.real_if_close(np.round(u))
print

np.real_if_close(np.round(dotm(P, D, la.inv(P)), 6))

from sympy import symbols, init_printing, roots, solve, eye
from sympy.matrices import Matrix

init_printing()

x = symbols('x')

M = Matrix([
    [1,  2, 1],
    [6, -1, 0],
    [-1,-2,-1]
])

M

# Find characteristic polynomial
poly = M.charpoly(x)
poly.as_poly()

# eigenvalues are the roots
roots(poly)

#  Provided function
def plot_matrix_transform(m):
    """Show the geometric effect of m on the standard unit vectors e1 and e2."""

    e1 = np.array([1,0])
    e2 = np.array([0,1])
    v1 = np.dot(m, e1)
    v2 = np.dot(m, e2)

    X = np.zeros((2,2))
    Y = np.zeros((2,2))
    pts = np.array([e1,e2,v1,v2])
    U = pts[:, 0]
    V = pts[:, 1]
    C = [0,1,0,1]

    xmin = min(-1, U.min())
    xmax = max(1, U.max())
    ymin = min(-1, V.min())
    ymax = max(-1, V.max())

    plt.figure(figsize=(6,6))
    plt.quiver(X, Y, U, V, C, angles='xy', scale_units='xy', scale=1)
    plt.axis([xmin, xmax, ymin, ymax]);

### Example usage
m = np.array([[1,2],[3,4]])
plot_matrix_transform(m)

# YOUR CODE HERE

A1 = np.diag([2,3])
A2 = np.array([[2,3],[1,1.5]])
A3 = np.array([[np.cos(1), -np.sin(1)], [np.sin(1), np.cos(1)]])
A4 = np.array([[np.cos(1), np.sin(1)], [np.sin(1), -np.cos(1)]])

print A1, la.det(A1)
print
print A2, la.det(A2)
print
print A3, la.det(A3)
print
print A4, la.det(A4)

print A3.dot(A3.T)

print A4.dot(A4.T)

# A diagnoal matrix simply scales the vectors
# This gives insight into what the eigendecomposition tells us
plot_matrix_transform(A1)

# A singluar matrix collapses one vector onto another
# The determinant is zero becasue the parallelogram area is zero
plot_matrix_transform(A2)

# An orthogoanl matrix preservees length and angle
# Hence the area is also preserved and the determinant is 1
# In 2D it is etiher a rotation (shown here)
plot_matrix_transform(A3)

# or a refelction
# The reflection does not preserve orietnation
# This is indicated by the determinatn being -1
plot_matrix_transform(A4)

# The tranpose of an orthogonal matrix is its inverse
plot_matrix_transform(A3.T)

def transform(u, v):
    """Retruns a matrix that converts u into e1 (1,0) and v into e2 (0,1)."""
    return la.inv(np.vstack([u, v]).T)

u = np.random.random(2)
v = np.random.random(2)

M = transform(u, v)
print u, M.dot(u)
print v, M.dot(v)

plt.figure(figsize=(12,8))
x = np.load('x.npy')
y = np.load('y.npy')
plt.plot(x, y, 'o')

### YOUR CODE HERE
p0 = np.poly1d(np.polyfit(x, y, 0))
p2 = np.poly1d(np.polyfit(x, y, 2))
p5 = np.poly1d(np.polyfit(x, y, 5))
p50 = np.poly1d(np.polyfit(x, y, 50))

plt.plot(x, p0(x), 'b:', linewidth=2)
plt.plot(x, p2(x), 'b-', linewidth=2)
plt.plot(x, p5(x), 'g--', linewidth=2)
plt.plot(x, p50(x), 'r-.', linewidth=2)
plt.legend(['Data', 'Constant', 'Qaudratic', 'Order 5', 'Order 50'], loc='best', fontsize=20);

%load_ext version_information

%version_information numpy, scipy