#include "mex.h"
#include "matlab.h"
#include "math.h"
#define LOG(a) ((a)>1e-20 ? log((a)):log(1e-20))
#include <afxwin.h>         // MFC core and standard components
#include <afxext.h>         // MFC extensions
#include <afxdisp.h>        // MFC Automation classes
#include <afxdtctl.h>		// MFC support for Internet Explorer 4 Common Controls
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	if((nrhs!=3)||(nlhs!=1))
	{
		mexPrintf("********************--------correl_dim----------*******************\n");
		mexPrintf("本函數計算相關積供后續程序計算關聯維\n");
		mexPrintf("c=correl_dim(data,n_eps,stepsize)\n");
		mexPrintf("輸入：data-------待處理數據 相點坐標維數*點數；\n");
		mexPrintf("      n_eps---搜尋徑點數 1*1; \n");
		mexPrintf("      stepsize---搜尋徑步長 1*1;   eps[l]=maxdist*2.0/pow(buchang,n_eps-l);\n");
		mexPrintf("輸出：c----------相關積結果 搜尋徑點數*1;\n");
		mexPrintf("********************--------correl_dim----------*******************\n");
		mexErrMsgTxt("你想怎么用亞");
	}
    
	int Mr0,Nr0;
	Mr0=mxGetM(prhs[0]);
	Nr0=mxGetN(prhs[0]);
	if(Nr0<102)
		mexErrMsgTxt("數據點數過少，不宜計算。");//延遲100個數據點，對于特征數據這個是不必要的。
	
	int n_eps=mxGetScalar(prhs[1]);	
	if(n_eps<1)
		mexErrMsgTxt("n_eps<1");
	
	double buchang=mxGetScalar(prhs[2]);
	
	double *y,*jieguo;
	plhs[0]=mxCreateDoubleMatrix(n_eps,1,mxREAL);
	y     =mxGetPr(prhs[0]);
	jieguo=mxGetPr(plhs[0]);	
	double *eps=(double*)calloc(n_eps,sizeof(double));	
	
	for(int l=0;l<n_eps;l++)
	{
		eps[l]=0;
		jieguo[l]=0;
	}

	double maxdist=0;
	for(int j=0;j<Nr0;j++)
	{
		double dist=0;
		double *a=y;
		double *b=y+j*Mr0;
		for(int mm=0;mm<Mr0;mm++)
			dist=dist+(a[mm]-b[mm])*(a[mm]-b[mm]);
		if(maxdist<dist)
			maxdist=dist;
	}
	if(maxdist<=0)
		mexErrMsgTxt("數都一樣");
	for(l=0;l<n_eps;l++)
		eps[l]=maxdist*2.0/pow(buchang,n_eps-l);
	for(j=0;j<Nr0-100;j++)
	{
		for(int k=j+100;k<Nr0;k++)
		{
			double dist=0;
			double *a=y+j*Mr0;
			double *b=y+k*Mr0;
			for(int mm=0;mm<Mr0;mm++)
				dist=dist+(a[mm]-b[mm])*(a[mm]-b[mm]);
			for(int l=0;l<n_eps;l++)
			{
				if(dist<=eps[l])
					jieguo[l]+=1;	
			}
		}
	}
	for(l=0;l<n_eps;l++)
	{
		jieguo[l]=LOG(jieguo[l])/(0.5*LOG(buchang));
	}
	free(eps);
}