from math import radians, degrees, cos, sin, acos
import math
from operator import add, sub, mul, neg
from collections import namedtuple

XYZ = namedtuple("xyz_vector", "x y z")
IVM = namedtuple("ivm_vector", "a b c d")

root2   = 2.0**0.5

        
class Qvector:
    """Quadray vector"""

    def __init__(self, arg):
        """Initialize a vector at an (x,y,z)"""
        self.coords = self.norm(arg)

    def __repr__(self):
        return repr(self.coords)

    def norm(self, arg):
        """Normalize such that 4-tuple all non-negative members."""
        return IVM(*tuple(map(sub, arg, [min(arg)] * 4))) 
    
    def norm0(self):
        """Normalize such that sum of 4-tuple members = 0"""
        q = self.coords
        return IVM(*tuple(map(sub, q, [sum(q)/4.0] * 4))) 

    @property
    def a(self):
        return self.coords.a

    @property
    def b(self):
        return self.coords.b

    @property
    def c(self):
        return self.coords.c

    @property
    def d(self):
        return self.coords.d
    
    def __eq__(self, other):
        return self.coords == other.coords
        
    def __lt__(self, other):
        return self.coords < other.coords

    def __gt__(self, other):
        return self.coords > other.coords
    
    def __hash__(self):
        return hash(self.coords)
    
    def __mul__(self, scalar):
        """Return vector (self) * scalar."""
        newcoords = [scalar * dim for dim in self.coords]
        return Qvector(newcoords)

    __rmul__ = __mul__ # allow scalar * vector

    def __truediv__(self,scalar):
        """Return vector (self) * 1/scalar"""        
        return self.__mul__(1.0/scalar)
    
    def __add__(self,v1):
        """Add a vector to this vector, return a vector""" 
        newcoords = tuple(map(add, v1.coords, self.coords))
        return Qvector(newcoords)
        
    def __sub__(self,v1):
        """Subtract vector from this vector, return a vector"""
        return self.__add__(-v1)
    
    def __neg__(self):      
        """Return a vector, the negative of this one."""
        return Qvector(tuple(map(neg, self.coords)))
                  
    def dot(self,v1):
        """Return the dot product of self with another vector.
        return a scalar
        
        >>> s1 = a.dot(b)/(a.length() * b.length())
        >>> degrees(acos(s1))
        109.47122063449069
        """
        return 0.5 * sum(map(mul, self.norm0(), v1.norm0()))

    def length(self):
        """Return this vector's length"""
        return self.dot(self) ** 0.5
        
    def cross(self,v1):
        """Return the cross product of self with another vector.
        return a Qvector"""
        A = Qvector((1,0,0,0))
        B = Qvector((0,1,0,0))
        C = Qvector((0,0,1,0))
        D = Qvector((0,0,0,1))
        a1,b1,c1,d1 = v1.coords
        a2,b2,c2,d2 = self.coords
        k= (2.0**0.5)/4.0
        sum =   (A*c1*d2 - A*d1*c2 - A*b1*d2 + A*b1*c2
               + A*b2*d1 - A*b2*c1 - B*c1*d2 + B*d1*c2 
               + b1*C*d2 - b1*D*c2 - b2*C*d1 + b2*D*c1 
               + a1*B*d2 - a1*B*c2 - a1*C*d2 + a1*D*c2
               + a1*b2*C - a1*b2*D - a2*B*d1 + a2*B*c1 
               + a2*C*d1 - a2*D*c1 - a2*b1*C + a2*b1*D)
        return k*sum

    def angle(self, v1):
        return self.xyz().angle(v1.xyz())
        
    def xyz(self):
        a,b,c,d     =  self.coords
        k           =  0.5/root2
        xyz         = (k * (a - b - c + d),
                       k * (a - b + c - d),
                       k * (a + b - c - d))
        return Vector(xyz)

"""
Explanation: STEM for Philosophers series<br/>Oregon Curriculum Network
Quadray Coordinates: Getting Started
What are quadray coordinates and how are they used in philosophy?  For more background on this question, read my Thinking Outside the Box: Language Games in Mathematics, on Medium
Lets start out with a stripped down XYZ Vector class that works pretty much as expected, in allowing for vector addition and subtraction, multiplication by a scalar.
End of explanation
"""


class Vector:

    def __init__(self, arg):
        """Initialize a vector at an (x,y,z)"""
        self.xyz = XYZ(*map(float,arg))

    def __repr__(self):
        return repr(self.xyz)
    
    @property
    def x(self):
        return self.xyz.x

    @property
    def y(self):
        return self.xyz.y

    @property
    def z(self):
        return self.xyz.z
    
    def __mul__(self, scalar):
        """Return vector (self) * scalar."""
        newcoords = [scalar * dim for dim in self.xyz]
        return type(self)(newcoords)

    __rmul__ = __mul__ # allow scalar * vector

    def __truediv__(self,scalar):
        """Return vector (self) * 1/scalar"""        
        return self.__mul__(1.0/scalar)
    
    def __add__(self,v1):
        """Add a vector to this vector, return a vector""" 
        newcoords = map(add, v1.xyz, self.xyz)
        return type(self)(newcoords)
        
    def __sub__(self,v1):
        """Subtract vector from this vector, return a vector"""
        return self.__add__(-v1)
    
    def __neg__(self):      
        """Return a vector, the negative of this one."""
        return type(self)(tuple(map(neg, self.xyz)))

    def unit(self):
        return self.__mul__(1.0/self.length())

    def dot(self,v1):
        """Return scalar dot product of this with another vector."""
        return sum(map(mul , v1.xyz, self.xyz))

    def cross(self,v1):
        """Return the vector cross product of this with another vector"""
        newcoords = (self.y * v1.z - self.z * v1.y, 
                     self.z * v1.x - self.x * v1.z,
                     self.x * v1.y - self.y * v1.x )
        return type(self)(newcoords)
    
    def length(self):
        """Return this vector's length"""
        return self.dot(self) ** 0.5

    def quadray(self):
        """return (a, b, c, d) quadray based on current (x, y, z)"""
        x, y, z = self.xyz
        k = 2/root2
        a = k * ((x >= 0)* ( x) + (y >= 0) * ( y) + (z >= 0) * ( z))
        b = k * ((x <  0)* (-x) + (y <  0) * (-y) + (z >= 0) * ( z))
        c = k * ((x <  0)* (-x) + (y >= 0) * ( y) + (z <  0) * (-z))
        d = k * ((x >= 0)* ( x) + (y <  0) * (-y) + (z <  0) * (-z))
        return Qvector((a, b, c, d))

"""
Explanation: Converting to xyz will not work yet, as the Vector class is not yet defined.  That's what's coming.
End of explanation
"""


octant0 = Vector((root2/4, root2/4, root2/4))
print(octant0.xyz)
q0 = octant0.quadray()
print(q0)

q0.length()

"""
Explanation: At the end is a method for outputting in quadray coordinates.  
Some design decisions were taken, conventions followed, in how the XYZ and IVM systems were overlaid. 
Lets not worry about that for now and just imagine a cube with edges sqrt(2)/2, one corner in each octant.  The face diagonals will have length 1 in this case.
For example, in the all positive octant (+ + +) we would have a point at (sqrt(2)/4, sqrt(2)/4, sqrt(2)/4).
End of explanation
"""


octant1 = Vector((-root2/4, root2/4, root2/4))  # neighboring octant
diff = octant0 - octant1
print("diff.quadray()", diff.quadray())
print("Length between adjacent corners: ", diff.length())

"""
Explanation: This might seem strange already. What appears to be a unit vector, has some irrational length. The cube below is flipped over somehow, but gives the idea.  Think of (1,0,0,0) as being in the all positive octant (+, +, +) of XYZ.
Another issue with this cube, in the context of the surrounding computations, is all edges are doubled, with cube face diagonals set at 2.  Divide every linear dimension through by 2 to get the computations here.
The body diagonal of a cube with unit face diagonas is $\sqrt{6}/2$ meaning each quadray, which goes from the cube center to a corner, has a length of $\sqrt(6)/4$. That's the lenght of q0 shown above.

<div style="text-align: center;">
Fig 1:  Divide Through by 2 for Dimensions
</div>
End of explanation
"""


diff.quadray().length()

a = Qvector((1,0,0,0))
b = Qvector((0,1,0,0))
(a-b).length()

"""
Explanation: This confirms the cube has the expected edge length.  This is not the length between two quadray tips, which is $1.0$, but between two adjacent corners of our cube, so $\sqrt{2}/2$.
End of explanation
"""


# add up three quadrays and negate their sum, to get the other Qray
a = Qvector((1,0,0,0))
c = Qvector((0,0,1,0))
d = Qvector((0,0,0,1))
v_sum = -(a + c + d)
print("Canonical representation:", v_sum)
print("Alternative expression:  ", v_sum.norm0())
print("v_sum length:            ", v_sum.length())

"""
Explanation: Half the cube's vertices will align with the four spokes of the caltrop (in blue). These correspond the the vertexes of an embedded tetrahedron of edges 2 (in red), or edges 1 if measuring in D units.

Note that the Qvector class comes with two ways to express a Qray in canonical lowest terms.  One way preserves the non-negative coordinate address for every point.  The other way assures that the 4-tuple coordinates sum to zero.  I'm using the former for all representations (repr) whereas the latter gets used in various internal computations.
End of explanation
"""


