Programação para as Ciências Experimentais 2008/9

Ludwig Krippahl, 2009

Programação para as Ciências Experimentais

2008/9

Teórica 12

Ludwig Krippahl, 2009 2

Na aula de hoje...

Minimização multidimensional Exemplo: estimativa do efeito de um

antibiótico no crescimento bacteriano. Trabalho Prático 2


Crescimento bacteriano

Equação de Verhulst:

dB/dt = cB – mB2




dB/dt = cB – mB2

Variação do número de organismos




dB/dt = cB – mB2

É o ritmo de crescimento vezes o número de organismos




dB/dt = cB – mB2

Menos a taxa de mortalidade vezes o quadrado desse número. A mortalidade

resulta da competição por recursos.





Problema:• Dado um conjunto de medições, ajustar os

parâmetros da equação



Problema:• Dado um conjunto de medições, ajustar os

parâmetros da equação

• Mas são dois parâmetros: crescimento e mortalidade. Precisamos de uma minimização a duas dimensões.


Minimização multidimensional

Método:• Minimizar uma variável de cada vez até

chegar a um ponto fixo, a menos da precisão desejada



Método :• Partir de um ponto inicial, um valor para cada

variável.

• Encontrar o mínimo de uma variável.

• Alterar o vector das variáveis

• Encontrar o mínimo da próxima.

• Repetir para todas, as vezes que for necessário.


Minimização multidimensionalPonto inicial



Mínimo de X



Mínimo de X

Mínimo de Y



Mínimo de X

Mínimo de Y

Novo mínimo de X



Vamos modificar a minfn• Para partir do ponto dado e não ser preciso

especificar os três pontos iniciais (é mais eficiente começar com 3 pontos juntos quando próximo do mínimo)

• Para procurar o mínimo de uma de várias variáveis.


Os 3 pontos iniciais

X1 é o ponto dado



X1

Xm próximo de X1



X1

Xm próximo de X1Desce?

Se não, troca



X1

Xm

X2 a 1.618*(Xm-X1)



X1

Xm

X2

Y2>Ym?

Não, continua:

X1=Xm

Xm=X2



X1

Xm

Y2>Ym?

Não, continua:

X1=Xm

Xm=X2



X1

Xm

Y2>Ym?

Não, continua:

X1=Xm

Xm=X2

X2



X1 Xm

Y2>Ym?

Sim.Devolve:

X1 Xm X2 Ym

(Ym para começar a minimização)

X2



function [x1,xm,x2,ym]=

mininicial(funcao,params,vars,indice,delta)





Os valores a devolver, os 3 pontos de x e o y do meio que precisamos para começar a minimização.





Nome da função a minimizar e os parâmetros constantes que precisamos para a avaliar. E.g.: os coeficientes do polinómio, os dados experimentais, etc





Um vector com os valores das N variáveis da função (a função tem várias dimensões)





O índice no vector da variável que estamos a minimizar agora (a função tem várias, mas só conseguimos lidar com uma de cada vez)





Tamanho do passo inicial para calcular o primeiro Xm a partir do X inicial, que será o valor que vem em vars(indice).



razaod=1.618;

x1=vars(indice);y1=feval(funcao,params,vars);

xm=x1+delta;vars(indice)=xm;ym=feval(funcao,params,vars);

Razão dourada



razaod=1.618;



Calcular x1 e y1.

Nota: a função cujo nome foi dado em funcao precisa de todo o vars, mas só alteramos o elemento indice



razaod=1.618;



Mesma coisa para xm e ym



if ym>y1 t=ym;ym=y1;y1=t;t=xm;xm=x1;x1=t;

endifx2=xm+razaod*(xm-x1);vars(indice)=x2;y2=feval(funcao,params,vars);

Se for “a subir” troca o x1 com o xm, e o y1 com o ym.





Calcula o x2 a uma distância de xm igual a 1.618.. vezes o intervalo (xm-x1).





Se x1>xm, continua para valores maiores.

Se foi trocado, x1-xm é negativo e segue para valores menores.



while y2<ymx1=xm;xm=x2;ym=y2;x2=xm+razaod*(xm-x1);vars(indice)=x2;y2=feval(funcao,params,vars);

endwhile

Enquanto continua “a descer”, avança com os pontos,


Minimização

Antes era:

function xm=minfn(func,params,x1,xm,x2,prec)

Agora é:

function

xm=minfnvec(func,params,vars,indice,prec)



function


Valor da variável que está a minimizar no mínimo da função considerando apenas esta variável.



function


Função (o nome, em string) e parâmetros constantes, como de costume.



function


Vector com os valores das várias variáveis no ponto inicial, de onde parte à procura do mínimo.



function


Índice da variável onde procurar o mínimo.



function


Precisão (tamanho do intervalo abaixo do qual consideramos ter encontrado o mínimo)



Esta função é igual à minfn, excepto:• Usa o mininicial para determinar os 3 pontos e

o valor de ym

[x1,xm,x2,ym]=mininicial(func,params,vars,indice,prec);

Nota: um bom valor para o delta é a precisão: começamos do intervalo mais pequeno.




o valor de ym

• Tem de atribuir o valor correcto a vars(indice) antes de chamar a função fornecida em func

xn=c1*xm+c2*x1;

vars(indice)=xn;

yn=feval(func,params,vars);




o valor de ym

• Tem que atribuir o valor correcto a vars(indice) antes de chamar a função fornecida em func

Mas minfnvec ainda só minimiza numa dimensão (a dimensão indicada em indice).



Precisamos de:

function xs=multimin(funcao,params,xs,prec)




Vector com os valores de todas as variáveis no mínimo da função




Nome da função a minimizar.




Parâmetros constantes...




Ponto inicial (valores de todas as variáveis da função de onde partir à procura do mínimo)




Precisão



while truexvs=xs;for f=1:length(xs)

xs(f)=minfnvec(funcao,params,xs,f,prec);endfordisp("Valores até agora:")disp(xs);fflush(stdout)if sum(abs(xs-xvs))<prec

breakendif

endwhileCiclo só termina no break.





breakendif

endwhile

Guarda os valores antigos dos xs





breakendif

endwhile

Minimiza em cada dimensão, actualizando o valor nos xs




xs(f)=minfnvec(funcao,params,xs,f,prec);endfordisp("Valores até agora:")disp(xs);fflush(stdout);if sum(abs(xs-xvs))<prec

breakendif

endwhile

Mostra o progresso do cálculo indicando os valores correntes.





breakendif

endwhile

Se o total da variação absoluta das variáveis é inferior à precisão, acabou.



De volta ao problema:

dB/dt = cB – mB2

• Integramos pelo método de Euler, com a função:

function mat=crescimento(cresc,mort,dt,qini,tfinal)




Matriz com os valores de tempo e número de bactérias em duas colunas




Taxa de crescimento




Taxa de mortalidade




Passo de integração




Quantidade inicial de organismos.




Tempo final.



mat=[0,qini];

B=qini;

for t=dt:dt:tfinal

dB=B*cresc-B^2*mort;

B=B+dB*dt;

mat=[mat;t,B];

endfor



Agora precisamos de calcular o erro do modelo aos dados experimentais. Análogo ao que fizemos para as reacções químicas, mas desta vez com duas variáveis.



function err=errocres(dados,vars)

mat=crescimento(vars(1),vars(2),10,0.1,400);

y=interpol(mat,dados(:,1));

err=sum((y-dados(:,2)).^2);

endfunction

Matriz com a simulação, 400 minutos, passo de 10 minutos.







endfunction

Nota: Quantidade em “kilobactérias”. Explicação adiante...







endfunction

Interpolar os valores simulados para os pontos dos dados experimentais.







endfunction

Erro quadrático....







endfunction

Para o erro quadrático médio usar mean ou dividir pelo total.



Para testar, simulamos dados com estes parâmetros: 10 pontos de 30 em 30 minutos. Nota: cada linha da matriz são 10 minutos.

vals=[0.040234,0.001877];

mat=crescimento(vals(1),vals(2),10,0.1,400);

dados=mat(3:3:30,:)



Escolhemos um ponto inicial diferente, e minimizamos:

xs=multimin("errocres",dados,[0.05,0.005],1e-4)



Simulamos com os parâmetros calculados e comparamos:

mat2=crescimento(xs(1),xs(2),10,0.1,400);hold offplot(dados(:,1),dados(:,2),"or");hold onplot(mat2(:,1),mat2(:,2));





Nota sobre as “kilobactérias”:• Com esta equação, se contarmos em unidades de uma

bactéria o parâmetro da mortalidade tem de ser mil vezes mais pequeno.

• Em geral, é melhor escolher as unidades de forma a que a função tenha uma escala semelhante nas várias dimensões. Desta forma a taxa de crescimento e de mortalidade têm apenas uma ordem de grandeza de diferença em vez de quatro.


Processar dados experimentais

Queremos estudar o efeito da meticilina no crescimento de uma bactéria.

Duas pessoas, Ana e Carlos, cresceram lotes da bactéria em meios com e sem meticilina, e contaram as colónias de amostras retiradas de 30 em 30 minutos.

A concentração inicial era de 100 bactérias por ml.



Os dados estão em 20 ficheiros 1.txt a 20.txtDados de crescimento

Meio:Normal

Preparador:Ana

25;0

59;1

...

296;40



Objectivo: ajustar o modelo de crescimento às 2 condições e comparar os parâmetros

• Ler os ficheiros para uma lista de estruturas

• Separar as medições por meio e/ou preparador, em matriz

• Calcular parâmetros.



Ler os ficheiros para uma lista de estruturas

function dados=leficheiros(num)


Ler os ficheiros

dados=[];for f=1:num

fid=fopen([num2str(f),".txt"],"r");reg.valores=[];while !feof(fid)

...endwhiledados(f)=reg;fclose(fid);

endfor

Percorre o número indicado de ficheiros numero.txt


Ler os ficheiros




endfor

Matriz para os valores neste ficheiro


Ler os ficheiros




endfor

Ler o ficheiro


Ler os ficheiros




endfor

Acrescenta registo à lista e fecha o ficheiro


Ler os ficheiross=fgetl(fid);if !isstr(s)

breakendifif findstr(s,"Meio:")!=[]

reg.meio=s(6:length(s));elseif findstr(s,"Preparador:")!=[]

reg.prep=s(12:length(s));elseif findstr(s,";")!=[]

m=split(s,";");reg.valores=[reg.valores;str2num(m(1,:)),str2num(m(2,:))];

endif;endwhile

Lê uma linha e testa se o resultado é string. Se não for é por ser -1, o que indica que não há linha para ler. Nesse caso termina o ciclo (pode haver linhas vazias no final do texto).



breakendifif !isempty(findstr(s,"Meio:"))

reg.meio=deblank(s(6:length(s)));elseif !isemoty(findstr(s,"Preparador:"))

reg.prep=deblank(s(12:length(s)));elseif !isempty(findstr(s,";"))


endif;endwhile

“Meio:” indica que se segue o meio (Normal ou Meticilina)







endif;endwhile

O preparador (Ana ou Carlos)







endif;endwhile

Um ; indica que é uma linha com os valores. Split, depois acrescenta à matriz.


Ler os ficheiros Exemplo:octave:25> l=leficheiros(20);octave:26> ll =( [1] = { meio = Normal prep = Ana valores = 34 0 ... 304 42 } [2] = ...


Organizar os dados

Queremos receber uma matriz tempo, contagens para cada meio e/ou preparador:

function mat=compiladados(lista,prep,meio)

Recebe a lista e dois strings com o preparador e meio (“” para qualquer um)


Organizar os dados

function mat=compiladados(lista,prep,meio)mat=[];for f=1:length(lista)

reg=lista(f);if (isempty(meio) || strcmp(meio,reg.meio)) && (isempty(prep)|| strcmp(prep,reg.prep))

mat=[mat;reg.valores];endif

endfor

Percorre a lista elemento a elemento.


Organizar os dados

function mat=compiladados(lista,prep,meio)mat=[];for f=1:length(lista)

reg=lista(f);if (isempty(meio) || strcmp(meio,reg.meio))

&& (isempty(prep)|| strcmp(prep,reg.prep)) mat=[mat;reg.valores];

endifendfor

Se é este o meio ou preparador, ou se “”, acrescenta.


Organizar os dados

Exemplo: todos os dados

l=leficheiros(20);

dados=compiladados(l,"","");

clearplot;

plot(dados(:,1),dados(:,2),"o");


Organizar os dados

Exemplo: todos os dados


Organizar os dados

Exemplos: separados por meio

dmet=compiladados(l,"","Meticilina");dsem=compiladados(l,"","Normal");clearplothold onplot(dmet(:,1),dmet(:,2),"or;Meticilina;");plot(dsem(:,1),dsem(:,2),"xg;Normal;");


Organizar os dados


Organizar os dados

Exemplos: separados por preparador

ana=compiladados(l,"Ana","");carlos=compiladados(l,"Carlos","");clearplothold onplot(ana(:,1),ana(:,2),"or;Ana;");plot(carlos(:,1),carlos(:,2),"xg;Carlos;");


Organizar os dados

Exemplos: separados por preparador


Ajustar o modelo

Separamos por meio:dmet=compiladados(l,"","Meticilina");

dsem=compiladados(l,"","Normal");

E minimizamos, a partir de uma estimativa inicial:

xs=multimin("errocres",dsem,[0.05,0.002],1e-4);

xm=multimin("errocres",dmet,[0.05,0.002],1e-4);


Ajustar o modelo

Simulamos com os parâmetros calculados:

sims=crescimento(xs(1),xs(2),10,0.1,400);

simm=crescimento(xm(1),xm(2),10,0.1,400);


Ajustar o modelo

E comparamos os dados com a simulação:

clearplothold onplot(dmet(:,1),dmet(:,2),"or");plot(dsem(:,1),dsem(:,2),"xg");plot(simm(:,1),simm(:,2),"-r;Meticilina;");plot(sims(:,1),sims(:,2),"-g;Normal;");


Ajustar o modelo

Compara-se no gráfico:


Exportar o resultado

Escrita formatada: fprintffprintf(id, formato, dados) Exemplo: escrever tabela em duas

colunas separadas por tab.fid=fopen("relatorio.txt","w");mat=compiladados(l,"","");fprintf(fid,"%i\t%i\r\n",mat’)fclose(fid)



fprintf(fid,"%i\t%i\r\n",mat’)

%i indica que é um número inteiro.




\t é o caracter tab.




\r\n indica uma nova linha.



% é um caracter especial na string de formatação, indica que o que se segue especifica o formato (ver fprint no manual)

\ é um caracter especial em qualquer string, usado para caracteres que não são visíveis (mudar de linha, tab, etc.)

Para mostrar escrever dois: \\, %%



Quando o fprintf ou o printf (para escrever no ecrã) recebem uma matriz percorrem todos os elementos da matriz aplicando a formatação na ordem indicada;

Atenção: percorre a primeira coluna toda, depois a segunda, etc.. Como temos dados em colunas, usar transposta.


Trabalho 2

O trabalho 2 é parecido.• Ler os ficheiros com dados e modelos.

• Ajustar as constantes.

• Exportar resultados


Trabalho 2

Problema• Medicamento no estômago (ME)

• Medicamento no sangue (MS)

ME MS Excretado

Duas constantes


Trabalho 2



d[ME]/dt = -[ME] ka

d[MS]/dt = [ME] ka – [MS] ke

[ ] mg/Kg


Trabalho 2



• As constantes podem depender• Da idade (<30 anos, >=30 anos)

• De tomar em jejum ou depois da refeição


Trabalho 2

Objectivos• Determinar as constantes para os 4 casos

(idade e quando toma)

• Determinar a dose necessária (mg/Kg) para manter uma concentração média de 50mg/Kg durante 6 dias tomando de 12 em 12 horas, conforme o caso.


Trabalho 2

O que fazer:• Ler os ficheiros (estruturas?)


Trabalho 2

Idade:49

Administrado em jejum

Horas [MS]

2.36 69

4.06 74

6.35 56

...

Horas: tempo da amostra

[MS]: mg/Kg


Trabalho 2


• Calcular as constantes que ajustam o modelo aos dados de cada ficheiro• Minimizar para as duas constantes: 18 valores

• Assumir ka dado, minimizar só uma: 15 valores


Trabalho 2


• Calcular as constantes que ajustam o modelo aos dados de cada ficheiro

• Calcular as doses necessárias para uma média de 50mg/Kg, 12/12h, 6 dias• Opcional, 2 valores

• Não explico como...


Trabalho 2


• Calcular as constantes que ajustam o modelo aos dados de cada ficheiro

• Calcular as doses necessárias para uma média de 50mg/Kg, 12/12h, 6 dias

• Gravar tudo para 4 ficheiros com os relatórios• Constantes, dosagem, dados experimentais


Próximas

Prática• Ajustar as curvas de crescimento

• 2 aulas para o trabalho

Teóricas• 20: Excel, revisões

• 27: Introdução à bioinformática

• 3: Estrutura do exame, revisões, dúvidas.


Dúvidas

Documents

Programação para as Ciências Experimentais 2008/9