請幫我弄清楚為什么我的代碼很慢以及加快速度的可能方法。我使用了矢量化、@inbounds、列主要索引、@floop 和預分配,所以我認為它會更快。我很茫然...
該代碼模擬細胞(如生物細胞)和突變細胞的隨機波。我正在使用 Euler-Murayama 方法傳播耦合(化學朗之萬)方程:

其中W/M分別表示野生型細胞和突變細胞的數量,K是承載能力(最大細胞數),i表示deme(位置),N(0,1)是正常隨機變數。
這是經過一段時間傳播后的波浪圖形:

我附上了下面代碼的重要部分,并盡我所能對其進行評論。
using Random, Distributions
using StatsBase
using Statistics
using FLoops
# CLE Parameters/Set-Up:
K = 100 # Carrying capacity (maximum number of cells)
M = 100 # Number of demes (locations)
T = 100_000_000 # number of time steps
dt = 1e-1 # time increment
g = Normal(0.0,sqrt(dt)) # normal distribution with ave = 0.0, std_dev = dt
r_w = 0.1 # Wild-type growth rate
r_m = 0.2 # Mutant growth rate
r_wm = [r_w, r_m]' # Growth rate vector (transposed)
N = 1 # number of independent processes (slow even when N = 1)
# initial wave (essentially a step-function of wild-types
# with 100 mutants at deme (location) 76)
state_init = Matrix(reshape(repeat([K, 0.0]',M 2),(M 2,2)))
state_init[M÷2 2:end,1] .= 0
state_init[76,2] = 100.0
state_init[1,:] .= [K,0]
state_init[end,:] .= [0,0]
state = deepcopy(state_init)
state_plus = zeros(size(state_init)) # state at demes i 1 instead of i (used for derivatives)
state_minus = zeros(size(state_init)) # state at demes i-1 instead of i (used for derivatives)
function sim!(state_init::Matrix{Float64}, state::Matrix{Float64},
T::Int64, dt::Float64, N::Int64, M::Int64, K::Int64,
hist_data::Array{Int64,3}, g::Normal{Float64})
@inbounds @floop for n in 1:N
state .= deepcopy(state_init) # initialize state
@inbounds for t in 1:T
state_plus .= circshift(state, -1) # make plus state
state_plus[1,:] .= [0,0] # fix boundary conditions
state_minus .= circshift(state, 1) # make minus state
state_minus[end,:] .= [0,0] # fix boundary conditions
state_shift = circshift(state, (0,1))
# make state where each deme has a vector
# (# mutants, # wild-types) instead of
# (# wild-types, # mutants)
######################################
# propagate state using Euler-Murayama method and
# restrict number of cells in a deme to be in the range [0,K]
# using clamp(). clamp() also prevents imaginary numbers from
# a negative number under the sqrt().
state .= clamp.(state .
dt .* (r_wm .* state .* (K .- state .- state_shift) .
K .* (state_plus .- 2.0 .* state . state_minus)) .
sqrt.(clamp.(
6 .* state .* (K .- state) .
(K .- 2.0 .* state) .* (state_plus .- 2.0 .* state .
state_minus) .- r_wm .* state .* (K .- state .- state_shift),
0.0, 1.0*K*K)) .*
rand(g,M 2), 0.0, 1.0*K)
######################################
end
end
end
sim!(state_init, state, T, dt, N, M, K, hist_data, g)
... the rest of the code is analysis and not the reason the code is slow.
uj5u.com熱心網友回復:
即使您正在預分配state_plusand state_minus,由于您使用的是circshift而不是circshift!. circshift分配新記憶體,不管它最終被分配給現有的預分配陣列。進行 3 億次分配肯定是代價高昂的!
嘗試
function sim!(state_init::Matrix{Float64}, state::Matrix{Float64},
T::Int64, dt::Float64, N::Int64, M::Int64, K::Int64,
hist_data::Array{Int64,3}, g::Normal{Float64})
@inbounds @floop for n in 1:N
state .= deepcopy(state_init) # initialize state
state_plus = zeros(size(state_init)) # state at demes i 1 instead of i (used for derivatives)
state_minus = zeros(size(state_init)) # state at demes i-1 instead of i (used for derivatives)
state_shift = zeros(size(state_init))
for t in 1:T
circshift!(state_plus, state, -1) # make plus state
state_plus[1,:] .= [0,0] # fix boundary conditions
circshift!(state_minus, state, 1) # make minus state
state_minus[end,:] .= [0,0] # fix boundary conditions
circshift!(state_shift, state, (0,1))
和
function sim!(state_init::Matrix{Float64}, state::Matrix{Float64},
T::Int64, dt::Float64, N::Int64, M::Int64, K::Int64,
hist_data::Array{Int64,3}, g::Normal{Float64})
state_plus = zeros(size(state_init)) # state at demes i 1 instead of i (used for derivatives)
state_minus = zeros(size(state_init)) # state at demes i-1 instead of i (used for derivatives)
state_shift = zeros(size(state_init))
@inbounds for n in 1:N
state .= deepcopy(state_init) # initialize state
for t in 1:T
circshift!(state_plus, state, -1) # make plus state
state_plus[1,:] .= [0,0] # fix boundary conditions
circshift!(state_minus, state, 1) # make minus state
state_minus[end,:] .= [0,0] # fix boundary conditions
circshift!(state_shift, state, (0,1))
第二個版本不使用. @floop_ 最好兩者都試一試,然后找出答案!state_plusNN
轉載請註明出處,本文鏈接:https://www.uj5u.com/gongcheng/529658.html
下一篇:IIS性能-首次加載aspx頁面
