flower-snarks.py

#!/usr/bin/python

DOC = """
Calculate that the number of CDCs of Flower snark $J_k$ is $c16^k \pm O(15^k)$.

We calculate the step matrix $M$, show that 16 is its eigenvalue by finding
its eigenvector. The eigenvector is guessed for the ratios of coefficients
of $M^10 x$ and $M^9 x$ where $x$ is the initial vector. This eigenvector
shares one non-zero coordinate with the final vector.
"""

import sys, os
sys.path.append(os.path.dirname(__file__) + "/..")

from graph_tools.all import *
from sage.all import identity_matrix, QQ, DiGraph, matrix

def make_P(u):
  coords = [ r for r, c in u.dict().keys() ]
  assert min(coords) == 1 and u[1, 0] == 1

  P = matrix(QQ, u.nrows(), u.nrows(), sparse=True)
  P[:, 0] = u
  P[0, 1] = 1
  for i in range(2, P.nrows()):
    P[i, i] = 1
    P[1, i] = -P[i, 0]

  return P

def deflate_eigenvector(M, u):
  coord = min(( r for r in range(u.nrows()) if u[r, 0] > 0 ), key=lambda r: u[r, 0])
  a = M[coord, :]
  return M - u * a / u[coord, 0]


def bound_spectral_radius(M, bound):
  def norm(A):
    """The matrix 1-norm. Max of the sum of columns.

    We do not use the norms from Sage because they are slow
    for sprase matrices.
    """
    sums = [ 0 ] * A.ncols()
    for (r, c), v in A.dict().items():
      sums[c] += abs(v)

    return max(sums)

  X = M
  p = 1

  while True:
    b = bound**p
    n = norm(X)
    print("  p=%i upper bound is %0.3lg" % (p, n**(1/p)))

    if n <= b: break

    X = X**2
    p *= 2

  print("  %i <= %i**%i" % (n, bound, p))


def test_eigenvalue(M, x):
  H = M - x*identity_matrix(QQ, M.ncols(), sparse=True)
  if H.is_singular():
    print("  Value %i is an eigenvalue of M with geometric multiplicity %i" %
      (x, H.ncols() - H.rank()))


def matrix_to_rev_graph(M):
  assert M.ncols() == M.nrows()
  return DiGraph([range(M.ncols()), M.dict().keys()],
    format="vertices_and_edges", loops=True)


if __name__ == "__main__":
  print(DOC)
  print("Stabilizing...")
  S = FlowerSnark.stabilize(CircuitDoubleCover, jordan=False)
  print("Done")
  M = S["step_matrix"]

  print("\nThe rank of the step matrix M: %s" % M.rank())
  upto = 16
  print("Testing if integers up to %i are eigenvalues:" % upto)
  for x in range(upto + 1):
    test_eigenvalue(M, x)

  fin = S["finalize"]

  print("\nCheking correctness of the matrix")
  N = 16
  A = S["initial_vector"]
  for i in range(3, N):
    nu1 = FlowerSnark.graph(i).eval(CircuitDoubleCover)
    A = M * A
    nu2 = fin * A
    assert nu1 == nu2
    c_approx = nu1 / 16**i
    print("  ν(J_%i) = %i          (c = %.3f)" % (i, nu1, c_approx))
  print("  Done")

  print("\nCalculating the eigenvector for eigenvalue 16")
  to_dict = lambda g: { b.boundary: b.value for b in g }
  gadget_to_dict = lambda i: to_dict(FlowerSnark.gadget(i).eval_gadget(CircuitDoubleCover))

  X10 = gadget_to_dict(10)
  X9 = gadget_to_dict(9)
  B = [ b for b in X10.keys() if b in X9 and X10[b] / X9[b] > 15 ]
  m = min( X10[b] for b in B )
  vec = { b: int(round(X10[b] / m)) for b in B }
  print("  Guessed eigenvector (with %i non-zero coordinates):" % len(vec))
  for b, v in sorted(vec.items()):
    print("    %i at %s" % (v, b))
  
  print("  Checking it...")
  nex = to_dict(FlowerSnark._next_gadget(FakeGadget(6,
    [ BoundaryValue(b, v) for b, v in vec.items() ])).eval_gadget(CircuitDoubleCover))

  assert nex.keys() == vec.keys()
  for b in vec:
    assert vec[b] * 16 == nex[b], "Oops at %b" % b

  assert set(S["variables"]).issuperset(vec.keys())
  index = { b: i for i, b in enumerate(S["variables"]) }
  vec_ = matrix(QQ, M.ncols(), 1, { (index[b], 0): v for b, v in vec.items() })
  assert M * vec_ == 16 * vec_

  print("  Done")

  print("  Boundaries of with the non-zero value in both the eigenvector and the final vector:")
  for i in range(vec_.nrows()):
    if vec_[i, 0] > 0 and fin[0, i] > 0:
      print("    %s" % (S["variables"][i],))

  print("Check that 16 is the largest eigenvalue")
  print("  Calculating deflation of the matrix")
  D = deflate_eigenvector(M, vec_)
  print("  Bounding the spectral radius by 15")
  bound_spectral_radius(D, 15)