14.000x Akselerasi atau Kemenangan Ilmu Komputer

Sebagai pengembang perangkat lunak ilmiah, saya melakukan banyak pemrograman. Dan kebanyakan orang di bidang ilmiah lainnya cenderung berpikir bahwa pemrograman "hanya" membuang kode dan menjalankannya. Saya memiliki hubungan kerja yang baik dengan banyak kolega, termasuk yang berasal dari negara lain ... Fisika, klimatologi, biologi, dll. Tetapi dalam hal pengembangan perangkat lunak, Anda mendapatkan kesan yang berbeda bahwa mereka berpikir: "Hei, apa yang bisa sulit?! Kami hanya menuliskan beberapa instruksi tentang apa yang harus dilakukan komputer, tekan tombol Run dan selesai - kami mendapatkan jawabannya! "



Masalahnya adalah, sangat mudah untuk menulis instruksi yang tidak berarti apa yang Anda pikirkan. Misalnya, program mungkin sepenuhnya menentang interpretasi komputer. Selain,secara harfiah tidak ada cara untuk mengetahui apakah suatu program akan berhenti sama sekali tanpa menjalankannya. Dan ada banyak sekali cara untuk memperlambat eksekusi program. Maksudku ... sangat pelan. Perlambatlah agar itu mengambil seluruh hidup Anda atau lebih. Ini paling sering terjadi dengan program yang ditulis oleh orang-orang tanpa pendidikan komputer, yaitu ilmuwan dari bidang lain. Tugas saya adalah memperbaiki program semacam itu.



Orang-orang tidak memahami bahwa ilmu komputer mengajarkan Anda teori komputasi, kompleksitas algoritme, komputasi (yaitu, apakah kita benar-benar dapat menghitung sesuatu? Terlalu sering kita menerima begitu saja bahwa kita bisa!) Ilmu komputer menyediakan pengetahuan, logika dan metode analisis untuk membantu menulis kode yang akan dijalankan dalam waktu singkat atau dengan penggunaan sumber daya minimal.



Izinkan saya menunjukkan kepada Anda contoh pengoptimalan besar dalam satu skrip sederhana yang ditulis oleh rekan saya.



Kami banyak meningkatkan skala dalam klimatologi. Kami mengambil pembacaan suhu dan curah hujan dari model iklim global skala besar dan memplotnya dengan jaringan lokal skala halus. Misalkan kisi global adalah 50x25 dan kisi lokalnya adalah 1000x500. Untuk setiap sel grid di grid lokal, kita ingin tahu sel grid mana di grid global yang bersangkutan.



Cara mudah untuk merepresentasikan soal adalah dengan memperkecil jarak antara L [n] dan G [n]. Ternyata pencarian seperti itu:

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

Sepertinya cukup sederhana. Namun, jika Anda melihat lebih dekat, kami melakukan banyak pekerjaan ekstra. Lihatlah algoritme dalam hal ukuran input.

    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)

Kode tersebut terlihat seperti ini:

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)

Sepertinya algoritme sederhana. Sangatlah mudah untuk menghitung jarak dan kemudian mencari nilai minimumnya. Namun dalam penerapan seperti itu, seiring dengan bertambahnya jumlah sel lokal, biaya komputasi kami meningkat berdasarkan produknya dengan jumlah sel di jaringan global. Data ANUSPLIN Kanada berisi 1068 × 510 sel (total 544680). Mari kita asumsikan GCM kita berisi 50x25 sel (total 1250). Jadi, biaya siklus dalam di "beberapa unit komputasi" adalah:

$$

$$(c_0 \cdot L \times G) + (c_1 \cdot L \times G) + (c_2 \cdot L \times G)$$

dimana anggota $$$c$Adalah konstanta yang sesuai dengan biaya penghitungan jarak antara dua titik, mencari titik minimum, dan mencari indeks larik. Faktanya, kami tidak (banyak) peduli tentang anggota konstan karena mereka tidak bergantung pada ukuran input. Jadi Anda bisa menambahkannya dan memperkirakan biaya komputasi;

$$

$$(c \cdot L \times G)$$

Jadi untuk rangkaian input ini, biaya kami adalah $$$544,680 \times 1,250 = 680,850,000$...



680 juta.



Sepertinya banyak ... meskipun, siapa yang tahu? Komputer itu cepat, bukan? Jika kita menjalankan implementasi yang naif, pada kenyataannya itu akan berjalan sedikit lebih cepat dari 1668 detik, yang kurang dari setengah jam.

> 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 

Tetapi apakah program harus berjalan selama 30 menit? Itulah masalahnya. Kami membandingkan dua kisi, yang masing-masing memiliki banyak struktur yang belum kami gunakan sama sekali. Misalnya, lintang dan bujur di kedua kisi diurutkan secara berurutan. Oleh karena itu, untuk mencari angka, Anda tidak perlu menelusuri seluruh larik. Anda dapat menggunakan algoritma halving - lihat titik di tengah dan putuskan bagian mana dari array yang kita inginkan. Kemudian mencari seluruh ruang akan mengambil basis dua logaritma dari seluruh ruang pencarian.



Struktur penting lainnya yang tidak kami gunakan adalah fakta bahwa garis lintang memiliki dimensi$$$x$, dan bujur - dalam dimensi $$$y$... Jadi, alih-alih melakukan operasi$$$x \times y$ karena kita bisa melakukannya $$$x + y$waktu. Ini adalah pengoptimalan yang sangat besar .



Seperti apa tampilan di pseudocode?

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

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

 2d      

Dalam kode:

## 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))
}

Biaya setiap pencarian pembagian adalah log dari ukuran input. Kali ini ukuran masukan kami dibagi menjadi ruang X dan Y, jadi kami akan menggunakan$$$G_x, G_y, L_x$ dan $$$L_y$ untuk Global, Lokal, X dan Y.

$$

$$cost = L_x \times log_2 G_x + L_y \times log_2 G_y + L_x \times L_y$$

Biayanya adalah 553.076. Nah, 553k terdengar jauh lebih baik daripada 680m. Bagaimana ini akan mempengaruhi runtime?

> 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 detik. Apa yang biasanya memakan waktu hampir setengah jam sekarang hanya membutuhkan sepersepuluh detik.

> 1668.868 / .117
[1] 14263.83

Soooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo, tapi bahkan saya terkejut dan terkesan oleh berapa banyak speedup ini. Kira-kira 14 ribu kali .



Sebelumnya, skrip tersebut sangat lambat sehingga menyimpan hasilnya ke disk untuk ditinjau secara manual oleh para ilmuwan sebelum melanjutkan. Sekarang semuanya dihitung dalam sekejap mata. Kami melakukan perhitungan seperti itu ratusan kali, sehingga kami menghemat waktu berhari-hari bahkan berminggu-minggu. Dan kami mendapat kesempatan untuk berinteraksi lebih baik dengan sistem, sehingga jam kerja para ilmuwan lebih menguntungkan ... mereka tidak duduk diam, menunggu akhir penghitungan. Semuanya siap sekaligus.



Saya harus menekankan bahwa peningkatan kinerja yang luar biasa ini tidak memerlukan pembelian sistem komputer besar, paralelisasi, atau peningkatan kompleksitas ... pada kenyataannya, kode algoritme yang lebih cepat bahkan lebih sederhana dan lebih fleksibel! Kemenangan lengkap ini dicapai hanya dengan membaca kode secara cermat dan memiliki pengetahuan tentang kompleksitas komputasi.