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

#define	MAX_STRLEN	 16
#define MAX_VARS	 32
#define MAX_OBSERVATIONS 365

typedef struct {

	float data[MAX_VARS];

} observation_t;

int  		nVariables = 0;
int  		nObservations = 0;

char 		VarNames[MAX_VARS][MAX_STRLEN];

observation_t	Observations[MAX_OBSERVATIONS][MAX_OBSERVATIONS]; 

float mincorrel = 0;

int bSumItems = 0;

float  GetAverage( observation_t *obs, int num, int column )
{

	int i;

	float val, mean = 0.0f;


	for ( i = 0 ; i < num ; i++ )
	{

		val = obs[ i ].data[column ];

		mean += val;

	}

	if ( i > 0 )
	{ 
		mean /= (float)( i );
	}
	
	
	return mean;
}


void PrintSamples( int level )
{

	int i, j;

	float val, sum;

	// printf("Data level %d \n\n", level);
	printf("\n");

	for( i = 0 ; i < nVariables ; i++ )
	{
		printf( "%s\t" , VarNames[i] );
	}
	printf("| SUM\n"); 

	for ( i = 0 ; i <= nVariables ; i++ )
	{
		printf("========");
	}
	printf("\n");

	for( i = 0 ; i < nObservations ; i++ )
	{
		printf( "%d\t", (int) Observations[ level ][ i ].data[0] );

		for( j = 1, sum = 0 ; j < nVariables ; j++ )
		{
			val = Observations[ level ][ i ].data[ j ] ;

			if ( val >= 0 ) printf( " " );
			
			printf( "%.2f \t", val );

			sum += val;
		}
		
		printf("| %.2f \n", sum );
	}
	
	for ( i = 0 ; i <= nVariables ; i++ )
	{
		printf("--------");
	}
	printf("\n%d\t", nObservations );
	for ( i = 1 , sum = 0 ; i < nVariables ; i++ )
	{
		val = GetAverage( Observations[0] , nObservations, i );
		if ( val < 0 )
			printf("%3.2f\t",  val);
		else
			printf(" %3.2f\t",  val);
		sum += val;
	}
	printf("| %.2f\n", sum );
}

int LoadSamples( char *filename )
{

	int i, dim;
	int lines, samples = 0;

	FILE *file;

	char buf[80];

	float *pData, sum = 0;


	file = fopen( filename, "ro" );

	for ( lines = 0 ; fgets( buf, 256, file ) ; lines++ )
	{
		if( lines == 0 )
		{
			nVariables = sscanf( buf, "%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s", 

				VarNames[0], VarNames[1], VarNames[2], VarNames[3], 
				VarNames[4], VarNames[5], VarNames[6], VarNames[7],
				VarNames[8], VarNames[9], VarNames[10], VarNames[11], 
				VarNames[12], VarNames[13], VarNames[14], VarNames[15] 	 	
			);
		}
		else
		{
			if( buf[0] != '#' )
			{			
				pData = (float *) &Observations[0][ nObservations + samples ] ;

				dim = sscanf( buf, "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",
					&pData[0], &pData[1], &pData[2], &pData[3], 
					&pData[4], &pData[5], &pData[6], &pData[7],
					&pData[8], &pData[9], &pData[10], &pData[11], 
					&pData[12], &pData[13], &pData[14], &pData[15]
				);
			
				if( dim > 0 )
				{
					samples++;					
				}
			}
		}
	}
	
	fclose( file );

	nObservations += samples;

	return samples;
}


float  GetVariance( observation_t *obs, int num, int mean, int column )
{

	int i;

	float val, var = 0.0f;


	for ( i = 0 ; i < num ; i++ )
	{

		val = ( obs[ i ].data[column ] - mean );

		var += val * val;

	}
		
	return var ;
}

float  GetMinMax( observation_t *obs, int num, float *minimum, float *maximum, int column )
{

	int i;

	float val, min = 1000000, max = -1000000;


	for ( i = 0 ; i < num ; i++ )
	{

		val = obs[ i ].data[column ];

		if ( val < min ) min = val;
		if ( val > max ) max = val; 
	}
		
	*minimum = min;
	*maximum = max;
}

float DeleteAverage( observation_t *obs, int num, int avg, int column )
{

	int i;

	for ( i = 0 ; i < num ; i++ )
	{
		obs[ i ].data[ column ] -= avg ;
	}

	return GetAverage( obs, num, column );
}


void ComputeLevel( int lag )
{

	int i, j;

	float diff;

	for ( i = 0 ; i < nObservations ; i++ )
	{
		Observations[ lag ][ i ].data[ 0 ] = Observations[ lag - 1 ][ i ].data[ 0 ] ; 
	}

	for( i = lag ; i < nObservations ; i++ )
	{
		for( j = 1 ; j < nVariables ; j++ )
		{
			diff  = Observations[ lag - 1 ][ i ].data[ j ];
			diff -= Observations[ lag - 1 ][ i - 1 ].data[ j ];
			Observations[ lag ][ i ].data[ j ] = diff; 
		}
	}
}
 
float GetCorrelation( int lag, int column )
{
	int t;

	float sum = 0.0f;

	float avg = GetAverage( Observations[0], nObservations, column );

	for( t = 0 ; t < nObservations - lag ; t++ )
	{
		sum += ( Observations[0][t].data[column] - avg ) * ( Observations[0][t+lag].data[column] - avg );
	}
	

	sum   /= ( nObservations - lag ) * sqrt ( GetVariance( Observations[0], nObservations, avg, column ) ) ; 

	return sum ; 
	
}

float GetAutoCorrelation( int lag, int column )
{
	int t;

	float num = 0.0f;
	float denum = 0.0f;

	float avg = GetAverage( Observations[0], nObservations, column );

	for( t = 0 ; t < nObservations - lag ; t++ )
	{
		num += ( Observations[0][t].data[column] - avg ) * ( Observations[0][t+lag].data[column] - avg );
	}
	
	for( t = 0 ; t < nObservations ; t++ )
	{
		denum += pow( Observations[0][t].data[column] - avg , 2.0f);
	}

	num   /= nObservations - lag; 
	denum /= nObservations - 1;

	return num / denum; 
	
}

float GetInterCorrelation( int lag, int column1, int column2 )
{
	int t;

	float num = 0.0f;
	float denum = 0.0f;


	float avg1 = GetAverage( Observations[0], nObservations, column1 );
	float avg2 = GetAverage( Observations[0], nObservations, column2 );

	for( t = 0 ; t < nObservations - lag ; t++ )
	{
		num += ( Observations[0][t].data[column1] - avg1 ) * ( Observations[0][t+lag].data[column2] - avg2 );
	}
	
	for( t = 0 ; t < nObservations ; t++ )
	{
		denum += (Observations[0][t].data[column1] - avg1) * (Observations[0][t].data[column2] - avg2); 
		
	}

	num   /= nObservations - lag; 
	denum /= nObservations - 1;

	return num / denum; 
	
}

int main( int argc,  char **argv )
{

	int i, j, k,l;

	float  avg, mean, moy, scale, idx, correl, sum, variance;
	float  min = 1000.0f, max = -1000.0f;

	printf("\n\nACF v0.1 by Cyclone \n\n");

	for( i = 1 ; i < argc; i++ )
	{
		if( argv[i][0] == '-' )
		{
			switch( argv[i][1] )
			{
				case 'a' : 
					sscanf( argv[i] + 2,  "%f", &mincorrel );
					printf("\n\nMinimum correlation set to %f \n", mincorrel); 
					break;

				case 's' :
					bSumItems = 1;
					break;
			} 
		}
		else
		{
			(void) LoadSamples( argv[i] );
		}	
	}

	PrintSamples(0);

/*	
	(void) DeleteAverage( Observations[0], nObservations, avg, 1);		
	printf("\nNew %s centered at %.2f \n", VarNames[ 1 ] , avg );	
	PrintSamples(0);
*/
	// printf("Compute lags : \n");
	for( i = 1 ; i < nObservations /* NUM_LAG */ ; i++ )
	{
		// printf("\nCompute the lag %d correlation \n", i);
		ComputeLevel(i);
		// PrintSamples(i);
		// printf(".");
	}
	// printf("\n");

	// INTRA CORRELATION
	for( j = 1 ; j < nVariables ; j++ )
	{
		printf("\n\nACF of the variable %s : ", VarNames[j] );

		avg = GetAverage( Observations[0], nObservations, j );	
		// printf("\nAverage of %s is %.2f \n\n", VarNames[ j ], avg );

		variance = GetVariance( Observations[0], nObservations, avg, j );
		// printf("\nVariance of %s is %.2f \n\n", VarNames[ j ],  sqrt(variance / nObservations));	

		GetMinMax( Observations[0], nObservations, &min, &max, j);
  
		printf("\tmin/avg/var/max = %.2f/%.2f/%.3f/%.2f \n\n", 
			min, avg, (float)(sqrt(variance)/nObservations), max); 

		for( i = 1 ; i < nObservations ; i++ )
		{
			correl = GetAutoCorrelation( i , j );

			// limit the display to "bigs" autocorrelation's coefs if wanted
			if( mincorrel != 0.0f )
			{
				if ( mincorrel < 0.0f )
				{
					if ( fabs(correl) < fabs(mincorrel) )	
						continue;
				}
				else
				{
					if ( correl < mincorrel )
						continue;
				}
			}

			// print the autocorrelation coef at this lag
			if( correl < 0 )
				printf("Lag = %3d Correl = %.5f ", i, correl ); 
			else
				printf("Lag = %3d Correl =  %.5f ", i, correl ); 
		
			// display it
			correl *= 100.0f;

			if( correl < 0 )
			{
				correl =  40.0f + correl; 

				for( l = 0 ; l < correl ; l++ )
				{
					printf( " " );
				}
				for( /* l = idx */ ; l < 40 ; l++ )
				{
					printf( "*" );
				}
				printf("|\n"); 		
			}
			else
			{
				for( l = 0 ; l < 40 ; l++ )
				{
					printf( " " );
				}
				printf( "|");
				for( l = 0 ; l < correl ; l++ )
				{
					printf( "*" );
				}
				printf("\n"); 		
			}
		} 
	}

	printf("\n\n\t\t====================================================================================================\n");

	// INTER CORRELATION
	for( j = 1 ; j < nVariables ; j++ )
	for( k = 1 ; k < nVariables ; k++ )
	{
		printf("\n\nIntercorrelation between variables %s and %s : \n\n", VarNames[j], VarNames[k] );

		avg = GetAverage( Observations[0], nObservations, j );	
		variance = GetVariance( Observations[0], nObservations, avg, j );
		GetMinMax( Observations[0], nObservations, &min, &max, j);
		printf("\tmin/avg/var/max of %s = %.2f/%.2f/%.3f/%.2f \n\n", 
			VarNames[j], min, avg, (float)(sqrt(variance)/nObservations), max); 

		avg = GetAverage( Observations[0], nObservations, k );	
		variance = GetVariance( Observations[0], nObservations, avg, k );
		GetMinMax( Observations[0], nObservations, &min, &max, k);
		printf("\tmin/avg/var/max of %s = %.2f/%.2f/%.3f/%.2f \n\n", 
			VarNames[k], min, avg, (float)(sqrt(variance)/nObservations), max); 

		for( i = 1 ; i < nObservations ; i++ )
		{
			correl = GetInterCorrelation( i , j, k);

			// limit the display to "bigs" autocorrelation's coefs if wanted
			if ( mincorrel != 0.0f )
			{
				if ( mincorrel < 0.0f )
				{
					if ( fabs(correl) < fabs(mincorrel) )	
						continue;
				}
				else
				{
					if ( correl < mincorrel )
						continue;
				}
			}

			// print the intercorrelation coef at this lag
			if( correl < 0 )
				printf("Lag = %3d Correl = %.5f ", i, correl ); 
			else
				printf("Lag = %3d Correl =  %.5f ", i, correl ); 
		
			// display it
			correl *= 100.0f;

			if( correl < 0 )
			{
				correl =  40.0f + correl; 

				for( l = 0 ; l < correl ; l++ )
				{
					printf( " " );
				}
				for( /* l = idx */ ; l < 40 ; l++ )
				{
					printf( "*" );
				}
				printf("|\n"); 		
			}
			else
			{
				for( l = 0 ; l < 40 ; l++ )
				{
					printf( " " );
				}
				printf( "|");
				for( l = 0 ; l < correl ; l++ )
				{
					printf( "*" );
				}
				printf("\n"); 		
			}
		} 
	}

	printf("\n\n");

	return 0;
}