Le problème est qu'il est incroyablement facile d'écrire des instructions qui ne signifient pas ce que vous pensez. Par exemple, un programme peut être complètement défiant de l'interprétation informatique. Outre,il n'y a littéralement aucun moyen de savoir si un programme se terminera du tout sans l'exécuter. Et il existe de très nombreuses façons de ralentir considérablement l'exécution d'un programme. Je veux dire ... vraiment ralentir. Ralentissez pour qu'il vous prenne toute votre vie ou plus. Cela se produit le plus souvent avec des programmes écrits par des personnes sans formation en informatique, c'est-à-dire des scientifiques d'autres domaines. Mon travail consiste à réparer ces programmes.
Les gens ne comprennent pas que l'informatique vous enseigne la théorie du calcul, la complexité des algorithmes, la calculabilité (c'est-à-dire, pouvons-nous vraiment calculer quelque chose? Trop souvent nous prenons pour acquis que nous le pouvons!) L'informatique fournit les connaissances, la logique et les méthodes d'analyse pour aider à écrire code qui s'exécutera en un minimum de temps ou avec une utilisation minimale des ressources.
Permettez-moi de vous montrer un exemple d'une énorme optimisation dans un script simple écrit par un de mes collègues.
Nous évoluons beaucoup en climatologie. Nous prenons des mesures de température et de précipitations à partir d'un modèle climatique mondial à grande échelle et les comparons à une grille locale à petite échelle. Disons que la grille globale est de 50x25 et la grille locale est de 1000x500. Pour chaque cellule de la grille de la grille locale, nous voulons savoir à quelle cellule de la grille globale elle correspond.
Une manière simple de représenter le problème est de minimiser la distance entre L [n] et G [n]. Il s'avère une telle recherche:
L[i]:
G[j]:
L[i] G[j]
L[i] * G
Cela semble assez simple. Cependant, si vous regardez de plus près, nous faisons beaucoup de travail supplémentaire. Regardez l'algorithme en termes de taille de l'entrée.
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)Le code ressemble à ceci:
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)Cela ressemble à un algorithme simple. Il est «facile» de calculer les distances puis de trouver le minimum. Mais dans une telle implémentation, à mesure que le nombre de cellules locales augmente, notre coût de calcul augmente de son produit avec le nombre de cellules dans la grille globale. Les données canadiennes ANUSPLIN contiennent 1068 × 510 cellules (544 680 au total). Disons que notre GCM contient 50x25 cellules (1250 au total). Ainsi, le coût du cycle interne en «quelques unités de calcul» est:
où les membres sont des constantes correspondant au coût du calcul de la distance entre deux points, de la recherche du point minimum et de la recherche de l'indice du tableau. En fait, nous ne nous soucions pas (beaucoup) des membres constants car ils ne dépendent pas de la taille de l'entrée. Vous pouvez donc simplement les additionner et estimer le coût de calcul;
Donc, pour cet ensemble d'entrées, notre coût est .
680 millions.
Cela semble beaucoup ... mais, qui sait? Les ordinateurs sont rapides, non? Si nous exécutons une implémentation naïve, elle fonctionnera en réalité un peu plus vite que 1668 secondes, soit un peu moins d'une demi-heure.
> 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 Mais le programme doit -il fonctionner pendant 30 minutes? C'est le problème. Nous comparons deux grilles, chacune contenant une tonne de structures que nous n'avons en aucun cas utilisées. Par exemple, les latitudes et les longitudes des deux grilles sont classées par ordre. Par conséquent, pour trouver un nombre, vous n'avez pas besoin de parcourir l'ensemble du tableau. Vous pouvez utiliser l'algorithme de réduction de moitié - regardez le point au milieu et décidez quelle moitié du tableau nous voulons. Ensuite, la recherche de l'espace entier prendra le logarithme de base deux de tout l'espace de recherche.
Une autre structure importante que nous n'avons pas utilisée est le fait que les latitudes vont en dimension et longitudes - en dimension . Ainsi, au lieu d'effectuer l'opération fois nous pouvons l'exécuter fois. C'est uneénormeoptimisation.
À quoi cela ressemble-t-il dans le pseudocode?
local[x]:
bisect_search(local[x], Global[x])
local[y]:
bisect_search(local[y], Global[y])
2d Dans du code:
## 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))
}Le coût de chaque recherche de bissection est le journal de la taille d'entrée. Cette fois, notre taille d'entrée est divisée en espaces X et Y, nous allons donc utiliser et pour Global, Local, X et Y.
Le coût est de 553 076. Eh bien, 553k sonne bien mieux que 680m. Comment cela affectera-t-il l'exécution?
> 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 seconde. Ce qui prenait auparavant près d'une demi-heure prend maintenant un peu plus d'un dixième de seconde.
> 1668.868 / .117
[1] 14263.83Soooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo, mais même moi , il surprend et impressionne combien l'accélération est. Environ 14 mille fois .
Auparavant, le script était si lent qu'il enregistrait le résultat sur le disque pour un examen manuel par des scientifiques avant de continuer. Désormais, tout est calculé en un clin d'œil. Nous effectuons de tels calculs des centaines de fois, ce qui nous permet de gagner des jours et même des semaines de temps de calcul. Et nous avons la possibilité de mieux interagir avec le système, afin que le temps de travail des scientifiques soit plus rentable ... ils ne restent pas les bras croisés, attendant la fin des calculs. Tout est prêt à la fois.
Je dois souligner que ces améliorations de performances épiques ne nécessitent pas l'achat de gros systèmes informatiques, la parallélisation ou une complexité accrue ... en fait, le code d'algorithme plus rapide est encore plus simple et plus polyvalent! Cette victoire complète a été obtenue simplement en lisant attentivement le code et en ayant une certaine connaissance de la complexité informatique.