class Perms:
"""
TITLE: Perms
AUTHOR: Bill McKeeman
MODIFIED:2012.01.8
PURPOSE:
(1) generate selected permutations of a set of input vectors
(2) generate vertex sets for selected poly(gon)(hedra)(tope)s
(3) generate edge sets for above
(4) generate face sets for above
(5) display 2-D projections for above
EXAMPLES for Perms:
P = Perms
print P.all([[1,2], [3,4]])
[[1,2], [2,1], [3,4], [4,3]]
print P.odd([7,8,9])
[[7, 9, 8], [8, 7, 9], [9, 8, 7]]
print P.even([[7,8,9]])
[[7, 8, 9], [9, 7, 8], [8, 9, 7]]
print cycles([[7,8,9]])
[[7, 8, 9], [8, 9, 7], [9, 7, 8]]
print signs([[1,2], [3,4]])
[[1, 2], [-1, 2], [1, -2], [-1, -2], [3, 4], [-3, 4], [3, -4], [-3, -4]]
print scale([[1,2]], 0.5)
[[0.5, 1]]
print unique([[1,2],[1,2],[1,3],[1,2]])
[[1,2], [1,3]]
v4cube = [[1,1,1,1]]
scale(signs(v4cube),0.5) # 16 vertices of unit 4-cube centered on origin
DESCRIPTION:
The vector [0,1,2,...n-1] has n factorial unique permutations,
of which half are odd and half are even. Class Perms computes and
these holds these permutations to use in generating various kinds
of permutations of a set of equal-length vectors M = [V0,V1,...].
The input vector elements must be comparable.
The kinds of permutations are all, even, odd, cycles and signs.
The kinds of permutations can be combined by using the result
of one call as input to another.
Repeated, or zero elements in the input will cause some vectors
to occur more than once in the result; the unique function removes
duplicates. A scaling factor can be applied to all the elements.
The idea of generating permutations blows up pretty quickly
limiting the usefulness of this class to short input vectors.
The vertex set for many geometric figures can be computed as
permutations of a small set of initial vertex coordinates.
"""
# len 0 1 2 internal tables
evens = [ [[ ]], [[0]], [[0,1]] ] # even permutations of ints
odds = [ [[ ]], [[ ]], [[1,0]] ] # odd permutations of ints
# Even and odd permutations of vectors [], [0], and [0,1]
# are held in class variables.
# When an input vector is longer than 2, function _extendTables
# increases the permutation tables to accommodate the new size.
def all(self, M): # input data vectors M = [V,V,...]
# print "all M =", M
n = len(M[0]) # length of an input vector
self._extendTables(n) # perhaps extend internal tables
res = self._mapData(self.evens[n]+self.odds[n], M) # apply all perms
return res
def even(self, M): # input data vectors
# print "even M =", M
n = len(M[0]) # length of an input vector
self._extendTables(n) # perhaps extend internal tables
res = self._mapData(self.evens[n], M) # apply even perms
return res
def odd(self, M): # input data vectors
# print "odd M =", M
n = len(M[0]) # length of an input vector
self._extendTables(n) # perhaps extend internal tables
res = self._mapData(self.odds[n], M) # apply odd perms
return res
def _mapData(self, ps, V): # data and perms
res = range(len(V)*len(ps)) # preallocate result matrix
ro = 0; # index of each result
for r in range(len(V)): # all input vectors
data = V[r] # rth input vector
for u in range(len(ps)): # each permutation
pm = ps[u] # a permutation
t = range(len(data)) # preallocate result vector
for v in range(len(data)): # each result element
t[v] = data[pm[v]] # permute element
res[ro] = t # permuted vector
ro = ro + 1 # ready for next
return res # element-by-element mapping
def _extendTables(self, n): # extend internal tables
lastn = len(self.evens)-1 # if needed
#print "lastn =", lastn
if n <= lastn: return # tables already OK
lastevens = self.evens[lastn] # [ [0,1] ]
lastodds = self.odds[lastn] # [ [1,0] ]
lastperms = lastevens+lastodds # [ [0,1], [1,0] ]
lp = len(lastperms) # 2
# print "lp =", lp
# preallocate for new evens and odds
m = factorial(lastn+1)//2 # 6//2 == 3
# print "m =", m
newevens = range(m) # [0,1,2]
ne = 0
newodds = newevens[:] # copy
no = 0
for ie in range(lp): # 0...1
te = lastperms[ie] # [0,1]
p = ie < lp/2 # evens are first half
for iee in range(len(te)+1)[::-1]: # start at end
newv = te[:iee] + [lastn] + te[iee:]
if p: # even
newevens[ne] = newv
ne = ne + 1
p = not p
else: # odd
newodds[no] = newv
no = no + 1
p = not p
self.evens = self.evens + [newevens] # extend even table
self.odds = self.odds + [newodds] # extend odd table
self._extendTables(n) # in case more is needed
def scale(M, x): # input data vectors
U = M[:] # copy
for vtx in U:
for e in range(len(vtx)):
vtx[e] = float(vtx[e]*x) # avoid symbolic results
return U
def cycles(M): # input data vectors
# print "cycle M =", M
n = len(M[0]) # length of an input vector
res = range(len(M)*n) # preallocate result
ir = 0
for v in M: # each input vector
for h in range(len(v)): # each cycle
res[ir] = v; # a result
ir = ir + 1
v = v[1:] + [v[0]] # move front to back
return res
def signs(M): # input data vectors
# print "signs M =", M
n = len(M[0]) # length of an input vector
res = M*2**n # preallocate result
ir = 0
for v in M: # each input vector
for k in range(2**n): # each sign change
ie = 0 # elements left to right
u = v[:]
while k > 0: # still something to flip
if k%2 == 1:
u[ie] = -u[ie] # flip sign
k = k//2 # discard bit
ie = ie + 1 # ready for next element
res[ir] = u # a result
ir = ir + 1 # ready for next result
return res
def unique(M): # input data vectors
# print "unique M =", M
U = M[:] # do not modify input
h = 0 # start lo
lim = len(U)
while h < lim: # lim may change
k = lim-1 # eliminate from hi
while k > h:
if U[h]==U[k]: # found duplicate
del U[k] # remove duplicate
lim = lim-1 # note shortened vector
k = k-1 # continue toward lo
h = h+1 # check next input
return U
print "finished compiling Perms"
finished compiling Perms finished compiling Perms |
# Test Perms
P = Perms()
print P.all([[0, 1, 2]])
print P.odd([[0, 1, 2]])
print P.even([[0, 1, 2]])
print signs([[0, 1, 2]])
print unique(signs([[0, 1, 2]]))
print cycles([[0, 1, 2]])
print unique(cycles(signs([[0, 1, 2]])))
print scale([[1, 1, 1]], 1/2)
print "finished testing Perms"
[[0, 1, 2], [2, 0, 1], [1, 2, 0], [0, 2, 1], [1, 0, 2], [2, 1, 0]] [[0, 2, 1], [1, 0, 2], [2, 1, 0]] [[0, 1, 2], [2, 0, 1], [1, 2, 0]] [[0, 1, 2], [0, 1, 2], [0, -1, 2], [0, -1, 2], [0, 1, -2], [0, 1, -2], [0, -1, -2], [0, -1, -2]] [[0, 1, 2], [0, -1, 2], [0, 1, -2], [0, -1, -2]] [[0, 1, 2], [1, 2, 0], [2, 0, 1]] [[0, 1, 2], [1, 2, 0], [2, 0, 1], [0, -1, 2], [-1, 2, 0], [2, 0, -1], [0, 1, -2], [1, -2, 0], [-2, 0, 1], [0, -1, -2], [-1, -2, 0], [-2, 0, -1]] [[0.5, 0.5, 0.5]] finished testing Perms [[0, 1, 2], [2, 0, 1], [1, 2, 0], [0, 2, 1], [1, 0, 2], [2, 1, 0]] [[0, 2, 1], [1, 0, 2], [2, 1, 0]] [[0, 1, 2], [2, 0, 1], [1, 2, 0]] [[0, 1, 2], [0, 1, 2], [0, -1, 2], [0, -1, 2], [0, 1, -2], [0, 1, -2], [0, -1, -2], [0, -1, -2]] [[0, 1, 2], [0, -1, 2], [0, 1, -2], [0, -1, -2]] [[0, 1, 2], [1, 2, 0], [2, 0, 1]] [[0, 1, 2], [1, 2, 0], [2, 0, 1], [0, -1, 2], [-1, 2, 0], [2, 0, -1], [0, 1, -2], [1, -2, 0], [-2, 0, 1], [0, -1, -2], [-1, -2, 0], [-2, 0, -1]] [[0.5, 0.5, 0.5]] finished testing Perms |
"""
TITLE: Vertices
AUTHOR: Bill McKeeman
MODIFIED:2012.01.08
PURPOSE: Generate selected vertex sets for polytopes.
METHOD: Most entries are named by the Schlafli symbols.
For example, if a polyhedron vertex is formed of 3 squares,
the symbol is s43 (In mathematics texts, the symbol is written {4,3}.)
The correspondence between familiar names and symbols is given
in the comments below.
REF:
http://www.mathworks.com/matlabcentral/fileexchange/18523-polytopes
http://en.wikipedia.org/wiki/Schlafli_symbol
http://en.wikipedia.org/wiki/Simplex
http://mathworld.wolfram.com/Polytope.html
http://en.wikipedia.org/wiki/Convex_regular_4-polytope
The Derived Polytopes in Euclidian N-Space, W.M. McKeeman, MS Thesis, GWU 1961
EXAMPLE:
print s43() # vertices of cube
"""
tau = float((sqrt(5)+1)/2) # golden ratio
def polygon(n): # any polygon in 2-D
a = 2*pi/n # Schlafli symbol {n}
return [[cos(a*i), sin(a*i)] for i in range(n)]
def s3(): # triangle {3}
n = 2
V0 = [1,0] # master vertex
y1 = sqrt(1-(-1/n)^2) # pythagoras
V1 = [-1/n, y1] # V0.V1 = -1/n
y2 = -sqrt(1-(-1/n)^2)
V2 = [-1/n, y2]
vs = [V0, V1, V2]
return vs
#return [[1,0], [-1/2, sqrt(3)/2], [-1/2, -sqrt(3)/2]]
#return cycles([[0,0,1]]) # in 3-D
def s33(): # tetrahedron {3,3}
# xi^2 + yi^2 + zi^2 = 1 # Pythagoras in 3-D
# xi*xj + yi*yj + zi*zj = -1/n # inner products unit vectors
n = 3
x0 = 1; y0 = 0; z0 = 0
V0 = [x0, y0, z0] # master vertex V0 V1 V2 V3
# x0 x1 x2 x3
x1 = -1/n # V0.V1 = -1/n y0 y1 y2 y3
x2 = -1/n # V0.V2 = -1/n z0 z1 z2 z3
x3 = -1/n # V0.V3 = -1/n
# V0 V1 V2 V3
z1 = 0 # K known K K K K
# U unknown 0 U U U
# 0 zero 0 0 U U
y1 = sqrt(1 - x1^2 - z1^2) # pythagoras V0 V1 V2 V3
V1 = [x1, y1, z1] # K K K K
# 0 K U U
# 0 0 U U
y2 = (-1/n - x1*x2)/y1 # V1.V2 = -1/n V0 V1 V2 V3
y3 = (-1/n - x1*x3)/y1 # V1.V3 = -1/n K K K K
# 0 K K K
# 0 0 U U
z2 = sqrt(1-x2^2-y2^2) # pythagoras V0 V1 V2 V3
z3 = -sqrt(1-x3^2-y3^2) # K K K K
V2 = [x2, y2, z2] # 0 K K K
V3 = [x3, y3, z3] # 0 0 K K
return [V0, V1, V2, V3] # tetrahedron
#return [[1, 0, 0],
# [-1/3, sqrt(8)/3, 0],
# [-1/3, -sqrt(2)/3, sqrt(6)/3],
# [-1/3, -sqrt(2)/3, -sqrt(6)/3]]
#return cycles([[0,0,0,1]]) # in 4-D
def s333(): # simplex {3,3,3}
# xi^2 + yi^2 + zi^2 + wi^2 = 1 # Pythagoras in 4-D
# xi*xj + yi*yj + zi*zj + wi*wj = -1/n # inner products unit vectors
n = 4
x0 = 1; y0 = 0; z0 = 0; w0 = 0;
V0 = [x0, y0, z0, w0] # master vertex
# V0 V1 V2 V3 V4
x1 = -1/n # V0.V1 = -1/n 1 K K K K
x2 = -1/n # V0.V2 = -1/n 0 y1 y2 y3 y4
x3 = -1/n # V0.V3 = -1/n 0 z1 z2 z3 z4
x4 = -1/n # V0.V4 = -1/n 0 w1 w2 w3 w4
# V0 V1 V2 V3 V4
z1 = 0; w1 = 0 # 1 K K K K
y1 = sqrt(1 - x1^2 - z1^2 - w1^2) # 0 K y2 y3 y4
V1 = [x1, y1, z1, w1] # 0 0 z2 z3 z4
# 0 0 w2 w3 w4
# V0 V1 V2 V3 V4
# 1 K K K K
y2 = (-1/n - x1*x2)/y1 # V1.V2 = -1/n 0 K K K K
y3 = (-1/n - x1*x3)/y1 # V1.V3 = -1/n 0 0 z2 z3 z4
y4 = (-1/n - x1*x4)/y1 # V1.V4 = -1/n 0 0 w2 w3 w4
w2 = 0 # V0 V1 V2 V3 V4
z2 = sqrt(1 - x2^2 - y2^2) # pythagorus 1 K K K K
V2 = [x2, y2, z2, w2] # 0 K K K K
# 0 0 K z3 z4
# 0 0 0 w3 w4
z3 = (-1/n - x2*x3 - y2+y3)/z2 # V2.V3 = -1/n V0 V1 V2 V3 V4
z4 = (-1/n - x2*x4 - y2+y4)/z2 # V2.V3 = -1/n 1 K K K K
# 0 K K K K
w3 = sqrt(1 - x3^2 - y3^2 - z3^2) # 0 0 K K K
w4 = -sqrt(1 - x4^2 - y4^2 - z4^2) # 0 0 0 w3 w4
V3 = [x3, y3, z3, w3]
V4 = [x4, y4, z4, w4]
return [V0, V1, V2, V3, V4] # 4-simplex
#return [[ 1, 0, 0, 0],
# [-1/4, sqrt(15)/4, 0, 0],
# [-1/4, -sqrt(15)/12, sqrt(5/6), 0],
# [-1/4, -sqrt(15)/12, -3*sqrt(5/6)/8, 5*sqrt(11/6)/8],
# [-1/4, -sqrt(15)/12, -3*sqrt(5/6)/8, -5*sqrt(11/6)/8]]
#return cycles([[0,0,0,0,1]]) # and so on in (n+1)-D
def segment(): # line segment
return signs([[1]])
def s4(): # square {4}
return signs([[1,1]])
def s43(): # cube {4,3}
return signs([[1,1,1]])
def s433(): # hypercube {4,3,3}
return signs([[1,1,1,1]]) # and so on
def s34(): # octahedron 3-cross
return unique(signs(cycles([[0,0,1]])))
def s334(): # 4-cross 16-cell
return unique(signs(cycles([[0,0,0,1]]))) # and so on
def s343(): # hyperdiamond 24-cell
return unique(signs(P.all([[0,0,1,1]])))
def s5(): # pentagon
return [[0,3*tau/2,-1/2],
[tau/2,tau+1/2,-1], [-tau/2,tau+1/2,-1],
[1/2,tau,-1-tau/2], [-1/2,tau,-1-tau/2]]
def s35(): # icosahedron
return unique(signs(cycles([[0,1,tau]])))
def s335(): # hyper icosahedron
v4i = unique(signs(cycles([[2,0,0,0]]))) + \
signs([[1,1,1,1]]) + \
unique(signs(P.even([[tau, 1, 1/tau, 0]])))
return v4i # last of the series
def s53(): # dodecahedron
return unique(signs(cycles([[0,1/tau,tau]])+[[1,1,1]]))
def s533(): # hyper dodecahedron 120-cell
vs1 = [[2, 2, 0, 0],
[float(sqrt(5)), 1, 1, 1],
[tau, tau, tau, tau^-2],
[tau^-1, tau^-1, tau^-1, tau^2]]
v1 = unique(P.all(vs1))
vs2 = [[tau^2, tau^-2, 1, 0],
[float(sqrt(5)), tau^-1, tau, 0],
[2, 1, tau, tau^-1]]
v2 = P.even(vs2)
v4d = unique(signs(v1+v2))
return v4d # last of the series
def buckyball():
d = lambda a,b: (a + b*tau)/2
bv = [[d(0,0),d(0,3),d(1,0)], [d(1,0),d(0,2),d(2,1)], [d(2,0),d(0,1),d(1,2)]]
return unique(signs(cycles(bv))) # bucky ball
def rotate(p, i0, i1, radians): # rotate about i0,i1 axes
s = sin(radians)
c = cos(radians)
res = p[:] # preallocate output
for i in range(len(p)):
vs = p[i][:] # copy point
x = vs[i0]
y = vs[i1]
vs[i0] = c*x + s*y # rotate pairs
vs[i1] = -s*x + c*y
res[i] = vs # record result
return res
print "finished compiling vertices"
finished compiling vertices finished compiling vertices |
# test Vertices
print s3()
print s33()
print s333()
#print rotate(s33(), 0, 1, 0.1)
#print s43()
vs = s335()
print "hyper icosahedron has ", len(vs), " vertices"
vs = s533()
print "hyper dodecahedron has ", len(vs), " vertices" # this takes awhile
print "finished testing vertices"
[[1, 0], [-1/2, 1/2*sqrt(3)], [-1/2, -1/2*sqrt(3)]] [[1, 0, 0], [-1/3, 2/3*sqrt(2), 0], [-1/3, -1/3*sqrt(2), sqrt(2/3)], [-1/3, -1/3*sqrt(2), -sqrt(2/3)]] [[1, 0, 0, 0], [-1/4, 1/4*sqrt(15), 0, 0], [-1/4, -1/12*sqrt(15), sqrt(5/6), 0], [-1/4, -1/12*sqrt(15), -3/8*sqrt(5/6), 5/8*sqrt(11/6)], [-1/4, -1/12*sqrt(15), -3/8*sqrt(5/6), -5/8*sqrt(11/6)]] hyper icosahedron has 120 vertices hyper dodecahedron has 600 vertices finished testing vertices [[1, 0], [-1/2, 1/2*sqrt(3)], [-1/2, -1/2*sqrt(3)]] [[1, 0, 0], [-1/3, 2/3*sqrt(2), 0], [-1/3, -1/3*sqrt(2), sqrt(2/3)], [-1/3, -1/3*sqrt(2), -sqrt(2/3)]] [[1, 0, 0, 0], [-1/4, 1/4*sqrt(15), 0, 0], [-1/4, -1/12*sqrt(15), sqrt(5/6), 0], [-1/4, -1/12*sqrt(15), -3/8*sqrt(5/6), 5/8*sqrt(11/6)], [-1/4, -1/12*sqrt(15), -3/8*sqrt(5/6), -5/8*sqrt(11/6)]] hyper icosahedron has 120 vertices hyper dodecahedron has 600 vertices finished testing vertices |
"""
TITLE: Edges & Edge Length
AUTHOR: Bill McKeeman
MODIFIED:2012.01.20
PURPOSE: compute edges given the vertices of a poly(gon)(hedron)(choron)
METHOD: The distance between pairs of vertices is minimum for edges.
Find the minimum and then collect the vertex pairs for the edges.
"""
def EdgeLength(vs):
nd = len(vs[0]) # number of dimensions
minedge = infinity # find edge length
for v0 in range(len(vs)-1):
for v1 in range(v0+1,len(vs)):
distsq = float(sum([(vs[v0][iv]-vs[v1][iv])**2 for iv in range(nd)]))
if distsq < minedge:
minedge = distsq
# print "edge length = ", minedge
return minedge
def Edges(vs):
nd = len(vs[0]) # number of dimensions
for v0 in range(len(vs)): # get rid of symbolic vals
for iv in range(nd):
vs[v0][iv] = float(vs[v0][iv])
minedge = EdgeLength(vs) # get edge length
edg = [] # no edges yet
for v0 in range(len(vs)-1): # try all pairs
for v1 in range(v0+1,len(vs)):
distsq = sum([(vs[v0][iv]-vs[v1][iv])**2 for iv in range(nd)])
# print "distsq = ", distsq
if abs(distsq - minedge) < minedge/2: # vertices are closest
edg += [[v0, v1]] # record new edge
# print "edge count = ", len(edg)
return edg
print "finished compiling Edges"
finished compiling Edges finished compiling Edges |
# test edges
v = Edges(s3())
print v
v = Edges(s33())
print v
v = Edges(s333())
print v
[[0, 1], [0, 2], [1, 2]] [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]] [[0, 1], [0, 2], [0, 3], [0, 4], [1, 2], [1, 3], [1, 4], [2, 3], [2, 4], [3, 4]] [[0, 1], [0, 2], [1, 2]] [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]] [[0, 1], [0, 2], [0, 3], [0, 4], [1, 2], [1, 3], [1, 4], [2, 3], [2, 4], [3, 4]] |
"""
FILE: lineareq.py
PURPOSE: solve nXn system of linear equations
AUTHOR: Bill McKeeman
MODIFIED:2012.01.27
METHOD: Mx = b
nXn matrix M + rhs vector b is nX(n+1) input
Subtract rows to LR triangular form.
Back-solve to get x.
"""
def cvmul(c, V): # scale vector V by constant c
return [c*V[iv] for iv in range(len(V))]
def vvsub(V0, V1): # vector V0 - vector V1
return [V0[iv]-V1[iv] for iv in range(len(V1))]
def dot(V0, V1):
return sum([V0[iv]*V1[iv] for iv in range(len(V1))])
def lr(M): # reduce to upper right triangle
# print "M =", M
R = M[:]
for ir in range(len(R)):
if R[ir][ir] == 0: # must have non zero pivot
pivrow = -1;
for iv in range(ir+1, len(R)): # find pivot
if R[iv][ir] != 0:
pivrow = iv
if pivrow != -1: # exchange equations
tmp = R[ir]
R[ir] = R[pivrow]
R[pivrow] = tmp
else:
return "unsolvable" # cannot solve
for iv in range(ir+1, len(R)): # zero below pivot
# print "before iv =", iv, "R[iv] =", R[iv]
R[iv] = vvsub(R[iv], cvmul(R[iv][ir]/R[ir][ir], R[ir]))
# print "after iv =", iv, "R[iv] =", R[iv]
return R
def back(lr): # back solve for x
# print "lr =", lr
n = len(lr)
res = range(n) # preallocate result
for ir in range(n-1,-1,-1):
t = sum(lr[ir][ix]*res[ix] for ix in range(ir+1,n))
if lr[ir][ir] == 0: return "unsolvable"
res[ir] = (lr[ir][n] - t)/lr[ir][ir]
return res
print "finished compiling lineareq"
finished compiling lineareq finished compiling lineareq |
"""
test lineareq
"""
print sum([1,2,3])
print cvmul(3, [1,2,3])
print vvsub([2,3], [1,1])
print dot([1,2], [4,5])
res = lr([[1, 1, 6], [1, -1, 4]])
print res
print back(res)
print "should be = ", [5,1]
res = lr([[0,1,1,1], [0,2,3,1], [0,3,4,1]])
print res
res = lr([[1,2,3,1], [1,2,3,1], [1,2,3,1]])
print res
res = lr([[1,1,1,4], [0,1,1,5], [2,2,2,6]])
print res
6 [3, 6, 9] [1, 2] 14 [[1, 1, 6], [0, -2, -2]] [5, 1] should be = [5, 1] unsolvable unsolvable unsolvable 6 [3, 6, 9] [1, 2] 14 [[1, 1, 6], [0, -2, -2]] [5, 1] should be = [5, 1] unsolvable unsolvable unsolvable |
"""
TITLE: Faces
AUTHOR: Bill McKeeman
MODIFIED:2012.01.20
PURPOSE: compute faces given the vertices of a poly(hedron)(choron)
METHOD: In 3-D, the vertices of a face all lie in the same plane,
x/a + y/b + z/c = 1 (which does not go through the origin)
and are connected by edges.
If the plane is parallel to an axis, some of a,b,c may be inf.
For each vertex x0,y0,z0
(1) find two connected vertices x1,y1,z1 and x2,y2,x2
(2) plug the 3 sets of vertex values into the linear equation,
(3) solve the resulting three equations for 1/a, 1/b and 1/c:
(a) discard the equation if it is a duplicate or passes through origin
(4) find all other vertices that solve the equation
(5) these vertices constitute a face or a slice.
(6) discard the slices (how?)
The result is a vector of a vector of vertices.
"""
def Faces(vs): # input is vertex set
nd = len(vs[0]) # number of dimensions
edg = Edges(vs)
# print "edges =", edg
fl = [-1 for iz in range(len(vs))]
zm = [fl[:] for iz in range(len(vs))] # nXn array of no connection
for ei in range(len(edg)): # connection matrix
v0 = edg[ei][0]
v1 = edg[ei][1]
zm[v0][v1] = ei # remember which edge
zm[v1][v0] = ei
# for ei in range(len(zm)): print ei, zm[ei]
trios = []
for st in range(len(vs)): # each starting vertex
nbr = []
for nb in range(len(vs)): # find neighbors
if zm[st][nb] >= 0: # st-nb is an edge
nbr += [nb]
trios += [nbr]
# print "trios ="
# for tr in range(len(trios)): print tr, trios[tr]
eqns = [[]]
fcs = []
for st in range(len(trios)): # each starting vtx
for tr in range(len(trios[st])-1): # pairs of connected edges
for te in range(tr+1, len(trios[st])):
# print "st =", st, "tr =", tr, "te =", te
eqns = [vs[st] + [1], vs[trios[st][tr]] + [1], vs[trios[st][te]] + [1]]
# print "eqns =", eqns
lrf = lr(eqns)
# print "lr =", lrf
if lrf != "unsolvable": # otherwise skip it
coef= back(lrf)
# print "coef=", coef
if coef != "unsolvable":
fac = [];
for fc in range(len(vs)):
vv = vs[fc]
sol = dot(vv, coef) - 1
# print "sol =", sol
if abs(sol) < 0.01: # effectively 0
fac += [fc]
# print "fac =", fac
new = true; # optimistic
for si in range(len(fcs)):
if fac == fcs[si]:
new = false; # seen before
if new:
fcs += [fac] # record new face
# print "fcs =", fcs
return fcs
print "finished compiling Faces"
finished compiling Faces finished compiling Faces |
# test Faces # for cube 0 1
fig = s43() # 2 3 4 5
fcs = Faces(fig) # 6 7
print "faces+slices =", len(fcs)
for ip in range(11):
for ir in range(len(fcs)):
if len(fcs[ir]) == ip:
print "cube face has", ip, "vertices: ", fcs[ir]
# does not work for tetrahedron (imbedded in 4-space)
print "tetrahedron has", len(Faces(s33())), "faces"
print "octahedron has", len(Faces(s34())), "faces"
print "cube has", len(Faces(s43())), "faces"
print "dodecahedron has", len(Faces(s53())), "faces"
print "icosahedron has", len(Faces(s35())), "faces and slices"
#print "hyper cube has", len(Faces(s433())), "faces"
print "bucky has", len(Faces(buckyball())), "faces"
faces+slices = 6 cube face has 4 vertices: [0, 1, 2, 3] cube face has 4 vertices: [0, 1, 4, 5] cube face has 4 vertices: [0, 2, 4, 6] cube face has 4 vertices: [1, 3, 5, 7] cube face has 4 vertices: [2, 3, 6, 7] cube face has 4 vertices: [4, 5, 6, 7] tetrahedron has 4 faces octahedron has 8 faces cube has 6 faces dodecahedron has 12 faces icosahedron has 32 faces and slices bucky has 32 faces faces+slices = 6 cube face has 4 vertices: [0, 1, 2, 3] cube face has 4 vertices: [0, 1, 4, 5] cube face has 4 vertices: [0, 2, 4, 6] cube face has 4 vertices: [1, 3, 5, 7] cube face has 4 vertices: [2, 3, 6, 7] cube face has 4 vertices: [4, 5, 6, 7] tetrahedron has 4 faces octahedron has 8 faces cube has 6 faces dodecahedron has 12 faces icosahedron has 32 faces and slices bucky has 32 faces |
"""
TITLE: plotpoly
AUTHOR: Bill McKeeman
MODIFIED:2012.01.08
PURPOSE: make graphics object to plot a poly(gon)(hedron)(choron)
"""
def plotpoly(vs): # input is vertex set
nd = len(vs[0]) # number of dimensions
edg = Edges(vs) # pairs of vertices
# print "edge count = ", len(edg)
dig = ['0','1','2','3','4','5','6','7','8','9','a','b','c','d','e','f']
hx = lambda(v): dig[v//16] + dig[v%16] # two-digit hex
cx = lambda(v): '#'+hx(v)+hx(255-v)+'00' # color
nvs = len(vs)
x = range(nvs) # preallocate
y = range(nvs)
z = range(nvs) # for 3-D depth of field
w = range(nvs) # for 4-D depth of field
c = 'blue' # default
maxz = 0
maxw = 0
for ip in range(nvs): # x,y coords of polytope
x[ip] = vs[ip][0]
y[ip] = vs[ip][1]
if nd> 2: # 3rd D
z[ip] = vs[ip][2]
#print "z[", ip, "]=", z[ip]
if z[ip] > maxz:
maxz = z[ip]
if nd > 3: # 4th D
w[ip] = vs[ip][3]
if w[ip] > maxw:
maxw = w[ip]
# print "maxz =", maxz
# print "maxw =", maxw
ed = edg[0]; # first edge
if nd> 2: # polyhedra and above
av = (z[ed[0]]+z[ed[1]])/2 # mid-point
depth = (av+maxz)/maxz/2 # 0 <= depth <= 1
else:
depth = 2 # polygon, no depth needed
if nd> 3: # polytope
av = (w[ed[0]]+w[ed[1]])/2 # mid-point
depth4 = (av+maxw)/maxw/2 # 0 <= depth4 <= 1
# print "depth4 = ", depth4
# print "color =", floor(depth4*255)
c = cx(floor(depth4*255))
L = line([(x[ed[0]],y[ed[0]]), (x[ed[1]],y[ed[1]])], thickness=depth+0.2, color=c)
for ei in range(1,len(edg)): # rest of edges
ed = edg[ei] # e.g. [1,2]
av = (z[ed[0]]+z[ed[1]])/2 # mid-point
if nd> 2:
depth = (av+maxz)/maxz/2 # positive
else:
depth = 2;
# print "depth =", depth
if nd > 3:
av = (w[ed[0]]+w[ed[1]])/2 # mid-point
depth4 = (av+maxw)/maxw/2
# print "depth4 =", depth4
# print "color =", floor(depth4*255)
c = cx(floor(depth4*255))
L += line([(x[ed[0]],y[ed[0]]), (x[ed[1]],y[ed[1]])], thickness=depth+0.2, color=c)
return L # graphics object
print "finished compiling plotpoly"
finished compiling plotpoly finished compiling plotpoly |
# looking for a pattern
V3 = [[ 1, 0],
[-1/2, sqrt(3/4)],
[-1/2, -sqrt(3/4)]]
V33 = [[ 1, 0, 0],
[-1/3, sqrt(8/9), 0],
[-1/3, -sqrt(2/9), sqrt(2/3)],
[-1/3, -sqrt(2/9), -sqrt(2/3)]]
V333 = [[ 1, 0, 0, 0],
[-1/4, sqrt(40/48), 0, 0],
[-1/4, -sqrt(5/48), sqrt(40/48), 0],
[-1/4, -sqrt(5/48), -sqrt(15/128), sqrt(275/384)],
[-1/4, -sqrt(5/48), -sqrt(15/128), -sqrt(275/384)]]
print "v333 =", V333
plotpoly(rotate(rotate(rotate(V333, 0,1, 0.1), 1,2, 0.2),0,3, 0.2))
v333 = [[1, 0, 0, 0], [-1/4, sqrt(5/6), 0, 0], [-1/4, -1/4*sqrt(5/3), sqrt(5/6), 0], [-1/4, -1/4*sqrt(5/3), -1/8*sqrt(15/2), 5/8*sqrt(11/6)], [-1/4, -1/4*sqrt(5/3), -1/8*sqrt(15/2), -5/8*sqrt(11/6)]] v333 = [[1, 0, 0, 0], [-1/4, sqrt(5/6), 0, 0], [-1/4, -1/4*sqrt(5/3), sqrt(5/6), 0], [-1/4, -1/4*sqrt(5/3), -1/8*sqrt(15/2), 5/8*sqrt(11/6)], [-1/4, -1/4*sqrt(5/3), -1/8*sqrt(15/2), -5/8*sqrt(11/6)]] |
# test plotpoly for triangle
fig = s3()
L = plotpoly(fig)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for tetrahedron
fig = s33()
fig = rotate(fig, 0,1, 0.2)
fig = rotate(fig, 0,2, 0.2)
fig = rotate(fig, 1,2, 0.3)
L = plotpoly(fig)
plot(L, aspect_ratio='equal')
|
|
# test plotpoly for hyper tetrahedron
fig = s333()
fig = rotate(fig, 0,1, 0.2)
fig = rotate(fig, 0,2, 0.3)
fig = rotate(fig, 1,2, 0.2)
fig = rotate(fig, 0,3, 0.1)
fig = rotate(fig, 1,3, 0.1)
L = plotpoly(fig)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for 4-D cube
fig= s433()
fig= rotate(fig, 0,1, 0.3) # an arbitrary 4-D rotation
fig= rotate(fig, 0,2, 0.5)
fig= rotate(fig, 1,3, 0.6)
L = plotpoly(fig)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for cube
cb = s43()
cb = rotate(cb, 0,1, 0.3)
cb = rotate(cb, 1,2, 0.5)
L = plotpoly(cb)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for octahedron
cb = s34()
cb = rotate(cb, 0,1, 0.3)
cb = rotate(cb, 1,2, 0.5)
L = plotpoly(cb)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for dodacahedron
cb = s53()
cb = rotate(cb, 0,1, 0.3)
cb = rotate(cb, 1,2, 0.5)
L = plotpoly(cb)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for icosacahedron
cb = s35()
cb = rotate(cb, 0,1, 0.4)
cb = rotate(cb, 0,2, 0.4)
L = plotpoly(cb)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for buckyball
bb = buckyball()
bb = rotate(bb, 0,1, 0.3)
bb = rotate(bb, 1,2, 0.5)
L = plotpoly(bb)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for 4-cross
vs = s334()
vs = rotate(vs, 0,1, 0.3)
vs = rotate(vs, 0,2, 0.5)
vs = rotate(vs, 1,3, 0.6)
L = plotpoly(vs)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for hyper diamond
vs = s343()
vs = rotate(vs, 0,1, 0.3)
vs = rotate(vs, 0,2, 0.5)
vs = rotate(vs, 1,3, 0.6)
L = plotpoly(vs)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for hyper icosahedron
vs = s335()
vs = rotate(vs, 0,1, 0.3)
vs = rotate(vs, 0,2, 0.5)
vs = rotate(vs, 1,3, 0.6)
L = plotpoly(vs)
plot(L, aspect_ratio='equal', axes=false)
|
|
# test plotpoly for hyper dodecahedron
# be patient: this takes awhile
vs = s533()
vs = rotate(vs, 0,1, 0.3)
vs = rotate(vs, 0,2, 0.4)
vs = rotate(vs, 0,3, 0.7)
vs = rotate(vs, 1,2, 0.2)
L = plotpoly(vs)
plot(L, aspect_ratio='equal', axes=false)
|
|
|
|