/*	NOTES:

		21-SEP-90 -- discovered that monitoring the D-norm value
		that it may jump down suddenly, this is due to the walking
		off of a reflection peak off the end, note that the output
		is shifted by length of the filter.  As this length gets
		longer, real reflection peaks can walk off the end.
		This occured for (2) at 42 filters.

		22-SEP-90 -- modified to accept JORDAN data.
		25-SEP-90 -- modified to accept multip[le scans/input
		26-SEP-90 -- modified plot_all to do specific scaling
		12-JUN-94 -- modified to work with C++
*/

// System Include Files
#include<math.h>
#include<values.h>
#include <sys\stat.h>

// Application Include Files
#include "winlib.h"
#include "paramete.h"
#include "meddprm.h"
#include "general.h"
#include "inverse.h"
#include "files.h"
#include "manipult.h"
#include "geofiles.h"
//#include "ddf_util.h"


//  Sets automatic scaling on all plots
#define scale_min 0
#define scale_max 0

//  L_Max is defined in paramete.h
#define N 1				// Number of Scans
#define M 250				// Number of Samples per Scan
//#define M 96


#define part_of_M 21

// L is the filter length which is determined by the maximization of the
// norm.  The number of filters is determined in a loop which begins with
	
#define L_begin 17	      // the beginning number of filters
#define L_max 	18     	      // the last number of filters <minus one>
#define L_inc 4               // the number of filters to increment by in the loop.


//#define beginning_scan 100
#define beginning_scan 1
#define ending_scan 616
//#define ending_scan 500

//#define output_starts_at L
#define output_starts_at 0

#define input_front_exponential_termination    FALSE
#define input_front_exponential_damping_factor  (int)(M * .05)
#define input_back_exponential_termination     FALSE
#define input_back_exponential_damping_factor   (int)(M * .5)
#define output_front_exponential_termination   FALSE
#define output_front_exponential_damping_factor (int)(M * .01)
#define output_back_exponential_termination    FALSE
#define output_back_exponential_damping_factor  (int)( (double)M * .5)



#define plot_nothing 0     	// define plot level
#define plot_only_end_results 1
#define plot_each_iteration 2
#define plot_algorithm_partial 3
#define plot_algorithm_tracking 4
#define plot_everything 5
#define debug1 6






	char input_filename[44];
	char output_filename[44];
	char reflect_filename[44];

	int plot_query = plot_nothing;
//	int plot_query = plot_algorithm_tracking;
//	int plot_query = plot_algorithm_partial;
//	int plot_query = plot_only_end_results;
	int save_error_file = TRUE;

	int data_source = realistic_simulation;
	unsigned long int scan_number;

	int L;

	double data[N][M];
	double dataX[N][M + L_Max];
        double reflections[M];
	double Y_optimum[N][M + L_Max];
	double Y[N][M + L_Max];
	double f[N][M + L_Max];
	double C[N][L_Max];
	double R[L_Max], W[L_Max][L_Max], F[L_Max];
	double D_norm_array[L_Max];

	double data_max = 0.0;



	void MEDD_deconvolution(HWND hWnd);



void medd(HWND hWnd)
{
        #include "common.h"
	#include <FCNTL.H>



	int i,j;
	int l;
        double data_local_max, data_local_min;
	double ave;
	double max_of_input_data;

	FILE *id_medded;
        int id_binary;

	if (chdir("c:\\jordan\\data") != 0)
	{
		sprintf(print_buffer,"Error, cannot change directory\nfrom medd");
		stop_screen(hWnd, print_buffer);

	}


//  Specify input filename

	reflect_filename[0] = 0;
    	strcpy( reflect_filename, "reflect.001" );

	input_filename[0] = 0;
//    	strcpy( input_filename, "jord34.003" );
// 	strcpy( input_filename, "jnk03.dat" );
 	strcpy( input_filename, "reflect.001" );

// Delete and re-create any output files of the same name

    	output_filename[0] = 0;
	strcpy( output_filename, "jord34.002" );

//	ddf_name[0] = 0;
//	strcpy( ddf_name, "jord34.ddf" );


//read_and_plot_data_file(hWnd,1, 10, M, output_filename);


	_fmode = O_BINARY;
	if (id_binary = creat( output_filename, S_IWRITE ) == -1)
	{
		sprintf(print_buffer,"Error, cannot delete binary file %s:\n %s\n%d\n", output_filename, sys_errlist[errno], id_binary);
		stop_screen(hWnd, print_buffer);
	}
////	close(id_binary);








// Begin reading in scans, processing then and saving the processed data
// to disk

/* if( read_ddf() == -1 )
    {
	sprintf(print_buffer,"Descriptor file not found");
	stop_screen(hWnd, print_buffer);
    }

	sprintf(print_buffer,"%s",hdr);
	stop_screen(hWnd, print_buffer);
	sprintf(print_buffer,"%s",list[0]);
	stop_screen(hWnd, print_buffer);
 */

for (scan_number = beginning_scan - 1; scan_number <= ending_scan - 1; scan_number ++)
{
//   data_source = simple_simulation;
//   data_source = real_simple_simulation;
//   data_source = realistic_simulation;
//   data_source = jordan_data;
//   data_source = gray_simulation;
   data_source = binary_simulation_1;



   if (scan_number == 1)
   {

//	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);


	if (data_source == 4)
	{
//		printf("\nenter starting scan number: ");
//		scanf("%d", &scan_number);
//

// May want unique scaling for this selection:
//
//		scale_min = 100;
//		scale_max = 150;

	}
   }                      

if (plot_query >= plot_everything)
{
sprintf(print_buffer,"Data Source Selected: %d\nscan #%d",data_source,scan_number);
stop_screen(hWnd, print_buffer);
}

   switch (data_source)
   {
	case jordan_data:
		get_array_from_disk(hWnd, scan_number, M, &data[0][0], input_filename, binary_file);
                break;
	case simple_simulation:
		get_data_file(hWnd,N,M,&data[0][0],data_source,scan_number);
		break;
	case real_simple_simulation:
		get_data_file(hWnd,N,M,&data[0][0],data_source,scan_number);
		break;
	case realistic_simulation:
		get_data_file(hWnd,N,M,&data[0][0],data_source,scan_number);
		break;
	case gray_simulation:
		get_data_file(hWnd,N,M,&data[0][0],data_source,scan_number);
		break;
	case binary_simulation_1:
		get_array_from_disk(hWnd, scan_number, M, &data[0][0], input_filename, binary_file);
		get_array_from_disk(hWnd, scan_number, M, &reflections[0], reflect_filename, binary_file);
		break;
	default:
		sprintf(print_buffer,"Data Source Undefined");
		stop_screen(hWnd, print_buffer);
	}



	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M; j ++)
		{
			dataX[i][j] = data[i][j];

//debug			dataX[i][j] = .6  + 0*exp((double)(j)/170.0 ) + .03 *sin((double)j/100.0 * 16.0 * pi * scan_number * 24);

		}

	}

	if (machine == tower && plot_query >= plot_only_end_results)
        {


		sprintf(print_buffer, "MEDD Input, x(i)");
		plot_all(hWnd,-N,M,&data[0][0],print_buffer,
			scale_min,scale_max);


	}



/*
	if (machine == tower && plot_query >= plot_only_end_results)
        {
		sprintf(print_buffer, "MEDD Input, x(i)");
		plot_all(hWnd,-N,M,&dataX[0][0],print_buffer,
			scale_min,scale_max);


	}
*/
//  Force a exponential taper onto the front end of the data (if requested)
	if (input_front_exponential_termination)
		gaussian_termination_of_array(M, &dataX[0][0], FRONT_END, input_front_exponential_damping_factor);

//  Force a exponential taper onto the front end of the data (if requested)
//  Now MEDD will produce no output if the last elements of the data array
//  are zero, so after damping any reflection signals from the end,
//  a slowly varying (expenential) offset is added to prevent MEDD
//  from giving a zero output.
	if (input_back_exponential_termination)
             {
		gaussian_termination_of_array(M, &dataX[0][0], BACK_END,  input_back_exponential_damping_factor );

		sprintf(print_buffer, "MEDD Input w/ taper, x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][0],print_buffer,
			scale_min,scale_max);

		max_of_input_data = maximum(
			absolute_value( maximum_of_array(M, &dataX[0][0]) ),
			absolute_value( minimum_of_array(M, &dataX[0][0]) ));


		for (i = 0; i < N; i ++)
//			for (j = 0; j < M; j ++)
			for (j = M-input_back_exponential_damping_factor; j < M; j ++)
	dataX[i][j] = dataX[i][j] + max_of_input_data * .1 * exp((double)(j)/(double)M - 1) ;
//	dataX[i][j] = dataX[i][j] + max_of_input_data * .01;

		sprintf(print_buffer, "MEDD Input w/ non-zero termination, x(i)");
	        if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][0],print_buffer,
			scale_min,scale_max);


	    }




//	MEDD_deconvolution(hWnd);



		sprintf(print_buffer, "MEDD Output before, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);


//  remove any offsets about the average
		ave = average_of_array(M, &dataX[0][output_starts_at]);


		for (j = 0; j < M+L; j ++)
		{
			dataX[0][j] = dataX[0][j] - ave;
		}
		sprintf(print_buffer, "MEDD Output ave'd, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);


//  Force a exponential taper onto the front and back end of the data (if requested)
	if (output_front_exponential_termination)
		gaussian_termination_of_array(M+output_starts_at, &dataX[0][output_starts_at], FRONT_END, output_front_exponential_damping_factor);

	if (output_back_exponential_termination)
		gaussian_termination_of_array(M+output_starts_at, &dataX[0][output_starts_at], BACK_END, output_back_exponential_damping_factor);

		sprintf(print_buffer, "MEDD Output damped'd, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);



//  Normalize data so the entire range fits from +1 to -1
		normalize((int) M, &dataX[0][output_starts_at]);

		sprintf(print_buffer, "MEDD Output norm'd, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);

goto zzz;
//  Run it thru an exponential amplifier to emphasize the spikes
		for (j = 0; j < M+L; j ++)
		{
			dataX[0][j] = pow(dataX[0][j], 2);
		}

		normalize((int) M, &dataX[0][output_starts_at]);

		sprintf(print_buffer, "MEDD Output exp'ed, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);

zzz:
// The polarity of Reflections and absorptions are reversed in the 
// Output, so correct by inverting the waveform
		for (j = 0; j < M+L; j ++)
		{
//			dataX[0][j] =-dataX[0][j];
//
			dataX[0][j] = absolute_value(dataX[0][j]);
		}
		sprintf(print_buffer, "MEDD Output abs, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);



		for (j = 0; j < M+L; j ++)
		{
			dataX[0][j] = dataX[0][j] - 1.0 ;
		}
		sprintf(print_buffer, "MEDD Output -1'ed, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);


//	save_array_to_disk(hWnd, scan_number, M+L, &dataX[0][0], output_filename, binary_file);
	save_array_to_disk(hWnd, scan_number, M, &dataX[0][output_starts_at], output_filename, binary_file);



// which do you really want to do?
//???	if (machine == tower && plot_query >= plot_everything)
//        {
//		sprintf(print_buffer, "MEDD Output, Y(i) = f * x(i)");
//		plot_all(hWnd,-N,M+L,&dataX[0][0],print_buffer,
//			-0.5,0.5);
//        }

	  if (machine == tower && plot_query >= plot_only_end_results)
	    if( data_source == binary_simulation_1)
	    {

		sprintf(print_buffer, "Reflection Series");
		plot_all(hWnd,-N,M,&reflections[0],print_buffer,
			scale_min,scale_max);

            }
		sprintf(print_buffer, "MEDD Output, Y(i) = f * x(i)");
		if (machine == tower && plot_query >= plot_only_end_results)
		  plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);

	if (machine == tower && plot_query >= plot_only_end_results)
        {
		sprintf(print_buffer, "MEDD Output, Y(i) = f * x(i)");
		plot_all(hWnd,-N,M,&dataX[0][output_starts_at],print_buffer, scale_min,scale_max);
	}
		
		data_local_max = maximum_of_array( (int) (M * .9), &dataX[0][9]);
		data_local_min = minimum_of_array( (int) (M * .9), &dataX[0][9]);
//		data_local_min = maximum(absolute_value(data_max_x),absolute_value(data_min_x));



//		sprintf(print_buffer, "Blow up of MEDD Output, Y(i) = f * x(i)");
//		plot_all(hWnd,-N,M+L,&dataX[0][0],print_buffer, data_local_min, data_local_max);


   }

	sprintf(print_buffer,"MEDD has Completed %d scans", scan_number );
	stop_screen(hWnd, print_buffer);


}


void MEDD_deconvolution(HWND hWnd)
{
	#include "common.h"

	int ickx,icky;

	int i,j,k,s,u;
	double expect_Y;
	double D_norm, D_norm_final;
	double medd_filter, medd_output;
	double optimum_D_norm, optimum_output;
	int optimum_L = L_begin;

	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)
				{
		     
//if (j+k < 0 || j+k > M + L_Max)
//{
//sprintf(print_buffer, "in C loop: i %d,j %d,j+k %d\n",i,j,j+k);
//sprintf("C: dataX[i][j + k] = dataX[%d][%d] = %f\n",i,j+k,dataX[i][j+k]);
//stop_screen(hWnd, print_buffer);
//}

				    C[i][k] =
				    C[i][k] + dataX[i][j] * dataX[i][j+k];
				}
			}

		}

	}


	for (k = 0; k < L; k ++)
	{
		R[k] = 0;
		for (i = 0; i < N; i ++)
		{
   			R[k] = R[k] + C[i][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 (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]);
//if (j-s < 0 || j-s > M + L_Max)
//{
//sprintf(print_buffer,"s,i,j-s,k,W[s][k],dataX[i][j-s],F[k]:\n %d %d %d %d %f, %f, %f\n",
//s,i,j-s,k,W[s][k],dataX[i][j-s],F[k]);
//stop_screen(hWnd, print_buffer);
//}

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




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

/*
for (icky = 0; icky < M; icky ++)
printf("before s-k bussiness:	dataX[0][%d] = %#g\n",
icky,dataX[0][icky]);
*/


					/* 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];
					}

//sprintf(print_buffer, "Y[%d][%d] = %#g\n",k,s,Y[k][s]);
//stop_screen(hWnd, print_buffer);
				}
			}
						/* 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];

//sprintf(print_buffer,"\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 = (absolute_value(Y[i][j]))/(expect_Y);
			else
				D_norm = 0;


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

			}


		}
	}

/*
// I don't know where this came from...???
// but it was commented out for as long as I can remember
 
	D_norm_array[L] = 0;
	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
			D_norm_array[L] = D_norm_array[L] + Y[i][j] * Y[i][j];
	}
*/


	D_norm_array[L] = D_norm_final;



/* 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];
	}

/*  Don't really need...
	if (machine == tower && plot_query >= plot_everything)
	{	sprintf(print_buffer, "MEDD filter values, f");
		plot_all(hWnd,-1,M,&f,print_buffer,scale_min,scale_max);
        }

		sprintf(print_buffer, "MEDD Input, x(i)");
//		plot_all(hWnd,-N,M,&data[0][0],print_buffer, scale_min,scale_max);
*/

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		f[0][j] = Y[0][j];
//sprintf(print_buffer, "Y[0][j] = %#g\n", Y[0][j]);
//stop_screen(hWnd, print_buffer);

	}


	if (plot_query >=  plot_algorithm_tracking)
        {
	  sprintf(print_buffer, "RESULTS for this Filter:\nD_norm: %f\nMEDD Output: %f\nCurrent number of Filters: %d\nError Factor: %#g\n",
	  D_norm_final, medd_output,L,D_norm_array[L]);
	  stop_screen(hWnd, print_buffer);

			sprintf(print_buffer, "MEDD output for this filter, f");
			plot_all(hWnd,-1,M+L,&f,print_buffer,scale_min, scale_max);


	}


///  Use this if unknown machine causes portability problems.
///  Also, can save a file of filter (L) information here.
//  
//	if (machine != unknown && plot_query >= plot_everything)
//	{
//		sprintf(message_text,"  Continue?");
//  		stop_screen(hWnd, message_text);
//
//	}
// ---------------------


		 if (L == L_begin) goto label1;

		 if (D_norm_array[L] > D_norm_array[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];
//sprintf(print_buffer, "Y_optimum[i][j] = %#g\n", Y_optimum[i][j]);
//stop_screen(hWnd, print_buffer);
				}
			}
		}
	}



// End of L loop -----------


	i = 0;

 if (machine == unknown || save_error_file)
		id_err = fopen("medderr.dat", "wt");

	for (L = L_begin; L < L_max; L = L + L_inc)
	{
		if (D_norm_array[L] != 0)
		{
			f[0][i] = log10(D_norm_array[L]);
                }
		else
		{
			sprintf(print_buffer, "Warning!  D_norm is zero.  Chech you data.\n");
			stop_screen(hWnd, print_buffer);
		}

		if (machine == unknown || save_error_file)
			fprintf(id_err, "%e\t%e\n", (double)(L), D_norm_array[L]);

		i++;
	}

	if (machine == unknown || save_error_file)
		fclose(id_err);


	if (machine == tower && plot_query >= plot_algorithm_partial)
	{
		sprintf(print_buffer, "Log of Error vs. filter number");
		plot_all(hWnd,-1,i,&f,print_buffer, scale_min, scale_max);
	}

if (plot_query >= plot_only_end_results)
{
	sprintf(print_buffer, "FINAL RESULTS:\nD_norm: %f\nMEDD Output: %f\nBest Number of Filters: %d\nBest Error Factor: %#g\n",
	optimum_D_norm,optimum_output,optimum_L,D_norm_array[optimum_L]);
	stop_screen(hWnd, print_buffer);
}
///  Use this if unknown machine causes portability problems.
///  Also, can save a file of filter (L) information here.
//  
//	if (machine != unknown && plot_query >= plot_everything)
//	{
//		sprintf(message_text,"  Continue?");
//  		stop_screen(hWnd, message_text);
//	}
// ---------------------


			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];

//sprintf(print_buffer, "dataX = %#g\n", dataX[i][j]);
//stop_screen(hWnd, print_buffer);

				}
			}
	L = optimum_L;
sprintf(message_text," End of Process.  Continue?");
//stop_screen(hWnd, message_text);

	return;
}