この記事は「R Advent Calendar 2017」の22日目です。

前書き

ingressという位置情報ゲームの 完全多重フィールド を探し出すプログラムをRでつくりました、というお話です。
ingressで完全多重を作りたい人以外に需要があるのかわかりませんが、何かが誰かの役に立つかもしれないので、とりあえずプログラムを公開して、その説明を書き残すことにしました。

書いてみるとそこそこの長文になってしまいました。目次的なものを書いておきます。

はじめに
ingressって何？完全多重って何？などなど

とりあえず動かす
サンプルデータを使って4重の完全多重を探す方法を説明しています。

ポータルの一覧とその座標を取得する方法
IITCプラグインを使ってポータルの緯度経度を取得する方法です。

完全多重が構成できるかを調べる方法
メイン処理のmlcf.Rの処理内容の説明です。Rの処理で工夫したことなども書いてあります。

完全多重を構築する手順と必要鍵数を考える方法
mlcf4.path.Rの処理内容の説明です。できるだけ短い距離でポータルを辿るルートを探しています。

均一５重の構成・手順を求める方法
5重の均一多重を探す方法を4重との違いを中心に説明しています。

Rのadvent calendarにエントリしたのに、あまりRについて語れませんでした。位置情報を扱うRのプログラムの一例の紹介だと思ってください。

はじめに

ingress紹介

Niantic社がポケモンGO以前から公開している位置情報ゲームです。
https://www.ingress.com/

プレイヤーはゲーム開始時に２つの陣営（青と緑）のどちらに属するかを選びます。
相手陣営のポータル（ポケモンGOのポケスポットに相当）に行って、そのポータルを攻撃して自陣営のものにします。
自陣営のポータル同士を「リンク」という仮想の線で結ぶことができます。
３点をリンクで結んでできる三角形を「コントロールフィールド」（CFやフィールドとも呼ばれる、以下フィールド）と呼びます。

で、世界各地に作られたフィールドの面積の陣営ごとの合計を日々競い合っています。

交叉するリンクは作れない、フィールドの中のポータルからはリンクを結ぶことができない、など制約がありますがそのあたりは割愛。

詳しくは「ingress 入門」で検索して出てくる記事を参照してください。
https://www.google.co.jp/search?q=ingress+入門

多重フィールド

ingressにはいろんな楽しみ方があって、その一つに多重フィールドがあります。
フィールドを包含するフィールドを作る、ということを繰り返していくと、複数のフィールドが重なった状態（多重フィールド）ができあがります。
普通に多重フィールドを作っていくと、一部の多重度が高く、その他の多重度が低い、という状態になります。（画面上で見ると、中心部の色が濃くて、周辺部に行くほど色が薄くなります）

完全多重フィールド

ある一定の手順でフィールドを作ると、フィールド内の多重度を均一にすることができて、これを 均一多重フィールド と呼びます。
さらに均一多重フィールドの内部に、使用されていないポータルがない状態を 完全多重フィールド と呼びます。

完全多重の説明は http://yugioh-hack.com/ingress-perfect-cf/ が詳しいです。

この完全多重を構成するポータルの配置を人力で探し出すのはかなり難しいです。
そもそも5重の完全多重自体が極めて少ないです。（都内近郊を見渡しても10パターン以下です）

仮に見つかったとしても完全多重フィールドを作るための手順を考えるのが、かなり面倒だったりします。

このプログラムで何ができるの？

ポータルの位置情報（緯度経度）を読み込んで、 完全多重フィールドのポータル構成を探し出します。見つかったポータル構成をgooglemapに表示したり、 完全多重フィールドの構築手順と必要な鍵数を計算 します。

なぜRで作ったの？

Rで作る必然性は全くなく、pythonやjavaやC++でもよかったのです。なんとなくRでした。

とりあえず動かす

ソースコードなど

ソースとサンプルデータをGitHubにおいています。
https://github.com/mktkbt/img_mlcf

ファイル	説明
mlcf.R	完全多重の構成を探すプログラム
mlcf4.path.R	完全多重の構築手順を求めるプログラム
parallel.R	マルチスレッド処理用
progress.R	進行状況の表示用
lnglatdist.R	緯度経度から距離を求める
defaultargs.R	デバッグ用に初期値を格納しておく
mlcf.html	実行結果をgooglemapに表示する際に使用するhtmlのテンプレート
data/hibiyakoen.txt	サンプルデータ（日比谷公園近辺のポータル情報）

mlcf.Rで完全多重を探して、mlcf4.path.Rで見つかった完全多重の構築手順を求めます。

Windows版RはUTF-8をうまく扱ってくれない

Windows版RがUTF-8をうまく扱ってくれないみたいなので、vagrantでlinux仮想環境を使いました。
一部使えない文字がありますが、ShiftJISを使うようにすれば（プログラムも修正すれば）動くかもしれません。（いいかげんですみません）

完全多重を探す

次のコマンドを実行します。

Rscript mlcf.R -4 -cpu4 -clear -f250 -n16 data/hibiyakoen

引数	説明
-4	完全4重を求めます完全5重を求めるときは-5を指定 -3や-6は指定できません
-cpu4	マルチスレッド処理のスレッド数(4)を指定
-clear	前回の実行結果（途中経過）を読み込まない
-f250	頂点間の距離を250メートル以下に限定
-n16	マルチスレッドを16処理単位ごとに実行
data/hibiyakoen	ポータル情報のファイル、拡張子は .txt 固定

実行結果のファイルを２個

ファイル	説明
data/hibiyakoen.mlcf4.map.html	見つかった多重フィールドをgooglemapに表示します
data/hibiyakoen.mlcf4.rdat	見つかった完全多重の構成情報、Rの変数をsave関数で書き出したもの

data/hibiyakoen.mlcf4.map.html をブラウザで開くと、見つかった完全多重の一覧（頂点のポータル名）が表示されます。クリックすると、そのフィールドの位置をgooglemapに表示します。

※ googlemap表示のhtmlは、自分のPCの解像度の合わせて要素サイズを固定で指定しています。見栄えが良くない場合は170行目あたりを編集してください。

完全多重を作る手順を求める

以下のコマンドを実行します。

Rscript mlcf4.path.R data/hibiyakoen 1

引数	説明
data/hibiyakoen	ポータル情報のファイル、拡張子は .txt 固定
N	見つかった完全多重のうちN個目の多重を作成するための情報を求めます

実行結果のファイル１個

ファイル	説明
data/hibiyakoen.N.mlcf4.path.txt	指定された完全多重を作成するための手順（ポータルをどの順番で巡るか、どうリンクを張るか）と必要な鍵数の情報 Nはmlcf4.path.Rの引数で指定したもの

data/hibiyakoen.N.mlcf4.path.txt をテキストエディタで開くと、以下を参照できます。

行番号	説明
1~16	使用するポータルの一覧
18~23	経路パターン（３個の頂点をどの順番で辿るか）括弧内の数字は移動距離です
25	IITCのドローデータ
28~128	経路パターンごとのポータルを辿る順番移動距離を短くするように計算しています
132~208	経路パターン1のリンク手順・必要鍵数・経路情報のIITCドローデータ
以下同様	経路パターン2~6のリンク手順・必要鍵数・経路情報のIITCドローデータ

リンク手順の内容

リンク手順は以下のような内容です。

21番この植樹の名前は？    41
     波 
     思い出ベンチ 
モミジバスズカケノキ  34
     波 
     21番この植樹の名前は？

「21番この植樹の名前は？」から「波」「思い出ベンチ」の順番でリンクを張ります。
次に「モミジバスズカケノキ」から「波」「21番この植樹の名前は？」の順番でリンクを張ります。

ポータルを辿る順番、リンクを張る順番が極めて重要です。この順番でリンクを張っていかないと完全多重が作れません。

ポータル名の後ろの数字は直前のポータルとの直線距離です。
「21番この植樹の名前は？」と「モミジバスズカケノキ」の直線距離は34メートルです。

経路パターンや経路パターンごとのポータルを辿る順番に <NG> が表示されていることがあります。
これについては後述します。

実行！

作戦準備は完了、あとは鍵を集めて多重を作るだけ。

処理内容の説明

ポータルの一覧とその座標を取得する方法

この章の内容はingressを知らない人には全くわからない内容ですね、すみません。

IITCプラグインで「Ingress Dual Map Exporter」を使います
https://gist.github.com/OllieTerrance/8547503

IITCで完全多重を調べたい範囲を表示して「IDM Export」で画面に表示されているポータルの一覧（ポータル名・緯度経度）をCSV形式で取得できます。
これをテキストエディタに張り付けてUTF-8で保存。
以降の処理の都合上、拡張子を .txt にする必要があります。
（サンプルデータの hibiyakoen.txt もこの方法で作成しています）

きわめて簡単なIITCの説明

ingressのポータル情報などをわかりやすく画面表示してくれるプログラムで、greasemonkey上で動きます。https://iitc.me/desktop/ を参照してください。

完全多重が構成できるかを調べる方法

メイン処理（mlcf.R）の解説です。

任意の３点の内側にあるポータルの数を調べる

完全多重フィールドの内側ポータル数は一定です。
内側ポータル数がこれよりも多くても少ないと完全多重を構成できないので、最初に任意の３点を選び、その内側のポータル数を数えることで、完全多重にならないパターンを除外します。

多重度	内側ポータル数
2	1
3	4
4	13
5	40

※ プログラムでは均一多重の構成も求められるよう、余分ポータル数を指定できるようにしています。

再帰的に完全多重の構成を確認

内側ポータル数が上記の数になっていれば、内側ポータルを１個選びます。
３頂点と選んだポータルで３個のフィールドを構成して、その３個のフィールド全てが（多重度が１少ない）完全多重になっているかを調べます。

計算量を削減するための工夫

計算量が多くなるので計算量を削減するための工夫をしています。

ポータルAを選ぶ
Aから2キロ以内のポータルを選ぶ
2の中からポータルBを選ぶ
2の中からポータルBから2キロ以内のポータルを選ぶ（残ったのはポータルABの両方から2キロ以内）
4のポータルについて線分ABとの位置関係を求めておく
4の中からポータルCを選ぶ
4のポータルについて線分AC・線分BCとの位置関係を求める
5と7を使って、4の中からABCの内部にあるポータルを選ぶ
8で選んだポータルの数が必要な内側ポータル数と同じでなければ、ABCでは完全多重を構成できない
8で選んだポータルで完全多重を構成できるかを詳細に調べる

8の「ABCの内部にあるポータルを選ぶ」は次の条件の論理和を求めています
* 線分ABに対してCと同じ側にある
* 線分ACに対してBと同じ側にある
* 線分BCに対してAと同じ側にある

計算量を減らして高速化しようとしたら、その判定処理や抽出処理に時間がかかって、計算量は減ったけど処理時間は増えた、などという試行錯誤を繰り返してます。高速化は難しいですね。

ポータル間の距離で対象を限定しています。この距離を長くすればもっと多くの完全多重の構成を見つけることができるのですが、計算量が増えて処理に時間がかかること、構築するときの移動距離が長くなること、見つかった構成は細長いものになって見栄えがよくないことから、距離の制限をいれています。

マルチスレッドで途中経過を画面表示

計算量が多いのでマルチスレッドで処理させています。
マルチスレッドで処理実行させると途中経過がわかりません。

doSNOWを使うと進行状況を表示できるようですが、それは使わずに自前でやってみました。

# 普通にマルチスレッド
foreach (i in 1:10000) %dopar% {
  # process here
}

# 進行状況を表示しながらマルチスレッド
c1.proc <- 0
chunk <- 100
while (c1.proc < 10000) {
  // 進行状況をここで表示
  c1.chunk <- (c1.proc + 1):min(c1.proc + chunk, 10000)
  foreach (c1 = c1.chunk) %dopar% {
    # process here
  }
}

doSNOWでのprogress barの表示する方法

https://stackoverflow.com/questions/5423760/how-do-you-create-a-progress-bar-when-using-the-foreach-function-in-r を参照してください。

緯度経度を使って２地点間の距離を求める

http://yamadarake.jp/trdi/report000001.html を使わせていただきました。

lnglat.dist <- function(x1, y1, x2, y2, a = 6378137.0, b = 6356752.314140) { 

  deg2rad <- function(x){ 
    x * pi / 180 
  } 

  dy <- deg2rad(y1 - y2) 
  dx <- deg2rad(x1 - x2) 
  my <- deg2rad((y1 + y2) / 2) 
  e2 <- (a^2 - b^2) / a^2 
  Mnum <- a * (1 - e2) 
  W <- sqrt(1 - e2 * sin(my)^2) 
  M <- Mnum / W^3 
  N <- a / W 
  d <- sqrt((dy * M)^2 + (dx * N * cos(my))^2) 
  return(d) 
}

portalsは2列（＝経度・緯度、列名は"lng"と"lat"）N行の行列で、次にように使っています。

# 距離を求める

lnglat.dist(portals[1,"lng"], portals[1,"lat"], portals[2,"lng"], portals[2,"lat"])

2地点間だけでなく複数地点との距離も求められます。

# x2・y2にベクトル（長さが同じ）を渡せば
# x1・y1から複数地点地点への距離をまとめて求められます

lnglat.dist(portals[1,"lng"], portals[1,"lat"], portals[,"lng"], portals[,"lat"])

ある地点が２地点を結ぶ線分のどちらにあるかを判定する

３地点（フィールド）の内側にあるポータルを探すために、線分と地点の位置関係を評価しています。

portalsは2列（＝経度・緯度、列名は"lng"と"lat"）N行の行列、
sele は判定対象かどうかのbooleanのベクトル（長さはportalsの行数と同じ）です。

同じ計算結果ですが、少しづつ高速化しています。

# 地道に計算する方法
  vertex1 <- portals[v1, ]
  vertex2 <- portals[v2, ]
  vx1 <- vertex2["lat"] - vertex1["lat"]
  vy1 <- vertex2["lng"] - vertex1["lng"]
  vx2 <- portals[sele, "lat"] - vertex1["lat"]
  vy2 <- portals[sele, "lng"] - vertex1["lng"]
  dirs <- sign(vx1 * vy2 - vy1 * vx2)

# 変数への代入をやめて少し高速化
  vertex1 <- portals[v1, ]
  vertex2 <- portals[v2, ]
  dirs <- 
    sign((vertex2["lat"] - vertex1["lat"]) * (portals[sele, "lng"] - vertex1["lng"]) -
         (vertex2["lng"] - vertex1["lng"]) * (portals[sele, "lat"] - vertex1["lat"]))

# 外積計算を使ってさらに高速化
  vertex1 <- portals[v1, ]
  vertex2 <- portals[v2, ]
  dirs <- 
    sign( 
      matrix(c((portals[sele, "lng"] - vertex1["lng"]), (portals[sele, "lat"] - vertex1["lat"])), byrow=F, nrow=sum(sele))
        %*% 
      matrix(c((vertex2["lat"] - vertex1["lat"]), (vertex1["lng"] - vertex2["lng"])), byrow=T, nrow=2)

# 複数の線分との位置関係をまとめて求めることでさらに高速化
  vertex1 <- portals[v1, ]
  dirs <- 
    sign( 
      matrix(c((portals[sele, "lng"] - vertex1["lng"]), (portals[sele, "lat"] - vertex1["lat"])), byrow=F, nrow=sum(sele))
        %*% 
      matrix(c((portals[,"lat"] - vertex1["lat"]), (vertex1["lng"] - portals[,"lng"])), byrow=T, nrow=2) )

外積を使って複数の線分に対しての位置関係をまとめて求めることで高速化を図りました。
v1を1とした場合、dirs[3,2] はポータル1とポータル2を結ぶ線分に対しての、ポータル3の位置を示しています。

完全多重の構築方法を考える方法

mlcf4.path.R の解説です。

３個の頂点(ABC)をどの順番で辿るかで６個の経路パターンがあります
- ABC
- ACB
- BAC
- BCA
- CAB
- CBA

各経路パターンごとに ①どの順番でポータルを辿るか 、 ②どこからどこにリンクを張るか と ③必要な鍵数 を計算しています。

どの順番でポータルを辿るか

全てのポータルを辿る経路のうち、できるだけ距離の短いものを探しています。

クラスタ分析してポータルをグループ化します。
クラスタ化の結果がいまいち（１つのクラスタにポータルが集中している、または、その逆で分散していてクラスタ数が多くなった）の場合は、グループ化の基準（クラスタ間距離の閾値）を変更して、グループ化をやり直します。

次に各クラスタの重心間の距離を使って、どういう順番でクラスタを回るのが効率がいいか（距離が短いか）を求めます。

最後にクラスタ内のポータルをどういう順番で回るのが効率が効率がいいかを求めます、

これで求められるのは短距離ではないですが、実際には直線で移動できるわけでもないので、そこそこ短い移動距離のパターンが見つかればＯＫとしています。

どこからどこにリンクを張るか

どのポータル間にリンクを張るかは決まっているので、辿る順番に従って、あとで回るポータルから、先に回ったポータルにリンクを張るようにします。

このままだと、フィールドの内部のポータルからはリンクを張れないという制約のために、リンクを張れないケースが出てきます。
この制約を回避して最適な手順を求めるようにするのはあきらめて、手作業でポータルを辿る順番を指定するようにしました（手抜きです）。

出力したファイル（.mlcf4.path.txt）を編集して、制約にひっかかるリンク（ <NG> が表示された行）を上に（早く行くように）移動してファイルを保存して、コマンドを実行すると（出力ファイルが存在する場合はそれを読み込んで）、リンク手順が制約にひっからないかをチェックして結果をあらためてファイルに書き出します。

必要な鍵数

どこからどこにリンクを張るかを計算した後で、リンク先の数を集計しています。

5重の経路の求め方

経路は４重に対して求めています。
完全５重の経路を求める場合は、５重に含まれる３つの４重の経路を求めて、それを人手でつなぐようにしてください。

均一５重の構成・手順を求める方法

均一５重パターンを求める方法を説明します

以下のコマンドを実行します。
Rscript mlcf.R -588888 -cpu4 -clear -n16 -f500 data/hibiyakoen

引数	説明
-588888	均一5重を求めます、余分ポータル数は8まで許容
-f500	頂点間の距離を500メートル以下に限定

完全5重だけを探すなら -500000 を指定します。
ただし、完全5重はほどんど見つからないので、長時間実行して多重が見つからないという結果になる可能性が高いです。
-533333 などにしておくと、時々均一多重が見つかってgooglemapに表示されていくので、実行していることが分かって安心です。

実行結果のファイル２個

ファイル	説明
data/hibiyakoen.mlcf5.map.html	見つかった多重フィールドをgooglemapに表示します
data/hibiyakoen.mlcf5.rdat	見つかった完全多重の構成情報、Rの変数をsave関数で書き出したもの

多重の作り方を求める

以下のコマンドを実行します（最後の引数を1,2,3にして3回コマンドを実行）

Rscript mlcf4.path.R data/hibiyakoen 1 1
Rscript mlcf4.path.R data/hibiyakoen 1 2
Rscript mlcf4.path.R data/hibiyakoen 1 3

引数	説明
data/hibiyakoen	ポータル情報のファイル
N	見つかった均一５重のうちN個目の多重を作成するための情報を求めます
M	均一５重に含まれる３つの均一多重のうちM個目の多重を作成するための情報を求めます 1,2,3を指定して３回実行します

実行結果のファイルを3個

ファイル	説明
data/hibiyakoen.N.1.mlcf5.proc.txt data/hibiyakoen.N.2.mlcf5.proc.txt data/hibiyakoen.N.3.mlcf5.proc.txt	指定された多重を作成するための手順と必要な鍵数の情報です

5重の手順と鍵数は、この3つの4重の手順を見ながら自分で考えてください。（手抜きです、すみません）
それでもゼロから考えるよりは楽なはずなので。

まとめ

ingressで多重作ってみたい人、よかったらお使いください
Rで位置情報扱ってみたい人、参考になるかならないかわかりませんが、ご自由にお使いください
core i5の4スレッドPCで、ちょっと広めのエリアで動かすと2～3日CPU100%で張り付きます、もうちょっと早くしたい。
これを使ってそこそこ広い範囲（例えば都心部全体）を複数に小分けして多重を探しています。それについてはまた書くかもしれません。

自分について

Rは初心者レベルです（パンダ来ました）
このプログラムは1年以上前に作って隠し持っていました
ingressは引退してます
QiitaもAdvent Calendarも初めてです

Rでingress - 完全多重を探すプログラムをRで作ってみた

前書き