14.000x aceleração ou vitória da ciência da computação

Como desenvolvedor de software científico, faço muita programação. E a maioria das pessoas em outros campos científicos tende a pensar que programar é "apenas" lançar código e executá-lo. Tenho boas relações de trabalho com muitos colegas, inclusive de outros países ... Física, climatologia, biologia, etc. Mas quando se trata de desenvolvimento de software, você tem a nítida impressão de que eles pensam: “Ei, o que poderia ser difícil ?! Apenas escrevemos algumas instruções sobre o que o computador deve fazer, pressionamos o botão "Executar" e pronto - já temos a resposta! "



O problema é que é incrivelmente fácil escrever instruções que não significam o que você pensa. Por exemplo, um programa pode desafiar completamente a interpretação do computador. Além disso,literalmente, não há como saber se um programa será encerrado sem executá-lo. E há muitas, muitas maneiras de retardar muito a execução de um programa. Quer dizer ... muito devagar. Diminua a velocidade de forma que leve toda a sua vida ou mais. Isso geralmente acontece com programas que são escritos por pessoas sem educação em computação, ou seja, cientistas de outras áreas. Meu trabalho é consertar esses programas.



As pessoas não entendem que a ciência da computação ensina a teoria da computação, a complexidade dos algoritmos, a computabilidade (ou seja, podemos realmente calcular alguma coisa? Freqüentemente, tomamos como certo que podemos!) código que será executado em um período mínimo de tempo ou com o uso mínimo de recursos.



Deixe-me mostrar um exemplo de uma grande otimização em um script simples escrito por um colega meu.



Escalamos muito em climatologia. Fazemos leituras de temperatura e precipitação de um modelo climático global em grande escala e as comparamos a uma grade local de escala fina. Digamos que a grade global seja 50x25 e a grade local 1000x500. Para cada célula da grade na grade local, queremos saber a qual célula da grade global ela corresponde.



Uma maneira simples de representar o problema é minimizar a distância entre L [n] e G [n]. Acontece que essa pesquisa:



    L[i]:
      G[j]:
        L[i]  G[j]
       L[i] * G
    


Parece bastante simples. No entanto, se você olhar de perto, faremos muito trabalho extra. Observe o algoritmo em termos do tamanho da entrada.



    L[i]:                           #  L 
      G[j]:                        #  L x G 
        L[i]  G[j]                 #  L x G 
       d[i*j]              #  G  L  (L x G)
   ,       #  G  L  (L x G)


O código é parecido com este:



obs.lon <- ncvar_get(nc.obs, 'lon')
obs.lat <- ncvar_get(nc.obs, 'lat')
n.lon <- length(obs.lon)
n.lat <- length(obs.lat)

obs.lats <- matrix(obs.lat, nrow=n.lon, ncol=n.lat, byrow=TRUE)
obs.lons <- matrix(obs.lon, nrow=n.lon, ncol=n.lat)
obs.time <- netcdf.calendar(nc.obs)

gcm.lon <- ncvar_get(nc.gcm, 'lon')-360
gcm.lat <- ncvar_get(nc.gcm, 'lat')
gcm.lats <- matrix(gcm.lat, ncol=length(gcm.lat), nrow=length(gcm.lon),
                   byrow=TRUE)
gcm.lons <- matrix(gcm.lon, ncol=length(gcm.lat), nrow=length(gcm.lon))
gcm.lons.lats <- cbind(c(gcm.lons), c(gcm.lats))

# Figure out which GCM grid boxes are associated with each fine-scale grid point
# Confine search to 10 deg. x 10 deg. neighbourhood

dxy <- 10
mdist <- function(x, y)
    apply(abs(sweep(data.matrix(y), 2, data.matrix(x), '-')), 1, sum)
nn <- list()
for (i in seq_along(obs.lons)) {
    if((i %% 500)==0) cat(i, '')
    gcm.lims <- ((gcm.lons.lats[,1] >= (obs.lons[i]-dxy)) &
                 (gcm.lons.lats[,1] <= (obs.lons[i]+dxy))) &
                ((gcm.lons.lats[,2] >= (obs.lats[i]-dxy)) &
                 (gcm.lons.lats[,2] <= (obs.lats[i]+dxy)))
    gcm.lims <- which(gcm.lims)
    nn.min <- which.min(mdist(c(obs.lons[i], obs.lats[i]),
                        gcm.lons.lats[gcm.lims,]))
    nn[[i]] <- gcm.lims[nn.min]
}
nn <- unlist(nn)


Parece um algoritmo simples. É “fácil” calcular as distâncias e depois encontrar o mínimo. Mas, em tal implementação, conforme o número de células locais cresce, nosso custo computacional aumenta por seu produto com o número de células na grade global. Os dados ANUSPLIN canadenses contêm 1068 × 510 células (544 680 no total). Digamos que nosso GCM contenha 50 x 25 células (1250 no total). Assim, o custo do ciclo interno em "algumas unidades computacionais" é:



(c0eu×G)+(c1eu×G)+(c2eu×G)



onde membros cSão constantes que correspondem ao custo de calcular a distância entre dois pontos, encontrar o ponto mínimo e encontrar o índice da matriz. Na verdade, não nos importamos (muito) com membros constantes porque eles não dependem do tamanho da entrada. Portanto, você pode apenas adicioná-los e estimar o custo de computação;



(ceu×G)



Portanto, para este conjunto de entradas, nosso custo é 544,680×1,250=680,850,000...



680 milhões.



Parece muito ... embora, quem sabe? Os computadores são rápidos, certo? Se executarmos uma implementação ingênua, na realidade ela será executada um pouco mais rápido do que 1668 segundos, o que é um pouco menos de meia hora.



> source('BCCA/naive.implementation.R')
500 1000 1500 2000 2500 3000 ... 543000 543500 544000 544500 [1] "Elapsed Time"
    user   system  elapsed 
1668.868    8.926 1681.728 


Mas o programa deve ser executado por 30 minutos? Esse é o problema. Estamos comparando duas grades, cada uma com uma tonelada de estruturas que não usamos de forma alguma. Por exemplo, as latitudes e longitudes em ambas as grades estão em ordem de classificação. Portanto, para encontrar um número, você não precisa percorrer todo o array. Você pode usar o algoritmo de redução pela metade - olhe para o ponto no meio e decida qual metade do array queremos. Em seguida, pesquisar em todo o espaço terá o logaritmo de base dois de todo o espaço de pesquisa.



Outra estrutura importante que não usamos é o fato de que as latitudes vão em dimensãox, e longitudes - na dimensão y... Assim, em vez de realizar a operaçãox×y já que podemos fazer isso x+yTempo. Esta é uma grande otimização.



Como fica em pseudocódigo?



  local[x]:
    bisect_search(local[x], Global[x])

  local[y]:
    bisect_search(local[y], Global[y])

 2d      


Em código:



## Perform a binary search on the *sorted* vector v
## Return the array index of the element closest to x
find.nearest <- function(x, v) {
    if (length(v) == 1) {
        return(1)
    }
    if (length(v) == 2) {
        return(which.min(abs(v - x)))
    }
    mid <- ceiling(length(v) / 2)
    if (x == v[mid]) {
        return(mid)
    } else if (x < v[mid]) {
        return(find.nearest(x, v[1:mid]))
    }
    else {
        return((mid - 1) + find.nearest(x, v[mid:length(v)]))
    }
}

regrid.one.dim <- function(coarse.points, fine.points) {
    return(sapply(fine.points, find.nearest, coarse.points))
}

## Take a fine scale (e.g. ANUSPLINE) grid of latitudes and longitudes
## and find the indicies that correspond to a coarse scale (e.g. a GCM) grid
## Since the search is essentially a minimizing distance in 2 dimensions
## We can actually search independently in each dimensions separately (which
## is a huge optimization, making the run time x + y instead of x * y) and
## then reconstruct the indices to create a full grid
regrid.coarse.to.fine <- function(coarse.lats, coarse.lons, fine.lats, fine.lons) {
    xi <- regrid.one.dim(gcm.lon, obs.lon)
    yi <- regrid.one.dim(gcm.lat, obs.lat)
    ## Two dimensional grid of indices
    xi <- matrix(xi, ncol=length(fine.lats), nrow=length(fine.lons), byrow=F)
    yi <- matrix(yi, ncol=length(fine.lats), nrow=length(fine.lons), byrow=T)
    return(list(xi=xi, yi=yi))
}


O custo de cada pesquisa de bisseção é log do tamanho da entrada. Desta vez, nosso tamanho de entrada é dividido em espaços X e Y, então usaremosGx,Gy,eux e euy para Global, Local, X e Y.



cost=eux×euog2Gx+euy×euog2Gy+eux×euy



O custo é de 553.076. Bem, 553k soa muito melhor do que 680 milhões. Como isso afetará o tempo de execução?



> ptm <- proc.time(); rv <- regrid.coarse.to.fine(gcm.lat, gcm.lon, obs.lat, obs.lon); print('Elapsed Time'); print(proc.time() - ptm)[1] "Elapsed Time"
   user  system elapsed 
  0.117   0.000   0.117 
> str(rv)
List of 2
 $ xi: num [1:1068, 1:510] 15 15 15 15 15 15 15 15 15 15 ...
 $ yi: num [1:1068, 1:510] 13 13 13 13 13 13 13 13 13 13 ...
> 


0,117 segundos. O que costumava levar quase meia hora agora leva pouco mais de um décimo de segundo.



> 1668.868 / .117
[1] 14263.83


Soooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo, mas ainda estou surpreso e impressionado com o quanto o aumento de velocidade é. Aproximadamente 14 mil vezes .



Anteriormente, o script era tão lento que salvava o resultado no disco para revisão manual pelos cientistas antes de continuar. Agora tudo é calculado em um piscar de olhos. Fazemos esses cálculos centenas de vezes, de modo que, no final, economizamos dias e até semanas de tempo de computação. E a gente tem a oportunidade de interagir melhor com o sistema, para que a jornada de trabalho dos cientistas seja mais lucrativa ... eles não ficam parados esperando o fim dos cálculos. Tudo está pronto de uma vez.



Devo enfatizar que essas melhorias épicas de desempenho não requerem a compra de nenhum grande sistema de computador, paralelização ou aumento da complexidade ... na verdade, o código do algoritmo mais rápido é ainda mais simples e versátil! Essa vitória completa foi alcançada simplesmente pela leitura cuidadosa do código e algum conhecimento da complexidade computacional.



All Articles