From 822fd0bc365abc2fe27a5bc1171f4a3bb2991d95 Mon Sep 17 00:00:00 2001 From: "David A. van Leeuwen" Date: Sat, 22 Nov 2014 19:51:57 +0100 Subject: [PATCH] =?UTF-8?q?Cosmetic=20changes,=20nx=20->=20n=E2=82=93=20in?= =?UTF-8?q?=20comments=20*=20->=20=C3=97=20(this=20is=20\times=20not=20x,?= =?UTF-8?q?=20boy=20are=20we=20pedantic)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/bayes.jl | 56 ++++++++++++++++++------------------- src/data.jl | 6 ++-- src/rand.jl | 6 ++-- src/recognizer.jl | 6 ++-- src/stats.jl | 70 +++++++++++++++++++++++------------------------ src/train.jl | 58 +++++++++++++++++++-------------------- 6 files changed, 101 insertions(+), 101 deletions(-) diff --git a/src/bayes.jl b/src/bayes.jl index 0f29b8b..6a93a27 100644 --- a/src/bayes.jl +++ b/src/bayes.jl @@ -17,10 +17,10 @@ function GMMprior{T<:FloatingPoint}(d::Int, alpha::T, beta::T) end Base.copy(p::GMMprior) = GMMprior(p.α0, p.β0, copy(p.m0), p.ν0, copy(p.W0)) -## initialize from a GMM and nx, the number of points used to train the GMM. +## initialize from a GMM and nₓ, the number of points used to train the GMM. function VGMM{T}(g::GMM{T}, π::GMMprior{T}) - nx = g.nx - N = g.w * nx + nₓ = g.nx + N = g.w * nₓ mx = g.μ if kind(g) == :diag S = covars(full(g)) @@ -60,8 +60,8 @@ function mstep{T}(π::GMMprior, N::Vector{T}, mx::Matrix{T}, S::Vector) α = π.α0 + N # ng, 10.58 ν = π.ν0 + N + 1 # ng, 10.63 β = π.β0 + N # ng, 10.60 - m = similar(mx) # ng * d - W = Array(eltype(FullCov{T}), ng) # ng * (d*d) + m = similar(mx) # ng × d + W = Array(eltype(FullCov{T}), ng) # ng × (d*d) d = size(mx,2) limit = √ eps(eltype(N)) invW0 = inv(π.W0) @@ -94,16 +94,16 @@ end ## log(ρ_nk) from 10.46, start with a very slow implementation, and optimize it further ## This is as the heart of VGMM training, it should be super efficient, possibly ## at the cost of readability. -## Removing the inner loop over nx, replacing it by a matmul led to a major speed increase +## Removing the inner loop over nₓ, replacing it by a matmul led to a major speed increase ## not leading to further speed increase for (Δ * vg.W[k]) .* Δ ## - c = Δ * chol(vg.W[k],:L), c .* c ## - Base.BLAS.symm('R', 'U', vg.ν[k], vg.W[k], Δ) .* Δ function logρ(vg::VGMM, x::Matrix, ex::Tuple) - (nx, d) = size(x) + (nₓ, d) = size(x) d == vg.d || error("dimension mismatch") ng = vg.n - EμΛ = similar(x, nx, ng) # nx * ng - Δ = similar(x) # nx * d + EμΛ = similar(x, nₓ, ng) # nₓ × ng + Δ = similar(x) # nₓ × d for k=1:ng ### d/vg.β[k] + vg.ν[k] * (x_i - m_k)' W_k (x_i = m_k) forall i ## Δ = (x_i - m_k)' W_k (x_i = m_k) @@ -118,7 +118,7 @@ end function rnk(vg::VGMM, x::Matrix, ex::Tuple) # ρ = exp(logρ(g, x)) # broadcast(/, ρ, sum(ρ, 2)) - logr = logρ(vg, x, ex) # nx * ng + logr = logρ(vg, x, ex) # nₓ × ng broadcast!(-, logr, logr, logsumexp(logr, 2)) ## we slowly replace logr by r, this is just cosmetic! r = logr @@ -138,29 +138,29 @@ function rnk(vg::VGMM, x::Matrix, ex::Tuple) end ## OK, also do stats. -## Like for the GMM, we return nx, (some value), zeroth, first, second order stats +## Like for the GMM, we return nₓ, (some value), zeroth, first, second order stats ## All return values can be accumulated, except r, which we need for ## lowerbound ElogpZπ and ElogqZ function stats{T}(vg::VGMM, x::Matrix{T}, ex::Tuple) ng = vg.n - (nx, d) = size(x) - if nx == 0 + (nₓ, d) = size(x) + if nₓ == 0 return 0, zero(RT), zeros(RT, ng), zeros(RT, ng, d), [zeros(RT, d,d) for k=1:ng] end - r, ElogpZπqZ = rnk(vg, x, ex) # nx * ng + r, ElogpZπqZ = rnk(vg, x, ex) # nₓ × ng N = vec(sum(r, 1)) # ng - F = r' * x # ng * d + F = r' * x # ng × d ## S_k = sum_i r_ik x_i x_i' ## Sm = x' * hcat([broadcast(*, r[:,k], x) for k=1:ng]...) - SS = similar(x, nx, d*ng) # nx*nd*ng mem, nx*nd*ng multiplications + SS = similar(x, nₓ, d*ng) # nₓ*nd*ng mem, nₓ*nd*ng multiplications for k=1:ng - for j=1:d for i=1:nx + for j=1:d for i=1:nₓ @inbounds SS[i,(k-1)*d+j] = r[i,k]*x[i,j] end end end - Sm = x' * SS # d * (d*ng) mem + Sm = x' * SS # d * (d * ng) mem S = Matrix{T}[Sm[:,(k-1)*d+1:k*d] for k=1:ng] - return nx, ElogpZπqZ, N, F, S + return nₓ, ElogpZπqZ, N, F, S end function stats(vg::VGMM, d::Data, ex::Tuple; parallel=false) @@ -203,7 +203,7 @@ function lowerbound(vg::VGMM, N::Vector, mx::Matrix, S::Vector, ## E[log p(x|Ζ,μ,Λ)] 10.71 Elogll = 0. for k in gaussians - Δ = mx[k,:] - m[k,:] # 1 * d + Δ = mx[k,:] - m[k,:] # 1 × d Wk = precision(W[k]) Elogll += 0.5N[k] * (ElogdetΛ[k] - d/β[k] - d*log(2π) - ν[k] * (trAB(S[k], Wk) + Δ * Wk ⋅ Δ)) # check chol efficiency @@ -212,7 +212,7 @@ function lowerbound(vg::VGMM, N::Vector, mx::Matrix, S::Vector, Elogpπ = lgamma(ng*α0) - ng*lgamma(α0) - (α0-1)sum(Elogπ) # E[log p(π)] 10.73 ElogpμΛ = ng*logB(W0,ν0) # E[log p(μ, Λ)] B.79 for k in gaussians - Δ = m[k,:] - m0' # 1 * d + Δ = m[k,:] - m0' # 1 × d Wk = precision(W[k]) ElogpμΛ += 0.5(d*log(β0/(2π)) + (ν0-d)ElogdetΛ[k] - d*β0/β[k] -β0*ν[k] * dot(Δ*Wk, Δ) - ν[k]*trAB(W0inv, Wk)) @@ -238,7 +238,7 @@ function emstep!(vg::VGMM, x::DataOrMatrix) ## E-like step Elogπ, ElogdetΛ = expectations(vg) ## N, mx, S, r = threestats(vg, x, (Elogπ, ElogdetΛ)) - nx, ElogZπqZ, N, F, S = stats(vg, x, (Elogπ, ElogdetΛ)) + nₓ, ElogZπqZ, N, F, S = stats(vg, x, (Elogπ, ElogdetΛ)) mx = F ./ N for k = 1:vg.n mk = mx[k,:] @@ -257,24 +257,24 @@ function emstep!(vg::VGMM, x::DataOrMatrix) L = lowerbound(vg, N, mx, S, Elogπ, ElogdetΛ, ElogZπqZ) ## and finally compute M-like step vg.α, vg.β, vg.m, vg.ν, vg.W = mstep(vg.π, N, mx, S) - L, nx + L, nₓ end ## This is called em!, but it is not really expectation maximization I think function em!(vg::VGMM, x; nIter=50) L = Float64[] - nx = 0 + nₓ = 0 for i=1:nIter - lb, nx = emstep!(vg, x) + lb, nₓ = emstep!(vg, x) push!(L, lb) if i>1 && isapprox(L[i], L[i-1], rtol=0) nIter=i break end - addhist!(vg, @sprintf("iteration %d, lowerbound %f", i, last(L)/nx/vg.d)) + addhist!(vg, @sprintf("iteration %d, lowerbound %f", i, last(L)/nₓ/vg.d)) end - L ./= nx*vg.d - addhist!(vg, @sprintf("%d variational Bayes EM-like iterations using %d data points, final lowerbound %f", nIter, nx, last(L))) + L ./= nₓ*vg.d + addhist!(vg, @sprintf("%d variational Bayes EM-like iterations using %d data points, final lowerbound %f", nIter, nₓ, last(L))) L end diff --git a/src/data.jl b/src/data.jl index 071cf3d..77cebd7 100644 --- a/src/data.jl +++ b/src/data.jl @@ -73,15 +73,15 @@ Base.done(x::Data, state::Int) = state == length(x) ## worker, so that file caching at the local machine is beneficial function dmap(f::Function, x::Data) if kind(x) == :file - nx = length(x) + nₓ = length(x) w = workers() nw = length(w) worker(i) = w[1 .+ ((i-1) % nw)] - results = cell(nx) + results = cell(nₓ) getnext(i) = x.list[i] load = x.API[:load] @sync begin - for i = 1:nx + for i = 1:nₓ @async begin next = getnext(i) results[i] = remotecall_fetch(worker(i), s->f(load(s)), next) diff --git a/src/rand.jl b/src/rand.jl index cdcb664..a7588a0 100644 --- a/src/rand.jl +++ b/src/rand.jl @@ -52,11 +52,11 @@ function Base.rand(gmm::GMM, n::Int) gmmkind = kind(gmm) for i=1:gmm.n ind = find(index.==i) - nx = length(ind) + nₓ = length(ind) if gmmkind == :diag - x[ind,:] = gmm.μ[i,:] .+ √gmm.Σ[i,:] .* randn(nx,gmm.d) + x[ind,:] = gmm.μ[i,:] .+ √gmm.Σ[i,:] .* randn(nₓ,gmm.d) elseif gmmkind == :full - x[ind,:] = rand(MvNormal(vec(gmm.μ[i,:]), covar(gmm.Σ[i])), nx)' + x[ind,:] = rand(MvNormal(vec(gmm.μ[i,:]), covar(gmm.Σ[i])), nₓ)' else error("Unknown kind") end diff --git a/src/recognizer.jl b/src/recognizer.jl index 39ba7a7..a53a72c 100644 --- a/src/recognizer.jl +++ b/src/recognizer.jl @@ -15,7 +15,7 @@ import Base.map ## Maximum A Posteriori adapt a gmm function map{T<:FloatingPoint}(gmm::GMM, x::Matrix{T}, r::Real=16.; means::Bool=true, weights::Bool=false, covars::Bool=false) - nx, ll, N, F, S = stats(gmm, x) + nₓ, ll, N, F, S = stats(gmm, x) α = N ./ (N+r) if weights w = α .* N / sum(N) + (1-α) .* gmm.w @@ -35,13 +35,13 @@ function map{T<:FloatingPoint}(gmm::GMM, x::Matrix{T}, r::Real=16.; means::Bool= Σ = gmm.Σ end hist = vcat(gmm.hist, History(@sprintf "MAP adapted with %d data points relevance %3.1f %s %s %s" size(x,1) r means ? "means" : "" weights ? "weights" : "" covars ? "covars" : "")) - return GMM(w, μ, Σ, hist, nx) + return GMM(w, μ, Σ, hist, nₓ) end ## compute a supervector from a MAP adapted utterance. function Base.vec(gmm::GMM, x::Matrix, r=16.) kind(gmm) == :diag || error("Sorry, can't compute MAP adapted supervector for full covariance matrix GMM yet") - nx, ll, N, F, S = stats(gmm, x) + nₓ, ll, N, F, S = stats(gmm, x) α = N ./ (N+r) Δμ = broadcast(*, α./N, F) - broadcast(*, α, gmm.μ) v = Δμ .* sqrt(broadcast(/, gmm.w, gmm.Σ)) diff --git a/src/stats.jl b/src/stats.jl index 1219990..96bbe85 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -7,7 +7,7 @@ function setmem(m::Float64) global mem=m end -## stats(gmm, x, order) computes zero, first, ... upto order (<=2) statistics of +## stats(gmm, x, order) computes zero, first, ... upto order (≤2) statistics of ## a feature file aligned to the gmm. The statistics are ordered (ng ## * d), as by the general rule for dimension order in types.jl. ## Note: these are _uncentered_ statistics. @@ -36,39 +36,39 @@ end ## diagonal covariance function stats{GT,T<:FloatingPoint}(gmm::GMM{GT,DiagCov{GT}}, x::Matrix{T}, order::Int) RT = promote_type(GT,T) - (nx, d) = size(x) + (nₓ, d) = size(x) ng = gmm.n gmm.d == d || error("dimension mismatch for data") - 0 <= order <= 2 || error("order out of range") - prec::Matrix{RT} = 1./gmm.Σ # ng * d - mp::Matrix{RT} = gmm.μ .* prec # mean*precision, ng * d + 1 ≤ order ≤ 2 || error("order out of range") + prec::Matrix{RT} = 1./gmm.Σ # ng × d + mp::Matrix{RT} = gmm.μ .* prec # mean*precision, ng × d ## note that we add exp(-sm2p/2) later to pxx for numerical stability - a::Matrix{RT} = gmm.w ./ (((2π)^(d/2)) * √ prod(gmm.Σ,2)) # ng * 1 - sm2p::Matrix{RT} = dot(mp, gmm.μ, 2) # sum over d mean^2*precision, ng * 1 - xx = x .* x # nx * d -## γ = broadcast(*, a', exp(x * mp' .- 0.5xx * prec')) # nx * ng, Likelihood per frame per Gaussian - γ = x * mp' # nx * ng, nx * d * ng multiplications + a::Matrix{RT} = gmm.w ./ ((2π)^(d/2) * √ prod(gmm.Σ,2)) # ng × 1 + sm2p::Matrix{RT} = dot(mp, gmm.μ, 2) # sum over d mean^2*precision, ng × 1 + xx = x .* x # nₓ × d +## γ = broadcast(*, a', exp(x * mp' .- 0.5xx * prec')) # nₓ × ng, Likelihood per frame per Gaussian + γ = x * mp' # nₓ × ng, nₓ * d * ng multiplications Base.BLAS.gemm!('N', 'T', -one(RT)/2, xx, prec, one(RT), γ) for j = 1:ng la = log(a[j]) - 0.5sm2p[j] - for i = 1:nx + for i = 1:nₓ @inbounds γ[i,j] += la end end for i = 1:length(γ) @inbounds γ[i] = exp(γ[i]) end - lpf=sum(γ,2) # nx * 1, Likelihood per frame - broadcast!(/, γ, γ, lpf .+ (lpf .== 0)) # nx * ng, posterior per frame per gaussian + lpf=sum(γ,2) # nₓ × 1, Likelihood per frame + broadcast!(/, γ, γ, lpf .+ (lpf .== 0)) # nₓ × ng, posterior per frame per gaussian ## zeroth order - N = vec(sum(γ, 1)) # ng * 1, vec() + N = vec(sum(γ, 1)) # ng, vec() ## first order - F = γ' * x # ng * d, Julia has efficient a' * b + F = γ' * x # ng × d, Julia has efficient a' * b llh = sum(log(lpf)) # total log likeliood if order==1 - return (nx, llh, N, F) + return (nₓ, llh, N, F) else ## second order - S = γ' * xx # ng * d - return (nx, llh, N, F, S) + S = γ' * xx # ng × d + return (nₓ, llh, N, F, S) end end @@ -76,30 +76,30 @@ end ## this is a `slow' implementation, based on posterior() function stats{GT,T<:FloatingPoint}(gmm::GMM{GT,FullCov{GT}}, x::Array{T,2}, order::Int) RT = promote_type(GT,T) - (nx, d) = size(x) + (nₓ, d) = size(x) ng = gmm.n gmm.d == d || error("dimension mismatch for data") - 0 <= order <= 2 || error("order out of range") - γ, ll = posterior(gmm, x) # nx * ng, both + 1 ≤ order ≤ 2 || error("order out of range") + γ, ll = posterior(gmm, x) # nₓ × ng, both llh = sum(logsumexp(ll .+ log(gmm.w)', 2)) ## zeroth order N = vec(sum(γ, 1)) ## first order F = γ' * x if order == 1 - return nx, llh, N, F + return nₓ, llh, N, F end ## S_k = Σ_i γ _ik x_i' * x S = Matrix{RT}[] - γx = Array(RT, nx, d) + γx = Array(RT, nₓ, d) @inbounds for k=1:ng - #broadcast!(*, γx, γ[:,k], x) # nx * d mults - for j = 1:d for i=1:nx + #broadcast!(*, γx, γ[:,k], x) # nₓ × d mults + for j = 1:d for i=1:nₓ γx[i,j] = γ[i,k]*x[i,j] end end - push!(S, x' * γx) # nx * d^2 mults + push!(S, x' * γx) # nₓ * d^2 mults end - return nx, llh, N, F, S + return nₓ, llh, N, F, S end @@ -120,17 +120,17 @@ end function stats{T<:FloatingPoint}(gmm::GMM, x::Matrix{T}; order::Int=2, parallel=false) parallel &= nworkers() > 1 ng = gmm.n - (nx, d) = size(x) + (nₓ, d) = size(x) if kind(gmm) == :diag - bytes = sizeof(T) * ((2d +2)ng + (d + ng + 1)nx) + bytes = sizeof(T) * ((2d +2)ng + (d + ng + 1)nₓ) elseif kind(gmm) == :full - bytes = sizeof(T) * ((d + d^2 + 5nx + nx*d)ng + (2d + 2)nx) + bytes = sizeof(T) * ((d + d^2 + 5nₓ + nₓ*d)ng + (2d + 2)nₓ) end blocks = iceil(bytes / (mem * (1<<30))) if parallel - blocks= min(nx, max(blocks, nworkers())) + blocks= min(nₓ, max(blocks, nworkers())) end - l = nx / blocks # chop array into smaller pieces xx + l = nₓ / blocks # chop array into smaller pieces xx xx = Matrix{T}[x[round(i*l+1):round((i+1)l),:] for i=0:(blocks-1)] if parallel r = pmap(x->stats(gmm, x, order), xx) @@ -165,9 +165,9 @@ end function csstats{T<:FloatingPoint}(gmm::GMM, x::DataOrMatrix{T}, order::Int=2) kind(gmm) == :diag || error("Can only do centered and scaled stats for diag covariance") if order==1 - nx, llh, N, F = stats(gmm, x, order) + nₓ, llh, N, F = stats(gmm, x, order) else - nx, llh, N, F, S = stats(gmm, x, order) + nₓ, llh, N, F, S = stats(gmm, x, order) end Nμ = N .* gmm.μ f = (F - Nμ) ./ gmm.Σ @@ -186,7 +186,7 @@ CSstats(gmm::GMM, x::DataOrMatrix) = CSstats(csstats(gmm, x, 1)) ## centered stats, but not scaled by UBM covariance ## check full covariance... function cstats{T<:FloatingPoint}(gmm::GMM, x::DataOrMatrix{T}, parallel=false) - nx, llh, N, F, S = stats(gmm, x, order=2, parallel=parallel) + nₓ, llh, N, F, S = stats(gmm, x, order=2, parallel=parallel) Nμ = N .* gmm.μ ## center the statistics gmmkind = kind(gmm) diff --git a/src/train.jl b/src/train.jl index 7456cdf..2c97832 100644 --- a/src/train.jl +++ b/src/train.jl @@ -40,11 +40,11 @@ GMM{T<:FloatingPoint}(n::Int, x::Vector{T}; method::Symbol=:kmeans, nInit::Int=5 ## initialize GMM using Clustering.kmeans (which uses a method similar to kmeans++) function GMMk{T}(n::Int, x::DataOrMatrix{T}; kind=:diag, nInit::Int=50, nIter::Int=10, sparse=0) - nx, d = size(x) - hist = [History(@sprintf("Initializing GMM, %d Gaussians %s covariance %d dimensions using %d data points", n, diag, d, nx))] + nₓ, d = size(x) + hist = [History(@sprintf("Initializing GMM, %d Gaussians %s covariance %d dimensions using %d data points", n, diag, d, nₓ))] ## subsample x to max 1000 points per mean nneeded = 1000*n - if nx < nneeded + if nₓ < nneeded if isa(x, Matrix) xx = x else @@ -52,7 +52,7 @@ function GMMk{T}(n::Int, x::DataOrMatrix{T}; kind=:diag, nInit::Int=50, nIter::I end else if isa(x, Matrix) - xx = x[sample(1:nx, nneeded, replace=false),:] + xx = x[sample(1:nₓ, nneeded, replace=false),:] else ## Data. Sample an equal amount from every entry in the list x. This reads in ## all data, and may require a lot of memory for very long lists. @@ -200,9 +200,9 @@ function em!(gmm::GMM, x::DataOrMatrix; nIter::Int = 10, varfloor::Float64=1e-3, gmmkind = kind(gmm) for i=1:nIter ## E-step - nx, ll[i], N, F, S = stats(gmm, x, parallel=true) + nₓ, ll[i], N, F, S = stats(gmm, x, parallel=true) ## M-step - gmm.w = N / nx + gmm.w = N / nₓ gmm.μ = F ./ N if gmmkind == :diag gmm.Σ = S ./ N - gmm.μ.^2 @@ -226,17 +226,17 @@ function em!(gmm::GMM, x::DataOrMatrix; nIter::Int = 10, varfloor::Float64=1e-3, error("Unknown kind") end addhist!(gmm, @sprintf("iteration %d, average log likelihood %f", - i, ll[i] / (nx*d))) + i, ll[i] / (nₓ*d))) end if nIter>0 - ll /= nx * d + ll /= nₓ * d finalll = ll[nIter] else finalll = avll(gmm, x) - nx = size(x,1) + nₓ = size(x,1) end - gmm.nx = nx - addhist!(gmm,@sprintf("EM with %d data points %d iterations avll %f\n%3.1f data points per parameter",nx,nIter,finalll,nx/nparams(gmm))) + gmm.nx = nₓ + addhist!(gmm,@sprintf("EM with %d data points %d iterations avll %f\n%3.1f data points per parameter",nₓ,nIter,finalll,nₓ/nparams(gmm))) ll end @@ -248,17 +248,17 @@ end function llpg{GT,T<:FloatingPoint}(gmm::GMM{GT,DiagCov{GT}}, x::Matrix{T}) RT = promote_type(GT,T) ## ng = gmm.n - (nx, d) = size(x) - prec::Matrix{RT} = 1./gmm.Σ # ng * d - mp = gmm.μ .* prec # mean*precision, ng * d + (nₓ, d) = size(x) + prec::Matrix{RT} = 1./gmm.Σ # ng × d + mp = gmm.μ .* prec # mean*precision, ng × d ## note that we add exp(-sm2p/2) later to pxx for numerical stability - normalization = 0.5 * (d * log(2π) .+ sum(log(gmm.Σ),2)) # ng * 1 - sm2p = sum(mp .* gmm.μ, 2) # sum over d mean^2*precision, ng * 1 + normalization = 0.5 * (d * log(2π) .+ sum(log(gmm.Σ),2)) # ng × 1 + sm2p = sum(mp .* gmm.μ, 2) # sum over d mean^2*precision, ng × 1 ## from here on data-dependent calculations - xx = x.^2 # nx * d - pxx = sm2p' .+ xx * prec' # nx * ng - mpx = x * mp' # nx * ng - # L = broadcast(*, a', exp(mpx-0.5pxx)) # nx * ng, Likelihood per frame per Gaussian + xx = x.^2 # nₓ × d + pxx = sm2p' .+ xx * prec' # nₓ × ng + mpx = x * mp' # nₓ × ng + # L = broadcast(*, a', exp(mpx-0.5pxx)) # nₓ × ng, Likelihood per frame per Gaussian mpx-0.5pxx .- normalization' end @@ -266,25 +266,25 @@ end ## compute Δ_i = (x_i - μ)' Λ (x_i - μ) ## Note: the return type of Δ should be the promote_type of x and μ/ciΣ function xμTΛxμ!(Δ::Matrix, x::Matrix, μ::Matrix, ciΣ::Triangular) - # broadcast!(-, Δ, x, μ) # size: nx * d, add ops: nx * d - (nx, d) = size(x) + # broadcast!(-, Δ, x, μ) # size: nₓ × d, add ops: nₓ * d + (nₓ, d) = size(x) @inbounds for j = 1:d μj = μ[j] - for i = 1:nx + for i = 1:nₓ Δ[i,j] = x[i,j] - μj end end - A_mul_Bc!(Δ, ciΣ) # size: nx * d, mult ops nx*d^2 + A_mul_Bc!(Δ, ciΣ) # size: nₓ × d, mult ops nₓ*d^2 end ## full covariance version of llpg() function llpg{GT,T<:FloatingPoint}(gmm::GMM{GT,FullCov{GT}}, x::Matrix{T}) RT = promote_type(GT,T) - (nx, d) = size(x) + (nₓ, d) = size(x) ng = gmm.n d==gmm.d || error ("Inconsistent size gmm and x") - ll = Array(RT, nx, ng) - Δ = Array(RT, nx, d) + ll = Array(RT, nₓ, ng) + Δ = Array(RT, nₓ, d) ## Σ's now are inverse choleski's, so logdet becomes -2sum(log(diag)) normalization = [0.5d*log(2π) - sum(log(diag((gmm.Σ[k])))) for k=1:ng] for k=1:ng @@ -311,9 +311,9 @@ import Distributions.posterior ## this function returns the posterior for component j: p_ij = p(j | gmm, x_i) ## TODO: This is a slow and memory-intensive implementation. It is better to ## use the approaches used in stats() -function posterior{GT,T<:FloatingPoint}(gmm::GMM{GT}, x::Matrix{T}) # nx * ng +function posterior{GT,T<:FloatingPoint}(gmm::GMM{GT}, x::Matrix{T}) # nₓ × ng RT = promote_type(GT,T) - (nx, d) = size(x) + (nₓ, d) = size(x) ng = gmm.n d==gmm.d || error("Inconsistent size gmm and x") ll = llpg(gmm, x)