Here’ s the MWE. Hopefully it won’t bring trouble to others’ laptop 
#minimum working example
using Distributions #Multinomial()
using Random #Multinomial()
using PoissonRandom #pois_rand
using LazyArrays #ApplyArray(vcat,)
#For test the function
Nmax=5.5e7
Nmin=5.5e6
D= gaussian_1d_localmig3(Nmax,Nmin,1,0.1,5e-5,5) #8 seconds
# plot(D[2]',yscale=:log10,ylims=(1e-3,1),legend=false)
# savefig("..\\photo\\oscipopu_without migration")
function gaussian_1d_localmig3(Nmax, Nmin,twoNmu,s0,mig,nDemes)
T=200
## Preallocating dataframes
migCell=Array{Array{Float64,1},1}(undef,nDemes) # stores how many migrants per deme
#one column is enough
migLR=Array{Array{Float64,2},1}(undef,nDemes) #left right migration amount
D=Array{Array{Float64,2},1}(undef,nDemes)#deme cell array - each with their own frequency array
x = zeros(1,T+1)
x[end,1] = 1 #Initial frequency array of the wt (no mutants present yet)
for i in 1:nDemes # populate deme cell array with initial wt-frequency of 1 @ T= 0
D[i] = x
end
#enable input s0, mig instead of population-scaled parameter twoNs, twoNmig
N=2(Nmax^(-1)+Nmin^(-1))^-1
twoNs=s0*N*2
twoNmig=mig*N*2
for t in 2:T+1 #Generations #T+1
println("generation############",t)
prd=12 #a period of 12 generation
NDeme=0.5(Nmax+Nmin)+0.5(Nmax+Nmin)*sin(2pi*t/prd) #oscillating population size
if NDeme==0 #usually in (9+12n) generation
NDeme=200 #to avoid casuing NANs after being dividing NDeme
end
NDeme=round(Int,NDeme)
println("oscillating NDeme ", NDeme)
## Population Parameters
mu = twoNmu/2/NDeme/nDemes # mutation coefficient
s0 = twoNs/2/NDeme/nDemes #selection coefficient, s_ to avoid name crash with s in Julia
mig = twoNmig/2/NDeme/nDemes # dispersal rate
for i in 1:nDemes
println("deme#####", i)
x=D[i]
nan=findall(isnan,x)
println("nan in x_beginning of Selection ",nan) #check the input at the beginning
if !isempty(nan)
println("x ",D[i][:,t-1])
end
## Filling up the frequency array
##Selection
##Selection, genetic drift
xwt = x[end,t-1] #wild type frequency at t-1
X = x[1:end-1,t-1] # mutant frequencies at t-1
ss=xwt*s0 #selection modification
xx=X+ss*X
xx[xx.<0].=0
xxK = 1-sum(xx) #Frequency of wild type (wt is the Kth allele)
if xxK < 0 # in case wild type frequency is negative
xxK = 0 # set wt frequency to 0
xx = xx/sum(xx) #resets all mutant frequencies to sum to 1
end
## Genetic drift
println("NDeme",NDeme)
n = rand(Multinomial(Int(NDeme),[xx;xxK])) #Multinomial random numbers
z=n/NDeme
x[:,t] = z #n outputed in individuals
D[i]=x
nan=findall(isnan,x)
println("nan in x_end of drift ",nan)
#[3:7,9], problem in frequency vector in 9th generation
if !isempty(nan)
println("x ",D[i][:,t])
end
##Mutation
xwt = x[end,t] #wild type frequency at t-1
m = pois_rand(NDeme*mu*xwt) #number of de novo mutants generated
if m != 0 # with >= mutant alleles generated
xN = zeros(m,T+1) #rows of new mutant alleles
for f in 1:nDemes #add xN to all demes
x = D[f]
x=ApplyArray(vcat,xN,x) #new mutant piled on top
if f==i #frequency of the deme that mutant emerged as 1/NDeme
x[1:m,t] .= 1/NDeme # 1 individual arising per mutant seeding into above array
x[end,t] = x[end,t] - m/NDeme
end
D[f]=x
end
end
# x=D[i]#x still in this module
nan=findall(isnan,x)
println("nan_end of mutation ",nan)
#[3:7,9], problem in frequency vector in 9th generation
if !isempty(nan)
println("x ",D[i][:,t])
end
end # for i in 1:nDemes
D = NDeme*D #frequency=>individuals
##Migration
## Determining the number of migrants leaving the demes
for i in 1:nDemes
n_i = D[i] #extract from cell array, multi-step=>better computational efficiency
n_B = n_i[:,t] #frequency column in t generation, i-th deme
println("n_B _before MIG",n_B)
MIG = mig*n_B #expected number of migrants individuals #column vec
migLV = pois_rand.(MIG) #actual number of migrants
#Before leaving in migration
#column vector
for g in 1:length(n_B) # poisson may generate>there is=>negative frequency
#"migCell[i] = poissrnd(MIG)' #actual number of migrants"--line111
if migLV[g] > n_B[g] # looks for negative values and sets to zero
println("migLV[g] ", migLV[g])
println("n_B[g] ", n_B[g])
migLV[g] = trunc(n_B[g]) #both float vector
end
end
n_lv = n_B - migLV #after leaving; both column vector
n_i[:,t] = n_lv # storing population after migrants leave in
# println("migration_leaving_size(n_lv)",size(n_lv))
D[i] = n_i # stores "frequency" array back in to deme cell array
migCell[i] = migLV # leaving array! not left-right array (migLR)
# println("migLV ", migLV)
nan=findall(isnan,migLV)
println("nan_end of migLV ",nan)#NaN in which column=>which generation
end
## Determining if migrants move left or right
for i in 1:nDemes
c = size(migCell[i])[1] #c=column number, registers size of migrant
#array to preallocate size of migLR_i array
migLR_i = Int.(zeros(c,2)) #column 1 amount going left, 2 going right
tot_mig = migCell[i] #total: left+right
# println("tot_mig ", " max " maximum(tot_mig)," type ",typeof(tot_mig))
tot_mig = Int.(tot_mig)
if i == 1 #first deme
# println("first_size",size(tot_mig))
migLR_i[:,2] = tot_mig # left end, only toward right
migLR[i] = migLR_i
elseif i == nDemes #last deme
# println("last_size", size(tot_mig))
migLR_i[:,1] = tot_mig # right end, only toward left
migLR[i] = migLR_i
else # middle demes
# migLR_i(1,:) = binornd(tot_mig,0.5);
# println("middle_size", size(tot_mig))
# println("mig")
# println("middle_size",size(tot_mig))
migLR_i[:,1] = rand.(Binomial.(tot_mig,0.5))#toward left
#two elemetwise dot is needed
migLR_i[:,2] = tot_mig - migLR_i[:,1] #toward right
migLR[i] = migLR_i
end
nan=findall(isnan,migLR[i])
println("nan position_end of migLR ",nan)#NaN in which column=>which generation
#[3:7,9], problem in frequency vector in 9th generation
end
#migLR: column1 column2
#leaving: toward left toward right
#entering: fr Right fr Left
## Accounting for migrants entering in the "frequency" array, D
for i in 1:nDemes
if i == 1 # for first deme
n_i = D[i] # n before accounting for entering migrants
n_B = n_i[:,t]
n_frR = migLR[i+1] # n entering array !ROW VECTOR
n_frR = n_frR[:,1] #only migration from right, columnn vector
#column vector
println("size(n_B)",size(n_B))
println("size(n_frR)",size(n_frR))
n_aft = n_B + n_frR
n_i[:,t]=n_aft
D[i] = n_i # stores back into deme cell array
elseif i == nDemes # for last deme
n_i = D[i]
n_B = n_i[:,t]
n_frL = migLR[i-1] # n entering array !ROW VECTOR
n_frL = n_frL[:,2] #only migration from left, columnn vector
println("size(n_B)",size(n_B))
println("size(n_frL)",size(n_frL))
n_aft = n_B + n_frL # adds entering array into existing array
n_i[:,t]=n_aft
D[i] = n_i # stores back into deme cell array
else # for the middle demes
n_i = D[i] # n before entering, COLUMN VECTOR
n_B = n_i[:,t]
# println("migLR_size",size(migLR))
n_frR = migLR[i+1] # n entering array from right deme !ROW VECTOR
n_frR = n_frR[:,1] # n entering (only left migration therefore #row 1)!COLUMN VECTOR
n_frL = migLR[i-1] # n entering array from left deme !ROW VECTOR
n_frL = n_frL[:,2] # n entering (only right migration therefore
#row 2)!COLUMN VECTOR
println("size(n_B)",size(n_B),
"size(n_frL)",size(n_frL),
"size(n_frR))",size(n_frR))
n_aft = n_B + n_frR + n_frL # adds entering array into existing array
n_i[:,t]=n_aft
D[i] = n_i # stores back into deme cell array
end
nan=findall(isnan,n_i)
println("nan position_just after migration ",nan)
println("n_i just after migration ", n_aft)
end
#individuals=>frequency
for i in 1:nDemes
x_i=D[i]#column in frequency except in this generation
n_B=x_i[:,t]
x_aft=n_B/sum(n_B)
println("sum(n_B) ",sum(n_B))
x_i[:,t]=x_aft
D[i]=x_i
#check NaNs contained at the end of migration or not
nan=findall(isnan,x_i)
println("x_i ", size(x_aft))#how much elements in this column
println("nan position_after individuals=>frequency ",nan)
#how much NaNs in this column
end
end #for t = 2:T+1
return D
end #function
