El problema es que es increíblemente fácil escribir instrucciones que no significan lo que piensas. Por ejemplo, un programa puede desafiar completamente la interpretación por computadora. Además,literalmente no hay forma de saber si un programa terminará sin ejecutarlo. Y hay muchas, muchas formas de ralentizar mucho la ejecución de un programa. Quiero decir ... realmente más lento. Reduzca la velocidad para que le lleve toda su vida o más completarlo. Esto sucede con mayor frecuencia con programas escritos por personas sin educación informática, es decir, científicos de otros campos. Mi trabajo es arreglar esos programas.
La gente no entiende que la informática te enseña la teoría de la computación, la complejidad de los algoritmos, la computabilidad (es decir, ¿podemos realmente calcular algo? ¡Con demasiada frecuencia damos por sentado que podemos!) código que se ejecutará en una cantidad mínima de tiempo o con un uso mínimo de recursos.
Permítanme mostrarles un ejemplo de una gran optimización en un simple script escrito por un colega mío.
Escalamos mucho en climatología. Tomamos lecturas de temperatura y precipitación de un modelo climático global a gran escala y las comparamos con una cuadrícula local de escala fina. Digamos que la cuadrícula global es 50x25 y la cuadrícula local es 1000x500. Para cada celda de la cuadrícula de la cuadrícula local, queremos saber a qué celda de la cuadrícula global corresponde.
Una forma fácil de representar el problema es minimizar la distancia entre L [n] y G [n]. Resulta tal búsqueda:
L[i]:
G[j]:
L[i] G[j]
L[i] * G
Parece bastante simple. Sin embargo, si miras de cerca, hacemos mucho trabajo adicional. Mire el algoritmo en términos del tamaño de la 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)El código se parece a esto:
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 un algoritmo simple. Es “fácil” calcular las distancias y luego encontrar el mínimo. Pero en tal implementación, a medida que crece el número de celdas locales, nuestro costo computacional aumenta en su producto con el número de celdas en la cuadrícula global. Los datos canadienses ANUSPLIN contienen 1068 × 510 celdas (544 680 en total). Digamos que nuestro GCM contiene 50x25 celdas (1250 en total). Por lo tanto, el costo del ciclo interno en "algunas unidades computacionales" es:
donde los miembros son constantes que corresponden al costo de calcular la distancia entre dos puntos, encontrar el punto mínimo y encontrar el índice de matriz. Realmente no nos importan (mucho) los miembros constantes, porque no dependen del tamaño de la entrada. De modo que puede sumarlos y estimar el costo de cálculo;
Entonces, para este conjunto de insumos, nuestro costo es .
680 millones.
Parece mucho ... aunque, ¿quién sabe? Las computadoras son rápidas, ¿verdad? Si ejecutamos una implementación ingenua, en realidad se ejecutará un poco más rápido de 1668 segundos, que es un poco menos de media 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 Pero, ¿ debería el programa funcionar durante 30 minutos? Ese es el problema. Estamos comparando dos cuadrículas, cada una de las cuales tiene una tonelada de estructuras que no hemos usado de ninguna manera. Por ejemplo, las latitudes y longitudes en ambas cuadrículas están ordenadas. Por lo tanto, para encontrar un número, no es necesario recorrer toda la matriz. Puede utilizar el algoritmo de reducción a la mitad: observe el punto en el medio y decida qué mitad de la matriz queremos. Luego, buscar en todo el espacio tomará el logaritmo base dos de todo el espacio de búsqueda.
Otra estructura importante que no usamos es el hecho de que las latitudes van en dimensión , y longitudes - en dimensión . Así, en lugar de realizar la operación veces podemos ejecutarlo veces. Esta es unagranoptimización.
¿Cómo se ve en pseudocódigo?
local[x]:
bisect_search(local[x], Global[x])
local[y]:
bisect_search(local[y], Global[y])
2d En codigo:
## 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))
}El costo de cada búsqueda de bisección es el registro del tamaño de entrada. Esta vez nuestro tamaño de entrada se divide en espacios X e Y, así que usaremos y para Global, Local, X e Y.
El costo es de 553 076. Bueno, 553k suena mucho mejor que 680m. ¿Cómo afectará esto al tiempo de ejecución?
> 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. Lo que solía tomar casi media hora ahora toma poco más de una décima de segundo.
> 1668.868 / .117
[1] 14263.83Taaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa, pero incluso yo estoy sorprendido e impresionado por lo mucho que es la aceleración. Aproximadamente 14 mil veces .
Anteriormente, el script era tan lento que guardaba el resultado en el disco para que los científicos lo revisaran manualmente antes de continuar. Ahora todo está calculado en un abrir y cerrar de ojos. Realizamos estos cálculos cientos de veces, de modo que al final nos ahorramos días e incluso semanas de tiempo de cálculo. Y tenemos la oportunidad de interactuar mejor con el sistema, para que las horas de trabajo de los científicos sean más rentables ... no se quedan inactivos, esperando el final de los cálculos. Todo está listo a la vez.
Debo enfatizar que estas mejoras de rendimiento épicas no requieren la compra de ningún sistema informático grande, paralelización o mayor complejidad ... de hecho, ¡el código de algoritmo más rápido es aún más simple y más versátil! Esta victoria completa se logra simplemente leyendo atentamente el código y teniendo algún conocimiento de la complejidad computacional.