# A simple simulation of evolutionary branching in the MacArthur-Levins resource competition model. Written by Hiroshi C. Ito 2020.04.22. ##

# copyright (C) 2020 Hiroshi C. Ito
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License Version 2 as 
# published by the Free Software Foundation.
# http://www.r-project.org/Licenses/

K0=10.0; ##環境収容力の最大値
sK=1.0; ##環境収容力の幅
sa=0.5; ##競争効果の幅
##sa=0.7; 
##sa=1.2;
ndiv=64; ##形質軸の刻み幅
edge_die=1e-7; ##最小個体数密度。これ以下は絶滅。
repeat_p=1000; ##変異体を出現させずに個体数動態をまわすステップ数（速く計算するため）
dt=0.1; ##時間刻み幅
mutation_rate=3.0e-4; ##突然変異率
t_end=1500; ##総タイムステップ数
set.seed(23); ##乱数の種。これをコメントアウトすれば毎回少し違う進化動態になる

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

flag_quick=FALSE;
 ##環境収容力の分布
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+1),nrow=ndiv,ncol=(t_end+1)); ##出力用行列（個体数分布の時系列データ）
t_out=0.0*(1:(t_end+1)); ##出力用の時刻データ
time=0.0;  ##時刻
n_out[,1]=n;  ##出力用行列に初期状態を格納

t_exec_start=proc.time();
##ここからメインループ
for(t in 1:t_end){

    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){
        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以下の表現型は絶滅
     
    n_out[,t+1]=n;  ##出力用行列に格納
    t_out[t+1]=time;  ##出力用の時刻配列に格納

    ##edge_die以上の表現型を1としてプロンプトに表示
    if(t%%5==1){
        cat("x: ",paste0(as.integer(n>edge_die),collapse=""));  
        cat("  time:", time,"\n");  
    }
}


cat("calculation time:\n");
print(proc.time()-t_exec_start);
##10列おきにデータを使ってプロット
list=which(seq(length(t_out))%%10==1);  
filled.contour(xx,t_out[list],log(n_out[,list]+edge_die),xlab="Trait x",ylab="Time",main=paste0("sa:",sa, ", sK:",sK));
##dev.print(pdf,"test.pdf");  ##このコメントアウトを外せばtest.pdfとして画像が保存される（linuxだけ？）。