This is Me Using LJ as a Code Stash

[astarcambiata]

#! /usr/bin/env python # coding: UTF-8 #HCoxeter.py import cmath, sys from math import * from spraque3D import * from poincomp import * def maketuple(strin): """Maketup takes a (somewhat predictably and strictly formatted) matrix as a string input, and strips the string to create a nested tuple matrix.""" t = () #strin = strin + ' ,' # commas are used to partition strin into a tuple strin = strin.replace('[','') strin = strin.replace(']','*%') strin = strin.replace(',','*') strin = strin.replace(';', '*%') while '%' in strin: tadd = strin[:strin.find('%')] t = t + (tadd,) strin = strin[strin.find('%') + 1:] for i in range(len(t)): a = () ti = t[i] while '*' in ti: ad = float(ti[:ti.find('*')]) a = a + (ad,) ti = ti[ti.find('*') + 1:] t = t[:i] + (a,) + t[i + 1:] return t G = () while len(G) != 4: G = raw_input('Input Gram matrix (must be 4x4): ') G = maketuple(G) if len(G) != 4: print "Invalid Gram matrix dimensions." for i in range(len(G)): if len(G[i]) != 4: print "Invalid Gram matrix dimensions." Q = GramtoVertNorm(G) print Q GenTet = HypeTet(Q[0]) print GenTet #! /usr/bin/env python # coding: UTF-8 #poincomp.py import cmath, sys from math import * from spraque3D import * def vspherinv3(v,r): """Inversion of a vector v in a 2-sphere of radius r in R^3""" l = veclengthSq3(v) if l == 0.0: return '∞' w = r*r/l return vecscale3(v,w) def Qhcntr(p1,p2,r): """Quaternion utilizing method to find center and radius of an h-line through p1 & p2 in a Poincaré ball of radius r""" pl1 = veclength3(p1) pl2 = veclength3(p2) if pl1 == r and pl2 == r: if vecdot3(p1,p2) == 0.0: return vecadd3(p1,p2) else: a1 = p1 a2 = p2 elif pl1 != r: a1 = vecscale3(vecadd3(p1,vspherinv3(p1,r)), 0.500) if pl2 == r: a2 = p2 else: a2 = vecscale3(vecadd3(p2,vspherinv3(p2,r)), 0.500) else: a1 = p1 a2 = vecscale3(vecadd3(p2,vspherinv3(p2,r)), 0.500) n = vecross(a1,a2) if veclengthSq3(n) == 0.0: return zerovec3() b1 = vecnormalize3(spatrot(0.500*pi, n, a1)) b2 = vecnormalize3(spatrot(-0.500*pi, n, a2)) t = (vecdot3(a1,b2) + vecdot3(a2,b1))/veclengthSq3(vecsub3(b2,b1)) c = vecscale3(vecadd3(vecadd3(a1,vecscale3(b1,t)), vecadd3(a2,vecscale3(b2,t))), 0.500) return (t,) + c def Hplncntr(p1,p2,p3,R): """Linear algebra method to find h-planes in the Poincaré model of H^3""" a1 = vecsub3(p1,p2) a2 = vecsub3(p1,p3) n = vecross(a1,a2) if n == zerovec3(): return zeroQ() elif vecdot3(p1,n) == 0.0: return (0.0,) + vecnormalize3(n) else: M = (p1, p2, p3) c = R*R v = vecscale3((veclengthSq3(p1) + c, veclengthSq3(p2) + c, veclengthSq3(p3) + c), 0.500) h = Matvecprd3(M3inverse(M), v) r = sqrt(veclengthSq3(h) - c) return (r,) + h def H3HtoH3C(p): """Sterographic projection from the hyperboloid (Minkowski space) model of H^3 to the Poincaré ball.""" c = 1.000/(p[0]+1.000) return (p[1]*c, p[2]*c, p[3]*c) def H3CtoH3H(v): """Inverse steroegraphic projection from the Poincaré ball to the hyperboloid model.""" s = veclengthSq(v) c = 2.000/(1.000-s) return (0.500*(1.000+s)*c, c*v[0], c*v[1], c*v[2]) def GramtoVertNorm(G): """From a Gram Matrix, returns the Vertices and Normals of a Coxeter tetrahedron in the hyperboloid model. Uses Jacobi rotation to do eigen decomposition of the real symmetric Gram matrix.""" U, L = Jacobirot(G,10) U, L = EigenSort(U,L) D = DetDiag(L) J = MinkMeTensor(4) V = MatrixProd(SqrtDiagM(L),Mtranspose(U)) W = Mscale(MatrixProd(J,MatrixProd(InvSqrtDiagM(L),U)),sqrt(abs(D))) return (V,W) def Hreflect(p,u): return Vecsub(p,Vecscale(u,2.0*Minkerprod4(p,u)/Indefmetric4(u))) def TetSphere(p1,p2,p3,p4): """From the vertices of a tetrahedron in the Poincaré model, computes the unique euclidean sphere coincident with all four vertices.""" r1 = veclengthSq3(p1) r2 = veclengthSq3(p2) r3 = veclengthSq3(p3) r4 = veclengthSq3(p4) M1 = M4Determinant(((p1[0],p1[1],p1[2],1.0), (p2[0],p2[1],p2[2],1.0), (p3[0],p3[1],p3[2],1.0), (p4[0],p4[1],p4[2],1.0))) M2 = M4Determinant(((r1,p1[1],p1[2],1.0), (r2,p2[1],p2[2],1.0), (r3,p3[1],p3[2],1.0), (r4,p4[1],p4[2],1.0))) M3 = M4Determinant(((r1,p1[0],p1[2],1.0), (r2,p2[0],p2[2],1.0), (r3,p3[0],p3[2],1.0), (r4,p4[0],p4[2],1.0))) M4 = M4Determinant(((r1,p1[0],p1[1],1.0), (r2,p2[0],p2[1],1.0), (r3,p3[0],p3[1],1.0), (r4,p4[0],p4[1],1.0))) M5 = M4Determinant(((r1,p1[0],p1[1],p1[2]), (r2,p2[0],p2[1],p2[2]), (r3,p3[0],p3[1],p3[2]), (r4,p4[0],p4[1],p4[2]))) c = vecscale3((M2,-M3,M4),1.0/(2.00*M1)) r = sqrt(veclengthSq3(c) - (M5/M1)) return (r,) + c def HypeTet(V): """Takes the vertices of the Tetrahedron (in Hyperboloid model) and determines the spheres necessary for the ray-tracer to use CSG to illustrate it. """ p1,p2,p3,p4 = Mtranspose(V) p1,p2,p3,p4 = H3HtoH3C(p1),H3HtoH3C(p2),H3HtoH3C(p3),H3HtoH3C(p4) return (TetSphere(p1,p2,p3,p4), Hplncntr(p2,p3,p4,1.0), Hplncntr(p1,p3,p4,1.0), Hplncntr(p1,p2,p4,1.0), Hplncntr(p1,p2,p3,1.0)) import cmath, sys from math import * def Zerovec(n): v = () for i in range(n): v = v + (0.0,) return v def zerovec3(): return (0.0,0.0,0.0) def Veccopy(v): return v def Vecinverse(v): u = () for i in range(len(v)): u = u + (v[i],) return u def vecinverse3(v): return (-v[0],-v[1],-v[2]) def Vecadd(u,v): a = () for i in range(len(u)): a = a + (u[i] + v[i],) return a def vecadd3(v1,v2): return (v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]) def Vecsub(u,v): a = () for i in range(len(u)): a = a + (u[i] - v[i],) return a def vecsub3(v1,v2): return (v1[0]-v2[0],v1[1]-v2[1],v1[2]-v2[2]) def Vecscale(v,s): u = () for i in range(len(v)): u = u + (v[i]*s,) return u def vecscale3(v,s): return (v[0]*s,v[1]*s,v[2]*s) def VeclengthSq(v): l = 0.0 for i in range(len(v)): l = l + v[i]*v[i] return l def veclengthSq3(v): return v[0]*v[0]+v[1]*v[1]+v[2]*v[2] def Veclength(v): l = 0.0 for i in range(len(v0)): l = l + v[i]*v[i] return sqrt(l) def veclength3(v): return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]) def Vecnormalize(v): l = Veclength(v) if l == 0: return Zerovec(n) s = 1.0/l return Vecscale(v,s) def vecnormalize3(v): l = veclength3(v) if l == 0: return (0.0,0.0,0.0) return (v[0]/l,v[1]/l,v[2]/l) def Vecdot(v1,v2): a = 0.0 for i in range(len(v1)): a = a + v1[i]*v2[i] return a def vecdot3(v1,v2): return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2] def vecross(v1,v2): return (v1[1]*v2[2]-v1[2]*v2[1],v1[2]*v2[0]-v1[0]*v2[2],v1[0]*v2[1]-v1[1]*v2[0]) def perpendicular3(v): if v[1]==0 and v[2]==0: return vecross(v,vecadd3(v,(0,1,0))) return vecross(v,vecadd3(v,(1,0,0))) def IdentityMatrix(n): M=() for i in range(n): a = () for j in range(n): a = (0.0,) + a a = a[:i] + (1.0,) + a[i+1:] M = M + (a,) return M def Idmatrix3(): return ((1.0,0.0,0.0), (0.0,1.0,0.0), (0.0,0.0,1.0)) def Zeromatrix(n): M=() for i in range(n): a = () for j in range(n): a = (0.0,) + a M = M + (a,) return M def zeromatrix3(): return (zerovec3(), zerovec3(), zerovec3()) def Diagonals(M): v = () # Extracts a vector of the main diagonal of a matrix. for i in range(len(M)): d = M[i] v = v + (d[i],) return v def DetDiag(M): D = 1.0 # M must be an upper triangular, lower triangular, or Diagonal Matrix! Returns the determinant. for i in range(len(M)): v = M[i] D = D*v[i] return D def DiagMat(M): D=() # Makes a diagonal Matrix with the same main diagonal as M. n = len(M) for i in range(n): a = () d = M[i] for j in range(n): a = (0.0,) + a a = a[:i] + (d[i],) + a[i+1:] D = D + (a,) return D def SqrtDiagM(M): D=() # Makes a diagonal matrix whose non-zero entries are square roots of the absolute values of entries on M's main diagonal. n = len(M) for i in range(n): a = () d = M[i] for j in range(n): a = (0.0,) + a a = a[:i] + (sqrt(abs(d[i])),) + a[i+1:] D = D + (a,) return D def InvDiagM(M): D=() # Computes the inverse of the diagonal matrix whose nonzero entries are taken from M's main diagonal. n = len(M) for i in range(n): a = () d = M[i] for j in range(n): a = (0.0,) + a a = a[:i] + (1.0/d[i],) + a[i+1:] D = D + (a,) return D def InvSqrtDiagM(M): D=() # Makes a diagonal matrix whose non-zero entries are square roots of the absolute values of the reciprocals of entries on M's main diagonal. n = len(M) for i in range(n): a = () d = M[i] for j in range(n): a = (0.0,) + a a = a[:i] + (1.0/sqrt(abs(d[i])),) + a[i+1:] D = D + (a,) return D def Dyad(u,v): M=() # Makes a tensor from the vectors u and v. for i in range(len(u)): a = () for j in range(len(v)): a = a[:j] + (u[i]*v[j],) M = M + (a,) return M def Matvecprod(M,v): u = () # Computes the matrix-vector product of M and v. for i in range(len(M)): u = u[:i] + (Vecdot(M[i],v),) return u def Matvecprd3(M,v): return (vecdot3(M[0],v), vecdot3(M[1],v), vecdot3(M[2],v)) def Matadd(M1,M2): M=() for i in range(len(M1)): M = M + (Vecadd(M1[i],M2[i]),) return M def Matadd3(M1,M2): u1, u2, u3 = M1 v1, v2, v3 = M2 return (vecadd(u1,v1), vecadd(u2,v2), vecadd(u3,v3)) def Matsub(M1,M2): M=() for i in range(len(M1)): M = M + (Vecsub(M1[i],M2[i]),) return M def Matsub3(M1,M2): u1, u2, u3 = M1 v1, v2, v3 = M2 return (vecsub(u1,v1), vecsub(u2,v2), vecsub(u3,v3)) def Mtranspose(M): A = () for i in range(len(M[1])): a=() for j in range(len(M)): v = M[j] a = a + (v[i],) A = A + (a,) return A def M3transpose(M): u1, u2, u3 = M return ((u1[0], u2[0], u3[0]), (u1[1], u2[1], u3[1]), (u1[2], u2[2], u3[2])) def MatRowSwap(M,i,k): i = i - 1 #Swaps ith and kth rows... Needs i< k to work! k = k - 1 a = M[i] b = M[k] M = M[:i] + (b,) + M[i+1:k] + (a,) + M[k+1:] return M def MatColSwap(M,j,k): j = j - 1 #Swaps ith and kth columns... Needs i< k to work! k = k - 1 A = Mtranspose(M) a = A[j] b = A[k] M = Mtranspose(M[:j] + (b,) + M[j+1:k] + (a,) + M[k+1:]) return M def Mscale(M,s): A = () for i in range(len(M)): A = A + (Vecscale(M[i],s),) return A def M3scale(M,s): u1, u2, u3 = M return (vecscale3(u1,s), vecscale3(u2,s), vecscale3(u3,s)) def MatrixProd(M1,M2): A, B = M1, Mtranspose(M2) M = () for i in range(len(A)): a = A[i] m = () for j in range(len(B)): b = B[j] m = m + (Vecdot(a,b),) M = M + (m,) return M def M3product(M1,M2): u1, u2, u3 = M1 v1, v2, v3 = M3transpose(M2) return ((vecdot3(u1,v1), vecdot3(u1,v2), vecdot3(u1,v3)), (vecdot3(u2,v1), vecdot3(u2,v2), vecdot3(u2,v3)), (vecdot3(u3,v1), vecdot3(u3,v2), vecdot3(u3,v3))) def M3Determinant(M): u1, u2, u3 = M v1 = vecross(u2, u3) return vecdot3(u1, v1) def M3inverse(M): u1, u2, u3 = M v1 = vecross(u2, u3) D = vecdot3(u1, v1) if D == 0.0: print 'Matrix M = ' + repr(M) + ' is singular.' return zeromatrix3() v2 = vecross(u3, u1) v3 = vecross(u1, u2) c = 1.0/D B = M3transpose((v1,v2,v3)) return M3scale(B,c) def SubMat(M,i,j): M = Mtranspose(M[:i] + M[i + 1:]) return Mtranspose(M[:j] + M[j + 1:]) def Cofactor4(M,i,j): return M3Determinant(SubMat(M,i,j))*((-1.0)**(i+j)) def M4Determinant(M): return Vecdot(M[1], (Cofactor4(M,1,1), Cofactor4(M,1,2), Cofactor4(M,1,3), Cofactor4(M,1,4))) def JacobiMatrix(n, p, q, theta): M=() c, s = cos(theta), sin(theta) for i in range(n): v = () for j in range(n): if i == j== p or i == j == q: v = v + (c,) elif i == q and j == p: v = v + (s,) elif i == p and q == j: v = v + (-s,) elif i == j: v = v + (1.0,) else: v = v + (0.0,) M = M + (v,) return M def Jacobirot(M,sweeps=2): n = len(M) J = E = P = IdentityMatrix(n) Q = M theta = 0.0 for w in range(sweeps): for i in range(n - 1): q = Q[i] for j in range(n-i-1): k = i + j + 1 p = Q[k] if q[i] == p[k]: theta = 0.25*pi else: x = q[k] theta = 0.5*atan(2.0*x/(q[i] - p[k])) P = JacobiMatrix(n, i, k, theta) Q = MatrixProd(MatrixProd(Mtranspose(P),Q),P) E = MatrixProd(E,P) G = MatrixProd(MatrixProd(Mtranspose(E),M),E) Q = Mscale(Matadd(Q,G), 0.5) return (E,Q) def EigenSort(E,D): for i in range(len(D)): v = D[i] k = 0 if v[i] < v[k] and i > k: D = MatColSwap(MatRowSwap(D, k + 1, i + 1), k + 1, i + 1) E = MatColSwap(E, k + 1, i + 1) k = k + 1 return (E,D) def zeroQ(): return (1.0,0.0,0.0,0.0) def MinkQuadform(v): a = v[0]*v[0] for i in range(len(v)-1): a = a - v[i+1]*v[i+1] return a def MinkBiliform(v1,v2): a = v1[0]*v2[0] for i in range(len(v1)-1): a = a - v1[i+1]*v2[i+1] return a def MinkMeTensor(n): M=() d = (-1.0,) for k in range(n - 1): d = d + (1.0,) for i in range(n): a = () for j in range(n): a = (0.0,) + a a = a[:i] + (d[i],) + a[i+1:] M = M + (a,) return M def Minkerprod4(v1,v2): return v1[1]*v2[1]+v1[2]*v2[2]+v1[3]*v2[3]-v1[0]*v2[0] def Indefmetric4(v): return v[1]*v[1]+v[2]*v[2]+v[3]*v[3]-v[0]*v[0] def copyQ(q): return (q[0],q[1],q[2],q[3]) def vectoquat(v): return (0.0,v[0],v[1],v[2]) def addQ(q1,q2): return (q1[0]+q2[0],q1[1]+q2[1],q1[2]+q2[2],q1[3]+q2[3]) def subQ(q1,q2): return (q1[0]-q2[0],q1[1]-q2[1],q1[2]-q2[2],q1[3]-q2[3]) def scaleQ(q,s): return (q[0]*s,q[1]*s,q[2]*s,q[3]*s) def magnitudeSqQ(q): return q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3] def magnitudeQ(q): return sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]) def conjugateQ(q): return (q[0],-q[1],-q[2],-q[3]) def multiplyQ(q1,q2): w1,x1,y1,z1 = q1[0],q1[1],q1[2],q1[3] w2,x2,y2,z2 = q2[0],q2[1],q2[2],q2[3] return (w1*w2 - x1*x2 - y1*y2 - z1*z2, w1*x2 + x1*w2 + y1*z2 - z1*y2, w1*y2 + y1*w2 + z1*x2 - x1*z2, w1*z2 + z1*w2 + x1*y2 - y1*x2) def sclrprdprtQ(q1,q2): v1=(q1[1],q1[2],q1[3]) v2=(q2[1],q2[2],q2[3]) return q1[0]*q2[0]-vecdot3(v1,v2) def vecprdprtQ(q1,q2): w1,x1,y1,z1 = q1[0],q1[1],q1[2],q1[3] w2,x2,y2,z2 = q2[0],q2[1],q2[2],q2[3] return (w1*x2 + x1*w2 + y1*z2 - z1*y2, w1*y2 + y1*w2 + z1*x2 - x1*z2, w1*z2 + z1*w2 + x1*y2 - y1*x2) def anticommprdprtQ(q1,q2): v1=(q1[1],q1[2],q1[3]) v2=(q2[1],q2[2],q2[3]) w = (0.0,) + vecross(v1,v2) return w def normalizeQ(q): m = magnitudeQ(q) if m==0: return (1.0,0.0,0.0,0.0) return (q[0]/m,q[1]/m,q[2]/m,q[3]/m) def inverseQ(q): m2 = magnitudeSqQ(q) return (q[0]/m2,-q[1]/m2,-q[2]/m2,-q[3]/m2) def dotQ(q1,q2): return q1[0]*q2[0] + q1[1]*q2[1] + q1[2]*q2[2] + q1[3]*q2[3] def fromAngleAxisQ(radians,x,y,z): radians/=2.0 s = sin(radians)/sqrt(x*x+y*y+z*z) return normalizeQ((cos(radians),x*s,y*s,z*s)) def toMatrixQ(q): w,x,y,z = q[0], q[1], q[2], q[3] xx = 2.0*x*x yy = 2.0*y*y zz = 2.0*z*z xy = 2.0*x*y zw = 2.0*z*w xz = 2.0*x*z yw = 2.0*y*w yz = 2.0*y*z xw = 2.0*x*w return ((1.0-yy-zz, xy-zw, xz+yw, 0.0), (xy+zw, 1.0-xx-zz, yz-xw, 0.0), (xz-yw, yz+xw, 1.0-xx-yy, 0.0), (0.0, 0.0, 0.0, 1.0)) def rotateVectorQ(q,v): qw, qx, qy, qz = q[0], q[1], q[2], q[3] x, y, z = v[0], v[1], v[2] ww = qw*qw xx = qx*qx yy = qy*qy zz = qz*qz wx = qw*qx wy = qw*qy wz = qw*qz xy = qx*qy xz = qx*qz yz = qy*qz return (ww*x + xx*x - yy*x - zz*x + 2*((xy-wz)*y + (xz+wy)*z), ww*y - xx*y + yy*y - zz*y + 2*((xy+wz)*x + (yz-wx)*z), ww*z - xx*z - yy*z + zz*z + 2*((xz-wy)*x + (yz+wx)*y)) def spatrot(theta,axis,vector): q = (0.0,) + vector v = vecnormalize(axis) P = (cos(theta/2.0),) + Vecscale(v,sin(theta/2.0)) S = conjugateQ(P) a = multiplyQ(P,q) return vecprdprtQ(a,S) def interpolateQ(q1, q2, s, shortest=True): ca = dotQ(q1,q2) if shortest and ca<0: ca = -ca neg_q2 = True else: neg_q2 = False o = acos(ca) so = sin(o) if (abs(so)<=1E-12): return copyQ(q1) a = sin(o*(1.0-s)) / so b = sin(o*s) / so if neg_q2: return subQ(scaleQ(q1,a),scaleQ(q2,b)) else: return addQ(scaleQ(q1,a),scaleQ(q2,b))