/* ex_dwgraphs.c  */

#include "def.h"
#include "macro.h"
#include <math.h>

static INT reverse_vector();
static INT number_to_2set();
static INT ausgabe_graph();

#define VERBOSE 0

int noprint=0; /* 1 = noprint 0 = print */
double globalenergy;

INT anzahlknoten;
INT anzahlkanten;
INT maxkanten;  
INT complement=0;  

main()
{
	int eingabe;
	anfang();
		again:
		printf("menu \n");
		printf(" (1) input vertices/edges/candidates \n");
		printf(" (3) input vertices/edges/candidates/database \n");
		printf(" (2) input vertices/all edges/250 candidates \n");
		printf(" (4) input vertices/edges between /candidates/database \n");
		printf(" (99)quit \n");
		scanf("%d",&eingabe);
		switch(eingabe)
			{
			case 1: onesetofdata(); goto again;
			case 3: onesetofdatadatabase(); goto again;
			case 4: moresetofdatadatabase(); goto again;
			case 2: data2(); goto again;
			case 99:break;
			}

	ende();
}


data2()
{
	OP a,b,c,d,e,f,g,h,aa;
	OP max,zz,anz,len,rep; 
	OP maxknoten;
	double sumofenergy;
	int i,j,k;
	FILE *fp;
	a=callocobject(); 
	b=callocobject();
	c=callocobject();
	d=callocobject(); 
	e=callocobject();
	f=callocobject();
	g=callocobject();
	h=callocobject();
	aa=callocobject();
	max=callocobject();
	anz=callocobject();
	len=callocobject();
	rep=callocobject();
	zz=callocobject();
	maxknoten=callocobject();
	printeingabe("number of vertices");
	scan(INTEGER,a);
	anzahlknoten=S_I_I(a);
	binom(a,cons_zwei,maxknoten);
	printeingabe("number of edges");
	//scan(INTEGER,aa);
	//printeingabe("number of random graphs to be generated");
for (M_I_I(1,aa); le(aa,maxknoten); inc(aa) )
{
	// scan(INTEGER,anz); 
	sumofenergy=0.0;
	m_i_i(250,anz);
	first_part_VECTOR(a,b);
	m_il_v(1L,h); 
	i=0L;
	do {
#if VERBOSE > 0
		printf("partition:                                  ");
		println(b);
#endif
		m_part_perm(b,c);
#if VERBOSE > 0
		printf("permutation in S_%d:                         ",S_I_I(a));
		println(c);
#endif
		m_perm_2setsperm(c,d);
#if VERBOSE > 0
		printf("induced permutation                         ");
		println(d);
#endif
		ordcon(b,e);
#if VERBOSE > 0
		printf("order of the conjugacy class                ");
		println(e);
#endif
		berechne_fixpunkt_anzahl(d,aa,f);
		mult(e,f,g);
#if VERBOSE > 0
		printf("probability of this conjugacy class         ");
		println(g);
#endif
		m_il_v(3L,S_V_I(h,i));
		copy(f,S_V_I(S_V_I(h,i),0L));
		copy(d,S_V_I(S_V_I(h,i),1L));
		copy(g,S_V_I(S_V_I(h,i),2L));
		inc(h);
		++i;
	} while (next_apply(b));
	dec(h);
	reverse_vector(h);
	copy(S_P_L(S_V_I(S_V_I(h,0L),1L)),len);
	k=i;
	for (i=1L;i<k;++i)
		add_apply(S_V_I(S_V_I(h,i-1L),2L),S_V_I(S_V_I(h,i),2L));
	copy(S_V_I(S_V_I(h,k-1L),2L),max);
	if (S_O_K(max)==INTEGER) t_int_longint(max,max);
#if VERBOSE > 0
	println(max);
#endif
	for (j=0L;j<S_I_I(anz);++j) {
		random_longint(zz,max);
#if VERBOSE > 0
		println(zz);
#endif
		for (i=0L;i<k;++i)
			if (le(zz,S_V_I(S_V_I(h,i),2L))) {
#if VERBOSE > 0
				printf("Class %d\n",i);
#endif
				break;
			}
		if (i==0L) {
#if VERBOSE > 0
			printf("id \n");
#endif
			fixed_point_id(len,aa,rep);
		}
		else {
#if VERBOSE > 0
			printf("non_id \n");
#endif
			fixed_point_nonid(S_V_I(S_V_I(h,i),1L),S_V_I(S_V_I(h,i),0L),aa,rep);
		}
		// println(rep);
		noprint=1;
		ausgabe_graph(rep);
		sumofenergy += globalenergy;
		noprint=0;
	}
	printf("number of edges ");print(aa);
	printf(" sumofenergy = %f average =%f \n",sumofenergy,sumofenergy/S_I_I(anz));
} /* loop over all edges */
	freeall(a);
	freeall(b);
	freeall(c);
	freeall(d);
	freeall(e);
	freeall(f);
	freeall(g);
	freeall(h);
	freeall(aa);
	freeall(max);
	freeall(anz);
	freeall(maxknoten);
	freeall(len);
	freeall(rep);
	freeall(zz);
	print_time(); 
}
onesetofdatadatabase()
	{
	onesetofdata_co(2);
	}
moresetofdatadatabase()
	{
	onesetofdata_co(3);
	}
onesetofdata()
	{
	onesetofdata_co(1);
	}
onesetofdata_co(para) INT para;
{
	OP a,b,c,d,e,f,g,h,aa,bb;
	OP max,zz,anz,len,rep; 
	int i,j,k;
	FILE *fp;
	a=callocobject(); 
	b=callocobject();
	c=callocobject();
	d=callocobject(); 
	e=callocobject();
	f=callocobject();
	g=callocobject();
	h=callocobject();
	aa=callocobject();
	bb=callocobject();
	max=callocobject();
	anz=callocobject();
	len=callocobject();
	rep=callocobject();
	zz=callocobject();
	printeingabe("number of vertices");
	scan(INTEGER,a);
	anzahlknoten=S_I_I(a);
	maxkanten = (((anzahlknoten % 2) == 0 )? (anzahlknoten/2)*(anzahlknoten-1) : (((anzahlknoten-1)/2)*anzahlknoten)) ;
	if (para == 3) {
		printeingabe("lower number of edges");
		scan(INTEGER,aa);
		printeingabe("upper number of edges");
		scan(INTEGER,bb);
		}
	else	{
		printeingabe("number of edges");
		scan(INTEGER,aa);
		}
	printeingabe("number of random graphs to be generated");
	scan(INTEGER,anz);
nochmal:
	anzahlkanten=S_I_I(aa);
	if (anzahlkanten > (maxkanten/2)) // besser complement
		{ anzahlkanten = maxkanten-anzahlkanten; complement = 1; 
		  m_i_i(anzahlkanten,aa);
		}
	else
		{ complement = 0; }
	first_part_VECTOR(a,b);
	m_il_v(1L,h); 
	i=0L;
	do {
#if VERBOSE > 0
		printf("partition:                                  ");
		println(b);
#endif
		m_part_perm(b,c);
#if VERBOSE > 0
		printf("permutation in S_%d:                         ",S_I_I(a));
		println(c);
#endif
		m_perm_2setsperm(c,d);
#if VERBOSE > 0
		printf("induced permutation                         ");
		println(d);
#endif
		ordcon(b,e);
#if VERBOSE > 0
		printf("order of the conjugacy class                ");
		println(e);
#endif
		berechne_fixpunkt_anzahl(d,aa,f);
		mult(e,f,g);
#if VERBOSE > 0
		printf("probability of this conjugacy class         ");
		println(g);
#endif
		m_il_v(3L,S_V_I(h,i));
		copy(f,S_V_I(S_V_I(h,i),0L));
		copy(d,S_V_I(S_V_I(h,i),1L));
		copy(g,S_V_I(S_V_I(h,i),2L));
		inc(h);
		++i;
	} while (next_apply(b));
	dec(h);
	reverse_vector(h);
	copy(S_P_L(S_V_I(S_V_I(h,0L),1L)),len);
	k=i;
	for (i=1L;i<k;++i)
		add_apply(S_V_I(S_V_I(h,i-1L),2L),S_V_I(S_V_I(h,i),2L));
	copy(S_V_I(S_V_I(h,k-1L),2L),max);
	if (S_O_K(max)==INTEGER) t_int_longint(max,max);
#if VERBOSE > 0
	println(max);
#endif
	for (j=0L;j<S_I_I(anz);++j) {
		random_longint(zz,max);
#if VERBOSE > 0
		println(zz);
#endif
		for (i=0L;i<k;++i)
			if (le(zz,S_V_I(S_V_I(h,i),2L))) {
#if VERBOSE > 0
				printf("Class %d\n",i);
#endif
				break;
			}
		if (i==0L) {
#if VERBOSE > 0
			printf("id \n");
#endif
			fixed_point_id(len,aa,rep);
		}
		else {
#if VERBOSE > 0
			printf("non_id \n");
#endif
			fixed_point_nonid(S_V_I(S_V_I(h,i),1L),S_V_I(S_V_I(h,i),0L),aa,rep);
		}
		println(rep);
		ausgabe_graph(rep);
		/*if ((para==2) || (para==3)) // into database
			{
			if (complement==1) 
				save_datensatz(anzahlknoten,maxkanten-anzahlkanten,globalenergy);
			else
				save_datensatz(anzahlknoten,anzahlkanten,globalenergy);
			}*/
	}

	if (para==3) {
		if (complement==1) m_i_i(maxkanten-anzahlkanten,aa);
		inc(aa); if (le(aa,bb)) goto nochmal;
		}
	freeall(a);
	freeall(b);
	freeall(c);
	freeall(d);
	freeall(e);
	freeall(f);
	freeall(g);
	freeall(h);
	freeall(aa);
	freeall(bb);
	freeall(max);
	freeall(anz);
	freeall(len);
	freeall(rep);
	freeall(zz);
	print_time(); 
}

static INT reverse_vector(a) OP a;
{
	OP hilf=callocobject();
	INT i,erg=OK;
	for (i=0L;2L*i<S_V_LI(a);++i) {
		erg+=copy(S_V_I(a,i),hilf);
		erg+=copy(S_V_I(a,S_V_LI(a)-1L-i),S_V_I(a,i));
		erg+=copy(hilf,S_V_I(a,S_V_LI(a)-1L-i));
	}
	erg+=freeall(hilf);
	ENDR("reverse_vector");
}

static INT number_to_2set(a,i,j) INT a,*i,*j;
{
	INT k,l=0L;
	for (k=1L;;++k) {
		l+=k;
		if (a<=l) break;
	}
	*j=k+1L;
	*i=a-l+k;
}

double energy(OP m) 
/* computes the energy of the adjacency matrix m */
	{
	double a[42][42];
	double h[42][42];
	double energy;
	int i,j;long n = anzahlknoten;
        if (n > 42) { printf("number of vertices > 42; sorry; exiting!\n"); exit(0); }
	for (i=0;i<42;i++)
	for (j=0;j<42;j++) a[i][j] = 0.0;

	for (i=0;i<anzahlknoten;i++)
	for (j=0;j<anzahlknoten;j++)
		a[i][j]= (double) S_M_IJI(m,i,j);

	diag_(&n,a,h);
	if (noprint==0)printf("eigenvalues \n");
	energy = 0.0;
	for (j=0;j<anzahlknoten;j++)
		{
		energy += fabs(a[j][j]);
		if (noprint==0) printf("%3.3f ",a[j][j]);
		}
	if (noprint==0)printf("\n");
	return energy;
	}



static INT ausgabe_graph(a) OP a;
{
	INT i,j,k;
	OP mat;
	mat = callocobject();
	m_ilih_nm(anzahlknoten,anzahlknoten,mat);
	
	for (i=0L;i<S_V_LI(a);++i)
		if (S_V_II(a,i)==1L) {
			number_to_2set(i+1L,&j,&k);
			if (noprint==0)printf("{%d,%d},",j,k);
			M_I_I(1,S_M_IJ(mat,j-1,k-1));
			M_I_I(1,S_M_IJ(mat,k-1,j-1));
		}
	if (complement==1) // invert rep
		{
		int ii,jj;
		for (ii= 0;ii<S_M_HI(mat);ii++)
		for (jj= 0;jj<S_M_LI(mat);jj++)
			{
			if (ii!=jj){
				if (S_M_IJI(mat,ii,jj) == 0) m_i_i(1,S_M_IJ(mat,ii,jj));
				else m_i_i(0,S_M_IJ(mat,ii,jj));
					}
			}
		}
	if (noprint==0)printf("\n");
	if (noprint==0)println(mat);

	/* now the energy AK 211100 */

	{
	double e;
	e = energy(mat);
	if (noprint==0)printf("energy = %3.3f\n",e);
	globalenergy=e;
	zeilenposition=0;
	}
}