lyap_exp_cambridge.cpp

本函数计算距离增长结果供拟合lyapunov指数 详情见说明

#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
//         mex lyap_exp_cambridge.cpp

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	if((nrhs!=5)||(nlhs!=1))
	{
		mexPrintf("********************--------lyap_exp----------*******************\n");
		mexPrintf("本函数计算距离增长结果供拟合lyapunov指数\n");
		mexPrintf("c=lyap_exp_cambridge(data,n_eps,stepsize,yanchi,t)\n");
		mexPrintf("输入：0---data-------待处理数据 相点坐标维数*点数；\n");
		mexPrintf("      1---n_eps------搜寻半径数 1*1; \n");
		mexPrintf("      2---stepsize---搜寻半径变化步长   1*1;\n");
		mexPrintf("      3---yanchi-----延迟点数   1*1;\n");
		mexPrintf("      4---t----------采样时间   1*1;\n");
		mexPrintf("输出：0---c----------结果       yanchi*n_eps;\n");
		mexPrintf("基于cambridge书63页公式5.2\n");
		mexPrintf("********************--------lyap_exp----------*******************\n");
		mexErrMsgTxt("你想怎么用亚");
	}
	    
	int Mr0,Nr0;
	Mr0=mxGetM(prhs[0]);
	Nr0=mxGetN(prhs[0]);
	int yanchi=mxGetScalar(prhs[3]);
	if(Nr0<100+yanchi)
		mexErrMsgTxt("数据点数过少，不宜计算。");//延迟100个数据点，对于特征数据这个是不必要的。
	int n_eps=mxGetScalar(prhs[1]);
	if(n_eps<1)
		mexErrMsgTxt("n_eps<1");
	double stepsize=mxGetScalar(prhs[2]);
	double t=mxGetScalar(prhs[4]);
	
	double *y,*jieguo;
	plhs[0]=mxCreateDoubleMatrix(yanchi,n_eps,mxREAL);
	y     =mxGetPr(prhs[0]);
	jieguo=mxGetPr(plhs[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("数都一样");
	
	int tiaochulai=0;
	double eps=maxdist/100;
	int *geshu=(int *)calloc(n_eps,sizeof(int));
	double *temps=(double *)calloc(yanchi*n_eps,sizeof(double));
	double *epsxulie=(double *)calloc(n_eps,sizeof(double));
	double *tempyanchi=(double *)calloc(yanchi,sizeof(double));

	int zuodaodijigedianle=0;
	do{
		for(int l=0;l<n_eps*yanchi;l++)
			jieguo[l]=0;	
		for(j=0;j<n_eps;j++)
			epsxulie[j]=eps/pow(stepsize,j);
		for(j=0;j<Nr0-yanchi;j++)
		{
			for(int k=0;k<n_eps;k++)
				geshu[k]=0;
			for(k=0;k<n_eps*yanchi;k++)
				temps[k]=0;

			for(k=0;k<Nr0-yanchi;k++)
			{
				if (abs(k-j)<10) 
					continue;
				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]);
				if(dist<epsxulie[0])
				{
					tempyanchi[0]=dist;
					for(int yan=1;yan<yanchi;yan++)
					{
						tempyanchi[yan]=0;
						double *a_yan=a+yan*Mr0;
						double *b_yan=b+yan*Mr0;
						for(int mm=0;mm<Mr0;mm++)
							tempyanchi[yan]=tempyanchi[yan]+(a_yan[mm]-b_yan[mm])*(a_yan[mm]-b_yan[mm]);
					}
					for(int nneps=0;nneps<n_eps;nneps++)
					{					
						if(epsxulie[nneps]>dist)
						{
							geshu[nneps]++;
							for(int yan1=0;yan1<yanchi;yan1++)
							{
								temps[nneps*yanchi+yan1]+=tempyanchi[yan1];
							}						
						}
					}				
				}
			}
			
			for(k=0;k<n_eps;k++)
			{
				if(geshu[k]>10)
				{
					for(int jj=0;jj<yanchi;jj++)					
						jieguo[k*yanchi+jj]+=LOG(temps[k*yanchi+jj]/geshu[k]);
				}
				if(geshu[k]>Nr0/10)
				{
					mexPrintf("有%d个近邻点。\n",geshu[k]);
					tiaochulai=1;
					break;
				}
			}
			if(tiaochulai)
			{
				tiaochulai=0;
				eps=eps/2.0;
				break;
			}
		}
		if (j==Nr0-yanchi) {
			zuodaodijigedianle=1;
		}
	} while(zuodaodijigedianle!=1);

	for(int qq=0;qq<n_eps*yanchi;qq++)
	{
		jieguo[qq]/=(2.0*t*(Nr0-yanchi));
	}
	free(geshu);
	free(temps);
	free(tempyanchi);
	free(epsxulie);
}