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

char *file_pesi = "rete_pesi.bin";
char *file_immagini = "mnist/train-images.idx3-ubyte";
char *file_label = "mnist/train-labels.idx1-ubyte";

//Test dataset
/* char *file_immagini = "mnist/t10k-images.idx3-ubyte";
char *file_label = "mnist/t10k-labels.idx1-ubyte"; */

#include "mnist/mnist_manager.h"
// #include "cifar_10/cifar10_manager.h"
//#include "xor_manager.h"

// Siccome il char è un byte che rappresenta il valore tra 0 e 255. Per evitare confusioni definisco il tipo "byte" come in Java
typedef unsigned char byte;

double LRE = 1.414;
double soglia_sigmoide = 0.5;

typedef struct
{
    double *pesi;
    double bias;
    int size;
} Percettrone;

typedef struct
{
    Percettrone *percettroni;
    int size;
} Layer;

typedef struct
{
    Layer *layers;
    int size;
} ReteNeurale;

double randomico();
double *get_double_from_bytes(Istanza);

Percettrone inzializza_percettrone(int);
Layer inizializza_layer(int, int);
ReteNeurale inizializza_rete_neurale_deep_feed_forward(int, int, int);
ReteNeurale inizializza_rete_neurale_feed_forward(int, int, int);


double funzione_attivazione_percettrone(Percettrone p, double*, int);
double **elabora_funzioni_attivazione(ReteNeurale, Istanza, int);
double *elabora_funzioni_attivazione_layer(Layer, double*, int);
double **discesa_gradiente(ReteNeurale, double **, double, int);
double *calcola_gradiente_output(Layer layer, double*, double, int);
double calcola_gradiente_disceso(ReteNeurale, int, int, double**, double);
void aggiorna_pesi(ReteNeurale, double**, double**, Istanza);
void correggi_pesi_percettrone(Percettrone*, double*, double);
//void correggi_pesi_percettrone_byte(Percettrone*, Istanza, double, int);

int previsione(double);
double sigmoide(double);
double derivata_sigmoide(double);
double relu(double);
double derivata_relu(double);

void salvaReteNeurale(const char *, ReteNeurale *);
ReteNeurale *caricaReteNeurale(const char *);

void debug(double);
void debug_matrice(ReteNeurale*, double**);
void debug_vettore(double*, int);

/*
    ################# INIZIALIZZAZIONI E METODI UTILI ################################
*/
// Questa funzione genera un valore reale random compreso nell'intervallo [-1, 1]
double randomico()
{
    // Genero numeri nell'intervallo [-1,1]
    return ((double)(rand() % 101 * 0.01 * 2.0) -1.0);
}

Percettrone inizializza_percettrone(int n_pesi)
{
    Percettrone p;
    p.pesi = (double *)malloc(sizeof(double) * n_pesi);
    for (int i = 0; i < n_pesi; i++) {
        p.pesi[i] = randomico();
    }

    p.bias = randomico();

    p.size = n_pesi;

    return p;
}

Layer inizializza_layer(int n_percettroni, int n_pesi)
{
    Layer layer;
    layer.percettroni = (Percettrone *)malloc(sizeof(Percettrone) * n_percettroni);

    for (int i = 0; i < n_percettroni; i++)
    {
        layer.percettroni[i] = inizializza_percettrone(n_pesi);
    }

    layer.size = n_percettroni;

    return layer;
}

ReteNeurale inizializza_rete_neurale_feed_forward(int numero_layers, int numero_percettroni_iniziali, int numero_input) {
    ReteNeurale r;
    r.layers = (Layer *)malloc(sizeof(Layer) * numero_layers);
    r.size = numero_layers;

    // Funzione esponenziale inversa layer 5
    for (int livello = 0; livello < numero_layers; livello++)
    {
        double esponente = (double)livello / (double)numero_layers;
        double frazione = (double)1 / (double)numero_percettroni_iniziali;

        int numero_percettroni_livello = (int)((double)numero_percettroni_iniziali * pow(frazione, esponente));
        if (livello == numero_layers - 1)
            numero_percettroni_livello = 1;

        // printf("esponente %f -> frazione: %f\n", esponente, frazione);
        printf("Layer %d -> percettroni: %d\n", livello, numero_percettroni_livello);

        if (livello == 0)
            r.layers[livello] = inizializza_layer(numero_percettroni_livello, numero_input);
        else
            r.layers[livello] = inizializza_layer(numero_percettroni_livello, r.layers[livello - 1].size);
    }

    return r;
}

ReteNeurale inizializza_rete_neurale_deep_feed_forward(int numero_layers, int numero_percettroni_iniziali, int numero_input)
{
    ReteNeurale r;
    r.layers = (Layer *)malloc(sizeof(Layer) * numero_layers);
    r.size = numero_layers;

    for (int livello = 0; livello < numero_layers; livello++)
    {
        if (livello == 0) {
            printf("Layer %d -> percettroni: %d\n", livello, numero_percettroni_iniziali);
            r.layers[livello] = inizializza_layer(numero_percettroni_iniziali, numero_input);
        }
        else if(livello < numero_layers -1) {
            printf("Layer %d -> percettroni: %d\n", livello, numero_percettroni_iniziali*2);
            r.layers[livello] = inizializza_layer(numero_percettroni_iniziali * 2, r.layers[livello - 1].size);
        }
        else {
            printf("Layer %d -> percettroni: %d\n", livello, 1);
            r.layers[livello] = inizializza_layer(1, r.layers[livello - 1].size);
        }
    }    

    return r;
}

double *get_double_from_bytes(Istanza istanza) {
    double *inputs = (double *)malloc(sizeof(double *) * istanza.size);
    for(int i = 0; i < istanza.size; i++) {
        inputs[i] = (double)istanza.dati[i];
    }

    return inputs;
}
/*
    ################# FUNZIONI ATTIVAZIONE ################################
*/

double sigmoide(double valore) {
    return 1.0 / (1.0 + exp(-valore));
}

double derivata_sigmoide(double valore) {
    return valore * (1.0 - valore);
}

double relu(double valore) {
    if(valore > 0)
        return valore;
    
    return 0;
}

double derivata_relu(double valore) {
    if(valore > 0)
        return 1;

    return 0;
}


//tipo_funzione 1: relu, 2:sigmoide

double **elabora_funzioni_attivazione(ReteNeurale rete, Istanza istanza, int tipo_funzione)
{
    //Inizializzo il vettore bidimensionale che dovrò ritornare
    double **funzioni = (double **)malloc(sizeof(double *) * rete.size);

    //Il primo layer devo farlo a parte perchè prende gli input dal dataset e non dal layer precedente
    funzioni[0] = elabora_funzioni_attivazione_layer(rete.layers[0], get_double_from_bytes(istanza), tipo_funzione);
    
    for(int indice_layer = 1; indice_layer < rete.size - 1; indice_layer ++) {
        funzioni[indice_layer] = elabora_funzioni_attivazione_layer(rete.layers[indice_layer], funzioni[indice_layer-1], tipo_funzione);   
    }

    funzioni[rete.size-1] = elabora_funzioni_attivazione_layer(rete.layers[rete.size-1], funzioni[rete.size-2], 2);

    return funzioni;
}



//DA NAN
double *elabora_funzioni_attivazione_layer(Layer layer, double *inputs, int tipo_funzione) {

    double *funzioni = (double *)malloc(sizeof(double) * layer.size);

    /* for(int indice_percettrone = 0; indice_percettrone < layer.size; indice_percettrone ++) {
        funzioni[indice_percettrone] = funzione_attivazione_percettrone(layer.percettroni[indice_percettrone], inputs, tipo_funzione);
    } */

    debug_vettore(funzioni, layer.percettroni->size);
    return funzioni;
}

//Tipo 1 = relu, tipo 2 = sigmoide
double funzione_attivazione_percettrone(Percettrone p, double *valori, int tipo)
{
    double sommatoria = 0.0;
    for (int i = 0; i < p.size; i++) {
        sommatoria += (valori[i] * p.pesi[i]);
        //printf("valore [%f] peso[%f] ", valori[i], p.pesi[i]);
    }

    sommatoria += p.bias;
    //printf(" sommatoria %f\n",sommatoria);

    if(tipo == 1)
        //funzione = relu(sommatoria);
        return relu(sommatoria);
    else
     //funzione = sigmoide(sommatoria);
     return sigmoide(sommatoria);
}

/*
    ################# RETROPROPAGAZIONE ################################
*/

//NON CALCOLA I GRADIENTI, PARTE DIRETTAMENTE NAN, VERIFICARE LA PRESENZA DI DATI IN FUNZIONI E LE FUNZIONI DI CALCOLO

double **discesa_gradiente(ReteNeurale rete, double **funzioni, double errore, int tipo_derivata) {

    double **gradienti = (double**)malloc(sizeof(double*) * rete.size);

    //Determino il gradiente di output a parte perchè non prende gradienti discesi dal livello superiore
    gradienti[rete.size-1] = calcola_gradiente_output(rete.layers[rete.size-1], funzioni[rete.size-1], errore, tipo_derivata);

    //Determino gli altri livelli
    for (int indice_layer = rete.size - 2; indice_layer >= 0; indice_layer--)
    {
        gradienti[indice_layer] = (double*)malloc(sizeof(double) * rete.layers[indice_layer].size);
        
        for (int indice_percettrone = 0; indice_percettrone < rete.layers[indice_layer].size; indice_percettrone++)
        {
            //In base al tipo di funzione scelto, avvio relu o sigmoidi
            double derivata_attivazione;

            if(tipo_derivata == 1)
                derivata_attivazione = derivata_relu(funzioni[indice_layer][indice_percettrone]);
            else
                derivata_attivazione = derivata_sigmoide(funzioni[indice_layer][indice_percettrone]);

            gradienti[indice_layer][indice_percettrone] = calcola_gradiente_disceso(rete, indice_layer + 1, indice_percettrone, gradienti, derivata_attivazione);
        }
    }

    //debug_matrice(&rete, gradienti);

    return gradienti;
}

 double *calcola_gradiente_output(Layer layer, double *funzioni, double errore, int tipo_derivata) {
    
    double *gradienti = (double*)malloc(sizeof(double) * layer.size);

    //Qui darò per scontato che ci sia un singolo percettrone in uscita, da aggiornare in caso di multipercettrone

    double derivata_funzione;

    if(tipo_derivata == 1)
        derivata_funzione = derivata_relu(funzioni[0]);
    else
        derivata_funzione = derivata_sigmoide(funzioni[0]);

    gradienti[0] = errore * derivata_funzione;

    //debug_vettore(gradienti, layer.size);

    return gradienti;
}

double calcola_gradiente_disceso(ReteNeurale rete, int livello, int indice_peso, double **gradienti, double derivata_attivazione)
{
    // printf("Qui ci arrivo\n");
    double sommatoria = 0.0;
    // printf("Layer %d: N_percettroni: %d\n", livello, rete.layers[livello].size);
    // Calcolo la sommatoria dei gradienti dei percettroni per i pesi
    for (int indice_percettrone = 0; indice_percettrone < rete.layers[livello].size; indice_percettrone++)
    {
        sommatoria += (gradienti[livello][indice_peso] * rete.layers[livello].percettroni[indice_percettrone].pesi[indice_peso]);
    }

    return sommatoria * derivata_attivazione;
}

/*
    ################# PREVISIONE E CORREZIONI ################################
*/

void aggiorna_pesi(ReteNeurale rete_neurale, double **gradienti, double **funzioni_attivazione, Istanza istanza) {

    //debug_matrice(&rete_neurale, gradienti);

    // Applico la correzione dal penultimo layer andando indietro fino al secondo (il primo si fa diverso)
    for (int indice_layer = rete_neurale.size - 1; indice_layer >= 0; indice_layer--)
    {
        // Applico la correzione a tutti i percettroni del layer dal primo a seguire
        for (int indice_percettrone = 0; indice_percettrone < rete_neurale.layers[indice_layer].size; indice_percettrone++)
        {
            // Devo prendere il gradiente del percettrone e moltiplicarlo con gli input associati ai pesi
            if (indice_layer > 0)
            {
                correggi_pesi_percettrone(&rete_neurale.layers[indice_layer].percettroni[indice_percettrone], funzioni_attivazione[indice_layer-1], gradienti[indice_layer][indice_percettrone]);
            }
            else
            {
                correggi_pesi_percettrone(&rete_neurale.layers[indice_layer].percettroni[indice_percettrone], get_double_from_bytes(istanza), gradienti[indice_layer][indice_percettrone]);
            }
        }
    }
}

int previsione(double valore)
{
    if (valore >= soglia_sigmoide)
        return 1;
    else
        return 0;
}

/*
    Al secondo giro diventa NAN
*/

void correggi_pesi_percettrone(Percettrone *p, double *input, double gradiente_percettrone)
{
    for (int indice_peso = 0; indice_peso < p->size; indice_peso++)
    {
        // Determino il gradiente del peso
        double gradiente_peso = gradiente_percettrone * input[indice_peso];

        // Modifico il peso
        p->pesi[indice_peso] = p->pesi[indice_peso] + (gradiente_peso * LRE);
    }

    p->bias = p->bias + (gradiente_percettrone * LRE);
    //printf("\n");
}

/*
    ################# IMPORT EXPORT ################################
*/


void salvaReteNeurale(const char *filename, ReteNeurale *rete)
{
    FILE *file = fopen(filename, "wb");
    if (!file)
    {
        perror("Errore nell'apertura del file");
        exit(EXIT_FAILURE);
    }

    // Scrivi il numero di layer
    fwrite(&rete->size, sizeof(int), 1, file);

    // Scrivi ogni layer
    for (int i = 0; i < rete->size; i++)
    {
        Layer *layer = &rete->layers[i];
        fwrite(&layer->size, sizeof(int), 1, file);

        // Scrivi ogni percettrone nel layer
        for (int j = 0; j < layer->size; j++)
        {
            Percettrone *perc = &layer->percettroni[j];
            fwrite(&perc->size, sizeof(int), 1, file);
            fwrite(perc->pesi, sizeof(double), perc->size, file);
            fwrite(&perc->bias, sizeof(double), 1, file);
        }
    }

    fclose(file);
}

ReteNeurale *caricaReteNeurale(const char *filename)
{
    FILE *file = fopen(filename, "rb");

    if (!file)
    {
        perror("Errore nell'apertura del file");
        return NULL;
    }

    ReteNeurale *rete = malloc(sizeof(ReteNeurale));
    if (!rete)
    {
        perror("Errore nell'allocazione della memoria");
        return NULL;
    }

    // Leggi il numero di layer
    fread(&rete->size, sizeof(int), 1, file);
    rete->layers = malloc(rete->size * sizeof(Layer));
    if (!rete->layers)
    {
        perror("Errore nell'allocazione della memoria");
        return NULL;
    }

    // Leggi ogni layer
    for (int i = 0; i < rete->size; i++)
    {
        Layer *layer = &rete->layers[i];
        fread(&layer->size, sizeof(int), 1, file);
        layer->percettroni = malloc(layer->size * sizeof(Percettrone));
        if (!layer->percettroni)
        {
            perror("Errore nell'allocazione della memoria");
            return NULL;
        }

        // Leggi ogni percettrone nel layer
        for (int j = 0; j < layer->size; j++)
        {
            Percettrone *perc = &layer->percettroni[j];
            fread(&perc->size, sizeof(int), 1, file);
            perc->pesi = malloc(perc->size * sizeof(double));
            if (!perc->pesi)
            {
                perror("Errore nell'allocazione della memoria");
                return NULL;
            }
            fread(perc->pesi, sizeof(double), perc->size, file);
            fread(&perc->bias, sizeof(double), 1, file);
        }
    }

    fclose(file);
    return rete;
}



void debug(double valore)
{
    printf("valore: %f\n", valore);
}

void debug_matrice(ReteNeurale *rete_neurale, double **matrice)
{
    for (int count = 0; count < rete_neurale->size; count++)
    {
        for (int count_2 = 0; count_2 < rete_neurale->layers[count].size; count_2++)
        {
            printf("[%d][%d]: %f\t", count, count_2, matrice[count][count_2]);
        }
        printf("\n");
    }
}

void debug_vettore(double *vettore, int size)
{
    for (int count_2 = 0; count_2 < size; count_2++)
    {
        printf("[%d]: %f\t", count_2, vettore[count_2]);
    }
    printf("\n");
}