gauss


Copyright 2016 Allen B. Downey

MIT License: https://opensource.org/licenses/MIT


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