#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>

FILE *FichDist;
FILE *Out;

float	**Dist, DistMax, *D, *Dens, eps=.001, DgMoy;
char	**Et, Nom[20], Fich[50]="Data\\",**T;

int		N, NbClas, MaxCar=0, gpl, DgMin, DgMax;
unsigned int N2;
int		*Part, *Klas, *Kard;

static int compar(const void *e1, const void *e2) 
{
    float u = *(float *)e1 - *(float *)e2;
    return u == 0 ? 0 : (u < 0 ? +1 : -1);   /* tri d‚croissant */
}

float LecDist()
{	int i,j;
	float Dmax;

	printf("Name of the distance file ");
	gets(Nom); strcat(Fich,Nom); strcat(Fich,".Dis");
//	printf("%s\n",Fich);
	FichDist = fopen(Fich,"r"); assert(FichDist != NULL);
	strcpy(Fich,"Data\\ClasDens."); strcat(Fich,Nom); strcat(Fich,".out");
	Out = fopen(Fich,"w"); assert(Out != NULL);
	fprintf(Out,"Distance File : %s\n\n",Nom);

	fscanf(FichDist,"%d",&N);
	Dist = malloc((N) * sizeof(float *));
	Et = malloc((N) * sizeof(char *));
	T = malloc((N) * sizeof(char *));
	assert(Dist != NULL && Et != NULL && T != NULL);

	Dmax=0.;
	for (i=0;i<N;i++)
	{ 	Dist[i] = malloc((N) * sizeof(float));
		Et[i] = malloc(20);
		T[i] = malloc ((N) * sizeof(char) );	/* T[i][j]=1 ssi arete (i,j) */
		assert(Dist[i] != NULL && Et[i] != NULL && T[i] != NULL);

		fscanf(FichDist,"%s",Et[i]);
        if (strlen(Et[i])>MaxCar) MaxCar=strlen(Et[i]);
		for (j=0;j<N;j++)
		{	fscanf(FichDist,"%f",&Dist[i][j]);
			if (i!=j && Dist[i][j]<eps) Dist[i][j]=eps;
			if (Dist[i][j]>Dmax) Dmax=Dist[i][j];
	}	} 	/* printf("Dmax = %.2f\n",Dmax); */
  	fclose(FichDist);
    MaxCar += 1; gpl=80/MaxCar;
  	return Dmax;
}

void Vue()
{   int i,j;

    for (i=0; i<N; i++)
    {   printf("%*s : ",-MaxCar,Et[i]);
        for (j=0; j<N; j++) printf("%d",T[i][j]); printf(" : ");
        printf(" %.3f\n",Dens[i]);
    }   printf("\n");
//    fflush(stdin); getchar();
    return;
}

void ClasOut(int typ)
{   int i, j, k;

    for (k=1; k<=NbClas; k++)
	{   printf("Class %d : ",k);
	    fprintf(Out,"Class %d\n",k);
	    j=0;
		for (i=0; i<N; i++)
			if (Part[i]==k)
			{	printf("%s ",Et[i]);
			    fprintf(Out,"%*s ",-MaxCar,Et[i]);
			    j++; if (j%gpl==0) fprintf(Out,"\n");
            }
		printf("\n"); fprintf(Out,"\n\n");
	}   printf("\n"); fprintf(Out,"\n\n");
	if (typ==1)
    {	for (i=0; i<N; i++)
    	{	printf("%*s : ",MaxCar,Et[i]);
    	    for (j=0; j<N; j++)
                if (T[i][j]==1) printf("%*s ",-MaxCar,Et[j]);
            printf(" : %.2f\n",Dens[i]);
        }
    	printf("\n");
   	}
	return;
}

void QualiPart()
/*	PdI = poids des distances intraclasses
	PdE = poids des distances interclasses
	SdG = Somme des Cocy plus grandes distances
	SdP = Somme des N2-Cocy plus petites distances
*/
{	unsigned int	i,j,k, clas, DifSym=0, Cocy=0;
	unsigned long	NbT=0, NbGT=0;
	float	PdI=0., PdE=0., SdG=0., SdP=0., taux, Ssup,Sinf,seuil, Dij;

		for (i=1; i<NbClas; i++)
			for (j=i+1; j<=NbClas; j++)
				Cocy += Kard[i]*Kard[j];
		seuil=D[Cocy-1];
		if (seuil-D[Cocy]>.01) { Sinf=D[Cocy]; Ssup=seuil; }
		else {	i=Cocy-2; 
				while (D[i]-seuil<.01 && i>0) i--;
				Ssup=D[i];
				j=Cocy+1;
				while (seuil-D[j]<.01 && j<N2-1) j++;
				Sinf=D[j];
				// for (k=i; k<=j; k++) printf("%.1f ",D[k]);
			 }	// printf("seuils : %.1f ü %.1f > %.1f\n",Ssup,seuil,Sinf);
		for (i=1; i<N; i++)
			for (j=0; j<i; j++)
			{	if (Part[i]==Part[j] && Dist[i][j]>=Ssup) DifSym += 1;
				if (Part[i]!=Part[j] && Dist[i][j]<=Sinf) DifSym += 1;
			}
		taux=1.0*(N2-DifSym)/N2;
		printf("Rate of agreements on pairs %.3f\n",taux);
		fprintf(Out,"Rate of agreements on pairs %.3f\n",taux);

		for (i=1; i<N; i++)
			for (j=0; j<i; j++)
				if (Part[i]==Part[j]) PdI += Dist[i][j]; else PdE += Dist[i][j];
		for (i=0; i<Cocy; i++) SdG += D[i]; // printf("SdG = %.1f ",SdG);
		for (i=Cocy; i<N2; i++) SdP += D[i]; // printf("SdP = %.1f ",SdP);
		taux=(SdP/PdI)*(PdE/SdG);
		printf("Rate of weights %.3f\n",taux);
		fprintf(Out,"Rate of weights %.3f\n",taux);

		taux=(PdE/Cocy)/(PdI/(N2-Cocy));
		printf("Ratio of average lengths %.2f\n",taux);
		fprintf(Out,"Ratio of average lengths %.2f\n",taux);
		
		if (N>1000) return;
		for (i=1; i<N; i++)
			for (j=0; j<i; j++)
			{	if (Part[i]!=Part[j]) continue;
				Dij=Dist[i][j];
				for (k=0; k<N; k++)
				{	if (Part[k]==Part[i]) continue;
					NbT++;
					if (Dij<=Dist[i][k] && Dij<=Dist[j][k]) NbGT++;
					// else printf("(%d %d %d) ",i,j,k); 
			}	}
		taux=1.0*NbGT/NbT;
		printf("Rate of well designed triples %.3f\n",taux);
		fprintf(Out,"Rate of well designed triples %.3f\n\n",taux);
 		return;
}

/*--------------------------------------------------------------*/

int Densite(int typ, float seuil)
//  Construit la table T d'incidence et calcule les degr‚s et la densit‚ selon la formule typ
{   int i, j, k, d, Na=0;
    float Dmin, Dmax, SD;
    
    DgMin=N; DgMax=-1;
    for (i=0; i<N; i++)
	{   for (j=0; j<N; j++) T[i][j]=0;
		d=0; SD=0.; 
        Dmin=seuil+eps; Dmax=0.;
		for (j=0; j<N; j++)
		{	if (Dist[i][j]<seuil && j!=i) 
            {   T[i][j]=1;
                d++; SD += Dist[i][j]; 
                if (Dist[i][j]<Dmin) Dmin=Dist[i][j]; 
                if (Dist[i][j]>Dmax) Dmax=Dist[i][j];
            }   
        }
        if  (d<DgMin) DgMin=d; else if (d>DgMax) DgMax=d;
        Na += d;
        switch(typ)
        {   case 1:
            if (d==0) Dens[i]=0; else Dens[i]=(seuil-Dmin)/seuil;
            break;
            case 2: 
            if (d==0) Dens[i]=0; else Dens[i]=(seuil-(SD/d))/seuil;
            break;
            case 3: 
            if (d==0) Dens[i]=0; else Dens[i]=(seuil-Dmax)/seuil;
        }
    }   Na=Na/2;
    return Na;
}

int Noyaux(float DeMoy)
{   int     i,j,k, flag, card, piv, cl, NbCl;

	for (i=0; i<N; i++) { Klas[i]=0; Part[i]=0; }
// Recherche des extrema locaux
	DeMoy += eps;
	for (i=0; i<N; i++)
	{	if (Part[i]!=0 || Dens[i]<DeMoy) continue;   
	    Part[i]=1; // a priori c'est un optimum local
	    for (j=0; j<N; j++)
	    {   if (T[i][j]==0) continue;
	        if (Dens[j]-Dens[i]>eps) { Part[i]=-1; break; } // ya un adjascent de densit‚ sup‚rieure
	        if (Dens[i]-Dens[j]>eps) Part[j]=-1;  // C'est j qui est domin‚
        }   
    }
//   	printf("Opt : "); for (i=0; i<N; i++) printf("%d ",Part[i]); fflush(stdin); getchar();

//  On num‚rote les composantes connexes des optima 
	NbClas=1;
	for (i=0; i<N; i++)
    {   if (Part[i]!=1) continue;
        NbClas++; Part[i]=NbClas;
   		card=1; Klas[0]=i;
	    k=0; j=0;
	    while (k<card)
		{	piv=Klas[k]; 
			for (j=0; j<N; j++)
			{	if (T[piv][j]==0 || Part[j]!=1) continue; 
				Klas[card]=j; card++; Part[j]=NbClas;
			}	k++; 
		}    
    }
    for (i=0; i<N; i++) if (Part[i]>1) Part[i]--; else Part[i]=0;
    NbClas--;
//   	printf("Cop : "); for (i=0; i<N; i++) printf("%d ",Part[i]); fflush(stdin); getchar();
    
// On met das les noyaux les ‚l‚ments qui sont adjacents … un seul maximum local
// et qui ont une densit‚ sup‚rieure ou ‚gale … Demoy
	for (i=0; i<N; i++)
	{	if (Part[i]>0 || Dens[i]<DeMoy) continue;
        flag=0; cl=0;
	    for (j=0; j<N; j++)
	        if (T[i][j]==1 && Part[j]>0) 
                if (flag==0) { cl=Part[j]; flag=1; }
                else if (Part[j]!=cl) { flag=2; break; }
        if (flag==1) Part[i]=cl;
    }     
//   	printf("Noy : "); for (i=0; i<N; i++) printf("%d ",Part[i]); printf("\n");
    NbCl=0;
    for (i=0; i<N; i++) if (Part[i]) NbCl++;
    printf("There will be %d classes and %d elements in the kernels\n\n",NbClas,NbCl);
//    fflush(stdin); getchar();
    return NbCl;
}

int Extension()
// Construit les extension des classes
{   int     i, j, k,kk, k1, k2, NbCl, NbNc=0, flag=1, max;
    float   Dmax;

//    On affecte ceux qui ne sont connect‚s qu'… une seule classe
    for (i=0; i<N; i++)
    {   if (Part[i]>0) continue;
        for (k=0; k<=NbClas; k++) Kard[k]=0;
        for (j=0; j<N; j++)
            if (T[i][j]>0) Kard[Part[j]]++;
//        printf("%*s : ",-MaxCar,Et[i]);
//        for (k=1; k<=NbClas; k++) printf("%2d ",Kard[k]);
        NbCl=0;        // Nb. de classement possibles
        for (k=1; k<=NbClas; k++) 
            if (Kard[k]>0) { NbCl++; kk=k; }
        if (NbCl==1) Part[i]=kk; else NbNc++; // printf("-> %d",kk); 
//        printf("\n");
    }   // printf("\n");
//   	printf("Ext : "); for (i=0; i<N; i++) printf("%d ",Part[i]); printf("\n");
    printf("Extended classes (Nb.of elements %d)\n\n",N-NbNc);
    fprintf(Out,"Extended classes (Nb.of elements %d)\n\n",N-NbNc);
    ClasOut(0); fflush(stdin); getchar();
    
//    On affecte … la classe majoritaire ou a la classe … distance moyenne minimum (valeur dans Dens)
//  printf("DistMax = %.2f\n",DistMax);
    while (flag>0 && NbNc>0)
    {   flag=0;
        for (i=0; i<N; i++)
        {   if (Part[i]>0) continue;
            for (k=0; k<=NbClas; k++) { Kard[k]=0; Dens[k]=0.; }
            for (j=0; j<N; j++)
                if (T[i][j]>0 && Part[j]>0) { Kard[Part[j]]++; Dens[Part[j]] += Dist[i][j]; }
            for (k=1; k<=NbClas; k++) 
                if (Kard[k]>0) Dens[k]=Dens[k]/Kard[k];
/*            printf("%s : ",Et[i]);
            for (k=1; k<=NbClas; k++) printf("%3d ",Kard[k]); printf(" ; ");
            for (k=1; k<=NbClas; k++) printf("%2.2f ",Dens[k]); printf("\n");
*/
            max=0; k1=0;
            for (k=1; k<=NbClas; k++)
                if (Kard[k]>max) { max=Kard[k]; k1=k; }
            if (max==0) continue;
            Dmax=DistMax+eps; k2=0;       // Classe majoritaire
            for (k=1; k<=NbClas; k++)
                if (Kard[k]>0 && Dens[k]<Dmax) { Dmax=Dens[k]; k2=k; }
            if (k1==k2 || Kard[k1]*2>Kard[k2]*3) { NbNc--; flag=1; Part[i]=k1; }
            else { NbNc--; flag=1; Part[i]=k2; }
        }
    }
//   	printf("Tot : "); for (i=0; i<N; i++) printf("%d ",Part[i]); printf("\n");
    printf("Final classes\n\n");
    fprintf(Out,"Final classes\n\n");
    ClasOut(0); 
    return NbNc;
}


main()
{	int 	i,j,k, typ=0, prc, Na, flag=1, NbNc, NbCl;
	float 	seuil, Dmin, DMaxMin, DeMin,DeMoy,DeMax;
	unsigned int coup;

	DistMax=LecDist();
	N2=N*(N-1)/2;
	Part = malloc((N) * sizeof(int));			/* # de classes pour chaque ‚l‚ment */
	Klas = malloc((N) * sizeof(int));			/* Liste des ‚l‚ments d'une classe */
	Kard = malloc((N) * sizeof(int));			/* Liste des ‚l‚ments d'une classe */
	Dens = malloc((N) * sizeof(float));         /* Densit‚ en chaque sommet du graphe seuil */
	D = malloc((N2) * sizeof(float));           /* Distances dans un tableau unidimentionnel */

	assert(Part != NULL && Klas != NULL && Kard != NULL && Dens != NULL && D != NULL);

    k=0; DMaxMin=0.;
	for (i=0; i<N-1; i++)
	{   Dmin=DistMax;
        for (j=0; j<i; j++) 
            if (Dist[i][j]<Dmin) Dmin=Dist[i][j];
    	for (j=i+1; j<N; j++)
        {	D[k++] = Dist[i][j];
            if (Dist[i][j]<Dmin) Dmin=Dist[i][j];
        }
        if (Dmin>DMaxMin) DMaxMin=Dmin;
    }
    printf("Maximum distance value %.2f\n",DistMax);
    printf("MaxMin distance value %.2f\n\n",DMaxMin);
    
	qsort(D, N2, sizeof(float), compar);
	
	while (flag)
    {	typ=0;
        while (typ!=1 && typ != 2)
        {	printf("Do you want to define the graph by\n");
            printf("    a percentage of edges (1) or\n");
            printf("    a threshold distance value (2)\n");
            printf("Your choice ");
        	scanf("%d",&typ); printf("\n");
       	}
        if (typ==1)
        {   printf("Pourcentage of edges in the threshold graph ");
            scanf("%d",&prc);
            coup=N2-(N2*prc/100); // printf("coup %d\n",coup);
        	seuil=D[coup-1]+eps; // printf("Seuil %.2f\n",seuil);
            printf("Distance threshold %.2f\n",seuil-eps);
    	}
    	else if (typ==2)
    	{   printf("Threshold distance value ");
    	    scanf("%f",&seuil); seuil += eps;
        }
        Na=Densite(2,seuil);
        printf("Nb. of edges %d (percentage %.3f)\n",Na,1.*Na/N2);
        DgMoy=2.*Na/N;
        printf("Degree : Minimum %d, Average %.1f, Maximum %d\n",DgMin,DgMoy,DgMax);
        fprintf(Out,"Pourcentage of edges in the threshold graph %d\n",prc);
        fprintf(Out,"Distance threshold %.2f\n",seuil-eps);
        fprintf(Out,"Nb. of edges %d (percentage %.3f)\n",Na,1.*Na/N2);
        fprintf(Out,"Degree : Minimum %d, Average %.1f, Maximum %d\n",DgMin,DgMoy,DgMax);
        if (N<=50) Vue();
    
        j=0;
        DeMin=Dens[0]+eps; DeMoy=0.; DeMax=0.;
        for (i=0; i<N; i++)
        {   if (Dens[i]>DeMax) DeMax=Dens[i];
            DeMoy+=Dens[i];
            if (Dens[i]<DeMin) DeMin=Dens[i];
        }   
       	DeMoy=DeMoy/N;
        NbCl=Noyaux(DeMoy);
        printf("Do you want to cluster this graph (0) or to test another one (1) ");
        fflush(stdin); scanf("%d",&flag); printf("\n");
    }
    fprintf(Out,"Nb. of classes %d\n\n",NbClas);
    printf("Kernels  (Nb. of elements %d)\n\n",NbCl); 
    fprintf(Out,"Kernels (Nb. of elements %d)\n\n",NbCl);
    ClasOut(0); fflush(stdin); getchar();

	NbNc=Extension();
	
    if (NbNc>0)
    {	printf("Nb. of unclustered elements : %d\n",NbNc);
    	fprintf(Out,"Nb. of unclustered elements : %d\n",NbNc);
    	k=0;
    	for (i=0; i<N; i++) 
           if (Part[i]==0) 
           {   fprintf(Out,"%*s ",MaxCar,Et[i]);
               k++; if (k%gpl==0)fprintf(Out,"\n");
           }
        printf("\n\n"); fprintf(Out,"\n\n");
    }
    else 
    {	fprintf(Out,"Quality criteria of the final partition\n\n");
    	printf("Quality criteria of the final partition\n");
    	for (k=0; k<=NbClas; k++) Kard[k]=0;
    	for (i=0; i<N; i++) Kard[Part[i]]++; 
    	QualiPart();
   	}
    
	printf("C'est fini ");
    fflush(stdin); getchar();

}