function [cap_carico,matagenti,agente,mondo,stats,coeff_gini] = sugarscape_fertile_ridotto()
%function [cap_carico,matagenti,agente,mondo,stats,coeff_gini] = sugarscape_fertile_ridotto();
%
% Questa versione dello Sugarscape é caratterizzata dal fatto che gli agenti si riproducono
% (genetica del sistema).
% Gli output sono nell'ordine:
% - la capacità portante del sistema, 
% - la matrice di 0 1 dove gli 1 rappresentano le celle dello Sugarscape occupate da un agente,
% - la struttura dati contenente le caratteristiche di ogni agente rimasto all'ultimo step,
% - la matrice che rappresenta lo Sugarscape,
% - la struttura contenente tutte le statistiche evolutive,
% - 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 contenente le varie statistiche
stats=[];

%=======
% Ciclo che fa evolvere il sistema
for step=1:prm.npassi
	step
	% Richiamo della funzione che muove gli agenti, li fa raccogliere e
	% riprodurrre
	[agente, matagenti, sugarscape, prm] = 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,prm] = statistiche(sugarscape,step,prm,agente,stats);
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);

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

% Funzione che crea lo sugarscape, con le due colline e tutti gli agenti
function [sugarscape, sugmax, agente, matagenti] = crea(prm)

% 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 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 della matrice sugmax (capacita' massima di carico del sistema)
for i=1:prm.nr
    for j=1:prm.nc
	% Creazione delle due colline che fanno parte del panorama
	collina1(i,j) = ceil(prm.altezza * exp(-(i-prm.x0) * (i-prm.x0)/sigmaxx) * exp(-(j-prm.y0) * (j-prm.y0)/sigmayy));
	collina2(i,j) = ceil(prm.altezza * exp(-(i-prm.x1) * (i-prm.x1)/sigmaxx) * exp(-(j-prm.y1) * (j-prm.y1)/sigmayy));
    end;
end;

% Combinazione delle due colline per generare il panorama completo
% rifarla con il max al posto del piu'
sugmax = collina1 + collina2 - 1;

% La matrice random viene livellata per rientrare nella 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

	% Definizione della struttura agente (eta', morte, metabolismo e vista)
	agente(i).eta = ceil(60*rand);
	agente(i).morte = 59 + ceil(41*rand);	
	agente(i).met = ceil(5*rand);
	agente(i).vista = ceil(6*rand);
	
	% Livello delle scorte pari a 10
	agente(i).scorte = 10;
	
	% Memoria delle scorte iniziali
	agente(i).scorte_iniziali = agente(i).scorte;
	
	% Definizione delle caratteristiche dell'agente necessarie per la genetica
	
	% Sesso (0 = femmina, 1 = maschio)
	if (rand < 0.5)
		agente(i).sesso = 0;
	else
		agente(i). sesso = 1;
	end
	
	% Inizio eta' fertile (intervallo [12,15] sia per maschi che per femmine)
	agente(i).fertile = 11 + ceil(4*rand);
	
	% Fine eta' fertile
	% (intervallo [40,50] per le femmine; intervallo [50,60] per i maschi)
	if (agente(i).sesso == 0)
		agente(i).menopausa = 29 + ceil(11*rand);
	else
		agente(i).menopausa = 39 + ceil(11*rand);
	end
	
	% 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

% Funzione che si occupa dello spostamento degli agenti
function [agente,matagenti,sugarscape,prm] = 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 for che fa guardare l'agente nelle quattro 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
		
	% 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;
		
		% Facciamo raccogliere all'agente lo zucchero presente 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;
		
		% Controlliamo se l'agente e' in eta' fertile e in possesso di una adeguata quantita' di
		% scorte per riprodursi.
		if ((agente(selezionato).scorte >= 0.4 * agente(selezionato).scorte_iniziali) && (agente(selezionato).eta >= agente(selezionato).fertile) && (agente(selezionato).eta <= agente(selezionato).menopausa))
		
			% Richiamo la funzione che calcola l'intorno dell'agente
			[intorno_1] = trova_intorno(matagenti,agente(selezionato),prm);
			
			% Individuo, nell'intorno, i partners di sesso opposto e fertili
			[agenti_sesso_opposto] = individua_altro_sesso(intorno_1,agente(selezionato).sesso,prm,agente);
			
			% Se il vettore agenti_sesso_opposto non e' vuoto, faccio riprodurre l'agente corrente
			if (length(agenti_sesso_opposto) > 0)
                		
				% Assegno valori random ad ogni partner (aggiungendo una colonna al vettore, che diventa cosi' una matrice)
				for i=1:length(agenti_sesso_opposto)
					agenti_sesso_opposto(i,2) = rand;
				end
				% Scelgo il partner con cui riprodurmi (il cui numero entra in b)
				[a,b] = max(agenti_sesso_opposto(:,2));
                		partner = agenti_sesso_opposto(b,1);
				
				% Calcolo l'intorno del secondo partner
				[intorno_2] = trova_intorno(matagenti,agente(partner),prm);
				
				% Individuo gli spazi liberi nell'intorno del genitore 1
				[spazi_liberi_genitore1] = individua_spazi_liberi(intorno_1);
				% Individuo gli spazi liberi nell'intorno del genitore 2
				[spazi_liberi_genitore2] = individua_spazi_liberi(intorno_2);
				
				% Controllo che le due matrici non siano vuote
				if ((length(spazi_liberi_genitore1) == 0) && (length(spazi_liberi_genitore2) == 0))
					spazi_liberi_genitori = [];
				elseif (length(spazi_liberi_genitore1) == 0)
					spazi_liberi_genitori = spazi_liberi_genitore2;
				elseif (length(spazi_liberi_genitore2) == 0)
					spazi_liberi_genitori = spazi_liberi_genitore1;
				else
					spazi_liberi_genitori = [spazi_liberi_genitore1; spazi_liberi_genitore2];
				end
				
				% Controllo l'esistenza di uno spazio libero dove far nascere il
				% figlio
				if (length(spazi_liberi_genitori) > 0)
				
					% Scelgo la cella dove far nascere l'agente
					% (in b vi e' il numero della riga scelta)
					[c,d] = max(spazi_liberi_genitori(:,3));
                                        					
					% Creazione del nuovo agente
					agente(prm.numagenti+1).x = spazi_liberi_genitori(d,1);
					agente(prm.numagenti+1).y = spazi_liberi_genitori(d,2);
					agente(prm.numagenti+1).eta = 0;
					agente(prm.numagenti+1).morte = 60 + ceil(40*rand);
					
					if (rand < 0.5)
						agente(prm.numagenti+1).sesso = 0;
					else
						agente(prm.numagenti+1).sesso = 1;
					end
					
					agente(prm.numagenti+1).fertile = 11 + ceil(4*rand);
					
					if (agente(i).sesso == 0)
						agente(prm.numagenti+1).menopausa = 29 + ceil(11*rand);
					else
						agente(prm.numagenti+1).menopausa = 39 + ceil(11*rand);
					end
	
					% Attribuzione di vista e metabolismo mediante crossover
					if (rand < 0.5)
						agente(prm.numagenti+1).met = agente(selezionato).met;
					else
						agente(prm.numagenti+1).met = agente(partner).met;
					end
					
					if (rand < 0.5)
						agente(prm.numagenti+1).vista = agente(selezionato).vista;
					else
						agente(prm.numagenti+1).vista = agente(partner).vista;
					end

					% Attribuzione delle scorte al nuovo agente (che inizia con 10)
					agente(prm.numagenti+1).scorte = 10;

					% Memoria delle scorte iniziali
					agente(prm.numagenti+1).scorte_iniziali = agente(prm.numagenti+1).scorte;
	
					% Aggiornamento di matagenti
					matagenti(agente(prm.numagenti+1).y,agente(prm.numagenti+1).x) = 1;

					% Aggiornamento di prm.numagenti
					prm.numagenti = prm.numagenti + 1;
					
					% Aggiornamento per le statistiche
					prm.nati = prm.nati + 1;
					
				end
			end
			
		end
	end
end

% Funzione che individua l'intorno dell'agente ricevuto in input
function [intorno] = trova_intorno(matagenti,agente_corrente,prm)

% Inizializza la matrice possibili_partners
intorno = [];

% Controlla il proprio intorno

% DESTRA
% Controlla che l'agente non sia sul bordo destro dello sugarscape
if (agente_corrente.x+1 <= prm.nc)
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x+1,agente_corrente.y,matagenti(agente_corrente.y,agente_corrente.x+1)];
end

% SINISTRA
% Controlla che l'agente non sia sul bordo sinistro dello sugarscape
if (agente_corrente.x-1 > 0)
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x-1,agente_corrente.y,matagenti(agente_corrente.y,agente_corrente.x-1)];
end

% ALTO
% Controlla che l'agente non sia sul bordo superiore dello sugarscape
if (agente_corrente.y+1 <= prm.nr)
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x,agente_corrente.y+1,matagenti(agente_corrente.y+1,agente_corrente.x)];
end

% BASSO
% Controlla che l'agente non sia sul bordo inferiore dello sugarscape
if (agente_corrente.y-1 > 0)
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x,agente_corrente.y-1,matagenti(agente_corrente.y-1,agente_corrente.x)];
end

% ALTO a DESTRA
% Controlla che l'agente non sia sul bordo destro e/o superiore dello sugarscape
if ((agente_corrente.y+1 <= prm.nr) && (agente_corrente.x+1 <  prm.nc))
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x+1,agente_corrente.y+1,matagenti(agente_corrente.y+1,agente_corrente.x+1)];
end

% ALTO e SINISTRA
% Controlla che l'agente non sia sul bordo sinistro e/o superiore dello sugarscape
if ((agente_corrente.y+1 <= prm.nr) && (agente_corrente.x-1 > 0))
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x-1,agente_corrente.y+1,matagenti(agente_corrente.y+1,agente_corrente.x-1)];
end


% BASSO A DESTRA
% Controlla che l'agente non sia sul bordo destro e/o inferiore dello sugarscape
if ((agente_corrente.y-1 > 0) && (agente_corrente.x+1 <= prm.nc))
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x+1,agente_corrente.y-1,matagenti(agente_corrente.y-1,agente_corrente.x+1)];
end


% BASSO A SINISTRA
% Controlla che l'agente non sia sul bordo sinistro e/o inferiore dello sugarscape
if ((agente_corrente.y-1 > 0) && (agente_corrente.x-1 > 0))
	% Aggiunge alla matrice intorno le varie coordinate e l'eventuale occupazione (0/1)
	intorno = [intorno; agente_corrente.x-1,agente_corrente.y-1,matagenti(agente_corrente.y-1,agente_corrente.x-1)];
end

% Funzione che restituisce gli agenti di sesso opposto, e fertili, presenti nell'intorno
function [agenti_sesso_opposto] = individua_altro_sesso(intorno,sesso_agente_corrente,prm,agente)

% Individua il numero di righe dell'array "intorno"
[m,n] = size(intorno);

% Inizializza il vettore "agenti_sesso_opposto"
agenti_sesso_opposto = [];

% Scorro tutti gli agenti
for j=1:prm.numagenti
	% Scorro tutto l'intorno
	for i=1:m
		% Controllo se l'agente e' uno di quelli presenti nell'intorno
		if (agente(j).x == intorno(i,1) && (agente(j).y == intorno(i,2)))
			% Nel caso in cui sia anche fertile, aggiungo il suo numero al vettore
			% agenti_sesso_opposto
			if ((agente(j).sesso ~= sesso_agente_corrente) && (agente(j).scorte >= 0.4 * agente(j).scorte_iniziali) && (agente(j).eta >= agente(j).fertile) && (agente(j).eta <= agente(j).menopausa))
				agenti_sesso_opposto = [agenti_sesso_opposto;j];
			end
		end
	end
end

% Funzione che restituisce le coordinate delle celle libere intorno all'agente ricevuto in input
function [spazi_liberi] = individua_spazi_liberi(intorno)

% Inizializzo la matrice "spazi_liberi"
spazi_liberi = [];

% Calcolo il numero di righe della matrice "intorno" (e inserisco questo valore nella variabile a)
[a,b] = size(intorno);

% Ciclo che scorre la matrice "intorno", eliminando le righe occupate (terza colonna uguale ad 1)
for i=a:-1:1
	% Controllo se la cella identificata dalla righe della matrice intorno e' occupata
	if (intorno(i,3) == 1)
		% In quel caso elimino la riga in questione
		intorno(i,:) = [];
	else
		% Nel caso in cui sia libera, aggiungo una quarta colonna con un valore random
		intorno(i,4) = rand;
	end
end

% Elimino la terza colonna (che ora e' sempre uguale a zero) dalla matrice "intorno"
intorno(:,3) = [];

% Eguaglio la variabile di output "spazi_liberi" ad "intorno"
spazi_liberi = intorno;

% 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

% Funzione che fa invecchiare gli agenti
function [agente] = invecchia(agente,prm)

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

% 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
	% Controllo se è fertile per aggionare la la statistica numero agenti fertili
	if ((agente(i).scorte >= 0.4 * agente(i).scorte_iniziali) && (agente(i).eta >= agente(i).fertile) && (agente(i).eta <= agente(i).menopausa))
        	prm.fertili=prm.fertili+1;
    	end
    
	if ((agente(i).scorte <= 0) || (agente(i).eta >= agente(i).morte))
		% 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;
		% Aggiorno le statistiche
		prm.morti = prm.morti + 1;
    	end
end

% Funzione che fa crescere lo sugarscape
function a = cresci(a)

% Aumenta di 0.5 il livello di zucchero di ogni cella dello Sugarscape
a=a+0.5;

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

% Calcolo delle statistiche
function [stats,prm] = statistiche(sugarscape,step,prm,agente,stats);

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

% Parametri relativi allo step appena concluso
stats(step).morti = prm.morti;
stats(step).nati = prm.nati;
stats(step).fertili = prm.fertili;
stats(step).sugarmedio= mean(mean(sugarscape));

% Riazzero questi parametri
prm.morti = 0;
prm.nati = 0;
prm.fertili = 0;

% 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);
    
	% 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).met_sd=0;
	stats(step).vista_sd=0;
end

% Funzione che visualizza il movimento degli agenti durante l'evoluzione del sistema
function visualizza_sugarscape(sugarscape,agente,prm,step)

% Disegno dello sugarscape
figure(1);
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_vista_sd = [];
vettore_met_sd = [];
vettore_scorta_min = [];
vettore_scorta_media = [];
vettore_scorta_max = [];
vettore_nati = [];
vettore_morti = [];
vettore_fertili = [];
vettore_sugarmedio = [];

% Converto la struttura in vettori
for i=1:prm.npassi
	vettore_numagenti(length(vettore_numagenti)+1) = stats(i).numagenti_vivi;
	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_nati(length(vettore_nati)+1) = stats(i).nati;
	vettore_morti(length(vettore_morti)+1) = stats(i).morti;
	vettore_fertili(length(vettore_fertili)+1) = stats(i).fertili;
	vettore_sugarmedio(length(vettore_sugarmedio)+1) = stats(i).sugarmedio;
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');
legend('Metabolismo medio','Vista media');
plot(vettore_metmedio,'r*');
plot(vettore_vista_media,'b*');

% Visualizzo il grafico con l'andamento delle deviazione standard del
% metabolismo medio
figure(10);
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(11);
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
figure(4);
hold on;
title('Andamento della ricchezza minima, media e massima degli agenti');
plot(vettore_scorta_min,'r');
plot(vettore_scorta_media,'b');
plot(vettore_scorta_max,'g');
legend('Ricchezza minima', 'Ricchezza media', 'Ricchezza massima');
plot(vettore_scorta_min,'r*');
plot(vettore_scorta_media,'bo');
plot(vettore_scorta_max,'gp');

% 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

% Visualizzo il grafico con l'andamento degli agenti nati e morti
figure(7);
hold on;
title('Andamento degli agenti nati e morti');
plot(vettore_nati,'r');
plot(vettore_morti,'b');
legend('Numero agenti nati', 'Numero agenti morti');

% Visualizzo il grafico con l'andamento degli agenti fertili
figure(8);
hold on;
title('Andamento del numero di agenti fertili');
plot(vettore_fertili,'g');

% Visualizzo il grafico con l'andamento della media dello sugarscape
figure(9);
hold on;
title('Andamento della media dello sugarscape');
plot(vettore_sugarmedio,'r');
plot(vettore_sugarmedio,'gp');

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

prm.morti = 0;
prm.nati = 0;
prm.fertili = 0;

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

%Funzione che scrive le statistiche 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\t%s\t%s\t%s\t%s\t%s\r\n','step','sug_medio','met_medio','met_sd','vista_media','vista_sd','scorta_media','scorta_min','scorta_max','numagenti_vivi','numagenti_morti','numagenti_nati','numagenti_fertili');

% 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\t%d\t%d\t%d\t%d\t%d\r\n',i,stats(i).sugarmedio,stats(i).met_medio,stats(i).met_sd,stats(i).vista_media,stats(i).vista_sd,stats(i).scorta_media,stats(i).scorta_min,stats(i).scorta_max,stats(i).numagenti_vivi,stats(i).morti,stats(i).nati,stats(i).fertili);
end