correction of update algorithm

src/CumulantsUpdates.jl

@@ -11,4 +11,5 @@ module CumulantsUpdates
   include("operations.jl")
 
   export momentupdat, cumulantsupdat, vecnorm, cumnorms, cumupdatnorms
+  export momentarray, updat, moms2cums!
 end

src/updates.jl

@@ -77,3 +77,10 @@ function cumulantsupdat(cum::Vector{SymmetricTensor{T}}, X::Matrix{T}, Xup::Matr
   moms2cums!(Mup)
   Mup
 end
+
+
+function updat(M::Vector{SymmetricTensor{T}}, X::Matrix{T}, Xup::Matrix{T}) where T <: AbstractFloat
+  @inbounds Mup = [momentupdat(M[i], X, Xup) for i in 1:length(M)]
+  moms2cums!(Mup)
+  Mup
+end

test/comptimes.jl

@@ -7,50 +7,78 @@ using SymmetricTensors
 import CumulantsUpdates: cumulants, moment
 
 """
-  compspeedups(fcalc::Function, fup::Function, m::Int, n::Int, t::Int, tup::Vector{Int})
+  momspeedups(fcalc::Function, fup::Function, m::Int, n::Int, t::Int, tup::Vector{Int})
 
-Returns Vector, a computional speedup of m'th statics update of multivariate data
+Returns Vector, a computional speedup of m'th moment update of n[i] variate data
 """
 
-function compspeedups(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, f::Function,
-                                                                        fup::Function,
-                                                                        b::Int)
+function momspeedups(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, b::Int)
   compt = zeros(length(tup), length(n))
   updt = copy(compt)
   X = randn(15, 10)
-  M = f(X[1:10,:], m, b)
-  fup(M, X[1:10,:], X[10:15,:])
+  M = moment(X[1:10,:], m, b)
+  momentupdat(M, X[1:10,:], X[10:15,:])
   for i in 1:length(n)
-    println("n = ", n[i])
-    println("compute")
+    println("moment calc n = ", n[i])
     X = randn(t, n[i])
     t1 = Float64(time_ns())
-    M = f(X, m, b)
-    calctime = Float64(time_ns())-t1
+    M = moment(X, m, b)
+    compt[:,i] = Float64(time_ns())-t1
     for j in 1:length(tup)
-      println("update ", fup)
-      println("tup = ", tup[j])
+      println("update tup = ", tup[j])
       Xup = rand(tup[j], n[i])
       t2 = Float64(time_ns())
-      cup = fup(M, X, Xup)
+      cup = momentupdat(M, X, Xup)
       updt[j,i] = Float64(time_ns()) - t2
-      compt[j,i] = calctime
     end
   end
   compt, updt
 end
 
+function cumspeedups(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, b::Int)
+  compt = zeros(length(tup), length(n))
+  ccomp = copy(compt)
+  updt = copy(compt)
+  X = randn(15, 10)
+  cumulants(X[1:10,:], m, b)
+  M = momentarray(X[1:10,:], m, b)
+  M1 = copy(M)
+  moms2cums!(M1)
+  updat(M, X[1:10,:], X[10:15,:])
+  for i in 1:length(n)
+    println("cumulants calc n = ", n[i])
+    X = randn(t, n[i])
+    t1 = Float64(time_ns())
+    cumulants(X, m, b)
+    t2 = Float64(time_ns())
+    M = momentarray(X, m, b)
+    M1 = copy(M)
+    moms2cums!(M1)
+    ccomp[:,i] = Float64(time_ns())-t2
+    compt[:,i] = t2-t1
+    for j in 1:length(tup)
+      println("update tup = ", tup[j])
+      Xup = rand(tup[j], n[i])
+      t2 = Float64(time_ns())
+      cup = updat(M, X, Xup)
+      updt[j,i] = Float64(time_ns()) - t2
+    end
+  end
+  compt, updt, ccomp
+end
+
 function savecomptime(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, b::Int, p::Int)
   filename = replace("res/$(m)_$(t)_$(n)_$(tup)_$(p).jld", "[", "")
   filename = replace(filename, "]", "")
   filename = replace(filename, " ", "")
   compt = Dict{String, Any}()
-  momtime = compspeedups(m, n, t, tup, moment, momentupdat, b)
-  cumtime = compspeedups(m, n, t, tup, cumulants, cumulantsupdat, b)
+  momtime = momspeedups(m, n, t, tup, b)
+  cumtime = cumspeedups(m, n, t, tup, b)
   push!(compt, "moment" => momtime[1])
   push!(compt, "moment update" => momtime[2])
   push!(compt, "cumulants" => cumtime[1])
   push!(compt, "cumulants updat" => cumtime[2])
+  push!(compt, "cumulants bench" => cumtime[3])
   push!(compt, "tm" => t./(2*tup))
   push!(compt, "t" => t)
   push!(compt, "n" => n)
@@ -79,7 +107,7 @@ function main(args)
         arg_type = Int
       "--nvar", "-n"
         nargs = '*'
-        default = [15, 20, 25, 30]
+        default = [20, 25]
         help = "n, numbers of marginal variables"
         arg_type = Int
       "--dats", "-t"

test/comptimesblocks.jl

@@ -9,7 +9,7 @@ using ArgParse
 
 function comptime(c, data::Matrix{Float64}, datup::Matrix{Float64})
   t = time_ns()
-  cumulantsupdat(c, data, datup)
+  updat(c, data, datup)
   Float64(time_ns()-t)/1.0e9
 end
 
@@ -19,9 +19,9 @@ function savect(tup::Vector{Int}, n::Int, m::Int, p::Int)
   comptimes = zeros(maxb-1, length(tup))
   println("max block size = ", maxb)
   data = randn(maximum(tup)+10, n)
-  cumulantsupdat(cumulants(data[1:10, 1:10], 4), data[1:10, 1:10], data[1:5, 1:10])
+  updat(momentarray(data[1:10, 1:10], 4), data[1:10, 1:10], data[1:5, 1:10])
   for b in 2:maxb
-    c = cumulants(data, m, b)
+    c = momentarray(data, m, b)
     for k in 1:length(tup)
       datup = randn(tup[k], n)
       comptimes[b-1, k] = comptime(c, data, datup)

test/runtests.jl

@@ -3,7 +3,7 @@ using Cumulants
 using SymmetricTensors
 using CumulantsUpdates
 using JLD
-import CumulantsUpdates: rep, momentarray, moms2cums!, cums2moms, cnorms, normperelement
+import CumulantsUpdates: rep, momentarray, moms2cums!, cums2moms, cnorms, normperelement, updat
 
 include("updatetest.jl")
-include("operationstest.jl")
+#include("operationstest.jl")

test/updatetest.jl

@@ -85,6 +85,8 @@ end
     @test convert(Array, cup1[2]) ≈ convert(Array, cprim[2])
     @test convert(Array, cup1[3]) ≈ convert(Array, cprim[3])
   end
+end
+#=
   @testset "test on wider data" begin
     @test convert(Array, cc[1]) ≈ convert(Array, cup[1])
     @test convert(Array, cc[2]) ≈ convert(Array, cup[2])
@@ -118,3 +120,52 @@ end
   y = 2*ones(15,4)
   @test_throws BoundsError cumulantsupdat(c1, x, y)
 end
+=#
+
+@testset "cumulants update" begin
+  m = momentarray(X, 5)
+  cc = cumulants(Xprim, 5)
+  cup = updat(m, X, Xup)
+  @testset "simple test" begin
+    x = ones(6, 2)
+    y = 2*ones(2,2)
+    m1 = momentarray(x, 3)
+    cup1 = updat(m1, x, y)
+    cprim = cumulants(vcat(x,y)[3:end,:],3)
+    @test convert(Array, cup1[1]) ≈ convert(Array, cprim[1])
+    @test convert(Array, cup1[2]) ≈ convert(Array, cprim[2])
+    @test convert(Array, cup1[3]) ≈ convert(Array, cprim[3])
+  end
+  @testset "test on wider data" begin
+    @test convert(Array, cc[1]) ≈ convert(Array, cup[1])
+    @test convert(Array, cc[2]) ≈ convert(Array, cup[2])
+    @test convert(Array, cc[3]) ≈ convert(Array, cup[3])
+    @test convert(Array, cc[4]) ≈ convert(Array, cup[4])
+    @test convert(Array, cc[5]) ≈ convert(Array, cup[5])
+  end
+  addprocs(2)
+  eval(Expr(:toplevel, :(@everywhere using CumulantsUpdates)))
+  @testset "multiprocessing cumulants update" begin
+    cupp = updat(m[1:4], X, Xup)
+    @test convert(Array, cc[1]) ≈ convert(Array, cupp[1])
+    @test convert(Array, cc[2]) ≈ convert(Array, cupp[2])
+    @test convert(Array, cc[3]) ≈ convert(Array, cupp[3])
+    @test convert(Array, cc[4]) ≈ convert(Array, cupp[4])
+  end
+end
+
+@testset "updat exceptions" begin
+  x = ones(10,4);
+  y = 2*ones(5,3);
+  c1 = momentarray(x, 4);
+  @test_throws DimensionMismatch updat(c1, x, y)
+  y = 2*ones(5,4)
+  @test_throws DimensionMismatch updat(c1, x[:, 1:3], y)
+  @test_throws MethodError updat([c1[1], c1[3], c1[4]], x, y)
+  @test_throws MethodError updat([c1[1], c1[2], c1[4], c1[4]], x, y)
+  @test_throws MethodError updat([c1[1], c1[2], c1[4], c1[3]], x, y)
+  @test_throws MethodError updat([c1[2], c1[3], c1[4]], x, y)
+  @test_throws MethodError updat([c1[2], c1[4]], x, y)
+  y = 2*ones(15,4)
+  @test_throws BoundsError updat(c1, x, y)
+end