Linear average time automorphism algorithm for random graphs.
A 1996 Master's paper from the Department of Physics at Pennsylvania State University, reproduced in full. It develops an algorithm that finds the automorphisms of a random graph in linear average time, and with them the Variety — Lee Smolin's measure of how well each point of a structure can be told apart from every other, offered as a quantitative handle on complexity in diffeomorphism-invariant theories such as gravity.
Master’s Paper · 1996
Luis Ferreira · Department of Physics, Pennsylvania State University
Abstract
An algorithm is developed and implemented, whose purpose is to determine the automorphisms of a random graph in linear average time and to calculate a related function of a graph called the Variety. For most part of the probability range the Algorithm is very efficient and so we are able to compute, in a couple of minutes, graphs with up to 1000 vertices. For a short range of probabilities close to zero (or on the symmetric position in relation to 0.5) it might take a much longer time, due, probably, to an increase in symmetry and disconnectivity of the graph. This range's width depends inversely on the number of vertices of the graph.
Problem
There is no know solution for the isomorphism problem for all graphs. It might be a NP-complete problem, meaning not only that it might not be solvable in polynomial time but also that it might belong to a class of problems which if one is proven solvable then all of them are [1, 2]. However there are strategies that can be implemented to solve special cases like the random graph automorphism.
Physical Motivation
The concept of variety, as defined in the following discussion, may be useful for the study of complexity in diffeomorphism invariant theories like gravitational and cosmological theories. It may also provide a quantitative tool to study self-organizing systems [3, 4, 5].
Introduction
A graph is a set of vertices and connections between pairs of vertices. A connection can be in either one of two states: on or off. In a random graph each connection has a probability p of being on (p is between 0 and 1).
The 1st order view of a vertex is the set of vertices to which it is connected.
The nth order view of a vertex is the sub-graph formed by all vertices within n steps from the original vertex and all the connections among them.
The Indifference of a vertex is the number of steps we must travel to obtain a view that is different from all the other views. If we cannot distinguish a vertex from another we say that they are isomorphic and then they have Indifference equal to infinite.
The Indifference of a graph is the sum of the indifference of all the vertices.
We want to calculate a function of a random graph called the Variety.
The Variety measures the difference between the view of each point and all the other points in the graph. It is the normalised sum of the inverses of all the Indifference.
Heuristically this Variety is a measure of "the ability to differentiate each point of the graph from all the others".
(The definitions of this section are in accordance with [3, 4, 5]).
Preliminaries
The program requires some input, namely the probability p that each connection is on and the number of runs (ensemble of random graphs).
The resulting Variety data is written in a file called "var.out" and on the screen as well.
It creates a random graph with probability p and stores it as an Adjacency Matrix.
An Adjacency Matrix is a matrix of zeros and ones whose rows and columns represent the vertices and whose entries represent the connections.
The program accepts a probability range between 0 and 1 but it treats the graphs symmetrically on each side of the value 0.5 (the variety will be the same if we exchange the 0 and 1 entries of the adjacency matrix).
The way the program works is by dividing iteratively the vertices among different baskets which correspond to different views.
Each iteration extracts the n th order view of each vertex (sub-graph X).
Once the program obtains a basket with a single vertex it calculates the variety of that vertex and adds it to the total variety. It proceeds until it has calculated the variety of the entire graph.
To test for isomorphism (same view) of different vertices we adapted an algorithm by Babai and Kucera [6].
The algorithm attributes a Canonical Label to each subgraph.
A Canonical Label is an assignment of a linear order to the set of vertices of a graph.
Two graphs are isomorphic if there exists an order preserving map between their respective Canonical Labels.
Canonical Label Algorithm
(in accordance with [6]).
A single vertex classification procedure (procedure A) has a very high probability of success after just two refinement steps (since we are dealing with random graphs).
An additional depth-first search (procedure B) yields a canonical labelling of all graphs.
V is the set of vertices.
A partial order is a an order of equivalence classes of vertices, C1 < C2 < ... < Cn.
A Canonical Vertex Classification is a partial order P.
Procedure A: Vertex classification procedure.
The input is a sub-graph X.
The first iteration of the classification procedure orders the vertices into equivalence classes according to their valences.
The second iteration refines this order. It partial-orders the elements of each equivalence class by attributing to each a word according to the number of its neighbours in each equivalence class: let x be a vertex; denote the number of neighbours of x in Ci by Ni(x); then x is characterised by a an ordered set of numbers (N1(x), N2(x),..., Nn(x)) which constitutes a word; we partial-order each equivalence class by the lexicographic order of the words corresponding to the vertices in the class.
Let U(X) denote the set of those vertices of X whose equivalence class is not a singleton.
If |U(X)| is bigger than a certain threshold value we jump to procedure C, which we expect to be more efficient at handling those graphs.
Procedure B: Brute force permutation procedure.
The input is a graph X together with a Canonical Vertex Classification P.
P has at most |U(X)|! refinements to a linear order.
To each of these linear orders L, compute the 0-1 array A(L) = (aij) where i belongs to U(X), j belongs to V, aij = 1 or 0 according to weather i and j are adjacent or not, and the rows and columns of the array are arranged in the order defined by L.
Let Lo be any of those L's for which A(Lo) is smaller or equal to A(L) for every L. Output Lo, end.
Any of two such Lo's differ in an automorphism of X, fixing all vertices outside U(X). So we say that Lo is a Canonical Label.
Procedure A is a breadth-first search. Procedure B is a brute force which is efficient if preceded by a very successful breadth-first search.
The following procedures constitute a depth-first search.
Procedure C:
The input is a graph X and a positive integer s.
N is the number of ways to choose s vertices out of the total number of vertices of the graph multiplied by the factorial of those s vertices (number of permutations).
This procedure generates all s-sequences of vertices of X and then applies procedure D to each of them.
If the pair (X, Si) is rejected for every i = 1,..., N, reject the pair (X, s). Else let Ai be the adjacency matrix corresponding to the relative Label with respect to Si, and find the lexicographically first among these matrices.
Output the corresponding Canonical Label.
Procedure D:
The input is a pair (X, S) where S is a sequence of vertices of X.
X/S is the complement of S in X.
To each vertex v belonging to X/S assign a binary word a(v) where aiv =1 or 0 according to weather xi and v are adjacent or not.
If these words are not all different, reject the pair (X, S).
Else order V by xi < xj if i < j; xi < v for all v belonging to X/S; and v < w for v, w belonging to X/S if a(v) < a(w).
Output (X, <).
The Canonical Label Algorithm is then:
1 - Apply A. If A rejects X, go to 3. Else
2 - Apply B to the output of A. Output the result, end.
3 - Set s = 1.
4 - Apply D to the pair (X, s). If accepted, output the result, end. Else
5 - Set s = s + 1 and go to 4.
We store all the Canonical Labels attributed to vertices contained in a basket and compare them, taking care to ensure that the vertex for which we extracted the sub-graph X1 is mapped to the corresponding vertex in the sub-graph X2.
Variety Algorithm
This is the entire algorithm:
1 - Create random graph.
2 - Set Variety = 0, n = 1, basket = V (entire graph).
3 - For each element of the basket extract the number of vertices of the n_th order view (n_valence of that vertex). If n >1 & for all the elements in the basket this n_valence is the same as the (n-1)_valence then go to the next basket.
4 - If n > 1 & last basket then return to next basket in 10 (n = n -1).
5 - Subdivide the basket into smaller baskets containing vertices with the same n_valence and order the resulting baskets. Go to first basket.
6 - If the basket has only a single vertex then Variety = Variety + 1/n and go to the next basket.
7 - If n > 1 & last basket then return to next basket in 10 (n = n - 1) else if n =1 & last basket end.
8 - Extract n_th order view of each vertex belonging to the basket and apply the Canonical Label Algorithm to the sub-graph.
9 - Separate vertices with the same Canonical Label into smaller baskets and order the resulting baskets. Go to first basket.
10 - If the basket has only a single vertex then Variety = Variety + 1/n and go to the next basket.
11 - If n>1 & last basket then go to next basket in 6 else if n =1 & last basket end.
12 - Set n = n +1 and go to 3.
Considerations
The algorithm was implemented in C and it required the standard math and input/output libraries.
The decision to apply either procedure B or procedure C is made according to the fraction of |U(X)| which belongs to a single class and the size of this class.
The relevant time consuming part of the algorithm resides in the calculation of the permutations of the vertex labels.
Certain strategies were followed in the implementation in order to make the algorithm more efficient. The most relevant of those is the following:
To avoid calculating all the permutations we may consider only those that decrease the position in the lexicographic order. Since the calculation of the permutations is itself ordered, this is easily accomplished.
Analysis
For most part of the probability range the Algorithm is very efficient, since the running time is dominated by Procedure A, and so we are able to compute, in a couple of minutes, graphs with up to 1000 vertices.
For a short range of probabilities close to zero (or on the symmetric position in relation to 0.5) it might take a much longer time, due, probably, to an increase in symmetry and disconnectivity of the graph.
This range's width depends inversely on the number of vertices of the graph.
References
- I. Ponomarenko, "Graph Algebras and the Graph Isomorphism Problem", AAECC 5, 277-286 (1994). https://doi.org/10.1007/BF01225642
- D. Corneil & R. Read, "The Graph Isomorphism Disease", Journal of Graph Theory, Vol. 1 (1977) 339-363. https://onlinelibrary.wiley.com/doi/abs/10.1002/jgt.3190010410
- L. Smolin, "Variety, complexity and cosmology", Preprint. https://arxiv.org/abs/hep-th/9203041
- L. Smolin, "Cosmology as a problem in critical phenomena", gr-qc/9505022 in "Complex systems and binary networks eds. R. Lopez-Pena, R. Capovilla, R. Garcia-Pelayo, H. Waelbroeck and F. Zertuche (Springer-Verlag, Berlin, 1995). https://arxiv.org/abs/gr-qc/9505022
- L. Smolin, "Space and Time in the Quantum Universe" in "Conceptual Problems of Quantum Gravity" ed. by A. Ashtekar and J. Stachel, (Birkhauser, Boston, 1991). https://philpapers.org/rec/SMOSAT
- L. Babai & L. Kucera, "Canonical Labelling of Graphs in linear average time," IEEE 2/79, 39-46. https://ieeexplore.ieee.org/document/4567999
Program
The implementation is a single C89 translation unit. It is reproduced below as it was written, with only three changes needed to make it compile under a modern toolchain: a missing #include <stdlib.h>, an explicit int return type on main(), and a corrected pointer dereference in the canonical-label routine.
Show the full C source (variety.c, ~1000 lines)
/* VARIETY */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define NMAX 1000
void make_graph();
void type_graph();
void quick_sort(int *b[], int left, int right);
void swap(int *b[], int i, int j);
int sub_mat(int vert);
int sub_mat_rep(int i, int j, int flag3);
void add_val(int vert);
void k_val(int ite, int lef, int rig);
void first_val();
void can_label(int ver,int iter, int pos);
void label(int nodes);
void lexico(int lef, int rig, int ki);
int divis(int i, int j, float mod);
void permut(int ver, int i, int j);
void lexico_can(int lef, int rig);
void sub_fval(int ver);
void procedure_d(int ver);
void s_sequence(int k, int level, int ver);
void permutation(int i, int j);
void procedure_c();
void same(int lef, int rig, int ki);
int a[NMAX][NMAX], val[NMAX], fval[NMAX], firval[NMAX], smat[NMAX][NMAX],
fsmat[NMAX][NMAX], firsmat[NMAX][NMAX];
int *b[NMAX], *aux[NMAX], lote[NMAX][3];
float p;
/* double p; */
double conv, variety;
int i, j, k;
int iseed, indifference; /* iseed=random seed */
char xs;
int largo, largo2, largomajor;
int *c[NMAX], ucol[NMAX], ulin[NMAX], *d[NMAX], *e[NMAX], *e2[NMAX], **f[NMAX],
*swp, *swdp, order[NMAX], ord[NMAX], sa[NMAX], sb[NMAX], si[NMAX], cs[NMAX],
*h[NMAX], nodes, u, v, class, canclass, s, o, jup, vertices, klass;
int *labp, *classp, *canclassp, *permutp, *labelp, *wordp, *canlabelp;
int clip[NMAX], *clp[NMAX];
int count, first;
int inic;
int outofhere;
/*------------------------------------------------------------------*/
int main()
{
int iter, lef, rig, one, r, runs;
FILE *file, *fopen();
file=fopen("var.out","a");
iseed=55577;
one=0;
printf("\n\nNumber of vertices = %d\n\n", NMAX);
fprintf(file, "\nNumber of vertices = %d\n", NMAX);
printf("\nThe data will be written in a file called <var.out>\npress return to continue");
xs=getchar();
printf("\nprobability? (number between 0 and 1; if you want to exit just write 0)\n");
scanf("%f", &p);
printf("probability = %f\n", p);
while(p!=0){
printf("\n number of trials? ");
scanf("%d", &runs);
fprintf(file, "\nProbability = %f , number of trials = %d\n\nVariety:\n\n", p, runs);
for(r=0; r<runs; r++){
largomajor=0;
indifference=0;
variety=.0;
conv=.0;
lef=0;
rig=NMAX-1;
printf("\ninitializing variables...\n");
for(i=0;i<NMAX;i++)
a[i][i]=0;
for(i=0;i<NMAX;i++)
b[i]=&val[i];
for(i=0; i<NMAX; i++){
for(j=1; j<NMAX; j++)
smat[i][j]=-1;
lote[i][2]=1;
lote[i][0]=lote[i][1]=0;
smat[i][0]=i;
}
printf("done\n");
srand(iseed);
printf("\ncalculating the adjacency matrix...\n");
make_graph();
printf("done\n");
if(p>.5){ /* transforms adjacency matrix into anti-adjacency matrix */
for(i=0; i<NMAX; i++)
for(j=0; j<NMAX; j++)
a[i][j]=1-a[i][j];
for(i=0;i<NMAX;i++)
a[i][i]=0;
}
first_val();
printf("\n\nnumber of vertices completed:\n");
k_val(1,lef,rig);
variety/=NMAX;
printf("\n\n VARIETY =%f indifference =%3d largomajor =%2d", variety,indifference,largomajor);
printf("\nfinished trial%3d\n", r+1);
fprintf(file, "%f\n", variety);
}
printf("\nprobability? (number between 0 and 1; if you want to stop just write 0) ");
scanf("%f", &p);
}
printf("\nTHE END\n\n");
}
/*------------------------------------------------------------------*/
/* makes a graph with NMAX nodes and prob p */
void make_graph()
{
int r;
double size;
size=32768;
iseed=rand();
srand(iseed);
for(j=0;j<NMAX;j++){
for(k=0; k<j; k++){
r=rand();
if((r/size) < p )
a[j][k]=1;
else
a[j][k]=0;
a[k][j]=a[j][k];
}
}
}
/*------------------------------------------------------------------*/
/* prints graph */
void type_graph()
{
printf("\n N = %d , p = %f , iseed = %d \n\n ",NMAX,p,iseed);
for(i=0; i<NMAX; i++)
printf("%2d", i);
printf("\n\n");
for(j=0;j<NMAX;j++)
{
printf("%2d ", j);
for(k=0;k<NMAX;k++)
{
printf("%2d",a[j][k]);
}
printf("\n");
}
printf("\n");
}
/*------------------------------------------------------------------*/
/* creates matrix with first_valences */
void first_val()
{
for(i=0; i<NMAX; i++)
firval[i]=0;
for(i=0; i<NMAX; i++)
for(j=0; j<NMAX; j++)
firsmat[i][j]=-1;
for(i=0; i<NMAX; i++){
/* firsmat[i][0]=i; */
k=0;
for(j=0; j<NMAX; j++)
if(a[i][j]==1){
firval[i]+=1;
firsmat[i][k]=j;
k++;
}
}
}
/*------------------------------------------------------------------*/
/* sorts val[left]...val[right] into increasing order */
void quick_sort(int *b[], int left, int right)
{
int ii, last;
if (left >= right)
return;
swap(b, left, (left+right)/2);
last=left;
for (ii=left+1; ii<=right; ii++)
if (*b[ii]<*b[left])
swap(b, ++last, ii);
swap(b, left, last);
quick_sort(b, left, last-1);
quick_sort(b, last+1, right);
}
/*------------------------------------------------------------------*/
/* swap: interchange v[i] and v[j] */
void swap(int *b[], int ii, int jj)
{
int *temp;
temp=b[ii];
b[ii]=b[jj];
b[jj]=temp;
}
/*------------------------------------------------------------------*/
/* extracts the sub-matrix with the k_valences */
int sub_mat(int vert)
{
int flg1,flg2,flag1,flag2,flag3,jmp;
vert=b[vert]-val;
jmp=0;
flg2=lote[vert][2];
flg1=lote[vert][1];
flag1=lote[vert][0];
k=flag2=flag3=flg2;
for(i=flag1; i<flg2; i++)
*aux[i]=smat[vert][i];
while(flg1<flg2){
i=smat[vert][flg1];
for(j=flag2; j<flag3; j++)
*aux[j]=smat[vert][j];
quick_sort(aux, flag1, flag3-1);
inic=flag1;
for(j=0; firsmat[i][j]!=-1; j++)
if(sub_mat_rep(i,j,flag3)==0){
smat[vert][k]=firsmat[i][j];
k++;
}
flg1++;
flag2=flag3;
flag3=k;
}
lote[vert][0]=lote[vert][1];
lote[vert][1]=lote[vert][2];
lote[vert][2]=k;
if(k==flg1)
jmp=1;
return jmp;
}
/*------------------------------------------------------------------*/
/* avoids repeating valences */
int sub_mat_rep(int i, int j, int flag3)
{
int rep, ii;
if(firsmat[i][j]>*aux[flag3-1])
return 0;
else if(firsmat[i][j]==*aux[flag3-1])
return 1;
else{
for(ii=inic; ii<flag3; ii++)
if(firsmat[i][j]<=*aux[ii]){
inic=ii;
break;
}
if(firsmat[i][j]==*aux[ii])
rep=1;
else rep=0;
return rep;
}
}
/*------------------------------------------------------------------*/
/* sums the nodes from the k_valences */
void add_val(int vert)
{
vert=b[vert]-val;
for(i=1; smat[vert][i]!=-1 && i<NMAX; i++);
val[vert]=i-1;
}
/*------------------------------------------------------------------*/
/* recursively calculates the k_valences, orders the vertices and calculates the variety */
void k_val(int ite, int lef, int rig)
{
int flg1, flg2, jump, vertex, fg1, fg2, ver;
jump=0;
printf("\n(calculating sub_matrix; iteration # %d\n",ite);
for(vertex=lef; vertex<=rig; vertex++){
for(i=0; i<NMAX; i++)
aux[i]=&smat[b[vertex]-val][i];
jump=sub_mat(vertex);
add_val(vertex);
}
printf(" ... done)");
quick_sort(b,lef,rig);
flg1=flg2=lef;
while(flg2<=rig){
printf("\n %d",flg2);
while(flg2<rig && *b[flg2]==*b[flg2+1])
flg2++;
if(flg1==flg2){
indifference+=ite;
conv=ite;
variety+=1.0/conv;
jump=0;
flg2++;
flg1=flg2;
}
else if(jump==1){ /* eliminates the zero valences */
jump=0;
flg2++;
flg1=flg2;
}
else{
nodes=val[b[flg1]-val]+1;
vertices=flg2-flg1+1;
canlabelp=(int *)calloc(nodes*nodes*vertices, sizeof(int));
for(ver=flg1; ver<=flg2; ver++){
labelp=(int *)calloc(nodes, sizeof(int));
classp=(int *)calloc(nodes, sizeof(int));
can_label(b[ver]-val,ite, ver); /* stores a canonical label in labelp */
/* convert canlabel into number -- DEAD CODE, commented out
p=0;
q=r=1;
for(j=nodes-1; j>=0; j--)
for(k=nodes-1; k>=0; k--){
if(q>nodes*nodes) goto out;
else if(r<128) canlabel+=pow(2,r)*a[*(labelp+j)][*(labelp+k)];
else{
r=0;
*(canlabelp+(ver-flg1)*(nodes*nodes)/127+p)=canlabel;
canlabel=0;
p++;
}
q++;
r++;
}
out:;
*/
for(j=0; j<nodes; j++)
for(k=0; k<nodes; k++)
*(canlabelp+(ver-flg1)*nodes*nodes+j*nodes+k)=a[*(labelp+j)][*(labelp+k)];
free(classp);
free(labelp);
}
canclass=0;
k=0;
canclassp=(int *)calloc(vertices, sizeof(int));
lexico_can(0,vertices-1);
free(canlabelp);
for(i=flg1; i<=flg2; i++)
*b[i]=*(canclassp+i-flg1);
free(canclassp);
quick_sort(b, flg1, flg2);
fg1=fg2=flg1;
while(fg2<=flg2){
while(fg2<flg2 && *b[fg2]==*b[fg2+1])
fg2++;
if(fg1==fg2){
indifference+=ite;
conv=ite;
variety+=1.0/conv;
fg2+=1;
fg1=fg2;
}
else{
k_val(ite+1,fg1,fg2);
fg2+=1;
fg1=fg2;
}
}
}
flg2+=1;
flg1=flg2;
}
return;
}
/*------------------------------------------------------------------*/
/*------------------------------------------------------------------*/
/* atributes a canonical label to a graph */
void can_label(int ver,int iter, int pos)
{
int lg;
sub_fval(ver);
labp=(int *)calloc(nodes*(nodes+1), sizeof(int));
label(nodes);
free(labp);
if(u<2){
for(i=0; i<nodes; i++)
*(labelp+i)=smat[ver][**f[i]];
}
else{
lg=log(nodes)/log(log(nodes));
if(largo>9 && largo>u*2.0/3.0){ /* (u > lg) 7! = 5040 , 8! = 40320 */
printf("\n working on %d th vertex , largo = %d , u = %d", pos, largo, u);
printf("\n TROUBLE: (10! = 3628800 trials) If our breadth-first procedure is not efficient enough there is a real possibility that this might take quite a long while because we might have to calculate all %d! permutations\n", u);
procedure_d(ver);
}
else{
if(largo > 8 || (largo>7 && largo2>1))
printf("\n\n This might take a while since we have to calculate something of the order of %d! * %d! (9! = 362880) , u = %d\n", largo, largo2, u);
permutp=(int *)calloc(u*nodes, sizeof(int));
swp=(int *)calloc((u-1)*2, sizeof(int));
for(i=0; i<u; i++)
for(j=0; j<nodes; j++)
*(permutp+i*nodes+j)=a[smat[ver][*d[i]]][smat[ver][**f[j]]];
for(i=0; i<nodes; i++)
*(labelp+i)=smat[ver][**f[i]];
permut(ver,0,0);
free(swp);
free(permutp);
}
}
}
/*------------------------------------------------------------------*/
/* labelling a graph */
void label(int nodes)
{
float mod;
mod=4.0;
for(i=0; i<nodes; i++)
for(j=0; j<nodes; j++)
*(labp+i*nodes+j)=0;
for(i=0; i<nodes; i++){
clp[i]=(labp+i*nodes);
*(labp+i*nodes)=fval[i];
}
quick_sort(clp,0,nodes-1);
klass=1;
clip[0]=klass;
for(i=0; i<nodes-1; i++){
if(*clp[i]!=*clp[i+1]) klass++;
clip[i+1]=klass;
}
for(i=0; i<nodes; i++)
*clp[i]=clip[i];
for(i=0; i<nodes; i++)
for(j=0; fsmat[i][j]!=-1; j++)
(*(labp+i*nodes+*(labp+fsmat[i][j]*nodes)))++;
/* we could do mod4 */
/* for(i=0; i<nodes; i++)
for(j=1; j<=klass; j++)
*(labp+i*nodes+j)=divis(i,j,mod);
*/
/* convert label into number -- DEAD CODE, commented out
p=q=0;
while(q<nodes){
for(i=128*2/mod-1; i>=0; i--){
if(q<nodes) label+=pow(mod,i)*(*(labip+q));
else break;
q++;
}
*(labp+p)=label;
p++;
}
*/
largo2=largo=0;
u=v=0;
class=0;
lexico(0, nodes-1, 0);
}
/*------------------------------------------------------------------*/
/* lexicographic order */
void lexico(int lef, int rig, int ki)
{
int flg1,flg2; /* *c[nodes], ucol[nodes], ulin[nodes] */
if(ki==0)
for(i=0; i<nodes; i++){
c[i]=(labp+i*nodes);
ulin[i]=-1;
ucol[i]=-1;
*(classp+i)=-1;
}
else
for(i=lef; i<=rig; i++)
c[i]++;
quick_sort(c,lef,rig);
flg1=flg2=lef;
while(flg2<=rig){
while(flg2<rig && *c[flg2]==*c[flg2+1])
flg2++;
if(flg1!=flg2 && ki<klass) lexico(flg1,flg2, ki+1);
else{
for(i=flg1; i<=flg2; i++){
if(flg2-flg1+1>largo){
largo2=largo;
largo=(flg2-flg1+1);
}
if(flg2-flg1+1>largomajor)
largomajor=(flg2-flg1+1);
ucol[i]=(c[i]-labp)/nodes;
*(classp+ucol[i])=class;
if(flg1!=flg2){
ulin[u]=ucol[i];
d[u]=&ulin[u];
f[i]=&e[u];
e[u]=&ucol[i];
u++;
}
else{
f[i]=&e2[v];
e2[v]=&ucol[i];
v++;
}
}
class++;
}
flg2+=1;
flg1=flg2;
}
}
/*------------------------------------------------------------------*/
/* mod function */
int divis(int i, int j, float mod)
{
int x;
x=*(labp+i*nodes+j)/mod;
x=(*(labp+i*nodes+j)/mod-x)*mod;
return x;
}
/*------------------------------------------------------------------*/
/* permutations: returns canonical label in *f[nodes] */
void permut(int ver, int i, int j)
{
int equal;
while(j<u-1){
while(j<u && *(classp+*d[i])==*(classp+*d[j])){
*(swp+i)=i;
*(swp+u-1+i)=j;
permut(ver,i+1,i+1);
j++;
}
return;
}
for(i=0; i<u-1; i++){
swap(d,*(swp+i),*(swp+u-1+i));
swap(e,*(swp+i),*(swp+u-1+i));
}
equal=1;
for(i=0; i<u; i++){
for(j=0; j<nodes; j++){
if(equal==1){
if(a[smat[ver][*d[i]]][smat[ver][**f[j]]]>*(permutp+i*nodes+j)) goto out;
else if(a[smat[ver][*d[i]]][smat[ver][**f[j]]]<*(permutp+i*nodes+j)) equal=0;
}
*(permutp+i*nodes+j)=a[smat[ver][*d[i]]][smat[ver][**f[j]]];
}
}
for(i=0; i<nodes; i++)
*(labelp+i)=smat[ver][**f[i]];
out:;
for(i=u-2; i>=0; i--){
swap(d,*(swp+i),*(swp+u-1+i));
swap(e,*(swp+i),*(swp+u-1+i));
}
return;
}
/*------------------------------------------------------------------*/
/* lexicographic order of canonical labels */
void lexico_can(int lef, int rig)
{
int flg1,flg2;
if(k==0)
for(i=0; i<vertices; i++)
c[i]=(canlabelp+i*nodes*nodes);
else
for(i=lef; i<=rig; i++)
c[i]+=1;
k++;
quick_sort(c,lef,rig);
flg1=flg2=lef;
while(flg2<=rig){
while(flg2<rig && *c[flg2]==*c[flg2+1])
flg2++;
if(flg1!=flg2 && k<nodes*nodes) lexico_can(flg1,flg2);
else{
for(i=flg1; i<=flg2; i++)
*(canclassp+(c[i]-canlabelp)/(nodes*nodes))=canclass;
canclass++;
}
flg2+=1;
flg1=flg2;
}
return;
}
/*------------------------------------------------------------------*/
/* creates matrix with first_valences of sub_graph */
void sub_fval(int ver)
{
for(i=0; i<nodes; i++)
fval[i]=0;
for(i=0; i<nodes; i++)
for(j=0; j<nodes; j++)
fsmat[i][j]=-1;
for(i=0; i<nodes; i++){
k=0;
for(j=0; fsmat[ver][j]!=-1 && j<nodes; j++)
if(a[smat[ver][i]][smat[ver][j]]==1){
fval[i]+=1;
fsmat[i][k]=j;
k++;
}
}
}
/*------------------------------------------------------------------*/
/*------------------------------------------------------------------*/
/* tests the difficult graphs */
/* outputs canonical label in *labelp[nodes] */
void procedure_d(int ver)
{
/* s=6*sqrt(nodes); */
/* s=sqrt(nodes); */
s=1;
o=0;
while(s<nodes){
count=0;
swdp=(int *)calloc(s, sizeof(int));
s_sequence(0,0,ver);
free(swdp);
if(o!=0)
break;
else
s++;
}
}
/*------------------------------------------------------------------*/
/* extracts s elements from smat[nodes] */
void s_sequence(int k, int level, int ver)
{
int ii, kk;
while(level<s){
for(ii=k; ii<=(level+nodes-s); ii++){
sa[level]=smat[ver][ii];
si[level+1]=ii;
s_sequence(ii+1, level+1, ver);
}
return;
}
si[0]=-1;
si[level+1]=nodes;
kk=0;
for(i=0; i<s+1; i++) /* vertices not in the s-sequence */
if(si[i]+1!=si[i+1])
for(j=si[i]+1; j<si[i+1]; j++){
cs[kk]=smat[ver][j];
kk++;
}
first=1;
permutation(0,0);
return;
}
/*------------------------------------------------------------------*/
/* permutes the s elements from smat[nodes] */
void permutation(int i, int j)
{
int temp, equal, ii, jj;
jup=0;
while(j<s-1){
while(j<s){
*(swdp+i)=i;
*(swdp+s-1+i)=j;
permutation(i+1, i+1);
j++;
}
return;
}
for(ii=0; ii<s; ii++)
sb[ii]=sa[ii];
for(ii=0; ii<s-1; ii++){
temp=sb[*(swdp+ii)]; /* swap(i,j) */
sb[*(swdp+ii)]=sb[*(swdp+s-1+ii)];
sb[*(swdp+s-1+ii)]=temp;
}
count++;
wordp=(int *)calloc((nodes-s)*s, sizeof(int));
procedure_c();
free(wordp);
if(jup==0){
if(first==1){
for(i=0; i<s; i++)
*(labelp+i)=sb[i];
for(i=s; i<nodes; i++)
*(labelp+i)=ord[i-s];
first=0;
}
equal=1;
for(j=0; j<nodes; j++)
for(k=0; k<nodes; k++) /* canonical label in *labelp */
if(k<s){
if(equal==1){
if(a[sb[j]][sb[k]]>a[*(labelp+j)][*(labelp+k)]) goto out;
else if (a[sb[j]][sb[k]]<a[*(labelp+j)][*(labelp+k)]) equal=0;
}
*(labelp+k)=sb[k];
}
else{
if(equal==1){
if(a[ord[j-s]][ord[k-s]]>a[*(labelp+j)][*(labelp+k)]) goto out;
else if (a[ord[j-s]][ord[k-s]]<a[*(labelp+j)][*(labelp+k)]) equal=0;
}
*(labelp+k)=ord[k-s];
}
out:;
o++;
}
else jup=0;
return;
}
/*------------------------------------------------------------------*/
/* word[i][j] */
void procedure_c() /* wordp */
{
for(i=0; i<nodes-s; i++){
for(j=0; j<s; j++){
*(wordp+i*s+j)=a[cs[i]][sb[j]];
}
}
outofhere=0;
same(0, nodes-s-1, 0);
}
/*------------------------------------------------------------------*/
/* checks if words are all different and orders them */
void same(int lef, int rig, int ki) /* ord */
{
int flg1,flg2; /* *h[nodes-s] */
if(outofhere==1)
return;
if(ki==0)
for(i=0; i<nodes-s; i++)
h[i]=(wordp+i*s);
else
for(i=lef; i<=rig; i++)
h[i]+=1;
quick_sort(h,lef,rig);
flg1=flg2=lef;
while(flg2<=rig){
while(flg2<rig && *h[flg2]==*h[flg2+1])
flg2++;
if(flg1!=flg2 && ki<s-1)
same(flg1,flg2,ki+1);
else{
for(i=flg1; i<=flg2; i++){
if(flg1!=flg2){
jup=1;
outofhere=1;
goto same_word; /* tenho que sair destes loops */
}
else{
ord[i]=cs[(h[i]-wordp)/s];
}
}
}
flg2+=1;
flg1=flg2;
}
same_word:;
}
/*------------------------------------------------------------------*/
/* THE END */Afterword
Thirty years later · 2026
The paper above was written in 1996. The program it describes is about a thousand lines of C, and roughly nine hundred of them exist to answer a single question: are these two subgraphs the same shape?
That question now has a library answer. What follows is the same computation, from the same definitions, in thirty-five lines of Python.
What the thousand lines were for
Procedures A, B, C and D are a canonical-labelling algorithm. They take a subgraph and impose a linear order on its vertices such that two subgraphs receive order-preserving labels exactly when they are isomorphic. Procedure A classifies vertices by valence and then by neighbour-class words; Procedure B brute-forces the permutations that A could not separate; C and D fall back to a depth-first search over vertex sequences when the graph is too symmetric for either.
None of that is about Variety. It is scaffolding, and every line of it was necessary in 1996 because there was no other way to get the answer.
Two library calls replace all of it. A Weisfeiler-Leman hash sorts the views into buckets cheaply, and VF2++ settles the rare collisions exactly. The Variety algorithm itself — the basket refinement, the n_valence saturation test, the accumulation of one over n — survives essentially unchanged, because that part was never the hard part.
Procedure A was ahead of its time.
The breadth-first vertex classification of Procedure A is the Weisfeiler-Leman algorithm, arrived at independently. In 2019 it was shown that the expressive power of every standard graph neural network is bounded above by exactly this procedure: if 1-WL cannot tell two graphs apart, no such network can either. A subroutine written here to make a 1996 thesis tractable turned out to be a fundamental limit on a class of machine learning models that did not yet exist.
The program
#!/usr/bin/env python3
"""Variety of a random graph — the 1996 algorithm, in 2026 Python.
This is the modern counterpart of variety.c. It computes the same quantity from
the same definitions:
nth order view of v the subgraph induced by every vertex within n steps
of v, together with all connections among them,
considered as a graph rooted at v
indifference of v the least n whose view is unique among all vertices,
or infinity if no such n exists
Variety (1/N) * sum over v of 1/indifference(v)
The 1996 program spends about a thousand lines on Procedures A-D, which exist
only to decide whether two views are isomorphic. That question now has a library
answer, so the whole canonical-labelling apparatus collapses into two calls:
a Weisfeiler-Leman hash to bucket views cheaply, and VF2++ to confirm exactly
within a bucket. What is left is the Variety algorithm itself.
pip install networkx
python3 variety.py --n 200 --p 0.05 --trials 3
Luis Ferreira, 1996 (C) / 2026 (Python).
"""
import argparse
from collections import defaultdict
import networkx as nx
def rooted_view(G, v, n):
"""The nth order view of v: the n-ball, induced, with v marked as the root.
Marking the root is what the thesis calls "taking care to ensure that the
vertex for which we extracted the sub-graph X1 is mapped to the
corresponding vertex in the sub-graph X2".
"""
ball = nx.single_source_shortest_path_length(G, v, cutoff=n)
H = G.subgraph(ball).copy()
nx.set_node_attributes(H, {u: int(u == v) for u in H}, 'root')
return H
def variety(G):
"""Return (Variety, {vertex: indifference}). Vertices absent from the dict
are the isomorphic ones, of infinite indifference, contributing nothing."""
N = G.number_of_nodes()
indifference = {}
previous_size = dict.fromkeys(G, -1)
for n in range(1, N + 1):
pending = [v for v in G if v not in indifference]
if not pending:
break
views = {v: rooted_view(G, v, n) for v in pending}
size = {v: views[v].number_of_nodes() for v in pending}
# The thesis's n_valence test: once no view is still growing, no deeper
# view can separate anything, and whoever is left is indifferent.
if all(size[v] == previous_size[v] for v in pending):
break
previous_size = size
# Bucket by a cheap certificate; only views that collide need VF2++.
buckets = defaultdict(list)
for v in pending:
buckets[nx.weisfeiler_lehman_graph_hash(views[v], node_attr='root')].append(v)
for group in buckets.values():
if len(group) == 1:
indifference[group[0]] = n
continue
# Sort the colliding views into isomorphism classes, comparing each
# against one representative per class rather than against every
# other view. In the symmetric graphs near p = 0 and p = 1 -- the
# range the thesis warns about -- there is only one class, so this
# is the difference between O(g) and O(g*g) exact tests.
classes = []
for v in group:
for cls in classes:
if nx.vf2pp_is_isomorphic(views[v], views[cls[0]],
node_label='root'):
cls[1] += 1
break
else:
classes.append([v, 1])
for representative, count in classes:
if count == 1:
indifference[representative] = n
return sum(1.0 / n for n in indifference.values()) / N, indifference
def main():
ap = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument('--n', type=int, default=1000, help='number of vertices')
ap.add_argument('--p', type=float, required=True, help='edge probability')
ap.add_argument('--trials', type=int, default=1, help='ensemble size')
ap.add_argument('--seed', type=int, default=55577, help='the 1996 seed')
args = ap.parse_args()
print(f'N = {args.n}, p = {args.p}, trials = {args.trials}')
for t in range(args.trials):
G = nx.erdos_renyi_graph(args.n, args.p, seed=args.seed + t)
v, ind = variety(G)
print(f' trial {t + 1}: Variety = {v:.6f}'
f' ({len(ind)}/{args.n} vertices distinguished)')
if __name__ == '__main__':
main()
What changed, and what did not
The 1996 program handled a thousand vertices in a couple of minutes on roughly 150 MHz hardware. This one does the same graph in about a second, and most of that gain is the machine rather than the method.
What genuinely changed is the failure mode. The old program had to degrade to Procedures C and D when a graph was too symmetric for the breadth-first pass, and that was where the running time could become unpleasant. The modern version pushes the same difficulty into VF2++, which is a far better search than anything that could reasonably have been written by hand into a master's thesis.
What did not change is the definition, or the shape of the algorithm, or the fact that the interesting behaviour lives in a narrow band of probabilities near the ends of the range, where graphs are symmetric and disconnected and every method available is slow.
A note on the definitions
One sentence in the Introduction does not survive re-reading. The Indifference of a graph is defined there as the sum of the indifference of all its vertices. But a vertex that cannot be distinguished from any other has infinite indifference, and a single such vertex makes the whole sum infinite. For every graph sparse enough or dense enough to carry a symmetry, the quantity is simply infinite, and it cannot rank anything against anything.
The program never computes it as defined. In the C listing above, the running total is incremented at exactly two places, and both sit inside the branch that fires only once a basket has been reduced to a single vertex — that is, only for vertices that have just been distinguished. The infinite ones are never added, because they are never reached. What the variable actually holds is the sum over the distinguishable vertices alone, and it is printed to the screen and then discarded: never written to the output file, never used in a calculation. The definition and the implementation are two different quantities, and the paper uses neither.
Variety has no such problem, but for a reason the paper leaves unsaid. It is the normalised sum of the inverses of the indifferences, and the inverse of infinity is zero, so an indistinguishable vertex contributes nothing instead of poisoning the total. That convention is the load-bearing assumption of the entire definition, and it appears nowhere in the text. It should have been stated.
The discarded sentence is not worthless, though. Read as a predicate rather than as a measure, it says something exact:
The Indifference of a graph is finite precisely when every vertex is distinguishable from every other, which is precisely when the graph has no non-trivial automorphism.
That is a sharper statement than Variety provides, and it is a yes-or-no question with a clean answer. Measured, all vertices become distinguishable once the mean number of connections per vertex exceeds the natural logarithm of the number of vertices: at a hundred vertices the threshold falls at five connections, at three hundred at six, against logarithms of 4.6 and 5.7. Below it, isolated and duplicated vertices survive and remain interchangeable. Above it, the graph is asymmetric and every point has an identity of its own.
The range of probabilities over which all points are distinguishable is therefore bounded below by roughly log N over N, and it widens as the graph grows. It is also exactly symmetric, which the graded version is not. Exchanging every connection for its absence — which is what moving from p to 1 minus p amounts to — leaves the automorphism group untouched, so the upper threshold sits the same distance below one. The paper claims this symmetry for Variety, where it does not quite hold, because indifference counts steps and complementation does not preserve distances. For the question of whether the points can be told apart at all, the symmetry is exact.
Where Variety stops measuring
The paper never plotted its own quantity. Doing so now, beside a measure of how complex the graph actually is, shows something the 1996 text does not discuss.
The complexity measure is the one the paper's own reasoning appeals to. There is exactly one graph with no connections and exactly one with all of them, while the number of labelled graphs at intermediate densities is enormous. The count of graphs with m edges out of M possible is M choose m, and its logarithm, per edge slot, is the binary entropy h(q). It vanishes at both ends and peaks at one half, which is the shape the intuition expects.
Variety does not have that shape. It rises, reaches one, and stays there.
Two things follow. The first is that Variety is bounded from above. Indifference is at least 1, so each term 1/indifference is at most 1, and the normalised sum cannot exceed 1. Once a graph is dense enough that every vertex is unique after a single step, Variety is pinned, and a graph at p = 0.05 and a graph at p = 0.5 receive the same score.
At a thousand vertices that ceiling arrives at about p = 0.04, where the complexity measure has covered less than a quarter of its range and is still climbing. Everything the entropy does between there and p = 0.5, Variety reports as identical.
The second is that Variety is a staircase rather than a curve. It is an average of reciprocals of small integers, so it moves in steps as vertices migrate from indifference three to two to one. The shelf at exactly one half, visible at a thousand vertices between p = 0.006 and p = 0.012, is the regime in which a sparse vertex's immediate neighbourhood is very nearly a star and so carries only its degree, which is too little to be unique among a thousand, while the second step reveals the neighbours' degrees and is more than enough.
None of this makes Variety wrong. It measures what it claims to measure, the ability to differentiate each point of the graph from all the others. But that ability saturates, and it saturates early, at a density of order log N over N. A quantity intended to track how much structure a universe contains will report that almost every graph is equally maximal.
Source, figures and a modern WL-based reimplementation: github.com/ferrelm/variety-graph
Comentários
Enviar um comentário