{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# \n", "\n", "\n", "\n", "Copyright 2016 Allen B. Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "import thinkstats2\n", "import thinkplot\n", "\n", "import pandas as pd\n", "import numpy as np\n", "\n", "from fractions import Fraction\n", "\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def scalar_product(x, y):\n", " x = np.asarray(x)\n", " y = np.asarray(y)\n", " return np.sum(x * y)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "32" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scalar_product([1,2,3], (4,5,6))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "12" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scalar_product([1,2,3], 2)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "12" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scalar_product([1,2,3], [2])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "operands could not be broadcast together with shapes (3,) (4,) \n" ] } ], "source": [ "try:\n", " scalar_product([1,2,3], (4,5,6,7))\n", "except ValueError as e:\n", " print(e)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class ArrayWrapper:\n", " def __init__(self, array):\n", " self.array = np.asarray(array)\n", " \n", " def __eq__(self, other):\n", " return np.array_equal(self.array, other.array)\n", " \n", " def __add__(self, other):\n", " return self.__class__(self.array + other.array)\n", " \n", " def __sub__(self, other):\n", " return self.__class__(self.array - other.array)\n", " \n", " def __str__(self):\n", " return str(self.array)\n", " \n", " def __repr__(self):\n", " return '%s(\\n%s)' % (self.__class__.__name__, str(self.array))\n", " \n", " def __len__(self):\n", " return len(self.array)\n", " \n", " def __getitem__(self, index):\n", " return self.array[index]\n", " \n", " def __setitem__(self, index, elt):\n", " self.array[index] = elt\n", " \n", " @property\n", " def t(self):\n", " return self.__class__(self.array.transpose())\n", " \n", "class Vector(ArrayWrapper):\n", " def __mul__(self, other):\n", " return scalar_product(self.array, other.array)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def random_array(*shape):\n", " return np.random.randint(1, 10, shape)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[7 2 1])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = Vector(random_array(3))\n", "x" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(7, 2, 1)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x[0], x[1], x[2]" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x[1] += 1" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7\n", "3\n", "1\n" ] } ], "source": [ "for elt in x:\n", " print(elt)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[7 3 1])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = Vector(x.array)\n", "y" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x == y" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[7 3 1])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.t" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x == x.t" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[1 8 6])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = Vector(random_array(3))\n", "y" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x == y" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[ 8 11 7])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x+y" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[ 6 -5 -5])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x-y" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "37" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x*y" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def mm_product(array1, array2):\n", " dtype = np.result_type(array1, array2)\n", " array = np.zeros((len(array1), len(array2)), dtype=dtype)\n", " for i, row1 in enumerate(array1):\n", " for j, row2 in enumerate(array2):\n", " array[i][j] = scalar_product(row1, row2)\n", " return array" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class Matrix(ArrayWrapper):\n", " \n", " def __mul__(self, other):\n", " return self.__class__(mm_product(self.array, other.t.array))\n", " \n", " def __truediv__(self, other):\n", " return self.__class__(np.linalg.solve(self.array, other.array.flat))" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[3 2 3]\n", " [6 8 3]\n", " [1 4 5]])" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = Matrix(random_array(3, 3))\n", "A" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(A)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3 2 3]\n", "[6 8 3]\n", "[1 4 5]\n" ] } ], "source": [ "for row in A:\n", " print(row)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[4 6 3]\n", " [4 5 5]\n", " [3 8 5]])" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = Matrix(random_array(3, 3))\n", "B" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[ 7 8 6]\n", " [10 13 8]\n", " [ 4 12 10]])" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A+B" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[-1 -4 0]\n", " [ 2 3 -2]\n", " [-2 -4 0]])" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A-B" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[ 29 52 34]\n", " [ 65 100 73]\n", " [ 35 66 48]])" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A*B" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 29, 52, 34],\n", " [ 65, 100, 73],\n", " [ 35, 66, 48]])" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A.array.dot(B.array)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[8 6 4])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = Vector(random_array(3))\n", "x" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[ 64 48 32]\n", " [136 102 68]\n", " [ 80 60 40]])" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A*x" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def mv_product(A, x):\n", " dtype = np.result_type(A, x)\n", " array = np.zeros(len(A), dtype=dtype)\n", " for i, row in enumerate(A):\n", " array[i] = scalar_product(row, x)\n", " return Vector(array)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Vector(\n", "[ 48 108 52])" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mv_product(A.array, x.array)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 48, 108, 52])" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A.array.dot(x.array)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[8]\n", " [9]\n", " [3]])" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = Matrix(random_array(3, 1))\n", "x" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x == x.t" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[154]])" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.t * x" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[64 72 24]\n", " [72 81 27]\n", " [24 27 9]])" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x * x.t" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[160]\n", " [180]\n", " [ 60]])" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x * x" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[ 51]\n", " [129]\n", " [ 59]])" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A * x" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 51],\n", " [129],\n", " [ 59]])" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A.array.dot(x.array)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[2]])" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scalar = Matrix([[2]])\n", "scalar" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scalar == scalar.t" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[4]])" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scalar * scalar" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[16]\n", " [18]\n", " [ 6]])" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x * scalar" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[16]\n", " [34]\n", " [20]])" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A * scalar" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Matrix(\n", "[[ 51]\n", " [129]\n", " [ 59]])" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = A * x\n", "b" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 51],\n", " [129],\n", " [ 59]])" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b.array" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 8.],\n", " [ 9.],\n", " [ 3.]])" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.solve(A.array, b.array)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 8. 9. 3.]\n" ] } ], "source": [ "print(A / b)" ] }, { "cell_type": "code", "execution_count": 64, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(3, 3)" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A.array.shape" ] }, { "cell_type": "code", "execution_count": 65, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(3, 1)" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b.array.shape" ] }, { "cell_type": "code", "execution_count": 66, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3 2 3 51]\n", " [6 8 3 129]\n", " [1 4 5 59]]\n" ] } ], "source": [ "m = np.hstack([A.array, b.array]).astype(Fraction)\n", "print(m)" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3 2 3 51]\n", " [3 6 0 78]\n", " [1 4 5 59]]\n" ] } ], "source": [ "m[1] -= m[0]\n", "print(m)" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[3, 2, 3],\n", " [3, 6, 0],\n", " [1, 4, 5]], dtype=object)" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m[:, :-1]" ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([51, 78, 59], dtype=object)" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m[:, -1]" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def solve_augmented(m):\n", " m = m.astype(float)\n", " return np.linalg.solve(m[:, :-1], m[:,-1])" ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 8. 9. 3.]\n" ] } ], "source": [ "print(solve_augmented(m))" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(array([3, 2, 3, 51], dtype=object),\n", " 3,\n", " 3,\n", " array([Fraction(3, 1), Fraction(2, 1), Fraction(3, 1), Fraction(51, 1)], dtype=object))" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "row1 = 0\n", "row2 = 1\n", "col = 0\n", "pivot = m[row1, col]\n", "victim = m[row2, col]\n", "m[row1], pivot, victim, m[row1] * Fraction(victim, pivot)" ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3 2 3 51]\n", " [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)]\n", " [1 4 5 59]]\n" ] } ], "source": [ "m[row2] -= m[row1] * Fraction(victim, pivot)\n", "print(m)" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def clobber(m, row1, row2, col):\n", " pivot = m[row1, col]\n", " victim = m[row2, col]\n", " m[row2] -= m[row1] * Fraction(victim, pivot)" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3 2 3 51]\n", " [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)]\n", " [Fraction(0, 1) Fraction(10, 3) Fraction(4, 1) Fraction(42, 1)]]\n" ] } ], "source": [ "clobber(m, 0, 2, 0)\n", "print(m)" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3 2 3 51]\n", " [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)]\n", " [Fraction(0, 1) Fraction(0, 1) Fraction(13, 2) Fraction(39, 2)]]\n" ] } ], "source": [ "clobber(m, 1, 2, 1)\n", "print(m)" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3 2 3 51]\n", " [Fraction(0, 1) Fraction(4, 1) Fraction(-3, 1) Fraction(27, 1)]\n", " [Fraction(0, 1) Fraction(0, 1) Fraction(1, 1) Fraction(3, 1)]]\n" ] } ], "source": [ "m[2] /= m[2,2]\n", "print(m)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }