

#include "stdafx.h"
#include "markov.h"


///////////////////////////////////Markov  constructor/destructor///////////////////////////////////////////
Markov::Markov()
{

}
Markov::Markov(const wchar_t *fname)
{
	
	//int classnum = get_classnum();
    FILE *fp = _wfopen(fname, L"rt");
    if (fp != 0) {
		if (fwscanf(fp, L"%d", &classnum) != 1){
            fclose(fp);
			return;
		}
		if (fwscanf(fp, L"%d", &som_size) != 1) {
            fclose(fp);
            return;
        }
	}
		
}

Markov::~Markov()
{
	// De-Allocate memory to prevent memory leak
	for (int i = 0; i < som_size; ++i){
		delete [] transitionMat[i];		
	}
    delete [] transitionMat;
	delete [] firstStateMat;
	delete [] bmuClassIndexArray;
	
	// De-Allocate memory to prevent memory leak
    for (int i = 0; i < dirsize; ++i)
	  delete [] directionMat[i];
    delete [] directionMat;
	delete [] firstAngleMat;
	delete [] angleClassIndexArray;
	
	delete [] som_dist_map;
	delete [] int_data_carrier;
}
///////////////////////////////////Markov  constructor/destructor///////////////////////////////////////////
//--////////////////////////////////////////////////////////////
void Markov::evaluate(SOM *som, PREC rec, const wchar_t *som_map_name)
{
	int trained_classes = get_classnum();
	int somsize = get_somsize();
	wprintf(L"Number of Trained Classes: %d\n", trained_classes);
	wprintf(L"--------------------------------\n");
	
	//for number of classes	
	int trans1 = getTransitionVector(som, rec); // fills up bmuClassIndexArray with pattern's bmu indeces
	int trans2 = getTransitionAngle(rec); // fills up angleClassIndexArray with pattern's angle categories
	
	//wprintf(L"t1:%d t2:%d\n", trans1, trans2);
	float mean_tProb =  0.0;
	float mean_aProb =  0.0;
	float mean_dProb =  0.0;
	float* maxOutValue =  new float[trained_classes];
	
	FILE *sdm = _wfopen(som_map_name, L"rt");
	som_dist_map = readMat(sdm); 
	
	float maxTProb =0.0;
	int maxTClass = 0;

	for (int c = 0; c < trained_classes; c++) {
		//Vres probability gia kathe class [trans]
		mean_tProb = getTransProb(som, rec, bmuClassIndexArray, trans1, c+1, somsize) /trans1;
		//Vres probability gia kathe class [direction]
		mean_aProb = getAngleProb(rec, angleClassIndexArray, trans2, c+1) /trans2;
		//Vres logo ML/L gia kathe class [levenshtein distance]
		mean_dProb = getLevenshteinProb(bmuClassIndexArray, trans1, c+1) ;

		//Return Max of Psom*Pof*Lev
		maxOutValue[c] = mean_tProb * mean_aProb * mean_dProb;
		if (c > 0){
			if (maxOutValue[c] > maxTProb){
				maxTProb = maxOutValue[c];
				maxTClass = c+1;
			}
			
		}
		wprintf(L"Class %d strength: %f \n", c+1, maxOutValue[c]);
		wprintf(L"*********************************\n");

	}
	wprintf(L"\n** Winning Class = %d with value %1.4f **\n\n", maxTClass, maxTProb);

}
/////////////////////
float Markov::getTransProb(SOM *som, PREC classrec, int *bmuArray, int transitions, int classid, int somsize)
{
	float cProb = 0;
	wchar_t filename[_MAX_PATH] = L"";
	wchar_t firstState_name[_MAX_PATH] = L"";
	wchar_t str[_MAX_PATH] = L"";

	wcscat(filename, L"mm_results/mm_trans_matrix_class_");
	wcscat(firstState_name, L"mm_results/mm_firststate_vector_class_");
	_itow(classid, str, 10);
	wcscat(filename, str);
	wcscat(firstState_name, str);

	FILE *ff = _wfopen(firstState_name, L"rt");
	wprintf(L"Reading /%s...\n", firstState_name);
	//int vec_size2 = readMat(ff);
	float *firstStateMat= readMat(ff);
	//float *firstStateMat = new float[vec_size2];
	//firstStateMat = fdata_carrier;

	FILE *fp = _wfopen(filename, L"rt");
	
	//int vec_size = readMat(fp);
	float *MMtransMat= readMat(fp);
	wprintf(L"Reading /%s ...\n", filename);
	
	
	// START PROCESSING
	int index = 1;
	int win_node;
	int sourcenode =0;
	int dim = somsize;
	int mat_pos =1;
	
	// Starting with first State
	float dist;
	float max_value = 0, new_value;
	int requested_node = bmuArray[0];
	// Get distances for each som node from requested node 
	// and multiply to their firstState probability
	for (int n = 0; n < somsize; n++){
		
		dist = som_dist_map[(requested_node-1)*dim + n];
		new_value = (1-dist) * firstStateMat[n];
		if (new_value > max_value) {
            max_value = new_value;
			index = n+1; //cause it starts from zero			
        }
	}
	win_node = index;
	//if (max_value>0.01)
		//wprintf(L"First Max index: %d  with value: %f\n", win_node, max_value);
	cProb = max_value;	
	
	//Next nodes
	index = 1;
	for (int pos = 1; pos < transitions; pos++){
		requested_node = bmuArray[pos];
		for (int n = 0; n < somsize; n++){
			
			dist = som_dist_map[(requested_node-1)*dim + n];
			new_value = (1-dist) * MMtransMat[(win_node-1)*dim + n];
			if (new_value > max_value) {
				max_value = new_value;
				index = n+1; //cause it starts from zero			
			}
		}
		win_node = index;
		//if (max_value>0.01)
			//wprintf(L"Max index: %d  with value: %f\n", win_node, max_value);

		cProb += max_value;	
	}
	// END OF PROCESSING
	
	return cProb;
}
////////////////////
float Markov::getAngleProb(PREC classrec, int *angleArray, int transitions, int classid)
{
	float cProb= 1.0;
	wchar_t filename[_MAX_PATH] = L"";
	wchar_t firstState_name[_MAX_PATH] = L"";
	wchar_t str[_MAX_PATH] = L"";

	wcscat(filename, L"mm_results/mm_dir_matrix_class_");
	wcscat(firstState_name, L"mm_results/mm_firstangle_vector_class_");
	_itow(classid, str, 10);
	wcscat(filename, str);
	wcscat(firstState_name, str);

	FILE *ff = _wfopen(firstState_name, L"rt");
	wprintf(L"Reading /%s...\n", firstState_name);
	//int vec_size2 = readMat(ff);
	float *firstAngleMat= readMat(ff);
	//float *firstAngleMat = new float[vec_size2];
	//firstAngleMat = fdata_carrier;

	FILE *fp = _wfopen(filename, L"rt");
	wprintf(L"Reading /%s...\n", filename);
	//int vec_size = readMat(fp);
	float *MMangleMat = readMat(fp);
	
	int dim = dirsize;
	// START PROCESSING
	// Starting with first State
	float dist;
	float max_value = 0, new_value;
	int index = 1;
	int requested_node = angleArray[0];
	// Get distances for each som node from requested node 
	// and multiply to their firstState probability
	for (int n = 0; n < dim; n++){
		
		dist = get_angle_distance(requested_node, n);
		new_value = (1-dist) * firstAngleMat[n];
		if (new_value > max_value) {
            max_value = new_value;
			index = n+1; //cause it starts from zero			
        }
	}
	int win_node = index;
	//if (max_value>0.01)
		//wprintf(L"First Max index: %d  with value: %f\n", win_node, max_value);
	cProb = max_value;	
	
	//Next nodes
	index = 1;
	for (int pos = 1; pos < transitions; pos++){
		requested_node = angleArray[pos];
		for (int n = 0; n < dim; n++){
			
			dist = get_angle_distance(requested_node, n);
			new_value = (1-dist) * MMangleMat[(win_node-1)*dim + n];
			if (new_value > max_value) {
				max_value = new_value;
				index = n+1; //cause it starts from zero			
			}
		}
		win_node = index;
		//if (max_value>0.01)
			//wprintf(L"Max index: %d  with value: %f\n", win_node, max_value);

		cProb += max_value;	
	}
	// END OF PROCESSING
	
	
	/*
	
	
	int counter = 0;
	int mat_pos =1;
	for (int pos = 0; pos < transitions-1; pos++){
		// Get Class transition vector
		mat_pos = (angleArray[pos]-1)*dim + (angleArray[pos+1] -1);
		//if (fdata_carrier[mat_pos] > 0)
		//cProb = cProb + fdata_carrier[mat_pos];
		cProb = cProb + MMangleMat[mat_pos];
		counter++;

	}// END OF PROCESSING
	if (counter>0)
		cProb = cProb/counter;
	*/
	return cProb;
}
////////////////////
float Markov::getLevenshteinProb(int *bmuArray, int transitions, int classid)
{
	float cProb = 0;
	wchar_t filename[_MAX_PATH] = L"";
	wchar_t meandist[_MAX_PATH] = L"";
	wchar_t str[_MAX_PATH] = L"";

	wcscat(filename, L"mm_results/lev_gmedian_class_");
	wcscat(meandist, L"mm_results/lev_meandist_class_");
	_itow(classid, str, 10);
	wcscat(filename, str);
	wcscat(meandist, str);

	FILE *fp = _wfopen(filename, L"rt");
	
	//int vec_size = readMat(fp);
	//float *LEVtransMat= readMat(fp);
	int gmedvec_size = readGMedian(fp);
	wprintf(L"Reading /%s ...\n", filename);
	int * LEVtransMat = int_data_carrier;
	
	// START PROCESSING
	int test_size = vec2str(bmuArray, transitions);
	wchar_t * test_str = new wchar_t[test_size];
	wcscpy(test_str, lev_string);

	int gm_size = vec2str(LEVtransMat, gmedvec_size);
	wchar_t * gmedian_str = new wchar_t[gm_size];
	wcscpy(gmedian_str, lev_string);
		

	FILE *md = _wfopen(meandist, L"rt");
	
	int distL;
	
	// Find Levenshtein distance L of test NodeVector towards GM
	distL = getLevenshteinDistance(test_str, test_size, gmedian_str, gm_size);
	// Read mean Levenshtein distance ML of all training NodeVectors towards their GM
	float distML = readMeanDist(fp);

	//wprintf(L"ML = %1.4f\n", distML);
	if (distL > 0)
		cProb =  distML/(float)distL;
	else
		cProb = 1;
	
	return cProb;
}

////////////////////////////////////////////////////////////////////////////////
int Markov::getTransitionVector(SOM *som, PREC classrec){
	
	int bmuIndex=0;
	int indexnode=0;
	int previousNode;

	bmuClassIndexArray = new int[(int)classrec->entries.size()];//(int *)malloc((int)classrec.entries.size() * sizeof(int));

	//for number of patterns/////////////////////////////////////////////////////////////////////////////
	for (int i = 0; i < (int)classrec->entries.size()-1; i++) {
			//find bmu for vector
			bmuIndex = som->getBMUindex(classrec->entries[i]->vec);

			//out result - bmu winning node
			if (i==0){
				bmuClassIndexArray[indexnode] = bmuIndex;
				indexnode++;
				// check if patterns get same index in a row
			}else {
				previousNode = bmuClassIndexArray[indexnode-1];
				if ( previousNode != bmuIndex ){
					bmuClassIndexArray[indexnode] = bmuIndex;
					indexnode++;
				}
			}
	//////////////////////////////////////////////////////////////
	}//end for (patterns)
	int vectorlength = indexnode;
	for (int i = 0; i < vectorlength; i++) {
		//wprintf(L"Pattern %d   ->   BMU: %d \n", i+1, bmuClassIndexArray[i]);
	}
	return vectorlength;
}
int Markov::getTransitionVector(SOM *som, PREC classrec, int irep){
	
	int bmuIndex=0;
	int indexnode=0;
	int previousNode;
	//wprintf(L"Rep = %d   \n", irep);


	bmuClassIndexArray = new int[(int)classrec->entries.size()];//(int *)malloc((int)classrec.entries.size() * sizeof(int));

	//for number of patterns/////////////////////////////////////////////////////////////////////////////
	for (int i = 0; i < (int)classrec->entries.size()-1; i++) {

		if ( classrec->entries[i]->rep == irep){
			//find bmu for vector
			bmuIndex = som->getBMUindex(classrec->entries[i]->vec);

			//out result - bmu winning node
			if (i==0){
				bmuClassIndexArray[indexnode] = bmuIndex;
				indexnode++;
				// check if patterns get same index in a row
			}else {
				previousNode = bmuClassIndexArray[indexnode-1];
				if ( previousNode != bmuIndex ){
					bmuClassIndexArray[indexnode] = bmuIndex;
					indexnode++;
				}
			}
		}
	//////////////////////////////////////////////////////////////
	}//end for (patterns)
	int vectorlength = indexnode;
	for (int i = 0; i < vectorlength; i++) {
		//wprintf(L"Pattern %d   ->   BMU: %d \n", i+1, bmuClassIndexArray[i]);
	}
	return vectorlength;
}

void Markov::getVectorTransitionMatrixPerClass(SOM *som, PREC rec, int reps)
{
    wchar_t str[_MAX_PATH] = L"";	
	wchar_t transitions_file[_MAX_PATH] = L"mm_results/";
	wchar_t gmedian_file[_MAX_PATH] = L"mm_results/";
	wchar_t meandist_file[_MAX_PATH] = L"mm_results/";
	int somsize = som->get_nodes_number();
	int classid = (int)rec->entries[0]->cls;
	_itow(classid, str, 10);
	wcscat(transitions_file, L"lev_trans_vector_class_");
	wcscat(gmedian_file, L"lev_gmedian_class_");
	wcscat(meandist_file, L"lev_meandist_class_");
	wcscat(transitions_file, str);
	wcscat(gmedian_file, str);
	wcscat(meandist_file, str);
	

	firstStateMat = new float[somsize]; //Pinakas Prwtis Katastasis
    	
	
	//Create and initialize Transition Matrix
	transitionMat = new float*[somsize];// Pinakas Metavasewn
	for (int s = 0; s < somsize; ++s)
		transitionMat[s] = new float[somsize];
	for (int i = 0; i < somsize; ++i){
		firstStateMat[i] = 0;
		for (int j = 0; j < somsize; ++j){
			transitionMat[i][j] = 0;
		}
	}
	int maxTransitions =1;
	int transitions;
	int startPos = 1;
	int ** temp_mat= new int*[reps]; 
	for (int s = 0; s < reps; ++s)
			temp_mat[s] = new int[100];
	// START REPETITIONS IF ANY
	for (int r = 1; r <= reps; r++){
		// Get Class transition vector	
		transitions = getTransitionVector(som, rec, r);
		if (transitions > maxTransitions)
			maxTransitions = transitions;
		int transitionstotal = transitions -1;	
		               
		temp_mat[r-1] = bmuClassIndexArray;

		int fromNode;
		int toNode;
		wprintf(L"Rep %d :: Transitions Total %d \n", r, transitionstotal);
		startPos = bmuClassIndexArray[0];
		//wprintf(L"FirstState Node %d \n", startPos);
		firstStateMat[startPos-1] = firstStateMat[startPos-1] + 1;
		for (int i = 1; i < transitions; i++) {
			fromNode = bmuClassIndexArray[i-1];
			toNode = bmuClassIndexArray[i];
			transitionMat[fromNode-1][toNode-1] = transitionMat[fromNode-1][toNode-1] + 1;
		}	
	}// END OF REPETITIONS

	// Run repetitions again to Get general median tables and save  
	FILE *tv = _wfopen(transitions_file, L"wt");
	int ** test_mat= new int*[reps]; 
	for (int s = 0; s < reps; ++s)
			test_mat[s] = new int[maxTransitions];
	//int * gmedian_mat = new int[maxTransitions];
	int minsod = 100, sod = 0, rep_index = 1;
	float meansod=0;
	for (int r = 1; r <= reps; r++){
		transitions = getTransitionVector(som, rec, r);

		//gmedian_mat 
		for (int i = 0; i < maxTransitions; i++) {
			if (i<transitions)
				test_mat[r-1][i] = temp_mat[r-1][i];
			else
				test_mat[r-1][i] = 0;
		}
		if (r == 1){	
			minsod = getSumOfDistance(som, rec, test_mat[r-1], maxTransitions, temp_mat, reps);	
			meansod = (float)minsod/reps;
		}

		sod = getSumOfDistance(som, rec, test_mat[r-1], maxTransitions, temp_mat, reps);
		
		if (sod < minsod){
			minsod = sod;
			rep_index = r;
			meansod = (float)minsod/reps;
		}
		//save all vectors
		if (tv != 0) { 
			// read table
			for (int i = 0; i < maxTransitions; i++) {
				fwprintf(tv, L"%d ", test_mat[r-1][i]);
			}
			fwprintf(tv, L"\n"); 
		}
	}
	
	fclose(tv);
	//save Winning (set median) bmuArray of all reps
	FILE *gm = _wfopen(gmedian_file, L"wt");
	if (gm != 0) { 
			// read table
			for (int i = 0; i < maxTransitions; i++) {
				fwprintf(gm, L"%d ", test_mat[rep_index-1][i]);
			}
			fwprintf(gm, L"\n"); 
		}
	fclose(gm);
	//save Mean Distance of gmedian to all reps vectors
	FILE *md = _wfopen(meandist_file, L"wt");
	if (md != 0) { 
		// read table
		fwprintf(md, L"%f\n", meansod);
	}
	fclose(md);
	
	// Divide to number of transitions
	for (int i = 0; i < somsize; ++i){
		float rowsum = 0.0;
		for (int j = 0; j < somsize; ++j){
				rowsum = rowsum + transitionMat[i][j];
				
		}
		if (rowsum>0){
			for (int j = 0; j < somsize; ++j){
				transitionMat[i][j] = transitionMat[i][j]/rowsum;					
			}
		}
	}

	for (int i = 0; i < somsize; ++i){
			firstStateMat[i] = firstStateMat[i]/reps;				
	}
	
	
	
	// Save
	wchar_t result_file[_MAX_PATH] = L"mm_results/";
	wchar_t firstState_file[_MAX_PATH] = L"mm_results/";
	
	
	wcscat(result_file, L"mm_trans_matrix_class_");
	wcscat(result_file, str);
	wcscat(firstState_file, L"mm_firststate_vector_class_");
	wcscat(firstState_file, str);
	

    saveMat(result_file, transitionMat, somsize);
	saveVec(firstState_file, firstStateMat, somsize);

	
}
/////////////////////////////////////////////////////////////////////////////
//Angle Direction
void Markov::getAngleTransitionMatrixPerClass(PREC rec, int reps)
{
    	
	firstAngleMat = new float[dirsize]; //Pinakas Prwtis Katastasis

	//Create and initialize Transition Matrix
	directionMat = new float*[dirsize];
	for (int s = 0; s < dirsize; ++s)
		directionMat[s] = new float[dirsize];
	for (int i = 0; i < dirsize; ++i){
		firstAngleMat[i] = 0;
		for (int j = 0; j < dirsize; ++j){
			directionMat[i][j] = 0;
		}
	}
	
	int classid = (int)rec->entries[0]->cls;
	int startPos;
	// START REPETITIONS IF ANY
	for (int r = 1; r <= reps; r++){
		int transitions = getTransitionAngle(rec, r);
		int transitionstotal = transitions -1;
		
		int fromNode;
		int toNode;
		//wprintf(L"FirstState Node %d \n", startPos);
		startPos = angleClassIndexArray[0];
		firstAngleMat[startPos-1] = firstAngleMat[startPos-1] + 1;
		//wprintf(L"Transitions Total %d \n", transitionstotal);
		for (int i = 1; i < transitions; i++) {
			fromNode = angleClassIndexArray[i-1];
			toNode = angleClassIndexArray[i];
			directionMat[fromNode-1][toNode-1] = directionMat[fromNode-1][toNode-1] + 1;
		}	
		
	}// END OF REPETITIONS
	// Display Matrix
	//print2DMarkovMatrix(directionMat, dirsize);
	
	// Divide to number of transitions
	for (int i = 0; i < dirsize; ++i){
		float rowsum = 0.0;
		for (int j = 0; j < dirsize; ++j){
				rowsum = rowsum + directionMat[i][j];
				
		}
		if (rowsum>0){
			for (int j = 0; j < dirsize; ++j){
				directionMat[i][j] = directionMat[i][j]/rowsum;					
			}
		}
	}
	// For for state angles
	for (int i = 0; i < dirsize; ++i){
			firstAngleMat[i] = firstAngleMat[i]/reps;				
	}
	
	// Save
	wchar_t result_file[_MAX_PATH] = L"mm_results/";
	wchar_t firstAngle_file[_MAX_PATH] = L"mm_results/";
	wchar_t str[_MAX_PATH] = L"";
	wchar_t str2[_MAX_PATH] = L"";
	
	_itow(classid, str, 10);
	_itow(classid, str2, 10);
	wcscat(result_file, L"mm_dir_matrix_class_");
	wcscat(result_file, str);
	wcscat(firstAngle_file, L"mm_firstangle_vector_class_");
	wcscat(firstAngle_file, str2);

    saveMat(result_file, directionMat, dirsize);
	saveVec(firstAngle_file, firstAngleMat, dirsize);
	

}

int Markov::getTransitionAngle(PREC rec){
	
	int angleIndex=0;
	int indexnode=0;
	int vectorlength=0;
	int previousNode=0;
	//wprintf(L"Items = %d   \n", (int)rec->entries.size());

	//float angleArray[] = {0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0};
	angleClassIndexArray = new int[(int)rec->entries.size()];//(int *)malloc((int)rec->entries.size() * sizeof(int));
	//int angleClassIndexArray2[30];
	//REC rec = *classrec;
	//for number of patterns
	for (int i = 1; i < (int)rec->entries.size(); i++) {
		//find bmu for vector
		float a = get_angle(rec->entries[i]->vec, rec->entries[i-1]->vec);
		int angleIndex = get_angle_index(a);
		//wprintf(L"angle = %.2f   index = %d\n", a, angleIndex);
		//out result - angle winning class
		if (indexnode==0){
			angleClassIndexArray[indexnode] = angleIndex;
			indexnode++;
		}
		else
		{
			previousNode = angleClassIndexArray[indexnode-1];
			if ( previousNode != angleIndex ){
				angleClassIndexArray[indexnode] = angleIndex;
				indexnode++;
			}
		}
		
		
	}//end for (patterns)
	
	vectorlength = indexnode;
	//for (int i = 0; i < vectorlength; i++) {
	//	wprintf(L"Pattern %d   ->   Angle Class: %d \n", i+1, angleClassIndexArray[i]);
	//}
	return vectorlength;
}
int Markov::getTransitionAngle(PREC rec, int irep){
	
	int angleIndex=0;
	int indexnode=0;
	int vectorlength=0;
	int previousNode=0;
	wprintf(L"Rep = %d   \n", irep);

	//float angleArray[] = {0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0};
	angleClassIndexArray = new int[(int)rec->entries.size()];//(int *)malloc((int)rec->entries.size() * sizeof(int));
	//int angleClassIndexArray2[30];
	//REC rec = *classrec;
	//for number of patterns
	for (int i = 1; i < (int)rec->entries.size(); i++) {
		
		if ( rec->entries[i]->rep == irep){
			//find angle between vectors
			float a = get_angle(rec->entries[i]->vec, rec->entries[i-1]->vec);
			int angleIndex = get_angle_index(a);
			//wprintf(L"angle = %.2f   index = %d\n", a, angleIndex);
			//out result - angle winning class
			if (indexnode==0){
				angleClassIndexArray[indexnode] = angleIndex;
				indexnode++;
			}
			else
			{
				previousNode = angleClassIndexArray[indexnode-1];
				if ( previousNode != angleIndex ){
					angleClassIndexArray[indexnode] = angleIndex;
					indexnode++;
				}
			}
		}
		
	}//end for (patterns)
	
	vectorlength = indexnode;
	for (int i = 0; i < vectorlength; i++) {
		//wprintf(L"Pattern %d   ->   Angle Class: %d \n", i+1, angleClassIndexArray[i]);
	}
	return vectorlength;
}
///////////////////
float Markov::get_angle(const float *vec1, const float *vec2){
	int dimension = 2;
	float distance = 0.0f;
	float x = 0.0f, y = 0.0f;
	// Method1: cos (theta) = dotprod(vec1, vec2)/(abs(vec1)*abs(vec2));
	//for (int w = 0; w < dimension; w++) {
           // distance += (vec1[w] * vec2[w]);
            //n1 += vec1[w] * vec1[w];
            //n2 += vec2[w] * vec2[w];
    //}	
	//float angle = acos(distance / (sqrt(n1)*sqrt(n2))); // in degrees
	
	// Method2:  angle = atan2(vec2.y, vec2.x) - atan2(vec1.y, v1.x)	
	//d1=atan2(vec1[1], vec1[0]);
	//d2=atan2(vec2[1], vec2[0]);
	//float angle = d2-d1;

	// Method3 (Gogo):  angle = atan2(vec2.y -vec1.y, vec2.x - vec1.x)	
	x = vec2[0] - vec1[0];
	y = vec2[1] - vec1[1];
	float angle = atan2(y,x)* 180 / PI;

	if (angle < 0)
		angle = 360 + angle;
    return angle;
}
float *Markov::get_angle3d(const float *vec1, const float *vec2){
	int dimension = 3;
	float x = 0.0f, y = 0.0f, z = 0.0f;
	float *angle = new float[2];

	// Method
	x = vec2[0] - vec1[0];
	y = vec2[1] - vec1[1];
	z = vec2[2] - vec1[2];
	angle[0] = atan2(y,x)* 180 / PI;
	angle[1] = atan2(z,y)* 180 / PI;
	angle[2] = atan2(x,z)* 180 / PI;
	for (int i = 0; i < dimension; i++){
		if (angle[i] < 0)
			angle[i] = 360 + angle[i];
	}
    return angle;
}
//////////////////
int Markov::get_angle_index(float angle){
	int angclass = 8;
	if (angle>=0 && angle<45) angclass = 1;
	else if (angle>=45 && angle<90) angclass = 2;
	else if (angle>=90 && angle<135) angclass = 3;
	else if (angle>=135 && angle<180) angclass = 4;
	else if (angle>=180 && angle<225) angclass = 5;
	else if (angle>=225 && angle<270) angclass = 6;
	else if (angle>=270 && angle<315) angclass = 7;
	else if (angle>=315 && angle<360) angclass = 8;

    return angclass;
}
int Markov::get_angle_index(float *angle){
	int angclass = 8;
/*	if (angle>=0 && angle<45) angclass = 1;
	else if (angle>=45 && angle<90) angclass = 2;
	else if (angle>=90 && angle<135) angclass = 3;
	else if (angle>=135 && angle<180) angclass = 4;
	else if (angle>=180 && angle<225) angclass = 5;
	else if (angle>=225 && angle<270) angclass = 6;
	else if (angle>=270 && angle<315) angclass = 7;
	else if (angle>=315 && angle<360) angclass = 8;
	*/

    return angclass;
}

float Markov::get_angle_distance(int node, int x_node){
	
	float dist = 0;
	int diff = abs(node - x_node);
	if (diff==0) dist = 1;
	else if (diff=1) dist = 1/2;
	else if (diff=2) dist = 1/4;
	else if (diff=3) dist = 1/8;
	else if (diff=4) dist = 1/16;
	else if (diff=5) dist = 1/32;
	else if (diff=6) dist = 1/64;
	else if (diff=7) dist = 1/128;

    return dist;
}
////// Utility functions /////////////////////////////////
int Markov::getSumOfDistance(SOM *som, PREC rec,int *test_mat, int size1, int ** temp_mat, int rows){
	int sum = 0;
	//wchar_t temp_str;
	
	//int size1 = (sizeof test_mat)/(sizeof test_mat[0]), 
	int size2, transitions;
	//wprintf(L"size1 = %d\n",size1);
	int size_s1  = vec2str(test_mat, size1 ), size_s2;
	
	wchar_t * str1 = new wchar_t[size_s1];
	wchar_t * str2;
	wcscpy(str1, lev_string);


	for (int i = 0; i < rows; i++){		
		transitions = getTransitionVector(som, rec, i+1);
		size2 = transitions;
		//wprintf(L"size2%d = %d\n",i,size2);
		size_s2 = vec2str(temp_mat[i], size2);
		str2 = new wchar_t[size_s2]; 
		wcscpy(str2, lev_string);

		sum += getLevenshteinDistance(str1, size_s1, str2, size_s2);
	}
	
	//wprintf(L"size s1 = %d, str1 = %s, size s2 = %d, str2 = %s   sum = %d\n", size_s1,str1,size_s2,str2, sum);
	
	return sum;
}
// Convert a vector to a string
int Markov::vec2str(int * vec, int matsize){
	wchar_t mystring[_MAX_PATH];
	wchar_t node[_MAX_PATH];
	
	for (int i = 0; i < matsize; i++){
		_itow(vec[i], node, 10);
		if (i==0) wcscpy(mystring, node);
		else
			wcscat(mystring, node);
		wcscat(mystring, L"_");
		//wprintf(L"String = %s , sz:%d \n", mystring, wcslen(mystring));
	}
	int index = 0;
	int strlength = wcslen(mystring);

	wcscpy(lev_string, mystring);
	return strlength;
}


/////////////
int Markov::getLevenshteinDistance(wchar_t *s1, int size_s1, wchar_t *s2, int size_s2)
 {
    int m = size_s1;
	int n = size_s2;
	 // for all i and j, d[i,j] will hold the Levenshtein distance between
    // the first i characters of s1 and the first j characters of s2,(t);
    // note that d has (m+1)x(n+1) values
    int ** d =new int*[m+1];
	for (int i = 0; i < m +1; ++i)
		d[i] = new int[n+1]; 

  
    for (int i = 0; i <= m; i++)
      d[i][0] = i; // the distance of any first string to an empty second string
    for (int j = 0; j <= n; j++)
	  d[0][j] = j; // the distance of any second string to an empty first string
  
    for (int j = 1; j <= n; j++)
    {
		for (int i = 1; i <= m; i++)
		{
			if (s1[i] == s2[j]){  
				d[i][j] = d[i-1][j-1] ;  // no operation required
			} 
			else
			{
			   //d[i-1, j] + 1,  // a deletion
			   //d[i, j-1] + 1,  // an insertion
			   //d[i-1, j-1] + 1 // a substitution
				int del = d[i-1][j] + 1;
				int ins = d[i][j-1] + 1; 
				int sub = d[i-1][j-1] + 1;
	
				d[i][j] = min (del, min(ins, sub));
			}
		}
	}
  
   return d[m][n];
 }
      
int Markov::readGMedian(FILE *fp)
{
	int tmp;
	vector<int> fvec;

	while (fwscanf(fp, L"%d", &tmp) == 1){
		fvec.push_back(tmp);
	}
	
	int vector_length = (int)fvec.size();
	// pass data to global carrier
	int_data_carrier = new int[vector_length];
	for (int i = 0; i < vector_length; i++){
		 int_data_carrier[i] = fvec[i];

	}
	fclose(fp);

	return vector_length;
}

float* Markov::readMat(FILE *fp)
{
			
	float tmp;
	vector<float> fvec;

	while (fwscanf(fp, L"%f", &tmp) == 1){
		fvec.push_back(tmp);
	}

	int vector_length = (int)fvec.size();
	float * fdata_carrier = new float[vector_length];
	for (int i = 0; i < vector_length; i++){
		 fdata_carrier[i] = fvec[i];

	}
	fclose(fp);
	
	return fdata_carrier;
}
float Markov::readMeanDist(FILE *fp)
{
			
	float tmp;
	float distML;

	while (fwscanf(fp, L"%f", &tmp) == 1){
		distML = tmp;
	}

	fclose(fp);
	
	return distML;
}

int Markov::saveMat(const wchar_t *fname, float **matrix, int m_size) const
{        
        FILE *fp = _wfopen(fname, L"wt");
		//int m_size = (sizeof matrix);
		wprintf(L"Saving matrix size %d x %d...\n", m_size, m_size);
        if (fp != 0) {
                
                //save transition Matrix
                for (int n = 0; n < m_size; n++) {
                      
                        for (int w = 0; w < m_size; w++)
                                fwprintf(fp, L"%.4f	", matrix[n][w]);
                        fwprintf(fp, L"\n");
                }

                fclose(fp);
                return 0;
        } else
                return -1;
}
int Markov::saveVec(const wchar_t *fname, float *fvector, int m_size) const
{        
        FILE *fp = _wfopen(fname, L"wt");
		//int m_size = (sizeof matrix);
		wprintf(L"Saving vector size %d ...\n", m_size);
        if (fp != 0) {
                
                //save transition Matrix
             
                        for (int w = 0; w < m_size; w++)
                                fwprintf(fp, L"%.4f	", fvector[w]);
                        fwprintf(fp, L"\n");


                fclose(fp);
                return 0;
        } else
                return -1;
}
int Markov::saveMarkov(const wchar_t *fname, int classnum, int somsize ) const
{        
        FILE *fp = _wfopen(fname, L"wt");
		//int m_size = (sizeof matrix);

        if (fp != 0) {
                
                //save class number
              fwprintf(fp, L"%d	", classnum);
              fwprintf(fp, L"\n");
			  fwprintf(fp, L"%d	", somsize);
              fwprintf(fp, L"\n");


                fclose(fp);
                return 0;
        } else
                return -1;
}

/////////////////////////////////////////////////////////////////////////////
void Markov::print2DMarkovMatrix(float **mat, int matsize){
	for (int i = 0; i < matsize; ++i){
		for (int j = 0; j < matsize; ++j){
			wprintf(L"%.2f ", mat[i][j]);
		}
		wprintf(L"\n");
	}
}
