install.packages("rgl")
source("rgl_demo3d.R")
library(rgl) ##rglをロード win_pos_x=100; ##ウィンドウの生成位置のx座標 win_pos_y=200; ##ウィンドウの生成位置のx座標 win_size=600; ##ウィンドウサイズ(ピクセル) flag_pause=0; ##0:一時停止しない, 1:一時停止 flag_halt=0; ##1で終了 dt=0.07; ##時間刻み幅 step=0; ##タイムステップ ncol=50; ##色パレットの色数 col_pal=topo.colors(ncol); ##色パレット zmax=2; ##色パレットの上限色に対応する値 zmin=-2; ##色パレットの下限色に対応する値 ##zの値から色を割り当てる関数 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]); } ##マウスクリックで一時停止と終了を制御する関数: ##マウスクリックのy座標がウインドウの上部から20%の範囲にある場合は ##終了(flag_halt=1) ##それ以外の場合は一時停止(flag_pauseを0⇒1 or 1⇒0) 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); ##-1〜1までのベクトルを生成 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の更新(ここにシミュレーション1ステップの実行関数をセットすればいい) z=sin(2*sqrt(xx^2+yy^2)+t*0.3)+cos(xx+1.3*yy+t); ##zの値を曲面表示して、そのオブジェクト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が1ならば終了する }
source("rgl_demo2d.R")
library(rgl) ##rglをロード win_pos_x=100; ##ウィンドウの生成位置のx座標 win_pos_y=200; ##ウィンドウの生成位置のx座標 win_size=600; ##ウィンドウサイズ(ピクセル) flag_pause=0; ##0:一時停止しない, 1:一時停止 flag_halt=0; ##1で終了 dt=0.02; ##時間刻み幅 step=0; ##タイムステップ zmax=2;##等高線の上限値 zmin=-2;##等高線の下限値 pwin_size=600;##2次元プロットの解像度(ピクセル) ##マウスクリックで一時停止と終了を制御する関数: ##マウスクリックのy座標がウインドウの上部から20%の範囲にある場合は ##終了(flag_halt=1) ##それ以外の場合は一時停止(flag_pauseを0⇒1 or 1⇒0) 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); ##視点をz軸に沿う下向き方向にセット view3d(theta=0,phi=0,fov=0,zoom=0.6); x=seq(-1,1,,64); ##-1〜1までのベクトルを生成 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);##等高線の太さを1に設定 ccol=rep("black",nlev);##等高線の色を黒にする ccol[8]="red";##8番目だけ赤くする clwd[8]=3.0;##8番目だけ太さを3にする ##塗り潰し等高線(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の更新(ここにシミュレーション1ステップの実行関数をセットすればいい) 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が1ならば終了する }
source("conrgl_demo2.R")
source("conrgl-0.12_with_demo2.R")
library(rgl); ##rglをロード source("conrgl-0.12.R"); ##conrgl-0.12.Rを読む cgl.plot_dim=3; ##最初に表示する描画の次元 cgl.xlim=c(-1,1); ##x軸の表示範囲 cgl.ylim=c(-1,1); ##y軸の表示範囲 cgl.zlim=c(-2,2); ##z軸の表示範囲 cgl.xlab=NULL; ##x軸のラベル cgl.ylab=NULL; ##y軸のラベル cgl.zlab=NULL; ##z軸のラベル 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; ##時刻を0にする x=seq(-1,1,,64); ##-1〜1までのベクトルを生成 xx=matrix(x, nrow=64, ncol=64, byrow=T); ##2次元格子のx座標 yy=matrix(x, nrow=64, ncol=64, byrow=F); ##2次元格子のy座標 z=sin(2*sqrt(xx^2+yy^2)+t*0.3)+cos(xx+1.3*yy+t); ##各位置のzの値を計算 cgl.plot(); ##初期状態の描画 while(1){ if(flag_pause==0){##一時停止でないなら以下を実行する step=step+1; t=dt*step; ##zの更新(ここにシミュレーション1ステップの実行関数をセットすればいい) 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が1ならば終了する }
source("conrgl_demo_branch1d.R")
source("conrgl-0.12_with_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); ##x軸の表示範囲 cgl.ylim=c(0,time_end); ##y軸の表示範囲 ##cgl.zlim=c(0,K0*2); ##z軸の表示範囲 cgl.zlim=c(-1,1); ##z軸の表示範囲 cgl.xlab="Trait x (niche position)"; ##x軸のラベル cgl.ylab="Time"; ##y軸のラベル cgl.zlab=NULL; ##z軸のラベル 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の位置に0でない個体数をセット 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);
install.packages("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); ##-1〜1までのベクトルを生成 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); ##移動の効果の計算に使う(反対側と連結) ##被食者の移動速度(xが大きい場所ほど速い) rate_prey=0.06*(5*(xx+1.1)); ##被食者の移動速度(yが大きい場所ほど速い) 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); ##x軸の表示範囲 cgl.ylim=c(-1,1); ##y軸の表示範囲 cgl.zlim=c(log(amp_prey),2*log(amp_prey*K)); ##z軸の表示範囲 cgl.xlab ="x";##x軸のラベル cgl.ylab ="y";##y軸のラベル cgl.zlab=NULL; ##z軸のラベル 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; ##時刻を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が1ならば終了する }
install.packages("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); ##-1〜1までのベクトルを生成 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から2までのベクトルを生成 ##被食者の移動速度(xが大きい場所ほど速い) rate_prey=matrix(0.08*buf, nrow=nn, ncol=nn, byrow=F); ##被食者の移動速度(yが大きい場所ほど速い) 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); ##x軸の表示範囲 cgl.ylim=c(-1,1); ##y軸の表示範囲 cgl.zlim=c(0,2); ##z軸の表示範囲 cgl.xlab = "x"; ##x軸のラベル cgl.ylab = "y"; ##y軸のラベル cgl.zlab=NULL; ##z軸のラベル nlevs=7; ##等高線の数(被食者) nlevs2=7; ##等高線の数(捕食者) mypal2=heat.colors; ## 2D描画で被食者を等高線で表示した場合の線色用のパレット clwd=rep(1.0,nlevs2); ##2D描画で被食者を等高線で表示した場合の線の太さ 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; ##時刻を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が1ならば終了する }