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

/*
#include</users/bmatso/medd/paramete.c>
#include</users/bmatso/old467x/screen.inc>
*/

#include<\tc\med\medlib\paramete.c>
#include<\tc\physics\screen.inc>


#define xt 0
#define rt 1



#define simple_simulation 1
#define real_simple_simulation 2
#define realistic_simulation 3
#define jordan_data 4
#define gray_simulation 5


/*	#define L_Max 50   hey dumkopf in paramete.c L = 50 (max poss L) */
#define N 1
#define M 128



#define L_inc 1
#define L_begin 4
#define L_max 7			/*minus one */




	int machine = xt;

	int data_source = realistic_simulation;

	int L;

	double data[N][M];
	double dataX[N][M + L_Max];

	double data_max = 0.0;





main()
{

	int i,j;

printf("simple simulation from modified Cabrelli ex. 1 (1)\n");
printf("real simple simulation from     Cabrelli ex. 1 (2)\n");
printf("realistic simulation of GRP siginal            (3)\n");
printf("jordan data, excerpt from a real data file     (4)\n");
printf("gray simulation; digitized data: Gray ex. 1    (5)\n");
printf("\nenter data source: ");
scanf("%d", &data_source);


	get_data_file(M,&dataX[0][0],data_source);

	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M; j ++)
		{
			data[i][j] = dataX[i][j];
/*			printf("dataX[%d][%d] = %f\n",j,i,dataX[j][i]);*/
		}

/*		stop_screen();*/

	}

	if (machine == xt)
		plot_all(0,M,&data[0][0],
			"MEDD Input, x(i)");


	MEDD_deconvolution();


	if (machine == xt)
		plot_all(0,M,&dataX[0][0],
			"MEDD Output, Y(i) = f * x(i)");

	if (machine == rt)
		file_xchange_data(M, &dataX[0][0], "meddout.dat");


	exit(0);

}


MEDD_deconvolution()
{

	int ickx,icky;

	double C[N][L_Max];
	int i,j,k,s,u;
	double R[L_Max], W[L_Max][L_Max], F[L_Max], Y[N][M + L_Max];
	double expect_Y;
	double D_norm, D_norm_final;
	double medd_filter, medd_output;
	double Y_optimum[N][M + L_Max], E[L_Max];
	double optimum_D_norm, optimum_output;
	int optimum_L = L_begin;
	double f[N][M + L_Max];

	FILE *id_err;



	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L_Max - 1; j ++)
		{
			Y_optimum[i][j] = 0;
		}

	}


	for (L = L_begin; L < L_max; L = L + L_inc)
	{

	D_norm_final = 0;

					/* calculate auto-corelation
					coefficients of input traces */


	for (k = 0; k < L; k ++)
	{

		for (i = 0; i < N; i ++)
		{
			C[i][k] = 0;
			for (j = 0; j < M; j ++)
			{
				if( j + k  >= 0 && j + k < M)
				{

/* printf("C: dataX[i][j]     = dataX[%d][%d] = %f\n",i,j,dataX[i][j]);*/
/* printf("C: dataX[i][j + k] = dataX[%d][%d] = %f\n",i,j+k,dataX[i][j+k]);*/

				    C[i][k] =
				    C[i][k] + dataX[i][j] * dataX[i][j+k];
				}
			}
			printf("C[%d][%d] = %f\n",i,k,C[i][k]);

		}

	}


	for (k = 0; k < L; k ++)
	{
		R[k] = 0;
		for (i = 0; i < N; i ++)
		{
			R[k] = R[k] + C[i] [k];

		 }
	printf("R[%d] = %f\n",k,R[k]);

	}







					/* inversion of autocorrelation
					matrix R[I-J] using
					Levinson Recursion:
					     -1
					W = R   (L x L) - matrix.	*/

	invert_matrix(L,L-1,&R[0],&R[1],&R[1],&W[0][0]);

/*
	for (ickx = 0; ickx < L; ickx++)
	{
		for(icky = 0; icky < L+1; icky++)
		{

	printf("W[%d][%d] = %f\n",ickx,icky,W[ickx][icky]);
		}

	}

	stop_screen();
*/


	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
		{

					/* calculation of [I] [J]
					filter, F	*/

			for (k = 0; k < L; k ++)
			{
				F[k] = 0;
				for (s = 0; s < L; s ++)
				{
					if ( j - s >= 0 && j - s < M)
					{
/* this assumes that j-s is terminative not cyclic*/


/*		printf("dataX[i][j-s] = dataX[%d][%d] = %#g\n",
		i,j-s,dataX[i][j-s]);*/

					F[k] = F[k] + dataX[i][j-s] * W[s][k];
					}

/*	printf("\t\ts,i,j-s,k,W[s][k],dataX[i][j-s],F[k]: %d %d %d %d %f, %f, %f\n",
	s,i,j-s,k,W[s][k],dataX[i][j-s],F[k]);*/

				}
/*				printf("\tFilter F[%d] = %#g\n",k,F[k]);*/
			}

/*stop_screen();
for (icky = 0; icky < M; icky ++)
		printf("before s-k bussiness:	dataX[0][%d] = %#g\n",
		icky,dataX[0][icky]);
stop_screen();*/
					/* calculation of the matrix Y
					corosponding to the output of the
					[I] [J] filter, F	*/
			for (k = 0; k < N; k ++)
			{
				for (s = 0; s < M+L-1; s ++)
				{
					Y[k] [s] = 0;
					for (u = 0; u < L; u ++)
					if (s - u + 1 >= 0 && s - u + 1 < M)
					{
/*printf("dataX[k][s-u+1] = dataX[%d][%d] = %f\n",k,s-u+1,dataX[k][s-u+1]);*/

						Y[k][s] =
					Y[k][s] + F[u] * dataX[k][s-u+1];
					}

/*				printf("\t\tY[%d][%d] = %#g\n",k,s,Y[k][s]);*/
				}
			}
						/* Calculation of ||Y|| */

/*
for (ickx = 0; ickx < M+L-1; ickx ++)
printf("\t\t\tY[%d][%d] = %#g\n",0,ickx,Y[0][ickx]);
*/

			expect_Y = 0;
			for (k = 0; k < N; k ++)
			{
				for (s = 0; s < M + L - 1; s ++)
					expect_Y =
					expect_Y + Y[k][s] * Y[k][s];

/*printf("\t\t\tk = %d, j = %d, M + L = %d, ||Y|| = %#g\n",
				k,j,M+L,expect_Y);*/
			 }
			expect_Y = sqrt(expect_Y);



						/* Calculation of the [I][J]
						coefficients of the [I][J]
						output normalized	*/

			if(expect_Y != 0.0)
				D_norm = (fabs(Y[i][j]))/(expect_Y);
			else
				D_norm = 0;

printf(".");
/*	printf("D_norm = %#g\n",D_norm);*/

			if (D_norm > D_norm_final)
			{
				D_norm_final = D_norm;
/*				medd_filter  = F[?];*/
				medd_output  = expect_Y;

			}


		}
	}


	printf("\n");

	E[L] = 0;
	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
			E[L] = E[L] + Y[i][j] * Y[i][j];
	}

	printf("E[%d] = %#g\n",L,E[L]);


/*	for (ickx = 0; ickx < M+L-1; ickx ++)
printf("Y[%d][%d] = %f\n",0,ickx,Y[0][ickx]);*/

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		if (j < L) f[0][j] = F[j];
	}

	if (machine == xt)
		plot_all(0,M,&f,"MEDD filter values, f");

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		f[0][j] = Y[0][j];
	}


	if (machine == xt)
		plot_all(0,M,&f,"MEDD output for this filter, f");


	printf("D_norm: %f\nmedd_output: %f\n",D_norm_final,medd_output);
	printf("number of filters: %d, error: %#g\n",L,E[L]);

	if (machine != rt) stop_screen();

		 if (L == L_begin) goto label1;

		 if (E[L] <= E[optimum_L])
		 {
label1:		 optimum_L = L;
		 optimum_D_norm = D_norm_final;
		 optimum_output = medd_output;
			for (i = 0; i < N; i ++)
			{
				for (j = 0; j < M + L - 1; j ++)
				Y_optimum[i][j] = Y[i][j];
			}
		}
	}


	i = 0;

	if (machine == rt)
		id_err = fopen("medderr.dat", "wt");

	for (L = L_begin; L < L_max; L = L + L_inc)
	{
		f[0][i] = log10(E[L]);
		if (machine == rt)
			fprintf(id_err, "%e\t%e\n", (double)(L), E[L]);
		i++;
	}

	if (machine == rt)
		fclose(id_err);

	if (machine == xt)
		plot_all(0,i,&f,"Log of Error vs. filter number");


	printf("FINAL RESULTS:\n");
	printf("D_norm: %f\nmedd_output: %f\n",optimum_D_norm,optimum_output);
	printf("Best error[# of filters] = E[%d] = %#g\n",
	optimum_L,E[optimum_L]);

	if (machine != rt) stop_screen();

			for (i = 0; i < N; i ++)
			{
				for (j = 0; j < M + optimum_L - 1; j ++)
				{
					dataX[i][j] = 0;
					dataX[i][j] = Y_optimum[i][j];
				}
			}


}

file_xchange_data(lngth, in_data, filename)
	int lngth;
	double *in_data;
	char *filename[30];
{
	int j;
	FILE *id;

printf("opened: %s, lngth %d\n", filename, lngth);

	id = fopen(filename, "wt");

	for (j = 0; j < lngth; j ++)
		fprintf(id, "%e\t%e\n", (double)(j), *(in_data + j));

	fclose(id);

}

int get_data_file(lngth, in_data, data_source)
	int lngth;
	double *in_data;
	int data_source;
{
	double k0 = 0;
	double *k = &k0;

	FILE *id_trace;
	FILE *id_reflect;
	FILE *id_convolved;

	char argv_rf[5][30] =

		{"med2_rf1.dat",
		"med2_rf2.dat",
		"med2_rf3.dat",
		"med2_rf4.dat",
		"med2_rf5.dat"};

	char argv_cv[5][30] =

		{"med2_cv1.dat",
		"med2_cv2.dat",
		"med2_cv3.dat",
		"med2_cv4.dat",
		"med2_cv5.dat"};

		char arg_rt_info[5][60] =
		{"Reflectivity Series for Example 1: Simple Simulation w/NR.",
		"Reflectivity Series for Example 2: Cabrelli w/NR .",
		"Reflectivity Series for Example 3: Realistic Simulatiow/NR",
		"Reflectivity Series for Jordan Data, Segment 1/wNR.",
		"Reflectivity Series for Gray (Example 5) w/ NR."};

	int ch;
	int j;
	int i;

	if (data_source == simple_simulation)
	{
		j = 0;
	}

	if (data_source == real_simple_simulation)
	{
		j = 1;
	}

	if (data_source == realistic_simulation)
	{
		j = 2;
	}

	if (data_source == jordan_data)
	{
		j = 3;
	}

	if (data_source == gray_simulation)
	{
		j = 4;
	}





	if ((id_reflect = fopen(&argv_rf[j][0],"rt")) == NULL)
	{
		printf("Cannot open input reflection file %d \n", j);
		printf("Error  %s \n",stderr);
		fclose(id_reflect);
		return(0);

	}

	i = 0;
	ch = 2;
	while (ch == 2 && i < lngth)
	{
		ch = fscanf(id_reflect,"%le\t%le", k, (in_data + i));
/*		printf("ch = %d, j = %d, data = %e \n",ch,j,*(in_data + i));*/
		i++;
	}


	fclose(id_reflect);

	if (machine == xt)
		plot_all(0,i,in_data,&arg_rt_info[j][0]);

	if ((id_convolved = fopen(&argv_cv[j][0],"rt")) == NULL)
	{
		printf("Cannot open input file %s \n");
		printf("Error  %s \n",stderr);
		fclose(id_convolved);
		return(0);

	}

	i = 0;
	ch = 2;
	while (ch == 2 && i < lngth)
	{
		ch = fscanf(id_convolved,"%le\t%le", k, (in_data + i));
/*		printf("ch = %d, j = %d, data = %e \n",ch,j,*(in_data + i));*/
		i++;
	}

end_of_file:

	fclose(id_convolved);

	printf("number of data points inputted: %u \n",i-1);

	return(i-1);



}

/*

plot_all(i, j, data, name)
	int i;
	int j;
	double *data;
	char *name;
{

}

*/