#include <math.h>
#include "mex.h"
#include <stdio.h>
#include <stdlib.h>
#include <matrix.h>

double ** TwoArrayAlloc(int r,int c)
{  //兩維內存動態分配函數
	double  *x, **y;
	int n;
	
	x=(double *)calloc(r*c,sizeof(double));
	if(x==NULL)
    {
		printf("\nFAilue in memory applying.");
		//AfxMessageBox("Failure in memory applying.");  
		exit(0);
    }
	y=(double **)calloc(r,sizeof(double *));
	for(n=0;n<=r-1;++n)
		y[n]=&x[c*n];
	return(y);
}

void TwoArrayFree(double **x)//2 維內存釋放函數
{
	free(x[0]);
	free(x);
}

double sum(double *array,int n)
{
	int i;
	double result=0.0;
	for (i=0;i<n;i++)
		result=result+array[i];
	return result;

}


void mexFunction (int nlhs, mxArray *plhs[],			// 輸出參數個數，及輸出參數數組
				  int nrhs, const mxArray *prhs[])	// 輸入參數個數，及輸入參數數組
{
	double *pdata,*pE1;
	int m1,tau1,N1,P1;
	
	if (nrhs!=4) mexErrMsgTxt("需要4個參數!");  //檢查輸入參數的個數
	
    // 取得輸入參數
    pdata = mxGetPr(prhs[0]);      // 時間序列（列向量）      
    tau1 = (int) *mxGetPr(prhs[1]);            // 最小嵌入維數 
	m1 = (int) *mxGetPr(prhs[2]);            // 最大嵌入維數 
    P1 = (int) *mxGetPr(prhs[3]);      // 時間延遲   
	N1 = mxGetM(prhs[0]);              // 序列長度
	
	
    // 為輸出變量分配內存空間
	plhs[0]=mxCreateDoubleMatrix(N1,1,mxREAL);
	
	
	
	// 取得輸出參數指針
	pE1 = mxGetPr(plhs[0]);
    
	
    // 調用 C 運算函數 (該函數名不能和本文件名重名)
	LYA_FUNCTION(pdata,tau1,m1,N1,P1,pE1);    
	

	
 

//  求最大、最小和平均相點距離
	int i,j,k,ii,tau=6,m=8,P=50,M,N,Loc_DK,old_Loc_DK;
	double **Y,*lmd;
   double max_d = 0,d;                                       //  %最大相點距離
    double min_d = 1.0e+100;                                 // %最小相點距離
    double avg_dd = 0;
   double dlt_eps,min_eps,max_eps,DK;
   double sum_lmd,DK1,old_DK,avg_d;
   double point_num,cos_sita,zjfwcs,dnew,DOT,CTH,*lambda_wolf;
   int min_point=1;
   int MAX_CISHU=5;
   	
   N=length;
    M=N-(m-1)*tau; 
   lambda_wolf=(double*)malloc(sizeof(double));
   lmd=(double*)malloc((M-2)*sizeof(double));
    //Y=reconstitution(data,N,m,tau);%相空間重構
                                     //%重構相空間中相點的個數
    Y=TwoArrayAlloc(m,M);
	for (j=0;j<M;j++)  //相空間重構
	{	for (i=0;i<m;i++) 
	        Y[i][j]=*(pdata+i*tau+j);
	}
	

    for (i=1;i<=(M-1);i++)   //i = 1 : (M-1)
	{   for (j=i+1;j<=M;j++)    //j = i+1 : M
		{   d = 0;
            for (k=1;k<=m;k++)     //k = 1 : m
			{   d = d + (Y[k-1][i-1]-Y[k-1][j-1])*(Y[k-1][i-1]-Y[k-1][j-1]);   }
            d = sqrt(d);
            if (max_d < d)
			{  max_d = d;}
            if (min_d > d)
			{  min_d = d;}
            avg_dd = avg_dd + d;
        }
    }
    avg_d = 2*avg_dd/(M*(M-1));               // %平均相點距離
    dlt_eps = (avg_d - min_d) * 0.02 ;       //  %若在min_eps～max_eps中找不到演化相點時，對max_eps的放寬幅度
    min_eps = min_d + dlt_eps / 2 ;          //  %演化相點與當前相點距離的最小限
    max_eps = min_d + 2 * dlt_eps  ;         //  %&&演化相點與當前相點距離的最大限
    
//     從P+1～M-1個相點中找與第一個相點最近的相點位置(Loc_DK)及其最短距離DK

    DK = 1.0e+100;                            // %第i個相點到其最近距離點的距離
    Loc_DK = 2;                              //  %第i個相點對應的最近距離點的下標
    for  (i=(P+1);i<=(M-1);i++)         // i = (P+1):(M-1)   %限制短暫分離，從點P+1開始搜索
	{   d = 0;
        for (k=1;k<=m;k++)       //k = 1 : m
		{   d = d + (Y[k-1][i-1]-Y[k-1][0])*(Y[k-1][i-1]-Y[k-1][0]);}
        
        d = sqrt(d);
        if ((d<DK)&&(d>min_eps)) 
		{  DK = d;
           Loc_DK = i;                  //Loc_DK = i;
        }
    }
//    以下計算各相點對應的李氏數保存到lmd()數組中
//     i 為相點序號，從1到(M-1)，也是i-1點的演化點；Loc_DK為相點i-1對應最短距離的相點位置，DK為其對應的最短距離
//     Loc_DK+1為Loc_DK的演化點，DK1為i點到Loc_DK+1點的距離，稱為演化距離
//     前i個log2（DK1/DK）的累計和用于求i點的lambda值
//	double sun_lmd,DK1,old_Loc_DK,old_DK

    sum_lmd = 0 ;                              //% 存放前i個log2（DK1/DK）的累計和
    for (i=2;i<=(M-1);i++)      //i = 2 : (M-1)           // % 計算演化距離      
	{   DK1 = 0.0;
        for  (k=1;k<=m;k++) //k = 1 : m
		{  // DK1 = DK1 + (Y[k][i]-Y[k][Loc_DK+1])*(Y[k][i]-Y[k][Loc_DK+1]);
		    DK1= DK1+(Y[k-1][i-1]-Y[k-1][Loc_DK])*(Y[k-1][i-1]-Y[k-1][Loc_DK]);
			 //DK1= DK1+(*(pdata+(k-1)*tau+i-1)-*(pdata+(k-1)*tau+i))*(*(pdata+(k-1)*tau+i-1)-*(pdata+(k-1)*tau+i));
        }
        DK1 = sqrt(DK1);
        old_Loc_DK = Loc_DK ;                  //% 保存原最近位置相點
        old_DK=DK;

//     計算前i個log2（DK1/DK）的累計和以及保存i點的李氏指數
        if ((DK1 != 0)&&( DK != 0))
        {  sum_lmd = sum_lmd + log(DK1/DK) /log(2);}
        lmd[i-2] = sum_lmd/(i-1);
//     以下尋找i點的最短距離：要求距離在指定距離范圍內盡量短，與DK1的角度最小
	//	double point_num,cos_sita,zjfwcs,dnew,DOT,CTH,cos_sita;
        point_num = 0  ; //% &&在指定距離范圍內找到的候選相點的個數
        cos_sita = 0  ; //%&&夾角余弦的比較初值 ——要求一定是銳角
        zjfwcs=0     ;//%&&增加范圍次數
         while (point_num == 0)
		 {   //% * 搜索相點
            for (j=1;j<=(M-1);j++)//j = 1 : (M-1)
			{ 
				if (fabs(j-i)<=(P-1))     //%&&候選點距當前點太近，跳過！
				{  continue;  }
                
                //%*計算候選點與當前點的距離
                dnew = 0;
                for (k=1;k<=m;k++)//k = 1 : m
				{ dnew = dnew + (Y[k-1][i-1]-Y[k-1][j-1])*(Y[k-1][i-1]-Y[k-1][j-1]);}
                  dnew = sqrt(dnew);
                
                if ((dnew < min_eps)||( dnew > max_eps ))  // %&&不在距離范圍，跳過！
				{continue;  }           
                
                               
                //%*計算夾角余弦及比較
                DOT = 0;
                for (k=1;k<=m;k++)//k = 1 : m
				{   DOT = DOT+(Y[k-1][i-1]-Y[k-1][j-1])*(Y[k-1][i-1]-Y[k-1][old_Loc_DK]);}
                CTH = DOT/(dnew*DK1);
                
                if (acos(CTH) > (3.14151926/4) )    // %&&不是小于45度的角，跳過！
				{continue;}
                
                if (CTH > cos_sita)  // %&&新夾角小于過去已找到的相點的夾角，保留
				{   cos_sita = CTH;
                    Loc_DK = j;
                    DK = dnew;
                }

                point_num = point_num +1;
                
			}        
        
            if (point_num <= min_point)
			{
				max_eps = max_eps + dlt_eps;
               zjfwcs =zjfwcs +1;
               if (zjfwcs > MAX_CISHU)    //%&&超過最大放寬次數，改找最近的點
			   {   DK = 1.0e+100;
                   for  (ii=1;ii<=(M-1);ii++) //ii = 1 : (M-1)
				   {  if (fabs(i-ii)<=(P-1))      //%&&候選點距當前點太近，跳過！
					  {continue;}    
                      d = 0;
                      for (k=1;k<=m;k++)//k = 1 : m
					  {   d = d + (Y[k-1][i-1]-Y[k-1][ii-1])*(Y[k-1][i-1]-Y[k-1][ii-1]);}
                      
                      d = sqrt(d);
        
                      if ((d<DK)&&(d>min_eps)) 
					  { DK = d;
                         Loc_DK = ii;
                      }
				   }
                   break; 
               }
               point_num = 0 ;     //%&&擴大距離范圍后重新搜索
               cos_sita = 0;
            }
		 }
	}

//%取平均得到最大李雅普諾夫指數
//lambda_wolf=sum(lmd)/length(lmd);
	*lambda_wolf=sum(lmd,M-2)/(M-2);

TwoArrayFree(Y);
free(lmd);
  
    
}