はじめに
「バドミントン戦略シミュレーションとの衝動遊戯：第２ゲーム」
moratoriamuo.hatenablog.com
で, 高速化を図ってRcppによる書き換えを試みたが上手くいかなかったことを書いたところ, ill-identifiedさんからコメントをいただき, コメントにしたがってRcppの書き換えたところ, 上手くいった. 以下ではRcppの書き換えと, それによる速度比較について書く.
Rcppによる書き換え
書き換えたコードが以下のものである.
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix frallycpp(Function f1, Function f0, int ngame, double reffect){
NumericMatrix mresult(6,0);
NumericVector vresult(6);
for(int game = 0; game < ngame; ++game){
NumericVector vcrally(0);
int p1point = 0;
int p0point = 0;
int player = 1;
while(p1point < 21 && p0point < 21){
int judge = 1;
int shot = 0;
int crally = 0;
while(judge==1){
double effect = (shot==0)? 0: reffect;
shot = (player==1)? as<int>(f1(shot)): as<int>(f0(shot));
double prob = (shot==0)? 0.9-effect: 0.7-effect;
judge = (R::runif(0, 1) < prob)? 1: 0;
player = (player==1)? 0: 1;
crally = crally + 1;
}
vcrally.insert(vcrally.end(),crally);
if(player==1){
p1point = p1point + 1;
}else{p0point = p0point + 1;}
}
vresult[0] = mean(vcrally);
vresult[1] = var(vcrally);
vresult[2] = max(vcrally);
vresult[3] = p1point;
vresult[4] = p0point;
vresult[5] = (p1point>p0point)? 1: 0;
mresult = cbind(mresult,vresult);
}
mresult = transpose(mresult);
return mresult;
}
前回詰まっていた2点のうち、一点は型の間違いによるもの、もう一点はRcppのrunif()とRのrunif()で返り値の型が違うことによるものだった. これらを修正し, 対戦する戦略関数とシミュレーション回数とrisky shotの効果を引数にとって、試合結果の行列を返す関数を作成した.
Rの関数についても同様に書き換えた.
frallyR <- function(fstr1,fstr0,ngame,reffect){
#### vresultをrbindするmatrixの作成 ####
mresult <- NULL
vresult <- NULL
#シミュレーション
for(i in 1:ngame){
#### game result 変数の初期値 ####
vcrally <- c() #rally count vector
p1point <- 0 #player1のポイント
p0point <- 0 #player0のポイント
player <- 1 #first serverはplayer1
#### rallyの定義とgameの定義 ####
#rallyをjudge==0になるまでのwhileループとする.
#gameをどちらかのplayerのpointが21になるまでのwhileループとする.
while(p1point < 21 & p0point < 21){
#ラリー前に初期化する
judge <- 1 #shotの成功失敗の判定(1=>成功,0=>失敗)
shot <- 0 #ラリー最初のサーブは成功率へのeffectは0なので相手のshotはsafetyとみなす.
crally <- 0 #rally count
while(judge==1){
#手番のplayerの戦略関数をfstrに代入
fstr <- ifelse(player==1,fstr1,fstr0)
#相手のshotからeffectを決定
effect <- ifelse(shot==0, 0, reffect)
#自分のshotの決定
shot <- fstr(shot)
prob <- ifelse(shot==0, 0.9-effect, 0.7-effect)
#shotの成功判定
judge <- ifelse(runif(1)<prob, 1, 0)
#手番playerの交代
player <- ifelse(player==1, 0, 1)
#rally count追加
crally <- crally + 1
}
#ラリーカウントベクトルの更新
vcrally <- c(vcrally, crally)
#ポイントの更新
if(player==1){
p1point <- p1point + 1
}else{p0point <- p0point + 1}
}
#試合結果(ラリー数平均,ラリー数分散,最大ラリー数p1point,p0point,gamewinplayer)
vresult <- c(mean(vcrally),var(vcrally),max(vcrally),p1point,p0point,ifelse(p1point>p0point,1,0))
mresult <- rbind(mresult,vresult)
}
return(mresult)
}
さらに, 速度比較の対象として並列計算するケースを用意するために, シミュレーションを繰り返さずに結果を返す関数も作成しておく.
fonerallyR <- function(fstr1,fstr0,reffect){
vresult <- NULL
#### game result 変数の初期値 ####
vcrally <- c() #rally count vector
p1point <- 0 #player1のポイント
p0point <- 0 #player0のポイント
player <- 1 #first serverはplayer1
#### rallyの定義とgameの定義 ####
#rallyをjudge==0になるまでのwhileループとする.
#gameをどちらかのplayerのpointが21になるまでのwhileループとする.
while(p1point < 21 & p0point < 21){
#ラリー前に初期化する
judge <- 1 #shotの成功失敗の判定(1=>成功,0=>失敗)
shot <- 0 #ラリー最初のサーブは成功率へのeffectは0なので相手のshotはsafetyとみなす.
crally <- 0 #rally count
while(judge==1){
#手番のplayerの戦略関数をfstrに代入
fstr <- ifelse(player==1,fstr1,fstr0)
#相手のshotからeffectを決定
effect <- ifelse(shot==0, 0, reffect)
#自分のshotの決定
shot <- fstr(shot)
prob <- ifelse(shot==0, 0.9-effect, 0.7-effect)
#shotの成功判定
judge <- ifelse(runif(1)<prob, 1, 0)
#手番playerの交代
player <- ifelse(player==1, 0, 1)
#rally count追加
crally <- crally + 1
}
#ラリーカウントベクトルの更新
vcrally <- c(vcrally, crally)
#ポイントの更新
if(player==1){
p1point <- p1point + 1
}else{p0point <- p0point + 1}
}
#試合結果(ラリー数平均,ラリー数分散,最大ラリー数p1point,p0point,gamewinplayer)
vresult <- c(mean(vcrally),var(vcrally),max(vcrally),p1point,p0point,ifelse(p1point>p0point,1,0))
return(vresult)
}
速度比較
必要なライブラリーの読み込みと戦略関数と上で定義した3つ関数の読み込みを行う.
library(Rcpp)
library(foreach)
library(doSNOW)
## Loading required package: iterators
## Loading required package: snow
source("strategy.R")
sourceCpp("frally.cpp")
source("frally.R")
source("fonerally.R")
単純にRで計算するケース, Rcppを使用して高速化するケース, foreachを使用して並列計算するケースの３つのケースで比較する. 設定としては, プレイヤー1がrisky戦略, プレイヤー0がsafety戦略をとり, シミュレーション回数は10000回, risky shotの効果の大きさを0.284とする.
単純にRで計算するケース
system.time(
mr <- frallyR(fstr_allrisky,fstr_allsafety,10000,0.284)
)
## user system elapsed
## 16.03 0.11 16.22
Rcppを使用して高速化するケース
system.time(
mcpp <- frallycpp(fstr_allrisky,fstr_allsafety,10000,0.284)
)
## user system elapsed
## 23.60 0.00 23.64
foreachを使用して並列計算するケース
cl <- makeCluster(4, type = "SOCK")
registerDoSNOW(cl)
system.time(
mrp <- foreach(i=1:10000, .combine = "rbind") %dopar% fonerallyR(fstr_allrisky,fstr_allsafety,0.284)
)
## user system elapsed
## 8.21 0.47 13.58
stopCluster(cl)
結果の確認
速度比較はできたが, 一応, 中で同様のことができているか確認する. 扱いやすいようにデータフレームの形式に変換する. mr以外の2つの行列についても同様.
dfr <- as.data.frame(mr)
colnames(dfr) <- c("mean","var","max","p1point","p2point","winp")
rownames(dfr) <- NULL
head(dfr)
## mean var max p1point p2point winp
## 1 3.175000 3.276282 9 19 21 0
## 2 3.973684 14.999289 19 17 21 0
## 3 4.121212 8.047348 13 21 12 1
## 4 3.684211 7.411095 11 21 17 1
## 5 3.390244 8.343902 12 20 21 0
## 6 2.825000 3.276282 9 19 21 0
dim(dfr)
## [1] 10000 6
colMeans(dfr[,c("winp","mean","var","p1point","p2point")])
## winp mean var p1point p2point
## 0.495700 3.328941 5.794026 18.440500 18.600500
head(dfcpp)
## mean var max p1point p2point winp
## 1 3.171429 11.969748 21 14 21 0
## 2 3.000000 3.600000 8 20 21 0
## 3 3.625000 5.596774 9 21 11 1
## 4 3.400000 4.041026 10 21 19 1
## 5 3.051282 5.418354 11 18 21 0
## 6 3.317073 5.671951 9 21 20 1
dim(dfcpp)
## [1] 10000 6
colMeans(dfcpp[,c("winp","mean","var","p1point","p2point")])
## winp mean var p1point p2point
## 0.495600 3.331092 5.813098 18.517200 18.557800
head(dfrp)
## mean var max p1point p2point winp
## 1 3.390244 4.693902 8 20 21 0
## 2 3.111111 3.473016 9 21 15 1
## 3 3.151515 5.570076 12 21 12 1
## 4 4.027778 14.027778 16 15 21 0
## 5 2.536585 2.004878 6 21 20 1
## 6 3.717949 10.628880 15 18 21 0
dim(dfrp)
## [1] 10000 6
colMeans(dfrp[,c("winp","mean","var","p1point","p2point")])
## winp mean var p1point p2point
## 0.488000 3.328949 5.789926 18.462100 18.596600
おそらく出力に問題はないと思われる.
おわりに
昨年学んだ並列計算で高速化できたのは良かったけど, Rcppで書き換えたのに、むしろ遅いってどういうことなの……. コードの書き方がおそらく非効率なのだろう. コードの最適化やC++について今後少しずつ学んでみようと思う.
<追記: 2019/01/02>
C++側で, Rで定義した関数を呼び出すところがどうしても遅くなるようで, 高速化をしたいのなら, それは避けるべきみたい. C++で定義し直してやるとよい.
加えて, cbindはサイズが大きくなると遅くなるので, mresultのmatrixを先に用意して代入していく方がおそらく速くなる.
Twitterの有識者(鍵アカウントなので名前の公開は避けます.)の方々のご指摘で気付きました. ありがとうございます.
<追記おわり>
Rcppとforeachという持ち札が増えたのは良かった. Rcppのコードで詰まっていたところについてのコメントは大変助かりました. ill-identifiedさん, ありがとうございました.
