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

Copyright (C) 2000 Peter Clote. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

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


/******************************************

	entropyPlot0.c

	P. Clote, 7 Sept 1998
	Program computes entropy (single and dinucleotide)
	This program assumes complete data -- if there is
	incomplete data (eg R,Y, etc), then there is no
	contribution to nucleotide/dinucleotide count.
	Nevertheless, since the denominator is window size
	W (for nucleotides) and W-1 (for dinucleotides),
	this leads to small values if many data items are
	incomplete.
	In a complicated fashion, this is corrected in
	entropy1.c


	WARNING - This runs on my Tecra laptop running Linux.
	You may get a runtime error on another machine, 
	because the main loop is non-terminating (while (1)). 
	In entropy1.c this was corrected.

	

******************************************/

#include <math.h>
#include <stdio.h>
#define W 100		/* window size */
#define lg(x) ( ((x) == 0) ? 0 : (log(x)/log(2)) )

main(int argc, char *argv[])
{
	FILE *in, *out1, *out2;    
		/* out1 is entropy plot for single nucleotides */
		/* out2 is entropy plot for dinucleotides      */
	int n=0,i,j,a=0,c=0,g=0,t=0;  	
		/* a,c,g,t nucleotides */
    	int num = 0,s=0,k=0,w=0,r=0,y=0; 
		/* num nucleotides in file */
	char ch;
	char win[W];	/* window */
	int first,second,index1,index2;
	double ent1, ent2;		/*** entropy1 and entropy2 */
	double fa,fc,fg,ft;		/*** frequency of a,c,g,t */
	int pair[4][4];  		/*** for dinucleotide counts */
	double freq[4][4];  		/*** for dinucleotide frequencies */
	void error(char *);		/** error handling */

	/**** Error handling for input file ****/	
	if (argc != 4)
		{
	     	fprintf(stderr,"%s in out1 out2 \n",argv[0]);
		exit(1);
		}
	/***** Initialization of I/O files and pair array **/
	in = fopen(argv[1],"r");
	out1 = fopen(argv[2],"w");
	out2 = fopen(argv[3],"w");
	for (i=0;i<4;i++)
		for (j=0;j<4;j++) pair[i][j] = 0;

	/** remove first line header from M. jannaschii genome file **/
	while ( ( ch = fgetc(in)) != '\n' ); /** remove first line **/		

	/* fill up window. Assume file contains at least W char */
	i=0;
	while (i<W){
		ch=fgetc(in);
		if (ch == '\n') continue;
		win[i++]=ch;	
	}

	/** process first 2 char in window  **/
	/** assume first 2 char are among A,C,G,T **/
	first = win[0]; num++; 
	switch ( first )
		{
		case 'A' : a++; index1 = 0; break;
		case 'C' : c++; index1 = 1; break;
		case 'G' : g++; index1 = 2; break;
		case 'T' : t++; index1 = 3; break;
		default : error("Bad Data"); break;
		}
	second = win[1]; num++;
	switch ( second )
		{
		case 'A' : a++; index2 = 0; break;
		case 'C' : c++; index2 = 1; break;
		case 'G' : g++; index2 = 2; break;
		case 'T' : t++; index2 = 3; break;
		default : error("Bad Data"); break;
		}
	if ( (index1 < 4) && (index2 < 4) ) 
		{
		pair[index1][index2]++;
		num++;
		}

	/** process next W-2 char in file to read first window  **/
	/** assumes EOF does not come before first W characters **/
	for (i=2;i<W;i++) 
		{
		ch = win[i];
		index1 = index2; /* move index down */
		switch ( ch )
			{
			case 'A' : a++; index2 = 0; break;
			case 'C' : c++; index2 = 1; break;
			case 'G' : g++; index2 = 2; break;
			case 'T' : t++; index2 = 3; break;
			case 'N' : n++; index2 = 4; break;
			case 'S' : s++; index2 = 5; break;
			case 'K' : k++; index2 = 6; break;
			case 'W' : w++; index2 = 7; break;
			case 'R' : r++; index2 = 8; break;
			case 'Y' : y++; index2 = 9; break;
			default : printf("%c\n",second); error("Bad Data");  break;
			}
		if ( (index1 < 4) && (index2 < 4) ) 
			{
			pair[index1][index2]++;
			num++;
			}
			
		} 

	/*** compute nucleotide frequencies in window **/
	fa = (double) a / (double) W;
	fc = (double) c / (double) W;
	fg = (double) g / (double) W;
	ft = (double) t / (double) W;

	ent1 = -( (fa*lg(fa)) + (fc*lg(fc)) + (fg*lg(fg)) + (ft*lg(ft)) );

	/*** compute dinucleotide frequencies and entropy2 in window **/
	ent2 = 0;
	for (i=0;i<4;i++)
		for(j=0;j<4;j++)
			{
			freq[i][j] = (double) pair[i][j] / (double) (W-1);
			ent2 -= freq[i][j]*lg(freq[i][j]);
			}

	
	/*** output entropy1 and entropy2 */
	fprintf(out1,"%.10lf \n",ent1);
	fprintf(out2,"%.10lf \n",ent2);
		

	/** continue processing **/	
	do
		{
		/* remove contribution for nucleotide freq for first char */
		switch ( win[0] )
			{
			case 'A' :  a--; index1 = 0; break;
			case 'C' :  c--; index1 = 1; break;
			case 'G' :  g--; index1 = 2; break;
			case 'T' :  t--; index1 = 3; break;
			default :  break;
			}
		switch ( win[1] )
			{
			case 'A' :  index2 = 0; break;
			case 'C' :  index2 = 1; break;
			case 'G' :  index2 = 2; break;
			case 'T' :  index2 = 3; break;
			default :  break;
			}

		/* remove contribution for nucleotide freq for first char */
		if ( (index1<4) && (index2<4) )
			pair[index1][index2]--;

		/* move things down in window */
		for (i=0;i<W-1;i++) win[i]=win[i+1];

		ch = fgetc(in);
		if ( ch == EOF ) break;
		if ( ch == '\n' ) 
			{
			ch = fgetc(in);
			if ( ch == EOF ) break;
			}
		/**** assumes not 2 '\n' in a row  **/
		win[W-1]=ch;

		switch ( win[W-2] )
			{
			case 'A' : index1 = 0; break;
			case 'C' : index1 = 1; break;
			case 'G' : index1 = 2; break;
			case 'T' : index1 = 3; break;
			default :  break;
			}

		switch ( ch )
			{
			case 'A' : a++; index2 = 0; break;
			case 'C' : c++; index2 = 1; break;
			case 'G' : g++; index2 = 2; break;
			case 'T' : t++; index2 = 3; break;
			case 'N' : n++; index2 = 4; break;
			case 'S' : s++; index2 = 5; break;
			case 'K' : k++; index2 = 6; break;
			case 'W' : w++; index2 = 7; break;
			case 'R' : r++; index2 = 8; break;
			case 'Y' : y++; index2 = 9; break;
			default : printf("%c\n",second); error("Bad Data");  break;
			}
		num++;
		if ( (index1 < 4) && (index2 < 4) ) 
			pair[index1][index2]++;

		/*** compute nucleotide frequencies in window **/
		fa = (double) a / (double) W;
		fc = (double) c / (double) W;
		fg = (double) g / (double) W;
		ft = (double) t / (double) W;
		ent1 = -( fa*lg(fa) + fc*lg(fc) + fg*lg(fg) + ft*lg(ft) );
	
		/*** compute dinucleotide frequencies and entropy2 in window **/
		ent2 = 0;
		for (i=0;i<4;i++)
			for(j=0;j<4;j++)
				{
				freq[i][j] = (double) pair[i][j] / (double) (W-1);
				ent2 -= freq[i][j]*lg(freq[i][j]);
				}
		
		/*** output entropy1 and entropy2 */
		fprintf(out1,"%.10lf \n",ent1);
		fprintf(out2,"%.10lf \n",ent2);
	
	} while ( 1 );		/** infinite while loop */

}


void error(char * s)
{
	fprintf(stderr,"%s\n",s);
	exit(1);
}