[R-br] Criação de função fortran em pacote (base)

[Benilton Carvalho]

Pegue uma xicara de Nescau ou um bule de cha' de paciencia (pq o email
ficou longo) e so' passe para a linha seguinte apenas depois de ter
compreendido a linha atual.

Tente o seguinte:

f = function(i){
   message("Calculo com i=", i)
   2^i
}

E execute-a no R:

> f(1)
Calculo com i=1
[1] 2
> f(2)
Calculo com i=2
[1] 4
> f(3)
Calculo com i=3
[1] 8

Entao, retornando `as definicoes de logica de programacao la' do
comeco do curso, a funcao f() e' um conjunto de 2 instrucoes. A
primeira instrucao eh escrever na tela a mensagem "Calculo com
i=<valor>". A segunda instrucao e' determinar 2^i.

Tudo certo ate' aqui?

Agora, para economizar um pouco o seu tempo, vc nao precisa digitar
f(1), f(2),f(3). Voce pode usar lapply() para automatizar um pouco.

resultado0 <- lapply(1:3, f)

que vai retornar o seguinte:

> resultado0 <- lapply(1:3, f)
Calculo com i=1
Calculo com i=2
Calculo com i=3


Pausa para pensar no que aconteceu: o lapply(), de modo bem
simplificado, criou internamente um loop indo de 1 ate' 3 e chamou a
funcao f() para esses argumentos... A traducao do lapply() nesse caso
eh:

resultado1 = vector('list', 3)
for (i in 1:3)
   resultado1[[i]] = f(i)
resultado1

Para confirmar que eu nao estou enganado, teste:

identical(resultado0, resultado1)

Ou seja, em cada passo, as duas instrucoes contidas na sua funcao f()
foram executadas. Alem disso, o resultado da **ultima** instrucao eh o
resultado que eh guardado pelo sistema.

Agora, com o multicore::mclapply(), o que acontece eh o seguinte:

library(multicore)
resultado2 = mclapply(1:3, f, mc.cores=2)

"Traduzindo" o que o mclapply faz eh o seguinte:

1) Conte quantos elementos foram passados no primeiro argumento do
mclapply()... Foram passados 1, 2, 3 ; portanto, 3 elementos;
2) Veja se foi especificado o numero de nucleos a serem usados (nesse caso 2);
3) Abra 2 novas sessoes (invisiveis para o usuario) do R (por isso que
mc.cores=2) e copie a definicao da funcao f() dentro de cada nova
sessao (note: copia a funcao *INTEIRA* nao apenas a primeira ou
segunda linha);
4) Divida 3 tarefas (executar f(1), f(2) e f(3) ) entre as 2 sessoes
disponiveis;
4.1) No meu caso, a sessao 1 ficou encarregada de executar f(1) e
f(2); ao passo que a sessao 2, ficou por executar f(3)
5) Coloque as sessoes para executar suas tarefas. Nesse momento:
5.1) Sessao 1: calcula f(1)
5.2) Sessao 2: calcula f(2)
5.3) Sessao 1: nota que ainda tem que calcular f(2), entao vai fazer isso;
5.4) Sessao 2: grita "first!" e fica esperando a sessao 1 acabar;
5.5) Sessao 1: termina;
5.6) Sessao 1 e Sessao 2: avisam para o processo central do R (que eh
o q vc tem na tela) que estao prontos pra mandar os resultados
(lembre-se, o resultado eh apenas o resultado da 2a instrucao) de
volta;
5.7) Processo Central do R: colecta os resultados das Sessoes e guarda
em resultado2

Confirme que os resultados sao identicos:

identical(resultado1, resultado2)

**OBS IMPORTANTE**: Volte no passo 5 acima... Observe que cada Sessao
tem o conhecimento inteiro da definicao da funcao f()... As sessoes
nao abrem a funcao nem desmembram seu conteudo para execucao. Cada
sessao executa a funcao f() inteira, do inicio ao fim. Portanto, cada
sessao executa, independentemente das outras sessoes, a linha
message() e tambem a linha 2^i que defini originalmente na funcao f().

Pronto, dito isso a *sua* funcao f(), vc define como:

function(i){
    g1 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m1)
    g2 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m2)
    g3 <- mvrnorm(n=1, mat.cokri[[1]][2602:5202],m3)
    g4 <- mvrnorm(n=1, mat.cokri[[1]][2602:5202],m4)
}

cada linha g1/.../g4 eh apenas uma instrucao dentro da funcao
inteira... E cada nucleo executara' a funcao completa (note que
retornarah apenas o resultado g4, por um possivel erro seu na
definicao da funcao).

Assim, ao chamar:

mclapply(1:1, SUAFUNCAO, mc.cores=4)

vc esta' dizendo para o R que quer calcular o valor de:

SUAFUNCAO(1) ## pois 1:1 eh igual a 1

Apesar de vc ter reservado 4 nucleos, apenas 1 sera usado, pq vc so'
tem 1 funcao a ser calculada... (note que eu falei *funcao* nao
*instrucao*). Esse 1 nucleo, calculara TODAS as instrucoes g1/.../g4
usando internamente i=1.

Se vc chamar:

mclapply(1:2, SUAFUNCAO, mc.cores=4)

vc esta' dizendo ao R que vc quer determinar:

SUAFUNCAO(1)
SUAFUNCAO(2)

Apesar de ter reservado 4 nucleos, ele soh precisa de 2 nucleos para
calcular... Dai' SUAFUNCAO() sera' executada em cada nucleo, usando 2
sessoes... cada sessao usara' respectiva- e internamente i=1 e i=2.
Novamente, TODAS as instrucoes serao executadas por um dado nucleo.

Se vc chamar:

mclapply(1:4, SUAFUNCAO, mc.cores=4)

vc esta' dizendo ao R que vc quer:

SUAFUNCAO(1)
SUAFUNCAO(2)
SUAFUNCAO(3)
SUAFUNCAO(4)

E como vc alocou 4 nucleos, entao todos os 4 serao usados... A mesma
estoria acima se aplica (cada nucleo executara todas as instrucoes
dentro de SUAFUNCAO, usando internamente os valores 1, 2, 3 e 4 para a
variavel i).

Pronto... terminado...

O problema maior: a sua funcao nao usa a variavel "i" internamente...
Isso quer dizer que:

SUAFUNCAO(1) e SUAFUNCAO(2) tem o mesmo valor (que eh o valor de g4).

Desse jeito, o que quer q vc faca apenas repetira a mesma tarefa
outras tantas vezes.

Para ter g1/g2/g3/g4 executadas em paralelo, vc deve idealmente
re-escrever a funcao de acordo com o que eu descrevi na minha primeira
resposta desse ciclo de mensagens. Apenas para fins ilustrativos,
mostro abaixo uma possivel solucao que deve te dar o resultado
esperado (nao posso garantir, afinal nao tenho um exemplo
reproduzivel, nao e'?):

pars = list(list(n=1, mu=mat.cokri[[1]][1:2601], Sigma=m1),
                list(n=1, mu=mat.cokri[[1]][1:2601], Sigma=m2),
                list(n=1, mu=mat.cokri[[1]][2602:5202], Sigma=m3),
                list(n=1, mu=mat.cokri[[1]][2602:5202], Sigma=m4))
f = function(thePar)
   do.call(mvrnorm, thePar)
mclapply(pars, f,mc.cores=4)

(mas, volto a repetir: faca o exercicio de implementar o que sugeri
desde o comeco)

benilton