/*
This program carries out the "nearest neighbor test" described in 
Hudson, R.R. 2000  "A new statistic for detecting genetic differentiation"
Genetics (in press).
Data consists of a matrix of pairwise differences between sampled sequences,
as described in readme. An example data file is sample.dat .  The program
can usually be compiled by gcc -o snn snn.c  . And then run by 

snn y 1000 20 27 <sample.dat >sample.out   .

On the command line,  20 is the sample size from
locality one, and 27 is the sample size from locality two.  1000 is
the number of permutations to carry out.  y indicates that the input
data matrix should be copied to the output.
*/


#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <ctype.h>




main( int argc, char * argv[] )
{
	int nperms, nlocs,  nsam ;
	double **dij, pval, snno,   permt(), **tdij ;
	int i, loc, wf,  *ni, tni[3], loc1, loc2, j, sti, stj, endi, endj ;
	 
	if( argc < 4 ) {
	  printf(
 "usage: snn  y/n(print data?) n_permutations n1 n2  \n"); 
		exit(1) ;
		}
	nperms = atoi( argv[2] ) ;
	nlocs = 2 ;
	if( (ni = (int *)malloc( (size_t) (nlocs+1)*sizeof(int) )) == NULL)
		perror( "malloc error1\n") ;
	nsam = 0 ;
	for( i=0; i<nlocs; i++){
	   ni[i] = atoi( argv[3+i] );
	   nsam += ni[i];
	   }

	if( (dij = (double **)malloc((size_t)nsam*sizeof( double * ) ) ) == NULL)
		perror( "malloc error2\n") ;
	for(i=0;i<nsam;i++){
		if( (dij[i] = (double *)malloc( (size_t)nsam*sizeof(double))) == NULL)
			perror( "malloc error3\n") ;
	   }
  

	getmatrix(nsam, dij, argv[1][0] );
	
	fprintf(stdout,"\nSample configuration: ");
	loc=0;
	while( loc<nlocs) fprintf(stdout,"%d  ",ni[loc++]);

	fprintf(stdout,"\nNumber of permutations: %d  \n",nperms);

	pval = permt(nperms,nsam,nlocs,ni,dij, &snno ) ;

	fprintf(stdout,"Observed values of Snn: ");
	fprintf(stdout,"  %lf ( p-value: %lf)\n\n", snno, pval);
	
}
	
	double
permt(nperms,nsam,nlocs,ni,dij,psnno )
   int nperms,nlocs,nsam ;
   int *ni ;
   double **dij, *psnno ;
  {
  	int i, count, count2,  loc, start ;
  	double  eps = 1.0e-8, snno, snnf(), snn  ;
  	int *ran_ind ;

	if( (ran_ind = (int *) malloc( (size_t)nsam*sizeof(int) )) == NULL)
			perror( "malloc error9\n") ;

	for(i=0;i<nsam;i++)  ran_ind[i]=i;
	snno = *psnno = snnf( nsam, nlocs, ni,dij,ran_ind);
	count = count2 =  0 ;	
while (nperms - count++) {
		scramb(ran_ind,nsam);
		
	snn  = snnf( nsam, nlocs, ni,dij,ran_ind);

	if( snn <= snno ) count2++ ;

	}
   free( ran_ind);
   return( (double)count2/(double)nperms );	
	
}


	int	   
getmatrix(nsam,dij,pflag)
	int nsam;
	double **dij;
	char pflag;
{
	int i, j, jstart, jend ;
	char gamnam[35], s[1000];
	FILE *pf ;
	
	pf = stdin; 

	i = 0 ;
	jend = -1 ;
	while( i< nsam-1 ){
	  jstart = jend +1  ;
	  while(1){
		fgets(s,1000, pf);
		if( pflag == 'y') fputs(s,stdout);
		if( s[0]=='O' && s[1]=='T' && s[2]=='U'&&s[3]=='s' ) {
			jend = countwords(s) + jstart -2 ;
			if( i==0) jstart++ ;
			break;
			}
		}	
	for( i=0; i<jend; i++){
		fscanf(pf," %s",gamnam);
		if( pflag == 'y') fprintf(stdout,"%s",gamnam);
		dij[i][i] = 0.0 ;
		if( jstart < i+1) jstart = i+1 ;
		for( j=jstart; j<= jend; j++) {
			fscanf(pf," %lf", dij[i]+j );
			if( pflag == 'y') fprintf(stdout," %g",dij[i][j] );
			 dij[j][i]=dij[i][j];
			 }
		if( pflag == 'y' ) fprintf(stdout,"\n");
		}
	}
}

        int
countwords(s)
        char *s;
{
        int wc, cc=0, out = 1, in=0 , state, c ;

        wc = 0 ;
        state = out;
        while( s[cc] != '\n' ){
                c = s[cc] ;
                if( isspace( c) ) state = out ;
                else if( state == out) {
                        state = in ;
                        wc++;
                        }
                cc++;
                }
                return( wc) ;
}

	int
scramb(ran_vect, n)
	int *ran_vect, n;
{
	double ran1();
	int i, j, temp ;
	
	for( i = n-1 ; i> 0 ; i--){
		j = ran1()*(i + 1) ;
		temp = ran_vect[j];
		ran_vect[j] = ran_vect[i] ;
		ran_vect[i] = temp ;
		}
}



	double
ran1()
{
  double drand48();
	return( (double)drand48()  );
}


double snnf(int nsam, int npop, int *config, int **dij, int *ran_ind)
{
        int i, ind, mall, nmall, nmwithin , pop, indstart ;
        double p ;

        p = 0. ;

        for(ind= pop=0; pop<npop; pop++)
          for( i=0, indstart= ind ; i<config[pop]; i++,ind++){
                mall = minofdij( nsam, dij[ ran_ind[ind]],ran_ind[ind], &nmall );
                nmwithin = howman(nsam, dij[ ran_ind[ind]],ran_ind, ind,mall, indstart, config[pop] );
                p += (double)nmwithin/(double)nmall ;
                }
         return( p /= nsam ) ;

}

/* counts number of nearest neighbors that are within subpopulation */

int  howman(int nsam, int *vec, int *ran_ind, int ind, int mall, int indstart, int nsubpop)
{
        int count, i ;

     count = 0 ;
        for( i=indstart; i<indstart + nsubpop; i++)
                if( (i != ind) && ( vec[ ran_ind[i] ] == mall ) ) count++;
return( count ) ;
}

/* finds minimum distance to neighbor  and counts number of such nearest neighbors */

        int
minofdij(int nsam, int *vec, int ind, int *pnmall)
{
        int i, mind, count ;

        if( ind != 0 ) mind = vec[0] ;
        else mind = vec[1] ;
        count = 0 ;
        for( i=0; i<nsam; i++)
                if( i != ind ){
                         if( vec[i] < mind ) { count =1; mind = vec[i] ;}
                         else if ( vec[i] == mind ) count++ ;
                         }
        *pnmall = count ;
        return( mind ) ;
}