/* Reads graphs in g6 code or multicode (optional) from stdin and directs them
options:
ix means: the indegree of every vertex may be at most x.
oy means: the outdegree of every vertex may be at most y.
S means: allow that for every pair of vertices x,y at most one of the edges x-->y and y-->x may be present. By default both of them may be present in the same graph.
T means: Output directed graphs in T-code. This is a simple ASCII output format. Every line contains one graph. First the number of vertices, then the number of directed edges and then the list of directed edges with the start first and the end then. E.g.: 3 2 0 1 2 1 means 3 vertices, 2 directed edges: 0-->1 and 2-->1
B means: Output the directed graphs in a binary code. Every item of the code is an unsigned char. The first unsigned char is the number nv of vertices. The vertices are numbered 1..nv Then the list of vertices x for which there is a directed edge 1->x follow. This list is ended by a 0. Then the list of outgoing neighbours of 2 follows -- again ended with a 0, etc. The code is complete with the 0 ending the list of outgoing neighbours of nv.
Z means: Output the directed graphs in digraph6 code. See formats.txt for a complete definition.
C means: Do really construct all the directed graphs in memory, but don't output them. This is not a big difference in case of restricted in- and outdegrees, because all that is done extra is that edges are directed instead of just keeping track of in- and out-degrees. This option is intended only for testing purposes to test also routines that are normally not used when counting. Things that would speed up the counting also in some cases of restricted in- and out-degrees -- like multiplying the possibilities of assigning directions to edges that can be assigned directions independent of each other (depending on the degrees of the endvertices and overlaps) -- are not included. In case of not restrictive bounds on the in- and out-degree it not really constructing the graphs can be considerably faster. In cases of restricted in- and out-degrees the only difference is that the graph isn't modified... The fact that in case of no output the graph is not modified is mainly to save time for the one case of waterclusters, where large numbers were determined. If large numbers (without output) for other cases shall be determined, one should think about adding the multiplication routines.
m read multicode
This program uses different labelling routines -- all based on the ideas of
G. Brinkmann, Generating water clusters and other directed graphs, Journal of Mathematical Chemistry 46, 1112--1121 (2009)
October 10, 2011: corrected error caused by overflow of 32bit int used as hashvalue.
Sep, 2012: PROCESS feature added by BDM.
Oct, 2017: digraph6 output added by BDM.
*/
/* PROCESS feature * * If PROCESS is defined, it must expand as the name of a procedure * with prototype like void PROCESS(FILE *f, graph *g, int n). This * procedure will be called for each output graph before it is * output (or before output is suppressed) with f being the output * file (possibly NULL). * It is an error if n > WORDSIZE. * If NOCONVERSE is also defined, the call is only made for one of * each digraph and its converse.
*/
/* SUMMARY feature * * If SUMMARY is defined, it must expand as the name of a procedure * with prototype void SUMMARY(void). It is called at the end after * the normal summary.
*/
BOOG edgelist[MAX_BOGEN+1]; /* de lijst van bogen */
BOOG edgelist_final[MAX_BOGEN+1]; /* de lijst van bogen die nadat er eerst nontriviale automorphismen waren die dan
verdwenen nog gericht moeten worden.*/
BOOG *laatstepositie; /* The last position in one of these lists. Global in order not to have to copy it every time */
/* all the next variables must be evaluated immediately when orbits are constructed. They are reused in the next
iteration */ unsignedchar *operations=NULL; int *root_op=NULL, size_root=0, blocklength, orbitblocklength[MAXN], size_operations=0, number_operations=0; /* these arrays will be dynamically allocated and extended. Operations will be an array with "blocks" of length "blocklength=tobedirected[vertex_in_orbit]+3". They represent an operation the following way: the first entry is the central vertex then the list of vertices from which only edges come in -- ended by INFTY_UCHAR, then the list of vertices to which only outgoing edges go -- ended by INFTY_UCHAR. Then the list of vertices to which incoming AND outgoing edges go -- ended by INFTY_UCHAR.
root has as many entries as there are blocks in operations. root[i]=j means that when computing the orbits of
operations and constructing a union-find-tree, block number j is the root of the tree containing vertex i. */
unsignedchar *remember_operations[MAXN]; //the nonequivalent operations that are stored for every iteration int remember_size[MAXN]; // remember_size[i] is the number of characters allocated for remember_operations[i]
/* OPTIONS */
boolean dummybool; int mingerichtdeg; int remaining_doubles=0; /* hoeveel bogen kunnen ten hoogste nog in allebei
richtingen gericht worden ? */ int watermaxedges; /* what is the theoretical maximum for the number of edges which can directed in a way that doesn't give a conflict with the in/out-degree
bounds. */ int watermaxdeg; /* maxindeg+maxoutdeg */
int is_gericht[MAX_BOGEN][MAX_BOGEN]={{0}}; /* is_gericht[i][j]==1 als {i,j} is al gericht -- anders 0 */ int virtual_gericht[MAX_BOGEN][MAX_BOGEN]={{0}}; /* virtual_gericht[i][j]==1 als het al vastgelegd is dat {i,j} i->j gericht moet worden, 2 als hij j->i gericht moet worden en 0 als het nog niet vastgelegd is. Vastgelegd betekent door de beperkingen van in- en outgraad moet deze boog zo gericht
worden. */ int positie[MAXN][MAXN]; /* the value of positie[i][j] is the position of edge {1,j} in edgelist in case of
nontrivial automorphisms of the underlying undirected graph. */
int virtual_indeg[MAXN], virtual_outdeg[MAXN]; /* de graden als de nog niet gerichte maar al vastgelegde bogen meegerekend worden -- dat mag niet voor kanoniciteit
gebruikt worden */ int nextstep_depth;
#define SWITCHPAR_ORBSIZE 5 // when does serial orbit labelling switch to parallel -- and stay there #define MAXPAR_ORBSIZE 9 // maximum 16, when changing, change following MAXPAROPS too #define MAXPAROPS 19683 // set to minimum 3^MAXPAR_ORBSIZE unsignedint parops[MAXPAROPS]; // in the case of no degree restrictions this could in fact be permanently filled in -- only depending // on the edge orbit size. But the time for filling it in is so small that it is simply not worth it. int number_parops; // an operation is defined as follows: if edge number i in kleinste_orbit must be directed kleinste_orbit[i][0]->kleinste_orbit[i][1], // bits 2i and 2i+1 form the number 0 (so are both 0). For kleinste_orbit[i][0]<-kleinste_orbit[i][1] they form 2, for // a double edge they form 1 ---- so 2-type is always the inverse #define SETOP(op,edge,type) ((op) = ((op) & ~(3<<((edge)<<1))) | (type)<<((edge)<<1)) #define GETTYPE(op,edge) (((op)>>((edge)<<1))&3)
/* for the following macros global variables _x_ and _y_ are used */ int _x_, _y_; #define BILDBOOG_UNDIR(startboog,bildboog,permnummer) \
{ _x_=generators[permnummer][(startboog)[0]]; _y_=generators[permnummer][(startboog)[1]]; \ if (_x_<_y_) { (bildboog)[0]=_x_; (bildboog)[1]=_y_; } else { (bildboog)[0]=_y_; (bildboog)[1]=_x_; }} #define BILDBOOG(startboog,bildboog,permnummer) \
{ (bildboog)[0]=generators[permnummer][(startboog)[0]]; (bildboog)[1]=generators[permnummer][(startboog)[1]]; } #define POSBILD(startboog,permnummer) \
(positie[generators[permnummer][(startboog)[0]]][generators[permnummer][(startboog)[1]]])
/* a macro that sets numbers a and b equivalent to a */ #define SETEQUIV(a,b) { _y_=(b); while ((_x_=number[_y_])!=_y_) { number[_y_]=(a); _y_=_x_; } number[_y_]=(a); } /* zorgt ervoor dat nummer[a] zo is dat nummer[nummer[a]]=nummer[a] */ #define UPDATENUMBER(a) { _y_=(a); while ((_x_=number[_y_])!=_y_) { number[_y_]=number[_x_]; _y_=number[_y_]; } number[a]=number[_y_]; }
/* the quality of a directed edge */ #define QUALITY(a,b) (((saturated[a]+saturated[b])<<6)+colour[nextstep_depth][a]) // upper bound for the quality if the in- or out-deg of one of the entries is changed by one #define QUALITY_P1(a,b) (((saturated[a]+saturated[b]+1)<<6)+colour[nextstep_depth][a]) // quality if the in- or out-deg of one of the entries are both changed by one #define QUALITY_P2(a,b) (((saturated[a]+saturated[b]+2)<<6)+colour[nextstep_depth][a])
int _marks[MAX_BOGEN], _markvalue=INT_MAX; #define RESETMARKS { int i; if (_markvalue<INT_MAX) _markvalue++; else { _markvalue=1; \ for (i=0;i<MAX_BOGEN;i++) _marks[i]=0; } } #define ISMARKED(i) (_marks[i]==_markvalue) #define UNMARKED(i) (_marks[i]!=_markvalue) #define MARK(i) {_marks[i]=_markvalue;}
longlongint onlevel[MAX_BOGEN]={0};
/* variabelen voor nauty: */ int lab[MAX_BOGEN][MAXN]={{0}}, ptn[MAX_BOGEN][MAXN]={{0}}, orbits[MAX_BOGEN], colour[MAX_BOGEN][MAXN]={{0}}; /* MAX_BOGEN is zeker te veel -- maar wat zou een goede bovengrens zijn ? */ /* in the arrays lab[i][] and ptn[i][] the partition is stored after having labelled i orbits
completely. lab[0][] and ptn[0][] are the partitions as computed for the original graph. */ int rememberorbits[MAXN][MAXN];
int bufferlab[MAXN], bufferptn[MAXN];
graph canong[MAXN]; /* make sure these are always evaluated at once, because they are reused whenever new
edges are directed */
graph workg[MAXN], staticg[MAXN], canong[MAXN];
graph bit_orbit[MAXN]; int deg[MAXN]={0}, outdeg[MAXN]={0}, indeg[MAXN]={0}; int tobedirected[MAXN]={0}; /* how many edges have still to be directed? This can also be used to detect whether an edge that is in both directions is a double edge or just not yet directed in the routines for nontrivial symmetry. It is a double edge if
and only if one of the endpoints is ready with ready defines as: */ #define READY(i) (tobedirected[i]==0) #define NOT_READY(i) (tobedirected[i]>0) int aantal_toppen, aantal_bogen, aantal_gerichte_bogen; int max_doubles; int double_free[MAXN]={0}; /* hoeveel dubbele bogen kan top [i] nog krijgen? */
longlongint addnumber=0; /* How much must the counter be increased if a graph is found. Just relevant in case not all are really constructed and coded.
Then this gives a speedup.*/
//unsigned int waterclusteruse=0UL, water_v_use=0UL;
/****OPTIONS******/ int maxindeg=MAXN, maxoutdeg=MAXN, maxdirectdeg=MAXN, nodegbound=1; int double_allowed=1, direct_output=0;
#define NEXTEL(a,b) (b)=FIRSTBIT((a)& BITMASK(b)) /*(b)++, b=FIRSTBIT((a)<<(b))+(b)*/ #define FORALLELEMENTS(x,y) for (y= FIRSTBIT(x); y<WORDSIZE; NEXTEL((x),(y))) #define FORALLELEMENTS_BOUND(x,y,b) for (y=FIRSTBIT((x)& BITMASK(b)) ; y<WORDSIZE; NEXTEL((x),(y))) /* Attention: bound is a lower bound and the first possible value is bound+1 ! */
#define WRITEUP() { aantal_gerichte_grafen+=addnumber; if (direct_output==1) {writeTcode(workg,aantal_toppen);\ if (addnumber==2) writeTcode_invers(workg,aantal_toppen);} elseif (direct_output==3) \
{writeBcode(workg,aantal_toppen);\ if (addnumber==2) writeBcode_invers(workg,aantal_toppen);} elseif (direct_output==4) \
{writeZcode(workg,aantal_toppen);\ if (addnumber==2) writeZcode_invers(workg,aantal_toppen);} } #define WRITEUP_COUNT() { aantal_gerichte_grafen+=addnumber; }
#ifdef PROCESS void callprocess(graph *g, int aantal_toppen)
{ int i,j; graph gconv[MAXN];
void writeoperation(unsignedchar op[])
{ int i,counter;
fprintf(stderr,"Operation: "); for (i=counter=0; counter<3; i++)
{ if (op[i]==INFTY_UCHAR) { fprintf(stderr," ||"); counter++; } else fprintf(stderr," %d",op[i]);
}
fprintf(stderr,"\n");
}
void writelist(int list[])
{ int i;
for (i=0;list[i]>=0;i++) fprintf(stderr,"%d ",list[i]); fprintf(stderr,"\n");
}
void writegraph(graph *g) // for graphs as they are used in the vertexorbit routines
{ int i,j; staticint counter=0;
fprintf(stderr,"---------------------------------------------------------\n");
fprintf(stderr,"Graph number %d with %d vertices\n",++counter,aantal_toppen); for (i=0; i<aantal_toppen;i++)
{ fprintf(stderr,"%d:",i);
FORALLELEMENTS(g[i],j)
{ fprintf(stderr," %d", j); if (READY(i) || READY(j))
{ if (ISELEMENT(g+j,i)) fprintf(stderr,"[d]"); else fprintf(stderr,"[s]");
} else fprintf(stderr,"[ng]");
}
fprintf(stderr,"\n");
}
fprintf(stderr,"---------------------------------------------------------\n");
}
void sammle_permutationen(int count, int perm[], int orbits[], int numorbits, int stabvertex, int n)
{
memcpy(generators+number_of_generators,perm,sizeof(int)*n);
void init_nauty_options() /* initialises the nauty variables */
{ /* tc_level = 0 is not the default in the most recent editions of nauty. However, it is better for very small
graphs so we set it here. */
int lese_multicode(unsignedchar**code, int *codelaenge, FILE *fil)
/* Liest den code und gibt EOF zurueck, wenn das Ende der Datei erreicht ist, 1 sonst. Der Speicher fuer den code wird alloziert, falls noetig,
was entschieden wird anhand der lokalen Variablen maxknotenzahl */
{ staticint maxknotenzahl= -1; int codel, gepuffert=0; int knotenzahl,a=0,b=0,nuller; unsignedchar ucharpuffer;
if ((knotenzahl=getc(fil))==EOF) return EOF; if (knotenzahl==0) { fprintf(stderr,"Umschaltung auf short noch nicht implementiert.\n"); exit(0);
}
nuller=0; codel=1; if (knotenzahl=='>') /* koennte ein header sein -- oder 'ne 62, also ausreichend fuer
unsigned char */
{ gepuffert=1;
a=getc(fil); if(a==0) nuller++;
b=getc(fil); if(b==0) nuller++; /* jetzt wurden 3 Zeichen gelesen */ if ((a=='>') && (b=='m')) /*garantiert header*/
{ while ((ucharpuffer=getc(fil)) != '<') {} /* noch zweimal: */ ucharpuffer=getc(fil); if (ucharpuffer!='<') { fprintf(stderr,"Problems with header -- single '<'\n"); exit(1); } if ((knotenzahl=getc(fil))==EOF) return EOF; /* kein graph drin */
} /* else kein header */
}
if (knotenzahl > maxknotenzahl)
{ if (*code) free(*code);
*code=(unsignedchar *)malloc((knotenzahl*(knotenzahl-1)/2+knotenzahl)*sizeof(unsignedchar)); if (code==NULL) { fprintf(stderr,"Do not get memory for code\n"); exit(0); }
maxknotenzahl=knotenzahl;
}
(*code)[0]=knotenzahl; if (gepuffert) { codel=3; (*code)[1]=a; (*code)[2]=b; }
while (nuller<knotenzahl-1)
{ if (((*code)[codel]=getc(fil))==0) nuller++;
codel++; }
*codelaenge=codel; return 1;
}
void usage(char name[])
{
fprintf(stderr,"\nUsage: %s [ix] [oy] [m] [T] [C] [B] [Z] [S]\n",name);
fprintf(stderr," Read undirected graphs and orient them in various ways.\n\n");
fprintf(stderr,"The option ix restricts the maximum indegree to x.\n");
fprintf(stderr,"The option oy restricts the maximum outdegree to y.\n");
fprintf(stderr,"The default maximum in- and out-degrees are unlimited. \n");
fprintf(stderr,"T means: Output directed graphs in T-code -- for details see header\n");
fprintf(stderr,"B means: Output directed graphs in binary code -- for details see header\n");
fprintf(stderr,"Z means: Output directed graphs in digraph6 code\n");
fprintf(stderr,"C means: Do really construct all the directed graphs in memory,\n" " but don't output them (default)\n");
fprintf(stderr,"S means that for each edge only one direction must be chosen -- not both.\n");
fprintf(stderr,"Default is that both are allowed\n");
fprintf(stderr," -- so the edge a-b can become a->b AND b->a in the same output graph.\n");
fprintf(stderr,"m means: read multicode instead of g6 code \n"); exit(1);
void fill_edgelist() /* writes the edges in g into the list in a lexicographic way. Is used if no bounds for in- and out-degree are given
or few edges are left. */
{ int i,j,end;
for (i=j=0; i<aantal_toppen; i++) /* j is de positie in edgelist */
{ indeg_free[i]=maxindeg-indeg[i]; outdeg_free[i]=maxoutdeg-outdeg[i];
FORALLELEMENTS_BOUND(workg[i],end,i)
{
edgelist[j][0]=i; edgelist[j][1]=end; j++;
}
}
void fill_edgelist_order() /* writes the edges in gg into the list in a way that for i in 1...numberedges we always have that the starting point of edge i has the minimum degree of all vertices in the graph you get when removing all edges i+1...aantal_bogen from gg. The global variable positie is not assigned any values.
*/
{ int last_positie,i, buffer, olddeg, newdeg, beste, start, end; int list[MAXN][MAXN], listlen[MAXN]={0}; /* list[i] contains the vertices with degree i at that moment */ int toppositie[MAXN], bufferdeg[MAXN]; /* toppositie[i] is de positie van top i in de lijst */ int buren[MAXN]; /* de buren van de top waaraan gewerkt wordt */
graph dummy[MAXN];
last_positie=aantal_bogen-1; while (last_positie>=0)
{ for (beste=1;listlen[beste]==0; beste++);
(listlen[beste])--; /* zo wordt hij ook verwijderd */
start=list[beste][listlen[beste]]; for (i=0; i<beste; i++) /* alle buren opslaan*/
{ end=buren[i]=FIRSTBIT(dummy[start]); DELELEMENT(dummy+start,end); }
for (i=0; i<beste; i++) /* alle buren */
{
end=buren[i];
edgelist[last_positie][0]=start; edgelist[last_positie][1]=end;
last_positie--;
DELELEMENT(dummy+end,start); /* de buur verhuizen: */
olddeg=bufferdeg[end];
(bufferdeg[end])--;
newdeg=bufferdeg[end]; /* uit de oude lijst verwijderen */ if (listlen[olddeg]==1) listlen[olddeg]=0; else { (listlen[olddeg])--;
buffer= list[olddeg][listlen[olddeg]];
list[olddeg][toppositie[end]]=buffer;
toppositie[buffer]=toppositie[end]; } /* tot de nieuwe toevoegen */ if (newdeg)
{ list[newdeg][listlen[newdeg]]=end;
toppositie[end]=listlen[newdeg];
(listlen[newdeg])++;
}
}
void trivlabels_nowrite_nodouble(BOOG *positie)
{ int start, end;
start=(*positie)[0]; end=(*positie)[1];
if (positie==laatstepositie)
{ if (indeg_free[start] && outdeg_free[end]) /* bogen kan end->start gericht worden */
{
WRITEUP_COUNT();
} if (indeg_free[end] && outdeg_free[start]) /* bogen kan start->end gericht worden */
{
WRITEUP_COUNT();
}
return;
}
/* else -- niet op laatste positie */
if (indeg_free[start] && outdeg_free[end]) /* bogen kan end->start gericht worden */
{
(indeg_free[start])--; (outdeg_free[end])--;
trivlabels_nowrite_nodouble(positie+1);
(indeg_free[start])++; (outdeg_free[end])++;
} if (indeg_free[end] && outdeg_free[start]) /* bogen kan start->end gericht worden */
{
(indeg_free[end])--; (outdeg_free[start])--;
trivlabels_nowrite_nodouble(positie+1);
(indeg_free[end])++; (outdeg_free[start])++;
} return;
}
void trivlabels_nowrite(BOOG *positie)
{ int start, end, counter;
if (!remaining_doubles) { trivlabels_nowrite_nodouble(positie); return; }
start=(*positie)[0]; end=(*positie)[1];
if (positie==laatstepositie)
{ if (indeg_free[start] && outdeg_free[end]) /* bogen kan end->start gericht worden */
{ counter=1;
WRITEUP_COUNT();
} else counter=0; if (indeg_free[end] && outdeg_free[start]) /* bogen kan start->end gericht worden */
{ counter++;
WRITEUP_COUNT();
} if (remaining_doubles && (counter==2)) WRITEUP_COUNT();
return;
}
/* else -- niet op laatste positie */
if (indeg_free[start] && outdeg_free[end]) /* bogen kan end->start gericht worden */
{ counter=1;
(indeg_free[start])--; (outdeg_free[end])--;
trivlabels_nowrite(positie+1);
(indeg_free[start])++; (outdeg_free[end])++;
} else counter=0; if (indeg_free[end] && outdeg_free[start]) /* bogen kan start->end gericht worden */
{ counter++;
(indeg_free[end])--; (outdeg_free[start])--;
trivlabels_nowrite(positie+1);
(indeg_free[end])++; (outdeg_free[start])++;
} if (remaining_doubles && (counter==2) && double_free[start] && double_free[end])
{
(indeg_free[end])--; (outdeg_free[start])--;
(indeg_free[start])--; (outdeg_free[end])--;
double_free[start]--; double_free[end]--;
remaining_doubles--;
trivlabels_nowrite(positie+1);
remaining_doubles++;
double_free[start]++; double_free[end]++;
(indeg_free[end])++; (outdeg_free[start])++;
(indeg_free[start])++; (outdeg_free[end])++;
} return;
}
void trivlabels(BOOG *positie)
{ int start, end, counter;
start=(*positie)[0]; end=(*positie)[1];
if (positie==laatstepositie)
{ if (indeg_free[start] && outdeg_free[end]) /* bogen kan end->start gericht worden */
{ counter=1; //(indeg_free[start])--; (outdeg_free[end])--;
DELELEMENT(workg+start,end);
MAYBEPROCESS;
WRITEUP();
ADDELEMENT(workg+start,end); //(indeg_free[start])++; (outdeg_free[end])++;
} else counter=0; if (indeg_free[end] && outdeg_free[start]) /* bogen kan start->end gericht worden */
{ counter++; //(indeg_free[end])--; (outdeg_free[start])--;
DELELEMENT(workg+end,start);
MAYBEPROCESS;
WRITEUP();
ADDELEMENT(workg+end,start); //(indeg_free[end])++; (outdeg_free[start])++;
} if (remaining_doubles && (counter==2)) { remaining_doubles--; MAYBEPROCESS; WRITEUP(); remaining_doubles++; } return;
}
/* else -- niet op laatste positie */
if (indeg_free[start] && outdeg_free[end]) /* bogen kan end->start gericht worden */
{ counter=1;
(indeg_free[start])--; (outdeg_free[end])--;
DELELEMENT(workg+start,end);
trivlabels(positie+1);
ADDELEMENT(workg+start,end);
(indeg_free[start])++; (outdeg_free[end])++;
} else counter=0; if (indeg_free[end] && outdeg_free[start]) /* bogen kan start->end gericht worden */
{ counter++;
(indeg_free[end])--; (outdeg_free[start])--;
DELELEMENT(workg+end,start);
trivlabels(positie+1);
ADDELEMENT(workg+end,start);
(indeg_free[end])++; (outdeg_free[start])++;
} if (remaining_doubles && (counter==2) && double_free[start] && double_free[end])
{
(indeg_free[end])--; (outdeg_free[start])--;
(indeg_free[start])--; (outdeg_free[end])--;
double_free[start]--; double_free[end]--;
remaining_doubles--;
trivlabels(positie+1);
remaining_doubles++;
double_free[start]++; double_free[end]++;
(indeg_free[end])++; (outdeg_free[start])++;
(indeg_free[start])++; (outdeg_free[end])++;
}
return;
}
void trivlabels_init(BOOG *positie)
{ int remember;
if (!direct_output)
{ if (nodegbound)
{ remember=addnumber; if (remaining_doubles) { for (;positie<=laatstepositie;positie++) addnumber*=3;} else { for (;positie<=laatstepositie;positie++) addnumber*=2;}
WRITEUP();
addnumber=remember; return;
} /* else */ if (!remaining_doubles) trivlabels_nowrite_nodouble(positie); else trivlabels_nowrite(positie); return;
} /* else */
trivlabels(positie); return;
}
aantal_grafen_met_triv_group++; if (nodegbound) fill_edgelist(); else fill_edgelist_order(); for (i=0;i<aantal_toppen;i++) double_free[i]=maxdirectdeg-deg[i];
laatstepositie=edgelist+aantal_bogen-1; if (maxoutdeg==maxindeg) /* then every valid graph for one direction is also valid with
the directions reversed */
{
addnumber=2;
start=edgelist[0][0]; end=edgelist[0][1]; /* is_gericht[][] isn't used in the directing routine */
(outdeg_free[start])--; (indeg_free[end])--; //virtual_outdeg[start]++; virtual_indeg[end]++;
DELELEMENT(workg+end,start);
aantal_gerichte_bogen=1;
trivlabels_init(edgelist+1);
aantal_gerichte_bogen=0;
ADDELEMENT(workg+end,start);
(outdeg_free[start])++; (indeg_free[end])++;
int compute_image_operation(unsignedchar image[], unsignedchar original[], int generator[]) // returns 1 if image and original differ and 0 otherwise
{ int in_image[MAXN], out_image[MAXN], double_image[MAXN]; int run, num_in, num_out, num_double, change;
change=0;
image[0]=generator[original[0]];
if (original[0]!=image[0]) change=1;
for (run=1, num_in=0; original[run]!=INFTY_UCHAR; run++)
{ in_image[num_in]=generator[original[run]]; num_in++; }
for (run++, num_out=0; original[run]!=INFTY_UCHAR; run++)
{ out_image[num_out]=generator[original[run]]; num_out++; }
if (num_in>1) sort_decreasing(in_image,num_in); if (num_out>1) sort_decreasing(out_image,num_out); if (num_double>1) sort_decreasing(double_image,num_double);
for (run=1 ; num_in; run++)
{ num_in--; image[run]=in_image[num_in]; if (original[run]!=image[run]) change=1;}
image[run]=INFTY_UCHAR; run++; for ( ; num_out; run++)
{ num_out--; image[run]=out_image[num_out]; if (original[run]!=image[run]) change=1;}
image[run]=INFTY_UCHAR; run++; for ( ; num_double; run++)
{ num_double--; image[run]=double_image[num_double]; if (original[run]!=image[run]) change=1;}
image[run]=INFTY_UCHAR; //for ( ;run<blocklength;run++) image[run]=INFTY_UCHAR;
void construct_operations_final(int list[], int decided[],unsignedchar buffer_op[], int positie, int numberin, int numberout)
{ int i, center; staticunsignedchar *buffer;
/* here about position "positie" of the operation is decided. If "positie==blocklength" we are done.
This function decides on double edges.
*/
if (size_operations<=((number_operations+1)*blocklength))
{ buffer=operations;
operations=malloc((size_t)(2*size_operations)); if (operations==NULL)
{ fprintf(stderr,"Can't allocate %d bytes for operations -- exiting\n",2*size_operations); exit(1); }
memcpy(operations,buffer,size_operations);
free(buffer);
size_operations *= 2;
}
center=buffer_op[0]; // double_free[center] must be tested inside the loop to make sure // that it only leads to rejection if all edges have already been assigned
if (double_free[center])// otherwise the rest was filled in in the previous step for (i=0;list[i]>=0;i++) if (!decided[i])
{ if ( double_free[center] && double_free[list[i]] &&
((outdeg[center]+numberout < maxoutdeg) && (indeg[list[i]]<maxindeg)) &&
((indeg[center]+numberin < maxindeg) && (outdeg[list[i]]<maxoutdeg)))
{
buffer_op[positie]=list[i];
positie++; numberout++; numberin++;
} elsereturn;
}
void construct_operations_out(int list[],int liststart, int decided[],unsignedchar buffer_op[], int positie, int numberin, int numberout, int lowerlimit_outdeg, int mindouble)
{ int i, center, buffer; staticint outdeg_larger;
/* here about position "positie" of the operation is decided. If "positie==blocklength" we are done.
This function decides on outgoing edges.
liststart ensures that no earlier -- smaller -- entries are written after larger ones unless there has been an INFTY_UCHAR to separate them.
*/
center=buffer_op[0];
if (numberout==0) // just test once
{ buffer=outdeg[center] + tobedirected[center] -numberin; if (buffer<lowerlimit_outdeg) return; else
{ if (buffer>lowerlimit_outdeg) outdeg_larger=1; else outdeg_larger=0; }
}
if (double_free[center]) // then there is still something to decide
{ if (outdeg[center]+numberout < maxoutdeg) for (i=liststart;list[i]>=0;i++) if (!decided[i])
{ if (indeg[list[i]]<maxindeg)
{ decided[i]=1;
buffer_op[positie]=list[i];
construct_operations_out(list,i+1,decided,buffer_op,positie+1,
numberin, numberout+1,lowerlimit_outdeg,mindouble);
decided[i]=0;
}
}
buffer_op[positie]=INFTY_UCHAR; if (outdeg_larger ||
(indeg[center]+outdeg[center]+(tobedirected[center]<<1)-numberin-numberout-deg[center]>=mindouble)) // looks complicated but is just the final number of double edges
construct_operations_final(list,decided,buffer_op,positie+1,
numberin, numberout); return;
} //else // all the undecided entries have to be outgoing
for (i=liststart;list[i]>=0;i++) if (!decided[i])
{ if ((outdeg[center]+numberout < maxoutdeg) && (indeg[list[i]]<maxindeg))
{
buffer_op[positie]=list[i];
positie++; numberout++;
} elsereturn;
}
buffer_op[positie]=INFTY_UCHAR;
construct_operations_final(list,decided,buffer_op,positie+1,
numberin, numberout); return;
void construct_operations_in(int list[],int liststart, int decided[],unsignedchar buffer_op[], int positie, int numberin, int numberout, int lowerlimit_outdeg, int mindouble)
{ int i, center;
/* here about position "positie" of the operation is decided. If "positie==blocklength" we are done.
This function decides on incoming edges.
liststart ensures that no earlier -- smaller -- entries are written after larger ones unless there has been an INFTY_UCHAR to separate them.
*/
center=buffer_op[0];
if(indeg[center]+numberin < maxindeg) for (i=liststart;list[i]>=0;i++) if (!decided[i])
{ if (outdeg[list[i]]<maxoutdeg)
{ decided[i]=1;
buffer_op[positie]=list[i];
construct_operations_in(list,i+1,decided,buffer_op,positie+1,
numberin+1, numberout,lowerlimit_outdeg,mindouble);
decided[i]=0;
}
}
buffer_op[positie]=INFTY_UCHAR;
construct_operations_out(list,0,decided,buffer_op,positie+1,
numberin, numberout,lowerlimit_outdeg,mindouble); return;
}
void construct_extensions(int still_open[], int orbit[], graph touched, int first_in_orbit, graph sameorbit)
{ int top, top2, j, end, list[MAXN], decided[MAXN], error, lowerlimit_outdeg, readylist[MAXN], *readyrun, dummy; int minout, do_double=0, i, mindouble; unsignedchar buffer[MAXN+4];
graph pre_free, tbd1, bufferg, freevertices; // tbd1 is the set of vertices with exactly one edge to be directed left // decided[i]=0 if and only if about list[i] is already decided //
tbd1=(graph)0;
pre_free= all & ~touched;
readyrun=readylist;
lowerlimit_outdeg=0;
// Since either all or no vertices in the orbit have edges going to untouched vertices, one could split // this function to save some tests. But for the moment it should better stay like this
// if a vertex v has at least two edges to be directed and one of them leads to a vertex in the same // orbit that will be ready afterwards, then directing the edges around v would not have a // reverse operation. So these operations must be avoided.
// if only one edge is left to be directed and the endpoint is in the same orbit is must be directed double // or outgoing.
if (first_in_orbit)
{ for (number_operations=0; (top=(*still_open))>=0; still_open++) // number_operations is global, so it must not be passed to the next functions //if (NOT_READY(top))
{ if (tobedirected[top]==1)
{ for (j= FIRSTBIT(workg[top]); READY(j) ; NEXTEL((workg[top]),(j))); if (ISELEMENT(&sameorbit,j) && (tobedirected[j]==1)) construct_operations_one(top, j, 0, 1, double_allowed); else construct_operations_one(top, j, 1, 1, double_allowed);
} else if (!(staticg[top] & tbd1))
{ for (end=0, j= FIRSTBIT(workg[top]); (j<WORDSIZE); NEXTEL((workg[top]),(j)))
{ if (NOT_READY(j)) { list[end]=j; decided[end]=0; end++; }} // the edge has already been directed if and only if the other top has already been finished
list[end]= -1;
buffer[0]= top;
construct_operations_in(list,0,decided,buffer,1,0,0,0,0);
}
} return;
}
// else: not first in orbit
minout=INT_MAX;
if (staticg[orbit[0]]& ~touched) // all operations will contain edges to free vertices
{ for (j=0; (top=orbit[j])>=0; j++)
{ if (NOT_READY(top))
{ //ADDELEMENT(&sameorbit,top); if (tobedirected[top]>=2) ADDELEMENT(&pre_free,top); // will stay free unless is chosen top else ADDELEMENT(&tbd1,top);
} else// that is: ready
{ *readyrun=top; readyrun++; } // minout would not be used
}
} else
{
do_double=double_allowed; for (j=0; (top=orbit[j])>=0; j++)
{ if (NOT_READY(top))
{ //ADDELEMENT(&sameorbit,top); if (tobedirected[top]>=2) ADDELEMENT(&pre_free,top); // will stay free unless is chosen top else ADDELEMENT(&tbd1,top);
} else// that is: ready
{ bufferg=(workg[top]&sameorbit);
*readyrun=top; readyrun++; if ((bufferg) && (outdeg[top]<minout)) { minout=outdeg[top]; do_double=double_allowed; } if (do_double && (outdeg[top]==minout))
{
FORALLELEMENTS(bufferg,i)
{ if (READY(i) && !ISELEMENT(workg+i,top)) // no double edge, so double edges inside the orbit would not be accepted
{ do_double=0; i=WORDSIZE-1; } // BDM changed "WORDSIZE" to "WORDSIZE-1"
}
}
}
}
}
*readyrun= -1;
for (number_operations=0; (top=(*still_open))>=0; still_open++) // number_operations is global, so it must not be passed to the next functions //if (NOT_READY(top))
{ if (tobedirected[top]==1) // already ready vertices can't have less...
{ for (j= FIRSTBIT(workg[top]); READY(j) ; NEXTEL((workg[top]),(j))); if (ISELEMENT(&sameorbit,j) && (tobedirected[j]==1))
{ bufferg=workg[j]&sameorbit; DELELEMENT(&bufferg,top); if (outdeg[top]<minout)
{ if ((bufferg== (graph)0) || (outdeg[j]>outdeg[top])) // otherwise the end would be better
{ if (outdeg[top]==minout-1) construct_operations_one(top, j, 0, 1, do_double); else construct_operations_one(top, j, 0, 1, double_allowed);
} else if ((outdeg[j]==outdeg[top]) && double_allowed)
{ if (outdeg[top]==minout-1) construct_operations_one(top, j, 0, 0, do_double); else construct_operations_one(top, j, 0, 0, double_allowed);
}
}
} else construct_operations_one(top, j, 1, 1, double_allowed);
} else
{ if (!(staticg[top] & tbd1))
{
error=0;
freevertices = pre_free | (tbd1 & ~staticg[top]);
DELELEMENT1(&freevertices,top); // now freevertices is the set of vertices that will still be free after the addition of top for (lowerlimit_outdeg=mindouble=0, readyrun=readylist; !error && (top2= *readyrun)>=0; readyrun++)
{ bufferg= freevertices & staticg[top2]; if ((dummy=POPCOUNT(bufferg)))// will be a candidate afterwards
{ if (dummy < tobedirected[top]) error=1; else
{ if (dummy == tobedirected[top])
{ if (outdeg[top2]>lowerlimit_outdeg)
{ lowerlimit_outdeg=outdeg[top2]; mindouble=indeg[top2]+outdeg[top2]-deg[top2]; } else if (double_allowed && (outdeg[top2]==lowerlimit_outdeg))
{ if (indeg[top2]+outdeg[top2]-deg[top2]>mindouble)
mindouble=indeg[top2]+outdeg[top2]-deg[top2];
}
}
}
}
} if (!error)
{ for (end=0, j= FIRSTBIT(workg[top]);
(j<WORDSIZE) ; NEXTEL((workg[top]),(j))) if (NOT_READY(j)) { list[end]=j; decided[end]=0; end++; }
// the edge has already been directed if and only if the other top has already been finished
list[end]= -1;
buffer[0]= top; if (!error) construct_operations_in(list,0,decided,buffer,1,0,0,lowerlimit_outdeg,mindouble);
}
}
}
} return;
}
int compare_op(unsignedchar *op1, unsignedchar *op2) // returns 0 if operations are the same, something negative if op1 // is lexicographically smaller and something positive else
{ int counter;
if (*op1 != *op2) return (*op1 - *op2); // can never be INFTY_UCHAR
op1++; op2++;
counter=2; while (*op1 == *op2)
{ if (*op1 == INFTY_UCHAR) { if (!counter) return 0; counter--; }
op1++; op2++;
} return (*op1 - *op2);
}
int compute_orbits() /* compute the orbits of the group described by the global variable generators[][] on the operations described in the global variable operations. Returns the number of orbits.
*/
{ int i,j,buffer,run, orbits, index; unsignedchar image[MAXN+4];
if (number_operations>size_root)
{ size_root=number_operations;
free(root_op);
root_op=malloc((size_t)size_root*sizeof(int));
} if (root_op==NULL)
{ fprintf(stderr,"Can't allocate %d items for root_op -- exiting.\n",size_root); exit(0); }
for (i=0; i<number_operations; i++) root_op[i]=i;
for (i=0, orbits=number_operations; i<number_operations; i++) for (j=0; j<number_of_generators; j++)
{ if (compute_image_operation(image, operations+(i*blocklength), generators[j]))
{
index=search_op(image); //unite: while ((buffer=root_op[index])!=index) index=buffer; // no path compression here
run=i; while ((buffer=root_op[run])!=run) { root_op[run]=index; run=buffer; } if (run!=index)
{ orbits--; // the roots were different -- two trees are united
root_op[run]=index; }
}
} return orbits;
}
int canonical(unsignedchar operation[], int vertexorbit[], int *newgroup, int orbitid, graph touched)
// checks whether the operation is canonical in vertexorbit. Returns 1 if yes, 0 otherwise. // In newgroup it is stored whether a new group has been computed for this (*newgroup=1) or not // (*newgroup=0)
// If no edges go to vertices outside the set of vertices in already "finished" orbits or to vertices in // this orbit with still undecided edges, the canonical operation is directing one edge {x,y} x->y or x<->y. // Edges completely in this orbit and obtained by operations center x end y are chosen as follows: // First the starting point x is chosen with minimum outdegree. // For edges with starting points with the same minimal outdegree, edges are preferred which are single. // If the tested operation is double and single operations exist, the operation is rejected. // Among the remaining edges those are chosen where y has minimum indegree. // Among the remaining those are chosen where x is in the orbit of the one with smallest canonical number // (after the operation). y is the neighbour with smallest indegree -- and among those with the same indegree // the one with smallest canonical number (in the same given numbering of course).
// Otherwise: // Define a vertex to be free if it is not in an orbit that was already chosen (we call that untouched) // or it is in such an orbit but not all edges to it are already directed (that can only be in the last orbit). // In this case the central vertex of the canonical operation is one that is chosen in steps: // Among all with edges to free vertices it is chosen as one with a minimum number of these. // Among those one with max outdegree is chosen. // Among those one with maximum number of double edges is chosen. // Finally a more complicated artificial criterion is tested -- not suitable for a look ahead during the construction // of possible operations // Among the remaining, one with smallest canonical label is chosen.
{ int i, j, k, min, top, buffer, candidatelist[MAXN], number_candidates, center;
graph free, candidates, bufferset, dummy, sameorbit, bit_startlist; int edgelistcounter, edgelist[2*MAXN*MAXN][2], orb[2*MAXN*MAXN], canoncenter, canonend, endvertex=0; int finished[MAXN], finishedptn[MAXN], testoutdeg; int startlist[MAXN], *startrun; // list of possible starts in case of no candidates int numberends, ends[MAXN][MAXN], endin=0, *run; int *colour; longlongint buffer2, k2;
free= all & (~touched); // in fact all ready vertices are candidates if one has a non-empty // intersection with this as they are all in the same orbit -- so have the same number of edges // to this set
sameorbit=bit_orbit[orbitid];
for (run=vertexorbit; (*run)>=0; run++) if (NOT_READY(*run)) ADDELEMENT(&free,(*run));
center=operation[0];
candidates= (graph)0;
number_candidates=0;
min=INFTY_UCHAR; for (i=0; (top=vertexorbit[i])>=0; i++) if (READY(top) && (bufferset=(staticg[top] & free)))
{
buffer=POPCOUNT(bufferset); if (buffer<min) { candidates= (graph)0; ADDELEMENT(&candidates,top);
candidatelist[0]=top; number_candidates=1; min=buffer; } else if (buffer==min)
{ ADDELEMENT(&candidates,top); candidatelist[number_candidates]=top; number_candidates++; }
} if (number_candidates && !ISELEMENT(&candidates,center))
{ return 0;
}
if (number_candidates==1) { *newgroup=0; return 1; }
testoutdeg=outdeg[center]; // First canonicity criterion: maximum outdegree if (number_candidates)
{ for (i=0; i<number_candidates; )
{ k= outdeg[candidatelist[i]]-testoutdeg; if (k>0)
{ return 0;} else
{ if (k<0) { DELELEMENT(&candidates,candidatelist[i]);
number_candidates--;
candidatelist[i]=candidatelist[number_candidates]; } else i++; // otherwise the new element has to be tested
}
}
if (number_candidates==1) { *newgroup=0; return 1; } // end number_candidates>0
if (remaining_doubles != max_doubles) // there are double edges
{ buffer=indeg[center]+outdeg[center]-deg[center]; // number of double edges starting at center for (i=0; i<number_candidates; )
{ j=candidatelist[i];
k= indeg[j]+outdeg[j]-deg[j]-buffer; if (k>0)
{ return 0;} else
{ if (k<0) { DELELEMENT(&candidates,candidatelist[i]);
number_candidates--;
candidatelist[i]=candidatelist[number_candidates]; } else i++; // otherwise the new element has to be tested
}
}
} if (number_candidates==1) { *newgroup=0; return 1; }
//OK-- last try. A colour that can not really be used to avoid unnecessary operation in advance
colour=rememberorbits[orbitid];
dummy=workg[center]&(free | sameorbit);
buffer2=1LL;
FORALLELEMENTS(dummy,i)
{ buffer2 *= (tobedirected[i]<<12)+(indeg[i]<<9)+(outdeg[i]<<6)+(colour[i]<<3)+deg[i]+1;
buffer2 = buffer2%63241LL;
}
if (k2>0LL)
{ return 0;} else
{ if (k2<0LL) { DELELEMENT(&candidates,candidatelist[i]);
number_candidates--;
candidatelist[i]=candidatelist[number_candidates]; } else i++; // otherwise the new element has to be tested
}
}
if (number_candidates==1) { *newgroup=0; return 1; }
} // end number candidates > 0
// OK -- we have to work. First the canonical form must be computed: // Problem: the difference between double edges and edges that have not yet been directed cannot be // detected in the datastructure graph. To this end we put "ready" vertices in an extra partition -- // this way double edges cannot be mapped on undirected edges. This is only necessary if there are double edges.
else// (number_candidates==0) // finished vertices in this orbit always have an internal edge!
{ startrun=startlist;
bit_startlist=(graph)0; if (operation[2]==INFTY_UCHAR) // a double edge was added and all single edges starting at same outdeg are better
{ endvertex= operation[3];
endin=indeg[endvertex]; for (i=0; (k=vertexorbit[i])>=0; i++) if (READY(k) && (workg[k] & sameorbit))// has _outgoing_ edge into same orbit // the end points are automatically ready and all in the same orbit -- otherwise there would have been a candidate
{ if (outdeg[k]<testoutdeg) return 0; if (outdeg[k]==testoutdeg)
{ numberends=0;
bufferset= workg[k]&sameorbit;
FORALLELEMENTS(bufferset,j)
{ if (!ISELEMENT(workg+j,k)) return 0; // else: also double edge if (indeg[j]<endin) return 0; if (indeg[j]==endin) { ends[k][numberends]=j; numberends++; }
} if (numberends)
{ ends[k][numberends]= -1;
*startrun=k; startrun++; ADDELEMENT(&bit_startlist,k);
}
}
}
} else// het is dus een single edge
{ endvertex= operation[2];
endin=indeg[endvertex]; for (i=0; (k=vertexorbit[i])>=0; i++) if (READY(k) && (workg[k] & sameorbit))// has _outgoing_ edge into same orbit // the end points are automatically ready and all in the same orbit -- otherwise there would have been a candidate
{ if (outdeg[k]<testoutdeg) return 0; if (outdeg[k]==testoutdeg)
{ numberends=0;
bufferset= workg[k]&sameorbit;
FORALLELEMENTS(bufferset,j)
{ if (!ISELEMENT(workg+j,k)) // also single
{ if (indeg[j]<endin) return 0; if (indeg[j]==endin) { ends[k][numberends]=j; numberends++; }
}
} if (numberends)
{ ends[k][numberends]= -1;
*startrun=k; startrun++; ADDELEMENT(&bit_startlist,k);
}
}
}
}
*startrun= -1;
dummy=workg[center] & sameorbit; if ((startlist[1]== -1) && (ends[center][1]== -1)) { *newgroup=0; return 1; } // maar 1 boog mogelijk
}
if (double_allowed) // distinguish single and double edges for nauty()
{ for (i=j=k=0; i<aantal_toppen;i++)
{ top=lab[orbitid][i]; if (READY(top)) { finished[k]=top; finishedptn[k]=1; k++; } else { bufferlab[j]=top; bufferptn[j]=1; j++; } if (ptn[orbitid][i]==0)
{ if (j) bufferptn[j-1]=0; if (k) finishedptn[k-1]=0;
}
} for (i=0;i<k;i++,j++) { bufferptn[j]=finishedptn[i]; bufferlab[j]=finished[i]; }
} else
{
memcpy(bufferlab,lab[orbitid],aantal_toppen*sizeof(int));
memcpy(bufferptn,ptn[orbitid],aantal_toppen*sizeof(int));
}
number_of_generators=0;
nauty(workg,bufferlab,bufferptn,NULL,orbits,&options_directed_canon,&stats,workspace,100*MAXN,1,aantal_toppen,canong);
*newgroup=1;
//nautyuse[number_candidates]++;
// Being a free vertex is invariant under automorphisms and the set of edges to free vertices depends only // on the central vertex. So the central vertex characterizes the whole operation and orbits of central // vertices uniquely correspond to orbits of operations.
if (number_candidates>1)
{ for (i=0; !ISELEMENT(&candidates,bufferlab[i]); i++); // zoekt candidaat met kleinste kanonische label
// only remaining case: number_candidates=0 -- that is: only one edge between 2 vertices of vertexorbit was added
// operation must be outgoing or double edge -- but this is guaranteed by the construction of operations
// for (i=0; !(READY(bufferlab[i]) && ISELEMENT(&sameorbit,bufferlab[i]) && (workg[bufferlab[i]] & sameorbit)) for (i=0; !ISELEMENT(&bit_startlist,bufferlab[i]); i++); // Now we have the smallest vertex in vertexorbit with an edge to another vertex in vertexorbit
canoncenter=bufferlab[i];
if (orbits[center]!=orbits[canoncenter]) { return 0; }
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
