#include <oaidl.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <comutil.h>
#include <ole2.h>
#include <cstdlib>

#define dllexport __declspec(dllexport)

extern "C"{

dllexport VARIANT _stdcall CheckAndRetVariant(VARIANT X)
	{
	VARIANT dvout;
	switch(X.vt)
		{
		case VT_DISPATCH:
			{
			EXCEPINFO excep;
			DISPPARAMS dispparams;
			unsigned int uiArgErr;
			DISPID dispidValue;     
			LPOLESTR XName = L"Value";
					
			X.pdispVal->GetIDsOfNames(IID_NULL, &XName,
					1, LOCALE_SYSTEM_DEFAULT, &dispidValue);
			// PropertyGet has no arguments
			dispparams.cArgs = 0;
			dispparams.cNamedArgs = 0;
			
			// Invoke PropertyGet and store values in vSource
			X.pdispVal->Invoke(dispidValue, IID_NULL, 
						LOCALE_SYSTEM_DEFAULT, DISPATCH_PROPERTYGET, 
						&dispparams, &dvout, &excep, &uiArgErr);
			X.pdispVal->Release();
			
			return dvout;
			}
			break;

		default:
			return X;
			break;

		}

	VariantClear(&dvout);
	VariantClear(&X);


	}

dllexport double _stdcall NormDist(double x, double Std)

{

double pi = 3.14159265358979;

double result = 1/(sqrt(2*pi)*Std)*exp(-(x*x)/(2*Std*Std));

return result;
	
}


dllexport VARIANT _stdcall KernelReg(VARIANT Y, VARIANT X, double h)
    
{

	X = CheckAndRetVariant(X);
	Y = CheckAndRetVariant(Y);
	VARIANT varvar;
		
    long ncols;
	long nrows;
	long ncolsY;
	long i, j, jx;
		
	ncols=(X.parray)->rgsabound[0].cElements;
	nrows=(X.parray)->rgsabound[1].cElements;
	ncolsY=(Y.parray)->rgsabound[0].cElements;

	double *ArrayX = new double [nrows*ncols];
	double *ArrayY = new double [nrows*ncolsY];
		
	for (i=0; i<nrows; ++i)
		{
		for (j=0; j<ncols; ++j)
			{
			long indices[] = {i+1,j+1};
			SafeArrayGetElement(X.parray, indices, &varvar);
			ArrayX[(i*ncols)+j] = varvar.dblVal;
			}
			for (jx=0; jx<ncolsY; jx++)
			{
			long indicesY[] = {i+1,jx+1};
			SafeArrayGetElement(Y.parray, indicesY, &varvar);
			ArrayY[(i*ncolsY)+jx] = varvar.dblVal;
			}
		}
	

	double prod;
	double *kernel = new double [nrows];
	double *kprod = new double [nrows];
	double *ye = new double [nrows];


		for(i=0; i<nrows; i++)
		{
		double w = 0;

	    for(j=0;j<nrows; j++)
			{
			prod =1;

			for (jx=0;jx<ncols;jx++)
				{
				prod =prod*NormDist(ArrayX[(i*ncols)+jx]
				- ArrayX[(j*ncols)+jx], h);
				}
			
			kernel[j]= prod;
			w = w + kernel[j];
			}
			double sumprod=0;
			for(j=0;j<nrows;j++)
				{
				kprod[j]=kernel[j]/w;
				
				sumprod= sumprod + kprod[j]*ArrayY[(j*ncolsY)];
				}
			
			
			ye[(i*ncolsY)]=sumprod;
		}


	VARIANT result;
	result.vt = VT_ARRAY | VT_VARIANT;

	SAFEARRAYBOUND sab[2];
	sab[0].lLbound = 1; sab[0].cElements = nrows;
	sab[1].lLbound = 1; sab[1].cElements = ncolsY;

	result.parray = SafeArrayCreate(VT_VARIANT, 2, sab);
	
	VARIANT tmp;
	tmp.vt = VT_R8;

	for(i=0; i<nrows; i++)
		{
		for(j=0;j<ncolsY;j++)
			{
			tmp.dblVal= ye[(i*ncolsY)+j];
			long indi [] = {i+1,j+1};
			SafeArrayPutElement(result.parray, indi, &tmp);
			}
		}
	
	return result;
 
	/* Deallocate the arrays */
	delete [] ArrayX;
	delete [] ArrayY;
	delete [] kernel;
	delete [] kprod;
	delete [] ye;

	VariantClear(&X);
	VariantClear(&Y);
	VariantClear(&varvar);
	VariantClear(&result);
	
}

dllexport VARIANT _stdcall KernelRegVali(VARIANT Yat, VARIANT Xat, 
							  double h, VARIANT X)
{
	long n,v,Ndat;
	long i, j, jx;
	VARIANT varvar;
	
	X =CheckAndRetVariant(X);
	Yat = CheckAndRetVariant(Yat);
	Xat = CheckAndRetVariant(Xat);

	v=(X.parray)->rgsabound[0].cElements;
	Ndat=(X.parray)->rgsabound[1].cElements;
	n=(Yat.parray)->rgsabound[1].cElements;
	
	double *ArrayX = new double [Ndat*v];
	double *dXat = new double [n*v];
	double *dYat = new double [n];
		
	for(i=0; i<n; i++)
		{
		for(j=0; j<v; j++)
			{
			long indices[] = {i+1, j+1};
			SafeArrayGetElement(Xat.parray, indices, &varvar);
			dXat[(i*v)+j] =varvar.dblVal;
			}
			long indicesY [] = {i+1,1};
			SafeArrayGetElement(Yat.parray, indicesY, &varvar);
			dYat[i] =varvar.dblVal;
		}
	
	for (i=0; i<Ndat; i++)
		{
		for (j=0; j<v; j++)
			{
			long indices[] = {i+1,j+1};
			SafeArrayGetElement(X.parray, indices, &varvar);
			ArrayX[(i*v)+j] = varvar.dblVal;
			}
		}

	double prod;
	double *kernel = new double [n];
	double *kprod = new double [n];
	double *ye = new double [Ndat];

		for(i=0; i<Ndat; i++)
		{
		double w = 0;

	    for(j=0;j<n; j++)
			{
			prod =1;

			for (jx=0;jx<v;jx++)
				{
				prod =prod*NormDist(ArrayX[(i*v)+jx]
				- dXat[(j*v)+jx], h);
				}
			
			kernel[j]= prod;
			w = w + kernel[j];
			}
			double sumprod=0;
			for(j=0;j<n;j++)
				{
				kprod[j]=kernel[j]/w;
				
				sumprod= sumprod + kprod[j]*dYat[j];
				}

			ye[i]=sumprod;
		}


	VARIANT result;
	result.vt = VT_ARRAY | VT_VARIANT;

	SAFEARRAYBOUND sab[2];
	sab[0].lLbound = 1; sab[0].cElements = Ndat;
	sab[1].lLbound = 1; sab[1].cElements = 1;

	result.parray = SafeArrayCreate(VT_VARIANT, 2, sab);
	
	VARIANT tmp;
	tmp.vt = VT_R8;

	for(i=0; i<Ndat; i++)
		{
		for(j=0;j<1;j++)
			{
			tmp.dblVal= ye[i+j];
			long indi [] = {i+1,j+1};
			SafeArrayPutElement(result.parray, indi, &tmp);
			}
		}
	
	return result;
 
	/* Deallocate the arrays */
	delete [] ArrayX;
	delete [] kernel;
	delete [] kprod;
	delete [] ye;
	delete [] dXat;
	delete [] dYat;
	
	VariantClear(&X);
	VariantClear(&varvar);
	VariantClear(&Yat);
	VariantClear(&Xat);
	VariantClear(&result);
	
}

dllexport VARIANT _stdcall KernelReg_Cross(VARIANT Y, VARIANT X, double h)
    
{
	X = CheckAndRetVariant(X);
	Y = CheckAndRetVariant(Y);
	VARIANT varvar;
		
    long ncols;
	long nrows;
	long ncolsY;
	long i, j, jx;
		
	ncols=(X.parray)->rgsabound[0].cElements;
	nrows=(X.parray)->rgsabound[1].cElements;
	ncolsY=(Y.parray)->rgsabound[0].cElements;

	double *ArrayX = new double [nrows*ncols];
	double *ArrayY = new double [nrows*ncolsY];

		
	for (i=0; i<nrows; i++)
		{
		for (j=0; j<ncols; j++)
			{
			long indices[] = {i+1,j+1};
			SafeArrayGetElement(X.parray, indices, &varvar);
			ArrayX[(i*ncols)+j] = varvar.dblVal;
			}
			for (jx=0; jx<ncolsY; jx++)
			{
			long indicesY[] = {i+1,jx+1};
			SafeArrayGetElement(Y.parray, indicesY, &varvar);
			ArrayY[(i*ncolsY)+jx] = varvar.dblVal;
			}
		}
	

	double *CrossX =new double [(nrows-1)*ncols];
	double *CrossY =new double [(nrows-1)*ncolsY];
	double prod;
	double *kernel = new double [nrows-1];
	double *kprod = new double [nrows-1];
	double *ye = new double [nrows];


		for(i=0; i<nrows; i++)
		{
		
		long b =-1;
		
	    for(j=0;j<nrows-1; j++)
			{
			if ( j != i )
				b=b+1;
			else 
				{
				b=b+2;
				}
	
			for (jx=0;jx<ncols;jx++)
				{
				CrossX[(j*ncols)+jx] =ArrayX[(b*ncols)+jx];
				}
			CrossY[(j*ncolsY)]=ArrayY[(b*ncolsY)];
			}

		double w = 0;

	    for(j=0;j<nrows-1; j++)
			{
			prod =1;

			for (jx=0;jx<ncols;jx++)
				{
				prod =prod*NormDist(ArrayX[(i*ncols)+jx]
				- CrossX[(j*ncols)+jx], h);
				}
			
			kernel[j]= prod;
			w = w + kernel[j];
			}
			double sumprod=0;
			for(j=0;j<nrows-1;j++)
				{
				kprod[j]=kernel[j]/w;
				
				sumprod= sumprod + kprod[j]*CrossY[(j*ncolsY)];
				}
			
			
			ye[(i*ncolsY)]=sumprod;
		}


	VARIANT result;
	result.vt = VT_ARRAY | VT_VARIANT;

	SAFEARRAYBOUND sab[2];
	sab[0].lLbound = 1; sab[0].cElements = nrows;
	sab[1].lLbound = 1; sab[1].cElements = ncolsY;

	result.parray = SafeArrayCreate(VT_VARIANT, 2, sab);
	
	VARIANT tmp;
	tmp.vt = VT_R8;

	for(i=0; i<nrows; i++)
		{
		for(j=0;j<ncolsY;j++)
			{
			tmp.dblVal= ye[(i*ncolsY)+j];
			long indi [] = {i+1,j+1};
			SafeArrayPutElement(result.parray, indi, &tmp);
			}
		}
	
	return result; 
	// Deallocate the arrays

	delete [] ArrayX;
	delete [] ArrayY;
	delete [] kernel;
	delete [] kprod;
	delete [] ye;
	delete [] CrossX;
	delete [] CrossY;

	VariantClear(&X);
	VariantClear(&Y);
	VariantClear(&varvar);
	VariantClear(&result);
	
}

}