orca.cpp

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <ctime>
#include <iostream>
#include <fstream>
#include <set>
#include <unordered_map>
#include <algorithm>
using namespace std;


typedef long long int64;
typedef pair<int,int> PII;
typedef struct { int first, second, third; } TIII;

struct PAIR {
	int a, b;
	PAIR(int a0, int b0) { a=min(a0,b0); b=max(a0,b0); }
};
bool operator<(const PAIR &x, const PAIR &y) {
	if (x.a==y.a) return x.b<y.b;
	else return x.a<y.a;
}
bool operator==(const PAIR &x, const PAIR &y) {
	return x.a==y.a && x.b==y.b;
}
struct hash_PAIR {
	size_t operator()(const PAIR &x) const {
		return (x.a<<8) ^ (x.b<<0);
	}
};

struct TRIPLE {
	int a, b, c;
	TRIPLE(int a0, int b0, int c0) {
		a=a0; b=b0; c=c0;
		if (a>b) swap(a,b);
		if (b>c) swap(b,c);
		if (a>b) swap(a,b);
	}
};
bool operator<(const TRIPLE &x, const TRIPLE &y) {
	if (x.a==y.a) {
		if (x.b==y.b) return x.c<y.c;
		else return x.b<y.b;
	} else return x.a<y.a;
}
bool operator==(const TRIPLE &x, const TRIPLE &y) {
	return x.a==y.a && x.b==y.b && x.c==y.c;
}
struct hash_TRIPLE {
	size_t operator()(const TRIPLE &x) const {
		return (x.a<<16) ^ (x.b<<8) ^ (x.c<<0);
	}
};

unordered_map<PAIR, int, hash_PAIR> common2;
unordered_map<TRIPLE, int, hash_TRIPLE> common3;
unordered_map<PAIR, int, hash_PAIR>::iterator common2_it;
unordered_map<TRIPLE, int, hash_TRIPLE>::iterator common3_it;

#define common3_get(x) (((common3_it=common3.find(x))!=common3.end())?(common3_it->second):0)
#define common2_get(x) (((common2_it=common2.find(x))!=common2.end())?(common2_it->second):0)

int n,m; // n = number of nodes, m = number of edges
int *deg; // degrees of individual nodes
PAIR *edges; // list of edges

int **adj; // adj[x] - adjacency list of node x
PII **inc; // inc[x] - incidence list of node x: (y, edge id)
bool adjacent_list(int x, int y) { return binary_search(adj[x],adj[x]+deg[x],y); }
int *adj_matrix; // compressed adjacency matrix
const int adj_chunk = 8*sizeof(int);
bool adjacent_matrix(int x, int y) { return adj_matrix[(x*n+y)/adj_chunk]&(1<<((x*n+y)%adj_chunk)); }
bool (*adjacent)(int,int);
int getEdgeId(int x, int y) { return inc[x][lower_bound(adj[x],adj[x]+deg[x],y)-adj[x]].second; }

int64 **orbit; // orbit[x][o] - how many times does node x participate in orbit o
int64 **eorbit; // eorbit[x][o] - how many times does node x participate in edge orbit o


/** count graphlets on max 4 nodes */
void count4() {
	clock_t startTime, endTime;
	startTime = clock();
	clock_t startTime_all, endTime_all;
	startTime_all = startTime;
	int frac,frac_prev;

	// precompute triangles that span over edges
	printf("stage 1 - precomputing common nodes\n");
	int *tri = (int*)calloc(m,sizeof(int));
	frac_prev=-1;
	for (int i=0;i<m;i++) {
		frac = 100LL*i/m;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		int x=edges[i].a, y=edges[i].b;
		for (int xi=0,yi=0; xi<deg[x] && yi<deg[y]; ) {
			if (adj[x][xi]==adj[y][yi]) { tri[i]++; xi++; yi++; }
			else if (adj[x][xi]<adj[y][yi]) { xi++; }
			else { yi++; }
		}
	}
	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// count full graphlets
	printf("stage 2 - counting full graphlets\n");
	int64 *C4 = (int64*)calloc(n,sizeof(int64));
	int *neigh = (int*)malloc(n*sizeof(int)), nn;
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		for (int nx=0;nx<deg[x];nx++) {
			int y=adj[x][nx];
			if (y >= x) break;
			nn=0;
			for (int ny=0;ny<deg[y];ny++) {
				int z=adj[y][ny];
				if (z >= y) break;
				if (adjacent(x,z)==0) continue;
				neigh[nn++]=z;
			}
			for (int i=0;i<nn;i++) {
				int z = neigh[i];
				for (int j=i+1;j<nn;j++) {
					int zz = neigh[j];
					if (adjacent(z,zz)) {
						C4[x]++; C4[y]++; C4[z]++; C4[zz]++;
					}
				}
			}
		}
	}
	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// set up a system of equations relating orbits for every node
	printf("stage 3 - building systems of equations\n");
	int *common = (int*)calloc(n,sizeof(int));
	int *common_list = (int*)malloc(n*sizeof(int)), nc=0;
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}

		int64 f_12_14=0, f_10_13=0;
		int64 f_13_14=0, f_11_13=0;
		int64 f_7_11=0, f_5_8=0;
		int64 f_6_9=0, f_9_12=0, f_4_8=0, f_8_12=0;
		int64 f_14=C4[x];

		for (int i=0;i<nc;i++) common[common_list[i]]=0;
		nc=0;

		orbit[x][0]=deg[x];
		// x - middle node
		for (int nx1=0;nx1<deg[x];nx1++) {
			int y=inc[x][nx1].first, ey=inc[x][nx1].second;
			for (int ny=0;ny<deg[y];ny++) {
				int z=inc[y][ny].first, ez=inc[y][ny].second;
				if (adjacent(x,z)) { // triangle
					if (z<y) {
						f_12_14 += tri[ez]-1;
						f_10_13 += (deg[y]-1-tri[ez])+(deg[z]-1-tri[ez]);
					}
				} else {
					if (common[z]==0) common_list[nc++]=z;
					common[z]++;
				}
			}
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int z=inc[x][nx2].first, ez=inc[x][nx2].second;
				if (adjacent(y,z)) { // triangle
					orbit[x][3]++;
					f_13_14 += (tri[ey]-1)+(tri[ez]-1);
					f_11_13 += (deg[x]-1-tri[ey])+(deg[x]-1-tri[ez]);
				} else { // path
					orbit[x][2]++;
					f_7_11 += (deg[x]-1-tri[ey]-1)+(deg[x]-1-tri[ez]-1);
					f_5_8 += (deg[y]-1-tri[ey])+(deg[z]-1-tri[ez]);
				}
			}
		}
		// x - side node
		for (int nx1=0;nx1<deg[x];nx1++) {
			int y=inc[x][nx1].first, ey=inc[x][nx1].second;
			for (int ny=0;ny<deg[y];ny++) {
				int z=inc[y][ny].first, ez=inc[y][ny].second;
				if (x==z) continue;
				if (!adjacent(x,z)) { // path
					orbit[x][1]++;
					f_6_9 += (deg[y]-1-tri[ey]-1);
					f_9_12 += tri[ez];
					f_4_8 += (deg[z]-1-tri[ez]);
					f_8_12 += (common[z]-1);
				}
			}
		}

		// solve system of equations
		orbit[x][14]=(f_14);
		orbit[x][13]=(f_13_14-6*f_14)/2;
		orbit[x][12]=(f_12_14-3*f_14);
		orbit[x][11]=(f_11_13-f_13_14+6*f_14)/2;
		orbit[x][10]=(f_10_13-f_13_14+6*f_14);
		orbit[x][9]=(f_9_12-2*f_12_14+6*f_14)/2;
		orbit[x][8]=(f_8_12-2*f_12_14+6*f_14)/2;
		orbit[x][7]=(f_13_14+f_7_11-f_11_13-6*f_14)/6;
		orbit[x][6]=(2*f_12_14+f_6_9-f_9_12-6*f_14)/2;
		orbit[x][5]=(2*f_12_14+f_5_8-f_8_12-6*f_14);
		orbit[x][4]=(2*f_12_14+f_4_8-f_8_12-6*f_14);
	}

	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);

	endTime_all = endTime;
	printf("total: %.2f\n", (double)(endTime_all-startTime_all)/CLOCKS_PER_SEC);
}


/** count edge orbits of graphlets on max 4 nodes */
void ecount4() {
	clock_t startTime, endTime;
	startTime = clock();
	clock_t startTime_all, endTime_all;
	startTime_all = startTime;
	int frac,frac_prev;

	// precompute triangles that span over edges
	printf("stage 1 - precomputing common nodes\n");
	int *tri = (int*)calloc(m,sizeof(int));
	frac_prev=-1;
	for (int i=0;i<m;i++) {
		frac = 100LL*i/m;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		int x=edges[i].a, y=edges[i].b;
		for (int xi=0,yi=0; xi<deg[x] && yi<deg[y]; ) {
			if (adj[x][xi]==adj[y][yi]) { tri[i]++; xi++; yi++; }
			else if (adj[x][xi]<adj[y][yi]) { xi++; }
			else { yi++; }
		}
	}
	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// count full graphlets
	printf("stage 2 - counting full graphlets\n");
	int64 *C4 = (int64*)calloc(m,sizeof(int64));
	int *neighx = (int*)malloc(n*sizeof(int)); // lookup table - edges to neighbors of x
	memset(neighx,-1,n*sizeof(int));
	int *neigh = (int*)malloc(n*sizeof(int)), nn; // lookup table - common neighbors of x and y
	PII *neigh_edges = (PII*)malloc(n*sizeof(PII)); // list of common neighbors of x and y
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}

		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			neighx[y]=xy;
		}
		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			if (y >= x) break;
			nn=0;
			for (int ny=0;ny<deg[y];ny++) {
				int z=inc[y][ny].first, yz=inc[y][ny].second;
				if (z >= y) break;
				if (neighx[z]==-1) continue;
				int xz=neighx[z];
				neigh[nn]=z;
				neigh_edges[nn]={xz, yz};
				nn++;
			}
			for (int i=0;i<nn;i++) {
				int z = neigh[i], xz = neigh_edges[i].first, yz = neigh_edges[i].second;
				for (int j=i+1;j<nn;j++) {
					int w = neigh[j], xw = neigh_edges[j].first, yw = neigh_edges[j].second;
					if (adjacent(z,w)) {
						C4[xy]++;
						C4[xz]++; C4[yz]++;
						C4[xw]++; C4[yw]++;
						// another iteration to count this last(smallest) edge instead of calling getEdgeId
						//int zw=getEdgeId(z,w); C4[zw]++;
					}
				}
			}
		}
		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			neighx[y]=-1;
		}
	}
	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// count full graphlets for the smallest edge
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		for (int nx=deg[x]-1;nx>=0;nx--) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			if (y <= x) break;
			nn=0;
			for (int ny=deg[y]-1;ny>=0;ny--) {
				int z=adj[y][ny];
				if (z <= y) break;
				if (adjacent(x,z)==0) continue;
				neigh[nn++]=z;
			}
			for (int i=0;i<nn;i++) {
				int z = neigh[i];
				for (int j=i+1;j<nn;j++) {
					int zz = neigh[j];
					if (adjacent(z,zz)) {
						C4[xy]++;
					}
				}
			}
		}
	}
	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// set up a system of equations relating orbits for every node
	printf("stage 3 - building systems of equations\n");
	int *common = (int*)calloc(n,sizeof(int));
	int *common_list = (int*)malloc(n*sizeof(int)), nc=0;
	frac_prev=-1;

	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}

		// common nodes of x and some other node
		for (int i=0;i<nc;i++) common[common_list[i]]=0;
		nc=0;
		for (int nx=0;nx<deg[x];nx++) {
			int y=adj[x][nx];
			for (int ny=0;ny<deg[y];ny++) {
				int z=adj[y][ny];
				if (z==x) continue;
				if (common[z]==0) common_list[nc++]=z;
				common[z]++;
			}
		}

		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			int e=xy;
			for (int n1=0;n1<deg[x];n1++) {
				int z=inc[x][n1].first, xz=inc[x][n1].second;
				if (z==y) continue;
				if (adjacent(y,z)) { // triangle
					if (x<y) {
						eorbit[e][1]++;
						eorbit[e][10] += tri[xy]-1;
						eorbit[e][7] += deg[z]-2;
					}
					eorbit[e][9] += tri[xz]-1;
					eorbit[e][8] += deg[x]-2;
				}
			}
			for (int n1=0;n1<deg[y];n1++) {
				int z=inc[y][n1].first, yz=inc[y][n1].second;
				if (z==x) continue;
				if (!adjacent(x,z)) { // path x-y-z
					eorbit[e][0]++;
					eorbit[e][6] += tri[yz];
					eorbit[e][5] += common[z]-1;
					eorbit[e][4] += deg[y]-2;
					eorbit[e][3] += deg[x]-1;
					eorbit[e][2] += deg[z]-1;
				}
			}
		}
	}
	// solve system of equations
	for (int e=0;e<m;e++) {
		eorbit[e][11]=C4[e];
		eorbit[e][10]=(eorbit[e][10]-2*eorbit[e][11])/2;
		eorbit[e][9]=(eorbit[e][9]-4*eorbit[e][11]);
		eorbit[e][8]=(eorbit[e][8]-eorbit[e][9]-4*eorbit[e][10]-4*eorbit[e][11]);
		eorbit[e][7]=(eorbit[e][7]-eorbit[e][9]-2*eorbit[e][11]);
		eorbit[e][6]=(eorbit[e][6]-eorbit[e][9])/2;
		eorbit[e][5]=(eorbit[e][5]-eorbit[e][9])/2;
		eorbit[e][4]=(eorbit[e][4]-2*eorbit[e][6]-eorbit[e][8]-eorbit[e][9])/2;
		eorbit[e][3]=(eorbit[e][3]-2*eorbit[e][5]-eorbit[e][8]-eorbit[e][9])/2;
		eorbit[e][2]=(eorbit[e][2]-2*eorbit[e][5]-2*eorbit[e][6]-eorbit[e][9]);
	}

	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);

	endTime_all = endTime;
	printf("total: %.2f\n", (double)(endTime_all-startTime_all)/CLOCKS_PER_SEC);
}


/** count graphlets on max 5 nodes */
void count5() {
	clock_t startTime, endTime;
	startTime = clock();
	clock_t startTime_all, endTime_all;
	startTime_all = startTime;
	int frac,frac_prev;

	// precompute common nodes
	printf("stage 1 - precomputing common nodes\n");
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		for (int n1=0;n1<deg[x];n1++) {
			int a=adj[x][n1];
			for (int n2=n1+1;n2<deg[x];n2++) {
				int b=adj[x][n2];
				PAIR ab=PAIR(a,b);
				common2[ab]++;
				for (int n3=n2+1;n3<deg[x];n3++) {
					int c=adj[x][n3];
					int st = adjacent(a,b)+adjacent(a,c)+adjacent(b,c);
					if (st<2) continue;
					TRIPLE abc=TRIPLE(a,b,c);
					common3[abc]++;
				}
			}
		}
	}
	// precompute triangles that span over edges
	int *tri = (int*)calloc(m,sizeof(int));
	for (int i=0;i<m;i++) {
		int x=edges[i].a, y=edges[i].b;
		for (int xi=0,yi=0; xi<deg[x] && yi<deg[y]; ) {
			if (adj[x][xi]==adj[y][yi]) { tri[i]++; xi++; yi++; }
			else if (adj[x][xi]<adj[y][yi]) { xi++; }
			else { yi++; }
		}
	}
	endTime = clock();
	printf("%.2f sec\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// count full graphlets
	printf("stage 2 - counting full graphlets\n");
	int64 *C5 = (int64*)calloc(n,sizeof(int64));
	int *neigh = (int*)malloc(n*sizeof(int)), nn;
	int *neigh2 = (int*)malloc(n*sizeof(int)), nn2;
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		for (int nx=0;nx<deg[x];nx++) {
			int y=adj[x][nx];
			if (y >= x) break;
			nn=0;
			for (int ny=0;ny<deg[y];ny++) {
				int z=adj[y][ny];
				if (z >= y) break;
				if (adjacent(x,z)) {
					neigh[nn++]=z;
				}
			}
			for (int i=0;i<nn;i++) {
				int z = neigh[i];
				nn2=0;
				for (int j=i+1;j<nn;j++) {
					int zz = neigh[j];
					if (adjacent(z,zz)) {
						neigh2[nn2++]=zz;
					}
				}
				for (int i2=0;i2<nn2;i2++) {
					int zz = neigh2[i2];
					for (int j2=i2+1;j2<nn2;j2++) {
						int zzz = neigh2[j2];
						if (adjacent(zz,zzz)) {
							C5[x]++; C5[y]++; C5[z]++; C5[zz]++; C5[zzz]++;
						}
					}
				}
			}
		}
	}
	endTime = clock();
	printf("%.2f sec\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	int *common_x = (int*)calloc(n,sizeof(int));
	int *common_x_list = (int*)malloc(n*sizeof(int)), ncx=0;
	int *common_a = (int*)calloc(n,sizeof(int));
	int *common_a_list = (int*)malloc(n*sizeof(int)), nca=0;

	// set up a system of equations relating orbit counts
	printf("stage 3 - building systems of equations\n");
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}

		for (int i=0;i<ncx;i++) common_x[common_x_list[i]]=0;
		ncx=0;

		// smaller graphlets
		orbit[x][0] = deg[x];
		for (int nx1=0;nx1<deg[x];nx1++) {
			int a=adj[x][nx1];
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int b=adj[x][nx2];
				if (adjacent(a,b)) orbit[x][3]++;
				else orbit[x][2]++;
			}
			for (int na=0;na<deg[a];na++) {
				int b=adj[a][na];
				if (b!=x && !adjacent(x,b)) {
					orbit[x][1]++;
					if (common_x[b]==0) common_x_list[ncx++]=b;
					common_x[b]++;
				}
			}
		}

		int64 f_71=0, f_70=0, f_67=0, f_66=0, f_58=0, f_57=0; // 14
		int64 f_69=0, f_68=0, f_64=0, f_61=0, f_60=0, f_55=0, f_48=0, f_42=0, f_41=0; // 13
		int64 f_65=0, f_63=0, f_59=0, f_54=0, f_47=0, f_46=0, f_40=0; // 12
		int64 f_62=0, f_53=0, f_51=0, f_50=0, f_49=0, f_38=0, f_37=0, f_36=0; // 8
		int64 f_44=0, f_33=0, f_30=0, f_26=0; // 11
		int64 f_52=0, f_43=0, f_32=0, f_29=0, f_25=0; // 10
		int64 f_56=0, f_45=0, f_39=0, f_31=0, f_28=0, f_24=0; // 9
		int64 f_35=0, f_34=0, f_27=0, f_18=0, f_16=0, f_15=0; // 4
		int64 f_17=0; // 5
		int64 f_22=0, f_20=0, f_19=0; // 6
		int64 f_23=0, f_21=0; // 7

		for (int nx1=0;nx1<deg[x];nx1++) {
			int a=inc[x][nx1].first, xa=inc[x][nx1].second;

			for (int i=0;i<nca;i++) common_a[common_a_list[i]]=0;
			nca=0;
			for (int na=0;na<deg[a];na++) {
				int b=adj[a][na];
				for (int nb=0;nb<deg[b];nb++) {
					int c=adj[b][nb];
					if (c==a || adjacent(a,c)) continue;
					if (common_a[c]==0) common_a_list[nca++]=c;
					common_a[c]++;
				}
			}

			// x = orbit-14 (tetrahedron)
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (!adjacent(a,b)) continue;
				for (int nx3=nx2+1;nx3<deg[x];nx3++) {
					int c=inc[x][nx3].first, xc=inc[x][nx3].second;
					if (!adjacent(a,c) || !adjacent(b,c)) continue;
					orbit[x][14]++;
					f_70 += common3_get(TRIPLE(a,b,c))-1;
					f_71 += (tri[xa]>2 && tri[xb]>2)?(common3_get(TRIPLE(x,a,b))-1):0;
					f_71 += (tri[xa]>2 && tri[xc]>2)?(common3_get(TRIPLE(x,a,c))-1):0;
					f_71 += (tri[xb]>2 && tri[xc]>2)?(common3_get(TRIPLE(x,b,c))-1):0;
					f_67 += tri[xa]-2+tri[xb]-2+tri[xc]-2;
					f_66 += common2_get(PAIR(a,b))-2;
					f_66 += common2_get(PAIR(a,c))-2;
					f_66 += common2_get(PAIR(b,c))-2;
					f_58 += deg[x]-3;
					f_57 += deg[a]-3+deg[b]-3+deg[c]-3;
				}
			}

			// x = orbit-13 (diamond)
			for (int nx2=0;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (!adjacent(a,b)) continue;
				for (int nx3=nx2+1;nx3<deg[x];nx3++) {
					int c=inc[x][nx3].first, xc=inc[x][nx3].second;
					if (!adjacent(a,c) || adjacent(b,c)) continue;
					orbit[x][13]++;
					f_69 += (tri[xb]>1 && tri[xc]>1)?(common3_get(TRIPLE(x,b,c))-1):0;
					f_68 += common3_get(TRIPLE(a,b,c))-1;
					f_64 += common2_get(PAIR(b,c))-2;
					f_61 += tri[xb]-1+tri[xc]-1;
					f_60 += common2_get(PAIR(a,b))-1;
					f_60 += common2_get(PAIR(a,c))-1;
					f_55 += tri[xa]-2;
					f_48 += deg[b]-2+deg[c]-2;
					f_42 += deg[x]-3;
					f_41 += deg[a]-3;
				}
			}

			// x = orbit-12 (diamond)
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (!adjacent(a,b)) continue;
				for (int na=0;na<deg[a];na++) {
					int c=inc[a][na].first, ac=inc[a][na].second;
					if (c==x || adjacent(x,c) || !adjacent(b,c)) continue;
					orbit[x][12]++;
					f_65 += (tri[ac]>1)?common3_get(TRIPLE(a,b,c)):0;
					f_63 += common_x[c]-2;
					f_59 += tri[ac]-1+common2_get(PAIR(b,c))-1;
					f_54 += common2_get(PAIR(a,b))-2;
					f_47 += deg[x]-2;
					f_46 += deg[c]-2;
					f_40 += deg[a]-3+deg[b]-3;
				}
			}

			// x = orbit-8 (cycle)
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (adjacent(a,b)) continue;
				for (int na=0;na<deg[a];na++) {
					int c=inc[a][na].first, ac=inc[a][na].second;
					if (c==x || adjacent(x,c) || !adjacent(b,c)) continue;
					orbit[x][8]++;
					f_62 += (tri[ac]>0)?common3_get(TRIPLE(a,b,c)):0;
					f_53 += tri[xa]+tri[xb];
					f_51 += tri[ac]+common2_get(PAIR(c,b));
					f_50 += common_x[c]-2;
					f_49 += common_a[b]-2;
					f_38 += deg[x]-2;
					f_37 += deg[a]-2+deg[b]-2;
					f_36 += deg[c]-2;
				}
			}

			// x = orbit-11 (paw)
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (!adjacent(a,b)) continue;
				for (int nx3=0;nx3<deg[x];nx3++) {
					int c=inc[x][nx3].first, xc=inc[x][nx3].second;
					if (c==a || c==b || adjacent(a,c) || adjacent(b,c)) continue;
					orbit[x][11]++;
					f_44 += tri[xc];
					f_33 += deg[x]-3;
					f_30 += deg[c]-1;
					f_26 += deg[a]-2+deg[b]-2;
				}
			}

			// x = orbit-10 (paw)
			for (int nx2=0;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (!adjacent(a,b)) continue;
				for (int nb=0;nb<deg[b];nb++) {
					int c=inc[b][nb].first, bc=inc[b][nb].second;
					if (c==x || c==a || adjacent(a,c) || adjacent(x,c)) continue;
					orbit[x][10]++;
					f_52 += common_a[c]-1;
					f_43 += tri[bc];
					f_32 += deg[b]-3;
					f_29 += deg[c]-1;
					f_25 += deg[a]-2;
				}
			}

			// x = orbit-9 (paw)
			for (int na1=0;na1<deg[a];na1++) {
				int b=inc[a][na1].first, ab=inc[a][na1].second;
				if (b==x || adjacent(x,b)) continue;
				for (int na2=na1+1;na2<deg[a];na2++) {
					int c=inc[a][na2].first, ac=inc[a][na2].second;
					if (c==x || !adjacent(b,c) || adjacent(x,c)) continue;
					orbit[x][9]++;
					f_56 += (tri[ab]>1 && tri[ac]>1)?common3_get(TRIPLE(a,b,c)):0;
					f_45 += common2_get(PAIR(b,c))-1;
					f_39 += tri[ab]-1+tri[ac]-1;
					f_31 += deg[a]-3;
					f_28 += deg[x]-1;
					f_24 += deg[b]-2+deg[c]-2;
				}
			}

			// x = orbit-4 (path)
			for (int na=0;na<deg[a];na++) {
				int b=inc[a][na].first, ab=inc[a][na].second;
				if (b==x || adjacent(x,b)) continue;
				for (int nb=0;nb<deg[b];nb++) {
					int c=inc[b][nb].first, bc=inc[b][nb].second;
					if (c==a || adjacent(a,c) || adjacent(x,c)) continue;
					orbit[x][4]++;
					f_35 += common_a[c]-1;
					f_34 += common_x[c];
					f_27 += tri[bc];
					f_18 += deg[b]-2;
					f_16 += deg[x]-1;
					f_15 += deg[c]-1;
				}
			}

			// x = orbit-5 (path)
			for (int nx2=0;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (b==a || adjacent(a,b)) continue;
				for (int nb=0;nb<deg[b];nb++) {
					int c=inc[b][nb].first, bc=inc[b][nb].second;
					if (c==x || adjacent(a,c) || adjacent(x,c)) continue;
					orbit[x][5]++;
					f_17 += deg[a]-1;
				}
			}

			// x = orbit-6 (claw)
			for (int na1=0;na1<deg[a];na1++) {
				int b=inc[a][na1].first, ab=inc[a][na1].second;
				if (b==x || adjacent(x,b)) continue;
				for (int na2=na1+1;na2<deg[a];na2++) {
					int c=inc[a][na2].first, ac=inc[a][na2].second;
					if (c==x || adjacent(x,c) || adjacent(b,c)) continue;
					orbit[x][6]++;
					f_22 += deg[a]-3;
					f_20 += deg[x]-1;
					f_19 += deg[b]-1+deg[c]-1;
				}
			}

			// x = orbit-7 (claw)
			for (int nx2=nx1+1;nx2<deg[x];nx2++) {
				int b=inc[x][nx2].first, xb=inc[x][nx2].second;
				if (adjacent(a,b)) continue;
				for (int nx3=nx2+1;nx3<deg[x];nx3++) {
					int c=inc[x][nx3].first, xc=inc[x][nx3].second;
					if (adjacent(a,c) || adjacent(b,c)) continue;
					orbit[x][7]++;
					f_23 += deg[x]-3;
					f_21 += deg[a]-1+deg[b]-1+deg[c]-1;
				}
			}
		}

		// solve equations
		orbit[x][72] = C5[x];
		orbit[x][71] = (f_71-12*orbit[x][72])/2;
		orbit[x][70] = (f_70-4*orbit[x][72]);
		orbit[x][69] = (f_69-2*orbit[x][71])/4;
		orbit[x][68] = (f_68-2*orbit[x][71]);
		orbit[x][67] = (f_67-12*orbit[x][72]-4*orbit[x][71]);
		orbit[x][66] = (f_66-12*orbit[x][72]-2*orbit[x][71]-3*orbit[x][70]);
		orbit[x][65] = (f_65-3*orbit[x][70])/2;
		orbit[x][64] = (f_64-2*orbit[x][71]-4*orbit[x][69]-1*orbit[x][68]);
		orbit[x][63] = (f_63-3*orbit[x][70]-2*orbit[x][68]);
		orbit[x][62] = (f_62-1*orbit[x][68])/2;
		orbit[x][61] = (f_61-4*orbit[x][71]-8*orbit[x][69]-2*orbit[x][67])/2;
		orbit[x][60] = (f_60-4*orbit[x][71]-2*orbit[x][68]-2*orbit[x][67]);
		orbit[x][59] = (f_59-6*orbit[x][70]-2*orbit[x][68]-4*orbit[x][65]);
		orbit[x][58] = (f_58-4*orbit[x][72]-2*orbit[x][71]-1*orbit[x][67]);
		orbit[x][57] = (f_57-12*orbit[x][72]-4*orbit[x][71]-3*orbit[x][70]-1*orbit[x][67]-2*orbit[x][66]);
		orbit[x][56] = (f_56-2*orbit[x][65])/3;
		orbit[x][55] = (f_55-2*orbit[x][71]-2*orbit[x][67])/3;
		orbit[x][54] = (f_54-3*orbit[x][70]-1*orbit[x][66]-2*orbit[x][65])/2;
		orbit[x][53] = (f_53-2*orbit[x][68]-2*orbit[x][64]-2*orbit[x][63]);
		orbit[x][52] = (f_52-2*orbit[x][66]-2*orbit[x][64]-1*orbit[x][59])/2;
		orbit[x][51] = (f_51-2*orbit[x][68]-2*orbit[x][63]-4*orbit[x][62]);
		orbit[x][50] = (f_50-1*orbit[x][68]-2*orbit[x][63])/3;
		orbit[x][49] = (f_49-1*orbit[x][68]-1*orbit[x][64]-2*orbit[x][62])/2;
		orbit[x][48] = (f_48-4*orbit[x][71]-8*orbit[x][69]-2*orbit[x][68]-2*orbit[x][67]-2*orbit[x][64]-2*orbit[x][61]-1*orbit[x][60]);
		orbit[x][47] = (f_47-3*orbit[x][70]-2*orbit[x][68]-1*orbit[x][66]-1*orbit[x][63]-1*orbit[x][60]);
		orbit[x][46] = (f_46-3*orbit[x][70]-2*orbit[x][68]-2*orbit[x][65]-1*orbit[x][63]-1*orbit[x][59]);
		orbit[x][45] = (f_45-2*orbit[x][65]-2*orbit[x][62]-3*orbit[x][56]);
		orbit[x][44] = (f_44-1*orbit[x][67]-2*orbit[x][61])/4;
		orbit[x][43] = (f_43-2*orbit[x][66]-1*orbit[x][60]-1*orbit[x][59])/2;
		orbit[x][42] = (f_42-2*orbit[x][71]-4*orbit[x][69]-2*orbit[x][67]-2*orbit[x][61]-3*orbit[x][55]);
		orbit[x][41] = (f_41-2*orbit[x][71]-1*orbit[x][68]-2*orbit[x][67]-1*orbit[x][60]-3*orbit[x][55]);
		orbit[x][40] = (f_40-6*orbit[x][70]-2*orbit[x][68]-2*orbit[x][66]-4*orbit[x][65]-1*orbit[x][60]-1*orbit[x][59]-4*orbit[x][54]);
		orbit[x][39] = (f_39-4*orbit[x][65]-1*orbit[x][59]-6*orbit[x][56])/2;
		orbit[x][38] = (f_38-1*orbit[x][68]-1*orbit[x][64]-2*orbit[x][63]-1*orbit[x][53]-3*orbit[x][50]);
		orbit[x][37] = (f_37-2*orbit[x][68]-2*orbit[x][64]-2*orbit[x][63]-4*orbit[x][62]-1*orbit[x][53]-1*orbit[x][51]-4*orbit[x][49]);
		orbit[x][36] = (f_36-1*orbit[x][68]-2*orbit[x][63]-2*orbit[x][62]-1*orbit[x][51]-3*orbit[x][50]);
		orbit[x][35] = (f_35-1*orbit[x][59]-2*orbit[x][52]-2*orbit[x][45])/2;
		orbit[x][34] = (f_34-1*orbit[x][59]-2*orbit[x][52]-1*orbit[x][51])/2;
		orbit[x][33] = (f_33-1*orbit[x][67]-2*orbit[x][61]-3*orbit[x][58]-4*orbit[x][44]-2*orbit[x][42])/2;
		orbit[x][32] = (f_32-2*orbit[x][66]-1*orbit[x][60]-1*orbit[x][59]-2*orbit[x][57]-2*orbit[x][43]-2*orbit[x][41]-1*orbit[x][40])/2;
		orbit[x][31] = (f_31-2*orbit[x][65]-1*orbit[x][59]-3*orbit[x][56]-1*orbit[x][43]-2*orbit[x][39]);
		orbit[x][30] = (f_30-1*orbit[x][67]-1*orbit[x][63]-2*orbit[x][61]-1*orbit[x][53]-4*orbit[x][44]);
		orbit[x][29] = (f_29-2*orbit[x][66]-2*orbit[x][64]-1*orbit[x][60]-1*orbit[x][59]-1*orbit[x][53]-2*orbit[x][52]-2*orbit[x][43]);
		orbit[x][28] = (f_28-2*orbit[x][65]-2*orbit[x][62]-1*orbit[x][59]-1*orbit[x][51]-1*orbit[x][43]);
		orbit[x][27] = (f_27-1*orbit[x][59]-1*orbit[x][51]-2*orbit[x][45])/2;
		orbit[x][26] = (f_26-2*orbit[x][67]-2*orbit[x][63]-2*orbit[x][61]-6*orbit[x][58]-1*orbit[x][53]-2*orbit[x][47]-2*orbit[x][42]);
		orbit[x][25] = (f_25-2*orbit[x][66]-2*orbit[x][64]-1*orbit[x][59]-2*orbit[x][57]-2*orbit[x][52]-1*orbit[x][48]-1*orbit[x][40])/2;
		orbit[x][24] = (f_24-4*orbit[x][65]-4*orbit[x][62]-1*orbit[x][59]-6*orbit[x][56]-1*orbit[x][51]-2*orbit[x][45]-2*orbit[x][39]);
		orbit[x][23] = (f_23-1*orbit[x][55]-1*orbit[x][42]-2*orbit[x][33])/4;
		orbit[x][22] = (f_22-2*orbit[x][54]-1*orbit[x][40]-1*orbit[x][39]-1*orbit[x][32]-2*orbit[x][31])/3;
		orbit[x][21] = (f_21-3*orbit[x][55]-3*orbit[x][50]-2*orbit[x][42]-2*orbit[x][38]-2*orbit[x][33]);
		orbit[x][20] = (f_20-2*orbit[x][54]-2*orbit[x][49]-1*orbit[x][40]-1*orbit[x][37]-1*orbit[x][32]);
		orbit[x][19] = (f_19-4*orbit[x][54]-4*orbit[x][49]-1*orbit[x][40]-2*orbit[x][39]-1*orbit[x][37]-2*orbit[x][35]-2*orbit[x][31]);
		orbit[x][18] = (f_18-1*orbit[x][59]-1*orbit[x][51]-2*orbit[x][46]-2*orbit[x][45]-2*orbit[x][36]-2*orbit[x][27]-1*orbit[x][24])/2;
		orbit[x][17] = (f_17-1*orbit[x][60]-1*orbit[x][53]-1*orbit[x][51]-1*orbit[x][48]-1*orbit[x][37]-2*orbit[x][34]-2*orbit[x][30])/2;
		orbit[x][16] = (f_16-1*orbit[x][59]-2*orbit[x][52]-1*orbit[x][51]-2*orbit[x][46]-2*orbit[x][36]-2*orbit[x][34]-1*orbit[x][29]);
		orbit[x][15] = (f_15-1*orbit[x][59]-2*orbit[x][52]-1*orbit[x][51]-2*orbit[x][45]-2*orbit[x][35]-2*orbit[x][34]-2*orbit[x][27]);
	}
	endTime = clock();
	printf("%.2f sec\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);

	endTime_all = endTime;
	printf("total: %.2f sec\n", (double)(endTime_all-startTime_all)/CLOCKS_PER_SEC);
}


/** count edge orbits of graphlets on max 5 nodes */
void ecount5() {
	clock_t startTime, endTime;
	startTime = clock();
	clock_t startTime_all, endTime_all;
	startTime_all = startTime;
	int frac,frac_prev;

	// precompute common nodes
	printf("stage 1 - precomputing common nodes\n");
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}
		for (int n1=0;n1<deg[x];n1++) {
			int a=adj[x][n1];
			for (int n2=n1+1;n2<deg[x];n2++) {
				int b=adj[x][n2];
				PAIR ab=PAIR(a,b);
				common2[ab]++;
				for (int n3=n2+1;n3<deg[x];n3++) {
					int c=adj[x][n3];
					int st = adjacent(a,b)+adjacent(a,c)+adjacent(b,c);
					if (st<2) continue;
					TRIPLE abc=TRIPLE(a,b,c);
					common3[abc]++;
				}
			}
		}
	}
	// precompute triangles that span over edges
	int *tri = (int*)calloc(m,sizeof(int));
	for (int i=0;i<m;i++) {
		int x=edges[i].a, y=edges[i].b;
		for (int xi=0,yi=0; xi<deg[x] && yi<deg[y]; ) {
			if (adj[x][xi]==adj[y][yi]) { tri[i]++; xi++; yi++; }
			else if (adj[x][xi]<adj[y][yi]) { xi++; }
			else { yi++; }
		}
	}
	endTime = clock();
	printf("%.2f sec\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// count full graphlets
	printf("stage 2 - counting full graphlets\n");
	int64 *C5 = (int64*)calloc(m,sizeof(int64));
	int *neighx = (int*)malloc(n*sizeof(int)); // lookup table - edges to neighbors of x
	memset(neighx,-1,n*sizeof(int));
	int *neigh = (int*)malloc(n*sizeof(int)), nn; // lookup table - common neighbors of x and y
	PII *neigh_edges = (PII*)malloc(n*sizeof(PII)); // list of common neighbors of x and y
	int *neigh2 = (int*)malloc(n*sizeof(int)), nn2;
	TIII *neigh2_edges = (TIII*)malloc(n*sizeof(TIII));
	frac_prev=-1;
	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}

		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			neighx[y]=xy;
		}
		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			if (y >= x) break;
			nn=0;
			for (int ny=0;ny<deg[y];ny++) {
				int z=inc[y][ny].first, yz=inc[y][ny].second;
				if (z >= y) break;
				if (neighx[z]==-1) continue;
				int xz=neighx[z];
				neigh[nn]=z;
				neigh_edges[nn]={xz, yz};
				nn++;
			}
			for (int i=0;i<nn;i++) {
				int z = neigh[i], xz = neigh_edges[i].first, yz = neigh_edges[i].second;
				nn2 = 0;
				for (int j=i+1;j<nn;j++) {
					int w = neigh[j], xw = neigh_edges[j].first, yw = neigh_edges[j].second;
					if (adjacent(z,w)) {
						neigh2[nn2]=w;
						int zw=getEdgeId(z,w);
						neigh2_edges[nn2]={xw,yw,zw};
						nn2++;
					}
				}
				for (int i2=0;i2<nn2;i2++) {
					int z2 = neigh2[i2];
					int z2x=neigh2_edges[i2].first, z2y=neigh2_edges[i2].second, z2z=neigh2_edges[i2].third;
					for (int j2=i2+1;j2<nn2;j2++) {
						int z3 = neigh2[j2];
						int z3x=neigh2_edges[j2].first, z3y=neigh2_edges[j2].second, z3z=neigh2_edges[j2].third;
						if (adjacent(z2,z3)) {
							int zid=getEdgeId(z2,z3);
							C5[xy]++; C5[xz]++; C5[yz]++;
							C5[z2x]++; C5[z2y]++; C5[z2z]++;
							C5[z3x]++; C5[z3y]++; C5[z3z]++;
							C5[zid]++;
						}
					}
				}
			}
		}
		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			neighx[y]=-1;
		}
	}
	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);
	startTime = endTime;

	// set up a system of equations relating orbits for every node
	printf("stage 3 - building systems of equations\n");
	int *common_x = (int*)calloc(n,sizeof(int));
	int *common_x_list = (int*)malloc(n*sizeof(int)), nc_x=0;
	int *common_y = (int*)calloc(n,sizeof(int));
	int *common_y_list = (int*)malloc(n*sizeof(int)), nc_y=0;
	frac_prev=-1;

	for (int x=0;x<n;x++) {
		frac = 100LL*x/n;
		if (frac!=frac_prev) {
			printf("%d%%\r",frac);
			frac_prev=frac;
		}

		// common nodes of x and some other node
		for (int i=0;i<nc_x;i++) common_x[common_x_list[i]]=0;
		nc_x=0;
		for (int nx=0;nx<deg[x];nx++) {
			int a=adj[x][nx];
			for (int na=0;na<deg[a];na++) {
				int z=adj[a][na];
				if (z==x) continue;
				if (common_x[z]==0) common_x_list[nc_x++]=z;
				common_x[z]++;
			}
		}

		for (int nx=0;nx<deg[x];nx++) {
			int y=inc[x][nx].first, xy=inc[x][nx].second;
			int e=xy;
			if (y>=x) break;

			// common nodes of y and some other node
			for (int i=0;i<nc_y;i++) common_y[common_y_list[i]]=0;
			nc_y=0;
			for (int ny=0;ny<deg[y];ny++) {
				int a=adj[y][ny];
				for (int na=0;na<deg[a];na++) {
					int z=adj[a][na];
					if (z==y) continue;
					if (common_y[z]==0) common_y_list[nc_y++]=z;
					common_y[z]++;
				}
			}

			int64 f_66=0, f_65=0, f_62=0, f_61=0, f_60=0, f_51=0, f_50=0; // 11
			int64 f_64=0, f_58=0, f_55=0, f_48=0, f_41=0, f_35=0; // 10
			int64 f_63=0, f_59=0, f_57=0, f_54=0, f_53=0, f_52=0, f_47=0, f_40=0, f_39=0, f_34=0, f_33=0; // 9
			int64 f_45=0, f_36=0, f_26=0, f_23=0, f_19=0; // 7
			int64 f_49=0, f_38=0, f_37=0, f_32=0, f_25=0, f_22=0, f_18=0; // 6
			int64 f_56=0, f_46=0, f_44=0, f_43=0, f_42=0, f_31=0, f_30=0; // 5
			int64 f_27=0, f_17=0, f_15=0; // 4
			int64 f_20=0, f_16=0, f_13=0; // 3
			int64 f_29=0, f_28=0, f_24=0, f_21=0, f_14=0, f_12=0; // 2

			// smaller (3-node) graphlets
			orbit[x][0] = deg[x];
			for (int nx1=0;nx1<deg[x];nx1++) {
				int z=adj[x][nx1];
				if (z==y) continue;
				if (adjacent(y,z)) eorbit[e][1]++;
				else eorbit[e][0]++;
			}
			for (int ny=0;ny<deg[y];ny++) {
				int z=adj[y][ny];
				if (z==x) continue;
				if (!adjacent(x,z)) eorbit[e][0]++;
			}

			// edge-orbit 11 = (14,14)
			for (int nx1=0;nx1<deg[x];nx1++) {
				int a=adj[x][nx1], xa=inc[x][nx1].second;
				if (a==y || !adjacent(y,a)) continue;
				for (int nx2=nx1+1;nx2<deg[x];nx2++) {
					int b=adj[x][nx2], xb=inc[x][nx2].second;
					if (b==y || !adjacent(y,b) || !adjacent(a,b)) continue;
					int ya=getEdgeId(y,a), yb=getEdgeId(y,b), ab=getEdgeId(a,b);
					eorbit[e][11]++;
					f_66 += common3_get(TRIPLE(x,y,a))-1;
					f_66 += common3_get(TRIPLE(x,y,b))-1;
					f_65 += common3_get(TRIPLE(a,b,x))-1;
					f_65 += common3_get(TRIPLE(a,b,y))-1;
					f_62 += tri[xy]-2;
					f_61 += (tri[xa]-2)+(tri[xb]-2)+(tri[ya]-2)+(tri[yb]-2);
					f_60 += tri[ab]-2;
					f_51 += (deg[x]-3)+(deg[y]-3);
					f_50 += (deg[a]-3)+(deg[b]-3);
				}
			}

			// edge-orbit 10 = (13,13)
			for (int nx1=0;nx1<deg[x];nx1++) {
				int a=adj[x][nx1], xa=inc[x][nx1].second;
				if (a==y || !adjacent(y,a)) continue;
				for (int nx2=nx1+1;nx2<deg[x];nx2++) {
					int b=adj[x][nx2], xb=inc[x][nx2].second;
					if (b==y || !adjacent(y,b) || adjacent(a,b)) continue;
					int ya=getEdgeId(y,a), yb=getEdgeId(y,b);
					eorbit[e][10]++;
					f_64 += common3_get(TRIPLE(a,b,x))-1;
					f_64 += common3_get(TRIPLE(a,b,y))-1;
					f_58 += common2_get(PAIR(a,b))-2;
					f_55 += (tri[xa]-1)+(tri[xb]-1)+(tri[ya]-1)+(tri[yb]-1);
					f_48 += tri[xy]-2;
					f_41 += (deg[a]-2)+(deg[b]-2);
					f_35 += (deg[x]-3)+(deg[y]-3);
				}
			}

			// edge-orbit 9 = (12,13)
			for (int nx=0;nx<deg[x];nx++) {
				int a=adj[x][nx], xa=inc[x][nx].second;
				if (a==y) continue;
				for (int ny=0;ny<deg[y];ny++) {
					int b=adj[y][ny], yb=inc[y][ny].second;
					if (b==x || !adjacent(a,b)) continue;
					int adj_ya=adjacent(y,a), adj_xb=adjacent(x,b);
					if (adj_ya+adj_xb!=1) continue;
					int ab=getEdgeId(a,b);
					eorbit[e][9]++;
					if (adj_xb) {
						int xb=getEdgeId(x,b);
						f_63 += common3_get(TRIPLE(a,b,y))-1;
						f_59 += common3_get(TRIPLE(a,b,x));
						f_57 += common_y[a]-2;
						f_54 += tri[yb]-1;
						f_53 += tri[xa]-1;
						f_47 += tri[xb]-2;
						f_40 += deg[y]-2;
						f_39 += deg[a]-2;
						f_34 += deg[x]-3;
						f_33 += deg[b]-3;
					} else if (adj_ya) {
						int ya=getEdgeId(y,a);
						f_63 += common3_get(TRIPLE(a,b,x))-1;
						f_59 += common3_get(TRIPLE(a,b,y));
						f_57 += common_x[b]-2;
						f_54 += tri[xa]-1;
						f_53 += tri[yb]-1;
						f_47 += tri[ya]-2;
						f_40 += deg[x]-2;
						f_39 += deg[b]-2;
						f_34 += deg[y]-3;
						f_33 += deg[a]-3;
					}
					f_52 += tri[ab]-1;
				}
			}

			// edge-orbit 8 = (10,11)
			for (int nx=0;nx<deg[x];nx++) {
				int a=adj[x][nx];
				if (a==y || !adjacent(y,a)) continue;
				for (int nx1=0;nx1<deg[x];nx1++) {
					int b=adj[x][nx1];
					if (b==y || b==a || adjacent(y,b) || adjacent(a,b)) continue;
					eorbit[e][8]++;
				}
				for (int ny1=0;ny1<deg[y];ny1++) {
					int b=adj[y][ny1];
					if (b==x || b==a || adjacent(x,b) || adjacent(a,b)) continue;
					eorbit[e][8]++;
				}
			}

			// edge-orbit 7 = (10,10)
			for (int nx=0;nx<deg[x];nx++) {
				int a=adj[x][nx];
				if (a==y || !adjacent(y,a)) continue;
				for (int na=0;na<deg[a];na++) {
					int b=adj[a][na], ab=inc[a][na].second;
					if (b==x || b==y || adjacent(x,b) || adjacent(y,b)) continue;
					eorbit[e][7]++;
					f_45 += common_x[b]-1;
					f_45 += common_y[b]-1;
					f_36 += tri[ab];
					f_26 += deg[a]-3;
					f_23 += deg[b]-1;
					f_19 += (deg[x]-2)+(deg[y]-2);
				}
			}

			// edge-orbit 6 = (9,11)
			for (int ny1=0;ny1<deg[y];ny1++) {
				int a=adj[y][ny1], ya=inc[y][ny1].second;
				if (a==x || adjacent(x,a)) continue;
				for (int ny2=ny1+1;ny2<deg[y];ny2++) {
					int b=adj[y][ny2], yb=inc[y][ny2].second;
					if (b==x || adjacent(x,b) || !adjacent(a,b)) continue;
					int ab=getEdgeId(a,b);
					eorbit[e][6]++;
					f_49 += common3_get(TRIPLE(y,a,b));
					f_38 += tri[ab]-1;
					f_37 += tri[xy];
					f_32 += (tri[ya]-1)+(tri[yb]-1);
					f_25 += deg[y]-3;
					f_22 += deg[x]-1;
					f_18 += (deg[a]-2)+(deg[b]-2);
				}
			}
			for (int nx1=0;nx1<deg[x];nx1++) {
				int a=adj[x][nx1], xa=inc[x][nx1].second;
				if (a==y || adjacent(y,a)) continue;
				for (int nx2=nx1+1;nx2<deg[x];nx2++) {
					int b=adj[x][nx2], xb=inc[x][nx2].second;
					if (b==y || adjacent(y,b) || !adjacent(a,b)) continue;
					int ab=getEdgeId(a,b);
					eorbit[e][6]++;
					f_49 += common3_get(TRIPLE(x,a,b));
					f_38 += tri[ab]-1;
					f_37 += tri[xy];
					f_32 += (tri[xa]-1)+(tri[xb]-1);
					f_25 += deg[x]-3;
					f_22 += deg[y]-1;
					f_18 += (deg[a]-2)+(deg[b]-2);
				}
			}

			// edge-orbit 5 = (8,8)
			for (int nx=0;nx<deg[x];nx++) {
				int a=adj[x][nx], xa=inc[x][nx].second;
				if (a==y || adjacent(y,a)) continue;
				for (int ny=0;ny<deg[y];ny++) {
					int b=adj[y][ny], yb=inc[y][ny].second;
					if (b==x || adjacent(x,b) || !adjacent(a,b)) continue;
					int ab=getEdgeId(a,b);
					eorbit[e][5]++;
					f_56 += common3_get(TRIPLE(x,a,b));
					f_56 += common3_get(TRIPLE(y,a,b));
					f_46 += tri[xy];
					f_44 += tri[xa]+tri[yb];
					f_43 += tri[ab];
					f_42 += common_x[b]-2;
					f_42 += common_y[a]-2;
					f_31 += (deg[x]-2)+(deg[y]-2);
					f_30 += (deg[a]-2)+(deg[b]-2);
				}
			}

			// edge-orbit 4 = (6,7)
			for (int ny1=0;ny1<deg[y];ny1++) {
				int a=adj[y][ny1];
				if (a==x || adjacent(x,a)) continue;
				for (int ny2=ny1+1;ny2<deg[y];ny2++) {
					int b=adj[y][ny2];
					if (b==x || adjacent(x,b) || adjacent(a,b)) continue;
					eorbit[e][4]++;
					f_27 += tri[xy];
					f_17 += deg[y]-3;
					f_15 += (deg[a]-1)+(deg[b]-1);
				}
			}
			for (int nx1=0;nx1<deg[x];nx1++) {
				int a=adj[x][nx1];
				if (a==y || adjacent(y,a)) continue;
				for (int nx2=nx1+1;nx2<deg[x];nx2++) {
					int b=adj[x][nx2];
					if (b==y || adjacent(y,b) || adjacent(a,b)) continue;
					eorbit[e][4]++;
					f_27 += tri[xy];
					f_17 += deg[x]-3;
					f_15 += (deg[a]-1)+(deg[b]-1);
				}
			}

			// edge-orbit 3 = (5,5)
			for (int nx=0;nx<deg[x];nx++) {
				int a=adj[x][nx];
				if (a==y || adjacent(y,a)) continue;
				for (int ny=0;ny<deg[y];ny++) {
					int b=adj[y][ny];
					if (b==x || adjacent(x,b) || adjacent(a,b)) continue;
					eorbit[e][3]++;
					f_20 += tri[xy];
					f_16 += (deg[x]-2)+(deg[y]-2);
					f_13 += (deg[a]-1)+(deg[b]-1);
				}
			}

			// edge-orbit 2 = (4,5)
			for (int ny=0;ny<deg[y];ny++) {
				int a=adj[y][ny];
				if (a==x || adjacent(x,a)) continue;
				for (int na=0;na<deg[a];na++) {
					int b=adj[a][na], ab=inc[a][na].second;
					if (b==y || adjacent(y,b) || adjacent(x,b)) continue;
					eorbit[e][2]++;
					f_29 += common_y[b]-1;
					f_28 += common_x[b];
					f_24 += tri[xy];
					f_21 += tri[ab];
					f_14 += deg[a]-2;
					f_12 += deg[b]-1;
				}
			}
			for (int nx=0;nx<deg[x];nx++) {
				int a=adj[x][nx];
				if (a==y || adjacent(y,a)) continue;
				for (int na=0;na<deg[a];na++) {
					int b=adj[a][na], ab=inc[a][na].second;
					if (b==x || adjacent(x,b) || adjacent(y,b)) continue;
					eorbit[e][2]++;
					f_29 += common_x[b]-1;
					f_28 += common_y[b];
					f_24 += tri[xy];
					f_21 += tri[ab];
					f_14 += deg[a]-2;
					f_12 += deg[b]-1;
				}
			}

			// solve system of equations
			eorbit[e][67]=C5[e];
			eorbit[e][66]=(f_66-6*eorbit[e][67])/2;
			eorbit[e][65]=(f_65-6*eorbit[e][67]);
			eorbit[e][64]=(f_64-2*eorbit[e][66]);
			eorbit[e][63]=(f_63-2*eorbit[e][65])/2;
			eorbit[e][62]=(f_62-2*eorbit[e][66]-3*eorbit[e][67]);
			eorbit[e][61]=(f_61-2*eorbit[e][65]-4*eorbit[e][66]-12*eorbit[e][67]);
			eorbit[e][60]=(f_60-1*eorbit[e][65]-3*eorbit[e][67]);
			eorbit[e][59]=(f_59-2*eorbit[e][65])/2;
			eorbit[e][58]=(f_58-1*eorbit[e][64]-1*eorbit[e][66]);
			eorbit[e][57]=(f_57-2*eorbit[e][63]-2*eorbit[e][64]-2*eorbit[e][65]);
			eorbit[e][56]=(f_56-2*eorbit[e][63])/2;
			eorbit[e][55]=(f_55-4*eorbit[e][62]-2*eorbit[e][64]-4*eorbit[e][66]);
			eorbit[e][54]=(f_54-1*eorbit[e][61]-2*eorbit[e][63]-2*eorbit[e][65])/2;
			eorbit[e][53]=(f_53-2*eorbit[e][59]-2*eorbit[e][64]-2*eorbit[e][65]);
			eorbit[e][52]=(f_52-2*eorbit[e][59]-2*eorbit[e][63]-2*eorbit[e][65]);
			eorbit[e][51]=(f_51-1*eorbit[e][61]-2*eorbit[e][62]-1*eorbit[e][65]-4*eorbit[e][66]-6*eorbit[e][67]);
			eorbit[e][50]=(f_50-2*eorbit[e][60]-1*eorbit[e][61]-2*eorbit[e][65]-2*eorbit[e][66]-6*eorbit[e][67]);
			eorbit[e][49]=(f_49-1*eorbit[e][59])/3;
			eorbit[e][48]=(f_48-2*eorbit[e][62]-1*eorbit[e][66])/3;
			eorbit[e][47]=(f_47-2*eorbit[e][59]-1*eorbit[e][61]-2*eorbit[e][65])/2;
			eorbit[e][46]=(f_46-1*eorbit[e][57]-1*eorbit[e][63]);
			eorbit[e][45]=(f_45-1*eorbit[e][52]-4*eorbit[e][58]-4*eorbit[e][60]);
			eorbit[e][44]=(f_44-2*eorbit[e][56]-1*eorbit[e][57]-2*eorbit[e][63]);
			eorbit[e][43]=(f_43-2*eorbit[e][56]-1*eorbit[e][63]);
			eorbit[e][42]=(f_42-2*eorbit[e][56]-1*eorbit[e][57]-2*eorbit[e][63])/2;
			eorbit[e][41]=(f_41-1*eorbit[e][55]-2*eorbit[e][58]-2*eorbit[e][62]-2*eorbit[e][64]-2*eorbit[e][66]);
			eorbit[e][40]=(f_40-2*eorbit[e][54]-1*eorbit[e][55]-1*eorbit[e][57]-1*eorbit[e][61]-2*eorbit[e][63]-2*eorbit[e][64]-2*eorbit[e][65]);
			eorbit[e][39]=(f_39-1*eorbit[e][52]-1*eorbit[e][53]-1*eorbit[e][57]-2*eorbit[e][59]-2*eorbit[e][63]-2*eorbit[e][64]-2*eorbit[e][65]);
			eorbit[e][38]=(f_38-3*eorbit[e][49]-1*eorbit[e][56]-1*eorbit[e][59]);
			eorbit[e][37]=(f_37-1*eorbit[e][53]-1*eorbit[e][59]);
			eorbit[e][36]=(f_36-1*eorbit[e][52]-2*eorbit[e][60])/2;
			eorbit[e][35]=(f_35-6*eorbit[e][48]-1*eorbit[e][55]-4*eorbit[e][62]-1*eorbit[e][64]-2*eorbit[e][66]);
			eorbit[e][34]=(f_34-2*eorbit[e][47]-1*eorbit[e][53]-1*eorbit[e][55]-2*eorbit[e][59]-1*eorbit[e][61]-2*eorbit[e][64]-2*eorbit[e][65]);
			eorbit[e][33]=(f_33-2*eorbit[e][47]-1*eorbit[e][52]-2*eorbit[e][54]-2*eorbit[e][59]-1*eorbit[e][61]-2*eorbit[e][63]-2*eorbit[e][65]);
			eorbit[e][32]=(f_32-6*eorbit[e][49]-1*eorbit[e][53]-2*eorbit[e][59])/2;
			eorbit[e][31]=(f_31-2*eorbit[e][42]-1*eorbit[e][44]-2*eorbit[e][46]-2*eorbit[e][56]-2*eorbit[e][57]-2*eorbit[e][63]);
			eorbit[e][30]=(f_30-2*eorbit[e][42]-2*eorbit[e][43]-1*eorbit[e][44]-4*eorbit[e][56]-1*eorbit[e][57]-2*eorbit[e][63]);
			eorbit[e][29]=(f_29-2*eorbit[e][38]-1*eorbit[e][45]-1*eorbit[e][52])/2;
			eorbit[e][28]=(f_28-2*eorbit[e][43]-1*eorbit[e][45]-1*eorbit[e][52])/2;
			eorbit[e][27]=(f_27-1*eorbit[e][34]-1*eorbit[e][47]);
			eorbit[e][26]=(f_26-1*eorbit[e][33]-2*eorbit[e][36]-1*eorbit[e][50]-1*eorbit[e][52]-2*eorbit[e][60])/2;
			eorbit[e][25]=(f_25-2*eorbit[e][32]-1*eorbit[e][37]-3*eorbit[e][49]-1*eorbit[e][53]-1*eorbit[e][59]);
			eorbit[e][24]=(f_24-1*eorbit[e][39]-1*eorbit[e][45]-1*eorbit[e][52]);
			eorbit[e][23]=(f_23-2*eorbit[e][36]-1*eorbit[e][45]-1*eorbit[e][52]-2*eorbit[e][58]-2*eorbit[e][60]);
			eorbit[e][22]=(f_22-1*eorbit[e][37]-1*eorbit[e][44]-1*eorbit[e][53]-1*eorbit[e][56]-1*eorbit[e][59]);
			eorbit[e][21]=(f_21-2*eorbit[e][38]-2*eorbit[e][43]-1*eorbit[e][52])/2;
			eorbit[e][20]=(f_20-1*eorbit[e][40]-1*eorbit[e][54]);
			eorbit[e][19]=(f_19-1*eorbit[e][33]-2*eorbit[e][41]-1*eorbit[e][45]-2*eorbit[e][50]-1*eorbit[e][52]-4*eorbit[e][58]-4*eorbit[e][60]);
			eorbit[e][18]=(f_18-2*eorbit[e][32]-2*eorbit[e][38]-1*eorbit[e][44]-6*eorbit[e][49]-1*eorbit[e][53]-2*eorbit[e][56]-2*eorbit[e][59]);
			eorbit[e][17]=(f_17-2*eorbit[e][25]-1*eorbit[e][27]-1*eorbit[e][32]-1*eorbit[e][34]-1*eorbit[e][47])/3;
			eorbit[e][16]=(f_16-2*eorbit[e][20]-2*eorbit[e][22]-1*eorbit[e][31]-2*eorbit[e][40]-1*eorbit[e][44]-2*eorbit[e][54])/2;
			eorbit[e][15]=(f_15-2*eorbit[e][25]-2*eorbit[e][29]-1*eorbit[e][31]-2*eorbit[e][32]-1*eorbit[e][34]-2*eorbit[e][42]-2*eorbit[e][47]);
			eorbit[e][14]=(f_14-1*eorbit[e][18]-2*eorbit[e][21]-1*eorbit[e][30]-2*eorbit[e][38]-1*eorbit[e][39]-2*eorbit[e][43]-1*eorbit[e][52])/2;
			eorbit[e][13]=(f_13-2*eorbit[e][22]-2*eorbit[e][28]-1*eorbit[e][31]-1*eorbit[e][40]-2*eorbit[e][44]-2*eorbit[e][54]);
			eorbit[e][12]=(f_12-2*eorbit[e][21]-2*eorbit[e][28]-2*eorbit[e][29]-2*eorbit[e][38]-2*eorbit[e][43]-1*eorbit[e][45]-1*eorbit[e][52]);
		}
	}

	endTime = clock();
	printf("%.2f\n", (double)(endTime-startTime)/CLOCKS_PER_SEC);

	endTime_all = endTime;
	printf("total: %.2f\n", (double)(endTime_all-startTime_all)/CLOCKS_PER_SEC);
}


fstream fin, fout; // input and output files
int GS=5;
string orbit_type;

int init(int argc, char *argv[]) {
	// open input, output files
	if (argc!=5) {
		cerr << "Incorrect number of arguments." << endl;
		cerr << "Usage: orca.exe [orbit type: node|edge] [graphlet size: 4/5] [graph - input file] [graphlets - output file]" << endl;
		return 0;
	}
	orbit_type = argv[1];
	if (orbit_type!="node" && orbit_type!="edge") {
		cerr << "Incorrect orbit type '" << argv[1] << "'. Should be 'node' or 'edge'." << endl;
		return 0;
	}
	sscanf(argv[2],"%d",&GS);
	if (GS!=4 && GS!=5) {
		cerr << "Incorrect graphlet size " << argv[2] << ". Should be 4 or 5." << endl;
		return 0;
	}
	fin.open(argv[3], fstream::in);
	fout.open(argv[4], fstream::out | fstream::binary);
	if (fin.fail()) {
		cerr << "Failed to open file " << argv[2] << endl;
		return 0;
	}
	if (fout.fail()) {
		cerr << "Failed to open file " << argv[3] << endl;
		return 0;
	}
	// read input graph
	fin >> n >> m;
	int d_max=0;
	edges = (PAIR*)malloc(m*sizeof(PAIR));
	deg = (int*)calloc(n,sizeof(int));
	for (int i=0;i<m;i++) {
		int a,b;
		fin >> a >> b;
		if (!(0<=a && a<n) || !(0<=b && b<n)) {
			cerr << "Node ids should be between 0 and n-1." << endl;
			return 0;
		}
		if (a==b) {
			cerr << "Self loops (edge from x to x) are not allowed." << endl;
			return 0;
		}
		deg[a]++; deg[b]++;
		edges[i]=PAIR(a,b);
	}
	for (int i=0;i<n;i++) d_max=max(d_max,deg[i]);
	printf("nodes: %d\n",n);
	printf("edges: %d\n",m);
	printf("max degree: %d\n",d_max);
	fin.close();
	if ((int)(set<PAIR>(edges,edges+m).size())!=m) {
		cerr << "Input file contains duplicate undirected edges." << endl;
		return 0;
	}
	// set up adjacency matrix if it's smaller than 100MB
	if ((int64)n*n < 100LL*1024*1024*8) {
		adjacent = adjacent_matrix;
		adj_matrix = (int*)calloc((n*n)/adj_chunk+1,sizeof(int));
		for (int i=0;i<m;i++) {
			int a=edges[i].a, b=edges[i].b;
			adj_matrix[(a*n+b)/adj_chunk]|=(1<<((a*n+b)%adj_chunk));
			adj_matrix[(b*n+a)/adj_chunk]|=(1<<((b*n+a)%adj_chunk));
		}
	} else {
		adjacent = adjacent_list;
	}
	// set up adjacency, incidence lists
	adj = (int**)malloc(n*sizeof(int*));
	for (int i=0;i<n;i++) adj[i] = (int*)malloc(deg[i]*sizeof(int));
	inc = (PII**)malloc(n*sizeof(PII*));
	for (int i=0;i<n;i++) inc[i] = (PII*)malloc(deg[i]*sizeof(PII));
	int *d = (int*)calloc(n,sizeof(int));
	for (int i=0;i<m;i++) {
		int a=edges[i].a, b=edges[i].b;
		adj[a][d[a]]=b; adj[b][d[b]]=a;
		inc[a][d[a]]=PII(b,i); inc[b][d[b]]=PII(a,i);
		d[a]++; d[b]++;
	}
	for (int i=0;i<n;i++) {
		sort(adj[i],adj[i]+deg[i]);
		sort(inc[i],inc[i]+deg[i]);
	}
	// initialize orbit counts
	orbit = (int64**)malloc(n*sizeof(int64*));
	for (int i=0;i<n;i++) orbit[i] = (int64*)calloc(73,sizeof(int64));
	// initialize edge orbit counts
	eorbit = (int64**)malloc(m*sizeof(int64*));
	for (int i=0;i<m;i++) eorbit[i] = (int64*)calloc(68,sizeof(int64));
	return 1;
}

void writeResults(int g=5) {
	int no[] = {0,0,1,4,15,73};
	for (int i=0;i<n;i++) {
		for (int j=0;j<no[g];j++) {
			if (j!=0) fout << " ";
			fout << orbit[i][j];
		}
		fout << endl;
	}
	fout.close();
}

void writeEdgeResults(int g=5) {
	int no[] = {0,0,0,2,12,68};
	for (int i=0;i<m;i++) {
		for (int j=0;j<no[g];j++) {
			if (j!=0) fout << " ";
			fout << eorbit[i][j];
		}
		fout << endl;
	}
	fout.close();
}

int main(int argc, char *argv[]) {
	if (!init(argc, argv)) {
		cerr << "Stopping!" << endl;
		return 0;
	}
	if (orbit_type=="node") {
		printf("Counting NODE orbits of graphlets on %d nodes.\n\n",GS);
		if (GS==4) count4();
		if (GS==5) count5();
		writeResults(GS);
	} else {
		printf("Counting EDGE orbits of graphlets on %d nodes.\n\n",GS);
		if (GS==4) ecount4();
		if (GS==5) ecount5();
		writeEdgeResults(GS);
	}
	return 0;
}