/*

bandwidth2.c
Kernel bandwidth selection for 1D data.

This is **sample code** for the technical report:

Alexandre, L.A, "A Solve-the-Equation Approach for Unidimensional Data Kernel Bandwidth Selection",
Technical Report ISBN 978-989-654-005-0, DI-UBI, Covilhã, Portugal, 2008.
Link: http://www.di.ubi.pt/~lfbaa/pubs/tecrep2008.pdf

Please cite this report if you use this code.

THERE IS NO GUARANTY OF ANY KIND WITH THIS SOFTWARE.
YOU RUN IT AT YOUR OWN RISK.

Author: Luis A. Alexandre
Date: 2007-12-27

*/

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

double dados[1000]; // lazy!

double k4f(double h, int n)
{
	int i,j;
	double out=0, aux, h2=h*h;
	for(i=0; i<n; i++)
		for(j=0; j<i-1; j++)
		{
			aux=pow(dados[i]-dados[j],2);
			out+=(pow(aux-6*h2,2)-24*h2*h2)*exp(-(aux/(4*h2)));
		}

	out/=(2.*pow(h,3));
	out+=3*n*h;
	return out;
}

// algoritmo recursivo - eq. 11
double hopt(double n)
{
	double S, e=0.0001, h0=1 ,h1=h0+e+0.1, k4;
	int i=0;

	while((fabs(h1-h0)>e) && (i<1000))
	{
		h0=h1;
		k4=k4f(h0,n);
		h1=pow(4*n*pow(h0,6)/k4,0.2);
		h1=(h0+h1)/2;  
		printf("%d h0=%.5lf h1=%.5lf dif=%.5f\n",i,h0,h1,fabs(h1-h0));
		i++;
	}
	return h1;
}

int main()
{
	int i=0;
	FILE *fp;
	double n;

	fp=fopen("buffaloSnowfall.dat","rt");
//	fp=fopen("OF107.dat","rt");
//	fp=fopen("swissForgedBankNotesBottom.dat","rt");
	if(fp==NULL)
	{
		printf("Error opening the file\n");
		exit(1);
	}

	// load data
	while(fscanf(fp,"%lf",&n)!=EOF && i<1000)
	{
		dados[i]=n;
		i++;
	}

	if(i==1000)
	{
		printf("Data has more than 1000 samples! Change the size of the data structure.\n");
		exit(1);
	}

	printf("n=%d hopt=%.4f\n",i,hopt(i));

	return 0;
}