advanced.sage 7.52 KB
Newer Older
Peter Korcsok's avatar
Peter Korcsok committed
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
def Partitions(List):
  """
    Generates all proper partitions of List.
  """
  if len(List)==0:
    yield []
    return
  if len(List)==1:
    return
  
  for A in Subsets( List[:-1] ):
    if len(A) > 0:
      M = List[:-1] 
      for a in A:
        M.remove(a) 
      AA = A.list() + [List[-1]] 
      for r in Partitions(M):
        yield r + [ AA ]
# Partitions


def Classes(Partition):
  """
    Returns a list determining the class
    where the partition belongs.
  """
  m = max([ max(p) for p in Partition])
  Res = range(m + 1)
  
  for P in Partition:
    m = min(P)
    for p in P:
      Res[p] = m
  
  return Res
# Classes


def IsMinPartition(Part, Gamma):
  """
    Returns True if given partition is
    lexicographically minimal in its class.
  """
  X = Classes(Part)
  
  for r in [ GammaPartition(g, Part) for g in Gamma]:
    x = Classes(r)
    if x < X:
      return False
  
  return True
# IsMinPartition


def GammaPartition(Gamma, Partition):
  """
    Applies permutation Gamma on Partition.
  """
  Res = Set([])
  
  for P in Partition:
    P_ = Set([Gamma[r]-1 for r in P])
    Res = Res.union(Set([P_]))
  
  return Res
# GammaPartition


def PartitionClasses(Partitions, Gamma):
  """
    Returns list of classes of partitions.
  """
  Res = []
  for Partition in Partitions:
    if IsMinPartition(Partition, Gamma):
      Res += [ Set([ GammaPartition(g, Partition) for g in Gamma ]) ]
  
  return Res
# PartitionClasses


def NTuples(n):
  """
    Generates all n-tuples of {1, 2, 3, 4}^n.
  """
  if n==0:
    yield []
    return
  
  for t in NTuples(n-1):
    for i in [1, 2, 3, 4]:
      yield t + [i]
# NTuples


def CountTerminalValues(n):
  """
    Counts all values of n terminals with
    - first item = 1,
    - sum of all items = 0.
  """
  if n<=1:
    return 0

  Count = 0
  for t in NTuples(n-2):
    T = [1] + t
    s = Integers(5)(0)
    for i in T:
      s += i
    if s!=0:
      Count += 1

  return Count
# CountTerminalValues


def TerminalValues(n):
  """
    Generates all values of n terminals with
    - first item = 1,
    - sum of all items = 0.
  """
  if n<=1:
    yield []
    return
  
  for t in NTuples(n-2):
    T = [1] + t
    s = Integers(5)(0)
    for i in T:
      s += i
    if s!=0:
      yield T + [-s]
# TerminalValues


def Compatibility(Partition, TerminalValues):
  """
    Returns 1 if the partition is compatible with
    values on terminals.
  """
  for Class in Partition:
    s = Integers(5)(0)
    for Index in Class:
      s += TerminalValues[Index]
    if s!=0:
      return 0

  return 1
# Compatibility


def Chi(PartitionClasses, TerminalValues):
  """
    Returns the vector of compatibility.
  """
  return [sum(Compatibility(P, TerminalValues) for P in Part) for Part in PartitionClasses]
# Chi


def IsFlow(TerminalValues):
  """
    Returns True if there exists a flow
    with given values on terminals.
  """
  val = [ Integers(5)(i) for i in range(1,5) ]
  for Val in TerminalValues:
    val = [ v + Val for v in val if v + Val > 0 ]
  
  return len(val) > 0
# IsFlow


def EncodeTerminalValues(Values):
  """
    Encodes values of terminals into one number.
  """
  if len(Values) == 1:
    return Integer(Values[-1] - 1)
  
  return 4 * EncodeTerminalValues(Values[:-1]) + Integer(Values[-1] - 1)
# EncodeTerminalValues


def GammaValues(Values, Gamma):
  """
    Applies permutation Gamma on Values.
  """
  Res = []
  for i in range(len(Values)):
    Res += [ Values[Gamma[i] - 1] ]
  
  s = 1/Integers(5)(Res[0])
  if s != 1:
    return [ i * s for i in Res]
  
  return Res
# GammaValues


def Generate(n, Buffer, Verbose, Save=True, Type=3):
  """
    Generates the matrices M_n and M'_n and
    counts their ranks.
    If some matrix is too long the rank is computed
    by parts with matrix buffer of given size.
    
    The parameter Verbose determines the number
    of rows after that a progress message is printed.
    The parameter Save determines whether the result
    matrices are saved on disk.
    The parameter Type determines which matrices
    are generated: 1 - M_n, 2 - M'_n, 3 - both.
  """
  from datetime import datetime
  
  print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "initiating..."
  
  G = DihedralGroup(n)
  Gamma = [ g.tuple() for g in G ]
  
  PartCls = PartitionClasses(Partitions(range(n)), Gamma)
  PartSml = [ range(n-2), [ n-2, n-1 ] ]
  
  Count = CountTerminalValues(n)
  Vals = TerminalValues(n)
  
  Used = []
  for i in range(4^(n-1)):
    Used += [ False ]
  
  Mn = []
  MMn = []
  rn = 0
  rrn = 0
  proc = 0
  verb = 0
  
  print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "processing terminal values..."
  
  for H in Vals:
    proc += 1
    verb += 1
    
    if verb >= Verbose:
      print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "processed", proc, "of", Count, "...", (100 * proc/Count).round(), "%"
      verb = 0
    
    if Used[EncodeTerminalValues(H)]:
      continue
    
    Used[EncodeTerminalValues(H)] = True
    
    if IsFlow(H):
      X = Chi(PartCls, H)
      for G in Gamma:
        Used[EncodeTerminalValues(GammaValues(H, G))] = True
      
      if Type & 1 == 1:
        Mn += [ X ]
        rn += 1
        if rn >= Buffer:
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "processed", proc, "of", Count, "...", (100 * proc/Count).round(), "%"
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "buffer for Mn overflow (", rn, ")..."
          M = Matrix(QQ, Mn)
          rn = M.rank()
          P, L, U = M.LU()
          Mn = [ U[i].list() for i in range(rn) ]
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "buffer for Mn completed (", rn, ")"
      
      if Type & 2 == 2:
        MMn += [ X ]
        rrn += 1
        if rrn >= Buffer:
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "processed", proc, "of", Count, "...", (100 * proc/Count).round(), "%"
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "buffer for MMn overflow (", rrn, ")..."
          M = Matrix(QQ, MMn)
          rrn = M.rank()
          P, L, U = M.LU()
          MMn = [ U[i].list() for i in range(rrn) ]
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "buffer for MMn completed (", rrn, ")"
      
      continue
    
    if (Compatibility(PartSml, H) == 1) and IsFlow(H[:-2]):
      X = Chi(PartCls, H)
      for G in Gamma:
        Used[EncodeTerminalValues(GammaValues(H, G))] = True
      
      if Type & 2 == 2:
        MMn += [ X ]
        rrn += 1
        if rrn >= Buffer:
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "processed", proc, "of", Count, "...", (100 * proc/Count).round(), "%"
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "buffer for MMn overflow (", rrn, ")..."
          M = Matrix(QQ, MMn)
          rrn = M.rank()
          P, L, U = M.LU()
          MMn = [ U[i].list() for i in range(rrn) ]
          print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "buffer for MMn completed (", rrn, ")"
  
  print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "processing terminal values completed"
  
  Mn  = Matrix(Mn)
  MMn = Matrix(MMn)
  
  if Save:
    print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "saving the matrices..."
    
    s = "saves/advanced/" + datetime.now().strftime('%Y-%m-%d_%H-%M-%S') + "_" + n.str() + "_"                   
    if Type & 1 == 1:
      save(Mn, s + "1")
    if Type & 2 == 2:
      save(MMn, s + "2")
  
  print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "counting the ranks..."
  
  rn  = Mn.rank()
  rrn = MMn.rank()
  
  print datetime.now().strftime('%Y-%m-%d %H:%M:%S'), "counting the ranks completed"
  
  return [ n, rn, rrn ]
# Generate