Skip to content

Commit

Permalink
Merge pull request #6 from ZKSI/optimization
Browse files Browse the repository at this point in the history 
Optimization
  • Loading branch information
@lpawela
lpawela authored Jun 21, 2017
2 parents e67588a + 5918759 commit ae7321a
Show file tree
Hide file tree
Showing 8 changed files with 234 additions and 193 deletions.
25 changes: 18 additions & 7 deletions src/cumulant.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,22 @@ julia> mom_el(M, (1,1), (1,1), 2)
julia> mom_el(M, (1,1), (2,2), 2)
37.0
```
"""

blockel{T <: AbstractFloat}(X::Matrix{T}, i::Tuple, j::Tuple, b::Int) =
mean(mapreduce(k::Int -> X[:,(j[k]-1)*b+i[k]], .*, 1:length(i)))
function blockel{T <: AbstractFloat}(data::Matrix{T}, mi::Tuple, mj::Tuple, b::Int)
ret = 0.
t = size(data, 1)
for l in 1:t
temp = 1.
for k in 1:length(mi)
@inbounds ind = (mj[k]-1)*b+mi[k]
@inbounds temp *= data[l,ind]
end
ret += temp
end
ret/t
end

"""
Expand Down Expand Up @@ -93,12 +105,10 @@ is a block size. Uses multicore parallel implementation via pmap()
function momentnc{T <: AbstractFloat}(x::Matrix{T}, m::Int, b::Int = 2)
t = size(x, 1)
f(z::Matrix{T}) = moment1c(z, m, b)
#k = max(2, nprocs()-1)
k = nprocs()
k = length(workers())
r = ceil(Int, t/k)
y = [x[ind2range(i, r, t), :] for i in 1:k]
ret = pmap(f, y)
#ret = pmap(z::Matrix{T} -> moment1c(z, m, b), y)
(r*sum(ret[1:(end-1)])+(t-(k-1)*r)*ret[end])/t
end

Expand All @@ -111,7 +121,7 @@ is a block size. Calls 1 core or multicore moment function.
"""

moment{T <: AbstractFloat}(X::Matrix{T}, m::Int, b::Int=2) =
(nprocs()==1)? moment1c(X, m, b): momentnc(X, m, b)
(length(workers())>1)? momentnc(X, m, b):moment1c(X, m, b)

# ---- following code is used to caclulate cumulants in SymmetricTensor form----
"""
Expand Down Expand Up @@ -234,14 +244,15 @@ julia> outprodblocks(IndexPart(Array{Int64,1}[[1,2],[3,4]],[2,2],4,2), blocks)
8.0 16.0
```
"""

function outprodblocks{T <: AbstractFloat}(inp::IndexPart,
blocks::Vector{Array{T}})
b = size(blocks[1], 1)
block = zeros(T, fill(b, inp.nind)...)
for i = 1:(b^inp.nind)
muli = ind2sub((fill(b, inp.nind)...), i)
@inbounds block[muli...] =
mapreduce(i -> blocks[i][muli[inp.part[i]]...], *, 1:inp.npart)
mapreduce(k -> blocks[k][muli[inp.part[k]]...], *, 1:inp.npart)
end
block
end
Expand Down
72 changes: 60 additions & 12 deletions src/mom2cum.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,35 @@
"""
outer(a::Vector{Float}, b::Vector{Float})
Return Vector{Float} , vectorsed outer/kroneker product o vectors a and b
Auxiliary function for rawmoment
"""
function outer{T <: AbstractFloat}(a::Vector{T}, b::Vector{T})
sa = size(a,1)
sb = size(b,1)
R = Vector{T}(sa*sb)
m = 1
for i = 1:sb, j = 1:sa
@inbounds R[m] = a[j]*b[i]
m += 1
end
return R
end

"""
updvec!(A::Vector{Float}, B::Vector{Float})
Returns updated Vector{Float} A, by adding elementwisely Vector{Float} B
Auxiliary function for rawmoment
"""
function updvec!{T<: AbstractFloat}(A::Vector{T}, B::Vector{T})
n = size(A, 1)
for i=1:n
@inbounds A[i] += B[i]
end
return A
end

"""
rawmoment(X::Matrix{T}, m::Int = 4)
Expand All @@ -7,23 +39,22 @@ pyramid structures and blocks
Returns Array{Float, m}, the m'th moment's tensor
"""

function rawmoment{T <: AbstractFloat}(X::Matrix{T}, m::Int = 4)
function rawmoment{T <: AbstractFloat}(X::Matrix{T}, m::Int = 4)
t,n = size(X)
if m == 1
return mean(X, 1)[1,:]
else
y = zeros(T, fill(n, m)...)
y = zeros(T, n^m)
z = T[1.]
for i in 1:t
@inbounds xi = X[i, :]
z = xi
for j in 2:m
@inbounds z = reshape(kron(xi', vec(z)), fill(n, j)...)
#@inbounds z = reshape(vec(z)*xi', fill(n, j)...)
for j in 1:m
z = outer(X[i, :], z)
end
y = y + z
updvec!(y, z)
z = T[1.]
end
end
y/t
reshape(y/t, fill(n, m)...)
end

"""
Expand All @@ -32,7 +63,7 @@ end
Returns [Array{Float, 1}, ..., Array{Float, k}] noncentral moment tensors of
order 1, ..., k
"""
raw_moments_upto_k{T <: AbstractFloat}(X::Matrix{T}, k::Int = 4) =
raw_moments_upto_k{T <: AbstractFloat}(X::Matrix{T}, k::Int = 4) =
[rawmoment(X, i) for i in 1:k]

"""
Expand Down Expand Up @@ -74,6 +105,23 @@ Returns Array{Float, n} the n'th cumulant tensor
"""
function onecumulant{T <: AbstractFloat}(ind::Tuple, raw::Vector{Array{T}},
spp::Vector, sppl::Vector{Vector{Int}}, dpp::Vector{Int})
ret = zero(T)
for i in 1:length(spp)
part = spp[i]
beln = length(part)
k = sppl[i]
temp = one(T)
for r in 1:beln
temp *= raw[k[r]][ind[part[r]]...]
end
ret += dpp[beln]*temp
end
ret
end


function onecumulant11{T <: AbstractFloat}(ind::Tuple, raw::Vector{Array{T}},
spp::Vector, sppl::Vector{Vector{Int}}, dpp::Vector{Int})
ret = zero(T)
for i in 1:length(spp)
Expand All @@ -88,7 +136,7 @@ end
"""
cumulatsfrommoments(x::Matrix{Float}, k::Int)
Returns a vector of 2, .., k dims Arrays{Float} of cumulant tensors
Returns a vector of 1,2, .., k dims Arrays{Float} of cumulant tensors
"""
mom2cums{T <: AbstractFloat}(x::Matrix{T}, k::Int) =
cumulants_from_moments(raw_moments_upto_k(x, k))[2:end]
cumulants_from_moments(raw_moments_upto_k(x, k))
114 changes: 52 additions & 62 deletions src/naivecumulants.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,8 @@ julia> momel(M, (1,1,1,1))
```
"""
momel{T <: AbstractFloat}(X::Matrix{T}, ind::Tuple) =
mean(mapreduce(i::Int -> X[:,ind[i]], .*, 1:length(ind)))

momel{T <: AbstractFloat}(X::Matrix{T}, multind::Tuple) = blockel(X, multind, multind, 0)

"""
Expand Down Expand Up @@ -50,21 +49,6 @@ end

# --- uses the naive method to calculate cumulants 2 - 6

"""
cumel(X::Matrix{Float}, ind::Tuple)
Returns Float an element of cumulant's tensor at ind multiindex
```jldoctest
julia> M = [[-0.88626 0.279571];[-0.704774 0.131896]];
julia> cumel(M, (1,1,1,1))
-0.80112701050308
```
"""
cumel{T<:AbstractFloat}(X::Matrix{T}, ind::Tuple) = momel(X, ind) + mixel(X, ind)

"""
mixel(X::Matrix{T}, ind::Tuple)
Expand All @@ -82,42 +66,49 @@ mixel(M, (1,1,1,1,1,1))
```
"""

function mixel{T<:AbstractFloat}(X::Matrix{T}, ind::Tuple)
A = X[:,ind[1]]
B = X[:,ind[2]]
C = X[:,ind[3]]
D = X[:,ind[4]]
if length(ind) == 4
return -mean(A.*B)*mean(C.*D) - mean(A.*C)*mean(B.*D) - mean(A.*D)*mean(B.*C)
elseif length(ind) == 5
E = X[:,ind[5]]
a = -mean(A.*B.*C)*mean(D.*E) - mean(A.*B.*D)*mean(C.*E) - mean(A.*B.*E)*mean(D.*C)
a -= mean(D.*B.*C)*mean(A.*E) + mean(E.*B.*C)*mean(D.*A) + mean(A.*D.*C)*mean(B.*E)
a -= mean(A.*E.*C)*mean(B.*D) + mean(D.*E.*C)*mean(A.*B) + mean(D.*B.*E)*mean(A.*C)
a -= mean(A.*D.*E)*mean(C.*B)
return a
elseif length(ind) == 6
E = X[:,ind[5]]
F = X[:,ind[6]]
a1 = -mean(A.*B.*C)*mean(D.*E.*F) - mean(A.*B.*D)*mean(C.*E.*F) - mean(A.*B.*E)*mean(C.*D.*F)
a1 -= mean(A.*B.*F)*mean(C.*D.*E) + mean(A.*C.*D)*mean(B.*E.*F) + mean(A.*C.*E)*mean(B.*D.*F)
a1 -= mean(A.*C.*F)*mean(B.*D.*E) + mean(A.*D.*E)*mean(B.*C.*F) + mean(A.*D.*F)*mean(B.*C.*E)
a1 -= mean(A.*E.*F)*mean(B.*C.*D)
a2 = -mean(A.*B.*C.*D)*mean(E.*F) - mean(A.*B.*C.*E)*mean(D.*F) - mean(A.*B.*C.*F)*mean(D.*E)
a2 -= mean(A.*B.*D.*E)*mean(C.*F) + mean(A.*B.*D.*F)*mean(C.*E) + mean(A.*B.*E.*F)*mean(C.*D)
a2 -= mean(A.*B)*mean(C.*D.*E.*F) + mean(A.*C.*D.*E)*mean(B.*F) + mean(A.*C.*D.*F)*mean(B.*E)
a2 -= mean(A.*C.*E.*F)*mean(B.*D) + mean(A.*C)*mean(B.*D.*E.*F) + mean(A.*D.*E.*F)*mean(B.*C)
a2 -= mean(A.*D)*mean(B.*C.*E.*F) + mean(A.*E)*mean(B.*C.*D.*F) + mean(A.*F)*mean(B.*C.*D.*E)
a3 = -mean(A.*B)*mean(C.*D)*mean(E.*F) - mean(A.*B)*mean(C.*E)*mean(D.*F)
a3 -= mean(A.*B)*mean(C.*F)*mean(D.*E) + mean(A.*C)*mean(B.*D)*mean(E.*F)
a3 -= mean(A.*C)*mean(B.*E)*mean(D.*F) + mean(A.*C)*mean(B.*F)*mean(D.*E)
a3 -= mean(A.*D)*mean(B.*C)*mean(E.*F) + mean(A.*E)*mean(B.*C)*mean(D.*F)
a3 -= mean(A.*F)*mean(B.*C)*mean(D.*E) + mean(A.*D)*mean(B.*E)*mean(C.*F)
a3 -= mean(A.*D)*mean(B.*F)*mean(C.*E) + mean(A.*E)*mean(B.*D)*mean(C.*F)
a3 -= mean(A.*F)*mean(B.*D)*mean(C.*E) + mean(A.*E)*mean(B.*F)*mean(C.*D)
a3 -= mean(A.*F)*mean(B.*E)*mean(C.*D)
return a1+a2-2*a3
function mixel{T<:AbstractFloat}(X::Matrix{T}, i::Tuple)
a = zero(T)
if length(i) == 4
a -= momel(X, (i[1],i[2]))*momel(X, (i[3],i[4]))
a -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[4])) + momel(X, (i[1],i[4]))*momel(X, (i[2],i[3]))
elseif length(i) == 5
a -= momel(X, (i[1],i[2],i[3]))*momel(X, (i[4],i[5])) + momel(X, (i[1],i[2],i[4]))*momel(X, (i[3],i[5]))
a -= momel(X, (i[1],i[2],i[5]))*momel(X, (i[3],i[4])) + momel(X, (i[2],i[3],i[4]))*momel(X, (i[1],i[5]))
a -= momel(X, (i[2],i[3],i[5]))*momel(X, (i[1],i[4])) + momel(X, (i[1],i[3],i[4]))*momel(X, (i[2],i[5]))
a -= momel(X, (i[1],i[3],i[5]))*momel(X, (i[2],i[4])) + momel(X, (i[3],i[4],i[5]))*momel(X, (i[1],i[2]))
a -= momel(X, (i[2],i[4],i[5]))*momel(X, (i[1],i[3])) + momel(X, (i[1],i[4],i[5]))*momel(X, (i[2],i[3]))
elseif length(i) == 6
a1 = -momel(X, (i[1],i[2],i[3]))*momel(X, (i[4],i[5], i[6])) - momel(X, (i[1],i[2],i[4]))*momel(X, (i[3],i[5], i[6]))
a1 -= momel(X, (i[1],i[2],i[5]))*momel(X, (i[3],i[4], i[6])) + momel(X, (i[1],i[2],i[6]))*momel(X, (i[3],i[4], i[5]))
a1 -= momel(X, (i[1],i[3],i[4]))*momel(X, (i[2],i[5], i[6])) + momel(X, (i[1],i[3],i[5]))*momel(X, (i[2],i[4], i[6]))
a1 -= momel(X, (i[1],i[3],i[6]))*momel(X, (i[2],i[4], i[5])) + momel(X, (i[1],i[4],i[5]))*momel(X, (i[2],i[3], i[6]))
a1 -= momel(X, (i[1],i[4],i[6]))*momel(X, (i[2],i[3], i[5])) + momel(X, (i[1],i[5],i[6]))*momel(X, (i[2],i[3], i[4]))
a2 = -momel(X, (i[1],i[2],i[3],i[4]))*momel(X, (i[5], i[6])) - momel(X, (i[1],i[2],i[3],i[5]))*momel(X, (i[4],i[6]))
a2 -= momel(X, (i[1],i[2],i[3],i[6]))*momel(X, (i[4], i[5])) + momel(X, (i[1],i[2],i[4],i[5]))*momel(X, (i[3], i[6]))
a2 -= momel(X, (i[1],i[2],i[4],i[6]))*momel(X, (i[3], i[5])) + momel(X, (i[1],i[2],i[5],i[6]))*momel(X, (i[3], i[4]))
a2 -= momel(X, (i[3],i[4],i[5],i[6]))*momel(X, (i[1], i[2])) + momel(X, (i[1],i[3],i[4],i[5]))*momel(X, (i[2], i[6]))
a2 -= momel(X, (i[1],i[3],i[4],i[6]))*momel(X, (i[2], i[5])) + momel(X, (i[1],i[3],i[5],i[6]))*momel(X, (i[2], i[4]))
a2 -= momel(X, (i[2],i[4],i[5],i[6]))*momel(X, (i[1], i[3])) + momel(X, (i[1],i[4],i[5],i[6]))*momel(X, (i[2], i[3]))
a2 -= momel(X, (i[2],i[3],i[5],i[6]))*momel(X, (i[1], i[4])) + momel(X, (i[2],i[3],i[4],i[6]))*momel(X, (i[1], i[5]))
a2 -= momel(X, (i[2],i[3],i[4],i[5]))*momel(X, (i[1], i[6]))
a3 = -momel(X, (i[1],i[2]))*momel(X, (i[3],i[4]))*momel(X, (i[5], i[6]))
a3 -= momel(X, (i[1],i[2]))*momel(X, (i[3],i[5]))*momel(X, (i[4], i[6]))
a3 -= momel(X, (i[1],i[2]))*momel(X, (i[3],i[6]))*momel(X, (i[4], i[5]))
a3 -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[4]))*momel(X, (i[5], i[6]))
a3 -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[5]))*momel(X, (i[4], i[6]))
a3 -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[6]))*momel(X, (i[4], i[5]))
a3 -= momel(X, (i[1],i[4]))*momel(X, (i[2],i[3]))*momel(X, (i[5], i[6]))
a3 -= momel(X, (i[1],i[5]))*momel(X, (i[2],i[3]))*momel(X, (i[4], i[6]))
a3 -= momel(X, (i[1],i[6]))*momel(X, (i[2],i[3]))*momel(X, (i[4], i[5]))
a3 -= momel(X, (i[1],i[4]))*momel(X, (i[2],i[5]))*momel(X, (i[3], i[6]))
a3 -= momel(X, (i[1],i[4]))*momel(X, (i[2],i[6]))*momel(X, (i[3], i[5]))
a3 -= momel(X, (i[1],i[5]))*momel(X, (i[2],i[4]))*momel(X, (i[3], i[6]))
a3 -= momel(X, (i[1],i[6]))*momel(X, (i[2],i[4]))*momel(X, (i[3], i[5]))
a3 -= momel(X, (i[1],i[5]))*momel(X, (i[2],i[6]))*momel(X, (i[3], i[4]))
a3 -= momel(X, (i[1],i[6]))*momel(X, (i[2],i[5]))*momel(X, (i[3], i[4]))
a += a1+a2-2*a3
end
a
end
"""
Expand All @@ -141,19 +132,18 @@ julia> naivecumulant(M, 3)
```
"""
function naivecumulant{T<:AbstractFloat}(X::Matrix{T}, m::Int = 4)
m < 7 || throw(AssertionError("naive implementation of $m cumulant not supported"))
if m == 1
return naivemoment(X,m)
end
X = X .- mean(X, 1)
n = size(X, 2)
cumulant = zeros(T, fill(n, m)...)
if m in [2,3]
for i = 1:(n^m)
ind = ind2sub((fill(n, m)...), i)
@inbounds cumulant[ind...] = momel(X, ind)
end
elseif m in [4,5,6]
ret = naivemoment(X,m)
if m in [4,5,6]
n = size(X, 2)
for i = 1:(n^m)
ind = ind2sub((fill(n, m)...), i)
@inbounds cumulant[ind...] = cumel(X, ind)
@inbounds ret[ind...] += mixel(X, ind)
end
end
cumulant
return ret
end
30 changes: 13 additions & 17 deletions src/pyramidcumulants.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,5 @@
#--- seminaive formula uses partitions and reccurence, but does not use blocks

"""
momentel(data::Matrix{T}, multind::Tuple)
Returns number, the single element of moment's tensor.
"""
momentel{T <: AbstractFloat}(data::Matrix{T}, multind::Tuple) =
mean(mapreduce(i -> data[:,multind[i]], .*, 1:length(multind)))

"""
part(n::Int)
Expand Down Expand Up @@ -59,7 +50,7 @@ julia> pyramidmoment(M, 3)
"""
function pyramidmoment{T<:AbstractFloat}(data::Matrix{T}, m::Int)
n = size(data,2)
ret = zeros(fill(n, m)...)
ret = zeros(T, fill(n, m)...)
for ind in indices(m, n)
@inbounds temp = momel(data, ind)
for per in collect(permutations([ind...]))
Expand All @@ -81,7 +72,7 @@ function mixedel{T<:AbstractFloat}(cum::Vector{Array{T}}, mulind::Tuple,
for k = 1:length(parts)
prod = 1.
for el in parts[k]
@inbounds prod*= cum[size(el,1)-1][mulind[el]...]
@inbounds prod*= cum[size(el,1)][mulind[el]...]
end
sum += prod
end
Expand All @@ -95,7 +86,7 @@ end
Returns Array{Float, m}, the mixed array (sum of products of lower cumulants)
"""
function mixedarr{T<:AbstractFloat}(cumulants::Vector{Array{T}}, m::Int)
n = size(cumulants[1], 2)
n = size(cumulants[1], 1)
sumofprod = zeros(fill(n, m)...)
parts = part(m)
for ind in indices(m, n)
Expand Down Expand Up @@ -128,12 +119,17 @@ julia> pyramidcumulants(M, 3)[2]
0.0 0.0
```
"""
function pyramidcumulants{T<:AbstractFloat}(X::Matrix{T}, m::Int)
function pyramidcumulants{T<:AbstractFloat}(X::Matrix{T}, m::Int = 4)
cumulants = Array{T}[]
push!(cumulants, pyramidmoment(X, 1))
X = X .- mean(X, 1)
cumulants = [Base.covm(X, 0, 1, false), pyramidmoment(X, 3)]
for i in 4:m
cumulant = pyramidmoment(X, i) - mixedarr(cumulants, i)
push!(cumulants, cumulant)
for i in 2:m
if i < 4
push!(cumulants, pyramidmoment(X, i))
else
cumulant = pyramidmoment(X, i) - mixedarr(cumulants, i)
push!(cumulants, cumulant)
end
end
cumulants
end
2 changes: 1 addition & 1 deletion test/comptimeproc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ end


function main()
plot(50000, 50, 4)
plot(100000, 52, 4)
end

main()
Loading

0 comments on commit ae7321a

Please sign in to comment.