EM Algorithm

EM_Algorithm.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector> 
#include <cstdlib>
#include <cmath>

#define MAX_TOK 10  
#define PI 3.14159265359

using namespace std;

// ---------------------------- Functions ------------------------------------------

// ======== Load File Function =================
string* StringSplit(string strOrigin, string strTok){
	int     cutAt;                            //자르는위치
	int     index     = 0;                    //문자열인덱스
	string* strResult = new string[MAX_TOK];  //결과return 할변수

	//strTok을찾을때까지반복
	while ((cutAt = strOrigin.find_first_of(strTok)) != strOrigin.npos)
	{
	    if (cutAt > 0)  //자르는위치가0보다크면(성공시)
	    {
	        strResult[index++] = strOrigin.substr(0, cutAt);  //결과배열에추가
	    }
	    strOrigin = strOrigin.substr(cutAt+1);  //원본은자른부분제외한나머지
	}

	if(strOrigin.length() > 0)  //원본이아직남았으면
	{
	    strResult[index++] = strOrigin.substr(0, cutAt);  //나머지를결과배열에추가
	}

	return strResult;  //결과return
/*
This function split string row to the tokenizers
*/
}

vector< vector<float> > LoadFile(string filename){
/* 
This function loads data and convert txt into float matrix 
*/
	string line;
	vector< vector<float> > Hans_Matrix; 
	vector<float> row; 
	ifstream MyStream (filename);

	if (MyStream.is_open()) {
    	while ( getline (MyStream,line) ){
			row.clear();
			string *token = StringSplit(line, "\t");
			// For each row in Matrix 
			for (int idx = 0; idx < 4; idx ++){
				float Str_to_Double = atof(token[idx].c_str()); // String Number 
				row.push_back(Str_to_Double);
			}
			Hans_Matrix.push_back(row);
    	}
    	MyStream.close();
  	}
  	return Hans_Matrix;
}

// ======== Performance Check Function =================

float Performance_Check(vector< vector<float> > My_Answer_R, vector< vector<float> > My_Answer_G, vector< vector<float> > My_Answer_B, 
	vector< vector<float> > R_Matrix, vector< vector<float> > G_Matrix, vector< vector<float> > B_Matrix){
	/* 
		Compute the Accuracy
	*/ 
	double score = 0.0; 
	float Accuracy; 
	bool isPresent (vector <vector<float> >, vector<float>);
	for (int idx = 0; idx < My_Answer_R.size(); idx++){
		vector<float> Current =  My_Answer_R[idx];
		if (isPresent(R_Matrix, Current) ){
			score ++; 
		}
	}

	for (int idx = 0; idx < My_Answer_G.size(); idx++){
		vector<float> Current =  My_Answer_G[idx];
		if (isPresent(G_Matrix, Current) ){
			score ++; 
		}
	}

	for (int idx = 0; idx < My_Answer_B.size(); idx++){
		vector<float> Current =  My_Answer_B[idx];
		if (isPresent(B_Matrix, Current) ){
			score ++; 
		}
	}

	Accuracy = score / 900; 
	return Accuracy;

}

bool isPresent(vector< vector<float> > Vector_Matrix, vector<float> Target){
/* 
	Check if Target Vector is in the Vector Matrix or Not 
*/
	bool result = false;
	for (int idx = 0; idx < Vector_Matrix.size(); idx ++){
		if (Vector_Matrix[idx] == Target){
			result = true;
			break;
		}
		else{
			result = false;
		}
	}

	return result; 
}

// ============================ Clustering Function  ===============

float Distance(vector<float> A, vector<float> B){
/* 
This function compute the distance
A.size = B.size
*/	
	float distance = 0; 
	for (int idx = 0; idx < A.size(); idx ++){
		distance += pow((A[idx] - B[idx]),2);
	}
	return sqrt(distance);

}

vector< vector< vector<float> > > K_Means_Clustering(vector< vector<float> > Train_Matrix, vector< vector<float> > Centroid_Matrix, int Num_Cluster){
	int K = Num_Cluster;
	vector< vector <vector<float> > > Cluster_Map; 

	for (int idx = 0; idx < K; idx ++){
		Cluster_Map.push_back(Centroid_Matrix);
	}
	
	for (int idx =0 ; idx < K; idx ++){
		Cluster_Map[idx].clear();
	}

	for (int idx = 0; idx < Train_Matrix.size(); idx++){
		vector<float> Current = Train_Matrix[idx];
		float min_distance = 100000.0; 
		float distance = 0.0;
		int cluster_idx = 100; 

		for (int temp = 0; temp < K; temp++){
			distance = Distance(Current, Centroid_Matrix[temp]);
			if (distance < min_distance){
				min_distance = distance; 
				cluster_idx = temp;  // Among 0,1,2 if K = 3 
			}
		}
		
		Cluster_Map[cluster_idx].push_back(Current);
		
	}
	return Cluster_Map;


}

// ========== Normal Distribution COmputation Function  ===============

float Power_Compute(float num, int K){
    float result = num;

    for (int idx = 1; idx < K; idx ++){
        result = result * num;
    }
    return result;
}

vector< vector<float> > Centroid_for_Mean(vector< vector<float> > Train_Matrix, int num_centroid){
	int K = num_centroid; 

	// =========== Initialization for each iteration (while goes from here) =============
	//		======== AT STEP 0 ============
	vector< vector <float> > Centroids;  // Declare Centriod 
	vector< vector <float> > Prev_Centroids;  // Declare Centriod 
	vector< vector <float> > Cluster_Sum ;  // Declare Cluster Sum 

	int num_iter = 0; 

	for (int idx = 0; idx < K; idx ++){
		int random_key = rand() % (Train_Matrix.size());
		Centroids.push_back(Train_Matrix[random_key]);
		Cluster_Sum.push_back(Train_Matrix[random_key]);
	}
	
	while (true) {
		// ========== After step 0, Cluster Sum is initialzied as centriod 
		Prev_Centroids = Centroids;

		// =========== Define Counter Cluster ========= 
		vector<int> Count_Cluster; 
		for (int temp = 0; temp < K; temp++){
			Count_Cluster.push_back(0);
		}

		// =========== Find all centroid for each cluster =============
		for (int idx = 0; idx < Train_Matrix.size(); idx++){
			
			// =============== Initialize for each data point ==========
			float distance; // declare distance 
			vector <float> Current = Train_Matrix[idx]; // current target point 
			float Min_Distance = 100000.0;  
			int Cluster_Idx = 100; 
			
			// =============== Measure distance from each centriods to the current point ==========
			// ---------- for each centroid
			for (int idx_cent = 0; idx_cent < Centroids.size(); idx_cent ++){
				// ------ measure distance from each centroid to the current point 
				distance = Distance(Centroids[idx_cent], Current);
				if (distance < Min_Distance){
					Cluster_Idx = idx_cent;
					Min_Distance = distance; 
				} // ------- then find the minimum distance and minimum centroid 
			}
			// =============== Find the new centroids ==========
			// OK, now we found the centroid index ==> Sum! 
			// Cluster_Sum[Cluster_Idx] + Current;
			vector<float> Sum_Cluster = Cluster_Sum[Cluster_Idx];
			Count_Cluster[Cluster_Idx] ++; 
			for (int idx_sum = 0; idx_sum < Sum_Cluster.size(); idx_sum++){
				Cluster_Sum[Cluster_Idx][idx_sum] = Sum_Cluster[idx_sum] + Current[idx_sum];
			}
		}
		// =================== Update Centriod =====================
		for (int temp = 0; temp < K; temp ++){
			for (int temp_eachpt = 0; temp_eachpt < Cluster_Sum[temp].size(); temp_eachpt++){
				Centroids[temp][temp_eachpt] = Cluster_Sum[temp][temp_eachpt] / Count_Cluster[temp];
			}
		}
		// ========= Initialize ClusterSum ========
		Cluster_Sum.clear();
		Cluster_Sum = Centroids;
		
		// ==================== Print ===============
		// cout << num_iter << endl; 
		// for (int temp = 0; temp < K; temp ++){
		// 	for (int temp_elem = 0; temp_elem < Centroids[temp].size(); temp_elem++){
		// 		cout << Centroids[temp][temp_elem] << " ";
		// 	}
		// 	cout << endl;
		// }
		// cout << endl; 
		num_iter ++; 

		// ================= Check Stopping Condition ================
		float distance_sum = 0.0; 
		for (int temp = 0; temp < K; temp ++){
			distance_sum += Distance(Centroids[temp], Prev_Centroids[temp]);
		}

		if (distance_sum < 0.000001 || num_iter > 100){
			break;
		}
	}
	return Centroids;
}

vector< vector< vector<float> > > Sample_Covariance_for_each_Cluster(vector< vector< vector<float> > > Cluster_Map,  
	vector< vector<float> > Initial_Sample_Means, int K, int dim, int N){

	vector< vector< vector<float> > > COV_MATRIX; 

	for (int idx = 0; idx < K; idx ++){
		COV_MATRIX.push_back(Initial_Sample_Means);
	}
	
	for (int idx = 0 ; idx < K; idx ++){
		COV_MATRIX[idx].clear();
	}
	// cout << COV_MATRIX.size() << endl;
	// Initialize 
	vector< vector<float> > Sum_Cov; 
	vector<float> Temp_for_Cov; 
	for (int idx = 0; idx < dim; idx ++ ){
		for (int idx2 = 0; idx2 < dim; idx2++){
			Temp_for_Cov.push_back(0);
		}
		Sum_Cov.push_back(Temp_for_Cov);
		Temp_for_Cov.clear();
	}

	// for (int temp = 0; temp < Sum_Cov.size(); temp ++){
	// 	for (int temp2 = 0; temp2 < Sum_Cov[temp].size(); temp2++){
	// 		cout << Sum_Cov[temp][temp2] << " ";
	// 	}
	// 	cout << endl; 
	// }
	// cout << " ============================== " << endl; 

	for (int Cluster_Idx = 0; Cluster_Idx < K; Cluster_Idx++ ){
		vector< vector<float> > Cluster_Data = Cluster_Map[Cluster_Idx];
		vector<float> Cluster_Mean = Initial_Sample_Means[Cluster_Idx];

		for (int each_point_idx = 0; each_point_idx < Cluster_Data.size(); each_point_idx ++){
			vector<float> Each_Point = Cluster_Data[each_point_idx];
			vector<float> Difference; 
			vector<float> Temp_for_Matrix; 
			vector< vector<float> > Difference_Matrix; 

			// Then measure the variance distance 
				// Each_Point - Initial_Sample_Means[Cluster_Idx]
			for (int each_elem = 0; each_elem < dim; each_elem ++){
				Difference.push_back( Cluster_Mean[each_elem] - Each_Point[each_elem]) ;
			}


			for (int each_elem = 0; each_elem < Each_Point.size(); each_elem ++){
				for (int each_elem2 = 0; each_elem2 < Each_Point.size(); each_elem2 ++){
					Temp_for_Matrix.push_back(Difference[each_elem] * Difference[each_elem2]);
				}
				Difference_Matrix.push_back(Temp_for_Matrix);
				Temp_for_Matrix.clear();
			}

			for (int temp = 0; temp < dim; temp ++){
				for (int temp2 = 0; temp2 < dim; temp2 ++){
					Sum_Cov[temp][temp2] += Difference_Matrix[temp][temp2];
				}
			}
			// Clear
			Difference.clear();
			Temp_for_Cov.clear();
			Difference_Matrix.clear();
		}

		// 1/(n-1)
		for (int temp = 0; temp < dim; temp ++){
			for (int temp2 = 0; temp2 < dim; temp2 ++){
				Sum_Cov[temp][temp2] = (Sum_Cov[temp][temp2] / (N-1));
			}
		}

		// COV_MATRIX.push_back(Sum_Cov);
		COV_MATRIX[Cluster_Idx] = Sum_Cov;
	}

	return COV_MATRIX;

}

float Normal_Distribution_Density(vector<float> data_point, vector<float> mean,  vector< vector<float> > Cov_Matrix){
	// Function Declarazation 
	float Abs_Determinant(vector< vector<float> >);
	float Power_Compute(float, int);
	vector< vector<float> > Inverse_Matrix( vector< vector<float> >);
	vector <vector <float> > Matrix_Multiplication(vector< vector<float> >, vector< vector<float> >);
	vector< vector<float> > Transpose(vector< vector<float> >);


	float dim = data_point.size();
	// float part1 = pow((1/float(sqrt(2*PI)), dim));
	float part1 = sqrt(2*PI);
	part1 = 1/float(part1);
	part1 = Power_Compute(part1, dim);
	float part2 = sqrt(1/ float(Abs_Determinant(Cov_Matrix)));
	if (Abs_Determinant(Cov_Matrix) == 0){
		cout << "HAHAHAHAHAHAHA" << endl;
	}
	float density; 

	vector<float> centered_data; 
	for (int idx = 0; idx < dim; idx++){
		float val = (data_point[idx] - mean[idx]);
		centered_data.push_back(val);
	}
	vector< vector<float> > Inverse_Covariance = Inverse_Matrix(Cov_Matrix);
	vector< vector<float> > Centered_data_Matrix; 
	Centered_data_Matrix.push_back(centered_data);
	vector< vector<float> > Transposed_Centered_data = Transpose(Centered_data_Matrix);
	vector< vector<float> > part3 = Matrix_Multiplication(Centered_data_Matrix, Inverse_Covariance);
	vector< vector<float> > part4;
	part4 = Matrix_Multiplication(part3, Transposed_Centered_data);
	float part5 = -0.5 * part4[0][0];
	part5 = exp(part5);
	density = part1 * part2 * part5;

	return density;

}

float Log_Normal_Density(vector<float> data_point, vector<float> mean,  vector< vector<float> > Cov_Matrix){
    float Abs_Determinant(vector< vector<float> >);
    float Power_Compute(float, int);
    vector< vector<float> > Inverse_Matrix( vector< vector<float> >);
    vector <vector <float> > Matrix_Multiplication(vector< vector<float> >, vector< vector<float> >);
    vector< vector<float> > Transpose(vector< vector<float> >);

    float dim = data_point.size();
    
    float part1 = -(dim / float(2)) * log (2*PI);
    float part2 = -0.5 * log(Abs_Determinant(Cov_Matrix));
    part1 = part1 + part2; // Normalized Condition 
    
    
    float density; 
    vector<float> centered_data; 
    for (int idx = 0; idx < dim; idx++){
        float val = (data_point[idx] - mean[idx]);
        centered_data.push_back(val);
    }
    vector< vector<float> > Inverse_Covariance = Inverse_Matrix(Cov_Matrix);
    vector< vector<float> > Centered_data_Matrix; 
    Centered_data_Matrix.push_back(centered_data);
    vector< vector<float> > Transposed_Centered_data = Transpose(Centered_data_Matrix);
    vector< vector<float> > part3 = Matrix_Multiplication(Centered_data_Matrix, Inverse_Covariance);
    vector< vector<float> > part4;
    part4 = Matrix_Multiplication(part3, Transposed_Centered_data);
    float part5 = -0.5 * part4[0][0];
    // part5 = exp(part5);

    density = part1 + part5;

    return density;




}

float Gaussian_Mixture_Density(vector<float> data_point, vector< vector<float> > Mean_Set, 
	vector< vector< vector<float> > > COV_Set, vector<float> Weight_Set  ){

	float Normal_Distribution_Density(vector<float>, vector<float>, vector< vector<float> > );
	int dim = data_point.size();
	int K = Mean_Set.size();

	float density_sum = 0.0; 
	float density; 
	float very_small = 0.0000001;

	for (int idx = 0; idx < K; idx ++){
		density = Normal_Distribution_Density(data_point, Mean_Set[idx], COV_Set[idx]);
		// if (density < very_small){
		// 	density = very_small;
		// }
		density = density * Weight_Set[idx];
		density_sum += density;
	}
	if (density_sum > 0) {
		return density_sum;
	}
	else{
		return very_small;
	}

}

float log_GMM(vector<float> Data_X, vector< vector<float> > Mean_Set, vector< vector< vector<float> > > COV_Set, vector<float> Weight_Set  ){

	float Log_Normal_Density(vector<float> , vector<float> ,  vector< vector<float> > );
	vector<float> each_data_per_cluster;
	float sum_density = 0;
	float each_density; 
	
	for (int cluster_idx = 0; cluster_idx < Mean_Set.size(); cluster_idx ++){
		float log_density_per_cluster = Log_Normal_Density(Data_X, Mean_Set[cluster_idx], COV_Set[cluster_idx]);
		log_density_per_cluster += log(Weight_Set[cluster_idx]);
		each_data_per_cluster.push_back(log_density_per_cluster);
	}

	for (int cluster_idx = 0; cluster_idx < Mean_Set.size(); cluster_idx ++){
		each_density = each_data_per_cluster[cluster_idx];
		sum_density += each_density;
	}
	

	return sum_density;

}


float Score_likelihood(vector< vector<float> > Data_Matrix, vector< vector<float> > Mean_Set, vector< vector< vector<float> > > COV_Set, vector<float> Weight_Set    ){
	float Gaussian_Mixture_Density(vector<float> , vector< vector<float> > , vector< vector< vector<float> > > , vector<float>   );

	float score_density_sum = 0; 
	float score_density; 

	for (int data_idx = 0; data_idx < Data_Matrix.size(); data_idx ++){
		score_density = Gaussian_Mixture_Density(Data_Matrix[data_idx], Mean_Set, COV_Set, Weight_Set);
		score_density_sum += score_density;
	}
	return score_density_sum;


}



// ======================== Matrix Operation ============================

void Print_Matrix(vector< vector<float> > MATRIX){
	for (int row = 0; row < MATRIX.size(); row++){
		for (int col = 0; col < MATRIX[row].size(); col++){
			cout << MATRIX[row][col] << " | ";
		}
		cout << endl;
	}
}

void Print_Vector(vector<float> A){
	for (int idx = 0; idx < A.size(); idx++){
		cout << A[idx] << " ";
	}
}

void Print_COV(vector< vector< vector<float> > > A){
	for (int mat_idx = 0; mat_idx < A.size(); mat_idx ++){
		cout << mat_idx+1 << " Cluster" << endl;
		vector< vector<float> > MATRIX = A[mat_idx];
		for (int row = 0; row < MATRIX.size(); row++){
			for (int col = 0; col < MATRIX[row].size(); col++){
				cout << MATRIX[row][col] << " ";
			}
			cout << endl;
		}
		cout << "---------------" << endl; 
	}
}

vector < vector<float> > eye (int n){
    vector< vector<float> > b;
    vector<float> hans; 
    
    for (int idx = 0; idx < n; idx++){
        for (int idx2 = 0; idx2 < n; idx2++){
            hans.push_back(0);
        }
        b.push_back(hans);
        hans.clear();
    }

    for (int row = 0; row < n; row++){
        for (int col = 0; col < n; col ++){
            if (row == col){
                b[row][col] = 1;
            }
            else{
                b[row][col] = 0;
            }
        }
    }

    return b;
}

vector< vector<float> > Inverse_Matrix( vector< vector<float> > a) {
    int n=0,m=0,i=0,j=0,p=0,q=0;
    float temp, tem, temp1, temp2, temp4, temp5; 
    n = a.size();
    vector< vector<float> > b = eye(n);

    for(i=0;i<n;i++) // n == Matrix Order 
    {
        temp=a[i][i];
        if(temp<0)
        temp=temp*(-1);
        p=i;
        for(j=i+1;j<n;j++)
        {
            if(a[j][i]<0)
                tem=a[j][i]*(-1);
            else
                tem=a[j][i];
            if(temp<0)
                temp=temp*(-1);
            if(tem>temp)
            {
                p=j;
                temp=a[j][i];
            }
        }

        //row exchange in both the matrix
        for(j=0;j<n;j++)
        {
            temp1=a[i][j];

            a[i][j]=a[p][j];

            a[p][j]=temp1;

            temp2=b[i][j];
            b[i][j]=b[p][j];
            b[p][j]=temp2;
        }
        //dividing the row by a[i][i]

        temp4=a[i][i];
        for(j=0;j<n;j++)
        {
            a[i][j]=(float)a[i][j]/temp4;
            b[i][j]=(float)b[i][j]/temp4;
        }
        //making other elements 0 in order to make the matrix a[][] an indentity matrix and obtaining a inverse b[][] matrix
        for(q=0;q<n;q++)
        {
            if(q==i)
                continue;
            temp5=a[q][i];
            for(j=0;j<n;j++)
            {
                a[q][j]=a[q][j]-(temp5*a[i][j]);
                b[q][j]=b[q][j]-(temp5*b[i][j]);
            }
        }
    }

    return b;

}

vector< vector<float> > Transpose(vector< vector<float> > A){
	// A = 1 * 4;
	vector< vector<float> > transposed_matrix; 
	vector<float> temp; 
	int row_size = A.size(); // = 1
	int col_size = A[0].size(); // = 4 

	for (int row = 0; row < col_size; row ++){
		for (int col = 0; col < row_size; col ++){
			// col will be 0
			temp.push_back(A[col][row]);
		}
		transposed_matrix.push_back(temp);
		temp.clear();
	}
	return transposed_matrix;
}

vector< vector<float> > Vector_to_Matrix(vector<float> A){
	vector< vector<float> > result; 
	vector<float> temp_result; 
	for (int idx = 0; idx < A.size(); idx++){
		temp_result.push_back(A[idx]);
		result.push_back(temp_result);
		temp_result.clear();
	}
	return result;
}

vector <vector <float> > Matrix_Multiplication(vector< vector<float> > A, vector< vector<float> > B){
    int N = A.size();
    int M = B[0].size();
    int common = A[0].size();
    vector< vector<float> > result;
    vector<float> temp; 

    // ============ Initializaiton ===============
    for (int row = 0; row < N; row ++){
        for (int col = 0; col < M; col ++){
            temp.push_back(0);
        }
        result.push_back(temp);
        temp.clear();
    }
    for (int row = 0; row < N; row ++){
        for (int col = 0; col < M; col ++){
            for (int inner = 0; inner < common; inner ++){
                result[row][col] += A[row][inner] * B[inner][col];
            }
        }
    }

    return result;
}

vector <vector<float> > Zeros(int K){
    vector<float> row_temp;
    vector< vector<float> > Zero_Matrix;

    for (int row = 0; row < K; row++){
        for (int col = 0; col < K; col ++){
            row_temp.push_back(0);
        }
        Zero_Matrix.push_back(row_temp);
        row_temp.clear();
    }
    return Zero_Matrix;
}

vector< vector<float> > pivotize(vector< vector<float> > m){
    int n;
    n = m.size();
    vector< vector<float> > ID= eye(n);
    vector<float> temp; 

    for (int idx = 0; idx < n; idx ++){
        int row = -100000000; 
        int max_idx = 0;
        for (int idx2 = idx; idx2 < n; idx2++){
            if (row < abs(m[idx2][idx])){
                row = abs(m[idx2][idx]);
                max_idx = idx2;
            }
        }


        if (idx != max_idx){
            temp = ID[idx];
            ID[idx] = ID[max_idx];
            ID[max_idx] = temp;
        }
    }
    return ID; 
}

vector <vector<float> > Input_Matrix(int n){
    vector< vector<float> > b;
    vector<float> hans; 
    float input;
    
    for (int idx = 0; idx < n; idx++){
        for (int idx2 = 0; idx2 < n; idx2++){
            cout << idx+1 << "," << idx2+1 << " : ";
            cin >> input; 
            hans.push_back(input);
        }
        b.push_back(hans);
        hans.clear();
    }
    return b;
}

float Abs_Determinant(vector< vector<float> > A){
    vector< vector< vector<float> > > LU_Decompose( vector< vector<float> >  );
    vector< vector< vector<float> > > LUP = LU_Decompose(A); 
    vector< vector<float> > L = LUP[0];
    vector< vector<float> > U = LUP[1];
    vector< vector<float> > P = LUP[2];    

    float L_Prod = 1;
    float U_Prod = 1; 
    float result; 

    for (int idx = 0; idx < L.size(); idx ++){
        L_Prod = L_Prod * L[idx][idx];
    }
    cout << endl; 
    for (int idx = 0; idx < U.size(); idx ++){
        U_Prod = U_Prod * U[idx][idx];
    }

    result = L_Prod * U_Prod;
    if (result < 0){
        result = result * (-1);
    }

    return result;
}

vector< vector< vector<float> > > LU_Decompose( vector< vector<float> > A ){
    int n = A.size();
    vector< vector< vector<float> > > LUP; 

    vector< vector<float> > L = Zeros(n);
    vector< vector<float> > U = Zeros(n);
    vector< vector<float> > P = pivotize(A);
    vector< vector<float> > A2 = Matrix_Multiplication(P,A);

    for (int j = 0; j < n; j++){
        L[j][j] = 1.0;
        for (int i = 0; i < j+1; i++){
            float s1 = 0; 
            for (int k = 0; k < i; k++){
                s1 += U[k][j] * L[i][k];
            }
            U[i][j] = A2[i][j] - s1;
        }
        for (int i = j; i < n; i ++){
            float s2 = 0; 
            for (int k = 0; k < j; k++){
                s2 += U[k][j] * L[i][k];
            }
            L[i][j] = (A2[i][j] - s2) / U[j][j];
        }
    }
    LUP.push_back(L);
    LUP.push_back(U);
    LUP.push_back(P);

    return LUP; 
}




// ---------------------------- Program  ------------------------------------------
int main () {
	/* ============== DATA LOAD AND GENERATING MATRIX =============== */ 
	string Train_Blue = "Wb_train.txt"; 
	string Train_Red= "Wr_train.txt";
	string Train_Green = "Wg_train.txt";

	string Test_Blue = "Wb_test.txt"; 
	string Test_Red= "Wr_test.txt";
	string Test_Green = "Wg_test.txt";
	

	vector< vector<float> > Train_Blue_Matrix = LoadFile(Train_Blue);
	vector< vector<float> > Train_Red_Matrix = LoadFile(Train_Red);
	vector< vector<float> > Train_Green_Matrix = LoadFile(Train_Green);
	vector< vector<float> > Test_Blue_Matrix = LoadFile(Test_Blue);
	vector< vector<float> > Test_Red_Matrix = LoadFile(Test_Red);
	vector< vector<float> > Test_Green_Matrix = LoadFile(Test_Green);

	vector< vector<float> > Train_Matrix; 
	vector< vector<float> > Test_Matrix;

	Train_Matrix.insert(Train_Matrix.end(), Train_Blue_Matrix.begin(), Train_Blue_Matrix.end());
	Train_Matrix.insert(Train_Matrix.end(), Train_Red_Matrix.begin(), Train_Red_Matrix.end());
	Train_Matrix.insert(Train_Matrix.end(), Train_Green_Matrix.begin(), Train_Green_Matrix.end());

	Test_Matrix.insert(Test_Matrix.end(), Test_Blue_Matrix.begin(), Test_Blue_Matrix.end());
	Test_Matrix.insert(Test_Matrix.end(), Test_Red_Matrix.begin(), Test_Red_Matrix.end());
	Test_Matrix.insert(Test_Matrix.end(), Test_Green_Matrix.begin(), Test_Green_Matrix.end());

	vector< vector<float> > My_Answer_Red; // index = 0
	vector< vector<float> > My_Answer_Green; // index = 1 
	vector< vector<float> > My_Answer_Blue; // index = 2


	/* ============== Compute Initial Sample Mean and Variance and Weight =============== */ 

	int K;
	int dim = 4; 
	int N = Train_Matrix.size();
	string color; 
	vector< vector<float> > Train_Data_Set;
	cout << "Which Color You want? (R,G,B)" << endl;
	cin >> color; 
	if (color == "B"){
		Train_Data_Set = Train_Blue_Matrix;
	}
	else if (color == "G"){
		Train_Data_Set = Train_Green_Matrix;	
	}

	else if (color == "R"){
		Train_Data_Set = Train_Red_Matrix;		
	}
	cout << "How many clusters you want?" << endl;
	cin >> K;  

	 

	vector< vector<float> > Initial_Sample_Means = Centroid_for_Mean(Train_Data_Set ,K);
	vector< vector< vector<float> > > Cluster_Map = K_Means_Clustering(Train_Data_Set, Initial_Sample_Means,K);

	

	vector< vector< vector<float> > > COV_MATRIX = Sample_Covariance_for_each_Cluster(Cluster_Map, Initial_Sample_Means, K, dim, N);

	
	

	vector<float> Weight; 
	for (int idx = 0; idx < K; idx ++){
		Weight.push_back(1/float(K));
	}



	int num_iter = 0; 




	
	// ==================== EM Algorithm ==================

	// 재료준비 
	vector< vector<float> > Sample_Means = Initial_Sample_Means;
	vector<float> Initial_Weight =  Weight;
	vector< vector< vector<float> > > Initial_COV_MATRIX = COV_MATRIX;
	vector< vector<float> > Cluster_Prob;
	

	// Initialize Cluster Probability MAtrix 
	for(int idx = 0; idx < Train_Data_Set.size(); idx ++){
		Cluster_Prob.push_back(Weight);
	}
	for (int idx = 0; idx < Train_Data_Set.size(); idx ++){
		Cluster_Prob[idx].clear();
	}

	vector<float> Cluster_Number; 
	

	// Initialize Cluster_Number Matrix 
	for (int idx = 0; idx < Cluster_Map.size(); idx ++){
		float each_cluster_element_number = Cluster_Map[idx].size();
		Cluster_Number.push_back(each_cluster_element_number);
	}


	float prev_likelihood = Score_likelihood(Train_Data_Set, Sample_Means, COV_MATRIX, Weight);
	float current_likelihood; 

	vector<float>sumsum; 
	for (int idx = 0; idx < Train_Data_Set[0].size(); idx++){
		sumsum.push_back(0);
	}
	for (int data_idx = 0; data_idx < Train_Data_Set.size(); data_idx ++){
		for (int each_idx = 0; each_idx < Train_Data_Set[0].size(); each_idx ++){
			sumsum[each_idx] += Train_Data_Set[data_idx][each_idx];
		}
	}

	Print_Vector(sumsum);


	while (true){
		//  ======================= E_STEP =========================
		current_likelihood = prev_likelihood;
		for (int data_idx = 0; data_idx < Train_Data_Set.size(); data_idx++ ){
			// Compute P(zj | Xi) = Xi 가 j 클러스터에 속할 확률 
			// Compute Cluster likelihood for each data points and for each clusters 
			vector<float> Data_X = Train_Data_Set[data_idx];
			vector<float> cluster_prob_for_each_data;

			// float log_GMM_density = log_GMM(Data_X, Sample_Means, COV_MATRIX, Weight);
			float GMM_density = Gaussian_Mixture_Density(Data_X, Sample_Means, COV_MATRIX, Weight);
			float Each_Data_Density_Per_Cluster;
			
			
			for (int cluster_idx = 0; cluster_idx < K; cluster_idx++){
				// for each cluster
				// float just_density = Normal_Distribution_Density(Data_X, Sample_Means[cluster_idx], COV_MATRIX[cluster_idx]);
				Each_Data_Density_Per_Cluster = Normal_Distribution_Density(Data_X, Sample_Means[cluster_idx], COV_MATRIX[cluster_idx]);
				Each_Data_Density_Per_Cluster *= Weight[cluster_idx];
				
				float Weight_Density = Each_Data_Density_Per_Cluster / GMM_density; // P(zj = 1 | X)
				// cout << cluster_idx << " | " << Each_Data_Density_Per_Cluster << " | " << Weight_Density;
				cluster_prob_for_each_data.push_back(Weight_Density);
				
			}
			// Print_Vector(cluster_prob_for_each_data);
			Cluster_Prob[data_idx] = (cluster_prob_for_each_data);
			cluster_prob_for_each_data.clear();
		}
		// Complete to make Cluster Probability Matrix 
		// each row represents data point 
		// each column represents each clusters 

		// M_STEP 
		for (int cluster_idx = 0; cluster_idx < K; cluster_idx++) 
		{
			float each_cluster_number = 0.0; // For counting how many data points belong to the each cluster 
			vector<float> each_cluster_mean = Sample_Means[cluster_idx]; // Previous Sample Mean 
			vector< vector<float> > each_cluster_cov = COV_MATRIX[cluster_idx]; // Previous Sample Covariance 
			vector< vector<float> > Sample_COV_MATRIX; // Compute Current Sample Covariance 
			vector< vector<float> > accumulate_Sample_COV_MATRIX; // Compute Current Sample Covariance 
			vector<float> accumulated_each_data;
			accumulated_each_data.clear();
			accumulate_Sample_COV_MATRIX.clear();
			// Print_Matrix(each_cluster_cov);
			
			for (int temp_idx = 0; temp_idx < Train_Data_Set[0].size(); temp_idx ++){
				accumulated_each_data.push_back(0);
			} // Initializing the current sample mean 

			vector<float> temp_sample_cov;
			for (int temp_idx1 = 0; temp_idx1 < each_cluster_cov.size(); temp_idx1++){
					for (int temp_idx2 = 0; temp_idx2 < each_cluster_cov[0].size(); temp_idx2++){
						temp_sample_cov.push_back(0);
					}
					accumulate_Sample_COV_MATRIX.push_back(temp_sample_cov);
					temp_sample_cov.clear();
				} // Initialize the current sample covariance
			// Print_Matrix(accumulate_Sample_COV_MATRIX);

			
			for (int data_idx = 0; data_idx < Train_Data_Set.size(); data_idx++){
				// for each data point, initialize the variable 
				vector<float> each_data = Train_Data_Set[data_idx]; // each data 
				vector<float> Centered_for_each_data ; // centered data to the sample mean 

				// the probability that the data point belonging to the cluster 
				float each_cluster_prob = Cluster_Prob[data_idx][cluster_idx];
				
				// Compute (P(z|X) * X)
				for (int temp_idx = 0; temp_idx < Train_Data_Set[0].size(); temp_idx++){
					// cout << each_cluster_prob << endl;
					each_data[temp_idx] = each_cluster_prob * each_data[temp_idx];
				}
				
				// Compute Uj 
				for(int update_idx = 0; update_idx < Train_Data_Set[0].size(); update_idx++){
					accumulated_each_data[update_idx] += each_data[update_idx];
				}

				// Compute Nj 
				each_cluster_number += each_cluster_prob; 
				// cout << each_cluster_number << endl;
				
			}
			// Complete to iterate for all data points. 

			// =========== Compute Nj, Uj, COVj, Wj ===========
			// cout <<each_cluster_number << endl;
			Cluster_Number[cluster_idx] = each_cluster_number; // Nj 
			// cout << endl;
			// cout << cluster_idx << " | ";
			// Print_Vector(accumulated_each_data); // For computing uj 
			// cout << endl; 
			// cout << each_cluster_number << " | ";
			// cout << each_cluster_number << endl;
			


			// ============ Uj ===============
			for (int each_point_idx = 0; each_point_idx < Train_Data_Set[0].size(); each_point_idx++){
				float each_mu = accumulated_each_data[each_point_idx] / each_cluster_number;
				accumulated_each_data[each_point_idx] = each_mu;
			}


			// ============ COVj ===============
			for (int temp_idx1 = 0; temp_idx1<Sample_COV_MATRIX.size(); temp_idx1++){
					for (int temp_idx2 = 0; temp_idx2 < Sample_COV_MATRIX[0].size(); temp_idx2++){
						accumulate_Sample_COV_MATRIX[temp_idx1][temp_idx2] = accumulate_Sample_COV_MATRIX[temp_idx1][temp_idx2] / each_cluster_number;
					}
				}
			Sample_Means[cluster_idx] = accumulated_each_data; // Update Uj
			vector<float> Centered;
			vector< vector<float> > Each_Sample_Cov;
			vector< vector<float> > Sum_Each_Sample_Cov = eye(Train_Data_Set[0].size()); 

			for (int data_idx = 0; data_idx < Train_Data_Set.size(); data_idx++){
				vector<float> Data_X = Train_Data_Set[data_idx];
				// Centering 
				Centered = Data_X; 
				for (int each_idx = 0; each_idx < Data_X.size(); each_idx++){
					Centered[each_idx] = Data_X[each_idx] - accumulated_each_data[each_idx];
				}
				Each_Sample_Cov = Matrix_Multiplication( Vector_to_Matrix(Centered), Transpose(Vector_to_Matrix(Centered)));
				for (int row_temp = 0; row_temp < Each_Sample_Cov.size(); row_temp ++){
					for (int col_temp = 0; col_temp < Each_Sample_Cov[0].size(); col_temp++){
						Each_Sample_Cov[row_temp][col_temp] *= Cluster_Prob[data_idx][cluster_idx];
					}
				}

				for (int row_temp = 0; row_temp < Each_Sample_Cov.size(); row_temp ++){
					for (int col_temp = 0; col_temp < Each_Sample_Cov[0].size(); col_temp++){
						Sum_Each_Sample_Cov[row_temp][col_temp] += Each_Sample_Cov[row_temp][col_temp];
					}
				}
			}

			for (int row_temp = 0; row_temp < Each_Sample_Cov.size(); row_temp ++){
					for (int col_temp = 0; col_temp < Each_Sample_Cov[0].size(); col_temp++){
						Sum_Each_Sample_Cov[row_temp][col_temp] /= Cluster_Number[cluster_idx];
					}
				}

			// Print_Matrix(Sample_Means);
			// cout << endl;
			// COV_MATRIX[cluster_idx] = accumulate_Sample_COV_MATRIX; // update COVj
			COV_MATRIX[cluster_idx] = Sum_Each_Sample_Cov; 
			Weight[cluster_idx] = each_cluster_number / float(Train_Data_Set.size());
		}
		// STOP CONDITION 
		current_likelihood = Score_likelihood(Train_Data_Set, Sample_Means, COV_MATRIX, Weight);
		if (abs(current_likelihood - prev_likelihood) < 0.00001 || num_iter > 10){
			cout << prev_likelihood << " " << current_likelihood << endl;
			break;
		}
		num_iter ++;
		// cout << current_likelihood << endl;

	}

		


		cout << "============================ Initial_Parameters ===================" << endl; 
		cout << " ------ Initial Sample Weights ---------" << endl;
		Print_Vector(Initial_Weight);
		cout << endl; cout << endl;
		cout << " ------ Initial Sample Means ---------" << endl;
		Print_Matrix(Initial_Sample_Means);
		cout << endl; 
		cout << " ------ Initial Sample COV ---------" << endl;
		Print_COV(Initial_COV_MATRIX);

		cout << endl;cout << endl; cout << endl;

		cout << "============================ After_EM_Parameters ===================" << endl; 
		cout << "numIter : " << num_iter << endl;
		cout << " ------ EM Sample Weights ---------" << endl;
		Print_Vector(Weight);
		cout << endl; cout << endl;
		cout << " ------ EM Sample Means ---------" << endl;
		Print_Matrix(Sample_Means);
		cout << endl; 
		cout << " ------ EM Sample COV ---------" << endl;
		Print_COV(COV_MATRIX);

		// cout << endl;
		// Print_Vector(Cluster_Number);
		// cout << endl; 
		// Print_Matrix(Cluster_Prob);
		// cout << endl;
		// cout << 3.62109e-23 + 2 << endl;

		// cout << Abs_Determinant(Inverse_Matrix (COV_MATRIX[0])) << endl;
		// Print_Matrix( Inverse_Matrix (COV_MATRIX[1]));
		// cout << Abs_Determinant(Inverse_Matrix (COV_MATRIX[1])) << endl;

		// for (int idx = 0; idx < Train_Data_Set.size(); idx++){
		// 	cout << Gaussian_Mixture_Density(Train_Data_Set[idx], Sample_Means, COV_MATRIX, Weight)  << endl;
		// }

		// for (int idx = 0; idx < Train_Data_Set.size(); idx++){
		// 	for (int cluster_idx = 0; cluster_idx < K; cluster_idx ++){
		// 	cout << Log_Normal_Density(Train_Data_Set[idx], Sample_Means[cluster_idx], COV_MATRIX[cluster_idx]) << " ";
		// 	}
		// 	cout << endl;
		// }

		
	 	// Print_Vector(log_GMM(Train_Data_Set, Sample_Means, COV_MATRIX, Weight));





	cout << endl; 
	return 0;
}