Autômato Celular para simulação da Esquistossomose em uma região.
Izabela Ralime / Manoel Rodrigues /Ronaldo Ramos
       25 de junho de 2008

(* 1- ESCOLHA DO DIRETÓRIO*)
(* ------------------------------------------------------------------------------------- *)

If[ChoiceDialog["ESTE É O CÓDIGO DE UM AUTÔMATO CELULAR QUE SIMULA A EXPANSÃO DA ESQUISTOSSOMOSE NUMA DETERMINADA ÁREA.

AUTORES:
        IZABELA RALIME
        MANOEL RODRIGUES
        RONALDO RAMOS

ATUALIZADO POR :
		GENIVALDO GUEIROS
		CAIO XAVIER
"],
dados=InputString["ENTRE COM O NOME DO DIRETÓRIO ONDE ESTÁ SALVO O BANCO DE DADOS DO LEVANTAMENTO DE CAMPO:

A pasta deve estar salva no diretório C:"]];
(* ------------------------------------------------------------------------------------- *)

(* 2- MATRIZES AUXILIARES*)

(* Fução que cria matrizes nulas de dimensão (m x n) e (m x n x 3) *)
(* ------------------------------------------------------------------------------------- *)
matrizAuxMolusco[A_]:=(Table[0,{i,Dimensions[A][[1]]},{j,Dimensions[A][[2]]},{k,3}])
matrizAux[A_]:=(Table[0,{i,Dimensions[A][[1]]},{j,Dimensions[A][[2]]}]);
(* ------------------------------------------------------------------------------------- *)

(* 3- IMPORTAÇÃO DOS DADOS*)

(* Iportação dos dados a serem utilizados na simulação*)
(* ------------------------------------------------------------------------------------- *)
populacaoSuscetivel=Import["C:/"<>dados<>"/populacaoSucetivel.txt","Table"];
populacaoInfectada=Import["C:/"<>dados<>"/populacaoInfectada.txt","Table"];
populacaoRecuperada=Import["C:/"<>dados<>"/populacaoRecuperada.txt","Table"];

criarmatrizMoluscoSadio[A_]:=(a=matrizAuxMolusco[A];
	For[i=1,i≤Dimensions[A][[1]],i++,
		a[[i]]=Import["C:/"<>dados<>"/populacaoMoluscoSadio/"<>IntegerString[i,10]<>".txt","Table"]];Return[a])

populacaoMoluscoSadio=criarmatrizMoluscoSadio[populacaoSuscetivel];

criarmatrizMoluscoInfectada[A_]:=(a=matrizAuxMolusco[A];  
	 For[i=1,i≤Dimensions[A][[1]],i++,
		a[[i]]=Import["C:/"<>dados<>"/populacaoMoluscoInfectada/"<>IntegerString[i,10]<>".txt","Table"]];Return[a])

populacaoMoluscoInfectada=criarmatrizMoluscoInfectada[populacaoSuscetivel];

dadosChuva=Import["C:/"<>dados<>"/dadosChuva.txt","Table"];
areaAlagadaMedia=Import["C:/"<>dados<>"/areaAlagadaMedia.txt","Table"];

criarmatrizConexaoAgua[A_]:=(a=matrizAuxMolusco[A];
		For[i=1,i≤Dimensions[A][[1]],i++,
			a[[i]]=Import["C:/"<>dados<>"/conexao/"<>IntegerString[i,10]<>".txt","Table"]];Return[a])

conexaoAgua=criarmatrizConexaoAgua[populacaoSuscetivel];
(* ------------------------------------------------------------------------------------- *)

(* 4- DRISTIBUIÇÃO DAS CHUVAS *)

(* Dados retirado da estação Goiana *)
(* ------------------------------------------------------------------------------------- *)

interpolacaoDados = Interpolation[dadosChuva];

chuvaAreaAlagada[A_,d_]:=(a=matrizAux[A];
		For[i=1,i≤Dimensions[A][[1]],i++,
			For[j=1,j≤Dimensions[A][[2]],j++,
				a[[i,j]]=A[[i,j]]*(interpolacaoDados [d]/200)]];Return[a])

Plot[interpolacaoDados[x], {x,1,365}, Epilog->Map[Point, dadosChuva],AxesOrigin->{1,0},AxesLabel->{dias,mm³ },ColorFunction->Function[{x,y},Hue[y]]]
(* ------------------------------------------------------------------------------------- *)

(* Grafico ilustrativo da distribuição das chuvas durante o ano *)

(* 5- CRESCIMENTO DOS MOLUSCOS *)

(* Função de Crescimento Populacional:Modelo de Verhulst *)
(* ------------------------------------------------------------------------------------- *)
n[t_,n0_,l_,k_]:=If[(n0+(l-n0) (E^(-k t)))==0,Return[0],(l n0)/(n0+(l-n0) (E^(-k t)))]
(* ------------------------------------------------------------------------------------- *)


(* função que retorna o número limite de molusco por celula *)
(* ------------------------------------------------------------------------------------- *)
limiteMolusco[i_,j_]:=(areaAlagada[[i,j]]*500);  (* Limite da população em função do ambiente *)

constanteMolusco=0.02; (* constante de crescimento *)
(* area alagavel de cada celula *)
(* ------------------------------------------------------------------------------------- *)


(* função de crescimento  do molusco *)
(* ------------------------------------------------------------------------------------- *)
(* entrar com a populacao de molusco sadio e infectado *)
crescimentoMolusco[A_]:=(a=matrizAuxMolusco[A];For[i=1,i≤Dimensions[A][[1]],i++,For[j=1,j≤Dimensions[A][[2]],j++,For[k=1,k≤Dimensions[A][[3]],k++,a[[i,j,k]]=n[1,A[[i,j,k]],limiteMolusco[i,j],constanteMolusco]]]];Return[N[a-A]])
(* ------------------------------------------------------------------------------------- *)


(* função que retorna a taxa de infecção dos moluscos *)
(* ------------------------------------------------------------------------------------- *)
taxaInfMolusco[i_,j_]:=(If[(300*taxaSaneamento*areaAlagada[[i,j]])==0,0,populacaoInfectada[[i,j]]/(300*taxaSaneamento*areaAlagada[[i,j]])])
(* ------------------------------------------------------------------------------------- *)


(* função de crescimento  do molusco infectado *)
(* ------------------------------------------------------------------------------------- *)
(* entra com populacaoMoluscoSadio e populacaoMoluscoInfectado *)
crescimentoMoluscoInfectado[A_,B_]:=(a=matrizAuxMolusco[A];b=A+B;
	For[i=1,i≤Dimensions[A][[1]],i++,
		For[j=1,j≤Dimensions[A][[2]],j++,
			For[k=1,k≤3,k++,
				If[B[[i,j,k]]≤b[[i,j,k]]/100,
					For[w=1,w≤A[[i,j,k]],w++,
						If[RandomReal[]≤taxaInfMolusco[i,j],a[[i,j,k]]=a[[i,j,k]]+1]]]]]];Return[a])
(* ------------------------------------------------------------------------------------- *)

(* 6- MOVIMENTAÇÃO DOS MOLUSCOS *)

(* função que simula o deslocamento do molusco pelas celulas *)
(* ------------------------------------------------------------------------------------- *)
taxaMovMoluscoSubir[d_]:=(1/(100*interpolacaoDados[d]))

taxaMovMoluscodescer[d_]:=(interpolacaoDados[d]/100000)

(* matriz que indica a conexão das celulas com suas vizinhas: 1->sobe o fluxo, 2->desce o fluxo, 3-> sobe o fluxo p/fora do mundo, 4->desce o fluxo p/fora do mundo *)
escolha1[i_,j_]:=(d=RandomChoice[Position[conexaoAgua[[i,j]],1]][[1]];Which[d==1,{x=i-1,y=j-1},d==2,{x=i-1,y=j},d==3,{x=i-1,y=j+1},d==4,{x=i,y=j-1},d==5,{x=i,y=j+1},d==6,{x=i+1,y=j-1},d==7,{x=i+1,y=j},d==8,{x=i+1,y=j+1}])

(* funções que determina qual deve ser a celula escolhida para descer o fluxo d'água  *)
escolha2[i_,j_]:=(d=RandomChoice[Position[conexaoAgua[[i,j]],2]][[1]];Which[d==1,{x=i-1,y=j-1},d==2,{x=i-1,y=j},d==3,{x=i-1,y=j+1},d==4,{x=i,y=j-1},d==5,{x=i,y=j+1},d==6,{x=i+1,y=j-1},d==7,{x=i+1,y=j},d==8,{x=i+1,y=j+1}])

passeioMolusco[A_,dias_]:=(a=matrizAuxMolusco[A];b=matrizAuxMolusco[A];
For[i=1,i≤Dimensions[A][[1]],i++,
	For[j=1,j≤Dimensions[A][[2]],j++,
		For[w=1,w≤Dimensions[A][[3]],w++,
			For[k=1,k≤A[[i,j,w]],k++,
				If[MemberQ[conexaoAgua[[i,j]],1]==True,
					If[RandomReal[]≤taxaMovMoluscoSubir[dias],aa=escolha1[i,j];a[[aa[[1]],aa[[2]],w]]=a[[aa[[1]],aa[[2]],w]]+1;
					b[[i,j,w]]=b[[i,j,w]]+1,
						If[MemberQ[conexaoAgua[[i,j]],2]==True,
							If[RandomReal[]≤taxaMovMoluscodescer[dias],aa=escolha2[i,j];a[[aa[[1]],aa[[2]],w]]=a[[aa[[1]],aa[[2]],w]]+1;
							b[[i,j,w]]=b[[i,j,w]]+1],
							If[MemberQ[conexaoAgua[[i,j]],4]==True,If[RandomReal[]≤taxaMovMoluscodescer[dias],b[[i,j,w]]=b[[i,j,w]]+1]]]],
								If[MemberQ[conexaoAgua[[i,j]],3]==True,
									If[RandomReal[]≤taxaMovMoluscoSubir[dias],b[[i,j,w]]=b[[i,j,w]]+1,
										If[MemberQ[conexaoAgua[[i,j]],2]==True,
											If[RandomReal[]≤taxaMovMoluscodescer[dias],aa=escolha2[i,j];
										      a[[aa[[1]],aa[[2]],w]]=a[[aa[[1]],aa[[2]],w]]+1;b[[i,j,w]]=b[[i,j,w]]+1]]],
													If[MemberQ[conexaoAgua[[i,j]],2]==True,
														If[RandomReal[]≤taxaMovMoluscodescer[dias],aa=escolha2[i,j];
														a[[aa[[1]],aa[[2]],w]]=a[[aa[[1]],aa[[2]],w]]+1;b[[i,j,w]]=b[[i,j,w]]+1],
															If[MemberQ[conexaoAgua[[i,j]],4]==True,
																If[RandomReal[]≤taxaMovMoluscodescer[dias],b[[i,j,w]]=b[[i,j,w]]+1]]]]]]]]];Return[{a,b}])
(* ------------------------------------------------------------------------------------- *)

(* 7- TAXAS E CONSTANTES RELATIVAS A POPULAÇÃO HUMANA *)

(* taxas relativas a população humana *)
(* ------------------------------------------------------------------------------------- *)
taxaSaneamento=0.4;

taxaCura=0.002;

taxaMovimentacao=0.8;

costanteInf=0.0003;

costanteReInf=0.000025;

0.00003
0.000025

(* função que retorna a taxa de infecção da celula *)
taxaInf[i_,j_]:=
	If[areaAlagada[[i,j]]*taxaSaneamento* costanteInf==0,0,(((populacaoMoluscoInfectada[[i,j,1]]*400+populacaoMoluscoInfectada[[i,j,2]]*4500+populacaoMoluscoInfectada[[i,j,3]]*2500))/areaAlagada[[i,j]]*taxaSaneamento* costanteInf)];


(* função que retorna a taxa de reinfecção da celula *)
taxaReInf[i_,j_]:=
	If[areaAlagada[[i,j]]*taxaSaneamento* costanteInf==0,0,(((populacaoMoluscoInfectada[[i,j,1]]*400+populacaoMoluscoInfectada[[i,j,2]]*4500+populacaoMoluscoInfectada[[i,j,3]]*2500))/areaAlagada[[i,j]]*taxaSaneamento* costanteReInf)];


(* glabrsta=4500, straminea=400, tenagophila=2500 *)
(* ------------------------------------------------------------------------------------- *)

(* 8- INFECÇÃO HUMANA *)



(* função que cria uma matriz com a soma dos humanos infectados apos o passeio *)
(* ------------------------------------------------------------------------------------- *)

(* função que cria uma matriz com a soma dos humanos infectados *)
(* ------------------------------------------------------------------------------------- *)
somaInfectado[A_]:=(a=matrizAux[A];
	For[i=1,i≤Dimensions[A][[1]],i++,
		For[j=1,j≤Dimensions[A][[2]],j++,
			For[k=1,k≤A[[i,j]],k++,
				If[RandomReal[]≤taxaInf[i,j],a[[i,j]]=a[[i,j]]+1]]]];
Return[a])
(* ------------------------------------------------------------------------------------- *)


(* função que cria uma matriz com a soma dos humanos curados *)
(* ------------------------------------------------------------------------------------- *)
somaCurados[A_]:=(a=matrizAux[A];
	For[i=1,i≤Dimensions[A][[1]],i++,
			For[j=1,j≤Dimensions[A][[2]],j++,
				For[k=1,k≤A[[i,j]],k++,
					If[RandomReal[]≤taxaCura,a[[i,j]]=a[[i,j]]+1]]]];
Return[a])
(* ------------------------------------------------------------------------------------- *)


(* função que cria uma matriz com a soma dos humanos Reinfectados *)
(* ------------------------------------------------------------------------------------- *)
somaReInfectado[A_]:=(a=matrizAux[A];
		For[i=1,i≤Dimensions[A][[1]],i++,
			For[j=1,j≤Dimensions[A][[2]],j++,
				For[k=1,k≤A[[i,j]],k++,
					If[RandomReal[]≤taxaReInf,a[[i,j]]=a[[i,j]]+1]]]];Return[a])
(* ------------------------------------------------------------------------------------- *)


(* função que provoca a movimentação pela vizinhança (8 vizinhos) *)
(* ------------------------------------------------------------------------------------- *)
passeio[i_,j_,A_]:=(If[i==1,{xmin=1,
xmax=i+1},If[i==Dimensions[A][[1]],{xmin=i-1,xmax=Dimensions[A][[1]]},{xmin=i-1,xmax=i+1}]];If[j==1,{ymin=1,
ymax=j+1},If[j==Dimensions[A][[2]],{ymin=j-1,ymax=Dimensions[A][[2]]},{ymin=j-1,ymax=j+1}]];xx={i,j};yy={i,j};While[ xx==yy,yy={RandomInteger[{xmin,xmax}],RandomInteger[{ymin,ymax}]}];Return[yy])


somaInfectadoPasseio[A_]:=(a=matrizAux[A];
	For[i=1,i≤Dimensions[A][[1]],i++,
		For[j=1,j≤Dimensions[A][[2]],j++,
			For[k=1,k≤A[[i,j]],k++,
				If[RandomReal[]≤taxaMovimentacao,b=passeio[i,j,A];
					If[RandomReal[]≤taxaInf[b[[1]],b[[2]]],a[[i,j]]=a[[i,j]]+1]]]]];
Return[a])
(* ------------------------------------------------------------------------------------- *)

(* função que cria uma matriz com a soma dos humanos Reinfectados apos o passeio *)
(* ------------------------------------------------------------------------------------- *)
somaReInfectadoPasseio[A_]:=(a=matrizAux[A];
	For[i=1,i≤Dimensions[A][[1]],i++,
		For[j=1,j≤Dimensions[A][[2]],j++,
			For[k=1,k≤A[[i,j]],k++,
				If[RandomReal[]≤taxaMovimentacao,b=passeio[i,j,A];
					If[RandomReal[]≤taxaReInf[b[[1]],b[[2]]],a[[i,j]]=a[[i,j]]+1]]]]];
Return[a])
(* ------------------------------------------------------------------------------------- *)

(* 9- NOVAS GERAÇÕES *)


(* função que calcula a razão infectatos por população total *)
(* ------------------------------------------------------------------------------------- *)
porcentagemPopulacao[A_,B_,C_]:=(a=matrizAux[A];For[i=1,i≤Dimensions[A][[1]],i++,For[j=1,j≤Dimensions[A][[2]],j++,If[A[[i,j]]+B[[i,j]]+C[[i,j]]==0,a[[i,j]]=0,a[[i,j]]=A[[i,j]]/(A[[i,j]]+B[[i,j]]+C[[i,j]])]]];Return[IntegerPart[10*N[a]]])
(* ------------------------------------------------------------------------------------- *)

(* função que calcula a proxima geração *)
(* ------------------------------------------------------------------------------------- *)
anos=0;
proximaGeracao[a_] := (If[a==1,Return[(porcentagemPopulacao[populacaoInfectada,populacaoSuscetivel,populacaoRecuperada])],areaAlagada=chuvaAreaAlagada[areaAlagadaMedia,a];populacaoMoluscoSadio=populacaoMoluscoSadio+crescimentoMolusco[populacaoMoluscoSadio]+crescimentoMolusco[populacaoMoluscoInfectada];moluscoInfectadoFuncao=crescimentoMoluscoInfectado[populacaoMoluscoSadio,populacaoMoluscoInfectada];populacaoMoluscoInfectada=populacaoMoluscoInfectada+moluscoInfectadoFuncao;populacaoMoluscoSadio=populacaoMoluscoSadio-moluscoInfectadoFuncao;moluscoSadioPasseio=passeioMolusco[populacaoMoluscoSadio,a];moluscoinfectadoPasseio=passeioMolusco[populacaoMoluscoInfectada,a];populacaoMoluscoSadio=populacaoMoluscoSadio+moluscoSadioPasseio[[1]]-moluscoSadioPasseio[[2]];populacaoMoluscoInfectada=populacaoMoluscoInfectada+moluscoinfectadoPasseio[[1]]-moluscoinfectadoPasseio[[2]];populacaoInfPasseio=somaInfectadoPasseio[populacaoSuscetivel];populacaoReInfPasseio=somaReInfectadoPasseio[populacaoRecuperada];populacaoSuscetivel=populacaoSuscetivel-populacaoInfPasseio;populacaoRecuperada=populacaoRecuperada-populacaoReInfPasseio;populacaoInfectada=populacaoInfectada+populacaoInfPasseio+populacaoReInfPasseio;infectados=somaInfectado[populacaoSuscetivel];reinfectados=somaReInfectado[populacaoRecuperada];
curados=somaCurados[populacaoInfectada];populacaoInfectada=populacaoInfectada+infectados+reinfectados-curados;populacaoSuscetivel=populacaoSuscetivel-infectados;populacaoRecuperada=populacaoRecuperada+curados-reinfectados;If[a==365,anos++];Return[porcentagemPopulacao[populacaoInfectada,populacaoSuscetivel,populacaoRecuperada]]])

(* 10- RELATÓRIO *)
(* função que gera um relatório por celula *)
relatorio[i_,j_]:=(Print[StringForm["população sucetível =          ``",populacaoSuscetivel[[i,j]]]]Print[StringForm["população recuperada =         ``",populacaoRecuperada[[i,j]]]]Print[StringForm["população infectada =          ``",populacaoInfectada[[i,j]]]]Print[StringForm["porcentagem infectada =        ``%",NumberForm[100*N[(populacaoInfectada[[i,j]]/(populacaoSuscetivel[[i,j]]+populacaoRecuperada[[i,j]]+populacaoInfectada[[i,j]]))],{2,2}]]]Print[StringForm["população molusco sadio =      ``",IntegerPart[populacaoMoluscoSadio][[i,j]]]]Print[StringForm["população molusco infectada =  ``",populacaoMoluscoInfectada[[i,j]]]]Print[StringForm["área alagada =                 ``m²",NumberForm[areaAlagada[[i,j]],{2,2}]]]Print[StringForm["taxa de infecção =             ``%",NumberForm[100*taxaInf[i,j],{2,2}]]]Print[StringForm["taxa de Reinfecção =           ``%",NumberForm[100*taxaReInf[i,j],{2,2}]]]Print[StringForm["taxa infc. molusco =           ``%",NumberForm[100*taxaInfMolusco[i,j],{2,2}]]];Return["Fim do relatório"])

(* 11- INTERFACE GRÁFICA *)

MostrarAreaInfectada[a_]:=Show[ArrayPlot[a,Mesh->True,ColorRules->{0-> LightCyan,1-> LightCyan,2->LightYellow,3->LightYellow,4->Yellow,5->Yellow,6->LightRed,7->LightRed,8->Red,9->Red,10->Black}]]

relatorio[3,2]

populacaoSuscetivel

populacaoInfectada

populacaoRecuperada

IntegerPart[populacaoMoluscoSadio]

populacaoMoluscoInfectada

anos