from __future__ import print_function, division
import thinkstats2
import thinkplot
import pandas as pd
import numpy as np
from fractions import Fraction
%matplotlib inline
def scalar_product(x, y):
x = np.asarray(x)
y = np.asarray(y)
return np.sum(x * y)
scalar_product([1,2,3], (4,5,6))
32
scalar_product([1,2,3], 2)
12
scalar_product([1,2,3], [2])
12
try:
scalar_product([1,2,3], (4,5,6,7))
except ValueError as e:
print(e)
operands could not be broadcast together with shapes (3,) (4,)
class ArrayWrapper:
def __init__(self, array):
self.array = np.asarray(array)
def __eq__(self, other):
return np.array_equal(self.array, other.array)
def __add__(self, other):
return self.__class__(self.array + other.array)
def __sub__(self, other):
return self.__class__(self.array - other.array)
def __str__(self):
return str(self.array)
def __repr__(self):
return '%s(\n%s)' % (self.__class__.__name__, str(self.array))
def __len__(self):
return len(self.array)
def __getitem__(self, index):
return self.array[index]
def __setitem__(self, index, elt):
self.array[index] = elt
@property
def t(self):
return self.__class__(self.array.transpose())
class Vector(ArrayWrapper):
def __mul__(self, other):
return scalar_product(self.array, other.array)
def random_array(*shape):
return np.random.randint(1, 10, shape)
x = Vector(random_array(3))
x
Vector( [7 2 1])
x[0], x[1], x[2]
(7, 2, 1)
x[1] += 1
for elt in x:
print(elt)
7 3 1
y = Vector(x.array)
y
Vector( [7 3 1])
x == y
True
x.t
Vector( [7 3 1])
x == x.t
True
y = Vector(random_array(3))
y
Vector( [1 8 6])
x == y
False
x+y
Vector( [ 8 11 7])
x-y
Vector( [ 6 -5 -5])
x*y
37
def mm_product(array1, array2):
dtype = np.result_type(array1, array2)
array = np.zeros((len(array1), len(array2)), dtype=dtype)
for i, row1 in enumerate(array1):
for j, row2 in enumerate(array2):
array[i][j] = scalar_product(row1, row2)
return array
class Matrix(ArrayWrapper):
def __mul__(self, other):
return self.__class__(mm_product(self.array, other.t.array))
def __truediv__(self, other):
return self.__class__(np.linalg.solve(self.array, other.array.flat))
A = Matrix(random_array(3, 3))
A
Matrix( [[3 2 3] [6 8 3] [1 4 5]])
len(A)
3
for row in A:
print(row)
[3 2 3] [6 8 3] [1 4 5]
B = Matrix(random_array(3, 3))
B
Matrix( [[4 6 3] [4 5 5] [3 8 5]])
A+B
Matrix( [[ 7 8 6] [10 13 8] [ 4 12 10]])
A-B
Matrix( [[-1 -4 0] [ 2 3 -2] [-2 -4 0]])
A*B
Matrix( [[ 29 52 34] [ 65 100 73] [ 35 66 48]])
A.array.dot(B.array)
array([[ 29, 52, 34], [ 65, 100, 73], [ 35, 66, 48]])
x = Vector(random_array(3))
x
Vector( [8 6 4])
A*x
Matrix( [[ 64 48 32] [136 102 68] [ 80 60 40]])
def mv_product(A, x):
dtype = np.result_type(A, x)
array = np.zeros(len(A), dtype=dtype)
for i, row in enumerate(A):
array[i] = scalar_product(row, x)
return Vector(array)
mv_product(A.array, x.array)
Vector( [ 48 108 52])
A.array.dot(x.array)
array([ 48, 108, 52])
x = Matrix(random_array(3, 1))
x
Matrix( [[8] [9] [3]])
x == x.t
False
x.t * x
Matrix( [[154]])
x * x.t
Matrix( [[64 72 24] [72 81 27] [24 27 9]])
x * x
Matrix( [[160] [180] [ 60]])
A * x
Matrix( [[ 51] [129] [ 59]])
A.array.dot(x.array)
array([[ 51], [129], [ 59]])
scalar = Matrix([[2]])
scalar
Matrix( [[2]])
scalar == scalar.t
True
scalar * scalar
Matrix( [[4]])
x * scalar
Matrix( [[16] [18] [ 6]])
A * scalar
Matrix( [[16] [34] [20]])
b = A * x
b
Matrix( [[ 51] [129] [ 59]])
b.array
array([[ 51], [129], [ 59]])
np.linalg.solve(A.array, b.array)
array([[ 8.], [ 9.], [ 3.]])
print(A / b)
[ 8. 9. 3.]
A.array.shape
(3, 3)
b.array.shape
(3, 1)
m = np.hstack([A.array, b.array]).astype(Fraction)
print(m)
[[3 2 3 51] [6 8 3 129] [1 4 5 59]]
m[1] -= m[0]
print(m)
[[3 2 3 51] [3 6 0 78] [1 4 5 59]]
m[:, :-1]
array([[3, 2, 3], [3, 6, 0], [1, 4, 5]], dtype=object)
m[:, -1]
array([51, 78, 59], dtype=object)
def solve_augmented(m):
m = m.astype(float)
return np.linalg.solve(m[:, :-1], m[:,-1])
print(solve_augmented(m))
[ 8. 9. 3.]
row1 = 0
row2 = 1
col = 0
pivot = m[row1, col]
victim = m[row2, col]
m[row1], pivot, victim, m[row1] * Fraction(victim, pivot)
(array([3, 2, 3, 51], dtype=object), 3, 3, array([Fraction(3, 1), Fraction(2, 1), Fraction(3, 1), Fraction(51, 1)], dtype=object))
m[row2] -= m[row1] * Fraction(victim, pivot)
print(m)
[[3 2 3 51] [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)] [1 4 5 59]]
def clobber(m, row1, row2, col):
pivot = m[row1, col]
victim = m[row2, col]
m[row2] -= m[row1] * Fraction(victim, pivot)
clobber(m, 0, 2, 0)
print(m)
[[3 2 3 51] [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)] [Fraction(0, 1) Fraction(10, 3) Fraction(4, 1) Fraction(42, 1)]]
clobber(m, 1, 2, 1)
print(m)
[[3 2 3 51] [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)] [Fraction(0, 1) Fraction(0, 1) Fraction(13, 2) Fraction(39, 2)]]
m[2] /= m[2,2]
print(m)
[[3 2 3 51] [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)] [Fraction(0, 1) Fraction(0, 1) Fraction(1, 1) Fraction(3, 1)]]