如何分析Julia内存分配和代码覆盖结果

Question

我正在用Gibbs抽样编写一个用于贝叶斯推断的包。由于这些方法通常在计算上很昂贵，所以我非常关心代码的性能。事实上，速度是我从Python转到Julia的原因。

在实现Dirichlet Process Model之后，我使用Coverage.jl和--track-allocation=user命令行选项分析了代码。

以下是覆盖结果

        - #=
        - DPM
        - 
        - Dirichlet Process Mixture Models
        - 
        - 25/08/2015
        - Adham Beyki, odinay@gmail.com
        - 
        - =#
        - 
        - type DPM{T}
        -   bayesian_component::T
        -   K::Int64
        -   aa::Float64
        -   a1::Float64
        -   a2::Float64
        -   K_hist::Vector{Int64}
        -   K_zz_dict::Dict{Int64, Vector{Int64}}
        - 
        -   DPM{T}(c::T, K::Int64, aa::Float64, a1::Float64, a2::Float64) = new(c, K, aa, a1, a2,
        -           Int64[], (Int64 => Vector{Int64})[])
        - end
        1 DPM{T}(c::T, K::Int64, aa::Real, a1::Real, a2::Real) = DPM{typeof(c)}(c, K, convert(Float64, aa),
        -           convert(Float64, a1), convert(Float64, a2))
        - 
        - function Base.show(io::IO, dpm::DPM)
        -   println(io, "Dirichlet Mixture Model with $(dpm.K) $(typeof(dpm.bayesian_component)) components")
        - end
        - 
        - function initialize_gibbs_sampler!(dpm::DPM, zz::Vector{Int64})
        -   # populates the cluster labels randomly
        1   zz[:] = rand(1:dpm.K, length(zz))
        - end
        - 
        - function DPM_sample_hyperparam(aa::Float64, a1::Float64, a2::Float64, K::Int64, NN::Int64, iters::Int64)
        - 
        -   # resampling concentration parameter based on Escobar and West 1995
      352   for n = 1:iters
     3504       eta = rand(Distributions.Beta(aa+1, NN))
     3504       rr = (a1+K-1) / (NN*(a2-log(NN)))
     3504       pi_eta = rr / (1+rr)
        - 
     3504       if rand() < pi_eta
        0           aa = rand(Distributions.Gamma(a1+K, 1/(a2-log(eta))))
        -       else
     3504           aa = rand(Distributions.Gamma(a1+K-1, 1/(a2-log(eta))))
        -       end
        -   end
      352   aa
        - end
        - 
        - function DPM_sample_pp{T1, T2}(
        -     bayesian_components::Vector{T1},
        -     xx::T2,
        -     nn::Vector{Float64},
        -     pp::Vector{Float64},
        -     aa::Float64)
        - 
  1760000   K = length(nn)
  1760000   @inbounds for kk = 1:K
 11384379     pp[kk] = log(nn[kk]) + logpredictive(bayesian_components[kk], xx)
        -   end
  1760000   pp[K+1] = log(aa) + logpredictive(bayesian_components[K+1], xx)
  1760000   normalize_pp!(pp, K+1)
  1760000   return sample(pp[1:K+1])
        - end
        - 
        - 
        - function collapsed_gibbs_sampler!{T1, T2}(
        -       dpm::DPM{T1},
        -       xx::Vector{T2},
        -       zz::Vector{Int64},
        -       n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64; max_clusters::Int64=100)
        - 
        - 
        2   NN = length(xx)                                          # number of data points
        2   nn = zeros(Float64, dpm.K)                       # count array
        2   n_iterations = n_burnins + (n_samples)*(n_lags+1)
        2   bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:dpm.K+1]
        2   dpm.K_hist = zeros(Int64, n_iterations)
        2   pp = zeros(Float64, max_clusters)
        - 
        2   tic()
        2   for ii = 1:NN
    10000       kk = zz[ii]
    10000       additem(bayesian_components[kk], xx[ii])
    10000       nn[kk] += 1
        -   end
        2   dpm.K_hist[1] = dpm.K
        2   elapsed_time = toq()
        - 
        2   for iteration = 1:n_iterations
        - 
      352       println("iteration: $iteration, KK: $(dpm.K), KK mode: $(indmax(hist(dpm.K_hist,
        -                       0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time")
        - 
      352       tic()
      352       @inbounds for ii = 1:NN
  1760000           kk = zz[ii]
  1760000           nn[kk] -= 1
  1760000           delitem(bayesian_components[kk], xx[ii])
        - 
        -           # remove the cluster if empty
  1760000           if nn[kk] == 0
      166               println("\tcomponent $kk has become inactive")
      166               splice!(nn, kk)
      166               splice!(bayesian_components, kk)
      166               dpm.K -= 1
        - 
        -               # shifting the labels one cluster back
   830166               idx = find(x -> x>kk, zz)
      166               zz[idx] -= 1
        -           end
        - 
  1760000           kk = DPM_sample_pp(bayesian_components, xx[ii], nn, pp, dpm.aa)
        - 
  1760000           if kk == dpm.K+1
      171               println("\tcomponent $kk activated.")
      171               push!(bayesian_components, deepcopy(dpm.bayesian_component))
      171               push!(nn, 0)
      171               dpm.K += 1
        -           end
        - 
  1760000           zz[ii] = kk
  1760000           nn[kk] += 1
  1760000           additem(bayesian_components[kk], xx[ii])
        -       end
        - 
      352       dpm.aa = DPM_sample_hyperparam(dpm.aa, dpm.a1, dpm.a2, dpm.K, NN, n_internals)
      352       dpm.K_hist[iteration] = dpm.K
      352       dpm.K_zz_dict[dpm.K] = deepcopy(zz)
      352       elapsed_time = toq()
        -   end
        - end
        - 
        - function truncated_gibbs_sampler{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, zz::Vector{Int64},
        -   n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64, K_truncation::Int64)
        - 
        -   NN = length(xx)                                             # number of data points
        -   nn = zeros(Int64, K_truncation)             # count array
        -   bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:K_truncation]
        -   n_iterations = n_burnins + (n_samples)*(n_lags+1)
        -   dpm.K_hist = zeros(Int64, n_iterations)
        -   states = (ASCIIString => Int64)[]
        -   n_states = 0
        - 
        -   tic()
        -   for ii = 1:NN
        -       kk = zz[ii]
        -       additem(bayesian_components[kk], xx[ii])
        -       nn[kk] += 1
        -   end
        -   dpm.K_hist[1] = dpm.K
        - 
        -   # constructing the sticks
        -   beta_VV = rand(Distributions.Beta(1.0, dpm.aa), K_truncation)
        -   beta_VV[end] = 1.0
        -   π = ones(Float64, K_truncation)
        -   π[2:end] = 1 - beta_VV[1:K_truncation-1]
        -   π = log(beta_VV) + log(cumprod(π))
        - 
        -   elapsed_time = toq()
        - 
        -   for iteration = 1:n_iterations
        - 
        -       println("iteration: $iteration, # active components: $(length(findn(nn)[1])), mode: $(indmax(hist(dpm.K_hist,
        -                       0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time \n", nn)
        - 
        -       tic()
        -       for ii = 1:NN
        -           kk = zz[ii]
        -           nn[kk] -= 1
        -           delitem(bayesian_components[kk], xx[ii])
        - 
        -           # resampling label
        -           pp = zeros(Float64, K_truncation)
        -           for kk = 1:K_truncation
        -               pp[kk] = π[kk] + logpredictive(bayesian_components[kk], xx[ii])
        -           end
        -           pp = exp(pp - maximum(pp))
        -           pp /= sum(pp)
        - 
        -           # sample from pp
        -           kk = sampleindex(pp)
        -           zz[ii] = kk
        -           nn[kk] += 1
        -           additem(bayesian_components[kk], xx[ii])
        - 
        -           for kk = 1:K_truncation-1
        -               gamma1 = 1 + nn[kk]
        -               gamma2 = dpm.aa + sum(nn[kk+1:end])
        -               beta_VV[kk] = rand(Distributions.Beta(gamma1, gamma2))
        -           end
        -           beta_VV[end] = 1.0
        -           π = ones(Float64, K_truncation)
        -           π[2:end] = 1 - beta_VV[1:K_truncation-1]
        -           π = log(beta_VV) + log(cumprod(π))
        - 
        -           # resampling concentration parameter based on Escobar and West 1995
        -           for internal_iters = 1:n_internals
        -                   eta = rand(Distributions.Beta(dpm.aa+1, NN))
        -               rr = (dpm.a1+dpm.K-1) / (NN*(dpm.a2-log(NN)))
        -               pi_eta = rr / (1+rr)
        - 
        -               if rand() < pi_eta
        -                   dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K, 1/(dpm.a2-log(eta))))
        -               else
        -                   dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K-1, 1/(dpm.a2-log(eta))))
        -               end
        -           end
        -       end
        - 
        -       nn_string = nn2string(nn)
        -       if !haskey(states, nn_string)
        -           n_states += 1
        -           states[nn_string] = n_states
        -       end
        -       dpm.K_hist[iteration] = states[nn_string]
        -       dpm.K_zz_dict[states[nn_string]] = deepcopy(zz)
        -       elapsed_time = toq()
        -   end
        -   return states
        - end
        - 
        - 
        - function posterior{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, K::Int64, K_truncation::Int64=0)
        2   n_components = 0
        1   if K_truncation == 0
        1       n_components = K
        -   else
        0       n_components = K_truncation
        -   end
        - 
        1   bayesian_components = [deepcopy(dpm.bayesian_component) for kk=1:n_components]
        1   zz = dpm.K_zz_dict[K]
        - 
        1   NN = length(xx)
        1   nn = zeros(Int64, n_components)
        - 
        1   for ii = 1:NN
     5000       kk = zz[ii]
     5000       additem(bayesian_components[kk], xx[ii])
     5000       nn[kk] += 1
        -   end
        - 
        1   return([posterior(bayesian_components[kk]) for kk=1:n_components], nn)
        - end
        - 

下面是内存分配：

        - #=
        - DPM
        - 
        - Dirichlet Process Mixture Models
        - 
        - 25/08/2015
        - Adham Beyki, odinay@gmail.com
        - 
        - =#
        - 
        - type DPM{T}
        -   bayesian_component::T
        -   K::Int64
        -   aa::Float64
        -   a1::Float64
        -   a2::Float64
        -   K_hist::Vector{Int64}
        -   K_zz_dict::Dict{Int64, Vector{Int64}}
        - 
        -   DPM{T}(c::T, K::Int64, aa::Float64, a1::Float64, a2::Float64) = new(c, K, aa, a1, a2,
        -           Int64[], (Int64 => Vector{Int64})[])
        - end
        0 DPM{T}(c::T, K::Int64, aa::Real, a1::Real, a2::Real) = DPM{typeof(c)}(c, K, convert(Float64, aa),
        -           convert(Float64, a1), convert(Float64, a2))
        - 
        - function Base.show(io::IO, dpm::DPM)
        -   println(io, "Dirichlet Mixture Model with $(dpm.K) $(typeof(dpm.bayesian_component)) components")
        - end
        - 
        - function initialize_gibbs_sampler!(dpm::DPM, zz::Vector{Int64})
        -   # populates the cluster labels randomly
        0   zz[:] = rand(1:dpm.K, length(zz))
        - end
        - 
        - function DPM_sample_hyperparam(aa::Float64, a1::Float64, a2::Float64, K::Int64, NN::Int64, iters::Int64)
        - 
        -   # resampling concentration parameter based on Escobar and West 1995
        0   for n = 1:iters
        0       eta = rand(Distributions.Beta(aa+1, NN))
        0       rr = (a1+K-1) / (NN*(a2-log(NN)))
        0       pi_eta = rr / (1+rr)
        - 
        0       if rand() < pi_eta
        0           aa = rand(Distributions.Gamma(a1+K, 1/(a2-log(eta))))
        -       else
        0           aa = rand(Distributions.Gamma(a1+K-1, 1/(a2-log(eta))))
        -       end
        -   end
        0   aa
        - end
        - 
        - function DPM_sample_pp{T1, T2}(
        -     bayesian_components::Vector{T1},
        -     xx::T2,
        -     nn::Vector{Float64},
        -     pp::Vector{Float64},
        -     aa::Float64)
        - 
        0   K = length(nn)
        0   @inbounds for kk = 1:K
        0     pp[kk] = log(nn[kk]) + logpredictive(bayesian_components[kk], xx)
        -   end
        0   pp[K+1] = log(aa) + logpredictive(bayesian_components[K+1], xx)
        0   normalize_pp!(pp, K+1)
        0   return sample(pp[1:K+1])
        - end
        - 
        - 
        - function collapsed_gibbs_sampler!{T1, T2}(
        -       dpm::DPM{T1},
        -       xx::Vector{T2},
        -       zz::Vector{Int64},
        -       n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64; max_clusters::Int64=100)
        - 
        - 
   191688   NN = length(xx)                                          # number of data points
       96   nn = zeros(Float64, dpm.K)                       # count array
        0   n_iterations = n_burnins + (n_samples)*(n_lags+1)
      384   bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:dpm.K+1]
     2864   dpm.K_hist = zeros(Int64, n_iterations)
      176   pp = zeros(Float64, max_clusters)
        - 
       48   tic()
        0   for ii = 1:NN
        0       kk = zz[ii]
        0       additem(bayesian_components[kk], xx[ii])
        0       nn[kk] += 1
        -   end
        0   dpm.K_hist[1] = dpm.K
        0   elapsed_time = toq()
        - 
        0   for iteration = 1:n_iterations
        - 
  5329296       println("iteration: $iteration, KK: $(dpm.K), KK mode: $(indmax(hist(dpm.K_hist,
        -                       0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time")
        - 
    16800       tic()
 28000000       @inbounds for ii = 1:NN
        0           kk = zz[ii]
        0           nn[kk] -= 1
        0           delitem(bayesian_components[kk], xx[ii])
        - 
        -           # remove the cluster if empty
        0           if nn[kk] == 0
   161880               println("\tcomponent $kk has become inactive")
        0               splice!(nn, kk)
        0               splice!(bayesian_components, kk)
        0               dpm.K -= 1
        - 
        -               # shifting the labels one cluster back
    69032               idx = find(x -> x>kk, zz)
    42944               zz[idx] -= 1
        -           end
        - 
        0           kk = DPM_sample_pp(bayesian_components, xx[ii], nn, pp, dpm.aa)
        - 
        0           if kk == dpm.K+1
   158976               println("\tcomponent $kk activated.")
    14144               push!(bayesian_components, deepcopy(dpm.bayesian_component))
     4872               push!(nn, 0)
        0               dpm.K += 1
        -           end
        - 
        0           zz[ii] = kk
        0           nn[kk] += 1
        0           additem(bayesian_components[kk], xx[ii])
        -       end
        - 
        0       dpm.aa = DPM_sample_hyperparam(dpm.aa, dpm.a1, dpm.a2, dpm.K, NN, n_internals)
        0       dpm.K_hist[iteration] = dpm.K
 14140000       dpm.K_zz_dict[dpm.K] = deepcopy(zz)
        0       elapsed_time = toq()
        -   end
        - end
        - 
        - function truncated_gibbs_sampler{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, zz::Vector{Int64},
        -   n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64, K_truncation::Int64)
        - 
        -   NN = length(xx)                                             # number of data points
        -   nn = zeros(Int64, K_truncation)             # count array
        -   bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:K_truncation]
        -   n_iterations = n_burnins + (n_samples)*(n_lags+1)
        -   dpm.K_hist = zeros(Int64, n_iterations)
        -   states = (ASCIIString => Int64)[]
        -   n_states = 0
        - 
        -   tic()
        -   for ii = 1:NN
        -       kk = zz[ii]
        -       additem(bayesian_components[kk], xx[ii])
        -       nn[kk] += 1
        -   end
        -   dpm.K_hist[1] = dpm.K
        - 
        -   # constructing the sticks
        -   beta_VV = rand(Distributions.Beta(1.0, dpm.aa), K_truncation)
        -   beta_VV[end] = 1.0
        -   π = ones(Float64, K_truncation)
        -   π[2:end] = 1 - beta_VV[1:K_truncation-1]
        -   π = log(beta_VV) + log(cumprod(π))
        - 
        -   elapsed_time = toq()
        - 
        -   for iteration = 1:n_iterations
        - 
        -       println("iteration: $iteration, # active components: $(length(findn(nn)[1])), mode: $(indmax(hist(dpm.K_hist,
        -                       0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time \n", nn)
        - 
        -       tic()
        -       for ii = 1:NN
        -           kk = zz[ii]
        -           nn[kk] -= 1
        -           delitem(bayesian_components[kk], xx[ii])
        - 
        -           # resampling label
        -           pp = zeros(Float64, K_truncation)
        -           for kk = 1:K_truncation
        -               pp[kk] = π[kk] + logpredictive(bayesian_components[kk], xx[ii])
        -           end
        -           pp = exp(pp - maximum(pp))
        -           pp /= sum(pp)
        - 
        -           # sample from pp
        -           kk = sampleindex(pp)
        -           zz[ii] = kk
        -           nn[kk] += 1
        -           additem(bayesian_components[kk], xx[ii])
        - 
        -           for kk = 1:K_truncation-1
        -               gamma1 = 1 + nn[kk]
        -               gamma2 = dpm.aa + sum(nn[kk+1:end])
        -               beta_VV[kk] = rand(Distributions.Beta(gamma1, gamma2))
        -           end
        -           beta_VV[end] = 1.0
        -           π = ones(Float64, K_truncation)
        -           π[2:end] = 1 - beta_VV[1:K_truncation-1]
        -           π = log(beta_VV) + log(cumprod(π))
        - 
        -           # resampling concentration parameter based on Escobar and West 1995
        -           for internal_iters = 1:n_internals
        -                   eta = rand(Distributions.Beta(dpm.aa+1, NN))
        -               rr = (dpm.a1+dpm.K-1) / (NN*(dpm.a2-log(NN)))
        -               pi_eta = rr / (1+rr)
        - 
        -               if rand() < pi_eta
        -                   dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K, 1/(dpm.a2-log(eta))))
        -               else
        -                   dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K-1, 1/(dpm.a2-log(eta))))
        -               end
        -           end
        -       end
        - 
        -       nn_string = nn2string(nn)
        -       if !haskey(states, nn_string)
        -           n_states += 1
        -           states[nn_string] = n_states
        -       end
        -       dpm.K_hist[iteration] = states[nn_string]
        -       dpm.K_zz_dict[states[nn_string]] = deepcopy(zz)
        -       elapsed_time = toq()
        -   end
        -   return states
        - end
        - 
        - 
        - function posterior{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, K::Int64, K_truncation::Int64=0)
        0   n_components = 0
        0   if K_truncation == 0
        0       n_components = K
        -   else
        0       n_components = K_truncation
        -   end
        - 
        0   bayesian_components = [deepcopy(dpm.bayesian_component) for kk=1:n_components]
        0   zz = dpm.K_zz_dict[K]
        - 
        0   NN = length(xx)
        0   nn = zeros(Int64, n_components)
        - 
        0   for ii = 1:NN
        0       kk = zz[ii]
        0       additem(bayesian_components[kk], xx[ii])
        0       nn[kk] += 1
        -   end
        - 
        0   return([posterior(bayesian_components[kk]) for kk=1:n_components], nn)
        - end
        - 

我似乎不明白的是，为什么一个简单的任务只运行两次，分配191688个单元(我假设这个单位是Bytes，我不确定)。

.cov：

    2   NN = length(xx)                  # number of data points

.mem：

   191688   NN = length(xx)              # number of data points

或者这个更糟：

cov：

  352       @inbounds for ii = 1:NN

mem：

  28000000      @inbounds for ii = 1:NN

Answer 1

简单地提到了in the docs，“在用户设置下，从REPL直接调用的任何函数的第一行将显示由于REPL代码本身发生的事件而进行的分配。”也可能相关：“更重要的是，JIT编译也增加了分配计数，因为大部分JIT编译器都是用Julia编写的(编译通常需要内存分配)。推荐的过程是通过执行所有要分析的命令强制编译，然后调用Profile.clear_malloc_data()来重置所有分配计数器。”

底线:第一行被指责为其他地方的分配，因为它是第一条重新开始报告分配的行。