-
Notifications
You must be signed in to change notification settings - Fork 0
/
FindGroups.py
75 lines (66 loc) · 2.21 KB
/
FindGroups.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import time, math
import matplotlib
import matplotlib.pyplot as plt
from Fractals.Sierpinski3d import *
from Fractals.Utilities import gen_all_params
import itertools
import matplotlib.animation as animation
import copy
matplotlib.use("Qt4Agg")
shifts = np.ones((8,3))
for i,cube in enumerate(range(8)):
if cube%2 == 1:
shifts[i,0] = -1
if (cube//2)%2 == 1:
shifts[i,1] = -1
if (cube//4)%2 == 1:
shifts[i,2] = -1
print(shifts)
mats = []
for i in range(48):
mats.append(O(i).rotation_matrix())
invs = np.zeros(48,dtype = int)
for i, Oi in enumerate(mats):
for j, Oj in enumerate(mats):
res = np.matmul(Oi,Oj)
if (res == np.eye(3)).all():
invs[i] = j
break
n = 2
for params in gen_all_params(n):
print(params)
frac = Sierpinski3d(params)
relevant_params = []
for p in params:
if p != -1:
relevant_params.append(p)
for di in range(48):
rotated_params = -np.ones(8,dtype=int)
print("\n\n Rotation: ", di)
for l, transform_index in enumerate(zip(frac.transforms,relevant_params)):
T, p = transform_index
Tp = copy.deepcopy(T)
Tp.matrix = 2*np.matmul(O(di).rotation_matrix(),T.matrix)
Tp.matrix = np.matmul(Tp.matrix,O(invs[di]).rotation_matrix())
#print("Matrix number" , l, ":\n ", Tp.matrix)
for j, mat in enumerate(mats):
if (mat == Tp.matrix).all():
#print("Element ", j)
break
Tp.shift = 2*np.matmul(O(di).rotation_matrix(),T.shift)
#print("Shift " , l, ":\n ", Tp.shift)
for k, shift in enumerate(shifts):
if (shift == Tp.shift).all():
#print("Subcube ", k)
break
rotated_params[k] = j
print(rotated_params)
rotated = Sierpinski3d(list(rotated_params))
ptcld, c = rotated.cubical_complex_midpoints(4 ,color = True)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(ptcld[:,0],ptcld[:,1],ptcld[:,2],c = c)
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
ax.set_zlim(-1,1)
plt.show()