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

char *file_pesi = "rete_mnist.bin";

#include "mnist/mnist_manager.h"

// #include "cifar-10/cifar10_manager.h";
// char *file_immagini = "cifar-10/data_batch_1.bin";
//  char *file_immagini = "cifar-10/data_batch_2.bin";
//  char *file_immagini = "cifar-10/data_batch_3.bin";
//  char *file_immagini = "cifar-10/data_batch_4.bin";
//  char *file_immagini = "cifar-10/data_batch_5.bin";
// char *file_immagini = "cifar-10/test_batch.bin";

// #include "xor_manager.h"

typedef unsigned char byte;

double LRE = 0.1;
double soglia_sigmoide = 0.5;
#define TOLLERANZA 97

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

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

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

/*
 * Genera un valore casuale nell'intervallo [-1, 1] per inizializzare pesi e bias.
 * Usa rand() modulo 101 (0-100) convertito in decimale (0.00-1.00),
 * poi scalato a [-1, 1] moltiplicando per 2 e sottraendo 1.
 */
double randomico();

/*
 * Inizializza un percettrone con n_pesi connessioni in entrata.
 * Assegna valori casuali a tutti i pesi e al bias usando randomico().
 * La dimensione n_pesi dipende dal layer precedente (o dagli input).
 */
Percettrone inzializza_percettrone(int);

/*
 * Crea una rete neurale con architettura dinamica:
 * - numero_input: dimensione dei dati in ingresso (es. 784 per MNIST 28x28)
 * - numero_layers: quanti strati nascosti + strato di output
 * - numero_percettroni_iniziali: percettroni nel primo layer
 * - numero_percettroni_finali: percettroni nell'output (10 per MNIST, 1 per XOR)
 *
 * I layer intermedi hanno numero decrescente di percettroni calcolato come:
 * numero_percettroni = numero_percettroni_iniziali * (layer_rimanenti / layer_totali)
 * Questo crea una piramide: molti neuroni in ingresso, pochi in uscita.
 */
ReteNeurale inizializza_rete_neurale(int, int, int, int);

/*
 * Crea un layer con n_percettroni, ognuno con n_pesi connessioni.
 * Ogni percettrone viene inizializzato indipendentemente con pesi casuali.
 */
Layer inizializza_layer(int, int);

/*
 * Funzione di attivazione sigmoide: trasforma la somma pesata in valore tra 0 e 1.
 * Calcola: sommatoria(input_i * peso_i) + bias, poi applica 1/(1+e^(-somma))
 * La sigmoide introduce non-linearità e "schiaccia" i valori estremi.
 */
double sigmoide(Percettrone p, double *);

/*
 * Derivata della sigmoide: necessaria per il calcolo del gradiente.
 * Formula: sigmoide(x) * (1 - sigmoide(x))
 * Indica quanto velocemente cambia l'output rispetto all'input.
 */
double derivata_sigmoide(double);

/*
 * Softmax: normalizza l'output dell'ultimo layer in probabilità (somma = 1).
 * Passaggi:
 * 1. Trova il valore massimo per stabilità numerica (evita overflow)
 * 2. Calcola e^(x_i - max) per ogni output
 * 3. Divide ogni valore per la somma totale
 * Utile per classificazione multi-classe: il percettrone con valore più alto
 * rappresenta la classe predetta con la probabilità più alta.
 */
void softmax(double *, int);

/*
 * Forward propagation: calcola gli output di tutti i layer partendo dall'input.
 *
 * Algoritmo:
 * 1. Converte i byte input (0-255) in double normalizzati (0.0-1.0) dividendo per 255
 * 2. Per il primo layer: applica sigmoide a ogni percettrone usando gli input
 * 3. Per i layer successivi: usa gli output del layer precedente come input
 *
 * Restituisce una matrice sigmoidi[layer][percettrone] con tutti i valori intermedi.
 * Questi valori servono poi per il calcolo dei gradienti (backpropagation).
 */
double **elabora_sigmoidi(ReteNeurale *, Istanza);

/*
 * Calcola i gradienti dell'errore per ogni percettrone usando backpropagation.
 *
 * Passaggi:
 * 1. GRADIENTI OUTPUT (ultimo layer):
 *    - Se il percettrone corrisponde alla classe corretta: errore = 1 - output
 *    - Altrimenti: errore = 0 - output
 *    - Moltiplica errore per derivata_sigmoide(output) per ottenere il gradiente
 *
 * 2. GRADIENTI NASCOSTI (layer precedenti):
 *    - Usa discesa_gradiente() per propagare l'errore all'indietro
 *    - Ogni neurone riceve la somma pesata dei gradienti dei neuroni successivi
 *
 * Il gradiente indica quanto e in che direzione modificare i pesi per ridurre l'errore.
 */
double **elabora_gradienti(ReteNeurale *, byte, double **);

/*
 * Propaga i gradienti dagli output verso i layer precedenti (backpropagation).
 *
 * Per ogni layer (dall'ultimo al primo):
 * - Per ogni percettrone del layer corrente:
 *   1. Calcola derivata_sigmoide dell'output del percettrone
 *   2. Calcola gradiente_disceso: somma pesata dei gradienti del layer successivo
 *      moltiplicati per i pesi delle connessioni
 *   3. Il gradiente finale = gradiente_disceso * derivata_sigmoide
 *
 * Questo è l'algoritmo della "catena" di derivate (chain rule) del calcolo differenziale.
 */
void discesa_gradiente(ReteNeurale *, double **, double **);

/*
 * Calcola quanto l'errore del layer success influenza questo percettrone.
 * Somma: gradiente[j] * peso_connessione[i][j] per ogni percettrone j del layer dopo.
 * Questo è il "peso" dell'errore che deve essere corretto a questo livello.
 */
double calcola_gradiente_disceso(ReteNeurale *, int, int, double **);
// double trova_gradiente_errore_softmax(double *, int , int);

/*
 * Aggiorna tutti i pesi della rete usando i gradienti calcolati (discesa del gradiente).
 *
 * Formula di aggiornamento: peso_nuovo = peso_vecchio + (gradiente * learning_rate)
 *
 * Per l'ultimo layer:
 * - gradiente_peso = gradiente_percettrone * output_layer_precedente
 * - Aggiorna pesi e bias
 *
 * Per i layer nascosti:
 * - Usa correggi_pesi_percettrone_double() per layer interni (input da double)
 * - Usa correggi_pesi_percettrone_byte() per primo layer (input da byte originale)
 *
 * Il learning rate (LRE) controlla la "velocità" di apprendimento.
 */
void aggiorna_pesi(ReteNeurale *, double **, double **, Istanza);

/*
 * Corregge i pesi di un percettrone usando input double (da layer precedenti).
 * Per ogni peso:
 * - gradiente_peso = gradiente_percettrone * input_corrispondente
 * - peso += gradiente_peso * learning_rate
 * - bias += gradiente_percettrone * learning_rate
 */
void correggi_pesi_percettrone_double(Percettrone *, int, double **, double);

/*
 * Corregge i pesi del primo layer usando input byte (dati originali).
 * Simile a correggi_pesi_percettrone_double ma converte i byte in double.
 * Necessario perché il primo layer riceve input non normalizzati dal dataset.
 */
void correggi_pesi_percettrone_byte(Percettrone *, Istanza, double, int);

/*
 * Predizione binaria: restituisce 1 se valore >= soglia_sigmoide (0.5), altrimenti 0.
 * Usata per classificazione binaria (es. XOR o singola categoria MNIST).
 */
/*
 * Predizione multi-classe: applica softmax all'output e restituisce l'indice
 * del percettrone con valore più alto (la classe predetta).
 * Usata per MNIST e CIFAR-10 (10 classi).
 */
int previsione_softmax(double *, int);

/*
 * Salva la rete neurale su file binario includendo:
 * - Numero di layer
 * - Per ogni layer: numero di percettroni
 * - Per ogni percettrone: numero di pesi, array pesi, bias
 *
 * Il formato binario permette di ricaricare il modello addestrato
 * senza dover ripetere l'allenamento.
 */
void salvaReteNeurale(const char *, ReteNeurale *);

/*
 * Carica una rete neurale da file binario.
 * Legge la struttura e alloca dinamicamente tutta la memoria necessaria.
 *
 * Restituisce NULL se il file non esiste o c'è un errore di allocazione.
 */
ReteNeurale *caricaReteNeurale(const char *);

/*
 * Stampa su stdout tutti i pesi e i bias della rete.
 * Utile per debugging e ispezione dello stato del modello.
 * Formato: Livello -> Percettrone -> Peso[i]: valore
 */
void stampa_pesi_rete(ReteNeurale *);

/*############# INIZIALIZZAZIONI #########################*/

double randomico()
{
    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(int numero_input, int numero_layers, int numero_percettroni_iniziali, int numero_percettroni_finali)
{
    srand(time(NULL));
    ReteNeurale r;
    r.layers = (Layer *)malloc(sizeof(Layer) * numero_layers);
    r.size = numero_layers;
    int somma_parametri = 0;

    for (int livello = 0; livello < numero_layers; livello++)
    {
        int numero_percettroni_livello = 1;

        if (livello == numero_layers - 1)
            numero_percettroni_livello = numero_percettroni_finali;
        else
        {
            double frazione = (double)(numero_layers - livello) / (double)numero_layers;
            numero_percettroni_livello = (int)((double)(numero_percettroni_iniziali * frazione));
        }

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

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

    printf("Rete neurale da %d parametri\n", somma_parametri);

    return r;
}

/*
    ################# PREVISIONI ################################
*/

double **elabora_gradienti(ReteNeurale *rete_neurale, byte output_corretto, double **sigmoidi)
{
    double **gradienti = (double **)malloc(sizeof(double *) * rete_neurale->size);

    for (int indice_layer = 0; indice_layer < rete_neurale->size; indice_layer++)
    {
        gradienti[indice_layer] = (double *)malloc(sizeof(double) * rete_neurale->layers[indice_layer].size);
    }

    for (int indice_percettrone = 0; indice_percettrone < rete_neurale->layers[rete_neurale->size - 1].size; indice_percettrone++)
    {
        double gradiente_errore = 0.0;
        if (indice_percettrone == output_corretto)
        {
            gradiente_errore = 1 - sigmoidi[rete_neurale->size - 1][indice_percettrone];
        }
        else
        {
            gradiente_errore = 0 - sigmoidi[rete_neurale->size - 1][indice_percettrone];
        }

        gradienti[rete_neurale->size - 1][indice_percettrone] = gradiente_errore * derivata_sigmoide(sigmoidi[rete_neurale->size - 1][indice_percettrone]);
    }

    discesa_gradiente(rete_neurale, sigmoidi, gradienti);

    return gradienti;
}

double sigmoide(Percettrone p, double *valori)
{
    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;

    // Sigmoide
    double risultato = 1.0 / (1.0 + exp(-sommatoria));
    // ReLU
    // double risultato = sommatoria > 0 ? sommatoria : 0.0f;

    // printf(" sommatoria %f -> %f\n",sommatoria, risultato);

    return risultato;
}

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

void softmax(double *input, int size)
{
    float max = input[0];

    for (int i = 1; i < size; i++)
        if (input[i] > max)
            max = input[i];

    float sum = 0.0f;
    for (int i = 0; i < size; i++)
    {
        input[i] = expf(input[i] - max);
        sum += input[i];
    }

    for (int i = 0; i < size; i++)
        input[i] /= sum;
}

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

int previsione_softmax(double *livello_percettroni, int size)
{

    softmax(livello_percettroni, size);
    int max = 0;

    for (int i = 1; i < size; i++)
    {
        if (livello_percettroni[i] > livello_percettroni[max])
            max = i;
    }

    return max;
}

void discesa_gradiente(ReteNeurale *rete, double **sigmoidi, double **gradienti)
{
    for (int indice_layer = rete->size - 2; indice_layer >= 0; indice_layer--)
    {
        for (int indice_percettrone = 0; indice_percettrone < rete->layers[indice_layer].size; indice_percettrone++)
        {
            double derivata_attivazione = derivata_sigmoide(sigmoidi[indice_layer][indice_percettrone]);
            double gradiente_disceso = calcola_gradiente_disceso(rete, indice_layer + 1, indice_percettrone, gradienti);

            gradienti[indice_layer][indice_percettrone] = gradiente_disceso * derivata_attivazione;
        }
    }
}

double calcola_gradiente_disceso(ReteNeurale *rete, int livello, int indice_peso, double **gradienti)
{
    double sommatoria = 0.0;
    for (int indice_percettrone = 0; indice_percettrone < rete->layers[livello].size; indice_percettrone++)
    {
        sommatoria += (gradienti[livello][indice_percettrone] * rete->layers[livello].percettroni[indice_percettrone].pesi[indice_peso]);
    }

    return sommatoria;
}

double **elabora_sigmoidi(ReteNeurale *rete, Istanza istanza)
{
    double **sigmoidi = (double **)malloc(sizeof(double *) * rete->size);
    double *inputs = (double *)malloc(sizeof(double) * N_INPUTS);
    for (int i = 0; i < N_INPUTS; i++)
    {
        inputs[i] = (double)istanza.dati[i] / 255.0;
    }

    sigmoidi[0] = (double *)malloc(sizeof(double) * rete->layers[0].size);
    for (int indice_percettrone = 0; indice_percettrone < rete->layers[0].size; indice_percettrone++)
    {
        sigmoidi[0][indice_percettrone] = sigmoide(rete->layers[0].percettroni[indice_percettrone], inputs);
    }

    free(inputs);

    for (int indice_layer = 1; indice_layer < rete->size; indice_layer++)
    {
        sigmoidi[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++)
        {
            sigmoidi[indice_layer][indice_percettrone] = sigmoide(rete->layers[indice_layer].percettroni[indice_percettrone], sigmoidi[indice_layer - 1]);
        }
    }

    return sigmoidi;
}

/*
    ################# CORREZIONI ################################
*/

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

    for (int indice_percettrone = 0; indice_percettrone < rete_neurale->layers[rete_neurale->size - 1].size; indice_percettrone++)
    {

        for (int indice_peso = 0; indice_peso < rete_neurale->layers[rete_neurale->size - 1].percettroni[indice_percettrone].size; indice_peso++)
        {
            double gradiente_peso = gradienti[rete_neurale->size - 1][indice_percettrone] * sigmoidi[rete_neurale->size - 2][indice_peso];
            rete_neurale->layers[rete_neurale->size - 1].percettroni[indice_percettrone].pesi[indice_peso] += gradiente_peso * LRE;
        }
        rete_neurale->layers[rete_neurale->size - 1].percettroni[indice_percettrone].bias += gradienti[rete_neurale->size - 1][indice_percettrone] * LRE;
    }

    for (int indice_layer = rete_neurale->size - 2; indice_layer >= 0; indice_layer--)
    {
        for (int indice_percettrone = 0; indice_percettrone < rete_neurale->layers[indice_layer].size; indice_percettrone++)
        {
            if (indice_layer != 0)
            {
                correggi_pesi_percettrone_double(&rete_neurale->layers[indice_layer].percettroni[indice_percettrone], indice_layer, sigmoidi, gradienti[indice_layer][indice_percettrone]);
            }
            else
            {
                correggi_pesi_percettrone_byte(&rete_neurale->layers[0].percettroni[indice_percettrone], istanza, gradienti[0][indice_percettrone], indice_percettrone);
            }
        }
    }
}

void correggi_pesi_percettrone_double(Percettrone *p, int layer, double **input, double gradiente_percettrone)
{

    for (int indice_peso = 0; indice_peso < p->size; indice_peso++)
    {
        double gradiente_peso = gradiente_percettrone * input[layer - 1][indice_peso];
        p->pesi[indice_peso] += (gradiente_peso * LRE);
    }

    p->bias += (gradiente_percettrone * LRE);
}

void correggi_pesi_percettrone_byte(Percettrone *p, Istanza input, double gradiente_percettrone, int indice_percettrone)
{
    for (int indice_peso = 0; indice_peso < p->size; indice_peso++)
    {
        double gradiente_peso = gradiente_percettrone * (double)input.dati[indice_peso];
        p->pesi[indice_peso] += (gradiente_peso * LRE);
    }

    p->bias += (gradiente_percettrone * LRE);
}

/*
    ################### ADDESTRAMENTO #########################
*/

char addestra(ReteNeurale *rete_neurale, Dataset set)
{
    int corrette = 0;
    for (int indice_set = 0; indice_set < set.size; indice_set++)
    {
        byte output_corretto = set.istanze[indice_set].classificazione;

        double **sigmoidi = elabora_sigmoidi(rete_neurale, set.istanze[indice_set]);

        double *sigmoidi_softmax = (double *)malloc(sizeof(double) * rete_neurale->layers[rete_neurale->size - 1].size);
        memcpy(sigmoidi_softmax, sigmoidi[rete_neurale->size - 1], rete_neurale->layers[rete_neurale->size - 1].size * sizeof(double));

        double **gradienti = elabora_gradienti(rete_neurale, output_corretto, sigmoidi);

        int previsto = previsione_softmax(sigmoidi_softmax, rete_neurale->layers[rete_neurale->size - 1].size);

        if (previsto == output_corretto)
        {
            corrette++;
        }

        aggiorna_pesi(rete_neurale, sigmoidi, gradienti, set.istanze[indice_set]);

        for (int i = 0; i < rete_neurale->size; i++)
            free(gradienti[i]);
        free(gradienti);

        for (int i = 0; i < rete_neurale->size; i++)
            free(sigmoidi[i]);
        free(sigmoidi);
        free(sigmoidi_softmax);
    }

    int percentuale = (int)((corrette * 100) / set.size);

    printf("Risposte corrette: %d%%\n", percentuale);

    if (percentuale >= TOLLERANZA)
        return 1;
    else
        return 0;
}

/*
    ################# 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);
    }

    fwrite(&rete->size, sizeof(int), 1, file);

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

        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;
    }

    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;
    }

    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;
        }

        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;
}

/*
    ################ STAMPE ############################
*/
void stampa_pesi_rete(ReteNeurale *rete)
{
    for (int indice_layer = 0; indice_layer < rete->size; indice_layer++)
    {
        printf("\nLivello %d", indice_layer);
        for (int indice_percettrone = 0; indice_percettrone < rete->layers[indice_layer].size; indice_percettrone++)
        {
            printf("\n\tPercettrone %d", indice_percettrone);
            for (int indice_peso = 0; indice_peso < rete->layers[indice_layer].percettroni[indice_percettrone].size; indice_peso++)
            {
                printf("\n\t\tPeso %d: %f", indice_peso, rete->layers[indice_layer].percettroni[indice_percettrone].pesi[indice_peso]);
            }
        }
    }
    printf("\n");
}