11-algorithms-ja.Rmd

# スクリプト、アルゴリズム、関数  {#algorithms}

```{r, include=FALSE}
source("code/before_script.R")
```

## 必須パッケージ  {- #prerequisites-11}

この章では、主に Base R を使用するため、必要なソフトウェアは最小限である。
これから開発するアルゴリズムの結果を確認するために **sf**\index{sf} パッケージだけを使用する。
Chapter \@ref(spatial-class) で紹介した地理クラスと、それを使ってさまざまな入力ファイル形式を表現する方法 (Chapter \@ref(read-write) 参照) について理解していることを前提にしている。

## イントロダクション  {#intro-algorithms}

Chapter \@ref(intro) は、ジオコンピュテーションは既存のツールを使うだけでなく、「共有可能な R スクリプトや関数の形で」新しいツールを開発することが重要であることを示した。
本章では、これらの再現性のあるコードの構成要素について学ぶ。
また、Chapter \@ref(gis) で使用されているような低レベルの幾何学的アルゴリズムも紹介する。
これを読めば、このようなアルゴリズムの仕組みを理解し、複数のデータセットに対して、多くの人が、何度も使えるようなコードを書くことができるようになるはずである。
本章だけでは、熟練したプログラマになることはできない。
プログラミングは難しく、十分な練習が必要である [@abelson_structure_1996] :

> プログラミングをそれ自体の知的活動として理解するためには、プログラミングに目を向けなければならないし、プログラムを読み、書かなければならない。

しかし、プログラミングを学ぶ強い理由がある。
この章では、プログラミングそのものを教えるわけではない。プログラミングについては、@wickham_advanced_2019、@gillespie_efficient_2016、@xiao_gis_2016 を推奨する。これらの書籍は R や他の言語について教えてくれる。また、地理データに焦点を当て、プログラミング能力を伸ばすための基礎を作ることができる。

本章は、再現性\index{さいげんせい@再現性}の重要性について、例を示しながら強調していきたい。
再現性の利点は、他の人があなたの研究を複製することを可能にするだけではない。
再現性のあるコードは、一度だけ実行されるように書かれたコードよりも、計算効率、スケーラビリティ (より大きなデータセットに対して実行するコードの能力)、適応やメンテナンスのしやすさなど、あらゆる面で優れていることが多いのである。

スクリプトは、再現可能な R コードの基礎であり、このトピックは、Section \@ref(scripts) でカバーされている。
アルゴリズムは、Section \@ref(geometric-algorithms) で説明されているように、一連のステップを使用して入力を変更し、その結果、出力を得るためのレシピである。
共有と再現を容易にするために、アルゴリズムを関数に配置することができる。
それが、Section \@ref(functions) のトピックである。
ポリゴンの重心\index{じゅうしん@重心}を求める例で、これらの概念を結びつけていこう。
Chapter \@ref(geometry-operations) で、重心の関数\index{じゅうしん@重心} `st_centroid()` をすでに紹介したが、この例は、一見単純な操作が比較的複雑なコードの結果であることを強調し、次の観察を保証する [@wise_gis_2001] 。

> 空間データの問題で最も興味深いのは、人間にとっては些細なことに見えることが、コンピュータにとっては驚くほど難しいということである。

この例は、@xiao_gis_2016 にならい、「世の中にあるものを複製するのではなく、世の中のものがどのように機能しているかを示す」という本章の第二の目的も反映している。

## スクリプト  {#scripts}

パッケージで配布される関数が R コードの構成要素だとすれば、スクリプトはそれらを論理的な順序でまとめる接着剤となる。
スクリプトとは、再現可能なワークフローを作り出す目的で、手動または **targets** などの自動化ツールで保存・実行される [@landau_targets_2021]。\index{さいげんせい@再現性}
プログラミングの初心者にとってスクリプトは敷居が高く聞こえるかもしれないが、単なるプレーンテキストファイルである。
スクリプトは、通常はその言語を表す拡張子で保存される。例えば、Python は `.py`、Rust は `.rs` である。
R スクリプトは一般に、`.R` 拡張子で保存され、実行内容を反映した名前が付けられる。
例として、この本のリポジトリの `code` フォルダに格納されているスクリプトファイル [`11-hello.R`](https://github.com/geocompx/geocompr/tree/main/code/11-hello.R) がある。
`11-hello.R` は、次の 2 行のコードが含まれているだけの簡単なスクリプトで、そのうち 1 行はコメントである。

```r
# Aim: provide a minimal R script
print("Hello geocompr")
```

このスクリプトの中身は、とりたててスゴいものではないが、「スクリプトは複雑である必要はない」という点を示している。
保存されたスクリプトは、`source()` を使って、その全体を呼び出したり、実行したりすることができる。
このコマンドの出力から、コメント行は無視され、`print()` コマンドが実行されることがわかる。

```{r 10-algorithms-1}
source("code/11-hello.R")
```

以下のように `bash` や `PowerShell` などのシステムコマンドラインシェルから R スクリプトを呼び出すこともできる。

```bash
Rscript code/11-hello.R
```

`RScript` 実行可能ファイルが [設定](https://www.reddit.com/r/Rlanguage/comments/zaovly/is_anybody_able_to_run_a_r_script_in_powershell/)されていれば、システムのシェルで `Hello geocompr` と表示される。
スクリプトファイルに何を書くべきで何を書くべきではないかについて厳密なルールがあるわけではない。壊れた再現性のないコードになることもよくあるので、テストが必要である。
有効な R を含まないコード行は、エラーを防ぐため、行頭に `#` を追加してコメントアウトする必要がある。`11-hello.R` スクリプトの 1 行目を参照。
守るべき基本的ルールもある。

- 順番に書く。映画の脚本と同じように、スクリプトも「設定」「データ処理」「結果保存」といった明確な順番が必要である (映画でいうところの「始まり」「中間」「終わり」にほぼ相当する)
- 他の人 (と未来の自分) が理解できるように、スクリプトにコメントを追加してみよう。
これは、例えば、RStudio\index{RStudio} で、「折りたたみ可能な」コードセクションの見出しを作成するショートカット `Ctrl+Shift+R` を使って行うことができる
- 特に、スクリプトは再現可能であるべきである。どんなコンピュータでも動作する自己完結型のスクリプトは、調子の良い日に自分のコンピュータでしか動作しないスクリプトよりも有用である。
これには、必要なパッケージを最初に添付し、データを永続的なソース (信頼できるウェブサイトなど) から読み込み、前のステップが実行されたことを確認することが含まれる。^[
前のステップは、コメントまたは `if(!exists("x")) source("x.R")` のような if 文で参照できる (オブジェクト `x` が見つからない場合、スクリプトファイル `x.R` を実行する)。
]

パッケージ化しないかいぎり R スクリプトで再現性を強制するのは難しいが、それを助けるツールはある。
RStudio\index{RStudio} は、デフォルトで R スクリプトを「コードチェック」し、不具合のあるコードに赤い波線を引く (下図参照)。
**repex** パッケージは、再現性のためのツールである。

```{r codecheck, echo=FALSE, fig.cap="RStudio でのコード確認の様子。この例は 11-centroid-alg.R スクリプトの 19 行目のカッコが閉じられていないことを示している。", fig.scap="Illustration of 'code checking' in RStudio."}
knitr::include_graphics("images/codecheck.png")
```

```{block2 spellcheck, type='rmdnote'}
再現性に役立つのは **reprex** パッケージである。
このパッケージの `reprex()` 関数は、R コードが再現可能かどうかを確認し、GitHub などのサイトで交流する際に使えるマークダウンを生成する。
詳細は reprex.tidyverse.org を参照。
```

\index{さいげんせい@再現性}
このセクションの内容は、あらゆるタイプの R スクリプトに適用できる。
ジオコンピュテーションのためのスクリプトで特に考慮すべき点は、GDAL など外部ライブラリへの依存が多い点である。実際、Chapter \@ref(read-write) のデータ入出力では GDAL をたくさん使用した。 
GIS ソフトウェアの依存は、Chaptger \@ref(gis) で解説したようにより多くの特別なジオアルゴリズムを実行するために必要になる。
地理データを扱うスクリプトは、入力データセットが特定のファイル形式であることを必要とする。
このような依存関係は、スクリプトのコメントとして、またはスクリプトの一部であるプロジェクトの他の場所でコメントするか、**renv** や Docker などで依存性として記述するべきである。

「保守的」プログラミング技術と適切なエラーメッセージは、要件が満たされていない時に依存性を確認し、ユーザと対話する時間を節約する。
R では `if ()` で表される If 文を使用し、特定の条件が揃った時のみにメッセージを送ったりコードを実行するようにコードを書く。
以下のコード例は、特定のファイルがない場合にユーザにメッセージを送る (訳註: メッセージに日本語を使わないことを推奨する。)。

```{r}
if (!file.exists("required_geo_data.gpkg")) {
  message("No file, required_geo_data.gpkg is missing!")
} 
```

このスクリプトが行う作業は、以下の再現例で示されている。このスクリプトは、`poly_mat` という、長さ 9 単位の辺を持つ正方形 (この意味は次のセクションで明らかになる) という前提条件のオブジェクトに対して動作する。
この例では、インターネットに接続していることを前提に、`source()` が URL (ここでは短縮版を使用) で動作することを示している。
そうでない場合は、[github.com/geocompx/geocompr](https://github.com/geocompx/geocompr) からダウンロードできる `geocompr` フォルダのルートディレクトリから R を実行していると仮定して、`source("code/11-centroid-alg.R")` で同じスクリプトを呼び出すことができる。

```{r 10-algorithms-2, eval=FALSE}
poly_mat = cbind(
  x = c(0, 9, 9, 0, 0),
  y = c(0, 0, 9, 9, 0)
)
# geocompr リポジトリの code/11-centroid-alg.R の短い URL
source("https://t.ly/0nzj")
```

```{r 10-algorithms-3, echo=FALSE, fig.show='hide'}
poly_mat = cbind(
  x = c(0, 9, 9, 0, 0),
  y = c(0, 0, 9, 9, 0)
)
if(curl::has_internet()) {
  source("https://raw.githubusercontent.com/geocompx/geocompr/main/code/11-centroid-alg.R")
} else {
  source("code/11-centroid-setup.R")
}
```

## 幾何学的アルゴリズム  {#geometric-algorithms}

アルゴリズム\index{あるごりずむ@アルゴリズム}は、コンピュータにおける料理のレシピに相当するものと理解することができる。
入力に対して実行されると、有用または美味しい出力が得られる指示の完全なセットである。
入力とは、料理においては小麦粉や砂糖などの食材で、アルゴリズムの場合はデータと入力パラメータとなる。
美味しいケーキはレシピの結果であるのと同様に、アルゴリズムの成功は環境や社会的に利点のある出力となりうる。
具体的なケーススタディに入る前に、アルゴリズムとスクリプト (Section \@ref(scripts) ) や関数 (Section \@ref(functions) で説明するように、アルゴリズムを一般化し、他でも使えるようにしたり簡単にする) の関係について簡単に説明する。

「アルゴリズム」\index{あるごりずむ@アルゴリズム}という言葉は、西暦825年のバグダッドで出版されたされた、初期の数学教科書に端を発している。
この本はラテン語に翻訳され、人気を博しただけでなく、著者の名字である [al-Khwārizmī](https://en.wikipedia.org/wiki/Muhammad_ibn_Musa_al-Khwarizmi) が「科学用語として不滅の名声を得て、Alchoarismi、Algorismi、そして最終的には algorithm になった」 [@bellos_alex_2011] 。
コンピュータ時代には、アルゴリズム\index{あるごりずむ@アルゴリズム}は、問題を解決する一連のステップを指し、その結果、あらかじめ定義された出力が得られる。
入力は、適切なデータ構造で正式に定義されなければならない [@wise_gis_2001]。
アルゴリズムは、多くの場合、コードで実装される前に、処理の目的を示すフローチャートや疑似コード\index{ぎじこーど@疑似コード}として開始される。
使い勝手をよくするために、一般的なアルゴリズムは関数内にパッケージ化されていることが多く、(関数のソースコードを見ない限り) 一部または全部の手順が隠されている場合がある (Section \@ref(functions) を参照)。

Chapter \@ref(gis) で見たような ジオアルゴリズム\index{あるごりずむ@アルゴリズム}は、地理的なデータを取り込み、一般的には地理的な結果を返すアルゴリズムである (同じものを表す別の用語として、<u>GIS アルゴリズム</u>、<u>幾何学的アルゴリズム</u>がある)。
簡単そうに聞こえるだろうが、このテーマは奥が深く、<u>計算幾何学</u>という学問分野全体がその研究に専念している [@berg_computational_2008]。このテーマに関する書籍も多数出版されている。
例えば、@orourke_computational_1998 は、再現可能で自由に利用できる C コードを用いて、徐々に難しくなる幾何学的アルゴリズムの範囲を紹介している。

幾何学的アルゴリズムの例としては、ポリゴンの重心\index{じゅうしん@重心}を求めるものがある。
重心\index{じゅうしん@重心}の計算には多くのアプローチがあり、中には特定のタイプの[空間データ](https://en.wikipedia.org/wiki/Centroid)に対してのみ機能するものもある。
本節では、視覚化しやすいアプローチを選択する。ポリゴンを多くの三角形に分割し、それぞれの重心\index{じゅうしん@重心}を求める。このアプローチは、他の重心アルゴリズムと並んで @kaiser_algorithms_1993 によって議論された [簡単な説明は @orourke_computational_1998]。
コードを書く前に、このアプローチをさらに個別のタスクに分解するのに役立つ (以降、ステップ 1 からステップ 4 と呼ぶが、これらは模式図や疑似コード\index{ぎじこーど@疑似コード} として提示すこともできる)。

1. ポリゴンを連続した三角形に分割する
2. 各三角形の重心\index{じゅうしん@重心}を求める
3. それぞれの三角形の面積を求める
4. 三角形の中心点の面積加重平均\index{じゅうしん@重心}を求める

一見、簡単そうに見えるが、言葉をコードに変換するには、入力に制約がある場合でも、試行錯誤を繰り返しながら作業を進める必要がある。
このアルゴリズムは、180°以上の内角を持たない<u>凸ポリゴン</u>に対してのみ動作し、星形は使用できない (パッケージの **decido** と **sfdct** は外部ライブラリを使用して非凸ポリゴンを三角測量できる。[geocompx.org](https://geocompx.org/) の [algorithm](https://geocompx.github.io/geocompkg/articles/algorithm.html) vignetteに示されている。)。 

ポリゴンの最も単純なデータ構造は、x と y の座標の行列で、各行はポリゴンの境界を順にたどる頂点を表し、最初と最後の行は同一である [@wise_gis_2001]。
今回は、*GIS Algorithms* [@xiao_gis_2016 Python コードは [github.com/gisalgs](https://github.com/gisalgs/geom) を参照] の例を参考に、Figure \@ref(fig:polymat) に示すように、5 つの頂点を持つポリゴンを Base R で作成する。

```{r centroid-setup, echo=FALSE, eval=FALSE}
# show where the data came from:
source("code/11-centroid-setup.R")
```

```{r 10-algorithms-4}
# ポリゴンを表現する行列を作成
x_coords = c(10, 20, 12, 0, 0, 10)
y_coords = c(0, 15, 20, 10, 0, 0)
poly_mat = cbind(x_coords, y_coords)
```

これで、例となるデータセットができたので、上記のステップ 1 に着手する準備が整った。
以下のコードでは、1 つの三角形 (`T1`) を作成して、この方法を示している。また、 [数式](https://math.stackexchange.com/a/1702606) $1/3(a + b + c)$ ($a$ から $c$ は三角形の頂点を表す座標) に基づいて重心\index{じゅうしん@重心}を計算するステップ 2 も示している。

```{r 10-algorithms-5}
# 原点を作成:
Origin = poly_mat[1, ]
# 三角形の行列を作成:
T1 = rbind(Origin, poly_mat[2:3, ], Origin) 
C1 = (T1[1,] + T1[2,] + T1[3,]) / 3
```

```{r, echo=FALSE}
# (Note: drop = FALSE preserves classes, resulting in a matrix):
C1_alternative = (T1[1, , drop = FALSE] + T1[2, , drop = FALSE] + T1[3, , drop = FALSE]) / 3
```

```{r polymat, echo=FALSE, fig.cap="ポリゴン重心計算問題。", fig.height="100", warning=FALSE}
# 最初のプロット。おそらく削除。
old_par = par(pty = "s") 
plot(poly_mat, cex = 3)
lines(poly_mat, lwd = 7)
lines(T1, col = "#fa8072", lwd = 2)
text(x = C1[1], y = C1[2], "C1", col = "#fa8072")
par(old_par)
```

ステップ 3 では、各三角形の面積を求めるので、大きな三角形の不釣り合いな影響を考慮した<u>加重平均</u>を計算する。 
三角形の面積を計算する公式は次の通りである [@kaiser_algorithms_1993]。

$$
\frac{A_x ( B_y − C_y ) + B_x ( C_y − A_y ) + C_x ( A_y − B_y )}{ 2 }
$$

ここで、$A$ から $C$ は三角形 `T1` の 3 点、$x$ と $y$ は x と y の次元を指す。
この式を、三角形の行列表現のデータを扱う R コードに翻訳すると、次のようになる (関数 `abs()` は、正の結果を保証する)。

```{r 10-algorithms-6}
# 行列 T1 で表される三角形の面積を計算
abs(T1[1, 1] * (T1[2, 2] - T1[3, 2]) +
    T1[2, 1] * (T1[3, 2] - T1[1, 2]) +
    T1[3, 1] * (T1[1, 2] - T1[2, 2])) / 2
```

このコードチャンクは正しい結果を出力する。^[
この結果は、以下の式で確認することができる (ベースが水平であると仮定している)。
面積は、ベース幅×高さの半分となる。$A = B * H / 2$.
この場合、$10 * 10 / 2 = 50$.
]
このコードは不格好で、別の三角行列で実行する場合、再入力しなければならない点が問題である。
より一般化するために、このコードを Section \@ref(functions) で関数に変換する方法を見てみよう。

ステップ 4 では、ステップ 2 と 3 を 1 つの三角形だけでなく、すべての三角形に対して行う必要がある (上の例)。
このため、ポリゴンを表すすべての三角形を作成するための<u>イテレーション</u> (繰り返し) が必要である。Figure \@ref(fig:polycent) に示す。
`lapply()`\index{るーぷしょり@ループ処理!lapply} と `vapply()`\index{るーぷしょり@ループ処理!vapply} が各三角形の反復処理に使われているのは、Base R で簡潔な解が得られるからである。^[
`vapply()`\index{るーぷしょり@ループ処理!vapply} のドキュメントについては `?lapply` を参照。イテレーションについては Chapter \@ref(location) を参照。
]

```{r 10-algorithms-7}
i = 2:(nrow(poly_mat) - 2)
T_all = lapply(i, function(x) {
  rbind(Origin, poly_mat[x:(x + 1), ], Origin)
})

C_list = lapply(T_all,  function(x) (x[1, ] + x[2, ] + x[3, ]) / 3)
C = do.call(rbind, C_list)

A = vapply(T_all, function(x) {
  abs(x[1, 1] * (x[2, 2] - x[3, 2]) +
        x[2, 1] * (x[3, 2] - x[1, 2]) +
        x[3, 1] * (x[1, 2] - x[2, 2]) ) / 2
  }, FUN.VALUE = double(1))
```

```{r polycent, fig.cap="複数の三角による繰り返し重心アルゴリズム。繰り返し 2 と 3 における X は、面積加重重心。", fig.scap="Iterative centroid algorithm with triangles.", echo=FALSE, fig.asp=0.3}
# idea: show animated version on web version
source("code/11-polycent.R")
```

これで、ステップ 4 を完了し、総面積を `sum(A)` で計算し、ポリゴンの重心\index{じゅうしん@重心}座標を  `weighted.mean(C [, 1] , A)` と  `weighted.mean(C [, 2] , A)` でポリゴンの重心座標を計算する (注意深い読者のための練習: これらのコマンドが動作することを確認してみよう)。
アルゴリズム\index{あるごりずむ@アルゴリズム}とスクリプトの関連性を示すために、このセクションの内容を凝縮して `11-centroid-alg.R` とした。
Section \@ref(scripts) の最後で、このスクリプトが正方形の重心\index{じゅうしん@重心}を計算する方法を見た。
アルゴリズムを<u>スクリプト</u>化することの素晴らしい点は、新しい `poly_mat` オブジェクト上で動作することである (`st_centroid()` を参照してこれらの結果を検証するには、以下の演習を参照)。

```{r 10-algorithms-8}
source("code/11-centroid-alg.R")
```

上記の例では、低レベルの地理的な操作は、base R で第一原理から開発することが<u>できる</u>ことが示されている。
また、すでに試行錯誤したソリューションが存在する場合、車輪の再発明をする価値はないことも示している。
もし、ポリゴンの重心\index{じゅうしん@重心}を求めるだけなら、`poly_mat` を **sf** オブジェクトとして表現し、代わりに既存の `sf::st_centroid()` 関数を使用する方が早かっただろう。
しかし、第一原理でアルゴリズムを書くことの大きな利点は、プロセスのすべての段階を理解できることであり、他の人のコードを使うときには保証されないことである。
さらに考慮すべきは性能である。R は C++\index{C++} のような低レベルの言語と比較すると数値計算が遅いことが多く、最適化が困難である (Section \@ref(software-for-geocomputation) 参照)。
新しい手法の開発を目的とするのであれば、計算効率を優先させるべきではない。
これは、「早すぎる最適化はプログラミングにおける諸悪の根源 (あるいは少なくともそのほとんど)」という言葉に集約される [@knuth_computer_1974]。

アルゴリズム\index{あるごりずむ@アルゴリズム}開発は大変な作業である。
このことは、Base R\index{R} を使った重心\index{じゅうしん@重心}アルゴリズム\index{あるごりずむ@アルゴリズム}の開発に費やした作業量から明らかである。このアルゴリズムは、実世界での応用が限られている問題に対する一つの、むしろ非効率的なアプローチに過ぎない。というのも、実際には凸ポリゴンは珍しいからである。
この経験は、GEOS\index{GEOS} や CGAL\index{CGAL} (計算幾何学アルゴリズムライブラリ, Computational Geometry Algorithms Library) など、高速に動作し、かつ幅広い入力ジオメトリタイプに対応する低レベル地理ライブラリの理解につながるはずである。
このようなライブラリのオープンソース化の大きな利点は、そのソースコード\index{そーすこーど@ソースコード}が、研究、理解、(技術と自信があれば) 改変のために容易に利用できることである。^[
CGAL\index{CGAL} の関数 `CGAL::centroid()` は、https://doc.cgal.org/latest/Kernel_23/group__centroid__grp.html で説明しているように、実際には 7 つのサブ関数で構成されており、幅広い入力データ型に対応できるようになっているが、私たちが作成したソリューションは特定の入力データ型にのみ対応する。
GEOS\index{GEOS} の関数 `Centroid::getCentroid()` の基礎となるソースコードは、https://github.com/libgeos/geos/search?q=getCentroid で見ることができる。
]

## 関数  {#functions}

アルゴリズムと同様に \index{あるごりずむ@アルゴリズム}、関数は入力を受け取り、出力を返す。
しかし、関数\index{かんすう@関数}は、「レシピ」そのものではなく、特定のプログラミング言語での実装を指している。
R では、関数\index{かんすう@関数}はそれ自体がオブジェクトであり、モジュール方式で作成したり結合したりすることができる。
例えば、重心\index{じゅうしん@重心}生成アルゴリズム\index{かんすう@関数}のステップ 2 を引き受ける関数を以下のように作成することができる。

```{r 10-algorithms-9}
t_centroid = function(x) {
  (x[1, ] + x[2, ] + x[3, ]) / 3
}
```

上記の例では、[関数](https://adv-r.hadley.nz/functions.html)の 2 つの重要な構成要素を示している。(1) 関数の中身 (*body*) は、中括弧内のコードで、関数が入力に対して何をするかを定義する。 (2) 引数 (*argument*) は 、関数が扱う引数のリスト。この場合は `x` である (3 番目の重要なコンポーネント、環境はこのセクションの範囲外)。
デフォルトでは、関数は最後に計算したオブジェクトを返す (`t_centroid()` の場合は重心\index{じゅうしん@重心}の座標)。^[
また、関数の本文に `return(output)` を追加することで、関数の出力を明示的に設定することができる。 `output` は返すべき結果である。
]

```{r 10-algorithms-10, eval=FALSE, echo=FALSE}
body(t_centroid)
formals(t_centroid)
environment(t_centroid)
```

この関数は、前のセクションのポリゴンの例から最初の三角形の面積を計算する以下のコマンドのように、渡した入力に対して動作する (Figure \@ref(fig:polycent) を参照)。

```{r 10-algorithms-11}
t_centroid(T1)
```

また、三角形の面積を計算する関数\index{かんすう@関数}を作成することができる。ここでは、`t_area()` と名付ける。

```{r 10-algorithms-12}
t_area = function(x) {
  abs(
    x[1, 1] * (x[2, 2] - x[3, 2]) +
    x[2, 1] * (x[3, 2] - x[1, 2]) +
    x[3, 1] * (x[1, 2] - x[2, 2])
  ) / 2
}
```

なお、この関数を作成した後は、1 行のコードで三角形の面積を計算できるようになり、冗長なコードの重複を避けることができる。
関数は、コードを一般化するためのメカニズムである。
新たに作成した関数\index{かんすう@関数} `t_area()` は、これまで使ってきた「三角行列」データ構造と同じ寸法を持つと仮定した任意のオブジェクト `x` を受け取り、その面積を返すもので、`T1` で図示すると次のようになる。

```{r 10-algorithms-13}
t_area(T1)
```

関数\index{かんすう@関数}を使って、高さ 1、底辺 3 の新しい三角行列の面積を求めることで、その一般化可能性を検証することができる。

```{r 10-algorithms-14}
t_new = cbind(x = c(0, 3, 3, 0),
              y = c(0, 0, 1, 0))
t_area(t_new)
```

関数の便利な点は、モジュール化されていることである。
出力が何であるかが分かっていれば、ある関数を別の関数の構成要素として利用することができる。
したがって、関数 `t_centroid()` と `t_area()` は、スクリプト `11-centroid-alg.R` というより大きな関数\index{かんすう@関数}のサブコンポーネントとして使うことができる。このスクリプトは、任意の凸ポリゴンの面積を計算する。
以下のコードでは、凸ポリゴンに対する `sf::st_centroid()` の動作を模倣する関数 `poly_centroid()` を作成している。^[
なお、作成した関数は、`lapply()`\index{るーぷしょり@ループ処理!lapply}  と `vapply()`\index{るーぷしょり@ループ処理!vapply} の関数呼び出しで繰り返し呼び出されている。
]

```{r 10-algorithms-15}
poly_centroid = function(poly_mat) {
  Origin = poly_mat[1, ] # create a point representing the origin
  i = 2:(nrow(poly_mat) - 2)
  T_all = lapply(i, function(x) {rbind(Origin, poly_mat[x:(x + 1), ], Origin)})
  C_list = lapply(T_all, t_centroid)
  C = do.call(rbind, C_list)
  A = vapply(T_all, t_area, FUN.VALUE = double(1))
  c(weighted.mean(C[, 1], A), weighted.mean(C[, 2], A))
}
```


```{r 10-algorithms-16, echo=FALSE, eval=FALSE}
# 少し複雑な、出力を伴う関数
poly_centroid = function(poly_mat, output = "matrix") {
  Origin = poly_mat[1, ] # 原点を作成
  i = 2:(nrow(poly_mat) - 2)
  T_all = T_all = lapply(i, function(x) {
    rbind(Origin, poly_mat[x:(x + 1), ], Origin)
  })
  C_list = lapply(T_all, t_centroid)
  C = do.call(rbind, C_list)
  A = vapply(T_all, t_area, FUN.VALUE = double(1))
  centroid_coords = c(weighted.mean(C[, 1], A), weighted.mean(C[, 2], A))
  if (output == "matrix") {
    return(centroid_coords)
  } else if (output == "area")
    return(sum(A))
}
```

```{r 10-algorithms-17}
poly_centroid(poly_mat)
```

`poly_centroid()` などの関数\index{かんすう@関数}はさらに拡張して、さまざまなタイプの出力を提供することができる。
例えば、結果をクラス `sfg` のオブジェクトとして返すには、結果を返す前に、「ラッパー」関数を用いて `poly_centroid()` の出力を変更することができる。

```{r 10-algorithms-18}
poly_centroid_sfg = function(x) {
  centroid_coords = poly_centroid(x)
  sf::st_point(centroid_coords)
}
```

以下のように、`sf::st_centroid()` からの出力と同じであることが確認できる。

```{r 10-algorithms-19}
poly_sfc = sf::st_polygon(list(poly_mat))
identical(poly_centroid_sfg(poly_mat), sf::st_centroid(poly_sfc))
```

## プログラミング  {#programming}

この章では、スクリプトからアルゴリズム\index{あるごりずむ@アルゴリズム}という厄介なトピックを経由して関数へと移った。
抽象的な議論だけでなく、具体的な問題を解決するために、それぞれの実用例を作成した。

- スクリプト `11-centroid-alg.R` を導入し、「ポリゴンマトリックス」で実行した。
- このスクリプトを動作させるための個々のステップは、アルゴリズム\index{あるごりずむ@アルゴリズム}、つまり計算レシピとして記述した。
- アルゴリズムを一般化するために、このアルゴリズムをモジュール関数に変換し、最終的にそれらを組み合わせて、前節の関数 `poly_centroid()` を作成した。

これらのステップは、簡単なことである。
しかし、プログラミングの技術とは、スクリプト、アルゴリズム、関数を<u>組み合わせ</u>て、性能の良い<u>システム</u>にすることなのである。
できたものは、堅牢で、皆が簡単に使えるものであるべきである。
この本を読んでいるほとんどの人がそうであるように、プログラミングの初心者であれば、前のセクションの結果を追って再現できることは、大きな達成感を得ることができるはずである。 
プログラミングができるようになるまでには、何時間もかけて熱心に勉強し、練習する必要がある。

新しいアルゴリズム\index{あるごりずむ@アルゴリズム}を効率的に実装する際の課題は、実運用で使用することを意図していない単純な関数を作成するために費やした作業量を考慮すると、はるかに遠い。現在の状態では、`poly_centroid()` はほとんどの (非凸) ポリゴンで失敗する!
ここから生じる疑問は、関数\index{かんすう@関数}をどのように一般化するかということである。
(1) 非凸ポリゴンを三角測量する方法を探す (geocompx.github.io/geocompkg/articles/ のオンライン記事 [Algorithms Extended](https://geocompx.github.io/geocompkg/articles/algorithm.html) で扱っている話題) と (2) 三角メッシュに依存しない他の重心のアルゴリズムを調べるという 2 つの選択肢がある。

もっと大きな疑問は、高性能なアルゴリズムがすでに実装され、`st_centroid()` のような関数にパッケージされているのに、ソリューションをプログラミングする価値があるのだろうか、ということである。
この具体的なケースにおける還元論的な答えは「価値はない」である。
広い意味で、プログラミングを学ぶことの利点を考えると、答えは「場合による」である。
プログラミングでは、あるメソッドを実装しようとすると、すでに誰かがその苦労をしていることに気づき、何時間も無駄にすることがよくある。
本章は、ジオアルゴリズムのプログラミングへの第一歩として捉えることができる。
しかし、一般化された解決策をプログラムする場合と、既存の高水準の解決策を利用する場合の教訓とも言える。
新しい関数を作るのが最善の場合もあれば、すでにある関数を使うのが良い場合もある。

「車輪の再発明をするな」という言葉は、他の人生の歩みと同様、いやそれ以上にプログラミングに当てはまる。
プロジェクトの最初に少し調査して考えることで、プログラミングの時間をどこに費やすのがベストなのかを決めることができる。
また、以下の 3 つの原則は、簡単なスクリプトであれ、何百もの関数で構成されるパッケージであれ、コードを書くときに労力を最大限に活用するのに役立つ。

1. [DRY](https://en.wikipedia.org/wiki/Don%27t_repeat_yourself) (don't repeat yourself): コードの繰り返しを最小限に抑え、より少ないコード行数で特定の問題を解決することを目指す。
この原則は、「R for Data Science」の「Functions」の章において、コードの繰り返しを減らすための関数の使用を参照して説明されている [@grolemund_r_2016]。
2. [KISS](https://en.wikipedia.org/wiki/KISS_principle) (keep it simple stupid): この原則は、複雑な解決策よりも単純な解決策を最初に試し、必要に応じて依存関係を使用し、スクリプトを簡潔に保つことを目指すことを示唆。 
この原則は、「ものごとはできるかぎりシンプルにすべきだ。しかし、シンプルすぎてもいけない。 」"things should be made as simple as possible, but no simpler" という[名言](https://www.nature.com/articles/d41586-018-05004-4) のコンピュータ版である。
3. Modularity: コードを明確に分割することで、メンテナンスが容易になる。
関数は、たった一つのことをするだけにして、それに専念するべきである。
もし関数が長くなりすぎたら、それを複数の小さな関数に分割し、それぞれを別の目的に再利用することを考えよう。

この章だけで、すぐに完璧な関数を作成できるようになることは保証していない。
しかし、この章の内容は、いつ挑戦するのが適切かを判断するのに役立つと確信している (問題を解決する既存の関数がない場合、プログラミングタスクが自分の能力の範囲内にある場合、そのソリューションの利点が開発にかかる時間コストを上回ると思われる場合)。
上記の原則と、上記の例を通しての実践的な経験を組み合わせることで、スクリプト、パッケージ作成、プログラミングのスキルを手に入れることができる。
プログラミングへの最初の一歩は時間がかかるが (以下の演習は急がないように)、長期的な見返りは大きくなるだろう。

## 演習  {#ex-algorithms}

```{r, echo=FALSE, results='asis'}
res = knitr::knit_child('_11-ex-ja.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE))
cat(res, sep = '\n')
```