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

//---------------------------------------------------------------------------
// 計算絕對值
double ABS(double x)
{
    double y;
    if (x>=0)
        y = x;
    else
        y = -x;
    return y;
}
//---------------------------------------------------------------------------

// 計算數組最大值
double MAX_VECTOR(double *p_vector,
                  int len_vector)
{
    int i;    
    double max_value = *p_vector;
    for (i=0; i<len_vector; i++)
    {   if (*(p_vector+i)>max_value)
	{
		max_value = *(p_vector+i);
	}
    }
    return max_value;
}




//求最短距離及向量
void min_dist(double *pdist,int n,int *place,double *pmindist,int i)
{  int k,j;    
   double min_value ;
   
   for (k=0;k<n;k++)
   {
	   if (*(pdist+k)!=0) 
	   {min_value=*(pdist+k);j=k+1;
	   break;}
   }

  for (k=0; k<n; k++)
  {   
	  if ((*(pdist+k)!=0)&&(*(pdist+k)<min_value))
	  {	 min_value = *(pdist+k); j=k+1; }
  }
  
  *(place+i)=j;
  *(pmindist+i)=min_value;

}



//-----------------------------------------------------------------------------
// 計算向量間的距離（最大范數）
double distance(double *Y1,double *Y2,int m1,int tau)
{  int i;
  double d,*pd,d_ij;
  pd = (double*)malloc(m1*sizeof(double));       // 聲明一個長度為 m 的 double 型數組
	for (i=1;i<=m1;i++)
	{
			d=*(Y1+(i-1)*tau)-*(Y2+(i-1)*tau);
			*(pd+i-1)=ABS(d);
	}
	d_ij=MAX_VECTOR(pd,m1);
    free(pd);
	return d_ij;
}


//---------------------------------------------------------------------------

// 定義 C 運算函數
void Cao_FUNCTION( double *pdata,      // 時間序列（列向量）
					int m_min,          // 最小嵌入維數       
					int m_max,          // 最大嵌入維數 
					int tau,            // 時間延遲
					int N ,  		     // 序列長度
					double *pE1,        //用于存放E1
					double *pE2)        //用于存放E2
					
{   int *place,i,j,m;
    double *pdist,*pmindist,*pE,*pEn;

   	pE = (double*)malloc((m_max-m_min+1)*sizeof(double)); 
	pEn = (double*)malloc((m_max-m_min+1)*sizeof(double)); 

 for (m=m_min;m<=m_max;m++)
 {    
     double ad=0,neard, dd=0;
      int n=N-m*tau,d;
      place = (int*) malloc(n*sizeof(int));    //用于存放重構后向量的標號
	  pmindist = (double*)malloc(n*sizeof(double)); 
	 // printf("共有 %d 步,正在計算第 %d 步!\n",m_max,m);
	 	  
    for (i=0;i<n;i++)         //計算第i個向量的距離最近的向量
	{  
		pdist = (double*)malloc(n*sizeof(double));   //存放第i個與每個向量的距離

       for (j=0;j<n;j++)
	   { *(pdist+j)=distance(pdata+i,pdata+j,m,tau);	}
	  
	    min_dist(pdist,n,place,pmindist,i);   //求第i個向量的最近距離及對應向量標號

      //--------------計算a(i,m)-------------------
      d=*(place+i);
      neard=*(pmindist+i);
      ad=ad+distance(pdata+i,pdata+d-1,m+1,tau)/neard;
	  dd=dd+ABS(*(pdata+i+m*tau)-*(pdata+d-1+m*tau));//求E*(d)
      free(pdist);
	}
	
      *(pE+m-m_min)=ad/n;   //求E(d)
      *(pEn+m-m_min)=dd/n; 	 //求E*(d)
      
      free(place);free(pmindist);
 }
//----------求E1(d)&E2(d)-------------------
  for (i=1;i<=m_max-m_min;i++)
  {
	   *(pE1+i-1)=*(pE+i)/(*(pE+i-1));
	   *(pE2+i-1)=*(pEn+i)/(*(pEn+i-1));
	   
  }
  free(pE);free(pEn);
}


//---------------------------------------------------------------------------
void 
mexFunction (int nlhs, mxArray *plhs[],			// 輸出參數個數，及輸出參數數組
			 int nrhs, const mxArray *prhs[])	// 輸入參數個數，及輸入參數數組
{
      double *pdata,*pE1,*pE2;
      int m_min,m_max,tau,N,nRows;

	  if (nrhs!=4) mexErrMsgTxt("需要4個參數!");  //檢查輸入參數的個數

    // 取得輸入參數
    pdata = mxGetPr(prhs[0]);      // 時間序列（列向量）      
    m_min = (int) *mxGetPr(prhs[1]);            // 最小嵌入維數 
	m_max = (int) *mxGetPr(prhs[2]);            // 最大嵌入維數 
    tau = (int) *mxGetPr(prhs[3]);      // 時間延遲   
	N = mxGetM(prhs[0]);              // 序列長度
   
	
    // 為輸出變量分配內存空間
	nRows=m_max-m_min;
	 plhs[0]= mxCreateDoubleMatrix(nRows,1,mxREAL); //用于存放E1
	 plhs[1]= mxCreateDoubleMatrix(nRows,1,mxREAL); //用于存放E2
	
	// 取得輸出參數指針
	pE1 = mxGetPr(plhs[0]);
    pE2 = mxGetPr(plhs[1]);

    // 調用 C 運算函數 (該函數名不能和本文件名重名)
     Cao_FUNCTION(pdata,m_min,m_max,tau,N,pE1,pE2);    
	
}