From 1815b6f67cc0510336571f14712719f5b753fac6 Mon Sep 17 00:00:00 2001 From: Sven Riwoldt Date: Sun, 5 Jan 2025 19:01:05 +0100 Subject: [PATCH] von Udacity --- von_Udacity/NOTES.md | 360 +++++++++++++++++++ von_Udacity/OBJECTS_USED_IN_TESTS.md | 165 +++++++++ von_Udacity/README.md | 8 + von_Udacity/hyperplane.py | 147 ++++++++ von_Udacity/line.py | 169 +++++++++ von_Udacity/linear_system.py | 520 +++++++++++++++++++++++++++ von_Udacity/plane.py | 171 +++++++++ von_Udacity/summary.py | 204 +++++++++++ von_Udacity/vector.py | 247 +++++++++++++ 9 files changed, 1991 insertions(+) create mode 100755 von_Udacity/NOTES.md create mode 100755 von_Udacity/OBJECTS_USED_IN_TESTS.md create mode 100755 von_Udacity/README.md create mode 100755 von_Udacity/hyperplane.py create mode 100755 von_Udacity/line.py create mode 100755 von_Udacity/linear_system.py create mode 100755 von_Udacity/plane.py create mode 100755 von_Udacity/summary.py create mode 100755 von_Udacity/vector.py diff --git a/von_Udacity/NOTES.md b/von_Udacity/NOTES.md new file mode 100755 index 0000000..ccd0bff --- /dev/null +++ b/von_Udacity/NOTES.md @@ -0,0 +1,360 @@ +### Vectors + + +A point represents a location, defined with x, y, z (2, -1, 0). +A vector represents a change in direction and it is also defined with x, y, z. +This time x, y and z represent the change in each of the coordinates. +A vectors doesn't have a specific location. + + +Vectors have magnitude and direction. The magnitude of a vector is calculated +using the pythagorian formula adapted to as many dimensions as the vector has. +e.g. 3 dimensions = `square-root(x**2 + y**2 + z**2)` + + +The Zero vector is the one that has all its coordinates equal to 0. The Zero +vector has no direction and no normalization. + +A unit vector is any vector in which it's magnitude is 1. + +The process of finding a unit vector from a given vector is called normalization. +There are 2 steps involved when normilizing a vector: +- A. Find magnitude. +- B. Multiply vector by 1/magnitude of vector. + +#### Dot product +The inner multiplication or dot product multiplication of vectors is a way of +multiplying 2 vectors and getting a scalar (a couple of references [here](http://betterexplained.com/articles/vector-calculus-understanding-the-dot-product/) and [here](http://mathworld.wolfram.com/InnerProduct.html)). +The dot product represents a scalar a number that applies the directional growth +of a vector into another one. + +It is calculated by multiplying their magnitudes times the multiplications of the +cos of their angle. + +`v1 * v2 = magnitude(v1) * magnitude(v2) * cos(v1-v2)` + +cos is bounded by -1 and +1 so: + +`v1 * v2 <= magnitude(v1) * magnitude(v2)` // Inequality of Cauchy-Schwartz. + +if `v1 * v2 == magnitude(v1) * magnitude(v2)` => Angle is 0deg (same angle). + +if `v1 * v2 == - (magnitude(v1) * magnitude(v2))` => Angle is 180deg +(opposite angle). + +if `v1 * v2 == 0 and v1 != 0 and v2 != 0` => Angle is 90deg + +`v * v => magnitude(v) * magnitude(v) * cos(0)` => `magnitude(v) = sqrt(v * v)` + +It can also be expressed as: + +`v1 * v2 = sum of multiplication of each of its coordinates` + +To find out the angle of 2 vectors: + +`arc-cosin ((v1*v2)/(magnitude(v1) * magnitude(2)))` => radians. + +`arc-cosin((normalize(v1) * normalize(v2)))`. + +#### Parallel and orthogonal vectors +2 vectors are parallel if one is a scalar multiple of another one. +e.g: v and 2v and 1/2v and Zero vector + +2 vectors are orthogonal if v1 * v2 is = 0. There are 2 possibilities. +A vector that is 90 deg with another vector or the Zero vector. + +Zero vector is both parallel and orthogonal to all other vectors. +Zero vector is the only one that is orthogonal to itself. + +Orthogonality is a tool for decomposing objects into combinations of simpler +objects in a structured way. + +#### Proyecting a vector into another vector +[Proyecting a vector into another vector](https://www.udacity.com/course/viewer#!/c-ud953/l-4374471116/m-4583493277): + +`v1 = v1-parallel (to base vector) + v1-orthogonal (to base vector)` + +`v1-parallel` is a cathetus, `v2-orthogonal` is the other cathetus and `v1` is the +hipotenuse. + + +`(*)` `cos angle = magnitude(v1-parallel) / magnitude(v1)` => + +`magnitude(v1-parallel) = cos angle * magnitude(v1)`. + + +Being v2 the base vector: + +`(v1 * v2) / (magnitude(v1) * magnitude(v2)) = cos(v1-v2)`. + +Substituting cos from `(*)`: + +`magnitude(v1-parallel) = magnitude(v1) * (v1 * v2) / (magnitude(v1) * magnitude(v2))` => + +(eliminating `magnitude(v1)` in right equation) => + +`magnitude(v1-parallel) = (v1 * v2) / magnitude(v2))` => + +`magnitude(v1-parallel) = v1 * normalize(v2)` + +v1 parallel has the same direction as v2. + +`magnitude(v1-parallel) * normalize(v2) = v1-parallel` => + +`(v1 * normalize(v2)) * normalize(v2) = v1-parallel` + + +#### Cross product +Cross product: It's only possible in 3-dimension vectors and it gives another +vector. + +Cross product of v and w is a vector that is orthogonal to both and its +`magintute = magnitude(v) * magnitude(w) * sin(angle v-w)`. + +Cross product of 2 parallel vectors is the 0 vector. +If v or w are the 0 vector the cross product will be the 0 vector. +Cross product is anticommutative, which means that the order matters, in this +case changing the order results on vectors having opposite directions. + +`cross product(v, w) = - cross product(v, w)` + +The area of the parallelogram created with v and w is the magnitude of the +cross product of v and w. + +The area of the triangle created by those 2 vectors v and w is half of +the area of the parallelogram + +For more information in cross products go [here](http://betterexplained.com/articles/cross-product/) + + +### Intersections + +Flat objects such as lines and planes are defined by linear equations: +- Can add and subtract vars and constants +- Can multiply a var by a constant + +`x + 2y = 1` // `y/2 - 2z = x` => Linear equations + +`x**2 - 1 = y` // `y/x = 3` => Non linear equations + +Finding intersections is solving systems of equations. We are basicaly converting a real world problems (linear equations) into a geometrical problem (n dimensional objects). + +#### Lines and systems with 2 dimensional objects +The most basic problem is finding intersection between lines. + +A line is defined with its basepoint (x0) and the direction of the vector(v). + +`x(t) = x0 + times_scalar(v, t)` => `x(t) = x0 + t * v` (parametrization of line) + +Parametrization of a line (dimensional object) is the act of conveting a system +of equations into a set of one dimension equations that all take the same +parameters in this case "t" (parameter). + +Note that x0 is a point, but it's treated as a vector in order to computer other +points (x(t)) + +A given line has infinetely many base points (basically each of its points). + +A given line has infinetely many directions (we can multiply the direction +vector by any postive scalar and it will give the same direction) + +`x(t) = [1, -2] + t*[3, 1]` + +for `t = 0` => `x = [1, -2]` => Base point + +for `t = 1` => `x = [1, -2] + [3, 1]` = `[3, -1]` => Another point on the line + +A generic way of representing a line: `Ax + By = k` (A, B not both 0) + +k => constant term + +Base point => `(0, k/B)` and `(k/A, 0)` + +if `k = 0` => + +`[A, B] (dot product) [x, y] = 0` => `[x, y]` is orthogonal to `[A, B]` => + +`[B, -A]` is orthogonal to `[A, B]` (that's how we figure out the direction vector) + +Direction vector => `[B, -A]` + +`Ax + By = k` => Gives us an implicit representation of the line such that for +a set of points (x, y): + +`[A, B] (dot product) [x, y] = k` + +A line is defined by its normal vector [A, B] and its constant term k. +An orthogonal or normal vector to the line is the one represented by [A, B]. + +2 lines are parallel if their normal vectors are parallel. + +If 2 lines are not parallel they'll have an intersection in a single point. + +If 2 lines are parallel there are 2 options: +- They never cross +- They are the same line, which means that they are conincidental and intersect +in infinite points + +So 2 lines could have: +- 1 point of intersection +- 0 points of intersection +- infinite points of intersection + +To solve the system: +- First check if 2 lines are parallel by checking if their normal vecotrs are parallel. +- If normal vectors are parallel, then check if they are the same vector. +- Make a line connecting a point of line 1 and line 2 and if that third +line is orthogonal to either normal vectors for line 1 and 2 => line 1 and 2 +are coincident or the same. +- If we have 2 lines that are not parallel, to know how they intersect: + + `Ax + By = k1` + + `Cx + Dy = k2` + +Either A or C need to not be 0, b/c if not the lines would both be horizontal +and these lines intersect with each other. + +if `A = 0` => swap A by C, B by D and k1 by k2 => A will never be 0 + +`(AD - BC)` is never 0 b/c that would mean that both lines are parallel + +`y = (-Ck1 + Ak2) / (AD - BC)` + +`x = (Dk1 - Bk2) / (AD - BC)` + + +If there are no intersections in a system => the system is "Inconsistent" + + +#### Planes and systems with 3 dimensional objects + +Three dimensional spaces are defined by linear equations of 3 variables +`Ax + By + Cz = k` => `[A B C] (dot product) [x y z] = k` + +`[A B C]` => Normal vector (similiar with lines) + +Changing k shifts the plane but it doesn't change the direction of the plane + +Same rules as with lines apply here. + +- Planes that are not parallel intersect on a line not a point. + +With 2 planes (a system of 2, 3-dimenstion equations) could have: +- No intersections (parallel planes) +- Intesect in a line +- Intersect in a whole plane + +With 3 planes (a system of 3, 3-dimenstion equations), the planes could intersect: +- In a line +- No solution +- A plane +- A single point + +The coefficients of a linear equation are the coordinates of a normal vector +to the object the equation defines + +In order to solve a system wtih 3 linear equations and 3 variables we need +a set of "tools" to manipulate its data (they need to be reversible): +- Swap equations +- Multiply each of the sides of the planes that is not 0 +- Add a multiple of one equation to another equation + +Gaussanian elimination: + +- Clear a variable in successive equations until there is only one equation with +one variable. +- At that point 3 eq with 3 variables have been converted to a system of 3 +equations, one with 3 variables, one with 2 variables and another one with +one variable (triangular form) + +The triangular form doesn't represent the same set of planes as the original +system, but it represents the same result. + +If when we are resolving the system and at any point we get the one of the sides +of an equation is 0 and the side is not 0 => Inconsistent system, the system has no +solution. + +If the system gets to a point in which we have 0 = 0, that means that the +last equation that we were soloving was redundant and gives no additional +information. + +The system in this case has been converted from a 3 plane system to a 2 plane +system, and now we are not looking for a single point of intersection but +for a line of intersection + +When the solution is not a point we need to parametrize the solution set. +We need to identify the pivot variables (lead variables), which are variables +in the leading term, when the system is in a triangular form. +The other variables are called the free variable. + +This free variables will become a parameter in the parametrization. So in +order to find the values of the pivotal variables for any given point in the +line we'll use the free variable (parameter) to compute it. + +Then we have to make the coefficient of pivot variable to be 1 + +Then we clear the pivot variable of the second equation in the first +equation leaving the system: + +`x + Az = B` + +`y + Cz = D` + +The leading variables have coefficient one and they only appear on one +equation each. The system is simplified as much as possible giving us the +"reduced row-echelon form" or rref. + +With the rref we are able to get a basepoint and a direction vector that define +the parametrization of the soluton for the system. + +`[x y z] = [(B - AZ) (D - Cz) z ]`=> `[B D 0] + z* [-A -C 1]` + +Base point `[B D 0]` + +Direction vector `[-A -C 1]` + +In order for a system to intersect in only one point it needs at least +3 equations, 3 planes. + +A system is inconsistent if we find `0 = k` during Gaussian elimination + +A consistent solution has a unique solution and each variable is a pivotal +or lead variable. + +Number of free variables is equal to the dimension of the solution set. + +----------- + +Hyperplane in n dimensions is defined as a single linear equation with n +variables. + +A hyperplane in n dimensions is an (n-1) dimensional object. // I don't +understand this concept yet. + +Gaussian elimination and parametrization through the rref works the same way +with planes as with hyperplanes. + +A 2 dimenstional plane in 3 dimensions. Single linear equation with 3 +dimensions + + +--------------- + +#### Bugs in the course + +[In video 23 of Intersections](https://www.udacity.com/course/viewer#!/c-ud953/l-4624329808/e-4897268655/m-4897268656) +it says that parametrization code is in the instructor notes, but it isn't. +It redirects [here](https://storage.googleapis.com/supplemental_media/udacityu/4897268656/linsys.py) +but this page has the linear system code, not the parametrization code. +My implementation of `Parametrization` which was partially copied from one of the +videos is at the end of [linear_system.py](https://github.com/omarrayward/Linear-Algebra-Refresher-Udacity/blob/master/linear_system.py) + +The same quiz video as has some typos in the second ecuation of the first system. +In the video it shows `-0.138x -0.138y + 0.244z = 0.319` but it should be +`-0.131x -0.131y + 0.244z = 0.319` + +The solution [video](https://www.udacity.com/course/viewer#!/c-ud953/l-4624329808/e-4897268655/m-4897268657) +uses the correct equation. + + diff --git a/von_Udacity/OBJECTS_USED_IN_TESTS.md b/von_Udacity/OBJECTS_USED_IN_TESTS.md new file mode 100755 index 0000000..089ad3e --- /dev/null +++ b/von_Udacity/OBJECTS_USED_IN_TESTS.md @@ -0,0 +1,165 @@ +###Vectors + +- Vectors used for test is video 4: +```python +v = Vector([8.218, -9.341]) +w = Vector([-1.129, 2.111]) +``` + +```python +v = Vector([7.119, 8.215]) +w = Vector([-8.223, 0.878]) +``` + +```python +v = Vector([1.671, -1.012, -0.318]) +``` + +- Vectors used for test is video 6: + +```python +v = Vector([-0.221, 7.437]) +``` +```python +v = Vector([8.813, -1.331, -6.247]) +``` +```python +v = Vector([5.581, -2.136]) +``` +```python +v = Vector([1.996, 3.108, -4.554]) +``` + +- Vectors used for test is video 8: +```python +v = Vector([7.887, 4.138]) +w = Vector([-8.802, 6.776]) +``` +```python +v = Vector([-5.955, -4.904, -1.874]) +w = Vector([-4.496, -8.755, 7.103]) +``` +```python +v = Vector([3.183, -7.627]) +w = Vector([-2.668, 5.319]) +``` +```python +v = Vector([7.35, 0.221, 5.188]) +w = Vector([2.751, 8.259, 3.985]) +``` + +- Vectors used for test is video 10: + +```python +v = Vector([-7.579, -7.88]) +w = Vector([22.737, 23.64]) +``` +```python +v = Vector([-2.029, 9.97, 4.172]) +w = Vector([-9.231, -6.639, -7.245]) +``` +```python +v = Vector([-2.328, -7.284, -1.214]) +w = Vector([-1.821, 1.072, -2.94]) +``` +```python +v = Vector([2.118, 4.827]) +w = Vector([0, 0]) +``` + + +- Vectors used for test is video 12: +```python +v = Vector([3.039, 1.879]) +w = Vector([0.825, 2.036]) +``` +```python +v = Vector([-9.88, -3.264, -8.159]) +w = Vector([-2.155, -9.353, -9.473]) +``` +```python +v = Vector([3.009, -6.172, 3.692, -2.51]) +w = Vector([6.404, -9.144, 2.759, 8.718]) +``` + + +- Vectors used for test is video 14: +```python +v1 = Vector([8.462, 7.893, -8.187]) +w1 = Vector([6.984, -5.975, 4.778]) +``` +```python +v2 = Vector([-8.987, -9.838, 5.031]) +w2 = Vector([-4.268, -1.861, -8.866]) +``` +```python +v3 = Vector([1.5, 9.547, 3.691]) +w3 = Vector([-6.007, 0.124, 5.772]) +``` + + +###Intersections + +- Lines used for test is video 4: +```python +line1 = Line(Vector([4.046, 2.836]), 1.21) +line2 = Line(Vector([10.115, 7.09]), 3.025) +``` +```python +line3 = Line(Vector([7.204, 3.182]), 8.68) +line4 = Line(Vector([8.172, 4.114]), 9.883) +``` +```python +line5 = Line(Vector([1.182, 5.562]), 6.744) +line6 = Line(Vector([1.773, 8.343]), 9.525) +``` + +- Planes used for test is video 7: +```python +plane1 = Plane(Vector([-0.412, 3.806, 0.728]), -3.46) +plane2 = Plane(Vector([1.03, -9.515, -1.82]), 8.65) +``` +```python +plane3 = Plane(Vector([2.611, 5.528, 0.283]), 4.6) +plane4 = Plane(Vector([7.715, 8.306, 5.342]), 3.76) +``` +```python +plane5 = Plane(Vector([-7.926, 8.625, -7.212]), -7.952) +plane6 = Plane(Vector([-2.642, 2.875, -2.404]), -2.443) +``` + +- Vectors used for test is video 22: +```python +p1 = Plane(Vector([5.862, 1.178, -10.366]), -8.15) +p2 = Plane(Vector([-2.931, -0.589, 5.183]), -4.075) +``` +```python +p1 = Plane(Vector([8.631, 5.112, -1.816]), -5.113) +p2 = Plane(Vector([4.315, 11.132, -5.27]), -6.775) +p3 = Plane(Vector([-2.158, 3.01, -1.727]), -0.831) +``` +```python +p1 = Plane(Vector([5.262, 2.739, -9.878]), -3.441) +p2 = Plane(Vector([5.111, 6.358, 7.638]), -2.152) +p3 = Plane(Vector([2.016, -9.924, -1.367]), -9.278) +p4 = Plane(Vector([2.167, -13.543, -18.883]), -10.567) +``` + +- Vectors used for test is video 23: +```python +# Note that in the video the second vector is written as: +# Vector([-0.138, -0.138, 0.244]), but bellow is the right implementation +p1 = Plane(Vector([0.786, 0.786, 0.588]), -0.714) +p2 = Plane(Vector([-0.131, -0.131, 0.244]), 0.319) +``` +```python +p1 = Plane(Vector([8.631, 5.112, -1.816]), -5.113) +p2 = Plane(Vector([4.315, 11.132, -5.27]), -6.775) +p3 = Plane(Vector([-2.158, 3.01, -1.727]), -0.831) +``` +```python +p1 = Plane(Vector([0.935, 1.76, -9.365]), -9.955) +p2 = Plane(Vector([0.187, 0.352, -1.873]), -1.991) +p3 = Plane(Vector([0.374, 0.704, -3.746]), -3.982) +p4 = Plane(Vector([-0.561, -1.056, 5.619]), 5.973) +``` diff --git a/von_Udacity/README.md b/von_Udacity/README.md new file mode 100755 index 0000000..31d927a --- /dev/null +++ b/von_Udacity/README.md @@ -0,0 +1,8 @@ +[Udacity Linear algebra refresher course](https://www.udacity.com/course/viewer#!/c-ud953) as a "pre-requisite" to do +[Udacity's Machine Learning nanodegree](https://www.udacity.com/course/machine-learning-engineer-nanodegree--nd009). + +- To avoid syntax errors when completing the assignments copy-paste the vectors, lines and planes in [OBJECTS_USED_IN_TESTS.md](https://github.com/omarrayward/Linear-Algebra-Refresher-Udacity/blob/master/OBJECTS_USED_IN_TESTS.md) used in them. + +- To follow the course explanation with my own notes go to [NOTES.md](https://github.com/omarrayward/Linear-Algebra-Refresher-Udacity/blob/master/NOTES.md) + +- I also found a couple of [bugs](https://github.com/omarrayward/Linear-Algebra-Refresher-Udacity/blob/master/NOTES.md#bugs-in-the-course) in the course. diff --git a/von_Udacity/hyperplane.py b/von_Udacity/hyperplane.py new file mode 100755 index 0000000..f7e479f --- /dev/null +++ b/von_Udacity/hyperplane.py @@ -0,0 +1,147 @@ +from decimal import Decimal, getcontext +from vector import Vector + +getcontext().prec = 30 + + +class MyDecimal(Decimal): + def is_near_zero(self, eps=1e-10): + return abs(self) < eps + + +class Hyperplane(object): + + NO_NONZERO_ELTS_FOUND_MSG = 'No nonzero elements found' + EITHER_DIM_OR_NORMAL_VEC_MUST_BE_PROVIDED_MSG = ( + 'Either the dimension of the hyperplane or the normal vector ' + 'must be provided') + + def __init__(self, dimension=None, normal_vector=None, constant_term=None): + if not dimension and not normal_vector: + raise Exception(self.EITHER_DIM_OR_NORMAL_VEC_MUST_BE_PROVIDED_MSG) + + elif not normal_vector: + self.dimension = dimension + all_zeros = ['0'] * self.dimension + normal_vector = Vector(all_zeros) + else: + self.dimension = normal_vector.dimension + + self.normal_vector = normal_vector + + if not constant_term: + constant_term = Decimal('0') + self.constant_term = Decimal(constant_term) + + self.set_basepoint() + + def set_basepoint(self): + try: + n = self.normal_vector + c = self.constant_term + basepoint_coords = ['0'] * self.dimension + + initial_index = Hyperplane.first_nonzero_index(n) + initial_coefficient = n[initial_index] + + basepoint_coords[initial_index] = c / initial_coefficient + self.basepoint = Vector(basepoint_coords) + + except Exception as e: + if str(e) == Hyperplane.NO_NONZERO_ELTS_FOUND_MSG: + self.basepoint = None + else: + raise e + + def __str__(self): + + num_decimal_places = 3 + + def write_coefficient(coefficient, is_initial_term=False): + coefficient = round(coefficient, num_decimal_places) + if coefficient % 1 == 0: + coefficient = int(coefficient) + + output = '' + + if coefficient < 0: + output += '-' + if coefficient > 0 and not is_initial_term: + output += '+' + + if not is_initial_term: + output += ' ' + + if abs(coefficient) != 1: + output += '{}'.format(abs(coefficient)) + + return output + + n = self.normal_vector + + try: + initial_index = Hyperplane.first_nonzero_index(n) + terms = [write_coefficient( + n[i], is_initial_term=(i == initial_index)) + + 'x_{}'.format(i + 1) + for i in range(self.dimension) + if round(n[i], num_decimal_places) != 0] + output = ' '.join(terms) + + except Exception as e: + if str(e) == self.NO_NONZERO_ELTS_FOUND_MSG: + output = '0' + else: + raise e + + constant = round(self.constant_term, num_decimal_places) + if constant % 1 == 0: + constant = int(constant) + output += ' = {}'.format(constant) + + return output + + def is_parallel(self, plane2): + return self.normal_vector.is_parallel(plane2.normal_vector) + + def __eq__(self, plane2): + if self.normal_vector.is_zero(): + if not plane2.normal_vector.is_zero(): + return False + + diff = self.constant_term - plane2.constant_term + return MyDecimal(diff).is_near_zero() + + elif plane2.normal_vector.is_zero(): + return False + + if not self.is_parallel(plane2): + return False + + basepoint_difference = self.basepoint.minus(plane2.basepoint) + return basepoint_difference.is_orthogonal(self.normal_vector) + + def __iter__(self): + self.current = 0 + return self + + def next(self): + if self.current >= len(self.normal_vector): + raise StopIteration + else: + current_value = self.normal_vector[self.current] + self.current += 1 + return current_value + + def __len__(self): + return len(self.normal_vector) + + def __getitem__(self, i): + return self.normal_vector[i] + + @staticmethod + def first_nonzero_index(iterable): + for k, item in enumerate(iterable): + if not MyDecimal(item).is_near_zero(): + return k + raise Exception(Hyperplane.NO_NONZERO_ELTS_FOUND_MSG) diff --git a/von_Udacity/line.py b/von_Udacity/line.py new file mode 100755 index 0000000..3d192b9 --- /dev/null +++ b/von_Udacity/line.py @@ -0,0 +1,169 @@ +from decimal import Decimal, getcontext +from vector import Vector + +getcontext().prec = 30 + + +class MyDecimal(Decimal): + def is_near_zero(self, eps=1e-10): + return abs(self) < eps + + +class Line(object): + + NO_NONZERO_ELTS_FOUND_MSG = 'No nonzero elements found' + + def __init__(self, normal_vector=None, constant_term=None): + self.dimension = 2 + + if not normal_vector: + all_zeros = ['0'] * self.dimension + normal_vector = Vector(all_zeros) + self.normal_vector = normal_vector + + if not constant_term: + constant_term = Decimal('0') + self.constant_term = Decimal(constant_term) + + self.set_basepoint() + + def set_basepoint(self): + try: + n = self.normal_vector + c = self.constant_term + basepoint_coords = ['0'] * self.dimension + + initial_index = Line.first_nonzero_index(n) + initial_coefficient = n[initial_index] + basepoint_coords[initial_index] = c / initial_coefficient + self.basepoint = Vector(basepoint_coords) + + except Exception as e: + if str(e) == Line.NO_NONZERO_ELTS_FOUND_MSG: + self.basepoint = None + else: + raise e + + def __str__(self): + + num_decimal_places = 3 + + def write_coefficient(coefficient, is_initial_term=False): + coefficient = round(coefficient, num_decimal_places) + if coefficient % 1 == 0: + coefficient = int(coefficient) + + output = '' + + if coefficient < 0: + output += '-' + if coefficient > 0 and not is_initial_term: + output += '+' + + if not is_initial_term: + output += ' ' + + if abs(coefficient) != 1: + output += '{}'.format(abs(coefficient)) + + return output + + n = self.normal_vector + + try: + initial_index = Line.first_nonzero_index(n) + terms = [write_coefficient(n[i], + is_initial_term=(i == initial_index)) + + 'x_{}'.format(i + 1) + for i in range(self.dimension) + if round(n[i], num_decimal_places) != 0] + output = ' '.join(terms) + + except Exception as e: + if str(e) == self.NO_NONZERO_ELTS_FOUND_MSG: + output = '0' + else: + raise e + + constant = round(self.constant_term, num_decimal_places) + if constant % 1 == 0: + constant = int(constant) + output += ' = {}'.format(constant) + + return output + + @staticmethod + def first_nonzero_index(iterable): + for k, item in enumerate(iterable): + if not MyDecimal(item).is_near_zero(): + return k + raise Exception(Line.NO_NONZERO_ELTS_FOUND_MSG) + + def is_parallel(self, line2): + return self.normal_vector.is_parallel(line2.normal_vector) + + def __eq__(self, line2): + if self.normal_vector.is_zero(): + if not line2.normal_vector.is_zero(): + return False + + diff = self.constant_term - line2.constant_term + return MyDecimal(diff).is_near_zero() + + elif line2.normal_vector.is_zero(): + return False + + if not self.is_parallel(line2): + return False + + basepoint_difference = self.basepoint.minus(line2.basepoint) + return basepoint_difference.is_orthogonal(self.normal_vector) + + def intersection(self, line2): + + a, b = self.normal_vector.coordinates + c, d = line2.normal_vector.coordinates + k1 = self.constant_term + k2 = line2.constant_term + denom = ((a * d) - (b * c)) + + if MyDecimal(denom).is_near_zero(): + if self == line2: + return self + else: + return None + + one_over_denom = Decimal('1') / ((a * d) - (b * c)) + x_num = (d * k1 - b * k2) + y_num = (-c * k1 + a * k2) + + return Vector([x_num, y_num]).times_scalar(one_over_denom) + + +# first system +# 4.046x + 2.836y = 1.21 +# 10.115x + 7.09y = 3.025 + +line1 = Line(Vector([4.046, 2.836]), 1.21) +line2 = Line(Vector([10.115, 7.09]), 3.025) + +print 'first system instersects in: {}'.format(line1.intersection(line2)) + + +# second system +# 7.204x + 3.182y = 8.68 +# 8.172x + 4.114y = 9.883 + +line3 = Line(Vector([7.204, 3.182]), 8.68) +line4 = Line(Vector([8.172, 4.114]), 9.883) + +print 'second system instersects in: {}'.format(line3.intersection(line4)) + +# third system +# 1.182x + 5.562y = 6.744 +# 1.773x + 8.343y = 9.525 + +line5 = Line(Vector([1.182, 5.562]), 6.744) +line6 = Line(Vector([1.773, 8.343]), 9.525) + +print 'third system instersects in: {}'.format(line5.intersection(line6)) diff --git a/von_Udacity/linear_system.py b/von_Udacity/linear_system.py new file mode 100755 index 0000000..cf1d0dd --- /dev/null +++ b/von_Udacity/linear_system.py @@ -0,0 +1,520 @@ +from copy import deepcopy +from decimal import Decimal, getcontext +from hyperplane import Hyperplane +from plane import Plane +from vector import Vector + +getcontext().prec = 30 + + +class MyDecimal(Decimal): + def is_near_zero(self, eps=1e-10): + return abs(self) < eps + + +def _get_new_plane(coefficient, plane): + new_normal_vector = plane.normal_vector.times_scalar(coefficient) + return Hyperplane(normal_vector=new_normal_vector, + constant_term=coefficient * plane.constant_term) + + +class LinearSystem(object): + + ALL_PLANES_MUST_BE_IN_SAME_DIM_MSG = ('All planes in the system should ' + 'live in the same dimension') + NO_SOLUTIONS_MSG = 'No solutions' + INF_SOLUTIONS_MSG = 'Infinitely many solutions' + + def __init__(self, planes): + try: + d = planes[0].dimension + for p in planes: + assert p.dimension == d + + self.planes = planes + self.dimension = d + + except AssertionError: + raise Exception(self.ALL_PLANES_MUST_BE_IN_SAME_DIM_MSG) + + def __len__(self): + return len(self.planes) + + def __getitem__(self, i): + return self.planes[i] + + def __setitem__(self, i, x): + try: + assert x.dimension == self.dimension + self.planes[i] = x + + except AssertionError: + raise Exception(self.ALL_PLANES_MUST_BE_IN_SAME_DIM_MSG) + + def __str__(self): + ret = 'Linear System:\n' + temp = ['Equation {}: {}'.format(i + 1, p) + for i, p in enumerate(self.planes)] + ret += '\n'.join(temp) + return ret + + def swap_rows(self, row1, row2): + self[row1], self[row2] = self[row2], self[row1] + + def multiply_coefficient_and_row(self, coefficient, row): + self[row] = _get_new_plane(coefficient, self[row]) + + def add_multiple_times_row_to_row(self, coefficient, row_to_add, + row_to_be_added_to): + recipient_plane = self[row_to_be_added_to] + new_plane = _get_new_plane(coefficient, self[row_to_add]) + new_normal_vector = \ + recipient_plane.normal_vector.plus(new_plane.normal_vector) + constant_term = new_plane.constant_term + recipient_plane.constant_term + self[row_to_be_added_to] = Hyperplane(normal_vector=new_normal_vector, + constant_term=constant_term) + + def indices_of_first_nonzero_terms_in_each_row(self): + num_equations = len(self) + + indices = [-1] * num_equations + + for i, p in enumerate(self.planes): + try: + p.first_nonzero_index(p.normal_vector) + indices[i] = p.first_nonzero_index(p.normal_vector) + except Exception as e: + if str(e) == Plane.NO_NONZERO_ELTS_FOUND_MSG: + continue + else: + raise e + + return indices + + def compute_triangular_form(self): + system = deepcopy(self) + + num_equations = len(system) + num_variables = system.dimension + + col = 0 + for row in range(num_equations): + while col < num_variables: + c = MyDecimal(system[row].normal_vector[col]) + if c.is_near_zero(): + swap_succeeded = system.did_swap_with_row_below(row, col) + if not swap_succeeded: + col += 1 + continue + + system.clear_coefficients_bellow(row, col) + col += 1 + break + + return system + + def did_swap_with_row_below(self, row, col): + num_equations = len(self) + + for k in range(row + 1, num_equations): + coefficient = MyDecimal(self[k].normal_vector[col]) + if not coefficient.is_near_zero(): + self.swap_rows(row, k) + return True + + return False + + def clear_coefficients_bellow(self, row, col): + num_equations = len(self) + beta = MyDecimal(self[row].normal_vector[col]) + + for row_to_be_added_to in range(row + 1, num_equations): + n = self[row_to_be_added_to].normal_vector + gamma = n[col] + alpha = -gamma / beta + self.add_multiple_times_row_to_row(alpha, row, row_to_be_added_to) + + def clear_coefficients_above(self, row, col): + for row_to_be_added_to in range(row)[::-1]: + n = self[row_to_be_added_to].normal_vector + alpha = -(n[col]) + self.add_multiple_times_row_to_row(alpha, row, row_to_be_added_to) + + def compute_rref(self): + tf = self.compute_triangular_form() + + num_equations = len(tf) + pivot_indices = tf.indices_of_first_nonzero_terms_in_each_row() + + for row in range(num_equations)[::-1]: + pivot_var = pivot_indices[row] + if pivot_var < 0: + continue + tf.scale_row_to_make_coefficient_equal_one(row, pivot_var) + tf.clear_coefficients_above(row, pivot_var) + + return tf + + def scale_row_to_make_coefficient_equal_one(self, row, col): + n = self[row].normal_vector + beta = Decimal('1.0') / n[col] + self.multiply_coefficient_and_row(beta, row) + + def do_gaussian_elimination(self): + rref = self.compute_rref() + + try: + rref.raise_excepion_if_contradictory_equation() + rref.raise_excepion_if_too_few_pivots() + except Exception as e: + return e.message + + num_variables = rref.dimension + solution_coordinates = [rref.planes[i].constant_term + for i in range(num_variables)] + + return Vector(solution_coordinates) + + def raise_excepion_if_contradictory_equation(self): + for plane in self.planes: + try: + plane.first_nonzero_index(plane.normal_vector) + + except Exception as e: + if str(e) == 'No nonzero elements found': + constant_term = MyDecimal(plane.constant_term) + if not constant_term.is_near_zero(): + raise Exception(self.NO_SOLUTIONS_MSG) + + else: + raise e + + def raise_excepion_if_too_few_pivots(self): + pivot_indices = self.indices_of_first_nonzero_terms_in_each_row() + num_pivots = sum([1 if index >= 0 else 0 for index in pivot_indices]) + num_variables = self.dimension + + if num_pivots < num_variables: + raise Exception(self.INF_SOLUTIONS_MSG) + + def compute_solution(self): + try: + return self.do_gaussian_elimination_and_parametrization() + + except Exception as e: + if str(e) == self.NO_SOLUTIONS_MSG: + return str(e) + else: + raise e + + def do_gaussian_elimination_and_parametrization(self): + rref = self.compute_rref() + rref.raise_excepion_if_contradictory_equation() + + direction_vectors = rref.extract_direction_vectors_for_parametrization() # NOQA + basepoint = rref.extract_basepoint_for_parametrization() + + return Parametrization(basepoint, direction_vectors) + + def extract_direction_vectors_for_parametrization(self): + num_variables = self.dimension + pivot_indices = self.indices_of_first_nonzero_terms_in_each_row() + free_variable_indices = set(range(num_variables)) - set(pivot_indices) + + direction_vectors = [] + + for free_var in free_variable_indices: + vector_coords = [0] * num_variables + vector_coords[free_var] = 1 + for index, plane in enumerate(self.planes): + pivot_var = pivot_indices[index] + if pivot_var < 0: + break + vector_coords[pivot_var] = -plane.normal_vector[free_var] + + direction_vectors.append(Vector(vector_coords)) + + return direction_vectors + + def extract_basepoint_for_parametrization(self): + num_variables = self.dimension + pivot_indices = self.indices_of_first_nonzero_terms_in_each_row() + + basepoint_coords = [0] * num_variables + + for index, plane in enumerate(self.planes): + pivot_var = pivot_indices[index] + if pivot_var < 0: + break + basepoint_coords[pivot_var] = plane.constant_term + + return Vector(basepoint_coords) + + +p0 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p1 = Plane(normal_vector=Vector(['0', '1', '0']), constant_term='2') +p2 = Plane(normal_vector=Vector(['1', '1', '-1']), constant_term='3') +p3 = Plane(normal_vector=Vector(['1', '0', '-2']), constant_term='2') + +s = LinearSystem([p0, p1, p2, p3]) +# Print initial system +# print s + +p0 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p1 = Plane(normal_vector=Vector(['0', '1', '0']), constant_term='2') +p2 = Plane(normal_vector=Vector(['1', '1', '-1']), constant_term='3') +p3 = Plane(normal_vector=Vector(['1', '0', '-2']), constant_term='2') + +s = LinearSystem([p0, p1, p2, p3]) +s.swap_rows(0, 1) +if not (s[0] == p1 and s[1] == p0 and s[2] == p2 and s[3] == p3): + print 'test case 1 failed' + +s.swap_rows(1, 3) +if not (s[0] == p1 and s[1] == p3 and s[2] == p2 and s[3] == p0): + print 'test case 2 failed' + +s.swap_rows(3, 1) +if not (s[0] == p1 and s[1] == p0 and s[2] == p2 and s[3] == p3): + print 'test case 3 failed' + +s.multiply_coefficient_and_row(1, 0) +if not (s[0] == p1 and s[1] == p0 and s[2] == p2 and s[3] == p3): + print 'test case 4 failed' + +s.multiply_coefficient_and_row(-1, 2) +new_s2 = Plane(normal_vector=Vector(['-1', '-1', '1']), constant_term='-3') +if not (s[0] == p1 and + s[1] == p0 and + s[2] == new_s2 and + s[3] == p3): + print 'test case 5 failed' + +s.multiply_coefficient_and_row(10, 1) +new_s1 = Plane(normal_vector=Vector(['10', '10', '10']), constant_term='10') +if not (s[0] == p1 and + s[1] == new_s1 and + s[2] == new_s2 and + s[3] == p3): + print 'test case 6 failed' + +s.add_multiple_times_row_to_row(0, 0, 1) +if not (s[0] == p1 and + s[1] == new_s1 and + s[2] == new_s2 and + s[3] == p3): + print 'test case 7 failed' + +s.add_multiple_times_row_to_row(1, 0, 1) +added_s1 = Plane(normal_vector=Vector(['10', '11', '10']), constant_term='12') +if not (s[0] == p1 and + s[1] == added_s1 and + s[2] == new_s2 and + s[3] == p3): + print 'test case 8 failed' + +s.add_multiple_times_row_to_row(-1, 1, 0) +new_s0 = Plane(normal_vector=Vector(['-10', '-10', '-10']), + constant_term='-10') +if not (s[0] == new_s0 and + s[1] == added_s1 and + s[2] == new_s2 and + s[3] == p3): + print 'test case 9 failed' + +p1 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['0', '1', '1']), constant_term='2') +s = LinearSystem([p1, p2]) +t = s.compute_triangular_form() +if not (t[0] == p1 and + t[1] == p2): + print 'test case 1 failed' + +p1 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='2') +s = LinearSystem([p1, p2]) +t = s.compute_triangular_form() +if not (t[0] == p1 and + t[1] == Plane(constant_term='1')): + print 'test case 2 failed' + +p1 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['0', '1', '0']), constant_term='2') +p3 = Plane(normal_vector=Vector(['1', '1', '-1']), constant_term='3') +p4 = Plane(normal_vector=Vector(['1', '0', '-2']), constant_term='2') +s = LinearSystem([p1, p2, p3, p4]) +t = s.compute_triangular_form() +if not (t[0] == p1 and + t[1] == p2 and + t[2] == Plane(normal_vector=Vector(['0', '0', '-2']), + constant_term='2') and + t[3] == Plane()): + print 'test case 3 failed' + +p1 = Plane(normal_vector=Vector(['0', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['1', '-1', '1']), constant_term='2') +p3 = Plane(normal_vector=Vector(['1', '2', '-5']), constant_term='3') +s = LinearSystem([p1, p2, p3]) +t = s.compute_triangular_form() +if not (t[0] == Plane(normal_vector=Vector(['1', '-1', '1']), + constant_term='2') and + t[1] == Plane(normal_vector=Vector(['0', '1', '1']), + constant_term='1') and + t[2] == Plane(normal_vector=Vector(['0', '0', '-9']), + constant_term='-2')): + print 'test case 4 failed' + + +# *************** + +p1 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['0', '1', '1']), constant_term='2') +s = LinearSystem([p1, p2]) +r = s.compute_rref() +if not (r[0] == Plane(normal_vector=Vector(['1', '0', '0']), + constant_term='-1') and + r[1] == p2): + print 'test case 1 failed' + +p1 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='2') +s = LinearSystem([p1, p2]) +r = s.compute_rref() +if not (r[0] == p1 and + r[1] == Plane(constant_term='1')): + print 'test case 2 failed' + +p1 = Plane(normal_vector=Vector(['1', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['0', '1', '0']), constant_term='2') +p3 = Plane(normal_vector=Vector(['1', '1', '-1']), constant_term='3') +p4 = Plane(normal_vector=Vector(['1', '0', '-2']), constant_term='2') +s = LinearSystem([p1, p2, p3, p4]) +r = s.compute_rref() +if not (r[0] == Plane(normal_vector=Vector(['1', '0', '0']), + constant_term='0') and + r[1] == p2 and + r[2] == Plane(normal_vector=Vector(['0', '0', '-2']), + constant_term='2') and + r[3] == Plane()): + print 'test case 3 failed' + +p1 = Plane(normal_vector=Vector(['0', '1', '1']), constant_term='1') +p2 = Plane(normal_vector=Vector(['1', '-1', '1']), constant_term='2') +p3 = Plane(normal_vector=Vector(['1', '2', '-5']), constant_term='3') +s = LinearSystem([p1, p2, p3]) +r = s.compute_rref() +if not (r[0] == Plane(normal_vector=Vector(['1', '0', '0']), + constant_term=Decimal('23') / Decimal('9')) and + r[1] == Plane(normal_vector=Vector(['0', '1', '0']), + constant_term=Decimal('7') / Decimal('9')) and + r[2] == Plane(normal_vector=Vector(['0', '0', '1']), + constant_term=Decimal('2') / Decimal('9'))): + print 'test case 4 failed' + + +# first system +p1 = Plane(Vector([5.862, 1.178, -10.366]), -8.15) +p2 = Plane(Vector([-2.931, -0.589, 5.183]), -4.075) +system1 = LinearSystem([p1, p2]) +print 'first system: {}'.format(system1.do_gaussian_elimination()) + + +# # second system +p1 = Plane(Vector([8.631, 5.112, -1.816]), -5.113) +p2 = Plane(Vector([4.315, 11.132, -5.27]), -6.775) +p3 = Plane(Vector([-2.158, 3.01, -1.727]), -0.831) +system2 = LinearSystem([p1, p2, p3]) +print 'second system: {}'.format(system2.do_gaussian_elimination()) + + +# third system +p1 = Plane(Vector([5.262, 2.739, -9.878]), -3.441) +p2 = Plane(Vector([5.111, 6.358, 7.638]), -2.152) +p3 = Plane(Vector([2.016, -9.924, -1.367]), -9.278) +p4 = Plane(Vector([2.167, -13.543, -18.883]), -10.567) +system3 = LinearSystem([p1, p2, p3, p4]) +print 'thrid system: {} '.format(system3.do_gaussian_elimination()) + + +class Parametrization(object): + + BASEPT_AND_DIR_VECTORS_MUST_BE_IN_SAME_DIM = ( + 'The basepoint and direction vectors should all live in the same ' + 'dimension') + + def __init__(self, basepoint, direction_vectors): + + self.basepoint = basepoint + self.direction_vectors = direction_vectors + self.dimension = self.basepoint.dimension + + try: + for v in direction_vectors: + assert v.dimension == self.dimension + + except AssertionError: + raise Exception(self.BASEPT_AND_DIR_VECTORS_MUST_BE_IN_SAME_DIM) + + def __str__(self): + + output = '' + for coord in range(self.dimension): + output += 'x_{} = {} '.format(coord + 1, + round(self.basepoint[coord], 3)) + for free_var, vector in enumerate(self.direction_vectors): + output += '+ {} t_{}'.format(round(vector[coord], 3), + free_var + 1) + output += '\n' + return output + + +p1 = Plane(normal_vector=Vector([0.786, 0.786, 0.588]), constant_term=-0.714) +p2 = Plane(normal_vector=Vector([-0.131, -0.131, 0.244]), constant_term=0.319) + +system = LinearSystem([p1, p2]) +print system.compute_solution() + + +p1 = Plane(Vector([8.631, 5.112, -1.816]), -5.113) +p2 = Plane(Vector([4.315, 11.132, -5.27]), -6.775) +p3 = Plane(Vector([-2.158, 3.01, -1.727]), -0.831) + +system = LinearSystem([p1, p2, p3]) +print system.compute_solution() + +p1 = Plane(Vector([0.935, 1.76, -9.365]), -9.955) +p2 = Plane(Vector([0.187, 0.352, -1.873]), -1.991) +p3 = Plane(Vector([0.374, 0.704, -3.746]), -3.982) +p4 = Plane(Vector([-0.561, -1.056, 5.619]), 5.973) + +print system.compute_solution() + + +# The systems bellow are just to test hyperplanes + +p1 = Hyperplane(normal_vector=Vector([0.786, 0.786]), constant_term=0.786) +p2 = Hyperplane(normal_vector=Vector([-0.131, -0.131]), constant_term=-0.131) + +system = LinearSystem([p1, p2]) +print system.compute_solution() + + +p1 = Hyperplane(normal_vector=Vector([2.102, 7.489, -0.786]), + constant_term=-5.113) +p2 = Hyperplane(normal_vector=Vector([-1.131, 8.318, -1.209]), + constant_term=-6.775) +p3 = Hyperplane(normal_vector=Vector([9.015, 5.873, -1.105]), + constant_term=-0.831) + +system = LinearSystem([p1, p2, p3]) +print system.compute_solution() + +p1 = Hyperplane(normal_vector=Vector([0.786, 0.786, 8.123, 1.111, -8.363]), + constant_term=-9.955) +p2 = Hyperplane(normal_vector=Vector([0.131, -0.131, 7.05, -2.813, 1.19]), + constant_term=-1.991) +p3 = Hyperplane(normal_vector=Vector([9.015, -5.873, -1.105, 2.013, -2.802]), + constant_term=-3.982) + +system = LinearSystem([p1, p2, p3]) +print system.compute_solution() diff --git a/von_Udacity/plane.py b/von_Udacity/plane.py new file mode 100755 index 0000000..62dee0b --- /dev/null +++ b/von_Udacity/plane.py @@ -0,0 +1,171 @@ +from decimal import Decimal, getcontext +from vector import Vector + +getcontext().prec = 30 + + +class MyDecimal(Decimal): + def is_near_zero(self, eps=1e-10): + return abs(self) < eps + + +class Plane(object): + + NO_NONZERO_ELTS_FOUND_MSG = 'No nonzero elements found' + + def __init__(self, normal_vector=None, constant_term=None): + self.dimension = 3 + + if not normal_vector: + all_zeros = ['0'] * self.dimension + normal_vector = Vector(all_zeros) + self.normal_vector = normal_vector + + if not constant_term: + constant_term = Decimal('0') + self.constant_term = Decimal(constant_term) + + self.set_basepoint() + + def set_basepoint(self): + try: + n = self.normal_vector + c = self.constant_term + basepoint_coords = ['0'] * self.dimension + + initial_index = Plane.first_nonzero_index(n) + initial_coefficient = n[initial_index] + + basepoint_coords[initial_index] = c / initial_coefficient + self.basepoint = Vector(basepoint_coords) + + except Exception as e: + if str(e) == Plane.NO_NONZERO_ELTS_FOUND_MSG: + self.basepoint = None + else: + raise e + + def __str__(self): + + num_decimal_places = 3 + + def write_coefficient(coefficient, is_initial_term=False): + coefficient = round(coefficient, num_decimal_places) + if coefficient % 1 == 0: + coefficient = int(coefficient) + + output = '' + + if coefficient < 0: + output += '-' + if coefficient > 0 and not is_initial_term: + output += '+' + + if not is_initial_term: + output += ' ' + + if abs(coefficient) != 1: + output += '{}'.format(abs(coefficient)) + + return output + + n = self.normal_vector + + try: + initial_index = Plane.first_nonzero_index(n) + terms = [write_coefficient( + n[i], is_initial_term=(i == initial_index)) + + 'x_{}'.format(i + 1) + for i in range(self.dimension) + if round(n[i], num_decimal_places) != 0] + output = ' '.join(terms) + + except Exception as e: + if str(e) == self.NO_NONZERO_ELTS_FOUND_MSG: + output = '0' + else: + raise e + + constant = round(self.constant_term, num_decimal_places) + if constant % 1 == 0: + constant = int(constant) + output += ' = {}'.format(constant) + + return output + + def is_parallel(self, plane2): + return self.normal_vector.is_parallel(plane2.normal_vector) + + def __eq__(self, plane2): + if self.normal_vector.is_zero(): + if not plane2.normal_vector.is_zero(): + return False + + diff = self.constant_term - plane2.constant_term + return MyDecimal(diff).is_near_zero() + + elif plane2.normal_vector.is_zero(): + return False + + if not self.is_parallel(plane2): + return False + + basepoint_difference = self.basepoint.minus(plane2.basepoint) + return basepoint_difference.is_orthogonal(self.normal_vector) + + def __iter__(self): + self.current = 0 + return self + + def next(self): + if self.current >= len(self.normal_vector): + raise StopIteration + else: + current_value = self.normal_vector[self.current] + self.current += 1 + return current_value + + def __len__(self): + return len(self.normal_vector) + + def __getitem__(self, i): + return self.normal_vector[i] + + @staticmethod + def first_nonzero_index(iterable): + for k, item in enumerate(iterable): + if not MyDecimal(item).is_near_zero(): + return k + raise Exception(Plane.NO_NONZERO_ELTS_FOUND_MSG) + + +if __name__ == '__main__': + # first system of planes: + # -0.412x + 3.806y + 0.728z = -3.46 + # 1.03x - 9.515y - 1.82z = 8.65 + + plane1 = Plane(Vector([-0.412, 3.806, 0.728]), -3.46) + plane2 = Plane(Vector([1.03, -9.515, -1.82]), 8.65) + + print '1 is parallel: {}'.format(plane1.is_parallel(plane2)) + print '1 is equal: {}'.format(plane1 == plane2) + + # second system of planes: + # 2.611x + 5.528y + 0.283z = 4.6 + # 7.715x + 8.306y + 5.342z = 3.76 + + plane3 = Plane(Vector([2.611, 5.528, 0.283]), 4.6) + plane4 = Plane(Vector([7.715, 8.306, 5.342]), 3.76) + + print '2 is parallel: {}'.format(plane3.is_parallel(plane4)) + print '2 is equal: {}'.format(plane3 == plane4) + + # third system of planes: + # -7.926x + 8.625y - 7.212z = -7.952 + # -2.642x + 2.875y - 2.404z = -2.443 + + plane5 = Plane(Vector([-7.926, 8.625, -7.212]), -7.952) + plane6 = Plane(Vector([-2.642, 2.875, -2.404]), -2.443) + + print '3 is parallel: {}'.format(plane5.is_parallel(plane6)) + print '3 is equal: {}'.format(plane5 == plane6) diff --git a/von_Udacity/summary.py b/von_Udacity/summary.py new file mode 100755 index 0000000..3946ee0 --- /dev/null +++ b/von_Udacity/summary.py @@ -0,0 +1,204 @@ +# Purpose: Combining the information learned from the course linear algebra into an interactive code which can be +# used to test equations easily. + +from vector import Vector +from line import Line +from linear_system import LinearSystem +from plane import Plane + + +def error_check_input(text, lower=-99999999999, upper=99999999999): + while True: + try: + operation = float(input(text)) + if operation >= lower and operation <= upper: + break + else: + print("Only enter one of the specified options") + except ValueError: + print("You may only enter a number") + return operation + + +# Function creates a vector using user input +def enter_vector(times=3): + coordinates = [] + for x in range(times): + coordinates.append(float(input("Enter coordinate: "))) + return Vector(coordinates) + + +def enter_vectors(vectors=2, choice=True, dimension=None): + if choice: + dimension = int(error_check_input("How many dimensions would you like your vectors to include(2-5): ", 2, 5)) + my_vectors = [] + for i in range(vectors): + print("Entering Vector #%s:" % str(i + 1)) + v = enter_vector(dimension) + my_vectors.append(v) + return my_vectors + + +def vector(): + print("Vector is a quantity with both magnitude and direction\n" + "Using vectors, we may perform mathematical operations including: \n" + "\t1) Vector addition\n" + "\t2) Vector Subtraction\n " + "\t3) Multiply the vector by a scalar\n" + "\t4) Find one vector's magnitude\n" + "\t5) Normalize a vector\n" + "\t6) Calculate dot product \n" + "\t7) Compute angle\n" + "\t8) Determine whether parallel\n" + "\t9) Determine whether orthogonal\n" + "\t10) Determine the projected vector\n" + "\t11) Determine the orthogonal vector\n" + "\t12) Determine the cross product\n" + "\t13) Find the area of the parallelogram formed\n" + "\t14) Find the area of the triangle formed\n") + option_to_execute = error_check_input("Pick one of the above operations to perform on vectors: ", 1, 14) + print("===================================================") + if option_to_execute == 1: + v, w = enter_vectors(2) + addition = v.plus(w) + print("Together the vectors form: %s" % addition) + elif option_to_execute == 2: + v, w = enter_vectors(2) + subtraction = v.minus(w) + print("Subtracting the vectors forms: %s" % subtraction) + elif option_to_execute == 3: + v = enter_vectors(1)[0] + scalar = error_check_input("Enter scalar: ") + multiplied = v.times_scalar(scalar) + print("Multiplying the vector by scalar returns: %s" % multiplied) + elif option_to_execute == 4: + v = enter_vectors(1)[0] + magnitude = v.magnitude() + print("The magnitude of the vector is: %s" % magnitude) + elif option_to_execute == 5: + v = enter_vectors(1)[0] + normalized = v.normalize() + print("The vector normalized is: %s" % normalized) + elif option_to_execute == 6: + v, w = enter_vectors(2) + dot_product = v.dot_product(w) + print('dot_product: {}'.format(round(dot_product, 3))) + elif option_to_execute == 7: + v, w = enter_vectors(2) + angle_degrees = v.get_angle_deg(w) + angle_radiants = v.get_angle_rad(w) + print("The first angle is:") + print("%s in radiant and %s in degrees " % (angle_degrees, angle_radiants)) + elif option_to_execute == 8: + v, w = enter_vectors(2) + is_parallel = v.is_parallel(w) + if is_parallel: + print("The vectors are parallel") + else: + print("The vectors aren't parallel") + elif option_to_execute == 9: + v, w = enter_vectors(2) + is_orthogonal = v.is_orthogonal(w) + if is_orthogonal: + print("The vectors are orthogonal") + else: + print("The vectors aren't orthogonal") + elif option_to_execute == 10: + v, w = enter_vectors(2) + projected_vector = v.get_projected_vector(w) + print("The projected vector is: %s" % projected_vector) + elif option_to_execute == 11: + v, w = enter_vectors(2) + orthogonal_vector = v.get_orthogonal_vector(w) + print("The orthogonal vector: %s" % orthogonal_vector) + elif option_to_execute == 12: + v, w = enter_vectors(2, False, 3) + cross_product = v.cross_product(w) + print("The cross product is: %s" % cross_product) + elif option_to_execute == 13: + v, w = enter_vectors(2, False, 3) + area_parallelogram = v.area_parallelogram(w) + print("The parallelogram formed area is: %s" % area_parallelogram) + elif option_to_execute == 14: + v, w = enter_vectors(2, False, 3) + area_triangle = v.area_triangle(w) + print("The triangle formed area is: %s" % area_triangle) + + +def enter_lines(times): + current_lines = [] + for i in range(times): + print("Enter the two-dimensional vector below:") + vec = enter_vector(2) + base_point = error_check_input("Enter the base point: ") + current_lines.append(Line(vec, base_point)) + return current_lines + + +def lines(): + print("A line can be defined by a basepoint and a direction vector. ") + line1, line2 = enter_lines(2) + print("The lines intersect in %s" % line1.intersection(line2)) + + +def linear_system(): + print("Linear systems use gaussian elimination to\n" + "\t1) Determine amount of possible solutions for planes\n" + "\t2) Solve a system of hyperplanes \n") + option_to_execute = error_check_input("Pick one of the above: ", 1, 2) + if option_to_execute == 1: + total_planes = int(error_check_input("Enter amount of planes(2/5): ", 2, 5)) + myPlanes = enter_plane(total_planes) + system1 = LinearSystem([i for i in myPlanes]) + print("System intersection: ", system1.do_gaussian_elimination()) + elif option_to_execute == 2: + total_planes = int(error_check_input("Enter amount of planes(2/5): ", 2, 5)) + myPlanes = enter_plane(total_planes) + system1 = LinearSystem([i for i in myPlanes]) + print(system1.compute_solution()) + + +def enter_plane(times): + planes = [] + for i in range(times): + print("Plane #%s" % str(i + 1)) + print("First, enter a 3 dimensional vector: ") + vector = enter_vector(3) + base_point = error_check_input("Now, enter base point: ") + planes.append(Plane(vector, base_point)) + return planes + + +def plane(): + print("A plane is similar to line, but in three dimensions\nLet's check whether two planes are parallel/equal: ") + p1, p2 = enter_plane(2) + if p1.is_parallel(p2): + print("The planes are parallel") + else: + print("The planes are not parallel") + if p1 == p2: + print("The planes are equal") + else: + print("The planes are not equal") + + +print("Welcome to the interactive math problem solving library. In this library you will find solutions to well-known" + "\nmath problems such as vectors, intersection, linear systems and planes.") +while True: + print("\n" + "\t1) Vectors\n" + "\t2) Intersection of lines\n" + "\t3) Linear Systems\n" + "\t4) Planes\n" + "\t5) Exit Program") + section_view = error_check_input("Choose the operation you would like to execute from the above: ", 1, 4) + if section_view == 1: + vector() + elif section_view == 2: + lines() + elif section_view == 3: + linear_system() + elif section_view == 4: + plane() + elif section_view == 5: + break diff --git a/von_Udacity/vector.py b/von_Udacity/vector.py new file mode 100755 index 0000000..0225600 --- /dev/null +++ b/von_Udacity/vector.py @@ -0,0 +1,247 @@ +from math import acos, sqrt, pi +from decimal import Decimal, getcontext + +getcontext().prec = 30 + + +class MyDecimal(Decimal): + def is_near_zero(self, eps=1e-10): + return abs(self) < eps + + +class Vector(object): + def __init__(self, coordinates): + try: + if not coordinates: + raise ValueError + self.coordinates = tuple([Decimal(c) for c in coordinates]) + self.dimension = len(coordinates) + + except ValueError: + raise ValueError('The coordinates must be nonempty') + + except TypeError: + raise TypeError('The coordinates must be an iterable') + + def __iter__(self): + self.current = 0 + return self + + def __next__(self): + if self.current >= len(self.coordinates): + raise StopIteration + else: + current_value = self.coordinates[self.current] + self.current += 1 + return current_value + + def next(self): + return self.__next__() + + def __len__(self): + return len(self.coordinates) + + def __getitem__(self, i): + return self.coordinates[i] + + def __str__(self): + return 'Vector: {}'.format([round(coord, 3) + for coord in self.coordinates]) + + def __eq__(self, v): + return self.coordinates == v.coordinates + + def is_zero(self): + return set(self.coordinates) == set([Decimal(0)]) + + def plus(self, other): + new_coordinates = [x + y for x, y in zip(self.coordinates, other.coordinates)] + return Vector(new_coordinates) + + def minus(self, other): + return Vector([coords[0] - coords[1] + for coords in zip(self.coordinates, other.coordinates)]) + + def times_scalar(self, factor): + return Vector([Decimal(factor) * coord for coord in self.coordinates]) + + def magnitude(self): + return Decimal(sqrt(sum([coord * coord + for coord in self.coordinates]))) + + def normalize(self): + try: + return self.times_scalar(Decimal('1.0') / self.magnitude()) + except ZeroDivisionError: + raise Exception('Cannot normalize the zero vector') + + def dot_product(self, other): + return sum(x * y for x, y in zip(self.coordinates, other.coordinates)) + + def get_angle_rad(self, other): + dot_prod = round(self.normalize().dot_product(other.normalize()), 3) + return acos(dot_prod) + + def get_angle_deg(self, other): + degrees_per_rad = 180. / pi + return degrees_per_rad * self.get_angle_rad(other) + + def is_parallel(self, other): + return (self.is_zero() or other.is_zero() or + self.get_angle_rad(other) in [0, pi]) + + def is_orthogonal(self, other): + return round(self.dot_product(other), 3) == 0 + + def get_projected_vector(self, other): + """ + Gets projection of vector v in b + """ + b_normalized = other.normalize() + return b_normalized.times_scalar(self.dot_product(b_normalized)) + + def get_orthogonal_vector(self, other): + return self.minus(self.get_projected_vector(other)) + + def cross_product(self, other): + [x1, y1, z1] = self.coordinates + [x2, y2, z2] = other.coordinates + x = (y1 * z2) - (y2 * z1) + y = -((x1 * z2) - (x2 * z1)) + z = (x1 * y2) - (x2 * y1) + return Vector([x, y, z]) + + def area_parallelogram(self, other): + return self.cross_product(other).magnitude() + + def area_triangle(self, other): + return self.cross_product(other).magnitude() / 2 + +if __name__ == '__main__': + v = Vector([8.218, -9.341]) + w = Vector([-1.129, 2.111]) + addition = v.plus(w) + print('addition: {}'.format(addition)) + + v = Vector([7.119, 8.215]) + w = Vector([-8.223, 0.878]) + subtraction = v.minus(w) + print('subtraction: {}'.format(subtraction)) + + v = Vector([1.671, -1.012, -0.318]) + multiplication = v.times_scalar(7.41) + print('multiplication: {}'.format(multiplication)) + + # ***************** + + v = Vector([-0.221, 7.437]) + first_magintude = v.magnitude() + print('first_magintude: {}'.format(round(first_magintude, 3))) + + v = Vector([8.813, -1.331, -6.247]) + second_magintude = v.magnitude() + print('second_magintude: {}'.format(round(second_magintude, 3))) + + v = Vector([5.581, -2.136]) + first_normalization = v.normalize() + print('first_normailization: {}'.format(first_normalization)) + + v = Vector([1.996, 3.108, -4.554]) + second_normalization = v.normalize() + print('second_normailization: {}'.format(second_normalization)) + + # ***************** + + v = Vector([7.887, 4.138]) + w = Vector([-8.802, 6.776]) + dot_product = v.dot_product(w) + print('first_dot_product: {}'.format(round(dot_product, 3))) + + v = Vector([-5.955, -4.904, -1.874]) + w = Vector([-4.496, -8.755, 7.103]) + dot_product = v.dot_product(w) + print('second_dot_product: {}'.format(round(dot_product, 3))) + + # ***************** + + v = Vector([3.183, -7.627]) + w = Vector([-2.668, 5.319]) + angle_rads = v.get_angle_rad(w) + print('first_angle_rads: {}'.format(angle_rads)) + + v = Vector([7.35, 0.221, 5.188]) + w = Vector([2.751, 8.259, 3.985]) + angle_degrees = v.get_angle_deg(w) + print('first_angle_rads: {}'.format(angle_degrees)) + + # ***************** + + v = Vector([-7.579, -7.88]) + w = Vector([22.737, 23.64]) + is_parallel = v.is_parallel(w) + is_orthogonal = v.is_orthogonal(w) + + print('1 parallel: {}, orthogonal: {}'.format(is_parallel, is_orthogonal)) + + v = Vector([-2.029, 9.97, 4.172]) + w = Vector([-9.231, -6.639, -7.245]) + is_parallel = v.is_parallel(w) + is_orthogonal = v.is_orthogonal(w) + + print('2 parallel: {}, orthogonal: {}'.format(is_parallel, is_orthogonal)) + + v = Vector([-2.328, -7.284, -1.214]) + w = Vector([-1.821, 1.072, -2.94]) + is_parallel = v.is_parallel(w) + is_orthogonal = v.is_orthogonal(w) + print('3 parallel: {}, orthogonal: {}'.format(is_parallel, is_orthogonal)) + + v = Vector([2.118, 4.827]) + w = Vector([0, 0]) + is_parallel = v.is_parallel(w) + is_orthogonal = v.is_orthogonal(w) + + print('4 parallel: {}, orthogonal: {}'.format(is_parallel, is_orthogonal)) + + # ***************** + + v = Vector([3.039, 1.879]) + w = Vector([0.825, 2.036]) + projected_vector = v.get_projected_vector(w) + + print('projected vector is: {}'.format(projected_vector)) + + v = Vector([-9.88, -3.264, -8.159]) + w = Vector([-2.155, -9.353, -9.473]) + orthogonal_vector = v.get_orthogonal_vector(w) + + print('orthogonal vector is: {}'.format(orthogonal_vector)) + + v = Vector([3.009, -6.172, 3.692, -2.51]) + w = Vector([6.404, -9.144, 2.759, 8.718]) + projected_vector = v.get_projected_vector(w) + orthogonal_vector = v.get_orthogonal_vector(w) + + print('second projected vector is: {}'.format(projected_vector)) + + print('second orthogonal vector is: {}'.format(orthogonal_vector)) + + # ***************** + + v1 = Vector([8.462, 7.893, -8.187]) + w1 = Vector([6.984, -5.975, 4.778]) + + v2 = Vector([-8.987, -9.838, 5.031]) + w2 = Vector([-4.268, -1.861, -8.866]) + + v3 = Vector([1.5, 9.547, 3.691]) + w3 = Vector([-6.007, 0.124, 5.772]) + + first_cross_product = v1.cross_product(w1) + print('cross product is: {}'.format(first_cross_product)) + + area_parallelogram = v2.area_parallelogram(w2) + print('area parallelogram is: {}'.format(round(area_parallelogram, 3))) + + area_triangle = v3.area_triangle(w3) + print('area triangle is: {}'.format(round(area_triangle, 3)))