Initial import.

This commit is contained in:
Kyle Isom 2022-03-07 23:47:39 -08:00
commit c0f599e563
9 changed files with 652 additions and 0 deletions

15
LICENSE Normal file
View File

@ -0,0 +1,15 @@
ISC License (ISC)
Copyright 2022 Kyle Isom <kyle@imap.cc>
Permission to use, copy, modify, and/or distribute this software for any
purpose with or without fee is hereby granted, provided that the above
copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

5
README.rst Normal file
View File

@ -0,0 +1,5 @@
kmath
=====
While learning maths, I've found it help to write code to explore the ideas. If I can implement them, it helps me understand them.

0
kmath/__init__.py Normal file
View File

308
kmath/linea.py Normal file
View File

@ -0,0 +1,308 @@
# -*- coding: utf-8 -*-
"""
``linea.vector``
Arbitrarily-sized vectors built on numpy's ndarray.
Components:
+ class Vector
+ exception NonConformantVectors
+ function: dot
+ function: angle
+ function: parallel
+ function: orthogonal
"""
# pylint: disable=C0103
import math
import numpy
EQUALITY_TOLERANCE = 0.001
class NonConformantVectors(Exception):
"""
A NonConformantVector is thrown when attempting to do operations
with vectors of differing sizes; generally, the size of the
/right-hand/ argument doesn't conform to the size of the /left-hand/
argument.
"""
def __init__(self, expected, actual):
"""
Initialize a new NonConformantVectors exception.
:param expected: length of the first vector
:param actual: length of the second vector
"""
self.expected = expected
self.actual = actual
Exception.__init__(self)
def __str__(self):
msg = "Expected the right-hand vector to have dimension {}, "
msg += "but it has a dimension of {}."
return msg.format(self.expected, self.actual)
class Vector:
"""
A vector is a one-dimensional vector of some arbitrary size. This size
is fixed and can't be changed later in the Vector's life.
"""
def __init__(self, a=None, *args):
"""
Initialise a vector, either using an iterable passed in or as a sequence
of values.
>>> print(Vector(1, 2, 3))
[1; 2; 3]
>>> print(Vector([1, 2, 3]))
[1; 2; 3]
"""
if len(args) > 0:
a = [a]
a.extend(args)
self.v = numpy.array(a)
else:
if isinstance(a, numpy.ndarray):
self.v = a
else:
self.v = numpy.array(a)
def __str__(self):
s = '['
for i in range(len(self.v)):
s += str(self.v[i])
if i < len(self.v) - 1:
s += '; '
s += ']'
return s
def __len__(self):
return len(self.v)
def __iter__(self):
return (x for x in self.v)
def __mul__(self, other):
o = list(map(lambda x: x * other, self.v))
return Vector(a=o)
def __rmul__(self, other):
return self * other
def __add__(self, other):
if len(self) != len(other):
raise NonConformantVectors(len(self), len(other))
return Vector(self.v + other.v)
def __radd__(self, other):
return self + other
def __sub__(self, other):
if len(self) != len(other):
raise NonConformantVectors(len(self), len(other))
return Vector(self.v - other.v)
def __eq__(self, other):
if not isinstance(other, Vector):
raise ValueError
if len(self) != len(other):
raise NonConformantVectors(len(self), len(other))
eq = numpy.isclose(self.v, other.v, EQUALITY_TOLERANCE)
if isinstance(eq, bool):
return eq
return eq.all()
def __repr__(self):
return 'Vector[{}]'.format(len(self))
def __getitem__(self, item):
return self.v[item]
def magnitude(self):
"""
Return the magnitude of the vector.
"""
return math.sqrt(sum(map(lambda x: x * x, self.v)))
def is_zero(self, tolerance=EQUALITY_TOLERANCE):
"""
Return True if the vector is a zero vector (within some tolerance).
"""
return numpy.isclose(self.magnitude(), 0, tolerance)
def unit(self):
"""
Return the unit vector of this vector. If this method is called on
a zero vector (i.e. is_zero returns True), a ValueError will be thrown.
"""
mag = self.magnitude()
if self.is_zero():
raise ValueError("cannot normalise the zero vector")
return self * (1 / mag)
def dot(self, other):
"""
Compute the dot product between this vector and the other vector.
"""
return dot(self, other)
def angle_with(self, other, in_degrees=False):
"""
Compute the angle between this vector and the other vector. The
default is to return the angle in radians. If in_degrees is True,
returns the angle in degrees.
"""
return angle(self, other, in_degrees=in_degrees)
def parallel_to(self, other):
"""Return True if the vector other is parallel to this vector."""
return parallel(self, other)
def orthogonal_to(self, other):
"""Return True if the vector other is orthogonal to this vector."""
return orthogonal(self, other)
def project_parallel(self, basis):
"""
Return the projection of this vector onto the given basis vector.
"""
unit_basis = basis.unit()
return Vector(dot(self, unit_basis) * unit_basis)
def project_orthogonal(self, basis):
"""
Compute the orthogonal projection of the vector from the given basis
vector.
"""
spar = self.project_parallel(basis)
return self - spar
# The dot (or inner) product determines the angle between two vectors.
def dot(v, w):
"""
Return the dot product of vectors v and w.
:rtype: int
"""
assert isinstance(v, Vector)
assert isinstance(v, Vector)
inner = sum([a * b for (a, b) in zip(v.v, w.v)])
# Cauchy-Schwartz inequality
assert abs(inner) <= (v.magnitude() * w.magnitude())
return inner
def angle(v, w, in_degrees=False, tolerance=EQUALITY_TOLERANCE):
"""
Return the angle between vectors v and w in radians. If in_degrees is
True, return the answer in degrees.
"""
try:
# check for floating point problems resulting in domain errors
inner = dot(v.unit(), w.unit())
if inner > 1:
inner = clamp_if_close(inner, 1.0, tolerance)
if inner < -1:
inner = clamp_if_close(inner, -1.0, tolerance)
except NonConformantVectors:
raise ValueError('Cannot determine the angle between the zero vector and another vector.')
except:
# Deliberately pass through any other exceptions.
raise
theta = math.acos(inner)
if in_degrees:
theta = r2d(theta)
return theta
def parallel(v, w, tolerance=EQUALITY_TOLERANCE):
"""Return True if vectors v and w are parallel."""
if len(v) != len(w):
raise NonConformantVectors(len(v), len(w))
if v.is_zero() or w.is_zero():
return True
# HERE LIES THE RUIN OF A DUMB
# I'd originally tried to code this by creating a list of v_i / w_i, and then reducing
# the list via comparison (e.g. are all the values in the list numpy.isclose?). I got it
# right, but... the video showed a better way. Lesson learned, think about things instead
# of blindly coding through.
theta = angle(v, w)
if numpy.isclose(theta, 0, tolerance):
return True
return numpy.isclose(theta, math.pi, tolerance)
def orthogonal(v, w, tolerance=EQUALITY_TOLERANCE):
"""
Return True if vectors v and w are orthogonal.
"""
if len(v) != len(w):
raise NonConformantVectors(len(v), len(w))
if v.is_zero() or w.is_zero():
return True
return numpy.isclose(dot(v, w), 0, tolerance)
def cross(v, w):
"""Return the cross product of the 3D vectors v and w."""
if len(v) != 3:
raise NonConformantVectors(3, len(v))
if len(w) != 3:
raise NonConformantVectors(3, len(w))
(x1, y1, z1) = v.v
(x2, y2, z2) = w.v
xc = (y1 * z2) - (y2 * z1)
# - (x1 * z2) - (x2 * z1)
yc = (x1 * z2) - (x2 * z1)
# x1 * y2 - x2 * y1
zc = (x1 * y2) - (x2 * y1)
return Vector(xc, -yc, zc)
def area_parallelogram(v, w):
"""Return the area of a parallelogram formed by Vectors v and w."""
u = cross(v, w)
return u.magnitude()
def area_triangle(v, w):
"""Return the area of a triangle formed by Vectors v and w."""
return 0.5 * area_parallelogram(v, w)
def r2d(rval):
"""
Convert the integer or floating point radian value to degrees. The
radian value must be between :math:`-2*\pi <= rval <= 2*\pi`.
>>> r2d(math.pi)
180.0
>>> r2d(3 * math.pi / 4)
135.0
>>> r2d(3 * math.pi)
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
"""
assert abs(rval) <= (2 * math.pi)
return rval * 180 / math.pi
def clamp_if_close(value, clamped, tolerance=EQUALITY_TOLERANCE):
"""
Return clamped if value is close to clamped (within tolerance), or return
value otherwise.
>>> clamp_if_close(0.99, 1.0, tolerance=0.1)
1.0
>>> clamp_if_close(0.99, 1.0, tolerance=0.001)
0.99
"""
if numpy.isclose(value, clamped, tolerance):
return clamped
return value

132
kmath/rat.py Normal file
View File

@ -0,0 +1,132 @@
import math
class Rational:
"""
Rational implements a rational number type.
"""
def __init__(self, n, d):
self.n = int(n)
self.d = int(d)
# hack to avoid DBZ
if d == 0:
self.n = 1
self.d = 0
def __float__(self) -> float:
return self.n / self.d
def __int__(self) -> int:
"""
The integer version of a Rational is the rounded float representation
of the number.
"""
return round(self.n / self.d)
def __add__(self, other):
lcm = int(math.lcm(self.d, other.d))
ks = lcm / self.d
ko = lcm / other.d
return Rational((self.n * ks) + (other.n * ko), lcm)
def __mul__(self, other):
return Rational(self.n * other.n, self.d * other.d)
def __truediv__(self, other):
return Rational(self.n * other.d, self.d * other.n)
def lcm(self, other):
return int(math.lcm(self.d, other.d))
def normalize(self, other):
lcm = self.lcm(other)
k = lcm / self.d
return self.scale(k)
def __normpair__(self, other):
rs = self.reduce()
ro = other.reduce()
ns = rs.normalize(ro)
no = ro.normalize(rs)
return ns, no
def reduce(self):
gcd = math.gcd(self.n, self.d)
return Rational(self.n / gcd, self.d / gcd)
def __eq__(self, other):
rs = self.reduce()
ro = other.reduce()
return (rs.n == ro.n) and (rs.d == ro.d)
def __gt__(self, other):
ns, no = self.__normpair__(other)
return ns.n > no.n
def __ge__(self, other):
ns, no = self.__normpair__(other)
return ns.n >= no.n
def __lt__(self, other):
ns, no = self.__normpair__(other)
return ns.n < no.n
def __le__(self, other):
ns, no = self.__normpair__(other)
return ns.n <= no.n
def scale(self, k):
return Rational(self.n * k, self.d * k)
def __sub__(self, other):
lcm = self.lcm(other)
ks = lcm / self.d
ko = lcm / other.d
return Rational((self.n * ks) - (other.n * ko), lcm)
def __pow__(self, power):
return Rational(self.n ** power, self.d ** power)
def __repr__(self):
return f'Rational({self.n}, {self.d})'
def Int(n) -> Rational:
"""
Int produces an integer version of a Rational (e.g. denominator is 1).
"""
return Rational(n, 1)
def ipow(n, m):
x = n
while m > 1:
x = x * n
return x
def expansion(v):
top = Int(1)
bottom = Int(2) + v
return top / bottom
def recurse(steps, start=None):
x = start
if x is None:
x = Rational(1, 2)
print(f'{x}')
while steps > 0:
x = expansion(x)
print(f'{x}')
steps -= 1
return x + Int(1)

2
requirements.txt Normal file
View File

@ -0,0 +1,2 @@
numpy
pytest

23
setup.py Normal file
View File

@ -0,0 +1,23 @@
# -*- coding: utf-8 -*-
# Learn more: https://github.com/kennethreitz/setup.py
from setuptools import setup, find_packages
with open('README.rst') as f:
readme = f.read()
with open('LICENSE') as f:
license_file = f.read()
setup(
name='kmath',
version='0.1.0',
description='Maths code for experimentation',
long_description=readme,
author='Kyle Isom',
author_email='kyle@imap.cc',
url='https://git.wntrmute.dev/kyle/pymath',
license=license_file,
packages=find_packages(exclude=('tests', 'docs'))
)

0
tests/__init__.py Normal file
View File

167
tests/linea.py Normal file
View File

@ -0,0 +1,167 @@
import numpy
import pytest
import kmath.linea as vec
def fequal(a, b):
eq = numpy.isclose(a, b, vec.EQUALITY_TOLERANCE)
if isinstance(eq, numpy.ndarray):
return eq.all()
return eq
def test_equality():
v1 = vec.Vector(1, 2, 3)
v2 = vec.Vector(1, 2)
v3 = vec.Vector(1, 2, 3)
v4 = vec.Vector(3, 4, 5)
assert (v1 == v3)
with pytest.raises(vec.NonConformantVectors):
assert v1 != v2
assert (v1 != v4)
# Video 4
def test_basic_operations():
v1 = vec.Vector(8.218, -9.341)
v2 = vec.Vector(-1.129, 2.111)
v3 = vec.Vector(7.089, -7.230)
assert (v1 + v2 == v3)
v1 = vec.Vector(7.119, 8.215)
v2 = vec.Vector(-8.223, 0.878)
v3 = vec.Vector(15.3420, 7.3370)
assert (v1 - v2 == v3)
v1 = vec.Vector(1.671, -1.012, -0.318)
k = 7.41
v2 = vec.Vector(12.3821, -7.4989, -2.3564)
assert (k * v1 == v2)
assert (v1 * k == v2)
# Video 6
def test_magnitude():
v1 = vec.Vector(-0.221, 7.437)
assert (fequal(v1.magnitude(), 7.4403))
v2 = vec.Vector(5.581, -2.136)
unit = vec.Vector(0.933935214087, -0.357442325262)
assert (v2.unit() == unit)
v3 = vec.Vector(8.813, -1.331, -6.247)
assert (fequal(v3.magnitude(), 10.8842))
v4 = vec.Vector(1.996, 3.108, -4.554)
unit = vec.Vector(.340401295943, 0.530043701298, -0.776647044953)
assert (v4.unit() == unit)
v5 = vec.Vector(0, 0, 0)
with pytest.raises(ValueError):
v5.unit()
# Video 8
def test_dot_product():
v1 = vec.Vector(1, 2, -1)
v2 = vec.Vector(3, 1, 0)
assert (fequal(vec.dot(v1, v2), 5))
assert (fequal(v1.dot(v2), 5))
v1 = vec.Vector(7.887, 4.138)
v2 = vec.Vector(-8.802, 6.776)
assert (fequal(vec.dot(v1, v2), -41.382))
v3 = vec.Vector(-5.955, -4.904, -1.874)
v4 = vec.Vector(-4.496, -8.755, 7.103)
assert (fequal(vec.dot(v3, v4), 56.3971))
v5 = vec.Vector(3.183, -7.627)
v6 = vec.Vector(-2.668, 5.319)
assert (fequal(vec.angle(v5, v6), 3.072))
v7 = vec.Vector(7.35, 0.221, 5.188)
v8 = vec.Vector(2.751, 8.259, 3.985)
theta = vec.angle(v7, v8, in_degrees=True)
assert (fequal(theta, 60.2758))
v9 = vec.Vector(0, 0)
with pytest.raises(ValueError):
v1.angle_with(v9)
# Video 10.
def test_parallel_orthogonal():
v1 = vec.Vector(-7.579, -7.88)
v2 = vec.Vector(22.737, 23.64)
assert (v1.parallel_to(v2))
assert (not v1.orthogonal_to(v2))
v3 = vec.Vector(-2.029, 9.97, 4.172)
v4 = vec.Vector(-9.231, -6.639, -7.245)
assert (not v3.parallel_to(v4))
assert (not v3.orthogonal_to(v4))
v5 = vec.Vector(-2.328, -7.284, -1.214)
v6 = vec.Vector(-1.821, 1.072, -2.94)
assert (not v5.parallel_to(v6))
assert (v5.orthogonal_to(v6))
v7 = vec.Vector(2.118, 4.827)
v8 = vec.Vector(0, 0)
assert v7.parallel_to(v8)
assert v7.orthogonal_to(v8)
# Video 12
def test_projection():
# sanity check
v1 = vec.Vector(1, 3)
v2 = vec.Vector(3, 3)
v3 = v1.project_parallel(v2)
v4 = v1.project_orthogonal(v2)
assert v3 + v4 == v1
v1 = vec.Vector(3.039, 1.879)
v2 = vec.Vector(0.825, 2.036)
v3 = vec.Vector(1.0826, 2.6717)
assert v1.project_parallel(v2) == v3
v4 = vec.Vector(-9.88, -3.264, -8.159)
v5 = vec.Vector(-2.155, -9.353, -9.473)
v6 = vec.Vector(-8.350, 3.376, -1.434)
assert v4.project_orthogonal(v5) == v6
v7 = vec.Vector(3.009, -6.172, 3.692, -2.510)
v8 = vec.Vector(6.404, -9.144, 2.759, 8.718)
v9 = vec.Vector(1.969, -2.811, 0.848, 2.680)
assert v7.project_parallel(v8) == v9
v10 = vec.Vector(1.040, -3.361, 2.844, -5.190)
assert v7.project_orthogonal(v8) == v10
assert (v9 + v10) == v7
def test_cross_product():
# sanity check
v1 = vec.Vector(5, 3, -2)
v2 = vec.Vector(-1, 0, 3)
v3 = vec.Vector(9, -13, 3)
assert vec.cross(v1, v2) == v3
v1 = vec.Vector(8.462, 7.893, -8.187)
v2 = vec.Vector(6.984, -5.975, 4.778)
v3 = vec.Vector(-11.205, -97.609, -105.685)
assert vec.cross(v1, v2) == v3
v4 = vec.Vector(-8.987, -9.838, 5.031)
v5 = vec.Vector(-4.268, -1.861, -8.866)
area = 142.122
assert fequal(vec.area_parallelogram(v4, v5), area)
v6 = vec.Vector(1.500, 9.547, 3.691)
v7 = vec.Vector(-6.007, 0.124, 5.772)
area = 42.565
assert fequal(vec.area_triangle(v6, v7), area)