0

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))
Goddard
  • 2,863
  • 31
  • 37
  • 1
    Questions seeking debugging help (**"why isn't this code working?"**) must include the desired behavior, *a specific problem or error* and *the shortest code necessary* to reproduce it **in the question itself**. Questions without **a clear problem statement** are not useful to other readers. See: [How to create a Minimal, Complete, and Verifiable Example](http://stackoverflow.com/help/mcve). Please don't dump a wall of code and expect people to troubleshoot the whole thing. Only post the **shortest code necessary** to replicate the problem. – MattDMo Oct 12 '16 at 23:45
  • 1
    Please see the bottom of my answer, I now have a solution. – OrderAndChaos Oct 13 '16 at 18:05

1 Answers1

2

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.

OrderAndChaos
  • 3,547
  • 2
  • 30
  • 57
  • I don't ever see the resulting value during debugging. How can I fix this? – Goddard Oct 13 '16 at 02:38
  • Litter the code with print statements, or set up a debugging environment in an IDE like PyCharm (which has a free version) and use breakpoints. https://www.jetbrains.com/help/pycharm/2016.1/debugging.html – OrderAndChaos Oct 13 '16 at 08:34
  • I posted it here to get help with the answer. I already knew where the problem was occurring. I am already using an IDE and know how to debug. I don't know where the problem is REALLY happening which is why I asked for help. – Goddard Oct 13 '16 at 16:28
  • Ha, sorry, fair enough. What are you expecting to get back from dot()? Do you think it be a rounding float error? – OrderAndChaos Oct 13 '16 at 16:33
  • How accurate do you need the result to be? – OrderAndChaos Oct 13 '16 at 16:35
  • So this only seems to happen in the first test is that right? – OrderAndChaos Oct 13 '16 at 16:37
  • I'm guessing you need it to be precise but the issue goes away if you reduce the precision to 16, so I'd imagine the problem as to the decimal place losing accuracy. You probably worked that out already. I'd suggest that you would be best off posting the question again but ask about maintaining accuracy with decimal places. – OrderAndChaos Oct 13 '16 at 16:58
  • @Goddard if this is indeed the problem, it's an easy fix: `angle_in_radians = acos(round(u1.dot(u2), 15))` – Mark Ransom Oct 13 '16 at 17:37
  • Hold on, it may be the `sqrt()` magintude. – OrderAndChaos Oct 13 '16 at 17:38
  • magnitude and normalise are returning float type objects, so It could be around there. – OrderAndChaos Oct 13 '16 at 17:39
  • I've just been reading up about Decimal and it should be accurate from what I can see, but it does look like a floating point error is occurring somewhere as you lose accuracy around 16 decimal places which is about how accurate floats are. – OrderAndChaos Oct 13 '16 at 17:42
  • Maybe this could help: http://stackoverflow.com/questions/10725522/arbitrary-precision-of-square-roots – OrderAndChaos Oct 13 '16 at 17:43
  • Let me know if you do get to the bottom of this, the round you suggest looks as though it would do the trick. I suspect the culprit is the sqrt() function though. – OrderAndChaos Oct 13 '16 at 17:57
  • So what does it actually do? – OrderAndChaos Oct 13 '16 at 18:05
  • Excellent, thanks. It is just vector math. I was trying to learn how to multiple vectors etc.. – Goddard Oct 14 '16 at 19:37