#include <iostream>
#include <math.h>
#include "bayes.h"
#include <cstdlib>
//計算均值向量U
Matrix getU(const Matrix &X)
{
int d = X.getCol();
int n = X.getRow();
Matrix U(d, 1);        //U=(u1,u2,...,ud)^T，
                //用 d x 1 的矩陣表示
for(int k=0; k<n; k++)    //k 代表第 k 個樣本
   for(int j=0; j<d; j++)   //j 代表樣本的第 j 個特征值
    U.set(j, 0, U.get(j,0)+X.get(k,j));
    //這里要注意對矩陣元素的訪問下標是從 0 開始的
return 1.0/n * U; //U 是所有樣本的均值，故加和后除以樣本總數
}

//計算協方差矩陣E
Matrix getE(const Matrix &X, const Matrix &U)
{
int d = X.getCol(); //X 的列數代表特征空間維數
int n = X.getRow(); //X 的每一行代表一個樣本
Matrix E(d, d);         //E 是 d x d 維矩陣，d 為每個樣本特征空間數
                 //故此處可以用 X 列數來初始化它
double *ar = new double[d];   //此數組用于從矩陣中提取行，創建新矩陣
                 //從而可以進行下面的矩陣運算
for(int k=0; k<n; k++){     //k 代表第 k 個樣本
   for(int j=0; j<d; j++)
    ar[j] = X.get(k, j);    //將 X 的一行元素值放入數組 ar 中

   Matrix Xk(d, 1, ar);     //利用 ar 創建矩陣
   E = E + (Xk - U) * trv(Xk - U); //(X-U) * (X-U)^T
}
delete[] ar; //釋放動態分配的數組空間
//E 是 n 個矩陣( (Xk - U)與其轉置的乘積 )的算術平均
return 1.0/n * E;
}

//貝葉斯判別函數
int bayesFun(const Matrix U[], const Matrix E[],
        const double Pw[],const Matrix &X, int c)
{
double testPw = 0;
for(int w=0; w<c; w++){
   if(Pw[w] <= 0){
    cout << "Wrong Pw[i]\n"; exit(0); }
   testPw += Pw[w];
   if(det(E[w]) == 0){
    cout << "Check input, E[" << w << "]=0!\n";
    exit(0); }
}
if(testPw != 1){
   cout << "Wrong Pw[i]\n"; exit(0); }

//對于貝葉斯判別函數中的 -d/2 * ln(2*PI)項因與 i 無關故可去掉
//另外(X-U)(E^-1)(X-U)^T 的計算結果雖然是一個數，但在這里
//其表示為一個一行一列的矩陣，因此需要對其求行列式值后才能
//賦給左邊的 double maxg
double maxg = -1.0/2 * det( trv(X-U[0]) * inv(E[0]) * (X-U[0]) )
         -1.0/2 * log(det(E[0])) + log(Pw[0]);
int code = 0; //用來返回類別代號

//求使得貝葉斯判別函數取得最大值類
for(int i=1; i<c; i++){
   double g =   -1.0/2 * det( trv(X-U[i]) * inv(E[i]) * (X-U[i]) )
         -1.0/2 * log(det(E[i])) + log(Pw[i]);
   if(g > maxg){
    maxg = g;
    code = i;
   }
}//for i

return code;
}//bayesFun()

//貝葉斯分類器
void bayesDepart(const Matrix X[], const double Pw[],
         const Matrix &x, int c, char* name[])
{
Matrix* U = new Matrix[c];
Matrix* E = new Matrix[c];

//先求出每個類的均值向量和協方差矩陣
for(int i=0; i<c; i++){
   U[i] = getU(X[i]);
   E[i] = getE(X[i], U[i]);
}
//調用貝葉斯判別函數，判斷X所屬類別，輸出類別名稱
cout << name[bayesFun(U,E,Pw,x,c)];
}