/*vdtwrecoge.c****************************************************************/

#include "vdtwrecoge.h"
#include "math.h"
double vdistance[DTWMAXNUM][DTWMAXNUM]; /*保存距離*/
double vdtwpath[DTWMAXNUM][DTWMAXNUM]; /*保存路徑*/

/*****************************************************************************/
/* VectorAdd(mvector *A,mvector *B,mvector *C)
/* 向量加法，結(jié)果存放在C中
/*****************************************************************************/
void VectorAdd(mvector *A,mvector *B,mvector *C)
{
	if(A->num != B->num)
		return;
	int dim = A->num;
	for(int i = 0; i < dim; i++)
	{
		C->data[i] = A->data[i] + B->data[i];
	}
	C->num = dim;
}
/*****************************************************************************/
/*VectorInit(mvector *A,int dim) 向量初始化
/* dim 維數(shù)
/******************************************************************************/

void VectorInit(mvector *A,int dim)
{
	if((dim > MAX_DIM)||(A == NULL))
		return ;
	for(int i = 0; i < dim; i++)
	{
		A->data[i] = 0.0;
	}
	A->num = dim;

}

/*****************************************************************************/
/*VectorMul(mvector *A,double mul) 向量數(shù)乘
/* 
/*****************************************************************************/
void VectorMul(mvector *A,double mul)
{
	for(int i = 0; i<MAX_DIM; i++)
	{
		A->data[i] *= mul;

	}

}
#define SQUARE(x) x*x
/*****************************************************************************/
/* VectorDistance,用于求兩個向量之間的距離，歐氏距離
/* A和B是待求距離的兩個向量
/*****************************************************************************/
double VectorDistance(mvector *A,mvector *B)
{
	double total = 0.0;
	double *vec1 = A->data;
	double *vec2 = B->data;
	double temp;
	int dim = A->num;
	if(A->num != B->num)
		return -1.0;
	for(int i = 0; i < dim; i++)
	{	
		temp = vec1[i] - vec2[i];
		total += SQUARE(temp);
	}
	total /= dim;
	total = sqrt(total);
	return total;
}
#undef SQUARE


/*****************************************************************************/
/* DTWDistance，求兩個數(shù)組之間的匹配距離
/* A,B分別為第一第二個數(shù)組，I，J為其數(shù)組長度，r為匹配窗口的大小
/* r的大小一般取為數(shù)組長度的1/10到1/30
/* 返回兩個數(shù)組之間的匹配距離,如果返回－1.0，表明數(shù)組長度太大了
/*****************************************************************************/
double VDTWDistanceFun(mvector *A,int I,mvector *B,int J,int r)
{
	int i,j;
	double dist;
	int istart,imax;
	double g1,g2,g3;
	int pathsig=1;/*路徑的標志*/
	int r2=r+ABS(I-J);/*匹配距離*/
	int dim = A->num;    //向量的維數(shù)
	if(A->num != B->num){
		return -1.0;
	}
	if((r2>I-1)||(r2>J-1)){
		return -1.0;
	}

	/*檢查參數(shù)的有效性*/
	if(I>DTWMAXNUM||J>DTWMAXNUM){
		//printf("Too big number\n");
		return -1.0;
	}
	
	/*進行一些必要的初始化*/
	for(i=0;i<I;i++){
		for(j=0;j<J;j++){
			vdtwpath[i][j]=0;
			vdistance[i][j]=DTWVERYBIG;
		}
	}
	
	/*動態(tài)規(guī)劃求最小距離*/
	/*這里我采用的路徑是 -------
	                          . |
	                        .   |
	                      .     |
	                    .       |
	 有可能有更好的前進路徑，我沒有細加論證，不過算法的數(shù)量級是一樣的
	 */
	vdistance[0][0]=(double)2*VectorDistance(&A[0],&B[0]);
	for(i=1;i<=r2;i++){
		vdistance[i][0]=vdistance[i-1][0]+VectorDistance(&A[i],&B[0]);
	}
	for(j=1;j<=r2;j++){
		vdistance[0][j]=vdistance[0][j-1]+VectorDistance(&A[0],&B[j]);
	}
	
	for(j=1;j<J;j++){
		istart=j-r2;
		if(j<=r2)
			istart=1;
		imax=j+r2;
		if(imax>=I)
			imax=I-1;

		
		for(i=istart;i<=imax;i++){
			g1=vdistance[i-1][j]+VectorDistance(&A[i],&B[j]);
			g2=vdistance[i-1][j-1]+2*VectorDistance(&A[i],&B[j]);
			g3=vdistance[i][j-1]+VectorDistance(&A[i],&B[j]);
			g2=MIN(g1,g2);
			g3=MIN(g2,g3);
			vdistance[i][j]=g3;
		}
	}
		
	dist=vdistance[I-1][J-1]/((double)(I+J));
	return dist;
}/*end DTWDistance*/

/*****************************************************************************/
/* DTWTemplate，進行建立模板的工作
/* 其中A為已經(jīng)建立好的模板，我們在以后加入訓(xùn)練樣本的時候，
/* 以已建立好的模板作為第一個參數(shù)，I為模板的長度，在這個模板中不再改變
/* B為新加入的訓(xùn)練樣本，J為B的長度，turn為訓(xùn)練的次數(shù)，在第一次
/* 用兩個數(shù)組建立模板時，r為1，這是出于權(quán)值的考慮
/* temp保存匹配最新訓(xùn)練后的模板，建議temp[DTWMAXNUM]，函數(shù)返回最新訓(xùn)練后模板的長度
/* 如果函數(shù)返回-1，表明訓(xùn)練樣本之間距離過大，需要重新選擇訓(xùn)練樣本，
/* tt為樣本之間距離的閾值，自行定義
/* rltdistance保存距離，第一次兩個數(shù)組建立模板時可以隨意賦予一個值，
/* 后面用前一次返回的值賦予該參數(shù)
/*****************************************************************************/
int VDTWTemplate(mvector *A,int I,mvector *B,int J,mvector *temp,int turn,double tt,double *rltdistance)
{
	double dist;
	int i,j;
	int pathsig=1;
	int dim = A->num;
	dist=VDTWDistanceFun(A,I,B,J,(int)(I/30));
	if(dist>tt){
		printf("\nSample doesn't match!\n");
		return -1;
	}
		
	if(turn==1)
		*rltdistance=dist;
	else{
		*rltdistance=((*rltdistance)*(turn-1)+dist)/turn;
	}
	/*尋找路徑,這里我采用了逆向搜索法*/
	i=I-1;
	j=J-1;
	while(j>=1||i>=1){
		double m;
		if(i>0&&j>0){
			m=MIN(MIN(vdistance[i-1][j],vdistance[i-1][j-1]),vdistance[i][j-1]);
			if(m==vdistance[i-1][j]){
				vdtwpath[i-1][j]=pathsig;
				i--;
			}
			else if(m==vdistance[i-1][j-1]){
				vdtwpath[i-1][j-1]=pathsig;
				i--;
				j--;
			}
			else{
				vdtwpath[i][j-1]=pathsig;
				j--;
			}
		}
		else if(i==0){
			vdtwpath[0][j-1]=pathsig;
			j--;
		}
		else{/*j==0*/
			vdtwpath[i-1][0]=pathsig;
			i--;
		}
	}
	vdtwpath[0][0]=pathsig;
	vdtwpath[I-1][J-1]=pathsig;
	
	/*建立模板*/
	for(i=0;i<I;i++){
		mvector ftemp;
		VectorInit(&ftemp,dim);
		int ntemp=0;
		for(j=0;j<J;j++){
			if(vdtwpath[i][j]==pathsig){
				VectorAdd(&ftemp,&B[j],&ftemp);
				ntemp++;
			}
		}
		VectorMul(&ftemp,1.0/(double)ntemp);
		//temp[i]=(A[i]*turn+ftemp)/((double)(turn+1));/*注意這里的權(quán)值*/
		VectorMul(&A[i],turn);
		VectorAdd(&ftemp,&A[i],&temp[i]);
		VectorMul(&temp[i],1.0/(double)(turn+1));
	}
	
	
	return I;/*返回模板的長度*/
}//end DTWTemplate