#include "mlp.h"


mlp::mlp(gsl_vector *neurons, activ_function *act_f, bool softmax): f(act_f), neurons(neurons), softmax(softmax)
{
	nInputs = (int)gsl_vector_get(neurons, 0);
	nOutputs = (int)gsl_vector_get(neurons, neurons->size - 1);
		
	W = (gsl_matrix**)malloc(sizeof(gsl_matrix*) * (neurons->size - 1));

	activationHidden = NULL;

	for(unsigned int n = 0; n < neurons->size - 1; n++)
	{
		W[n] = gsl_matrix_alloc((int)gsl_vector_get(neurons, n + 1), (int)gsl_vector_get(neurons, n));

		for(unsigned int i = 0; i < W[n]->size1; i++)
		{
			for(unsigned int j = 0; j < W[n]->size2; j++)
			{
				gsl_matrix_set(W[n], i, j, (rand() / (double)RAND_MAX) - 0.5);
			}
		}
	}
}

mlp::~mlp(void)
{
	for(unsigned int n = 0; n < neurons->size - 1; n++)
	{
		gsl_matrix_free(W[n]);
	}
	gsl_vector_free(neurons);
	if(activationHidden != NULL){
		gsl_matrix_free(activationHidden);
		activationHidden = NULL;
	}

	delete &f;
}

gsl_matrix *mlp::runNetwork(gsl_matrix *inputs){
	gsl_matrix *h_new;
	gsl_matrix *h_prev = gsl_matrix_alloc(inputs->size1, inputs->size2);
	gsl_matrix_memcpy(h_prev, inputs);

	gsl_vector *v;
	for(unsigned int n = 0; n < neurons->size - 1; n++){
		h_new = gsl_matrix_alloc(gsl_vector_get(neurons, n + 1), inputs->size2);
		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, W[n], h_prev, 0.0, h_new);

		// Softmax na vystupe
		if(softmax && (n == (neurons->size - 2))){
			f->applySoftmax(h_new);
		}else{
			f->apply(h_new);
		}

		//bias
		if(n != neurons->size - 2){ // test na (ne)vystupnu vrstvu
			v = gsl_vector_alloc(inputs->size2);
			gsl_vector_set_all(v, -1);
			gsl_matrix_set_row(h_new, h_new->size1 - 1, v);
			gsl_vector_free(v);
		}

		gsl_matrix_free(h_prev);
		h_prev = h_new;
		
		// odkladanie skrytej aktivacie
		if(n == 0) //prva skryta vrstva
		{	
			if(activationHidden != NULL){
				gsl_matrix_free(activationHidden);
				activationHidden = NULL;
			}

			activationHidden = gsl_matrix_alloc(h_new->size1, h_new->size2);
			gsl_matrix_memcpy(activationHidden, h_new);
		}
		//koniec odkladania
	}
	return h_prev;
}

gsl_vector *mlp::runNetwork(gsl_vector *input){
	gsl_matrix *inputs = gsl_matrix_alloc(input->size, 1);
	gsl_matrix_set_col(inputs, 0, input);
	
	gsl_matrix *outputs = runNetwork(inputs);
	gsl_vector *result = gsl_vector_alloc(outputs->size1);
	gsl_matrix_get_col(result, outputs,0);
	gsl_matrix_free(outputs);
	gsl_matrix_free(inputs);

	return result;
}
void mlp::saveNetwork(char *path){
	for(unsigned int i = 0; i < neurons->size - 1; i++){

		mkdir(path);

		stringstream filename;
		filename << path; 
		filename <<"W_";
		filename << i;
		filename << ".txt";

		char *str = (char*)malloc(sizeof(char) * filename.str().length() + 1);
		strcpy(str, filename.str().c_str());
	
		FILE *f = fopen(str, "w");
		gsl_matrix_fprintf (f, W[i], "%e");
		fclose(f);
	
		free(str);
	}
}

void mlp::loadNetwork(char *path){
	for(unsigned int i = 0; i < neurons->size - 1; i++){

		stringstream filename;
		filename << path; 
		filename <<"W_";
		filename << i;
		filename << ".txt";

		char *str = (char*)malloc(sizeof(char) * filename.str().length() + 1);
		strcpy(str, filename.str().c_str());

		FILE *f = fopen(str, "r");
		gsl_matrix_fscanf(f, W[i]);
		fclose(f);

		free(str);
	}
}
void mlp::trainNetwork(gsl_matrix *inputs, gsl_matrix *desired, double alpha){

	gsl_vector **h = (gsl_vector**)malloc(sizeof(gsl_vector*) * neurons->size);
	gsl_vector **delta = (gsl_vector**)malloc(sizeof(gsl_vector*) * neurons->size);
//print_matrix(inputs);	
	//alokovanie aktivacii a delt siete
	for(unsigned int n = 0; n < neurons->size; n++){
		h[n] = gsl_vector_alloc(gsl_vector_get(neurons, n));
		delta[n] = gsl_vector_alloc(gsl_vector_get(neurons, n));
	}
	//gsl_vector *d = gsl_vector_alloc(desired->size1);
	
	for(unsigned int step = 0; step < inputs->size2; step++){
		//gsl_matrix_get_col(d, desired, step);
		gsl_matrix_get_col(h[0], inputs, step);
		//presirenie vstupu
		for(unsigned int n = 0; n < neurons->size - 1; n++){
			//h[n + 1] = gsl_vector_alloc(gsl_vector_get(neurons, n + 1));
			//print_matrix(W[n]);

			gsl_blas_dgemv(CblasNoTrans, 1.0, W[n], h[n], 0.0, h[n + 1]);
			
			// Softmax na vystupe
			if(softmax && (n == (neurons->size - 2))){
				f->applySoftmax(h[n + 1]);
			}else{
				f->apply(h[n + 1]);
			}


			if(n != neurons->size - 2)
				gsl_vector_set(h[n + 1], h[n+1]->size - 1, -1);

			//print_vector(h[n+1]);
		}

		//delta[neurons->size - 1] = gsl_vector_alloc(gsl_vector_get(neurons, neurons->size - 1));


		/*		
		gsl_matrix_get_col(delta[neurons->size - 1], neurons, step);
		gsl_vector_sub(delta[neurons->size - 1], h[neurons->size - 1]);
		f->apply(h[neurons->size - 1]);
		gsl_vector_mul(delta[neurons->size - 1], h[neurons->size - 1]);
	*/
		// vystupna delta
		for(unsigned int i = 0; i < delta[neurons->size - 1]->size; i++){
			double y = gsl_vector_get(h[neurons->size - 1], i);
			if(!softmax){
				gsl_vector_set(delta[neurons->size - 1], i, (gsl_matrix_get(desired, i, step) - y) * f->actFunction_deriv(y)); 
			}else{
				gsl_vector_set(delta[neurons->size - 1], i, (gsl_matrix_get(desired, i, step) - y));
			}
		}

		for(unsigned int n = neurons->size - 2; n > 0 ; n--){
			//delta[n] = gsl_vector_alloc(gsl_vector_get(neurons, n));
			// ostatna delta
			for(unsigned int k = 0; k < delta[n]->size; k++){
				double delta_k = 0;
				double y = gsl_vector_get(h[n], k);
				for(unsigned int i = 0; i < delta[n + 1]->size; i++){
					delta_k += gsl_vector_get(delta[n + 1], i) * gsl_matrix_get(W[n], i, k);
					gsl_vector_set(delta[n], k, delta_k * f->actFunction_deriv(y)); 
				}
			}
		}

		//uprava vystupnych vah
		for(unsigned int i = 0; i < h[neurons->size - 1]->size; i++){
			for(unsigned int k = 0; k < h[neurons->size - 2]->size; k++){
				gsl_matrix_set(W[neurons->size - 2], i, k, gsl_matrix_get(W[neurons->size - 2], i, k) + alpha * gsl_vector_get(delta[neurons->size - 1], i) * gsl_vector_get(h[neurons->size - 2], k));
			}
		}

		for(unsigned int n = neurons->size - 2; n > 0 ; n--){
			//delta[n] = gsl_vector_alloc(gsl_vector_get(neurons, n));

			//uprava vah
			for(unsigned int k = 0; k < h[n]->size - 1; k++){
				for(unsigned int j = 0; j < h[n - 1]->size; j++){
					gsl_matrix_set(W[n - 1], k, j, gsl_matrix_get(W[n - 1], k, j) + alpha * gsl_vector_get(delta[n], k) * gsl_vector_get(h[n - 1], j));
				}
			}
		}
	}

	//dealokovanie aktivacii a delt siete
	for(unsigned int n = 0; n < neurons->size - 1; n++){
		gsl_vector_free(h[n]);
		gsl_vector_free(delta[n]);
	}

}
void mlp::trainNetwork(gsl_vector *input, gsl_vector *desired, double alpha){
	gsl_matrix *inputs = gsl_matrix_alloc(input->size, 1);
	gsl_matrix *des = gsl_matrix_alloc(desired->size, 1);
	gsl_matrix_set_col(inputs, 0, input);
	gsl_matrix_set_col(des, 0, desired);

	trainNetwork(inputs, des, alpha);
	
	gsl_matrix_free(inputs);
	gsl_matrix_free(des);

}
double mlp::compute_error(gsl_matrix *desired, gsl_matrix *output){
	// RMSE error sqrt(sum((d-x).^2) / n);
	if(desired->size1 != output->size1 || desired->size2 != desired->size2)
	{
		printf("WARNING: Nesuhlasia rozmery matic!");
		return 0;
	}
	double error = 0;

	for(unsigned int i = 0; i < desired->size1; i++){
		for(unsigned int j = 0; j < desired->size2; j++){
			double d = gsl_matrix_get(desired, i, j);
			double y = gsl_matrix_get(output, i, j);
			error += (d - y) * (d - y);
		}
	}
	return sqrt(error / ((double)desired->size1 * (double)desired->size2));
}