# gauss

``````

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)

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

``````