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

#include "ann.h"
void print_arr(double* arr,int len){
	int i;
	printf("[");
	for(i=0;i<len;i++)
		printf(" %f",arr[i]);
	printf(" ]");
}

void print_arr_arr(double** arr, int rows, int cols){
	int i,j;
	for(i=0;i<rows;i++){
		printf("[");
		for(j=0;j<cols;j++)
			printf(" %f",arr[i][j]);
		printf(" ]\n");
	}
}

void error(double* output, double* correct, double* error, int len){
	int i;
	for(i=0;i<len;i++)
		error[i] = correct[i] - output[i];
}

double threshold(double x){
	if(x>0.0) return 1.0;
	return -1.0;
}

double sigmoid(double x){
	return 1.0 / (1.0 + exp(-x));
}


double sigderiv(double x){
	return sigmoid(x)*(1.0 - sigmoid(x));
}

/*
 * The Perceptron
 */

void init_perceptron(unsigned int inputs, unsigned int outputs, struct perceptron *p){
	int i,j;
	p->num_input = inputs;
	p->num_output = outputs;
	p->weights = (double**)(malloc(inputs*sizeof(double*)));
	for(i=0;i<inputs;i++)
		p->weights[i] = (double*)(malloc(outputs*sizeof(double)));
	for(i=0;i<inputs;i++)
		for(j=0;j<outputs;j++)
			p->weights[i][j] = (2*drand48())-1;
}

void free_perceptron(struct perceptron *p){
	int i;
	for(i=0;i<p->num_input;i++)
		free(p->weights[i]);
	free(p->weights);
}

void run_perceptron(struct perceptron *p, double* input, double* output){
	int i,j;
	for(j=0;j<p->num_output;j++)
		output[j] = 0.0;
	for(i=0;i<p->num_input;i++)
		for(j=0;j<p->num_output;j++)
			output[j] += p->weights[i][j] * input[i];
	for(j=0;j<p->num_output;j++)
		output[j] = threshold(output[j]);
}

void update_perceptron(struct perceptron *p, double *input, double *output,
						double *correct_output, double nu){
	int i,j;
	run_perceptron(p,input,output);
	error(output,correct_output,output,p->num_output);
	for(i=0;i<p->num_input;i++)
		for(j=0;j<p->num_output;j++)
			p->weights[i][j] = p->weights[i][j] + nu * output[j] * input[i];
}

void train_perceptron(struct perceptron *p, double **inputs, double **correct_outputs,
						double nu, int num_examples){
	int i;
	double* output = (double *)(malloc(p->num_output*sizeof(double)));
	for(i=0;i<num_examples;i++){
		update_perceptron(p,inputs[i],output,correct_outputs[i],nu);
	}
	free(output);
}

/*
 * Two layer feed-forward ANN
 */

void init_two_layer_NN(unsigned int inputs, unsigned int outputs, unsigned int hidden,
						struct two_layer_NN *p){
	int i,j,k;
	p->num_input = inputs;
	p->num_output = outputs;
	p->num_hidden = hidden;
	
	p->tmp = (double*)(malloc(hidden*sizeof(double)));
	
	p->input_hidden_weights = (double**)(malloc(inputs*sizeof(double*)));
	for(i=0;i<inputs;i++)
		p->input_hidden_weights[i] = (double*)(malloc(hidden*sizeof(double)));
	p->hidden_output_weights = (double**)(malloc(hidden*sizeof(double*)));
	for(i=0;i<hidden;i++)
		p->hidden_output_weights[i] = (double*)(malloc(outputs*sizeof(double)));
	for(i=0;i<inputs;i++)
		for(j=0;j<hidden;j++)
			p->input_hidden_weights[i][j] = (2*drand48())-1;
	for(j=0;j<hidden;j++)
		for(k=0;k<outputs;k++)
			p->hidden_output_weights[j][k] = (2*drand48())-1;
}

void free_two_layer_NN(struct two_layer_NN *p){
	int i;
	for(i=0;i<p->num_input;i++)
		free(p->input_hidden_weights[i]);
	free(p->input_hidden_weights);
	for(i=0;i<p->num_hidden;i++)
		free(p->hidden_output_weights[i]);
	free(p->hidden_output_weights);
	free(p->tmp);
}

void run_two_layer_NN(struct two_layer_NN *p, double *input, double *output){
	int i,j,k;
	for(k=0;k<p->num_output;k++)
		output[k] = 0.0;
	for(j=0;j<p->num_hidden;j++)
		p->tmp[j] = 0.0;
	//input to hidden layer
	for(i=0;i<p->num_input;i++)
		for(j=0;j<p->num_hidden;j++)
			p->tmp[j] += p->input_hidden_weights[i][j] * input[i];
	//for(j=0;j<p->num_hidden;j++)
	//	p->tmp[j] = sigmoid(p->tmp[j]);
	for(j=0;j<p->num_hidden;j++)
		for(k=0;k<p->num_output;k++)
			output[k] += p->hidden_output_weights[j][k] * sigmoid(p->tmp[j]);
	for(k=0;k<p->num_output;k++)
		output[k] = sigmoid(output[k]);
}

void update_two_layer_NN(struct two_layer_NN* p, double* input, double* output,
							double* correct_output, double nu){
	int i,j,k;
	double* err = (double*)(malloc(p->num_output*sizeof(double)));
	run_two_layer_NN(p,input,output);
	error(output,correct_output,err,p->num_output);
	for(j=0;j<p->num_hidden;j++){
		for(k=0;k<p->num_output;k++)
			p->hidden_output_weights[j][k]
				= p->hidden_output_weights[j][k] 
				+ nu * sigmoid(p->tmp[j]) * err[k] * output[k]*(1.0 - output[k]);
	}
	for(i=0;i<p->num_input;i++){
		for(j=0;j<p->num_hidden;j++){
			double sum=0.0;
			for(k=0;k<p->num_output;k++){
				sum += p->hidden_output_weights[j][k] * err[k] * output[k]*(1.0 - output[k]);
			}
			//printf("*** %f ***\n",sigmoid(p->tmp[j]));
			sum *= sigderiv(p->tmp[j]);
			p->input_hidden_weights[i][j]
				= p->input_hidden_weights[i][j]
				+ nu * input[i] * sigmoid(p->tmp[j]) * sum;
		}
	}
	free(err);
}

void train_two_layer_NN(struct two_layer_NN* p, double **inputs, double** correct_outputs,
						double nu, int num_examples){
	int i;
	double* output = (double *)(malloc(p->num_output*sizeof(double)));
	for(i=0;i<num_examples;i++){
		update_two_layer_NN(p,inputs[i],output,correct_outputs[i],nu);
	}
	free(output);
}