Sage code for brute force generation of uniform permutation networks

unsorting.sage

def swap_sequences_inlined(n, swap_vars):
  f = factorial(n)

  def extend_swap_sequence(seq, equivalence_classes, distrib, swap_vars):
    if len(distrib) * (2^len(swap_vars)) < f:
      # We can't hit all permutations, so fail fast
      return

    if not swap_vars:
      yield (seq, [weight - 1/f for weight in distrib.values()])
      return

    for i in range(len(equivalence_classes)):
      for j in range(i, len(equivalence_classes)):
        if i == j and len(equivalence_classes[i]) == 1:
          continue

        nclasses = list(equivalence_classes)
        if i == j:
          sw = equivalence_classes[i][:2]
          if len(equivalence_classes[i]) > 2:
            nclasses.append(equivalence_classes[i][2:])
            nclasses[i] = sw
        else:
          sw = (equivalence_classes[i][0], equivalence_classes[j][0])
          if len(equivalence_classes[i]) > 1:
            nclasses.append(equivalence_classes[i][1:])
            nclasses[i] = equivalence_classes[i][:1]
          if len(equivalence_classes[j]) > 1:
            nclasses.append(equivalence_classes[j][1:])
            nclasses[j] = equivalence_classes[j][:1]

        if seq:
          # Impose canonical order on independent swaps to reduce redundancies
          if sw < seq[-1] and len(set(sw + seq[-1])) == 4:
            continue

          # Consecutive swaps of the same pair of indices are redundant
          if sw == seq[-1]:
            continue

        ndistrib = dict()
        for pi, weight in distrib.items():
          npi = Permutation(sw).right_action_product(pi)
          ndistrib[npi] = weight * swap_vars[0] + (ndistrib[npi] if npi in ndistrib else 0)
          ndistrib[pi] = weight * (1 - swap_vars[0]) + (ndistrib[pi] if pi in ndistrib else 0)

        for extension in extend_swap_sequence(seq + [sw], nclasses, ndistrib, swap_vars[1:]):
          yield extension

  return extend_swap_sequence([], [tuple(range(1, n+1))], { Permutation((n,)): 1 }, swap_vars)


def combinatorial_prime_ideal_generation(n, distrib, swap_vars):
  # At least n-1 of the swap_vars are equal to 1/2.
  all_prime_ideals = set()
  for halves in Subsets(swap_vars, n-1):
    expanded = [p.subs({v:1/2 for v in halves}) for p in distrib] + [2*v-1 for v in halves]
    # Deduplicate, but for some reason we need to convert to a list rather than a set
    expanded = list(set(expanded))
    all_prime_ideals.update(Ideal(expanded).minimal_associated_primes())

  return [prime for prime in all_prime_ideals if not prime.is_trivial()]


for n, k in [(4,5), (4,6), (5,9), (5,10)]:
  print(n, k)
  # Ensure that the variables are in order and we can decode the output
  gd = PolynomialRing(QQ, "p", k, order='lex').gens_dict()
  swap_vars = [gd[varname] for varname in sorted(gd)]

  for swap_seq, distrib in swap_sequences_inlined(n, swap_vars):
    alarm(20)
    try:
      I = ideal(distrib)
      primes = [prime for prime in I.minimal_associated_primes() if not prime.is_trivial()]
    except AlarmInterrupt:
      primes = combinatorial_prime_ideal_generation(n, distrib, swap_vars)

    cancel_alarm()

    if primes:
      print(swap_seq, primes)

  print("---")