function [cap_carico,matagenti,agente,mondo,stats] = sugarscape()
% function [cap_carico,matagenti,agente,mondo,stats] = sugarscape()
%--------------------------------------------------------------------------
% Versione base dello Sugarscape: gli agenti non invecchiano e
% non si riproducono
% Gli output sono nell'ordine:
% - la cpacità 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

% Chiude le eventuali finestre grafiche aperte
close all;
% Richiamo la funzione che prende i parametri di input dall'esterno
prm = leggi_parametri;

% Richiamo la funzione che crea tutte le condizioni iniziali
[sugarscape, sugmax, agente, matagenti] = crea(prm);

%visualizzo le condizioni iniziali
visualizza_sugarscape(sugarscape,agente,prm,0);
pause;

% Inizializzo il vettore stats che conterrà le statistiche finali del sistema
stats=[];

mondo=sugarscape;

% CICLO CHE FA EVOLVERE IL SISTEMA
%*********************************
% Ad ogni passo (dato in input) vengono mossi prm.numagenti agenti
% Questi, dopo aver scelto dove muoversi, raccolgono lo zucchero presente nella cella in cui si muovono,
% poi consumano il proprio metabolismo e successivamente compare l'angelo della morte che li elimina se
% rimasti senza scorte. Lo sugarscape in seguito cresce e si controlla che le celle non superino la
% capacita' massima consentita

for step=1:prm.npassi
    if (prm.numagenti~=0) %il sistema evolve finchè rimangono agenti
        % Richiamo della funzione che muove e fa raccogliere gli agenti
        [agente, matagenti,sugarscape] = muovi(sugarscape,prm,agente,matagenti);
        % Richiamo della funzione che fa "raccogliere" gli agenti
        % Richiamo della funzione che fa "consumare" gli agenti
        [agente] = consuma(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);
        % La funzione "controlla" avà in uscita la matrice sugarscape corretta
        sugarscape = controlla(sugarscape,sugmax);
        % La funzione "visualizza" mostra lo sugarscape con sopra gli agenti
        visualizza_sugarscape(sugarscape,agente,prm,step);
        % Calcolo ciclico delle statistiche
        if (prm.numagenti~=0) %ricontrollo che gli agenti non siano nulli dopo angelo della morte
            stats = statistiche(sugarscape,step,prm,agente,stats);
        end
     else 
       break;
    end
end

% 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);

% Visualizzazione dei grafici delle statistiche
visualizza_statistiche(stats,prm,agente);

%==========================================================================
%CREAZIONE 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
%********************
sugarscape=4*rand([prm.nr,prm.nc]);

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

% La funzione ceil() serve per convertire i valori reali di b in numeri interi
% (essendo la quantita' delle singole celle un numero intero tra 0 e 4)
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 nrxnc composta da elementi inizialmente tutti uguali a zero
% quando una cella e' occupata, il suo valore diventa uno
matagenti = zeros(prm.nr,prm.nc);

for i=1:prm.numagenti
    % definisco la struttura agente (scorte, metabolismo e vista)
	agente(i).scorte=10;
  	agente(i).met=ceil(5*rand);
   	agente(i).vista=ceil(6*rand);
   	% definisco le coordinate dell'agente
	agente(i).x=ceil(prm.nc*rand);
	agente(i).y=ceil(prm.nr*rand);
	% controllo che le coordinate scelte casualmente siano libere
        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, aggiorno matagenti
	% e non calcolo altre coordinate 
	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 tutti gli agenti
for j=1:prm.numagenti
%while contatore_delle_mosse~=prm.numagenti

	% Seleziono in maniera casuale l'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 non sia vuota la matrice delle destinazioni
	if length(poss_destinazioni) ~= 0
        % Ora abbiamo a disposizione il vettore delle destinazioni possibili
        % Dobbiamo scegliere 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
        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 all'agente la quantita' di zucchero corrispondente al suo metabolismo
function [agente] = consuma(agente,prm)

for i=1:prm.numagenti
	% Toglie dalla scorta dell'agente la quantita' pari al suo metabolismo
	agente(i).scorte = agente(i).scorte - agente(i).met;
end

%==========================================================================
%ANGELO DELLA MORTE
%==========================================================================
% Funzione che uccide gli agenti che non hanno sufficienti scorte per sopravvivere
function [agente,matagenti,prm] = ammazza(agente,matagenti,prm)

for i=prm.numagenti:-1:1
    if agente(i).scorte <= 0
        % Elimino l'agente dalla matrice "matagenti"
		matagenti(agente(i).y,agente(i).x) = 0;
		% Elimino l'agente di indice i
		agente(i) = [];
		% Abbasso di uno il numero di agenti nella struttura prm
		prm.numagenti = prm.numagenti - 1;
        
	end
end

%==========================================================================
%CRESCITA SUGARSCAPE
%==========================================================================
% Funzione che fa crescere lo 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 crescita 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;

% 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;
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);


%==========================================================================
%VISUlIZZAZIONE RISULTATI
%==========================================================================

% Funzione che visualizza i vari grafici durante l'evoluzione del sistema
%**************************************************************************
function visualizza_sugarscape(sugarscape,agente,prm,step)

% Visualizzazione degli agenti e del loro posizionamento sullo sugarscape
% Disegno dello sugarscape
figure(1);
imagesc(sugarscape);
hold on;
eval(['title (''sugarscape step ' ,int2str(step),''')']);
% Disegno gli agenti
if(prm.numagenti~=0)
    for i=1:prm.numagenti
        plot(agente(i).x,agente(i).y,'yd');
    end
end   
pause(0.1);

% Funzione che visualizza le varie statistiche
%***************************************************************************
function 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 = [];

% Converto la struttura in vettori
for i=1:length(stats)
	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_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;	
    vettore_scorta_min(length(vettore_scorta_min)+1) = stats(i).scorta_min;

end

% Visualizzo il grafico con l'andamento della popolazione degli agenti
figure(2);
hold on;
title('Andamento della popolazione degli agenti');
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 medio del metabolismo e della 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 minimo, medio e massimo della ricchezza e con la mediana
figure(4);
hold on;
title('Andamento minimo, medio e massimo della ricchezza 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    

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

%========================
%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
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