ABMEv_runworld.jl 8.16 KB
 Victor Boussange committed Jan 17, 2020 1 2 # for now we consider that competition is local within an array  Victor Boussange committed Mar 30, 2020 3 """  4  update_rates_std!(world,p::Dict,t::Float64)  Victor Boussange committed Mar 30, 2020 5 6 7 8 This standard updates takes - competition kernels of the form α(x,y) and - carrying capacity of the form K(x) """  Victor Boussange committed Mar 30, 2020 9 function update_rates_std!(world,p::Dict,t::Float64)  Victor Boussange committed Mar 30, 2020 10  α = p["alpha"];K=p["K"];  Victor Boussange committed Jan 17, 2020 11 12 13  traits = get_x.(world) # traits = get_xhist.(world) N = length(traits)  14  D = SharedArray{Float64}(N)  Victor Boussange committed Jan 17, 2020 15 16 17  # Here you should do a shared array to compute in parallel @sync @distributed for i in 1:(N-1) for j in i+1:N  Victor Boussange committed Mar 30, 2020 18  C = α(traits[i],traits[j])  19 20  D[i] += C D[j] += C  Victor Boussange committed Jan 17, 2020 21 22 23 24  end end # Here we can do it in parallel as well for (i,a) in enumerate(world)  25  a.d = D[i]  Victor Boussange committed Mar 30, 2020 26  a.b = K(traits[i])  Victor Boussange committed Jan 17, 2020 27 28 29 30 31 32  end end function update_rates_graph!(world,C,p::Dict,t::Float64) for e in edges(p["g"]) # agents on e  33  aidx_e = findall(a -> get_x(a,1)==e,world)  Victor Boussange committed Jan 17, 2020 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158  na = length(aidx_e) for i in aidx_e world[i].d = na^2 world[i].b = 1 end end end function update_rates_2D!(world,C,p::Dict,t::Float64) # Competition matrix traits = get_x.(world) # traits = get_xhist.(world) N = length(traits) # C = SharedArray{Float64}((N,N)) # Here you should do a shared array to compute in parallel @sync @distributed for i in 1:(N-1) C[i,i] = 1. for j in i+1:N # be careful, for xhist what is return is an array hence no need of [] C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"]) C[j,i] = C[i,j] end end C[N,N] = 1. # Here we can do it in parallel as well for (i,a) in enumerate(world) # we only update death rate # a.d = sum(C[i,:]) / K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1]) a.d = sum(C[i,:]) # /!| not ideal to assign at every time step the birth rate that is constant a.b = KK(traits[i][1:1], p["K0"], p["n_K"], p["sigma_K"][1:1], split_merge_move(t), -split_merge_move(t))/K(traits[i][2:2], p["K0"], p["n_K"], p["sigma_K"][2:2]) end end function update_rates_grad2D!(world,C,p::Dict,t::Float64) # Competition matrix traits = get_x.(world) # traits = get_xhist.(world) N = length(traits) # C = SharedArray{Float64}((N,N)) # Here you should do a shared array to compute in parallel @sync @distributed for i in 1:(N-1) C[i,i] = 1. for j in i+1:N # be careful, for xhist what is return is an array hence no need of [] C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"]) C[j,i] = C[i,j] end end C[N,N] = 1. # Here we can do it in parallel as well for (i,a) in enumerate(world) # we only update death rate a.d = sum(C[i,:]) a.b = K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1]) end end function update_rates_mountain!(world,C,p::Dict,t::Float64) # Competition matrix traits = get_x.(world) # traits = get_xhist.(world) N = length(traits) # C = SharedArray{Float64}((N,N)) # Here you should do a shared array to compute in parallel @sync @distributed for i in 1:(N-1) C[i,i] = 1. for j in i+1:N # be careful, for xhist what is return is an array hence no need of [] C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"]) C[j,i] = C[i,j] end end C[N,N] = 1. # Here we can do it in parallel as well for (i,a) in enumerate(world) # we only update death rate a.d = sum(C[i,:]) a.b = KK(traits[i][1:1], p["K0"], p["n_K"], p["sigma_K"][1:1], split_move(t), -split_move(t))/K(traits[i][2:2], p["K0"], p["n_K"], p["sigma_K"][2:2], μ=(tin(traits[i][1],-1.,0.)*(traits[i][1]+1) - tin(traits[i][1],0.,1.)*(traits[i][1]-1) ) * split_move(t)) end end function update_rates_std_split!(world,C,p::Dict,t::Float64) # Competition matrix traits = get_x.(world) # traits = get_xhist.(world) N = length(traits) # C = SharedArray{Float64}((N,N)) # Here you should do a shared array to compute in parallel @sync @distributed for i in 1:(N-1) C[i,i] = 1. for j in i+1:N C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"]) C[j,i] = C[i,j] end end C[N,N] = 1. # Here we can do it in parallel as well for (i,a) in enumerate(world) # we only update death rate # this could be imptrove since \alpha is symmetric, by using a symmetric matrix # a.d = sum(C[i,:]) / KK(traits[i][:,end],p["K0"],p["n_K"],p["sigma_K"],split_merge_move(t),-split_merge_move(t)) a.d = sum(C[i,:]) # /!| not ideal to assign at every time step the birth rate that is constant a.b = KK(traits[i][:,end],p["K0"],p["n_K"],p["sigma_K"],split_move(t),-split_move(t)) end end """  Victor committed Apr 02, 2020 159  function runWorld_store_WF(p,world0;mode="std")  Victor Boussange committed Jan 17, 2020 160 161 Wright Fisher process. Returns an array worldall with every agents. """  Victor committed Apr 02, 2020 162 function runWorld_store_WF(p,world0;mode="std")  Victor Boussange committed Jan 17, 2020 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183  worldall = repeat(world0,inner = (1,length(1:Int(p["tend"])))) N=length(world0); newworld = copy.(world0) if mode == "std" update_rates! = update_rates_std! elseif mode == "2D" update_rates! = update_rates_2D! elseif mode == "grad2D" update_rates! = update_rates_grad2D! elseif mode == "mountain" update_rates! = update_rates_mountain! elseif mode == "split" update_rates! = update_rates_std_split! elseif mode == "graph" update_rates! = update_rates_graph! else error("Mode $mode is not recognized") end for i in 1:(Int(p["tend"])-1) # we have to take views, otherwise it does not affect worldall world = @view worldall[:,i];newworld = @view worldall[:,i+1];  Victor committed Apr 02, 2020 184  updateWorld_WF!(world,newworld,p,update_rates!,Float64(i))  Victor Boussange committed Jan 17, 2020 185  end  Victor Boussange committed Jan 20, 2020 186  return worldall,collect(0:Int(p["tend"]-1))  Victor Boussange committed Jan 17, 2020 187 188 end """  Victor committed Apr 02, 2020 189  function runWorld_store_G(p,world0)  Victor Boussange committed Jan 17, 2020 190 Gillepsie process. Returns a tuple worldall,tspanarray  Victor committed Apr 02, 2020 191 192 If specified p["dt_saving"] determines the time step to save the simulation. If not only last time step is saved  Victor Boussange committed Jan 17, 2020 193 """  Victor committed Apr 02, 2020 194 function runWorld_store_G(p,world0)  Victor Boussange committed Jan 17, 2020 195 196  # we store the value of world every 100 changes by default tspan = zeros(1)  197  i = 1;j=0;dt = 1.  Victor Boussange committed Jan 17, 2020 198  N=length(world0);  Victor committed Apr 02, 2020 199 200 201  tspanarray = zeros(1); dt_saving = haskey(p,"dt_saving") ? p["dt_saving"] : p["tend"] + 1. # Array that stores the agents  Victor Boussange committed Jan 17, 2020 202  worldall = reshape(copy.(world0),N,1)  Victor Boussange committed Apr 01, 2020 203  # we instantiate C as the biggest size it can take  Victor Boussange committed Mar 30, 2020 204 205 206 207 208 209 210 211 212 213 214  update_rates_std!(skipmissing(world0),p,0.) while tspan[i]= dt_saving  Victor Boussange committed Jan 17, 2020 216 217  @info "saving world @ t =$(tspan[i])/ $(p["tend"])" j+=1;sw = size(worldall,2);  218 219  # we use <= because we save at the end of the wile loop if sw <= j  Victor Boussange committed Jan 17, 2020 220 221 222 223  # we double the size of worldall worldall = hcat(worldall,Array{Missing}(missing,N,sw)) end worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));  Victor Boussange committed Jan 17, 2020 224 225  push!(tspanarray,tspan[i]) end  Victor committed Apr 02, 2020 226  dt = updateWorld_G!(world0,p,update_rates_std!,tspan)  227  push!(tspan, tspan[end] + dt)  Victor Boussange committed Jan 17, 2020 228 229  i += 1 end  230 231 232 233  # Saving laste time step worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0))); push!(tspanarray,tspan[i]) @info "simulation stopped at t=$(tspanarray[end])"  Victor Boussange committed Jan 17, 2020 234 235  return worldall[:,1:j],tspanarray end