Merge pull request #10 from ZKSI/optimization

Optimization

src/cumulant.jl

@@ -76,7 +76,7 @@ function usebl(bind::Tuple, n::Int, b::Int, nbar::Int)
 end
 
 """
-    momentn1(X::Matrix}, m::Int, b::Int)
+    momentn1c(X::Matrix{Float}, m::Int, b::Int)
 
 Returns: SymmetricTensor{Float, m}, a tensor of the m'th moment of X, where b
 is a block size. Uses 1 core implementation
@@ -121,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) =
-  (length(workers())>1)? momentnc(X, m, b):moment1c(X, m, b)
+  (nworkers()>1)? momentnc(X, m, b):moment1c(X, m, b)
 
 # ---- following code is used to caclulate cumulants in SymmetricTensor form----
 """

src/mom2cum.jl

@@ -1,26 +1,25 @@
 """
-  outer(a::Vector{Float}, b::Vector{Float})
+  outer!(z::Vector{Float}, a::Vector{Float}, b::Vector{Float})
 
-Return Vector{Float} , vectorsed outer/kroneker product o vectors a and b
+Return z - 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})
+function outer!{T <: AbstractFloat}(z::Vector{T}, 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]
+    for i = 1:sb for j = 1:sa
+        @inbounds z[m] = a[j]*b[i]
         m += 1
     end
-    return R
+    end
 end
 
 """
   updvec!(A::Vector{Float}, B::Vector{Float})
 
 Returns updated Vector{Float} A, by adding elementwisely Vector{Float} B
-Auxiliary function for rawmoment 
+Auxiliary function for rawmoment
 """
 function updvec!{T<: AbstractFloat}(A::Vector{T}, B::Vector{T})
   n = size(A, 1)
@@ -44,14 +43,15 @@ function rawmoment{T <: AbstractFloat}(X::Matrix{T}, m::Int = 4)
   if m == 1
     return mean(X, 1)[1,:]
   else
+    z = [map(i -> zeros(T, n^i), 1:m)...]
     y = zeros(T, n^m)
-    z = T[1.]
     for i in 1:t
-      for j in 1:m
-        z = outer(X[i, :], z)
+      xi = X[i, :]
+      z[1] = xi
+      for j in 2:m
+        outer!(z[j], xi, z[j-1])
       end
-      updvec!(y, z)
-      z = T[1.]
+      updvec!(y, z[m])
     end
   end
   reshape(y/t, fill(n, m)...)
@@ -121,18 +121,6 @@ function onecumulant{T <: AbstractFloat}(ind::Tuple, raw::Vector{Array{T}},
 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)
-    part = spp[i]
-    beln = length(part)
-    k = sppl[i]
-    ret += dpp[beln]*mapreduce(i->raw[k[i]][ind[part[i]]...], *, 1:beln)
-  end
-  ret
-end
-
 """
   cumulatsfrommoments(x::Matrix{Float}, k::Int)
 

test/comptimeblocks.jl

@@ -0,0 +1,63 @@
+#!/usr/bin/env julia
+
+using Cumulants
+using JLD
+using ArgParse
+
+
+
+function comptime(data::Matrix{Float64}, ccalc::Function, m::Int, b::Int)
+  ccalc(data[1:4, 1:4], m, 2)
+  t = time_ns()
+  ccalc(data, m, b)
+  Float64(time_ns()-t)/1.0e9
+end
+
+
+function savect(t::Vector{Int}, n::Int, m::Int)
+  maxb = round(Int, sqrt(n))
+  comptimes = zeros(maxb-1, length(t))
+  for k in 1:length(t)
+    data = randn(t[k], n)
+    for i in 2:maxb
+      comptimes[i-1, k] = comptime(data, cumulants, m, i)
+    end
+  end
+  filename = replace("res2/$(m)_$(t)_$(n)_nblocks.jld", "[", "")
+  filename = replace(filename, "]", "")
+  compt = Dict{String, Any}("cumulants"=> comptimes)
+  push!(compt, "t" => t)
+  push!(compt, "n" => n)
+  push!(compt, "m" => m)
+  push!(compt, "x" => "block size")
+  push!(compt, "block size" => [collect(2:maxb)...])
+  push!(compt, "functions" => [["cumulants"]])
+  save(filename, compt)
+end
+
+
+function main(args)
+  s = ArgParseSettings("description")
+  @add_arg_table s begin
+      "--order", "-m"
+        help = "m, the order of cumulant, ndims of cumulant's tensor"
+        default = 4
+        arg_type = Int
+      "--nvar", "-n"
+        default = 50
+        help = "n, numbers of marginal variables"
+        arg_type = Int
+      "--dats", "-t"
+        help = "t, numbers of data records"
+        nargs = '*'
+        default = [10000, 15000]
+        arg_type = Int
+    end
+  parsed_args = parse_args(s)
+  m = parsed_args["order"]
+  n = parsed_args["nvar"]
+  t = parsed_args["dats"]
+  savect(t::Vector{Int}, n::Int, m::Int)
+end
+
+main(ARGS)

test/comptimeproc.jl

@@ -1,12 +1,10 @@
 #!/usr/bin/env julia
 
 using Cumulants
-using PyCall
-@pyimport matplotlib as mpl
-mpl.use("Agg")
-using PyPlot
-mpl.rc("text", usetex=true)
-mpl.rc("font", family="serif", size = 12)
+using JLD
+using ArgParse
+
+
 
 function comptime(data::Matrix{Float64}, ccalc::Function, m::Int, b::Int)
   ccalc(data[1:4, 1:4], m, b)
@@ -16,37 +14,64 @@ function comptime(data::Matrix{Float64}, ccalc::Function, m::Int, b::Int)
 end
 
 
-function comptimesonprocs(t::Int, n::Int, m::Int, p::Int = 12)
+function comptimesonprocs(t::Int, n::Int, m::Int, p::Int)
   data = randn(t, n)
-  time = Float64[]
-  prc = Float64[]
+  times = zeros(p)
   for i in 1:p
-    rmprocs(procs()[2:end])
-    addprocs(i-1)
-    println(nprocs())
+    addprocs(i)
+    println(nworkers())
     @everywhere using Cumulants
-    push!(time, comptime(data, moment, m, 3))
-    push!(prc, 1.*nprocs())
+    times[i] = comptime(data, moment, m, 3)
+    rmprocs(workers())
   end
-  time, prc
+  times
 end
 
 
-
-function plot(t::Int, n::Int, m::Int)
-  a, c = comptimesonprocs(t,n,m)
-  b = a[1]./a[1:end]
-  fig, ax = subplots(figsize = (4.6, 4.6))
-  ax[:plot](c, b, "--x", label= "m = $m, t = $t, n = $n")
-  ax[:set_ylabel]("speedup of computional time")
-  ax[:set_xlabel]("core numbers")
-  ax[:legend](fontsize = 12, loc = 4, ncol = 1)
-  fig[:savefig]("test$n$t.eps")
+function savect(t::Int, n::Int, m::Int, maxprocs::Int)
+  comptimes = zeros(maxprocs)
+  comptimes = comptimesonprocs(t,n,m,maxprocs)
+  onec = fill(comptimes[1], maxprocs)
+  filename = replace("res2/$(m)_$(t)_$(n)_nprocs.jld", "[", "")
+  filename = replace(filename, "]", "")
+  compt = Dict{String, Any}("cumulants1c"=> onec, "cumulantsnc"=> comptimes)
+  push!(compt, "t" => [t])
+  push!(compt, "n" => n)
+  push!(compt, "m" => m)
+  push!(compt, "x" => "procs")
+  push!(compt, "procs" => collect(1:maxprocs))
+  push!(compt, "functions" => [["cumulants1c", "cumulantsnc"]])
+  save(filename, compt)
 end
 
 
-function main()
-  plot(100000, 52, 4)
+function main(args)
+  s = ArgParseSettings("description")
+  @add_arg_table s begin
+      "--order", "-m"
+        help = "m, the order of cumulant, ndims of cumulant's tensor"
+        default = 4
+        arg_type = Int
+      "--nvar", "-n"
+        default = 40
+        help = "n, numbers of marginal variables"
+        arg_type = Int
+      "--dats", "-t"
+        help = "t, numbers of data records"
+        #nargs = '*'
+        default = 100000
+        arg_type = Int
+      "--maxprocs", "-p"
+        help = "maximal number of procs"
+        default = 4
+        arg_type = Int
+    end
+  parsed_args = parse_args(s)
+  m = parsed_args["order"]
+  n = parsed_args["nvar"]
+  t = parsed_args["dats"]
+  p = parsed_args["maxprocs"]
+  savect(t::Int, n::Int, m::Int, p::Int)
 end
 
-main()
+main(ARGS)

test/comptimes.jl

@@ -1,11 +1,8 @@
 #!/usr/bin/env julia
 
 using Cumulants
-using Distributions
 using JLD
 using ArgParse
-import SymmetricTensors: indices
-import Cumulants: rawmoment, mom2cums, moment
 
 
 """
@@ -45,17 +42,17 @@ end
 Save a file in jld format of the computional times of moment, naivemoment, rawmoment
 
 """
-function savecomptime(m::Int, t::Vector{Int}, n::Vector{Int}, cache::Bool)
+function savecomptime(m::Int, t::Vector{Int}, n::Vector{Int})
   filename = replace("res2/"*string(m)*string(t)*string(n)*".jld", "[", "_")
   filename = replace(filename, "]", "")
-  if !(isfile(filename) & cache)
-    fs = [moment, naivemoment, rawmoment, cumulants, mom2cums, naivecumulant]
-    compt = Dict{String, Any}("$f"[11:end] => comtimes(m, t, n, f) for f in fs)
-    push!(compt, "t" => t)
-    push!(compt, "n" => n)
-    push!(compt, "m" => m)
-    save(filename, compt)
-  end
+  fs = [moment, naivemoment, cumulants, naivecumulant]
+  compt = Dict{String, Any}("$f"[11:end] => comtimes(m, t, n, f) for f in fs)
+  push!(compt, "t" => t)
+  push!(compt, "n" => n)
+  push!(compt, "m" => m)
+  push!(compt, "x" => "n")
+  push!(compt, "functions" => [["naivemoment", "moment"], ["naivecumulant", "cumulants"]])
+  save(filename, compt)
 end
 
 """
@@ -81,18 +78,12 @@ function main(args)
         nargs = '*'
         default = [4000]
         arg_type = Int
-      "--cache", "-c"
-        help = "indicates if computional times should be saved in a file or read
-          from a file"
-        default = true
-        arg_type = Bool
     end
   parsed_args = parse_args(s)
   m = parsed_args["order"]
   n = parsed_args["nvar"]
   t = parsed_args["dats"]
-  cache = parsed_args["cache"]
-  savecomptime(m, t, n, cache)
+  savecomptime(m, t, n)
 end
 
 main(ARGS)

test/plotcomptimes.jl

@@ -4,39 +4,44 @@ mpl.use("Agg")
 using PyPlot
 using JLD
 
+function singleplot(filename::String, name::String, compare::String = "")
+    d = load("res2/$filename"*".jld")
+  if compare == ""
+    comptimes = d[name]
+    ylab = "computional time [s]"
+  else
+    comptimes = d[name]./d[compare]
+    ylab = "speedup"
+  end
+  x = d["x"]
+  t = d["t"]
+  m = d["m"]
+  mpl.rc("text", usetex=true)
+  mpl.rc("font", family="serif", size = 8)
+  fig, ax = subplots(figsize = (3, 2.3))
+  for i in 1:size(comptimes, 2)
+    tt = t[i]
+    ax[:plot](d[x], comptimes[:,i], "--x", label= "t = $tt")
+  end
+  PyPlot.ylabel(ylab, labelpad = -1)
+  PyPlot.xlabel(x, labelpad = -3)
+  ax[:legend](fontsize = 7, loc = 2, ncol = 1)
+  fig[:savefig]("res2/"*name*filename*".eps")
+end
+
 
 """
   pltspeedup(comptimes::Array{Float}, m::Int, n::Vector{Int}, T::Vector{Int}, label::String)
 
 Returns a figure in .eps format of the computional speedup of cumulants function
+
 """
 
 function pltspeedup(filename::String)
-  d = load("res2/$filename")
-  singleplot(d, "naivemoment", "moment")
-  singleplot(d, "rawmoment", "moment")
-  singleplot(d, "naivecumulant", "cumulants")
-  singleplot(d, "mom2cums", "cumulants")
-end
-
-function singleplot(d::Dict, name::String, compare::String)
-  comptimes = d[name]./d[compare]
-  t = d["t"]
-  n = d["n"]
-  m = d["m"]
-  mpl.rc("text", usetex=true)
-  mpl.rc("font", family="serif", size = 8)
-  fig, ax = subplots(figsize = (3, 2.3))
-  for i in 1:size(comptimes, 2)
-    tt = t[i]
-    ax[:plot](n, comptimes[:,i], "--x", label= "m = $m, t = $tt")
+  d = load("res2/$filename"*".jld")
+  for f in d["functions"]
+    singleplot(filename::String, f...)
   end
-  PyPlot.ylabel("speedup of computional time", labelpad = -1)
-  PyPlot.xlabel("m", labelpad = -3)
-  ax[:legend](fontsize = 8, loc = 2, ncol = 1)
-  name = replace("$name$m$t$n", "[", "_")
-  name = replace(name, "]", "")
-  fig[:savefig]("res2/"*name*".eps")
 end
 
-pltspeedup("4_50000_10,15,20,25,30,35.jld")
+pltspeedup("4_10000,15000_50_nblocks")