#include <stdio.h> 
/*#include <iostream.h>
*/
#include <stddef.h> 
/*#include <fstream.h>
*/
#include <stdlib.h>
#include <math.h> 

#define TFACTR 0.9
#define NR_END 1
#define FREE_ARG char*
#define NRANSI
#define IAran1 16807
#define IMran1 2147483647
#define AMran1 (1.0/IMran1)
#define IQran1 127773
#define IRran1 2836   
#define NDIVran1 (1+(IMran1-1)/32)
#define EPSran1 1.2e-7   
#define RNMXran1 (1.0-EPSran1)

/*int nunits=139;
int nclst=8;
int nMethod=9;
*/

void nrerror(char error_text[]);
float *vector(long nl, long nh);
int *ivector(long nl, long nh);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_matrix(float **m, long nrl,long nrh,long ncl,long nch);
float **matrix(long nrl,long nrh,long ncl,long nch);
float ALEN(float xdata[], float ydata[], int ngene, int ncity);
int irbit1(unsigned long *iseed);
int metrop(float de, float t);
float ran1(long *idum);
float revcst(float x[], float y[], int iorder[], int ngene,
             int ncity, int n[]);
void reverse(int iorder[], int ncity, int n[]);
void reversey(float y[], int iorder[], int ncity, int nunits, int n[]);
float trncst(float x[], float y[], int iorder[],
             int ngene, int ncity, int n[]);
void trnspt(int iorder[], int ncity, int n[]);
void trnspty(float y[], int iorder[], int ncity, int nunits, int n[]);
float anneal(float x[], float y[], int iorder[], int ngene, int ncity);

clustComp(ans, cldata, nunits, nmethod, ncl)
   float *ans, *cldata;
   int *nunits,*nmethod,*ncl;   
{
   int i,j,k,l;
   int *iorder;
   float *xdat,*ydat,**dataMat;
   iorder=ivector(1, *ncl);
   xdat=vector(1, *nunits);
   ydat=vector(1, *nunits);
   dataMat=matrix(1, *nunits,1, *nmethod);
   for(i=1; i<= *ncl; i++){
      iorder[i] = i;
   }
   l=1;
   for(i=1;i<= *nmethod;i++){
      for(j=1;j<= *nunits;j++){
         dataMat[j][i]=cldata[l-1];
         l+=1;
      }
   }
   l=0;
   for(i=1;i<= *nmethod;i++){
      for(k=1;k<= *nunits;k++) xdat[k]=dataMat[k][i];
      for(j=i+1;j<= *nmethod;j++){
         for(k=1;k<= *nunits;k++) ydat[k]=dataMat[k][j];
         ans[l]=anneal(xdat, ydat, iorder, *nunits, *ncl);
         l+=1; 
      }
   }
   free_ivector(iorder,1, *ncl);
   free_vector(xdat,1, *nunits);
   free_vector(ydat,1, *nunits);
   free_matrix(dataMat,1, *nunits,1, *nmethod);
}

 
void nrerror(char error_text[])
/*numberical receipes standard srroe handler*/
{
   fprintf(stderr, "Numerical Receipes run-time error..\n");
   fprintf(stderr, "%s\n", error_text);
   fprintf(stderr,"...now exiting to system...\n");
   exit(1);
}

float **matrix(long nrl,long nrh,long ncl,long nch)
/* allocate a matrix with subscript range m[nrl...nrh][ncl...nch] 
*/
{   
   long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
   float **m;
   void nrerror(char error_text[]);
   /* allocate pointers to rows */
   m=(float **) malloc((size_t)((nrow+1)*sizeof(float*)));
   if (!m) nrerror("allocation failure 1 in matrix()");
   m += 1;
   m -= nrl;
   /* allocate rows and set pointers to them */
   m[nrl]=(float *) malloc((size_t)((nrow*ncol+1)*sizeof(float)));
   if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
   m[nrl] += 1;
   m[nrl] -= ncl;
   for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
   /* return pointer to array of pointers to rows */
   return m;
}

void free_matrix(float **m, long nrl,long nrh,long ncl,long nch)
/* free a matrix allocated by dmatrix() */
{
   free((FREE_ARG) (m[nrl]+ncl-1));
   free((FREE_ARG) (m+nrl-1));
} 

float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl .. nh] */
{   
   float *v;
   v = (float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
   if (!v) nrerror("allocation failure in vector()");
   return v-nl+NR_END;
}

void free_vector(float *v, long nl, long nh)

{
   free((FREE_ARG) (v+nl-NR_END));
}   
   

/*calculate negative kappa value based on two cluster analyses */

float ALEN(float xdata[], float ydata[], int ngene, int ncity)

/* ncity is number of groups in cluster analysis; */
/*ngene is number of gene */
/* xdata[] is first cluster analysis */
/* ydata[] is the second cluster analysis */
{
   int j, i;
   float *sumi, *sumj;
   float sum1=0.0, sum2 =0.0, negkappa;
   for (i = 1; i <= ngene; i++){
      if (xdata[i] == ydata[i]) sum1 = sum1 +  1; 
   }
/* marginal number of gene for method 1 xdata */
   sumi = vector(1, ncity); 
/* marginal number of gene for method 2 ydata */
   sumj = vector(1, ncity); 
   for(i = 1; i <= ncity; i++){
      sumi[i] = 0.0;
      for(j = 1; j <= ngene; j++){
         if (xdata[j] == i) sumi[i] = sumi[i] +  1;
      }
   }
   for(i = 1; i <= ncity; i++){
      sumj[i] = 0.0;
      for(j = 1; j <= ngene; j++){
         if (ydata[j] == i) sumj[i] = sumj[i] +  1;
      }
   }   
   for(i = 1; i <= ncity; i++){
      sum2 += sumi[i]*sumj[i];
   }   
/*negative kappa value*/
   negkappa = (sum2 - ngene*sum1)/(ngene*ngene - sum2); 
   free_vector(sumj,1,ncity);
   free_vector(sumi,1,ncity);
   return negkappa;
}
   
float anneal(float x[], float y[], int iorder[], int ngene, int ncity)
{
	int ans,nover,nlimit,i1,i2;
	int i,j,k,nsucc,nn,idec;
	static int n[7];
	long idum;
	unsigned long iseed;
	float de,t;
        static float negkap;
	nover=100*ncity;
	nlimit=10*ncity;
	negkap=0.0;
	t=0.5;
        /* calculate initial negative kappa value */ 
        negkap = ALEN(x, y, ngene, ncity);
/*        printf("Initial negative kappa is %f\n", negkap);*/
        idum = -1;
	iseed=111;
	for (j=1;j<=100;j++) {
/*           printf("negative kappa is %f\n", negkap);
*/		nsucc=0;
		for (k=1;k<=nover;k++) {
/*           printf("get here %f\n", negkap);*/
			do {
				n[1]=1+(int) (ncity*ran1(&idum));
/*           printf("get here %i\n", nn);*/
				n[2]=1+(int) ((ncity-1)*ran1(&idum));
				if (n[2] >= n[1]) ++n[2];
				nn=1+((n[1]-n[2]+ncity-1) % ncity);
			} while (nn<3);
			idec=irbit1(&iseed);
/*           printf("get here %f\n", negkap);*/
			if (idec == 0) {
				n[3]=n[2]+(int) (abs(nn-2)*ran1(&idum))+1;
				n[3]=1+((n[3]-1) % ncity);
				de=trncst(x,y,iorder,ngene,ncity,n);
				ans=metrop(de,t);
				if (ans) {
                                        ++nsucc;
					negkap += de;
					trnspty(y,iorder,ncity,ngene,n);
				}
			} else {
				de=revcst(x,y,iorder,ngene,ncity,n);
				ans=metrop(de,t);
				if (ans) {
					++nsucc;
					negkap += de;
					reversey(y,iorder,ncity,ngene,n);
				}
			}
			if (nsucc >= nlimit) break;
		}
/*		printf("\n %s %10.6f %s %12.6f \n","T =",t,
			"	 Negative Kappa  =",negkap);
		printf("Successful Moves: %6d\n",nsucc);
*/		t *= TFACTR;
		if (nsucc == 0){
/*                   cout << negkap << endl;
*/                   return -negkap;
                }
      
        }
        return -negkap;
}
#undef TFACTR

float revcst(float x[], float y[], int iorder[],
             int ngene,  int ncity, int n[])
{
        float de;
	int i, j;
        int *jorder;
        float *z;
        jorder=ivector(1,ncity);
        z=vector(1,ngene); 
        for (i = 1; i <= ncity; i++) jorder[i] = iorder[i];
        /* original negative value */
        de = -ALEN(x, y, ngene, ncity);   
        reverse(jorder, ncity, n);         
        for (i = 1; i <= ngene; i++) {
           for (j=1;j<=ncity;j++) {
              if ((int)(y[i]) == iorder[j]) z[i] = jorder[j]; 
           }
        }
        /* difference b/w original and rearranged negative kappa value */
	de += ALEN(x, z, ngene, ncity);
	free_ivector(jorder,1,ncity);
        free_vector(z,1,ngene);
        return de;
}

void reverse(int iorder[], int ncity, int n[])
{
        int nn,j,k,l,itmp;
 
        nn=(1+((n[2]-n[1]+ncity) % ncity))/2;
        for (j=1;j<=nn;j++) {
                k=1 + ((n[1]+j-2) % ncity);
                l=1 + ((n[2]-j+ncity) % ncity);
                itmp=iorder[k];
                iorder[k]=iorder[l];
                iorder[l]=itmp;
        }
}

void reversey(float y[], int iorder[], int ncity,int nunits, int n[])
{
        int i,j;
        int *jorder;                         
        float *z;
        jorder=ivector(1,ncity);
        z=vector(1,nunits);
        for (i = 1; i <= ncity; i++) {
            jorder[i] = iorder[i];
        }
        for (i = 1; i <= nunits; i++) {
            z[i] = y[i];
         }
        reverse(iorder, ncity, n);
        for (i = 1; i <= nunits; i++) {
           for (j=1;j<=ncity;j++) {
              if ((int)(z[i]) == jorder[j]) y[i] = iorder[j];
           }                    
         }
        free_ivector(jorder,1,ncity);
        free_vector(z,1,nunits);
}

float trncst(float x[], float y[], int iorder[], 
             int ngene, int ncity, int n[])
{
        float de;
	int i, j;
        int *jorder;
        float *z;
        jorder=ivector(1,ncity);
        z=vector(1,ngene);
        for (i = 1; i <= ncity; i++) {
            jorder[i] = iorder[i];
        }
        /* original negative kappa value */
        de = -ALEN(x, y, ngene, ncity);
        trnspt(jorder, ncity, n);   
        for (i = 1; i <= ngene; i++) {
           for (j=1;j<=ncity;j++) {
              if ((int)(y[i]) == iorder[j]) z[i] = jorder[j];
           }
        }
        /* difference b/w original and rearranged negative value */
        de += ALEN(x, z, ngene, ncity);
        free_ivector(jorder,1,ncity);
        free_vector(z,1,ngene);
	return de;
}


void trnspt(int iorder[], int ncity, int n[])
{
	int m1,m2,m3,nn,j,jj,*jorder;
        int i;
        int *ivector(long nl, long nh);        
        void free_ivector(int *v, long nl, long nh);   
        jorder=ivector(1,ncity);
        n[4]=1 + (n[3] % ncity);
        n[5]=1 + ((n[1]+ncity-2) % ncity);
        n[6]=1 + (n[2] % ncity);
	m1=1 + ((n[2]-n[1]+ncity) % ncity);
	m2=1 + ((n[5]-n[4]+ncity) % ncity);
	m3=1 + ((n[3]-n[6]+ncity) % ncity);
	nn=1;
	for (j=1;j<=m1;j++) {
		jj=1 + ((j+n[1]-2) % ncity);
		jorder[nn++]=iorder[jj];
	}
	for (j=1;j<=m2;j++) {
		jj=1+((j+n[4]-2) % ncity);
		jorder[nn++]=iorder[jj];
	}
	for (j=1;j<=m3;j++) {
		jj=1 + ((j+n[6]-2) % ncity);
		jorder[nn++]=iorder[jj];
	}
	for (j=1;j<=ncity;j++)
           iorder[j]=jorder[j];
	free_ivector(jorder,1,ncity);
}


void trnspty(float y[], int iorder[], int ncity, int nunits, int n[])
{
        int i,j;                  
        int *jorder;
        float *z;        
        jorder=ivector(1,ncity);
        z=vector(1,nunits);
        for (i = 1; i <= ncity; i++) {
            jorder[i] = iorder[i];
        }
        for (i = 1; i <= nunits; i++) {
            z[i] = y[i];
         }                         
        trnspt(iorder, ncity, n);
        for (i = 1; i <= nunits; i++) {
           for (j=1;j<=ncity;j++) {
              if ((int)(z[i]) == jorder[j]) y[i] = iorder[j];
           }
        } 
        free_ivector(jorder,1,ncity);
        free_vector(z,1,nunits);
}


/* allocate a int vector with subscript range v[nl .. nh] */

int *ivector(long nl, long nh)
{
   int *v;
   v = (int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
   if (!v) nrerror("allocation failure in ivector()");
   return v-nl+NR_END;

}

void free_ivector(int *v, long nl, long nh)

{
   free((FREE_ARG) (v+nl-NR_END));
}

int metrop(float de, float t)
{
   float ran1(long *idum);
   static long gljdum = 1;
   return de < 0.0 || ran1(&gljdum) < exp(-de/t);
}

float ran1(long *idum)
{
   const int NTABran1=32;
   int j;
   long k;
   static long iy=0;
   static long iv[32];
   float temp;
   if (*idum <= 0 || !iy) {
      if (-(*idum)<1) *idum=1;
      else *idum = -(*idum);
      for(j=NTABran1+7;j>=0;j--) {
         k=(*idum)/IQran1;
         *idum=IAran1*(*idum-k*IQran1)-IRran1*k;
         if (*idum<0) *idum += IMran1;
         if (j<NTABran1) iv[j] = *idum;
      }
      iy=iv[0];
   }
   k=(*idum)/IQran1;
   *idum=IAran1*(*idum-k*IQran1)-IRran1*k;
   if (*idum<0) *idum += IMran1;
   j=iy/NDIVran1;
   iy=iv[j];
   iv[j] = *idum;
   if ((temp=AMran1*iy)>RNMXran1) return RNMXran1;
   else return temp;
}

int irbit1(unsigned long *iseed)
{
   unsigned long newbit;
   newbit = (*iseed >> 17) & 1
      ^(*iseed >> 4) & 1
      ^(*iseed >> 1) & 1
      ^(*iseed & 1);
   *iseed = (*iseed << 1) | newbit;
   return (int) newbit;
}