xldlas/xldlas/filters.c

/*      -------------------------------------------------------------------
	xldlas -- A Program for Statistics

	Copyright (C) 1996, 1997 Thor Sigvaldason

	This file contains all the filtering routines

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

#include "xldlas.h"

extern void say_status(char the_status[XLDLASMAX_INPUT]);
extern void inhibit_input();
extern void reenable_input();
extern void simple_line_output(char which_routine[XLDLASMAX_INPUT], char the_output[XLDLASMAX_INPUT]);


void do_detrend_filter(int f_begin, int f_end, int which_one)
{
	int i, j, k, length, max;
	float intercept, slope, running, temp1, temp2, coefs[3], mean;
	char a_string[XLDLASMAX_INPUT];
	double *bigmatrix;
	double *matrix;
	worksize = 0;
	mean = 0.0;
	if(f_end > data_matrix[which_one].obs) f_end = data_matrix[which_one].obs;
	for(i = f_begin - 1; i < f_end; i++)
	{
		if(*(fvector[which_one] + i) != missing_value)
		{
			working[worksize] = *(fvector[which_one] + i);
			mean = mean + working[worksize];
			worksize++;
		}
	}
	if(worksize < 2)
	{
		sprintf(a_string,"%s does not have enough observations", data_matrix[which_one].name);
		fl_show_alert(a_string, "","",TRUE);
		return;
	}
	mean = mean / worksize;
	length = worksize;
	bigmatrix = (double *) malloc ((length * 3) * sizeof (double));
	if(!bigmatrix)
	{
		fl_show_alert("Not Enough Memory to Build Detrend Matrix","","",TRUE);
		return;
	}
	for(i = 0; i < length; i++)
	{
		*(bigmatrix + i) = 1;
		*(bigmatrix + i + length) = i+1;
		*(bigmatrix + i + 2 * length) = working[i];
	}
	matrix = (double *) malloc (6 * sizeof (double));
	if(!matrix)
	{
		fl_show_alert("Not Enough Memory to Build Detrend Sub-Matrix","","",TRUE);
		free(bigmatrix);
		return;
	}
	for(i = 0; i < 2; i++)
	{
		for(j = 0; j <= 2; j++)
		{
			running = 0.0;
			for(k = 0; k < length; k++)
			{
				running = running + *(bigmatrix + k + (i * length)) * *(bigmatrix + k + (j * length));
			}
			*(matrix + i + (j * 2)) = running;
		}
	}
	for(i = 0; i < 2; i++)
	{
		max = i;
		for(j = i + 1; j < 2; j++)
		{
			temp1 = *(matrix + j + (i * 2));
			temp2 = *(matrix + max + (i * 2));
			if(fabs(temp1) > fabs(temp2)) max = j;
			for(k = 0; k <= 2; k++)
			{
				running = *(matrix + i + (k * 2));
				*(matrix + i + (k * 2)) = *(matrix + max + (k * 2));
				*(matrix + max + (k * 2)) = running;
			}
		}
		for(j = i + 1; j < 2; j++)
		{
			for(k = 2; k >= 0; k--)
			{
				temp1 = *(matrix + i + (i * 2));
				if(temp1 == 0)
				{
					sprintf(a_string,"%s results in a singular matrix",data_matrix[which_one].name);
					fl_show_alert(a_string,"Sorry, there's nothing I can do","",TRUE);
					free(bigmatrix);
					free(matrix);
					return;
				}
				*(matrix + j + (k * 2)) = *(matrix + j + (k * 2)) - *(matrix + i + (k * 2)) * *(matrix + j + (i * 2)) / *(matrix + i + (i * 2));
			}
		}
	}
	for(i = 0; i < 3; i++)
	{
		coefs[i] = 0.0;
	}
	for(i = 1; i >= 0; i--)
	{
		running = 0.0;
		for(j = i + 1; j <= 2; j++)
		{
			running = running + *(matrix + i + (j * 2)) * coefs[j];
			temp1 = *(matrix + i + (i * 2));
			if(temp1 == 0)
			{
				sprintf(a_string,"%s results in a singular matrix",data_matrix[which_one].name);
				fl_show_alert(a_string,"Sorry, there's nothing I can do","",TRUE);
				free(bigmatrix);
				free(matrix);
				return;
			}
			coefs[i]= ( *(matrix + i + (2 * 2)) - running ) / *(matrix + i + (i * 2));
		}
	}
	intercept = coefs[0];
	slope = coefs[1];
	for(i = f_begin - 1; i < f_end; i++)
	{
		if(*(fvector[which_one] + i) != missing_value)
		{
			*(fvector[which_one] + i) = mean + *(fvector[which_one] + i) - (intercept + ((i - f_begin) * slope));
		}
	}
	oktoquit = FALSE;
	sprintf(a_string,"%s detrended (intercept = %f, slope %f)", data_matrix[which_one].name, intercept, slope);
	simple_line_output("filter",a_string);
	free(bigmatrix);
	free(matrix);
}

void do_outliers_filter(int f_begin, int f_end, int which_one)
{
	int i, numb_removed, nrealobs;
	float mean, variance, scratch;
	char a_string[XLDLASMAX_INPUT];
	worksize = 0;
	mean = 0.0;
	nrealobs = 0;
	numb_removed = 0;
	if(f_end > data_matrix[which_one].obs) f_end = data_matrix[which_one].obs;
	for(i = f_begin - 1; i < f_end; i++)
	{
		working[worksize] = *(fvector[which_one] + i);
		if(working[worksize] != missing_value)
		{
			mean = mean + working[worksize];
			nrealobs++;
		}
		worksize++;
	}
	if(nrealobs < 2)
	{
		sprintf(a_string,"%s does not have enough observations", data_matrix[which_one].name);
		fl_show_alert(a_string, "","",TRUE);
		return;
	}
	mean = mean / nrealobs;
	scratch = 0.0;
	variance = 0.0;
	for(i = 0; i < worksize; i++)
	{
		if(working[i] != missing_value)
		{
			scratch = scratch + (working[i] - mean);
			variance = variance + pow((working[i] - mean),2);
		}
	}
	variance = (variance - scratch * scratch / nrealobs) / nrealobs;
	variance = pow(variance,0.5);
	for(i = 0; i < worksize; i++)
	{
		if(working[i] != missing_value)
		{
			if(working[i] < mean - 3.0 * variance || working[i] > mean + 3.0 * variance)
			{
				numb_removed++;
				working[i] = missing_value;
			}
		}
	}
	if(numb_removed > 0)
	{
		for(i = 0; i < worksize; i++)
		{
			*(fvector[which_one] + (f_begin - 1) + i) = working[i];
		}
		oktoquit = FALSE;
	}
	sprintf(a_string,"Outliers removed from %s (%d value(s) changed)", data_matrix[which_one].name, numb_removed);
	simple_line_output("filter",a_string);
}

void do_zo_filter(int f_begin, int f_end, int which_one)
{
	int i;
	float smallest, largest;
	char a_string[XLDLASMAX_INPUT];
	smallest = missing_value;
	largest = missing_value;
	if(f_end > data_matrix[which_one].obs) f_end = data_matrix[which_one].obs;
	for(i = f_begin - 1; i < f_end; i++)
	{
		if(*(fvector[which_one] + (f_begin -1) + i) != missing_value)
		{
			if(smallest == missing_value || smallest > *(fvector[which_one] + (f_begin - 1) + i))
			{
				smallest = *(fvector[which_one] + (f_begin - 1) + i);
			}
			if(largest == missing_value || largest < *(fvector[which_one] + (f_begin - 1) + i))
			{
				largest = *(fvector[which_one] + (f_begin - 1) + i);
			}
		}
	}
	if(smallest == missing_value)
	{
		fl_show_alert("No observations in selected range","","",TRUE);
		return;
	}
	if(smallest == largest)
	{
		fl_show_alert("No variation in selected range","(I'm leaving the data unchanged)","",TRUE);
		return;
	}
	for(i = f_begin -1; i < f_end; i++)
	{
		if(*(fvector[which_one] + (f_begin -1) + i) != missing_value)
		{
			*(fvector[which_one] + (f_begin - 1) + i) = ( *(fvector[which_one] + (f_begin - 1) + i) - smallest ) / ( largest - smallest);
		}
	}
	sprintf(a_string,"Values of %s (obs %i to %i) transformed to lie between 0.0 and 1.0", data_matrix[which_one].name, f_begin, f_end);
	simple_line_output("filter",a_string);
	oktoquit=FALSE;
}

void do_filter(FL_OBJECT *obj, long arg)
{
	int numb_filters, to_filter[MAX_VARS], i;
	numb_filters = 0;
	for(i = 0; i < numb_variables; i++)
	{
		if(fl_isselected_browser_line(filter_browser, i+1))
		{
			to_filter[numb_filters] = i;
			numb_filters++;
		}
	}
	if(numb_filters == 0)
	{
		fl_show_alert("No Variable(s) Selected","","",TRUE);
		return;
	}
	all_start = fl_get_counter_value(filter_from_counter);
	all_stop = fl_get_counter_value(filter_to_counter);
	if(all_stop < all_start)
	{
		fl_show_alert("From Value greater than To Value","","",TRUE);
		return;
	}
	say_status("Filtering...");
	if(fl_get_button(filter_outliers_button) == TRUE)
	{
		for(i = 0; i < numb_filters; i++)
		{
			do_outliers_filter(all_start, all_stop, to_filter[i]);
		}
	}
	if(fl_get_button(filter_detrend_button) == TRUE)
	{
		for(i = 0; i < numb_filters; i++)
		{
			do_detrend_filter(all_start, all_stop, to_filter[i]);
		}
	}
	if(fl_get_button(filter_zo_button) == TRUE)
	{
		for(i = 0; i < numb_filters; i++)
		{
			do_zo_filter(all_start, all_stop, to_filter[i]);
		}
	}
	say_status("Waiting for filter specs...");
}

void start_filters(int filt_numb)
{
	int i, largest;
	inhibit_input();
	largest = 0;
	fl_clear_browser(filter_browser);
	for(i = 0; i < numb_variables; i++)
	{
		if(largest < data_matrix[i].obs) largest = data_matrix[i].obs;
		fl_addto_browser(filter_browser,data_matrix[i].name);
	}
	fl_set_counter_value(filter_from_counter, 1);
	fl_set_counter_bounds(filter_from_counter, 1, largest);
	fl_set_counter_value(filter_to_counter, largest);
	fl_set_counter_bounds(filter_to_counter, 1, largest);
	say_status("Waiting for filtering specs...");
	if(window_geometry[XLDLAS_FILTER][0] != -1)
	{
		fl_set_form_geometry(filter_window,
					window_geometry[XLDLAS_FILTER][0],
					window_geometry[XLDLAS_FILTER][1],
					window_geometry[XLDLAS_FILTER][2],
					window_geometry[XLDLAS_FILTER][3]);
	}
	fl_show_form(filter_window,FL_PLACE_FREE,FL_FULLBORDER,"Do Some Filtering");
	fl_set_button(filter_outliers_button, 0);
	fl_set_button(filter_detrend_button, 0);
	fl_set_button(filter_zo_button, 0);
	if(filt_numb == 1) fl_set_button(filter_outliers_button, 1);
	if(filt_numb == 2) fl_set_button(filter_detrend_button, 1);
	if(filt_numb == 3) fl_set_button(filter_zo_button, 1);
}

void done_filter(FL_OBJECT *obj, long arg)
{
	window_geometry[XLDLAS_FILTER][0] = obj->form->x;
	window_geometry[XLDLAS_FILTER][1] = obj->form->y;
	window_geometry[XLDLAS_FILTER][2] = obj->form->w;
	window_geometry[XLDLAS_FILTER][3] = obj->form->h;
	fl_hide_form(filter_window);
	reenable_input();
	say_status("Ready");
}

void filter_routines(FL_OBJECT *menu, long user_data)
{
	int choice;
	choice = fl_get_menu(menu);
	start_filters(choice);
}