Rでシミュレーション制御

Rでシミュレーション制御！（工事中）

R言語とそのパッケージrgl （OpenGLをRで使うための関数群）でシミュレーションと可視化をマウス制御するプログラムを書いてみました。自由に改変・改良して遊んで下さい。（説明が足りてないと思うので、詳しい解説を追加する予定です。）もしミスがあったら教えて頂けると助かります。windowsで文字化けしないようにプログラムファイルの文字コードをShift-JISにしました。githubに置いた方はutf-8のまま。

３次元プロット
２次元プロット
2Dと3Dを統合+簡易メニュー
使用例１：１次元形質空間における進化的分岐
使用例２：２次元空間における捕食者と被食者の個体数動態
使用例３：２次元反応拡散モデル
windowsとmacでは動作が多少不安定かも。①「plotR」のボタンは押さない、②クリック連打をしない、③メニューによる描画の変更は一時停止してから行う、でエラー停止を減らせると思います。（僕のmacではrglウインドウとRデフォルトの描画ウインドウを両方開くと、次にrglを使う時にRが落ちます。）
チラつきを抑えるための仕組みについては「３次元プロット」のプログラムのメインループ（最後の部分）を見て下さい。

Linux,Windows,Macでの使用を想定して書いています。（主にLinux: Ubuntu 18.04, MacOs: Mojave, Windows: win7とwin10, R-3.4.4〜4.0, rgl-0.100.54で動作確認）
３次元プロットはrglの描画関数をそのまま使えばOK
２次元プロットは試行錯誤の末、Rで描いた２次元プロットの画像をrglのx-y平面に張り込む方法にしました。

ソースコードはgithubに置いたファイルの方が色分けされてて見やすいかも
実行するためにはまずrglをインストールして下さい
```
install.packages("rgl")
```

参考になるサイト

rglについて：

Heavy Watal: rglの解説

Rの関数の定義や振る舞い：

３次元プロットの場合

マウス左ボタンで回転（デフォルトのまま）
マウス右ボタンで一時停止（ウインドウの上端のバーの下辺りでクリックすると終了）
rglのウインドウを閉じるコマンドはrgl.close()
実行方法：プログラムファイル"rgl_demo3d.R"を保存して、そのディレクトリ（フォルダ）でRを起動して以下のコマンドを入力すれば計算が実行されグラフが表示されます（日本語のコメントが文字化けする場合は下に貼った中身をコピペして下さい）。windowsの場合はRを起動して「Rソースコードの読み込み」のような選択⇒rgl_demo3d.Rでも実行できるはず。
```
source("rgl_demo3d.R")
```

プログラムの中身は以下の通り。

library(rgl)  ##rglをロード

win_pos_x=100; ##ウィンドウの生成位置のｘ座標
win_pos_y=200; ##ウィンドウの生成位置のｘ座標
win_size=600; ##ウィンドウサイズ（ピクセル）

flag_pause=0; ##０：一時停止しない, １：一時停止
flag_halt=0; ##１で終了

dt=0.07; ##時間刻み幅
step=0; ##タイムステップ

ncol=50; ##色パレットの色数
col_pal=topo.colors(ncol); ##色パレット
zmax=2; ##色パレットの上限色に対応する値
zmin=-2; ##色パレットの下限色に対応する値

##ｚの値から色を割り当てる関数
calc_col <- function(z,zmin,zmax,ncol,col_pal){
    cv=as.integer(ncol*(z-zmin)/(zmax-zmin));
    cv[which(cv<1)]=1;
    cv[which(cv>ncol)]=ncol;
    return(col_pal[cv]);
}

##マウスクリックで一時停止と終了を制御する関数：
##マウスクリックのｙ座標がウインドウの上部から２０％の範囲にある場合は
##終了(flag_halt=1)
##それ以外の場合は一時停止(flag_pauseを０⇒１ or １⇒０)
mouse_pause <- function(mpos_x,mpos_y){
    if(mpos_y < win_size*0.2){
        flag_halt <<- flag_halt*0+1;
        cat("ended\n");
    }
    else{
        flag_pause <<- (flag_pause+1)%%2;
        cat("pause:",flag_pause,"\n");
    }
}

##rglのウインドウを生成
open3d(windowRect=c(win_pos_x, win_pos_y, win_pos_x+win_size, win_pos_y+win_size));

##マウス右ボタンのクリックでmouse_pose関数が実行されるようにセット
rgl.setMouseCallbacks(button=2, begin=mouse_pause,dev=rgl.cur(),subscene = currentSubscene3d(rgl.cur()));

##枠や視点の向きをセット
decorate3d(xlim=c(-1,1),ylim=c(-1,1),zlim=c(-2,2),box=F,axes=T,aspect=TRUE);
    par3d(userMatrix=(rotationMatrix(-0.3*pi,1,0,0) %*% rotationMatrix(0.1*pi,0,0,1)))


x=seq(-1,1,,64); ##-１〜１までのベクトルを生成
y=x;
xx=matrix(x, nrow=64, ncol=64, byrow=T); ##2次元格子のx座標
yy=matrix(x, nrow=64, ncol=64, byrow=F); ##2次元格子のy座標

##メインループ
while(1){
    if(flag_pause==0){##一時停止でないなら以下を実行する
        step=step+1;
        t=dt*step;
	##ｚの更新（ここにシミュレーション１ステップの実行関数をセットすればいい）	       
        z=sin(2*sqrt(xx^2+yy^2)+t*0.3)+cos(xx+1.3*yy+t);

        ##ｚの値を曲面表示して、そのオブジェクトidをb1に格納
        b1=terrain3d(x,y,z,color=calc_col(z,zmin,zmax,ncol,col_pal));

        ##最初以外は古い方のオブジェクトを削除   
        if(step>1)rgl.pop(id=b0);
        b0=b1; ##b1に格納していたオブジェクトidをb0に格納       
    }
    
    Sys.sleep(0.001); ##0.001秒待機
    if(flag_halt==1)return(0); ## flag_haltが１ならば終了する     
}

２次元プロットの場合

マウス左ボタンでrglウインドウに描画したものをR標準のウインドウに再描画（動画ではRのウインドウが小さいですが、サイズは自由に選べます）
マウス右ボタンで一時停止（ウインドウの上端のバーの下辺りでクリックすると終了）
通常のRの２次描画関数は一通り使えるはず
rglのウインドウを閉じるコマンドはrgl.close()
実行方法：プログラムファイル"rgl_demo2d.R"（９５行くらい）を保存して、
```
source("rgl_demo2d.R")
```

プログラムの中身は以下の通り。

library(rgl)  ##rglをロード

win_pos_x=100; ##ウィンドウの生成位置のｘ座標
win_pos_y=200; ##ウィンドウの生成位置のｘ座標
win_size=600; ##ウィンドウサイズ（ピクセル）

flag_pause=0; ##０：一時停止しない, １：一時停止
flag_halt=0; ##１で終了

dt=0.02; ##時間刻み幅
step=0; ##タイムステップ
zmax=2;##等高線の上限値
zmin=-2;##等高線の下限値

pwin_size=600;##2次元プロットの解像度（ピクセル）

##マウスクリックで一時停止と終了を制御する関数：
##マウスクリックのｙ座標がウインドウの上部から２０％の範囲にある場合は
##終了(flag_halt=1)
##それ以外の場合は一時停止(flag_pauseを０⇒１ or １⇒０)
mouse_pause <- function(mpos_x,mpos_y){
    if(mpos_y < win_size*0.2){
        flag_halt <<- flag_halt*0+1;
        cat("ended\n");
    }
    else{
        flag_pause <<- (flag_pause+1)%%2;
        cat("pause:",flag_pause,"\n");
    }
}

##Rデフォルトのウインドウにプロットする関数
mouse_pfunc <- function(mpos_x,mpos_y){
    pfunc(x,x,z);
    cat("plotted in R window\n");
}

##rglのウインドウを生成
open3d(windowRect=c(win_pos_x, win_pos_y, win_pos_x+win_size, win_pos_y+win_size));

##マウス右ボタンのクリックでmouse_pose関数が実行されるようにセット
rgl.setMouseCallbacks(button=2, begin=mouse_pause,dev=rgl.cur(),subscene = currentSubscene3d(rgl.cur()));

##マウス左ボタンのクリックでmouse_pfunc関数が実行されるようにセット
rgl.setMouseCallbacks(button=1, begin=mouse_pfunc,dev=rgl.cur(),subscene = currentSubscene3d(rgl.cur()));

##枠なしの設定
decorate3d(xlim=c(-1,1),ylim=c(-1,1),zlim=c(-2,2),box=FALSE,axes=FALSE,aspect=TRUE,xlab=NULL,ylab=NULL,zlab=NULL);

##視点をｚ軸に沿う下向き方向にセット
view3d(theta=0,phi=0,fov=0,zoom=0.6);

x=seq(-1,1,,64); ##-１〜１までのベクトルを生成
y=x;
xx=matrix(x, nrow=64, ncol=64, byrow=T); ##2次元格子のx座標
yy=matrix(x, nrow=64, ncol=64, byrow=F); ##2次元格子のy座標

##2次元プロット関数（Rの描画関数は大体使えるはず）
pfunc <-function(x,y,z){
    nlev=17;
    clwd=rep(1.0,nlev);##等高線の太さを１に設定
    ccol=rep("black",nlev);##等高線の色を黒にする
    ccol[8]="red";##8番目だけ赤くする
    clwd[8]=3.0;##8番目だけ太さを３にする
    ##塗り潰し等高線(filled.contour)と等高線(contour)を合わせてプロット
    filled.contour(x,x,z,col=topo.colors(20),xlab="x",ylab="y",cex.lab=2.0,levels=seq(zmin,zmax,,nlev),line=1.5,main=paste0("step: ",step),
                   plot.axes = contour(x,x,z, levels =seq(zmin,zmax,,nlev), drawlabels = TRUE, axes = FALSE, frame.plot = FFALSE, add = TRUE,labcex=1.3,col=ccol,lwd=clwd)
                   )
   
}

##メインループ
while(1){
    if(flag_pause==0){##一時停止でないなら以下を実行する
        step=step+1;
        t=dt*step;
        
        ##ｚの更新（ここにシミュレーション１ステップの実行関数をセットすればいい）
        z=sin(2*sqrt(xx^2+yy^2)+t*0.3)+cos(xx+1.3*yy+t);
        
        ##pfuncでプロットしたグラフをpng画像にしてx-y平面に張り込んで
        ##そのオブジェクトidをb1に格納
        b1=show2d(pfunc(x,x,z),width=pwin_size,height=pwin_size);

        ##最初以外は古い方のオブジェクトを削除   
        if(step > 1)rgl.pop(id=b0);    
        b0=b1;##b1に格納していたオブジェクトidをb0に格納       
        
    }
    
    Sys.sleep(0.001);##0.001秒待機
    if(flag_halt==1)return(0);## flag_haltが１ならば終了する     
}

2Dと3Dを統合（簡易メニュー付き）

マウス制御関数を"conrgl-0.12.R"というファイルにまとめました
マウス右クリックでメニューの各項目を実行できます
マウス左で回転。マウス中央で拡大/縮小
rglのウインドウを閉じるコマンドはrgl.close()
「command」をクリックすると、プロンプトが表示されるので、コマンドが使えます。Enterキーだけ入力するとコマンドモード終了。
自分用のメニュー項目を簡単に追加できます
主にLinuxで動作確認してます。windowsとmacでは動作が多少不安定っぽい。①「plotR」のボタンは押さない、②クリック連打をしない、③メニューによる描画の変更は一時停止してから行う、でエラー停止を減らせると思います。エラーを抑制するために、2D描画でシミュレーションをまわす場合はRのウインドウを閉じるようにしてます(plotR_only_pause=1を0にすれば解除)。
使い方：conrgl-0.12.Rとデモプログラム conrgl_demo2.Rを保存して同じディレクトリ（フォルダ）に置いて、
```
source("conrgl_demo2.R")
```
「conrgl-0.12.Rが見つからない」的なエラーが出る場合はconrgl-0.12.Rを同梱したデモファイルconrgl-0.12_with_demo2.Rを保存して
```
source("conrgl-0.12_with_demo2.R")
```
(windowsの場合はRを起動して「Rソースコードの読み込み」のような選択⇒conrgl-0.12_with_demo2.Rでも実行できるはず)
このデモではメニューに2D描画の等高線を増やすボタンと減らす項目を追加しています

conrgl_demo2.Rの中身は以下の通り。

library(rgl); ##rglをロード
source("conrgl-0.12.R"); ##conrgl-0.12.Rを読む
 
cgl.plot_dim=3; ##最初に表示する描画の次元
cgl.xlim=c(-1,1); ##ｘ軸の表示範囲
cgl.ylim=c(-1,1); ##ｙ軸の表示範囲
cgl.zlim=c(-2,2); ##ｚ軸の表示範囲
cgl.xlab=NULL; ##ｘ軸のラベル
cgl.ylab=NULL; ##ｙ軸のラベル
cgl.zlab=NULL; ##ｚ軸のラベル

flag_show_menu_state=0; ##メニュー各項目の状態の表示/非表示
dt=0.07; ##時間刻み幅
step=0; ##時間ステップ
nlevs=17; ##等高線の数
zmax=2; ##等高線の最大値
zmin=-2; ##等高線の最小値

 ##等高線を増やすボタンに対応する関数
action_contour_inc <- function(mst){
    nlevs <<- nlevs+1;
    if((flag_pause==1)+(flag_halt==1) > 0)try(cgl.plot());
    Sys.sleep(0.1);
}

 ##等高線を減らすボタンに対応する関数
action_contour_dec <- function(mst){
    nlevs <<- max(nlevs-1,1);
    if((flag_pause==1)+(flag_halt==1) > 0)try(cgl.plot());
    Sys.sleep(0.1);
}

 ##等高線を増やすボタン要素を作成
item_contour_inc=list(state=1,n=1,label=c("cont+"),color=c("dark green"));
 ##等高線を減らすボタン要素を作成
item_contour_dec=list(state=1,n=1,label=c("cont-"),color=c("dark green"));

 ##メニューに自分用のボタン要素を追加
menu=c(menu,list(contour_inc=item_contour_inc,
                 contour_dec=item_contour_dec));
action=c(action, c(action_contour_inc,
                   action_contour_dec));

nmenu=length(menu); ##メニューの数

##メニュー順序の入れ替え例
menu_order=seq(nmenu);
##menu_order=c(1,2,3,4,6,5,7,8);
action=action[menu_order]; 
menu=menu[menu_order]; 

##2DプロットとRウインドウでの描画を行う関数(cgl.plot()の中で実行される)
##ここを改変すれば自分用の描画ができる
cgl.plotR <- function(){
    clwd=rep(1.0,20);
    ccol=rep("black",20);
    ccol[8]="red";
    clwd[8]=3.0;
    filled.contour(x,x,z,col=mypal(nlevs-1),xlab="x",ylab="y",cex.lab=2.0,levels=seq(zmin,zmax,,nlevs),line=1.5,main=paste0("step: ",step, "  nlevs:",nlevs),
                   plot.axes = contour(x,x,z, levels =seq(zmin,zmax,,nlevs), drawlabels = TRUE, axes = FALSE, frame.plot = FFALSE, add = TRUE,labcex=1.1,col=ccol,lwd=clwd)
                   )
   
}

##3Dプロットを行う関数(cgl.plot()の中で実行される)
##注意！！この関数はcgl.plotRと異なり、描画関数の戻り値を返すように設定する!
##ここを改変すれば自分用の描画ができる
cgl.plotRgl <- function(){
                cv=as.integer(ncol*(z-zmin)/(zmax-zmin));
                cv[which(cv < 1)]=1;
                cv[which(cv > ncol)]=ncol;
               
                return(terrain3d(x,x,z,color=col_pal[cv]));
}

cgl.init(); ##初期化のための関数

t=0.0; ##時刻を０にする
x=seq(-1,1,,64); ##-１〜１までのベクトルを生成
xx=matrix(x, nrow=64, ncol=64, byrow=T); ##2次元格子のｘ座標
yy=matrix(x, nrow=64, ncol=64, byrow=F); ##2次元格子のｙ座標
z=sin(2*sqrt(xx^2+yy^2)+t*0.3)+cos(xx+1.3*yy+t); ##各位置のｚの値を計算

cgl.plot(); ##初期状態の描画

while(1){
    if(flag_pause==0){##一時停止でないなら以下を実行する
        step=step+1;
        t=dt*step;
        ##ｚの更新（ここにシミュレーション１ステップの実行関数をセットすればいい）	
        z=sin(2*sqrt(xx^2+yy^2)+t*0.3)+cos(xx+1.3*yy+t);
        cgl.plot(); ##描画
        if(step%%10==0)cat("step: ",step,"\n");
 
    }
    
    Sys.sleep(0.001); ##0.001秒待機
    if(flag_halt==1)return(0); ## flag_haltが１ならば終了する    
}

使用例１：１次元空間における進化的分岐

「RでAdaptive Dynamics」の１次元形質空間における進化的分岐のプログラムに適用しました。
実行方法：conrgl-0.12.Rとデモプログラム conrgl_demo_branch1d.Rを保存して同じディレクトリ（フォルダ）に置いて、
```
source("conrgl_demo_branch1d.R")
```
「conrgl-0.12.Rが見つからない」的なエラーが出る場合はconrgl-0.12.Rを同梱したデモファイルconrgl-0.12_with_demo_branch1d.Rを保存して
```
source("conrgl-0.12_with_demo_branch1d.R")
```
(windowsの場合はRを起動して「Rソースコードの読み込み」のような選択⇒conrgl-0.12_with_demo_branch1d.Rでも実行できるはず)

conrgl_demo_branch1d.Rの中身は以下の通り。

library(rgl); ##rglをロード
source("conrgl-0.12.R"); ##conrgl-0.12.Rを読む

K0=10.0; ##環境収容力の最大値
sK=1.0; ##環境収容力の幅
sa=0.5; ##競争効果の幅
##sa=0.7; 
##sa=1.2;
ndiv=64; ##形質軸の刻み幅
edge_die=1e-4; ##最小個体数密度。これ以下は絶滅。
repeat_p=1000; ##変異体を出現させずに個体数動態をまわすステップ数（速く計算するため）
dt=0.2; ##時間刻み幅
mutation_rate=3.0e-4; ##突然変異率

time_end=150000;##終了時刻
t_end=as.integer(time_end/(dt*repeat_p)); ##総タイムステップ数
set.seed(23); ##乱数の種。これをコメントアウトすれば毎回少し違う進化動態になる

xx=seq(-1,1,length=ndiv);  ##とり得る形質値のセット
x=xx;
n=xx*0.0;  ##形質軸上の個体数分布
ancestor=c(ndiv/8); ##初期の表現型を指定。
flag_quick=T;
out_interval=5;
show_interval=1;
wins=500;

amp_n=1000;
zcol_max=log(amp_n*K0); ##等高線の最大値
zcol_min=log(amp_n); ##等高線の最小値
cgl.plot_dim=3; ##最初に表示する描画の次元
cgl.xlim=c(-1,1); ##ｘ軸の表示範囲
cgl.ylim=c(0,time_end); ##ｙ軸の表示範囲
##cgl.zlim=c(0,K0*2); ##ｚ軸の表示範囲
cgl.zlim=c(-1,1); ##ｚ軸の表示範囲
cgl.xlab="Trait x (niche position)"; ##ｘ軸のラベル
cgl.ylab="Time"; ##ｙ軸のラベル
cgl.zlab=NULL; ##ｚ軸のラベル
cgl.phi=-0.1*pi;
cgl.theta=0.0*pi

nlevs=20; ##等高線の数（被食者）
nlevs2=7; ##等高線の数（捕食者）

Rxlab = "Trate x"; ##x軸のラベル
Rylab = "Rime"; ##y軸のラベル


 ##2DプロットとRウインドウでの描画を行う関数
cgl.plotR <- function(){
    zcol=log(n_out+edge_die);

    image(xx,t_out,zcol,useRaster=F,col=mypal(nlevs),xlab="Trait x",ylab="Time",cex.lab=1.3,main=paste0("Time: ",t_out[out_count]));
    
    lines(x=xx,y=time_end*0.5*land_out[,out_count]+1.0,col="cyan",lwd=2);
    lines(x=xx,y=time_end*0.01*log((n+edge_die)/edge_die)+1.0,col="red",lwd=2);
    
}

##3Dプロットを行う関数（描画関数の戻り値を返すように設定する!）
cgl.nb1=2;##オブジェクトの数
cgl.plotRgl <- function(){
    idd=out_count;
    if(idd<5)idd=5;
    
    ##zcol=log(amp_n*(n_out[,1:idd]+1));
    zcol=land_out;
    zcol_min=-1;
    zcol_max=1;
    cv=as.integer(ncol*(zcol-zcol_min)/(zcol_max-zcol_min));
    cv[which(cv < 1)]=1;
    cv[which(cv > ncol)]=ncol;
    cols=col_pal[cv];
    cols[which(n_out[,1:idd]>=edge_die)]="#FF5500";
    bid=terrain3d(xx,t_out[1:idd],(0.02*n_out+land_out)[,1:idd],color=cols);
    material3d(alpha=c(0.5));
    bid=c(bid,terrain3d(xx,t_out,0*land_out,color="blue"));
    material3d(alpha=c(1.0));
        return(bid);
}

 ##環境収容力の分布
K <- function(x){
    return (K0*exp(-x^2/(2.0*sK^2)));
}

 ##競争効果の関数
alpha <- function (x0,x1){
    return (exp(-(x0-x1)^2/(2.0*sa^2)));
}

 ##変異型x1の適応度
fitness <-function(x1,xx,n){
    return(1.0-sum(alpha(xx,x1)*n)/K(x1));
}

 ##突然変異の関数
mutate <- function(xx,n){
    for(i in 1:length(xx)){
        if(rbinom(1,size=1,prob=mutation_rate*n[i]*repeat_p*dt)>0){
            if(rbinom(1,size=1,prob=0.5)>0){
                if((i < ndiv))n[i+1]=n[i+1]+edge_die*2.0;
            }else{
                if((i > 1))n[i-1]=n[i-1]+edge_die*2.0;
            }
        }
    }
    return(n);
}

n[ancestor]=K(xx[ancestor])*0.5;  ##ancestorの位置に０でない個体数をセット
n_out=matrix(ndiv*(t_end/out_interval+1),nrow=ndiv,ncol=(t_end/out_interval+1))*0.0; ##出力用行列（個体数分布の時系列データ）
land_out=n_out*0.0;##適応度地形を格納する行列
##t_out=dt*(0:(t_end/out_interval)); ##出力用の時刻データ
t_out=dt*out_interval*repeat_p*(0:(t_end/out_interval)); ##出力用の時刻データ
time=0.0;  ##時刻
n_out[,1]=n;  ##出力用行列に初期状態を格納
t_exec_start=proc.time();
out_count=1;
##初期設定
cgl.init();
##プロットしておく
cgl.plot();
t=0;
while(1){
    if(t  == t_end+1)flag_halt <<- 1;
    if(flag_halt==1){
                return(0);
}
    if(flag_pause ==0){
        t = t+1;
        n=mutate(xx,n);  ##突然変異   
        mask_n=(n>edge_die);  ##存在する（最小個体数密度を上回る）表現型に目印をつける
        xxb=xx[mask_n];  ##存在する表現型だけにする（計算を速くするため）
        nb=n[mask_n];  ##存在する表現型だけの個体数密度
        dnb=nb*0.0;
        fit=dnb*0.0;

        if(flag_quick==TRUE){## exp関数を呼ぶ回数を減らす小細工
            len_nb=length(nb);
            alphab=matrix(len_nb*len_nb,nrow=len_nb,ncol=len_nb);
            for(i in 1:len_nb){
                for(j in 1:len_nb)alphab[i,j]=alpha(xxb[i],xxb[j]);
            }
            Kb=K(xxb);
        }
        
        for(k in 1:repeat_p){            
            if(flag_quick==TRUE){
                fit=(1.0-colSums(alphab*nb)/Kb); ##適応度の計算
                nb=nb+dt*nb*fit;  ##個体数密度を変化させる
            }
            else{
                for(i in 1:length(fit))fit[i]=fitness(xxb[i],xxb,nb);  ##適応度の計算 
                dnb=nb*fit;  ##増加率の計算
                nb=nb+dt*dnb; ##個体数密度を変化させる
            }
            
            time=time+dt;  ##時刻を更新
        }
        n[mask_n]=nb;  ##全ての表現型の個体数のベクトルに戻す
        n[which(n<=edge_die)]=0.0;  ##edge_die以下の表現型は絶滅
        
                
        if(t%% out_interval ==0){
            n_out[,out_count]=n;  ##出力用行列に格納
            out_count=out_count+1;
            cat("  time:", time,"\n");

            buf=rep(0,ndiv);
            for(i in 1:length(land_out[,1])) buf[i] = fitness(xx[i],xxb,nb);
            land_out[,out_count]=buf;
            
            if(out_count%% show_interval==0){

                cgl.plot();
                
            }
        }

    }
    Sys.sleep(0.001);
}


cat("calculation time:\n");
print(proc.time()-t_exec_start);

使用例２：２次元空間における捕食者と被食者の個体数動態

環境勾配：
- x座標が大きい場所⇒被食者の移動が速い
- y座標が大きい場所⇒捕食者の移動が速い
表示スタイル：
- 捕食者密度は高：黄色、低い：青色
- 被食者密度は３Dでは高さ、２Dでは赤い等高線
一番下の緑の項目を選ぶと捕食者と被食者の表示スタイルを入れ替えることができます。
conrgl-0.12_with_demo_pred_prey4.Rを保存して
```
install.packages("matrixStats")
```
で"matrixStats"をインストールしたら
```
source("conrgl-0.12_with_demo_pred_prey4.R")
```

メイン部分は以下の通り。

library(matrixStats);

 ##捕食者と被食者の描画スタイルを入れ替えるボタン"switch_pp"
item_switch_pp=list(state=1,n=2,label=c("prey:height\npred:color","pred:height\nprey:color"),color=c("darkgreen","darkgreen"));
 ##switch_ppの動作
action_switch_pp <- function(mst){ 
    if((flag_pause==1)+(flag_halt==1) > 0)cgl.plot();
}

##メニューにswitch_ppを追加
menu=c(menu,list(switch_pp=item_switch_pp));
action=c(action, c(action_switch_pp));
nmenu=length(menu); ##メニューの数

nn=128; ##格子点の数
view_interval=5; ##
r=4.0; ##被食者の内的自然増加率
K=20.0; ##被食者の環境収容力
e_rate=2.0; ##捕食効率
M_eat=5.0; ##最大捕食量
d=1; ##捕食者の死亡率
eff=0.8; ##栄養効率
dt=0.04; ##時間刻み幅
delta=0.5; ##格子間の距離

x=seq(-1,1,,nn); ##-１〜１までのベクトルを生成
xx=matrix(x, nrow=nn, ncol=nn, byrow=F); ##2次元格子のx座標
yy=matrix(x, nrow=nn, ncol=nn, byrow=T); ##2次元格子のy座標
##格子状の捕食者密度、被食者密度をセット
pred=xx*0.0;
prey=xx*0.0;
pred[10,10]=10.0;
prey[10,10]=0.1*K;
prey=prey+runif(length(prey))*0.001;

difmask=c(1,seq(nn),nn);##移動の効果の計算に使う（吸収壁）
##difmask=c(nn,seq(nn),1); ##移動の効果の計算に使う（反対側と連結）
##被食者の移動速度（ｘが大きい場所ほど速い）
rate_prey=0.06*(5*(xx+1.1));
##被食者の移動速度（ｙが大きい場所ほど速い）
rate_pred=0.03*(5*(yy+1.1));

amp_prey=1000; ##logスケール表示の調整パラメータ
amp_pred=1000; ##logスケール表示の調整パラメータ

zcol_max=log(amp_pred*K); ##等高線の最大値
zcol_min=log(amp_pred); ##等高線の最小値
cgl.plot_dim=3; ##最初に表示する描画の次元
cgl.xlim=c(-1,1); ##ｘ軸の表示範囲
cgl.ylim=c(-1,1); ##ｙ軸の表示範囲
cgl.zlim=c(log(amp_prey),2*log(amp_prey*K)); ##ｚ軸の表示範囲
cgl.xlab ="x";##ｘ軸のラベル
cgl.ylab ="y";##ｙ軸のラベル
cgl.zlab=NULL; ##ｚ軸のラベル
nlevs=7; ##等高線の数（被食者）
nlevs2=7; ##等高線の数（捕食者）

mypal2=heat.colors; ## current pallete
clwd=rep(1.0,nlevs2);
clwd[nlevs/2:nlevs2]=2.0;
Rxlab ="x (遅い<-- 被食者の移動速度 -->速い)"; ##Rウインドウでのx軸のラベル
Rylab ="y (遅い<-- 捕食者の移動速度 -->速い)"; ##Rウインドウでのy軸のラベル
 ##2DプロットとRウインドウでの描画を行う関数
cgl.plotR <- function(){
    if(menu$switch_pp$state==1){
        z <<- log(amp_prey*(prey+1));zcol <<- log(amp_pred*(pred+1));
    }
    else{
        z <<- log(amp_pred*(pred+1));zcol <<- log(amp_prey*(prey+1));
    }
    filled.contour(x,x,zcol,col=mypal(nlevs-1),xlab=Rxlab,ylab=Rylab, cex.lab=2.0,levels=seq(cgl.zlim[1],0.5*cgl.zlim[2],,nlevs),line=1.5,main=paste0("step: ",step, "  nlevs:",nlevs),
                   plot.axes = contour(x,x,z, levels =seq(zcol_min,zcol_max,,nlevs2), drawlabels = FALSE, axes = FALSE, frame.plot = FFALSE, add = TRUE,labcex=1.1,col=mypal2(nlevs2),lwd=clwd)
                   )   
}

 ##3Dプロットを行う関数（描画関数の戻り値を返すように設定する!）
cgl.plotRgl <- function(){
    if(menu$switch_pp$state==1){
        z <<- log(amp_prey*(prey+1));zcol <<- log(amp_pred*(pred+1));
    }
    else{
        z <<- log(amp_pred*(pred+1));zcol <<- log(amp_prey*(prey+1));
    }
        cv=as.integer(ncol*(zcol-zcol_min)/(zcol_max-zcol_min));
        cv[which(cv < 1)]=1;
        cv[which(cv > ncol)]=ncol;               
        return(terrain3d(x,x,z,color=col_pal[cv]));
}

step=0; ##時間ステップ
t=0.0; ##時刻を０にする
cgl.init(); ##初期化のための関数
cgl.plot(); ##初期状態の描画

while(1){
    if(flag_pause==0){##一時停止でないなら以下を実行する
        for(k in 1:view_interval){
            ##相互作用による変化率
        fresp = e_rate*prey/(1+prey/M_eat);
        dprey = r*prey*(1 - prey/K) - pred*fresp;
        dpred = eff*pred*fresp - d*pred;
            ##移動量
        preydif =  (colDiffs(colDiffs(prey[difmask,]))+rowDiffs(rowDiffs(prey[,difmask])))/(delta*delta);
        preddif =  (colDiffs(colDiffs(pred[difmask,]))+rowDiffs(rowDiffs(pred[,difmask])))/(delta*delta);
            ##個体数の変化
        prey = prey + dt*(dprey + rate_prey*preydif);
        pred = pred + dt*(dpred + rate_pred*preddif);
                
        step=step+1;
        t=dt*step;
        if(step%%100==0)cat("step: ",step,"\n");
        }
       
     
        cgl.plot(); ##描画         
    }
    
    Sys.sleep(0.001); ##0.001秒待機
    if(flag_halt==1)return(0); ## flag_haltが１ならば終了する    
}

使用例３：２次元反応拡散モデル

２成分の反応拡散モデル「Gray-Scottモデル」です（熱帯魚の模様とかを説明できるモデル）。ここでは１つ目の成分を捕食者、２つ目を捕食者と呼びます。

環境勾配：
- x座標が大きい場所⇒被食者の移動が速い
- y座標が大きい場所⇒捕食者の移動が速い
表示スタイル：
- 捕食者密度は高：黄色、低い：青色
- 被食者密度は３Dでは高さ、２Dでは赤い等高線
一番下の緑の項目を選ぶと捕食者と被食者の表示スタイルを入れ替えることができます。
conrgl-0.12_with_demo_rd1.Rを保存して,
```
install.packages("matrixStats")
```
で"matrixStats"をインストールしたら
```
source("conrgl-0.12_with_demo_rd1.R")
```

メイン部分は以下の通り（使用例２と大体同じ）。

library(matrixStats);##行列演算のパッケージをロード

 ##捕食者と被食者の描画スタイルを入れ替えるボタン"switch_pp"
item_switch_pp=list(state=1,n=2,label=c("prey:height\npred:color","pred:height\nprey:color"),color=c("darkgreen","darkgreen"));
 ##switch_ppの動作
action_switch_pp <- function(mst){ 
    if((flag_pause==1)+(flag_halt==1) > 0)cgl.plot();
}

##メニューにswitch_ppを追加
menu=c(menu,list(switch_pp=item_switch_pp));
action=c(action, c(action_switch_pp));
nmenu=length(menu); ##メニューの数

nn=100; ##格子点の数
view_interval=100; ##描画間隔
a = 0.023; ##被食者の世代交代速度
b = 0.055; ##a+bが捕食者の死亡率

dt=0.2; ##時間刻み幅
x=seq(-1,1,,nn); ##-１〜１までのベクトルを生成
yy=matrix(x, nrow=nn, ncol=nn, byrow=T); ##2次元格子のy座標
xx=matrix(x, nrow=nn, ncol=nn, byrow=F); ##2次元格子のx座標

prey=xx*0.0+1.0; ##被食者密度
pred=xx*0.0+0.001; ##捕食者密度
pred[nn/2,nn/3]=0.5; ##ここだけ0.5にしとく

difmask=c(1,seq(nn),nn); ##移動の効果の計算に使う（吸収壁）
##difmask=c(nn,seq(nn),1); ##移動の効果の計算に使う（反対側と連結）
buf=seq(0.3,2.0,,nn); ##0.3から２までのベクトルを生成
 ##被食者の移動速度（ｘが大きい場所ほど速い）
rate_prey=matrix(0.08*buf, nrow=nn, ncol=nn, byrow=F);
##被食者の移動速度（ｙが大きい場所ほど速い）
rate_pred=matrix(0.046*buf, nrow=nn, ncol=nn, byrow=T);

zcol_max=1; ##等高線の最大値
zcol_min=0; ##等高線の最小値
cgl.plot_dim=3; ##最初に表示する描画の次元
cgl.xlim=c(-1,1); ##ｘ軸の表示範囲
cgl.ylim=c(-1,1); ##ｙ軸の表示範囲
cgl.zlim=c(0,2); ##ｚ軸の表示範囲
cgl.xlab = "x"; ##ｘ軸のラベル
cgl.ylab = "y"; ##ｙ軸のラベル
cgl.zlab=NULL; ##ｚ軸のラベル
nlevs=7; ##等高線の数（被食者）
nlevs2=7; ##等高線の数（捕食者）

mypal2=heat.colors; ## ２D描画で被食者を等高線で表示した場合の線色用のパレット
clwd=rep(1.0,nlevs2); ##２D描画で被食者を等高線で表示した場合の線の太さ
clwd[nlevs/2:nlevs2]=2.0; ##上半分を太くしとく
Rxlab ="x (slow<-- Predator migration -->fast)"; ##Rウインドウでのx軸のラベル
Rylab ="y (slow<-- Prey migration -->fast)"; ##Rウインドウでのy軸のラベル
		    
 ##2DプロットとRウインドウでの描画を行う関数
cgl.plotR <- function(){
    if(menu$switch_pp$state==1){z <<- prey;zcol <<- pred;}
    else{z <<- pred;zcol <<- prey;}
    filled.contour(x,x,zcol,col=mypal(nlevs-1),xlab=Rxlab,ylab=Rylab, cex.lab=2.0,levels=seq(cgl.zlim[1],0.5*cgl.zlim[2],,nlevs),line=1.5,main=paste0("step: ",step, "  nlevs:",nlevs),
                   plot.axes = contour(x,x,z, levels =seq(zcol_min,zcol_max,,nlevs2), drawlabels = FALSE, axes = FALSE, frame.plot = FFALSE, add = TRUE,labcex=1.1,col=mypal2(nlevs2),lwd=clwd)
                   )   
}

 ##3Dプロットを行う関数（描画関数の戻り値を返すように設定する!）
cgl.plotRgl <- function(){
    if(menu$switch_pp$state==1){z <<- prey;zcol <<- pred;}
    else{z <<- pred;zcol <<- prey;}
    cv=as.integer(ncol*(zcol-zcol_min)/(zcol_max-zcol_min));
    cv[which(cv < 1)]=1;
    cv[which(cv > ncol)]=ncol;               
    return(terrain3d(x,x,z,color=col_pal[cv]));
}

step=0; ##時間ステップ
t=0.0; ##時刻を０にする
cgl.init(); ##初期化のための関数
cgl.plot(); ##初期状態の描画

while(1){
    if(flag_pause==0){##一時停止でないなら以下を実行する
        for(k in 1:view_interval){
            ##相互作用による変化率              
        dprey = -1.0*prey*pred*pred + a*(1.0-prey);
        dpred = prey*pred*pred - (a+b)*pred;
            ##移動量
        preydif =  (colDiffs(colDiffs(prey[difmask,]))+rowDiffs(rowDiffs(prey[,difmask])));
        preddif =  (colDiffs(colDiffs(pred[difmask,]))+rowDiffs(rowDiffs(pred[,difmask])));
         ##個体数の変化
        prey = prey + dt*(dprey + rate_prey*preydif);
        pred = pred + dt*(dpred + rate_pred*preddif);
    
        step=step+1;
        t=dt*step;
        if(step%%1000==0)cat("step: ",step,"\n");
        }

        cgl.plot(); ##描画         
    }
    
    Sys.sleep(0.001); ##0.001秒待機
    if(flag_halt==1)return(0); ## flag_haltが１ならば終了する    
}