In [1]:
# use magic to install all pylab functions (allows easy definition of n-dimensional empty arrays)
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [23]:
# import libraries
import pulp as lp
import numpy as np
import networkx as nx
from networkx.drawing.nx_pydot import write_dot
In [3]:
# Set params
order = 9
degree = 4
In [4]:
# Assumed params
l = 1
m = 2
In [5]:
# only interested in feasibility, so trivial objective
graphprob = lp.LpProblem(name="SRG({},{},{},{})_problem".format(order,degree,l,m), sense=lp.LpMinimize)
graphprob += 0, "feasibility_obj"
In [6]:
# define variables to indicate if vertices i,j are neighbours
edge_variables = empty((order,order), dtype='object')
for i in range(order):
    for j in range(order):
        name = 'A_{},{}'.format(str(i + 1), str(j + 1))
        var = lp.LpVariable(name, cat='Binary')
        edge_variables[i,j] = var
In [7]:
# define variables to indicate when vertex k is a mutual neighbour of vertices i<j
mutual_variables=empty((order,order,order), dtype='object')
for i in range(order):
    for j in range(order):
        for k in range(order):
            name = 'M_{},{},{}'.format(i + 1, j + 1, k + 1)
            var = lp.LpVariable(name, cat='Binary')
            mutual_variables[i,j,k] = var
In [8]:
# enforce symmetry and non self-adjacency
for i in range(order):
    graphprob += edge_variables[i,i] == 0, "vertex {} not self adjacent".format(i + 1)
    
for i in range(order):
    for j in range(order):
        if i != j:
            graphprob += edge_variables[i,j] == edge_variables[j,i], "{},{} symmetry".format(i + 1, j + 1)
In [9]:
# enforce vertex degree
for i in range(order):
    graphprob += lp.lpSum(edge_variables[i,0:order]) == degree, "vertex {} degree {}".format(i,degree)
In [10]:
# Enforce mutual neighbour conditions (special case of l = m - 1 here!)

for i in range(order):
    for j in range(order):
        if i != j:
            graphprob += (edge_variables[i,j] + lp.lpSum(mutual_variables[i,j,0:order])) == 2, "Neighbour conditions on {},{}".format(i + 1, j + 1)
In [11]:
# Make sure mutual neighbour indicators are reflected in adjacency structure!

# indicator set requires edges
for i in range(order):
    for j in range(order):
        if i != j:
            for k in range(order):
                graphprob+= edge_variables[i,k] >= mutual_variables[i,j,k]
                graphprob+= edge_variables[k,i] >= mutual_variables[i,j,k]
                
                
# mutual neighbouring requires indicator
for i in range(order):
    for j in range(order):
        if i != j:
            for k in range(order):
                graphprob+= edge_variables[i,k] + edge_variables[j,k] - mutual_variables[i,j,k] <= 1
In [12]:
# Symmetry condition on mutual neighbouring

for i in range(order):
    for j in range(order):
            if i != j:
                for k in range(order):
                    graphprob += mutual_variables[i,j,k] == mutual_variables[j,i,k]
In [13]:
len(graphprob.constraints)
Out[13]:
2754
In [14]:
len(graphprob.variables())
Out[14]:
729
In [15]:
# run the solver
%time graphprob.solve()
Wall time: 185 ms
Out[15]:
1
In [16]:
graphprob.status
Out[16]:
1
In [17]:
# populate a matrix in accordance with the edge list found
A=np.zeros((order,order))
for i in range(order):
    for j in range(order):
        if edge_variables[i,j].value()>0:
            A[i,j]=1
In [18]:
# Pieces for verification
Jv=np.ones((order,order))
Iv=np.eye(order)
In [19]:
# verify degree condition
AJ=np.dot(A,Jv)
np.all(AJ==degree)
Out[19]:
True
In [20]:
# verify quadratic condition
cond2=np.dot(A,A)+(m-l)*A+(m-degree)*Iv
np.all(cond2==m)
Out[20]:
True
In [21]:
# take a look at our result
G = nx.from_numpy_matrix(A) 
nx.draw_circular(G, with_labels=False)
C:\Users\Graeme\Anaconda3\lib\site-packages\networkx\drawing\nx_pylab.py:611: MatplotlibDeprecationWarning: isinstance(..., numbers.Number)
  if cb.is_numlike(alpha):
In [24]:
write_dot(G,'Paley9412.dot')
In [ ]: