#include "rbf_cloud.h"

/***********生命函數****************/
void write_matrix(char* textname, double * matrix [],int i,int j);
void write_matrix(char* textname, int * matrix [],int i,int j);
void drop_Matrix_Space();


/************函數實現*********************/
void center(double * input){
	double t=100;
	for(int i=1;i<=7;i++){
		t=100;
		for(int j=1;j<=A_C[i-1];j++)
			if(fabs(Ex[i][j]-input[i])<t){
				c_p[i]=j;
				t=fabs(Ex[i][j]-input[i]);
			}
		}
}

// 寫到文件
void write_matrix(char* textname, double * matrix [],int i,int j){
	fp=fopen(textname,"w");
	if(fp==NULL){
		cout<<"cannot open file"<<endl;
		exit(0);
	}
	for(int ii=1;ii<=i;ii++)
	{
		for(int jj =1;jj<=j;jj++)
		{
			fprintf(fp,"%lf",matrix[ii][jj]);
			fputc(' ',fp);
		}
		fputc('\n',fp);
	}
	fclose(fp);

}
void write_matrix(char* textname, int * matrix [],int i,int j){
	fp=fopen(textname,"w");
	if(fp==NULL){
		cout<<"cannot open file"<<endl;
		exit(0);
	}
	for(int ii=1;ii<=i;ii++)
	{
		for(int jj =1;jj<=j;jj++)
		{
			fprintf(fp,"%d",matrix[ii][jj]);
			fputc(' ',fp);
		}
		fputc('\n',fp);
	}
	fclose(fp);

}
void write_matrix(char* textname, vector<double> matrix){
	fp=fopen(textname,"w");
	if(fp==NULL){
		cout<<"cannot open file"<<endl;
		exit(0);
	}

	for(int j =1;j<=matrix.size();j++)
	{
		fprintf(fp,"%lf",matrix[j]);
		fputc(' ',fp);
	}
	fputc('\n',fp);

	fclose(fp);

}

//將已經用matlab生成好的中心讀入 輸入為study.txt 中心輸出到c[][];
void make_center()
{
	int i =1,j=1;
	double data=0.0;
	// open the center.txt
	fp=fopen("center.txt","r");	
	if(fp==NULL){
		cout<<"cannot open file"<<endl;
		exit(0);
	}
	//initial center
	for(i;i<=7;i++){
		for(j=1;j<=A_C[i-1];j++){
			fscanf(fp,"%lf",&data);
			Ex[i][j]=data;
			fgetc(fp);
		}

		for(j=1;j<=A_C[i-1];j++){
			fscanf(fp,"%lf",&data);
			En[i][j]=data;
			fgetc(fp);
		}

		for(j=1;j<=A_C[i-1];j++){
			fscanf(fp,"%lf",&data);
			He[i][j]=data;
			fgetc(fp);
		}

		for(j=1;j<=A_C[i-1];j++){
			fscanf(fp,"%lf",&data);
			Peak[i][j]=data/10000;
			fgetc(fp);
		}
		
	}
	// 形成中心層
	i = 1;
	for(int i1=1; i1<=A_C[0]; i1++)
		for(int i2=1; i2<=A_C[1]; i2++)
			for(int i3=1; i3<=A_C[2]; i3++)
				for(int i4=1; i4<=A_C[3]; i4++)
					for(int i5=1; i5<=A_C[4]; i5++)
						for(int i6=1; i6<=A_C[5]; i6++)
							for(int i7=1; i7<=A_C[6]; i7++)
							{
								C[i][1]=i1;
								C[i][2]=i2;
								C[i][3]=i3;
								C[i][4]=i4;
								C[i][5]=i5;
								C[i][6]=i6;
								C[i][7]=i7;
								i++;
							}
		
	//	中心形成完畢
	fclose(fp);
}

//生成隨機數
double _random(void) 
{ 
	int a; 
	double r; 

	a=rand()%32767; 
	r=(a+0.00)/32767.00; 

	return r; 

} 
//生成正態隨機數
double _sta(double mu,double sigma) 
{ 
	double r1,r2,r; 
			r1=_random(); 
			r2=_random(); 
			r=sqrt(-2*log(r1))*cos(2*PI*r2)*sigma+mu;
	return r;
} 

// 計算樣本的輸出向量
void compute_2(){
	//模糊度k=5
	int k =5;
	double t=-100; int position =0;
	for(int i=1;i<C_BUF;i++){
		//為每一維生成五個正態隨機數，期望是En[],方差是He[]
		for(int j =1;j<=7; j++){
			for(int ii=1;ii<=5;ii++){
				Enk[j][ii]=_sta(En[j][C[i][j]],He[j][C[i][j]]);
			}
		}
		//計算在這個中心上的輸出值
		double sum = 0.0;
		u[s_input][i] = 0.0;
		for(k=1;k<=5;k++){
			sum = 0.0;
			for(j=1;j<=7; j++){
				sum+=Peak[j][C[i][j]]*((sample[j]-Ex[j][C[i][j]])*(sample[j]-Ex[j][C[i][j]])/(Enk[j][k])*(Enk[j][k]));
			}
			sum/=5;
			u[s_input][i]+=exp(-sum/2);
		}
//		u[s_input][i]/=5;
		if(u[s_input][i]>t){
			t=u[s_input][i];
			position = i;
		}
	}
	
	cout<<"最大輸出："<<t<<endl;
	cout<<"最近中心 ";
	for(int j=1; j<=7;j++)
	{
		cout<<C[position][j]<<" ";
	}
	cout<<endl;
}

// 計算隱藏節點的輸出矩陣（在訓練時使用）
void compute_1(){
	//從訓練文件讀入一個樣本放入sample1[]中。
	fp=fopen("study.txt","r");
	if(fp==NULL){
		cout<<"cannot open file"<<endl;
		exit(0);
	}
	s_input = 1;
	double temp;
	while(!feof(fp)){
		for(int i=1; i<=7; i++){
			fscanf(fp,"%lf",&temp);
			sample[i]=temp;
			fgetc(fp);
		}
		center(sample);
		//計算此樣本的輸出向量
		compute_2();
		//讀取下一條樣本數據
		fscanf(fp,"%d",output1[s_input]+1);
		while(fgetc(fp)!='\n'&&!feof(fp));
		s_input ++;
	}
	s_input--;
	fclose(fp);

}


double randnumber()//-0.1~0.1的隨機函數 
{
	return (sin(rand())/10) ;
}

//用LMS法求權值
void compute_3(){

	
	//初始化
	double e = 0.0;
	for(int j=0; j<C_BUF; j++)
		w.push_back(randnumber());
	double temp = 0.0;
	//對每一個輸入向量的中心輸出進行權值逼近
	for(int i=1;i<=s_input;i++){
		temp = 0.0;
		for(int ii=1; ii<C_BUF; ii++){
			if(fabs(u[i][ii])<1e-6)
				;
			else temp+=w[ii]*u[i][ii];
		}
		e = output1[i][1]-temp;
		cout<<e<<endl;
		for(int j=1; j<C_BUF; j++){
			if(fabs(u[i][j])<1e-6)
				;
			else w[j] = w[j]+2*miu*e*u[i][j];
		}
	}
		write_matrix("quanzhi.txt",w);
	
}

//開辟矩陣空間
void creat_Matrix_Space(){

	Ex= new double*[N_BUF];
	for(int i=0;i<N_BUF;i++)
		Ex[i] = new double[5];
	En= new double*[N_BUF];
	for( i=0;i<N_BUF;i++)
		En[i] = new double[5];
	He= new double*[N_BUF];
	for( i=0;i<N_BUF;i++)
		He[i] = new double[5];
	Peak= new double*[N_BUF];
	for( i=0;i<N_BUF;i++)
		Peak[i] = new double[5];
	Enk = new double *[N_BUF];
	for( i=0;i<N_BUF;i++)
		Enk[i] = new double[6];
	C= new int*[C_BUF];
	for( i=0;i<C_BUF;i++)
		C[i]=new int[N_BUF];
	/*
	w= new double*[C_BUF];
	for(i=0;i<C_BUF;i++)
		w[i] = new double[O_BUF];
		*/

	sample = new double [8+1];

	output1= new int*[S_BUF];
	for(i=0;i<S_BUF;i++)
		output1[i]= new int[O_BUF];
	output2= new int*[S_BUF];
	for(i=0;i<S_BUF;i++)
		output2[i]= new int[O_BUF];
	u= new double*[max_node];
	for(i=0;i<max_node;i++)
		u[i]= new double[C_BUF];
}

//回收空間
void drop_Matrix_Space(){

	delete []Ex;
	delete []En;
	delete []Peak;
	delete []He;
	delete []C;
	delete []sample;
	delete []u;
	delete []output1;
}

// 調用MATLAB畫結果狀態圖

void plot(){
	//混合編程
	Engine *ep; //定義Matlab 引擎指針　 
    if (!(ep=engOpen(NULL))) //測試是否啟動 Matlab 引擎成功。 
    {
		cout <<"Can't start Matlab engine!" <<endl;
		exit(1);
    }
	//定義實數數組。  
	mxArray *xx = mxCreateDoubleMatrix(1,T_BUF-1, mxREAL);  
	mxArray *tt = mxCreateDoubleMatrix(1,T_BUF-1, mxREAL); 
	mxArray *oo = mxCreateDoubleMatrix(1,T_BUF-1, mxREAL);

	double m_i[T_BUF-1], m_o[T_BUF-1], o[T_BUF-1];

	for(int i=0;i<T_BUF-1;i++){
		m_i[i]=matrix_in[i+1][1];
		m_o[i]=matrix_out[i+1][1];
		o[i]=output2[i+1][1];
	}

	memcpy(mxGetPr(xx), m_i, (T_BUF-1)*sizeof(double));// 將數組 x 復制到 mxarray 數組 xx 中。 
	memcpy(mxGetPr(tt), m_o, (T_BUF-1)*sizeof(double));// 將數組 y 復制到 mxarray 數組 yy 中。 
	memcpy(mxGetPr(oo), o, (T_BUF-1)*sizeof(double));// 將數組 z 復制到 mxarray 數組 zz 中。

	engPutVariable(ep, "xx",xx);// 將 mxArray 數組 xx 寫入到 Matlab 工作空間，命名為 xx 
	engPutVariable(ep, "oo",oo);// 將 mxArray 數組 yy 寫入到 Matlab 工作空間，命名為 yy
	engPutVariable(ep, "tt",tt);// 將 mxArray 數組 zz 寫入到 Matlab 工作空間，命名為 zz 

    // 向 Matlab 引擎發送畫圖命令。 plot 為 Matlab 的畫圖函數
	engEvalString(ep, "title('測試結果圖')");
	engEvalString(ep, "xlabel('X坐標'),ylabel('Y坐標');hold on");
	engEvalString(ep, "plot(xx, oo,'r',xx, tt,'d');");

	mxDestroyArray(xx); // 銷毀 mxArray 數組 xx 和 yy 。 
	mxDestroyArray(tt);
	mxDestroyArray(oo);

	cout<<"繪圖完畢，按任一鍵退出!" <<endl;
	cin.get();
	engClose(ep); // 關閉 Matlab 引擎。 
	return;
}

//測試

void rbf_test()
{
	FILE * fp_test=NULL;
	fp_test=fopen("test.txt","r");
	if(!fp_test)
	{
		cout<<"cannot open file"<<endl;
		exit(0);
	}

	int in=1;

	// 開辟輸入輸出空間
	matrix_in = new double *[T_BUF];
	for(int i =1; i<T_BUF; i++)
		matrix_in[i]=new double[N_BUF];
	matrix_out= new double *[T_BUF];
	for( i =1; i<T_BUF; i++)
		matrix_out[i]=new double[O_BUF];


	while (!feof(fp_test))
	{	// 讀入一個實驗數據
		for(int j=1; j<N_BUF; j++)
		{
			fscanf(fp_test,"%lf",input+j);
			matrix_in[in][j] = input[j];
		}

		for( j=1;j<O_BUF;j++){
			fscanf(fp_test,"%d",output2[in]+j);
		}
		while(fgetc(fp_test)!='\n'&&!feof(fp_test));

		// 計算隱層節點輸出

		for(int i=1;i<C_BUF;i++){
			//為每一維生成五個正態隨機數，期望是En[],方差是He[]
			
			for(int j =1;j<=7; j++){
				for(int ii=1;ii<=5;ii++){
					int c = C[i][j];
					Enk[j][ii]=_sta(En[j][c],He[j][c]);
				}
			}
			//計算在這個中心上的輸出值
			double sum = 0.0;
			u[in][i] = 0.0;
			for(int k=1;k<=5;k++){
				sum = 0.0;
				for(j=1;j<=7; j++){
					int c = C[i][j];
					sum+=Peak[j][c]*((matrix_in[in][j]-Ex[j][c])*(matrix_in[in][j]-Ex[j][c])/(Enk[j][k]*Enk[j][k]));
				}
				u[in][i]+=exp(-sum/2);
			}
			u[in][i]/=5;
		}
		//*****************************

		double e = 0.0;
		double result =0.0;
		
		for(int jj=1;jj<C_BUF;jj++)
			result+=w[jj]*u[in][jj];

		matrix_out[in][1] = result;
		
		in++;
	}
	// 畫圖
	plot();

	delete [] matrix_in;
	delete [] matrix_out;
	fclose(fp_test);
}

void init(){
	w.clear();
	w.push_back(0);
	//讀入權值
	fp=fopen("quanzhi.txt","r");
	if(!fp)
	{
		cout<<"cannot open file"<<endl;
		exit(0);
	}
	int i =1;
	double temp;
	while(!feof(fp)){
		fscanf(fp,"%lf",&temp);
		w.push_back(temp);
		fgetc(fp);

		i++;
	}

}

void main(){
	
	creat_Matrix_Space();
	
	//用聚集法形成中心
	make_center();
	compute_1();
	compute_3();

	//訓練結束
	

	//開始試驗
	//先建立中心
	//make_center();
	init();
	rbf_test();
	//測試結束
	drop_Matrix_Space();
}