問題は、あなたが考えていることを意味しない指示を書くのは信じられないほど簡単だということです。たとえば、プログラムはコンピュータの解釈に完全に反抗する場合があります。その上、プログラムを実行せずにプログラムが終了するかどうかを文字通り判断する方法はありません。そして、プログラムの実行を大幅に遅くする方法はたくさんあります。つまり...本当に遅くなります。それがあなたの生涯以上かかるようにそれを遅くしてください。これは、コンピュータ教育を受けていない人々、つまり他の分野の科学者によって書かれたプログラムで最も頻繁に発生します。私の仕事はそのようなプログラムを修正することです。
コンピュータサイエンスが計算の理論、アルゴリズムの複雑さ、計算可能性を教えてくれることを人々は理解していません（つまり、実際に何かを計算できるのでしょうか？私たちができることを当然のことと思っていることが多すぎます！）最小限の時間または最小限のリソース使用で実行されるコード。
私の同僚が書いた1つの簡単なスクリプトでの大規模な最適化の例を示しましょう。
私たちは気候学で多くのスケーリングを行います。大規模な地球規模の気候モデルから温度と降水量の読み取り値を取得し、それらを細かいスケールのローカルグリッドに対してプロットします。グローバルグリッドが50x25で、ローカルグリッドが1000x500であるとします。ローカルグリッドの各グリッドセルについて、グローバルグリッドのどのグリッドセルに対応するかを知りたいと思います。
問題を表す簡単な方法は、L [n]とG [n]の間の距離を最小化することです。そのような検索が判明しました：
L[i]:
G[j]:
L[i] G[j]
L[i] * G
簡単そうです。ただし、よく見ると、余分な作業がたくさんあります。入力のサイズの観点からアルゴリズムを見てください。
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)コードは次のようになります。
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)単純なアルゴリズムのように見えます。距離を計算して最小値を見つけるのは「簡単」です。しかし、このような実装では、ローカルセルの数が増えると、グローバルグリッド内のセルの数との積によって計算コストが増加します。カナダのANUSPLINデータには、1,068×510セル（合計544,680）が含まれています。GCMに50x25セル（合計1250）が含まれているとしましょう。したがって、「いくつかの計算単位」の内部サイクルのコストは次のとおりです。
メンバーがいる場所 2点間の距離を計算し、最小点を見つけ、配列インデックスを見つけるコストに対応する定数です。実際、定数メンバーは入力のサイズに依存しないため、（あまり）気にしません。したがって、それらを合計して計算コストを見積もることができます。
したがって、この一連の入力の場合、コストは ..。
6億8000万。
それはたくさんのようです...しかし、誰が知っていますか？コンピューターは速いですよね？単純な実装を実行すると、実際には1668秒より少し速く実行されます。これは30分弱です。
> 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 しかし、プログラムは30分間実行する必要がありますか？それが問題です。2つのメッシュを比較しています。各メッシュには、まったく使用していない構造がたくさんあります。たとえば、両方のグリッドの緯度と経度は並べ替えられた順序になっています。したがって、番号を見つけるために、配列全体を調べる必要はありません。半分にするアルゴリズムを使用できます。中央のポイントを見て、配列の半分を決定します。次に、スペース全体を検索すると、検索スペース全体の2つの対数になります。
私たちが使用しなかったもう1つの重要な構造は、緯度が次元に入るという事実です。、および経度-次元内 ..。したがって、操作を実行する代わりに できるから 時間。これは大きな最適化です。
疑似コードではどのように見えますか？
local[x]:
bisect_search(local[x], Global[x])
local[y]:
bisect_search(local[y], Global[y])
2d コード内：
## 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))
}各二分検索のコストは、入力サイズのログです。今回は入力サイズがXスペースとYスペースに分割されているので、 そして グローバル、ローカル、XおよびYの場合。
コストは553,076です。まあ、553kは680mよりずっといいですね。これはランタイムにどのように影響しますか？
> 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秒。以前は30分近くかかっていましたが、今では10分の1秒強かかります。
> 1668.868 / .117
[1] 14263.83Soooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo、それでも私は高速化がどのくらいで驚いと感心します。約14000回。
以前は、スクリプトが非常に遅いため、続行する前に科学者が手動で確認できるように結果をディスクに保存していました。今、すべてが瞬く間に計算されます。このような計算は何百回も実行されるため、最終的には数日から数週間の計算時間を節約できます。そして、私たちはシステムとよりよく相互作用する機会を得て、科学者の労働時間がより有益になるようにします...彼らは計算の終了を待ってアイドル状態になりません。すべてが一度に準備ができています。
これらの壮大なパフォーマンスの向上は、大規模なコンピューターシステムの購入、並列化、または複雑さの増大を必要としないことを強調する必要があります...実際、より高速なアルゴリズムコードはさらに単純で、より用途が広いです！この完全な勝利は、コードを注意深く読み、計算の複雑さについてある程度の知識を持っているだけで達成されました。