function [cap_carico,matagenti,agente,mondo,stats, coeff_gini] = sugarscape_sostagenti()
% function [cap_carico,matagenti,agente,mondo,stats,coeff_gini] = sugarscape_sostagenti()
%
% Questa versione dello Sugarscape é caratterizzata dal fatto che gli agenti sono dotati di una
% certa età, ad ogni step invecchiano, fino a giungere alla loro età di morte. Una volta morti,
% indipendentemente dal fatto che ciò abbia avuto luogo a causa delle scorte o dell'età, vengono
% rimpiazzati da nuovi agenti.
% Gli output sono nell'ordine:
% - la capacità portante del sistema, 
% - la matrice di 0 1 dove gli 1 rappresentano celle dello sugarscape occupate da un agente,
% - la struttura dati contenneti le caratteristiche di ogni agente rimasto all'ultimo step,
% - la matrice che rappresenta sugarscape 
% - la struttura contenente tutte le statistiche
% - il coefficiente di gini.

    % Chiude le eventuali finestre grafiche aperte
    close all;

    % Richiamo la funzione che prende i parametri dall'esterno
    prm = leggi_parametri;

    % Richiamo la funzione che crea lo sugarscape e gli agenti
    [sugarscape, sugmax, agente, matagenti] = crea(prm);

    % Inizializzo il vettore stats
    stats=[];
    
    %visualizzo le condizioni iniziali
    visualizza_sugarscape(sugarscape,agente,prm,0);
    pause;

    % Ciclo che fa evolvere il sistema
    %**************************************************************************
    for step=1:prm.npassi
        step
        % Richiamo della funzione che fa muovere e raccogliere gli agenti
        [agente, matagenti, sugarscape] = muovi(sugarscape,prm,agente,matagenti);	
        % Richiamo della funzione che fa consumare gli agenti
        [agente] = consuma(agente,prm);
        % Richiamo della funzione che fa invecchiare gli agenti
        [agente] = invecchia(agente,prm);
        % Richiamo della funzione che fa apparire l'angelo della morte
        [agente,matagenti,prm] = ammazza(agente,matagenti,prm);
        % Richiamo della funzione che fa crescere lo sguarscape
        sugarscape = cresci(sugarscape);
        % Richiamo della funzione che controlla che le celle dello sugarscape non superino il limite
        % massimo di zucchero contenuto al loro interno
        sugarscape = controlla(sugarscape,sugmax);
        % Richiamo della visualizzazione che mostra lo sugarscape e gli agenti nel corso
        % dell'evoluzione
        %visualizza_sugarscape(sugarscape,agente,prm,step);
        % Richiamo della funzione che calcola ciclicamente le varie statistiche
        stats = statistiche(sugarscape,step,prm,agente,stats);
    end

    %visualizzo le condizioni finali
    visualizza_sugarscape(sugarscape,agente,prm,prm.npassi);
    pause;
    
    % Salvo il numero di agenti finali nella variabile cap_carico
    cap_carico = prm.numagenti;

    % Memorizzo tutte le statistiche finali in un file   
    scrivi_parametri(stats);

    % Istruzione per mandare in output a matlab lo sugarscape
    mondo=sugarscape;

    % Visualizzazione dei grafici delle statistiche (con il coeff. di Gini in output)
    [coeff_gini] = visualizza_statistiche(stats,prm,agente);
    
%==========================================================================
%CREAZIONE DELLE CONDIZIONI INIZIALI
%==========================================================================
function [sugarscape, sugmax, agente, matagenti] = crea(prm)
% Con questa funzione creo 
    %la matrice sugarscape(la condizione di partenza del mondo)
    %la matrice sugmax(la matrice che indica la capacitò massima delle
    %varie celle dello sugaerscape
    %i vari agenti con le loro proprietà in una struttura dati
    %la matrice matagenti che indica quali coordinate occupano gli agenti
    
    %CREAZIONE SUGARSCAPE
    %*********************
    % Creazione della matrice sugarscape (condizione di partenza del mondo)
    sugarscape=4*rand([prm.nr,prm.nc]);

    % La funzione ceil() serve per convertire valori reali in numeri interi
    % (essendo la quantita' di zucchero contenuta nelle singole celle un numero intero compreso tra 0 e 4)
    sugarscape = ceil(sugarscape);

    % Elevamento al quadrato dei coefficienti sigmax e sigmay
    sigmaxx=prm.sigmax*prm.sigmax;
    sigmayy=prm.sigmay*prm.sigmay;
    
    %CREAZIONE SUGMAX
    %******************
    % Elevo al quadrato i coefficienti sigmax e sigmay
    sigmaxx=prm.sigmax*prm.sigmax;
    sigmayy=prm.sigmay*prm.sigmay;
    %Uso l'equazione di due collinette per calcolare la capacità di ogni cella
    for i=1:prm.nr
        for j=1:prm.nc
        % Arrotondo separatamente con ceil le due colinnete  per
        % averle distaccate nello sugarscape
         b(i,j)=prm.altezza*exp(-(i-prm.x0)*(i-prm.x0)/sigmaxx)*exp(-(j-prm.y0)*(j-prm.y0)/sigmayy);
         c(i,j)=prm.altezza*exp(-(i-prm.x1)*(i-prm.x1)/sigmaxx)*exp(-(j-prm.y1)*(j-prm.y1)/sigmayy);
        end;
    b=ceil(b);
    c=ceil(c);
    sugmax=max(b,c);
    end;
    
    sugarscape = ceil(sugarscape);
    % Primo controllo: "livella" la matrice random con la capacita' massima del sistema
    sugarscape = min(sugarscape,sugmax);

    % CREAZIONE DEGLI AGENTI
    % matagenti e' una matrice di dimensione nrXnc, composta da elementi inizialmente tutti uguali a 0.
    % Quando una cella e' occupata da un agente, il suo valore diventa uno.
    matagenti = zeros(prm.nr,prm.nc);

    for i=1:prm.numagenti

        % definisco la struttura agente (eta', morte, scorte, metabolismo e vista)	
        agente(i).eta = ceil(59*rand);
        agente(i).morte = 59 + ceil(41*rand);
        agente(i).scorte = 10;
        agente(i).met = ceil(5*rand);
        agente(i).vista = ceil(6*rand);

        % Definizione delle coordinate iniziali dell'agente
        agente(i).x = ceil(prm.nc*rand);
        agente(i).y = ceil(prm.nr*rand);
        % Ri-calcolo delle coordinate, fino a quando non se ne trova una non occupata
            while (matagenti(agente(i).y,agente(i).x) ~= 0)
            agente(i).x = ceil(prm.nc*rand);
            agente(i).y = ceil(prm.nr*rand);
            end
        % Nel caso in cui le coordinate siano libere, viene aggiornata la matrice matagenti
        matagenti(agente(i).y,agente(i).x) = 1;
    end
    

%==================================================================
%MOVIMENTO e RACCOLTA AGENTI
%===================================================================
function [agente,matagenti,sugarscape] = muovi(sugarscape,prm,agente,matagenti);

    % Ciclo for che muove gli agenti (non tutti: il ciclo termina dopo prm.numagenti movimenti,
    % indipendentemente dal fatto che alcuni agenti si siano mossi piu' volte o siano rimasti fermi.
    for j=1:prm.numagenti

        % Selezione casuale dell'agente che si dovra' muovere  
        selezionato = ceil(prm.numagenti*rand);

        % La matrice poss_destinazioni conterra' le possibili celle in cui l'agente potra' spostarsi
        % (matrice tre colonne: la prima contiene la quantita' di zucchero della cella di destinazione, la seconda la sua coordinata x e la terza la coordinata sull'asse y)

        % Azzera/inizializza la matrice poss_destinazioni
        poss_destinazioni = [];

         % CICLO CHE FA GUARDARE L'AGENTE NELLE 4 DIREZIONI
         %*************************************************
        for i=1:agente(selezionato).vista
            % DESTINAZIONI AL DI SOPRA DELL'AGENTE
            % Controllo che la destinazione esista
            if agente(selezionato).y+i <= prm.nr
                % Controllo che la destinazione sia vuota
                if matagenti(agente(selezionato).y+i,agente(selezionato).x) == 0
                    % Se la destinazione esiste ed e' vuota, allora scrive una nuova riga nella matrice possibili_destinazioni
                    [nr_poss_destinazioni,var2] = size(poss_destinazioni);
                    poss_destinazioni(nr_poss_destinazioni+1,1)=sugarscape(agente(selezionato).y+i,agente(selezionato).x)-0.01*i+0.00001*rand;
                    poss_destinazioni(nr_poss_destinazioni+1,2)=agente(selezionato).x;
                    poss_destinazioni(nr_poss_destinazioni+1,3)=agente(selezionato).y+i;
                end
            end
            % DESTINAZIONI AL DI SOTTO DELL'AGENTE
            % Controllo che la destinazione esista
            if agente(selezionato).y-i > 0
                % Controllo che la destinazione sia vuota  
                if matagenti(agente(selezionato).y-i,agente(selezionato).x) == 0
                    % Se la destinazione esiste ed e' vuota, allora scrive una nuova riga nella matrice possibili_destinazioni
                    [nr_poss_destinazioni,var2] = size(poss_destinazioni);
                    poss_destinazioni(nr_poss_destinazioni+1,1) = sugarscape(agente(selezionato).y-i,(agente(selezionato).x))-0.01*i+0.00001*rand;
                    poss_destinazioni(nr_poss_destinazioni+1,2) = agente(selezionato).x;
                    poss_destinazioni(nr_poss_destinazioni+1,3) = agente(selezionato).y-i;
                end
            end            
            % DESTINAZIONI A DESTRA DELL'AGENTE
            % Controllo che la destinazione esista
            if agente(selezionato).x+i <= prm.nc
                % Controllo che la destinazione sia vuota 
                if matagenti(agente(selezionato).y,agente(selezionato).x+i) == 0
                    % Se la destinazione esiste ed e' vuota, allora scrive una nuova riga nella matrice possibili_destinazioni
                    [nr_poss_destinazioni,var2] = size(poss_destinazioni);
                    poss_destinazioni(nr_poss_destinazioni+1,1) = sugarscape(agente(selezionato).y,agente(selezionato).x+i)-0.01*i+0.00001*rand;
                    poss_destinazioni(nr_poss_destinazioni+1,2) = agente(selezionato).x+i;
                    poss_destinazioni(nr_poss_destinazioni+1,3) = agente(selezionato).y;
                end
            end	    
            % DESTINAZIONI A SINISTRA DELL'AGENTE
            % Controllo che la destinazione esista
            if agente(selezionato).x-i > 0
                % Controllo che la destinazione sia vuota 
                if matagenti(agente(selezionato).y,agente(selezionato).x-i) == 0
                    % Se la destinazione esiste ed e' vuota, allora scrive una nuova riga nella matrice possibili_destinazioni
                    [nr_poss_destinazioni,var2] = size(poss_destinazioni);
                    poss_destinazioni(nr_poss_destinazioni+1,1) = sugarscape(agente(selezionato).y,agente(selezionato).x-i)-0.01*i+0.00001*rand;
                    poss_destinazioni(nr_poss_destinazioni+1,2) = agente(selezionato).x-i;
                    poss_destinazioni(nr_poss_destinazioni+1,3) = agente(selezionato).y;
                end
            end
        end
        
        %SCELTA DESTINAZIONE
        %********************
        % Prima di muovere controlliamo che la matrice delle destinazioni possibili non sia vuota
        if (length(poss_destinazioni) ~= 0)
        % Ora abbiamo a disposizione il vettore delle destinazioni possibili

            % Scegliamo in quale cella muoverci
            [a,b] = max(poss_destinazioni(:,1));
            destinazione_scelta.x = poss_destinazioni(b,2);
            destinazione_scelta.y = poss_destinazioni(b,3);

            % Muoviamo l'agente nella casella scelta
            matagenti(agente(selezionato).y,agente(selezionato).x) = 0;
            matagenti(destinazione_scelta.y,destinazione_scelta.x) = 1;
            agente(selezionato).x = destinazione_scelta.x; 
            agente(selezionato).y = destinazione_scelta.y;
               
            %RACCOLTA ZUCCHERO
            % Facciamo raccogliere all'agente lo zucchero presenta nella cella in cui
            % si e'appena spostato
            agente(selezionato).scorte = agente(selezionato).scorte + sugarscape(agente(selezionato).y,agente(selezionato).x);
            sugarscape(agente(selezionato).y,agente(selezionato).x) = 0;
        end
    end

%==========================================================================
%CONSUMO ZUCCHERO
%==========================================================================
% Funzione che fa consumare agli agenti la quantita' di scorte corrispondenti al loro metabolismo
function [agente] = consuma(agente,prm)

    for i=1:prm.numagenti
        agente(i).scorte = agente(i).scorte - agente(i).met;
    end

%==========================================================================
%INVECCHIAMENTO AGENTI
%==========================================================================
function [agente] = invecchia(agente,prm)

    for i=1:prm.numagenti
        agente(i).eta = agente(i).eta + 1;
    end

%==========================================================================
%ANGELO DELLA MORTE
%==========================================================================
% Funzione che uccide gli agenti che non hanno sufficienti scorte per sopravvivere
% o che hanno raggiunto l'eta' di morte e li sostituisce con agenti nuovi
function [agente,matagenti,prm] = ammazza(agente,matagenti,prm)

for i=prm.numagenti:-1:1
	% Controlla se l'agente ha finito le scorte o ha raggiunto l'eta' di morte
	if ((agente(i).scorte <= 0) || (agente(i).eta >= agente(i).morte))
		% Nel caso in cui l'agente muore, viene sostituito da uno con eta' uguale a 0...
		agente(i).eta = 0;
		% ... e con gli altri parametri random
		agente(i).morte = 59 + ceil(41*rand);
		agente(i).scorte = 10;
		agente(i).met = ceil(5*rand);
		agente(i).vista = ceil(6*rand);
	end
end

%==========================================================================
%CRESCITA SUGARSCAPE
%==========================================================================
function a = cresci(a)

    % Aumenta di 1 il livello di zucchero di ogni cella all'interno della matrice A (sugarscape)
    a=a+1;

%==========================================================================
%LIVELLAMENTO SUGARSCAPE
%==========================================================================
function a = controlla(a,b)

    % se a è più piccolo della capacità di carico della cella, allora ok...
    % se a è più grande, dobbiamo ridurlo alla max capacità del terreno
    % possiamo quindi fare un controllo sul minimo delle matrici
    a = min(a,b);
    
%==========================================================================
% STATISTICHE FINALI
%==========================================================================
function stats = statistiche(sugarscape,step,prm,agente,stats);

    % Memorizzo il numero di agenti nella struttura stats
    stats(step).numagenti = prm.numagenti;

    % Nel caso in cui vi siano ancora agenti in vita, calcolo le varie statistiche
    if (prm.numagenti ~= 0)

        % Inizializzo le somme ed il vettore scorte
        somma_met = 0;
        somma_vis = 0;
        somma_scorte = 0;
        vettore_scorte = [];

        % Calcolo le somme dei valori di metabolismo, vista e scorte
        for i=1:prm.numagenti
            somma_met = somma_met + agente(i).met;
            somma_vis = somma_vis + agente(i).vista;
            vettore_scorte(length(vettore_scorte)+1) = agente(i).scorte;
            %memorizzo metabolismo e vista in un vettore per poter poi
            %calcolare la deviazioni standard dello step
            vettore_vista(i)=agente(i).vista;
            vettore_met(i)=agente(i).met; 
        end

        % Calcolo le varie statistiche
        stats(step).met_medio = somma_met/prm.numagenti;
        stats(step).vista_media = somma_vis/prm.numagenti;
        stats(step).scorta_media = sum(vettore_scorte)/prm.numagenti;
        stats(step).scorta_min = min(vettore_scorte);
        stats(step).scorta_max = max(vettore_scorte);
        stats(step).scorta_mediana = median(vettore_scorte);
        %calcolo le deviazioni standard
        met_diff=vettore_met-stats(step).met_medio;
        vista_diff=vettore_vista-stats(step).vista_media;
        square_met_diff=met_diff.^2;
        square_vista_diff=vista_diff.^2;
        stats(step).met_sd=sum(square_met_diff)/prm.numagenti;
        stats(step).vista_sd=sum(square_vista_diff)/prm.numagenti;

    else
        % Nel caso in cui non sia rimasto in vita alcun agente, inseriamo 0 in tutte le statistiche
        stats(step).met_medio = 0;
        stats(step).vista_media = 0;
        stats(step).scorta_media = 0;
        stats(step).scorta_min = 0;
        stats(step).scorta_max = 0;
        stats(step).scorta_mediana = 0;
        stats(step).met_sd=0;
        stats(step).vista_sd=0;

    end

%==========================================================================
%VISUALIZZAZIONE RISULTATI
%==========================================================================
% Funzione che visualizza i vari grafici durante l'evoluzione del sistema
%**************************************************************************
function visualizza_sugarscape(sugarscape,agente,prm,step)
   
    % Disegno dello sugarscape
    figure(1);
    %se è l'utlimo step la metto in un'altra finestra
    if(step==prm.npassi)
        figure(prm.npassi);
    end
    imagesc(sugarscape);
    hold on;
    eval(['title (''sugarscape step ' ,int2str(step),''')']);

    % Disegno degli agenti
    for i=1:prm.numagenti
        plot(agente(i).x,agente(i).y,'yd');
    end
    pause(0.15);

% Funzione che visualizza le varie statistiche
%**************************************************************************
function [coeff_gini] = visualizza_statistiche(stats,prm,agente)

    % Inizializzo i vari vettori
    vettore_numagenti = [];
    vettore_metmedio = [];
    vettore_vista_media = [];
    vettore_scorta_min = [];
    vettore_scorta_media = [];
    vettore_scorta_max = [];
    vettore_scorta_mediana = [];
    vettore_vista_sd = [];
    vettore_met_sd = [];

    % Converto la struttura in vettori
    for i=1:prm.npassi
        vettore_numagenti(length(vettore_numagenti)+1) = stats(i).numagenti;
        vettore_metmedio(length(vettore_metmedio)+1) = stats(i).met_medio;
        vettore_vista_media(length(vettore_vista_media)+1) = stats(i).vista_media;
        vettore_vista_sd(length(vettore_vista_sd)+1) = stats(i).vista_sd;
        vettore_met_sd(length(vettore_met_sd)+1) = stats(i).met_sd;
        vettore_scorta_min(length(vettore_scorta_min)+1) = stats(i).scorta_min;
        vettore_scorta_media(length(vettore_scorta_media)+1) = stats(i).scorta_media;
        vettore_scorta_max(length(vettore_scorta_max)+1) = stats(i).scorta_max;
        vettore_scorta_mediana(length(vettore_scorta_mediana)+1) = stats(i).scorta_mediana;
    end

    % Visualizzo il grafico con l'andamento della popolazione degli agenti
    figure(2);
    hold on;
    title('Andamento del numero di agenti nel corso della evoluzione');
    plot(vettore_numagenti);
    plot(vettore_numagenti,'*');

    % Visualizzo il grafico con l'andamento delle medie del metabolismo e della vista degli agenti nel tempo
    figure(3);
    hold on;
    title('Andamento delle medie di metabolismo e vista degli agenti');
    plot(vettore_metmedio,'r');
    plot(vettore_vista_media,'b');
    plot(vettore_metmedio,'r*');
    plot(vettore_vista_media,'b*');
    legend('Metabolismo medio','Vista media');
    
    % Visualizzo il grafico con l'andamento delle deviazione standard del
    % metabolismo medio
    figure(7);
    hold on;
    title('Andamento del metabolismo medio e della sua deviazione standard');
    errorbar(vettore_metmedio, vettore_met_sd, 'r');

    % Visualizzo il grafico con l'andamento delle deviazione standard della
    % vista media
    figure(8);
    hold on;
    title('Andamento della vista media e della sua deviazione standard');
    errorbar(vettore_vista_media,vettore_vista_sd, 'b');

    % Visualizzo il grafico con l'andamento minimo, medio e massimo della ricchezza e con la mediana
    figure(4);
    hold on;
    title('Andamento della ricchezza minima, media, mediana e massima degli agenti');
    plot(vettore_scorta_min,'r');
    plot(vettore_scorta_media,'b');
    plot(vettore_scorta_max,'g');
    plot(vettore_scorta_mediana','k');
    legend('Ricchezza minima', 'Ricchezza media', 'Ricchezza massima', 'Mediana');
    plot(vettore_scorta_min,'r*');
    plot(vettore_scorta_media,'bo');
    plot(vettore_scorta_max,'gp');
    plot(vettore_scorta_mediana','k+');

    % metto nel vettore scorta_fin le scorte di tutti gli agenti all'ultimo step
    % per visualizzarne la distribuzione
    vettore_scorta_fin=[];
    for k=1:prm.numagenti
        vettore_scorta_fin=[vettore_scorta_fin, agente(k).scorte];
    end    

    % Visualizza l'istogramma della ricchezza 
    figure(5);
    hold on;
    title('Distribuzione finale della ricchezza');
    hist(vettore_scorta_fin);

    % Calcolo del coefficiente di Gini
    %**********************************

    % Controlla che siano rimasti agenti nel sistema
    if (prm.numagenti > 1)
        % Calcolo della retta di uguaglianza
        for i=1:prm.numagenti
            retta_uguaglianza(i) = vettore_scorta_media(prm.npassi);
        end

        % Calcolo dell'area sottesa alla retta di uguaglianza
        area_retta_uguaglianza = cumsum(retta_uguaglianza);

        % Calcolo dell'area sottesa alla curva di Lorenz
        area_curva_lorenz = cumsum(sort(vettore_scorta_fin));

        % Calcolo dell'area compresa tra la retta di uguaglianza e la curva di Lorenz
        area_compresa = area_retta_uguaglianza - area_curva_lorenz;

        % Calcolo del coefficiente di Gini
        coeff_gini = area_compresa/area_retta_uguaglianza;

        % Disegno del grafico
        figure(6);
        title('Retta di uguaglianza del reddito e curva di Lorenz');
        hold on;
        plot(area_retta_uguaglianza,'r');
        plot(area_curva_lorenz);
        legend('Retta di uguaglianza del reddito', 'Curva di Lorenz');
    else
        % Nel caso in cui non siano rimasti agenti nel sistema imposta il Coefficiente di Gini = 0
        coeff_gini = 0;
    end

%========================
%INPUT
%========================
% Funzione che va a prendere i parametri del programma da un file esterno
function prm = leggi_parametri

    % Apre, in sola lettura, un file di dati esterno
    fid=fopen('input_sugarscape.dat','r');

    % Richiama i parametri dal file di testo che li contiene
    prm.nr = fscanf(fid,'%g',1);
    prm.nc = fscanf(fid,'%g',1);
    prm.npassi = fscanf(fid,'%g',1);
    prm.numagenti = fscanf(fid,'%g',1);
    prm.x0 = fscanf(fid,'%g',1);
    prm.y0 = fscanf(fid,'%g',1);
    prm.x1 = fscanf(fid,'%g',1);
    prm.y1 = fscanf(fid,'%g',1);
    prm.sigmax = fscanf(fid,'%g',1);
    prm.sigmay = fscanf(fid,'%g',1);
    prm.altezza = fscanf(fid,'%g',1);

    % Chiude il file di dati esterno aperto in precedenza
    fclose(fid);

%========================
%OUTPUT
%========================
%Funzione che scrive le statisctiche finali in un file
function scrivi_parametri(stats)

    % Apre, in scrittura, un file di dati esterno
    fid=fopen('output_sugarscape.dat','w');

    % Scrive l'intestazione del file
    fprintf(fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\r\n','step','met_medio','vista_media','scorta_media','scorta_min','scorta_mediana','scorta_max','numagenti');

    % Ciclo che scorre la struttura statistiche e ne memorizza il contenuto all'interno del file
    for i=1:length(stats)
        fprintf(fid, '%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\r\n',i,stats(i).met_medio,stats(i).vista_media,stats(i).scorta_media,stats(i).scorta_min,stats(i).scorta_mediana,stats(i).scorta_max,stats(i).numagenti);
    end