#include<stdio.h>
#include<conio.h>
#include<stdlib.h>
#include<math.h>
#include<windows.h>

#define maxiter 80000

//-------------------------------數據類型定義-----------------------
//矩陣結構定義

#define m 10304
#define n 1			//測試時為1，訓練時為7
#define k 5

#define k1 1
#define k2 1000000

double A[m][n];//原始矩陣
double W[m][k];//左乘矩陣
double H[k][n];//右乘矩陣

double H_[k][1];

//-------------------------------隨機矩陣生成函數-------------------

void matrixrand(int tag)
{
#define RANDMAX 1000000000//可以運行前設置，也可=A中數據的最大值+1~……
	int i,j;
	if(tag==0)//tag為0時為左乘矩陣的初始化
		for(i=0;i<m;i++)
			for(j=0;j<k;j++)
				W[i][j]=double(rand()%RANDMAX*1.0/RANDMAX);
	else
		if(tag==1)//tag為1時為右乘矩陣的初始化
			for(i=0;i<k;i++)
				for(j=0;j<n;j++)
					H[i][j]=double(rand()%RANDMAX*1.0/RANDMAX);
#undef RANDMAX
}

//-------------------------------信息讀入函數-----------------------

bool inforead(char *filename)
{
	int i,j;
	FILE *fp;
	fp=fopen(filename,"r");
	if(fp==NULL)
		return FALSE;
	for(i=0;i<n;i++)
		for(j=0;j<m;j++)
			fscanf(fp,"%lf",&A[j][i]);
	fclose(fp);
	return TRUE;
}

//-------------------------------正規化函數-------------------------

void normalization(void )
{
	int i,j;
	double max,min;
	for(i=0;i<m;i++)
	{
		max=min=A[i][0];
		for(j=0;j<n;j++)
		{
			if(A[i][j]<min)
				min=A[i][j];
			else
				if(max<A[i][j])
					max=A[i][j];
		}
		for(j=0;j<n;j++)
			A[i][j]=(A[i][j]-min)*1.0/(max-min);
	}
}

//-------------------------------歸一化H----------------------------

void NormalizeH(void )
{
	int i,j;
	double sum;
	//H歸一化
	for(i=0;i<k;i++)
	{
		sum=0;
		for(j=0;j<n;j++)
			sum+=H[i][j];
		for(j=0;j<n;j++)
			H[i][j]/=sum;
	}
}

//-------------------------------歸一化W----------------------------

void NormalizeW(void )
{
	int i,j;
	double sum;
	//W歸一化
	for(i=0;i<k;i++)
	{
		sum=0;
		for(j=0;j<m;j++)
			sum+=W[j][i];
		for(j=0;j<m;j++)
			W[j][i]/=sum;
	}
}

//-------------------------------計算H的迭代公式--------------------

void updateH(void )
{
	int i,j,con;
	double WT[k][m];//存儲W的轉置
	double Temp_1[k][n];//存儲WT與A的積
	double Temp_2[k][k];//存儲WT與W的積
	double Temp_3[k][n];//存儲Temp_2與H的積+1e-6
	double Temp_4[k][n];//存儲H.*Temp_1
	//轉置操作W
	for(i=0;i<k;i++)
		for(j=0;j<m;j++)
			WT[i][j]=W[j][i];
	//求WT與A的積
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
		{
			Temp_1[i][j]=0;
			for(con=0;con<m;con++)
				Temp_1[i][j]+=WT[i][con]*A[con][j];
		}
	//求WT與W的積
	for(i=0;i<k;i++)
		for(j=0;j<k;j++)
		{
			Temp_2[i][j]=0;
			for(con=0;con<m;con++)
				Temp_2[i][j]+=WT[i][con]*W[con][j];
		}
	//求Temp_2與H的積
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
		{
			Temp_3[i][j]=0;
			for(con=0;con<k;con++)
				Temp_3[i][j]+=Temp_2[i][con]*H[con][j];
			Temp_3[i][j]+=1e-9;
		}
	//求H.*Temp_1
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
			Temp_4[i][j]=H[i][j]*Temp_1[i][j];
	//求一趟循環后的H
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
			H[i][j]=Temp_4[i][j]/Temp_3[i][j];
}

//-------------------------------計算W的迭代公式--------------------

void updateW(void )
{
	int i,j,con;
	double HT[n][k];//存儲H的轉置
	double Temp_01[m][k];//存儲A與HT的積
	double Temp_02[k][k];//存儲H與HT的積
	double Temp_03[m][k];//存儲W與Temp_02的積+1e-6
	double Temp_04[m][k];//存儲W.*Temp_01
	//轉置操作H
	for(i=0;i<n;i++)
		for(j=0;j<k;j++)
			HT[i][j]=H[j][i];
	//求A與HT的積
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
		{
			Temp_01[i][j]=0;
			for(con=0;con<n;con++)
				Temp_01[i][j]+=A[i][con]*HT[con][j];
		}
	//求H與HT的積
	for(i=0;i<k;i++)
		for(j=0;j<k;j++)
		{
			Temp_02[i][j]=0;
			for(con=0;con<n;con++)
				Temp_02[i][j]+=H[i][con]*HT[con][j];
		}
	//求W與Temp_02的積
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
		{
			Temp_03[i][j]=0;
			for(con=0;con<k;con++)
				Temp_03[i][j]+=W[i][con]*Temp_02[con][j];
			Temp_03[i][j]+=1e-9;
		}
	//求W.*Temp_01
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
			Temp_04[i][j]=W[i][j]*Temp_01[i][j];
	//求一趟循環后的W
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
			W[i][j]=Temp_04[i][j]/Temp_03[i][j];
}

//-------------------------------原始算法函數-----------------------

void origNMF(void )
{
	int control=0;
	while(control<maxiter)
	{
		updateH();
		NormalizeH();//歸一化H
		updateW();
//		NormalizeW();//歸一化W
		control++;
	}
}

//-------------------------------判優函數---------------------------

double g_estimate(void )
{
	int i,con;
	double gestimate=0;
	double AWH[m][n];//存儲W*H所生成的矩陣
	for(i=0;i<m;i++)
	{
		AWH[i][0]=0;
		for(con=0;con<k;con++)
			AWH[i][0]+=W[i][con]*H[con][0];
	}
	for(i=0;i<m;i++)
		gestimate+=k1*pow((A[i][0]-AWH[i][0]),2);
	for(i=0;i<k;i++)
		gestimate+=k2*pow((H[i][0]-H_[i][0]),2);
	return gestimate;
}

//-------------------------------測試函數---------------------------

bool test(void )
{
	int min;
	double tempmin;
	double Cgestimate[10];
	/////////////////////////////////////////////////////////////////

	char Wname[10][10]={"W1.txt","W2.txt","W3.txt","W4.txt","W5.txt","W6.txt","W7.txt","W8.txt","W9.txt","W10.txt"};
	char hname[10][10]={"h1.txt","h2.txt","h3.txt","h4.txt","h5.txt","h6.txt","h7.txt","h8.txt","h9.txt","h10.txt"};

	/////////////////////////////////////////////////////////////////
	FILE *fp,*fpresult,*fptemp;
	if(!(fp=fopen("測試.txt","r")))
		return false;
	if(!(fpresult=fopen("結果.txt","w")))
		return false;

	for(int i=0;i<40;i++)
	{
		for(int j=0;j<m;j++)
			fscanf(fp,"%lf",&A[j][0]);
		int ii,jj;

		//計算對應的目標函數
		for(int testtemp=0;testtemp<10;testtemp++)
		{
			fptemp=fopen(Wname[testtemp],"r");
			for(ii=0;ii<m;ii++)
				for(jj=0;jj<k;jj++)
					fscanf(fptemp,"%lf",&W[ii][jj]);
			fclose(fptemp);

			fptemp=fopen(hname[testtemp],"r");
			for(ii=0;ii<k;ii++)
				fscanf(fptemp,"%lf",&H_[ii][0]);
			fclose(fptemp);

			matrixrand(1);//初始化H_矩陣

			for(j=0;j<200;j++)
				updateH();

			Cgestimate[testtemp]=g_estimate();
		}
		
		//找出Cgestimate中最小的值
		tempmin=Cgestimate[0];
		min=0;
		for(int t=1;t<10;t++)
			if(Cgestimate[t]<tempmin)
			{
				min=t;
				tempmin=Cgestimate[t];
			}

		fprintf(fpresult,"%d\n",min);
	}
	fclose(fp);
	fclose(fpresult);
	return true;
}

//-------------------------------主函數-----------------------------

void main()
{
//                  測試

	
	test();

	
/*//					訓練
	bool bo;
	bo=inforead("訓練2.txt");		//具體訓練時需要更改該處的名稱
	if(bo==FALSE)
		printf("文件打開失敗！！！\n");
	//矩陣的初始化
	matrixrand(0);
	matrixrand(1);

	origNMF();

	int i,j;
	FILE *fp;
	fp=fopen("W2.txt","w");			//具體訓練時需要更改該處的名稱
	for(i=0;i<m;i++)
	{
		for(j=0;j<k;j++)
			fprintf(fp,"%lf\t",W[i][j]);
		fprintf(fp,"\n");
	}
	fclose(fp);

	double Havg;
	fp=fopen("h2.txt","w");			//具體訓練時需要更改該處的名稱
	for(i=0;i<k;i++)
	{
		Havg=0;
		for(j=0;j<n;j++)
			Havg+=(H[i][j]*1.0/n);
		fprintf(fp,"%lf\n",Havg);
	}
	fclose(fp);
*/
	printf("OK!\n");
	getch();
}