Cosmetic changes, nx -> nₓ in comments * -> × (this is \times not x, …

…boy are we pedantic)

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
 

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)

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

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.Σ))

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,70 +36,70 @@ 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
 
 ## Full covariance
 ## 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)