In [1]:

    

from __future__ import print_function, division

import thinkstats2
import thinkplot

import pandas as pd
import numpy as np

from fractions import Fraction

%matplotlib inline


    

In [2]:

    

def scalar_product(x, y):
    x = np.asarray(x)
    y = np.asarray(y)
    return np.sum(x * y)


    

In [3]:

    

scalar_product([1,2,3], (4,5,6))


    

    
Out[3]:





32

In [4]:

    

scalar_product([1,2,3], 2)


    

    
Out[4]:





12

In [5]:

    

scalar_product([1,2,3], [2])


    

    
Out[5]:





12

In [6]:

    

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,) 

In [7]:

    

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)


    

In [8]:

    

def random_array(*shape):
    return np.random.randint(1, 10, shape)


    

In [9]:

    

x = Vector(random_array(3))
x


    

    
Out[9]:





Vector(
[7 2 1])

In [10]:

    

x[0], x[1], x[2]


    

    
Out[10]:





(7, 2, 1)

In [11]:

    

x[1] += 1


    

In [12]:

    

for elt in x:
    print(elt)


    

    




7
3
1

In [13]:

    

y = Vector(x.array)
y


    

    
Out[13]:





Vector(
[7 3 1])

In [14]:

    

x == y


    

    
Out[14]:





True

In [15]:

    

x.t


    

    
Out[15]:





Vector(
[7 3 1])

In [16]:

    

x == x.t


    

    
Out[16]:





True

In [17]:

    

y = Vector(random_array(3))
y


    

    
Out[17]:





Vector(
[1 8 6])

In [18]:

    

x == y


    

    
Out[18]:





False

In [19]:

    

x+y


    

    
Out[19]:





Vector(
[ 8 11  7])

In [20]:

    

x-y


    

    
Out[20]:





Vector(
[ 6 -5 -5])

In [21]:

    

x*y


    

    
Out[21]:





37

In [31]:

    

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


    

In [32]:

    

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))


    

In [33]:

    

A = Matrix(random_array(3, 3))
A


    

    
Out[33]:





Matrix(
[[3 2 3]
 [6 8 3]
 [1 4 5]])

In [34]:

    

len(A)


    

    
Out[34]:





3

In [35]:

    

for row in A:
    print(row)


    

    




[3 2 3]
[6 8 3]
[1 4 5]

In [36]:

    

B = Matrix(random_array(3, 3))
B


    

    
Out[36]:





Matrix(
[[4 6 3]
 [4 5 5]
 [3 8 5]])

In [37]:

    

A+B


    

    
Out[37]:





Matrix(
[[ 7  8  6]
 [10 13  8]
 [ 4 12 10]])

In [38]:

    

A-B


    

    
Out[38]:





Matrix(
[[-1 -4  0]
 [ 2  3 -2]
 [-2 -4  0]])

In [39]:

    

A*B


    

    
Out[39]:





Matrix(
[[ 29  52  34]
 [ 65 100  73]
 [ 35  66  48]])

In [40]:

    

A.array.dot(B.array)


    

    
Out[40]:





array([[ 29,  52,  34],
       [ 65, 100,  73],
       [ 35,  66,  48]])

In [41]:

    

x = Vector(random_array(3))
x


    

    
Out[41]:





Vector(
[8 6 4])

In [42]:

    

A*x


    

    
Out[42]:





Matrix(
[[ 64  48  32]
 [136 102  68]
 [ 80  60  40]])

In [45]:

    

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)


    

In [46]:

    

mv_product(A.array, x.array)


    

    
Out[46]:





Vector(
[ 48 108  52])

In [47]:

    

A.array.dot(x.array)


    

    
Out[47]:





array([ 48, 108,  52])

In [48]:

    

x = Matrix(random_array(3, 1))
x


    

    
Out[48]:





Matrix(
[[8]
 [9]
 [3]])

In [49]:

    

x == x.t


    

    
Out[49]:





False

In [50]:

    

x.t * x


    

    
Out[50]:





Matrix(
[[154]])

In [51]:

    

x * x.t


    

    
Out[51]:





Matrix(
[[64 72 24]
 [72 81 27]
 [24 27  9]])

In [52]:

    

x * x


    

    
Out[52]:





Matrix(
[[160]
 [180]
 [ 60]])

In [53]:

    

A * x


    

    
Out[53]:





Matrix(
[[ 51]
 [129]
 [ 59]])

In [54]:

    

A.array.dot(x.array)


    

    
Out[54]:





array([[ 51],
       [129],
       [ 59]])

In [55]:

    

scalar = Matrix([[2]])
scalar


    

    
Out[55]:





Matrix(
[[2]])

In [56]:

    

scalar == scalar.t


    

    
Out[56]:





True

In [57]:

    

scalar * scalar


    

    
Out[57]:





Matrix(
[[4]])

In [58]:

    

x * scalar


    

    
Out[58]:





Matrix(
[[16]
 [18]
 [ 6]])

In [59]:

    

A * scalar


    

    
Out[59]:





Matrix(
[[16]
 [34]
 [20]])

In [60]:

    

b = A * x
b


    

    
Out[60]:





Matrix(
[[ 51]
 [129]
 [ 59]])

In [61]:

    

b.array


    

    
Out[61]:





array([[ 51],
       [129],
       [ 59]])

In [62]:

    

np.linalg.solve(A.array, b.array)


    

    
Out[62]:





array([[ 8.],
       [ 9.],
       [ 3.]])

In [63]:

    

print(A / b)


    

    




[ 8.  9.  3.]

In [64]:

    

A.array.shape


    

    
Out[64]:





(3, 3)

In [65]:

    

b.array.shape


    

    
Out[65]:





(3, 1)

In [66]:

    

m = np.hstack([A.array, b.array]).astype(Fraction)
print(m)


    

    




[[3 2 3 51]
 [6 8 3 129]
 [1 4 5 59]]

In [67]:

    

m[1] -= m[0]
print(m)


    

    




[[3 2 3 51]
 [3 6 0 78]
 [1 4 5 59]]

In [68]:

    

m[:, :-1]


    

    
Out[68]:





array([[3, 2, 3],
       [3, 6, 0],
       [1, 4, 5]], dtype=object)

In [69]:

    

m[:, -1]


    

    
Out[69]:





array([51, 78, 59], dtype=object)

In [70]:

    

def solve_augmented(m):
    m = m.astype(float)
    return np.linalg.solve(m[:, :-1], m[:,-1])


    

In [71]:

    

print(solve_augmented(m))


    

    




[ 8.  9.  3.]

In [72]:

    

row1 = 0
row2 = 1
col = 0
pivot = m[row1, col]
victim = m[row2, col]
m[row1], pivot, victim, m[row1] * Fraction(victim, pivot)


    

    
Out[72]:





(array([3, 2, 3, 51], dtype=object),
 3,
 3,
 array([Fraction(3, 1), Fraction(2, 1), Fraction(3, 1), Fraction(51, 1)], dtype=object))

In [73]:

    

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]]

In [74]:

    

def clobber(m, row1, row2, col):
    pivot = m[row1, col]
    victim = m[row2, col]
    m[row2] -= m[row1] * Fraction(victim, pivot)


    

In [75]:

    

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)]]

In [76]:

    

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)]]

In [77]:

    

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)]]

In [ ]: