advanced.sage

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