CODEWITHSUNDEEP
pythonmathvector
Goddard
Goddard

Reputation: 3059

acos getting error math domain error

Trying to figure out why I am getting an error. My numbers are between -1 and 1, but still errors.

ValueError: math domain error

Any ideas?

Thanks

from math import sqrt, acos, pi
from decimal import Decimal, getcontext

getcontext().prec = 30


class Vector(object):
    CANNOT_NORMALIZE_ZERO_VECTOR_MSG = 'Cannot normalize the zero vector'

    def __init__(self, coordinates):
        try:
            if not coordinates:
                raise ValueError
            self.coordinates = tuple([Decimal(x) for x in coordinates])
            self.dimension = len(self.coordinates)

        except ValueError:
            raise ValueError('The coordinates must be nonempty')

        except TypeError:
            raise TypeError('The coordinates must be an iterable')

    def __str__(self):
        return 'Vector: {}'.format(self.coordinates)

    def __eq__(self, v):
        return self.coordinates == v.coordinates

    def magnitude(self):
        coordinates_squared = [x ** 2 for x in self.coordinates]
        return sqrt(sum(coordinates_squared))

    def normalized(self):
        try:
            magnitude = self.magnitude()
            return self.times_scalar(Decimal(1.0 / magnitude))

        except ZeroDivisionError:
            raise Exception('Cannot normalize the zero vector')

    def plus(self, v):
        new_coordinates = [x + y for x, y in zip(self.coordinates, v.coordinates)]
        return Vector(new_coordinates)

    def minus(self, v):
        new_coordinates = [x - y for x, y in zip(self.coordinates, v.coordinates)]
        return Vector(new_coordinates)

    def times_scalar(self, c):
        new_coordinates = [Decimal(c) * x for x in self.coordinates]
        return Vector(new_coordinates)

    def dot(self, v):
        return sum([x * y for x, y in zip(self.coordinates, v.coordinates)])

    def angle_with(self, v, in_degrees=False):
        try:
            u1 = self.normalized()
            u2 = v.normalized()
            angle_in_radians = acos(u1.dot(u2))

            if in_degrees:
                degrees_per_radian = 180. / pi
                return angle_in_radians * degrees_per_radian
            else:
                return angle_in_radians

        except Exception as e:
            if str(e) == self.CANNOT_NORMALIZE_ZERO_VECTOR_MSG:
                raise Exception('Cannot comput an angle with a zero vector')
            else:
                raise e

    def is_orthogonal_to(self, v, tolerance=1e-10):
        return abs(self.dot(v)) < tolerance

    def is_parallel_to(self, v):
        return self.is_zero() or v.is_zero() or self.angle_with(v) == 0 or self.angle_with(v) == pi

    def is_zero(self, tolerance=1e-10):
        return self.magnitude() < tolerance


print('first pair...')
v = Vector(['-7.579', '-7.88'])
w = Vector(['22.737', '23.64'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

print('second pair...')
v = Vector(['-2.029', '9.97', '4.172'])
w = Vector(['-9.231', '-6.639', '-7.245'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

print('third pair...')
v = Vector(['-2.328', '-7.284', '-1.214'])
w = Vector(['-1.821', '1.072', '-2.94'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

print('fourth pair...')
v = Vector(['2.118', '4.827'])
w = Vector(['0', '0'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

Upvotes: 0

Views: 2019

Answers (1)

OrderAndChaos
OrderAndChaos

Reputation: 3860

Could it be that u1.dot(u2) equals -1.00000000000000018058942747512

print(u2)
print(u1.dot(u2))
angle_in_radians = acos(u1.dot(u2))

This is around line 60

Update, with further tests:

getcontext().prec = 16

......

def dot(self, v):
    print(self.coordinates, v.coordinates)
    print("asf")
    result = 0
    for x, y in zip(self.coordinates, v.coordinates):
        print("=================")
        print("x: ", x)
        print("y: ", y)
        print("x*y: ", x*y)
        result += (x*y)
        print("=================")
    print("Result: ", result)
    print(sum([x * y for x, y in zip(self.coordinates, v.coordinates)]))
    return sum([x * y for x, y in zip(self.coordinates, v.coordinates)])

Results in:

=================
x:  -0.6932074151971374
y:  0.6932074151971375
x*y:  -0.4805365204842965
=================
=================
x:  -0.7207381490636552
y:  0.7207381490636553
x*y:  -0.5194634795157037
=================
Result:  -1.000000000000000
-1.000000000000000

But with:

getcontext().prec = 30

The decimal begins to drift.

=================
x:  -0.693207415197137377521618972764
y:  0.693207415197137482701372768190
x*y:  -0.480536520484296481693529594664
=================
=================
x:  -0.720738149063655170190045851086
y:  0.720738149063655279547013776664
x*y:  -0.519463479515703698895897880460
=================
Result:  -1.00000000000000018058942747512

Which leaves the result less than -1 breaking the acos() function.

After finding the floats were out, I looked through your code I noticed a couple of functions that return floats. The culprit is the sqrt() function which doesn't have a high enough accuracy.

def magnitude(self):
    coordinates_squared = [x ** 2 for x in self.coordinates]
    return Decimal(sum(coordinates_squared)).sqrt()

def normalized(self):
    try:
        magnitude = self.magnitude()
        return self.times_scalar(Decimal(1.0) / magnitude)

Using the Decimal(x).sqrt() function will fix your issue. You'll then need to update the normalized() function a bit too.

Upvotes: 2

Related Questions