我正在编写一个使用 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 个单位(我假设单位是字节，但我不确定)。
.cov:
2 NN = length(xx) # number of data points
.mem:
191688 NN = length(xx) # number of data points
或者这个更糟:
冠状病毒:
352 @inbounds for ii = 1:NN
内存:
28000000 @inbounds for ii = 1:NN
最佳答案
简单提一下答案in the docs , “在用户设置下，直接从 REPL 调用的任何函数的第一行都将显示由于 REPL 代码本身发生的事件而分配。”也可能相关:“更重要的是，JIT 编译还增加了分配计数，因为 Julia 的大部分编译器都是用 Julia 编写的(编译通常需要内存分配)。推荐的过程是通过执行所有你想要的命令来强制编译分析，然后调用 Profile.clear_malloc_data() 重置所有分配计数器。"
底线:第一行被指责为其他地方发生的分配，因为它是重新开始报告分配的第一行。
关于memory - 如何分析 Julia 内存分配和代码覆盖结果，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/32520295/
我正在学习如何使用Nokogiri，根据这段代码我遇到了一些问题:require'rubygems'require'mechanize'post_agent=WWW::Mechanize.newpost_page=post_agent.get('http://www.vbulletin.org/forum/showthread.php?t=230708')puts"\nabsolutepathwithtbodygivesnil"putspost_page.parser.xpath('/html/body/div/div/div/div/div/table/tbody/tr/td/div
总的来说，我对ruby还比较陌生，我正在为我正在创建的对象编写一些rspec测试用例。许多测试用例都非常基础，我只是想确保正确填充和返回值。我想知道是否有办法使用循环结构来执行此操作。不必为我要测试的每个方法都设置一个assertEquals。例如:describeitem,"TestingtheItem"doit"willhaveanullvaluetostart"doitem=Item.new#HereIcoulddotheitem.name.shouldbe_nil#thenIcoulddoitem.category.shouldbe_nilendend但我想要一些方法来使用
作为我的Rails应用程序的一部分，我编写了一个小导入程序，它从我们的LDAP系统中吸取数据并将其塞入一个用户表中。不幸的是，与LDAP相关的代码在遍历我们的32K用户时泄漏了大量内存，我一直无法弄清楚如何解决这个问题。这个问题似乎在某种程度上与LDAP库有关，因为当我删除对LDAP内容的调用时，内存使用情况会很好地稳定下来。此外，不断增加的对象是Net::BER::BerIdentifiedString和Net::BER::BerIdentifiedArray，它们都是LDAP库的一部分。当我运行导入时，内存使用量最终达到超过1GB的峰值。如果问题存在，我需要找到一些方法来更正我的代
关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题？更新问题，以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点？我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等)，但我确实想创建一个输出文件。
给定这段代码defcreate@upgrades=User.update_all(["role=?","upgraded"],:id=>params[:upgrade])redirect_toadmin_upgrades_path,:notice=>"Successfullyupgradeduser."end我如何在该操作中实际验证它们是否已保存或未重定向到适当的页面和消息？ 最佳答案 在Rails3中，update_all不返回任何有意义的信息，除了已更新的记录数(这可能取决于您的DBMS是否返回该信息)。http://ar.ru
我在我的项目目录中完成了compasscreate.和compassinitrails。几个问题:我已将我的.sass文件放在public/stylesheets中。这是放置它们的正确位置吗？当我运行compasswatch时，它不会自动编译这些.sass文件。我必须手动指定文件:compasswatchpublic/stylesheets/myfile.sass等。如何让它自动运行？文件ie.css、print.css和screen.css已放在stylesheets/compiled。如何在编译后不让它们重新出现的情况下删除它们？我自己编译的.sass文件编译成compiled/t
我正在寻找执行以下操作的正确语法(在Perl、Shell或Ruby中):#variabletoaccessthedatalinesappendedasafileEND_OF_SCRIPT_MARKERrawdatastartshereanditcontinues. 最佳答案 Perl用__DATA__做这个:#!/usr/bin/perlusestrict;usewarnings;while(){print;}__DATA__Texttoprintgoeshere 关于ruby-如何将脚
如何在buildr项目中使用Ruby？我在很多不同的项目中使用过Ruby、JRuby、Java和Clojure。我目前正在使用我的标准Ruby开发一个模拟应用程序，我想尝试使用Clojure后端(我确实喜欢功能代码)以及JRubygui和测试套件。我还可以看到在未来的不同项目中使用Scala作为后端。我想我要为我的项目尝试一下buildr(http://buildr.apache.org/)，但我注意到buildr似乎没有设置为在项目中使用JRuby代码本身!这看起来有点傻，因为该工具旨在统一通用的JVM语言并且是在ruby中构建的。除了将输出的jar包含在一个独特的、仅限ruby
Rackup通过Rack的默认处理程序成功运行任何Rack应用程序。例如:classRackAppdefcall(environment)['200',{'Content-Type'=>'text/html'},["Helloworld"]]endendrunRackApp.new但是当最后一行更改为使用Rack的内置CGI处理程序时，rackup给出“NoMethodErrorat/undefinedmethod`call'fornil:NilClass”:Rack::Handler::CGI.runRackApp.newRack的其他内置处理程序也提出了同样的反对意见。例如Rack
在rails源中:https://github.com/rails/rails/blob/master/activesupport/lib/active_support/lazy_load_hooks.rb可以看到以下内容@load_hooks=Hash.new{|h,k|h[k]=[]}在IRB中，它只是初始化一个空哈希。和做有什么区别@load_hooks=Hash.new 最佳答案 查看rubydocumentationforHashnew→new_hashclicktotogglesourcenew(obj)→new_has