/*Fisher.C
*recognition of handwritten numerals with Fishier Linear Classifier
*/
/*
*說明：一個數字16X16像素，每個像素是灰度值，
*歸一化為-1到1的雙精度浮點數
*
*struct Sample
*{
*	double data[RAWDATADIM];	//一個點的數據為8 byte
*	double feature[FEADIM];		//特征向量
*	int trueClass;	//真實類別
*	int classifiedClass[10]; //被識別的分類，對應位置1
*	int isClassified; // 1 for true and -1 for false
*};
*例如 6 對應的classifiedClass[10]為：-1-1-1-1-1-11-1-1-1
*
*/

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

#define TRAININGSNUM	 7291	//訓練樣本數
#define TESTINGSNUM		 2007	//測試樣本數
#define RAWDATADIM		 256	//每個數字的像數點數
#define FEADIM			 60	//提取的特征維數

#define SUCCESSFUL		 0
#define FAIL			 -1

struct Sample
{
	double data[RAWDATADIM];	//一個點的數據為8 byte
	double feature[FEADIM];		//特征向量
	int trueClass;	//真實類別
	int classifiedClass[10]; //被識別的分類，對應位置1
	int isClassified; // 1 for true and -1 for false
};

struct Sample trainingSample[TRAININGSNUM];	//訓練樣本
struct Sample testingSample[TESTINGSNUM];	//測試樣本

int points[FEADIM][2];	//特征點
double dirToMap[FEADIM][1];	//映射向量
double meanC1, meanC2;	//第一類，第二類映射后的均值
int nC1, nC2;	//兩類的訓練樣本數

/*******************************************************************
 *兩個矩陣相乘
 *參數：A×B＝C
 *      A[m,n]，B[n,k]，C[m,k]
 */
void brmul(a, b, m, n, k, c)
int m, n, k;
double a[], b[], c[];
{
	int i, j, l, u;
	for (i=0; i<=m-1; i++)
		for (j=0; j<=k-1; j++)
		{ 
			u = i*k+j; 
			c[u] = 0.0;
			for (l=0; l<=n-1; l++)
				c[u] = c[u] + a[i*n+l]*b[l*k+j];
		}
		
	return;
}

/****************************************************************
 *矩陣求逆
 *參數：A：矩陣
 *      n：矩陣地維數
 */
int brinv(a, n)
int n;
double a[];
{ 
	int *is, *js, i, j, k, l, u, v;
	double d, p;
	is = (int *) malloc(n*sizeof(int));
	js = (int *) malloc(n*sizeof(int));

	for (k=0; k<=n-1; k++)
	{ 
		d=0.0;
		for (i=k; i<=n-1; i++)
			for (j=k; j<=n-1; j++)
			{ 
				l=i*n+j; 
				p=fabs(a[l]);
				if (p>d) 
				{ 
					d=p; 
					is[k]=i; 
					js[k]=j;
				}
			}

		if (d+1.0==1.0)
		{ 
			free(is); 
			free(js); 
			printf("err**not inv\n");
			return(0);
		}
		if (is[k]!=k)
			for (j=0; j<=n-1; j++)
			{ 
				u=k*n+j; 
				v=is[k]*n+j;
				p=a[u]; 
				a[u]=a[v]; 
				a[v]=p;
			}
        
		if (js[k]!=k)
			for (i=0; i<=n-1; i++)
			{ 
				u=i*n+k; 
				v=i*n+js[k];
				p=a[u]; 
				a[u]=a[v]; 
				a[v]=p;
			}
				
		l=k*n+k;
		a[l]=1.0/a[l];
		for (j=0; j<=n-1; j++)
			if (j!=k)
			{ 
				u=k*n+j; 
				a[u]=a[u]*a[l];
			}
				
		for (i=0; i<=n-1; i++)
			if (i!=k)
				for (j=0; j<=n-1; j++)
					if (j!=k)
					{ 
						u=i*n+j;
						a[u]=a[u]-a[i*n+k]*a[k*n+j];
					}
        
		for (i=0; i<=n-1; i++)
			if (i!=k)
			{ 
				u=i*n+k; 
				a[u]=-a[u]*a[l];
			}
	}
    
	for (k=n-1; k>=0; k--)
	{ 
		if (js[k]!=k)
			for (j=0; j<=n-1; j++)
			{ 
				u=k*n+j; 
				v=js[k]*n+j;
				p=a[u]; 
				a[u]=a[v]; 
				a[v]=p;
			}
        
		if (is[k]!=k)
			for (i=0; i<=n-1; i++)
			{ 
				u=i*n+k; 
				v=i*n+is[k];
				p=a[u]; 
				a[u]=a[v]; 
				a[v]=p;
			}
	}
	
	free(is); 
	free(js);
	return(1);
}

/*
*retrieveSample (s, n, dfname, cfname)

*funtion to retrieve data from binary data file
*s: 樣本數組的頭指針，trainingSample or testingSample
*n: number of samples
*dfname: file name of sample data
*cfname: file name of information on true class
*
*/
int retrieveSample (s, n, dfname, cfname)
struct Sample *s;
int n;
char *dfname, *cfname;
{
	FILE *fdf, *fcf;
	
	int sclass[10];

	int i, j;
	//open datafile
	if	((fdf = fopen (dfname, "rb")) == NULL)
	{
		printf ("Can't open file %s\n", dfname);
		return FAIL;
	}
	
	//open classfile
	if ((fcf = fopen (cfname, "rb")) == NULL)
	{
		printf ("Can't open file %s\n", cfname);
		return FAIL;
	}
	
	for (i = 0; i < n; i++)
	{
		fread (s[i].data, sizeof (double), RAWDATADIM, fdf);//retrieve gray data
		fread (sclass, sizeof (int), 10, fcf);
		for (j = 0; j < 10; j++)
		{
			(*(s + i)).classifiedClass[j] = -1; //initiate -1
			if (sclass[j] == 1)
				(*(s + i)).trueClass = j;  //load trueClass
		}
		(*(s + i)).isClassified = -1;//have not been classified

	}

	fclose (fdf);
	fclose (fcf);
	return SUCCESSFUL;
}

/*
int selPoints (serial1,serial2)
選取特征點
serial1:第一類序號
serial2:第二類序號
*/

int selPoints (serial1,serial2)

int serial1;
int serial2;

{
	double sumC1[RAWDATADIM], sumC2[RAWDATADIM],diff[RAWDATADIM];
	int i,j;
	int n1, n2;

	double uBound, lBound;
	int index;

	n1 = n2 = 0;
	for (i = 0; i < RAWDATADIM; i++)
	{
		sumC1[i] = 0;
		sumC2[i] = 0;
	}

	for (i = 0; i < TRAININGSNUM; i++)
	{
		if (trainingSample[i].trueClass == serial1)
		{
			n1++;
			for (j = 0; j < RAWDATADIM; j++)
				sumC1[j] += trainingSample[i].data[j];
		}
		else if(trainingSample[i].trueClass == serial2) 
		{
			n2++;
			for (j =0; j < RAWDATADIM; j++)
				sumC2[j] += trainingSample[i].data[j];
		}
	}

	if (n1 == 0 || n2 == 0)
		return FAIL;

	for (i = 0; i < RAWDATADIM; i++)
	{
		sumC1[i] = sumC1[i]/n1;
		sumC2[i] = sumC2[i]/n2;
		diff[i] = fabs(sumC1[i] - sumC2[i]);
	}
	
	//下面尋找第一類和第二類偏差最小的點為特征點

	lBound = 0; uBound = 3;
	//lBound指示當前一輪的最大值，uBound指示前一輪的最大值
	for (j = 0; j < FEADIM; j++)
	{
		for (i = 0; i < RAWDATADIM; i++)
			if (lBound < diff[i] && uBound > diff[i])
			{
				lBound = diff[i];//find the Maximum
				index = i;
			}
		uBound = lBound;
		lBound = 0;
		points[j][0] = index / ((int)sqrt (RAWDATADIM));
		points[j][1] = index % ((int)sqrt (RAWDATADIM));
	}
	return SUCCESSFUL;
}


/*
void genFeature (s,n)
特征提取
s:樣本結構struct Sample的指針， trainingSample or testingSample
n:number of samples
*/
void genFeature (s,n)
struct Sample *s;
int n;

{
	int i, j, t;

	for (i = 0; i < n; i++)
	{
		for (j = 0; j < FEADIM; j++)
		{
			t = points[j][0] * (int)sqrt(RAWDATADIM) + points[j][1];
			s[i].feature[j] = s[i].data[t];
		}
	}
	return;
}


/*
int Fisher (serial1,serial2)

serial1第一類的序號，serial2第二類的序號
1、返回投影向量double dirToMap[FEADIM][1]，把多維空間映射到一維空間
2、計算meanC1,meanC2，即兩類的均值
3、計算nC1,nC2，兩類的樣本個數
*/
int Fisher (serial1, serial2)
int serial1;
int serial2;

{
	struct Sample *s;
	double meanVector1[FEADIM], meanVector2[FEADIM], diffVector[FEADIM];
	double sumC1[FEADIM], sumC2[FEADIM];

	double SW[FEADIM][FEADIM];//樣本總類內離散度矩陣
	double tempMtrx[FEADIM][FEADIM];
	double tempV1[FEADIM][1];
	double tempV2[1][FEADIM];
	double tempV[1][1];

	int i, j, k;
	int flag;

	s = trainingSample;

	nC1 = nC2 = 0;//初始化
	
	for (j = 0; j < FEADIM; j++)
	{
		sumC1[j] = 0;
		sumC2[j] = 0;
		for (k = 0; k < FEADIM; k++)
			SW[j][k] = 0;
	}

	for (i = 0; i < TRAININGSNUM; i++)
	{
		if (s[i].trueClass == serial1)
		{
			++nC1;//第一類樣本數增1
			for (j = 0; j < FEADIM; j++)
			{
				sumC1[j] += s[i].feature[j]; 
			}
		}
		if (s[i].trueClass == serial2) 
		{
			++nC2;//第二類樣本數增1
			for (j = 0; j < FEADIM; j++)
			{
				sumC2[j] += s[i].feature[j]; 
			}
		}
	}

	if (nC1 == 0 || nC2 == 0)//某一類沒有樣本，返回錯誤
	{
		return FAIL;
	}

	for (j = 0; j < FEADIM; j++)
	{
		meanVector1[j] = sumC1[j]/nC1;
		meanVector2[j] = sumC2[j]/nC2;
		diffVector[j] = meanVector1[j] - meanVector2[j];
	}

	for (i = 0; i < TRAININGSNUM; i++)
	{
		if (s[i].trueClass == serial1 || s[i].trueClass == serial2)
		{
			if (s[i].trueClass == serial1)
				for (j = 0; j < FEADIM; j++)
				{
					tempV1[j][0]  = s[i].feature[j] - meanVector1[j];
					tempV2[0][j] = tempV1[j][0];
				}
			if (s[i].trueClass == serial2) 
				for (j = 0; j < FEADIM; j++)
				{
					tempV1[j][0] = s[i].feature[j] - meanVector2[j];
					tempV2[0][j] = tempV1[j][0];
				}

		//向量求積，得一個樣本的離散度矩陣
		//	tempMtrx[FEADIM][FEADIM] = tempV1 * tempV2 
		
			brmul(tempV1,tempV2,FEADIM,1,FEADIM,tempMtrx);
		
			for (j = 0; j < FEADIM; j++)
				for (k = 0; k < FEADIM; k++)
					SW[j][k] =SW[j][k] + tempMtrx[j][k];
		}

	}

	//歸一化SW矩陣，以免求逆后數字太小而引起舍入誤差的增加
	for (j = 0; j <FEADIM; j++)
		for (k = 0; k < FEADIM; k++)
			SW[j][k] = SW[j][k] / (nC1 + nC2);


	//下面求dirToMap, meanC1, meanC2
	//dirToMap = SW- * (meanVector1 - meanVector2)
	//meanC1 = T(dirToMap) * meanVector1
	//meanC2 = T(dirToMap) * meanVector2
	for (j = 0; j < FEADIM; j++)
		tempV1[j][0] = diffVector[j];
	
	//矩陣求逆
	flag = brinv(SW,FEADIM);
	if (flag != 1)
		return FAIL;

	brmul(SW, tempV1, FEADIM, FEADIM, 1, dirToMap);//得到dirToMap

	for (j = 0; j < FEADIM; j++)
	{
		tempV2[0][j] = dirToMap[j][0];//轉置
		tempV1[j][0] = meanVector1[j];
	}
	brmul(tempV2, tempV1, 1, FEADIM, 1, tempV);
	meanC1 = tempV [0][0];//得到meanC1

	for (j = 0; j < FEADIM; j++)
	{
		tempV1[j][0] = meanVector2[j];
	}
	brmul(tempV2, tempV1, 1, FEADIM, 1, tempV);
	meanC2 = tempV[0][0];//得到meanC2

	return SUCCESSFUL;
}


/*
classifier (s)
s: certain testing sample
classifying using
dirToMap ,meanC1,meanC2, feature
分類時把多類問題轉化為兩類問題，可能性最大的類作為判決結果
如：0類和1類比較，可能性大的再和3類比較，依次類推
*/
int classifier (s)
struct Sample *s;
{
	
	int serial1;
	int serial2;

	int flag;//成功與否的標識
	double critical;//判決臨界值
	double tempV[1][1];
	double tempV1[FEADIM][1],tempV2[1][FEADIM];
	int count, j;
	
	serial1 =0;//第一類初始化為0

	for (count = 1; count < 10; count++)
	{
		serial2 = count;

		//選出特征點
		flag = selPoints (serial1, serial2);
		if (flag == FAIL)
			return FAIL;

		//提取特征
		genFeature (trainingSample, TRAININGSNUM);
		genFeature (testingSample, TESTINGSNUM);

		//計算分類器
		if (Fisher (serial1, serial2) != SUCCESSFUL)
			return FAIL;


		//開始判別

		if (nC1 ==0 || nC2 ==0)
			return FAIL;

		critical = (meanC1 * nC1 + meanC2 * nC2) / (nC1 + nC2);

		for (j = 0; j < FEADIM; j++)
		{
			tempV2[0][j] = dirToMap[j][0];//轉置
			tempV1[j][0] = (*s).feature[j];
		}

		brmul(tempV2, tempV1, 1, FEADIM, 1, tempV);
	
		if (tempV[0][0] < critical)
			serial1 = serial2;
	}
				
	(*s).isClassified = 1;
	(*s).classifiedClass[serial1] = 1;

	return SUCCESSFUL;
}


/*
void ftestingSample (resultf)
char *resultf：結果輸出文件名，
調用classifier ()，作判決，輸出結果
*/
void ftestingSample (resultf)

char *resultf;
{
	struct Sample *s;
	FILE *rf;
	
	int flag;
	int i, j, k, tempClass;
	int rightNum, errorNum;
	double ratio;
		
	s = testingSample;

	if	((rf = fopen (resultf, "a")) == NULL)
	{
		printf ("Can't open file %s\n", resultf);
		return;
	}

	errorNum = 0;
	for (i = 899; i < TESTINGSNUM; i++)
	{	
		flag = classifier (s + i);//判第i個樣本分類
		if (flag == FAIL)
			return;
			
		fprintf (rf, "Sample %d (%d):\n", i, s[i].trueClass);
		fprintf (rf, "Classified to :");
		
		printf ("Sample %d (%d):\n", i, s[i].trueClass);
		printf ("Classified to :");

		k = 0;

		for (j = 0; j < 10; j++)
		{
			if (s[i].classifiedClass[j] == 1)
			{
				k++;
				tempClass = j;
				fprintf (rf, "%d ", j);
				printf ("%d ", j);
			}
		}
		
		if (k !=1 || tempClass != s[i].trueClass)
		{
			fprintf (rf, "\nError\n\n");
			printf ("Error\n");
			errorNum++;
		}	
		else
		{
			fprintf (rf, "\n\n\n");
			printf ("\n");
		}	
	}
	
	//統計			
	rightNum = TESTINGSNUM - errorNum;
	ratio =  (double)rightNum/(double)TESTINGSNUM;
	fprintf (rf, "\n\nStatistics\n");
	fprintf (rf, "------------------------------------------\n\n");
	fprintf (rf,"Right :%d		Error:%d	Ratio:%f\n",
		rightNum, errorNum, ratio);

	fclose (rf);
	return;
}


void main (void)
{
	char *dfname1 = "D:\\usps\\trainData.bina";
	char *cfname1 = "D:\\usps\\trainDataResult.bina";
	char *dfname2 = "D:\\usps\\testData.bina";
	char *cfname2 = "D:\\usps\\testDataResult.bina";
	char *resultf = "D:\\usps\\result.txt";

	int flag;

	//提取數據
	flag = retrieveSample (
		trainingSample, TRAININGSNUM, dfname1, cfname1);
	if (flag != SUCCESSFUL) return;
	
	flag = retrieveSample (
		testingSample, TESTINGSNUM, dfname2, cfname2);
	if (flag != SUCCESSFUL) return;

	//測試樣本
	ftestingSample (resultf);
}