text:
using LinearAlgebra
using Rimu
using Test

@testset "Interface basics" begin
    @test eltype(StyleUnknown{String}()) == String
    @test StochasticStyle(['a', 'b']) == StyleUnknown{Char}()

    vector = [1, 2, 3]
    deposit!(vector, 1, 1, 1 => 1)
    @test vector == [2, 2, 3]
    @test storage(vector) ≡ vector
    @test localpart(vector) ≡ vector
    zero!(vector)
    @test vector == [0, 0, 0]

    ham = [1 0 0; 2 3 0; 5 6 0]
    @test offdiagonals(ham, 1) == [2 => 2, 3 => 5]
    @test offdiagonals(ham, 2) == [3 => 6]
    @test offdiagonals(ham, 3) == []

    @test num_offdiagonals(ham, 1) == 2
    @test num_offdiagonals(ham, 2) == 1
    @test num_offdiagonals(ham, 3) == 0

    @test get_offdiagonal(ham, 1, 1) == (2 => 2)
    @test get_offdiagonal(ham, 1, 2) == (3 => 5)
    @test get_offdiagonal(ham, 2, 1) == (3 => 6)

    @test starting_address(ham) == 3

    @test LOStructure(ham) == AdjointKnown()
    @test has_adjoint(ham)
end
@testset "lomc! with matrix" begin
    ham = [1 1 2 3 2;
           2 0 2 2 3;
           0 0 0 3 2;
           0 0 1 1 2;
           0 1 0 1 0]
    vector = ones(5)
    post_step = ProjectedEnergy(ham, vector)
    df, st = lomc!(ham, vector; laststep=10_000, post_step)
    eigs = eigen(ham)
    @test df.shift[end] ≈ eigs.values[1] rtol=0.01
    @test df.hproj[end] / df.vproj[end] ≈ eigs.values[1] rtol=0.01
    @test normalize(st.replicas[1].v) ≈ eigs.vectors[:, 1] rtol=0.01
end


text:
using PyPlot

# debug
function say(sth)
  println(sth)
end

XI100 = 124.342
XI80 = 101.879

function test_a(;verbose=true)
  n = 10000
  bins = 101

  U = rand(n)

  (range, counts) = hist(U, 0:1/bins:1)

  D = sum((counts - n/bins).^2/(n/bins))

  verbose ? say("n=$n, bins=$bins, D=$(round(D*1000)/1000) $(D>XI100?"> $XI100, Reject H0, bad random":"< $XI100, Fail to reject H0, good random")") : nothing

  D
end

# dist
function chi_square(x, nu)
  x.^(nu/2-1).*exp(-x/2)./(2^(nu/2)*gamma(nu/2))
end

function test_c(;verbose=true)
  n = 10000
  k0 = 9
  bins = k0^2
  ind(x) = floor(x*9)

  U = rand(n,2)

  Data = [Float64(k0 * ind(U[i,1]) + ind(U[i,2])) for i=1:n]

  (range, counts) = hist(Data[:], -1:80)

  D = sum((counts - n/bins).^2/(n/bins))

  verbose ? say("n=$n, k0=$k0, bins=$bins, D=$(round(D*1000)/1000) $(D>XI80?"> $XI80, Reject H0, bad random":"< $XI80, Fail to reject H0, good random")") : nothing

  D
end

ext: jl

max_stars_repo_name: innerlee/testRand

max_issues_repo_name: innerlee/testRand

max_forks_repo_name: innerlee/testRand


text:
module MyGraphModule

export MyGraph, MyGraphOut, addNode!, addArcNoDup!, allArcs, outArcs, inArcs, shortestPath, getArc

mutable struct MyArc
    from::Int32
    to::Int32
    cost::Float64
end

mutable struct MyGraph
    # outnodes[i] are the nodes that can be reached using outgoing arcs from node i
    # outCost[i][j] is the cost of the arc (i,outnodes[i][j])
    outArcs::Vector{Vector{Int32}}
    # inNodes[i] are the nodes that can be reached using incoming arcs to node i
    # inCost[i][j] is the cost of the arc (inNodes[i][j],i)
    inArcs::Vector{Vector{Int32}}
    arcs::Vector{MyArc}
    arcDict::Dict{Pair{Int32,Int32},Int64}
end
MyGraph() = MyGraph(Vector{Vector{Int32}}(), Vector{Vector{Int32}}(), Vector{MyArc}(), Dict{Pair{Int32,Int32},Int64}())
MyGraph(n::Int) = MyGraph([Vector{Int32}() for i = 1:n], [Vector{Int32}() for i = 1:n], Vector{MyArc}(), Dict{Pair{Int32,Int32},Int64}())

# returns the pair (|V|,|A|)
function size(G::MyGraph)
    return length(G.outArcs), length(G.arcs)
end

function getArc(g, arcId)
    return g.arcs[arcId]
end

function addNode!(G::MyGraph)
    push!(G.outArcs, Vector{Int32}())
    push!(G.inArcs, Vector{Int32}())
    # return the index of the new node
    return convert(Int32, length(G.outArcs))
end

# adds an arc. Does not check for duplication
function addArcNoDup!(G::MyGraph,from::Int32,to::Int32, cost::Float64)
    push!(G.arcs, MyArc(from,to,cost))
    push!(G.outArcs[from], length(G.arcs))
    push!(G.inArcs[to], length(G.arcs))
    G.arcDict[Pair(from,to)] = length(G.arcs)
end

function cost(G::MyGraph,from::Int32,to::Int32)
    if haskey(G.arcDict, Pair(from,to))
        return G.arcs[G.arcDict[Pair(from,to)]].cost
    else
        throw("cost(G,from,to): Arc not found")
    end
end

function allArcs(G::MyGraph)
    return G.arcs
end

mutable struct ArcIterator
    g::MyGraph
    arcIds::Vector{Int32}
end

function outArcs(g::MyGraph, i::Int32)
    if i > length(g.outArcs)
        throw("Node is out of range")
    end
    return ArcIterator(g, g.outArcs[i])
end

function inArcs(g::MyGraph, i::Int32)
    if i > length(g.inArcs)
        throw("Node is out of range")
    end
    return ArcIterator(g, g.inArcs[i])
end

function Base.iterate(ai::ArcIterator)
    if length(ai.arcIds) == 0
        return nothing
    else
        arc = ai.g.arcs[ai.arcIds[1]]
    end
    return arc,2
end

function Base.iterate(ai::ArcIterator, nextIndex)
    if nextIndex > length(ai.arcIds)
        return nothing
    else
        arc = ai.g.arcs[ai.arcIds[nextIndex]]
        return arc, nextIndex+1
    end
end

mutable struct PromisingNode
    prev::Int32
    cost::Float64
end

function getPath(node, prevNode)
    path = Vector{Int32}()
    push!(path, node)
    while (prevNode[node] != 0)
        node = prevNode[node]
        push!(path, node)
    end
    return reverse(path)
end

function shortestPath(g::MyGraph, nodeOrder::Vector{Int32}, nodeBonus::Vector{Float64}, source::Int32, sink::Int32, maxNPaths::Int64, eps::Float64)
    n,m = size(g)
    cost = typemax(Float64)*ones(Float64,n)
    prevNode = zeros(Int32,n)
    cost[source] = -nodeBonus[source]
    for i in nodeOrder
        if i != source && i!= sink
            nodeCost = typemax(Float64)
            prev = 0
            for arc in inArcs(g,i)
                if cost[arc.from] + arc.cost < nodeCost
                    nodeCost = cost[arc.from] + arc.cost
                    prev = arc.from
                end
            end
            cost[i] = nodeCost-nodeBonus[i]
            prevNode[i] = prev
            #println("$i: best cost: $(cost[i]), prev: $prev")
        end
    end
    # we've reached the sink node. Here we may want to return multiple paths with negative cost
    bestPaths = Vector{PromisingNode}()
    for arc in inArcs(g,sink)
        if cost[arc.from] + arc.cost - nodeBonus[sink] < eps
            push!(bestPaths, PromisingNode(arc.from, cost[arc.from] + arc.cost - nodeBonus[sink]))
        end
    end
    sort!(bestPaths, by = x -> x.cost)
    index = 1
    paths = Vector{Vector{Int32}}()
    costs = Vector{Float64}()
    while index <= maxNPaths && index <= length(bestPaths)
        push!(paths, getPath(bestPaths[index].prev,prevNode))
        push!(paths[index], sink)
        push!(costs, bestPaths[index].cost)
        index += 1
    end
    #println("paths returned: ", paths)
    #println("costs returned: ", costs)
    return paths, costs
end

###############################################################################
###############################################################################
###############################################################################

mutable struct MyGraphOut
    # outnodes[i] is a vector containing the nodes that can be reached using outgoing arcs from node i
    # outCost[i][j] is the cost of the arc (i,outnodes[i][j])
    outNodes::Vector{Vector{Int32}}
    outCost::Vector{Vector{Float64}}
end

#MyGraphOut() = MyGraphOut(Vector{Vector{Int32}}(), Vector{Vector{Float64}())
MyGraphOut() = MyGraphOut([],[])

# helper for iterator functions. Finds the next arc in our ordering
function nextArc(G::MyGraphOut,(fromNode,index)::Tuple{Int64,Int64})
    while (fromNode <= length(G.outNodes)) && (index > length(G.outNodes[fromNode]))
        index = 1
        fromNode += 1
    end
    return (fromNode,index)
end

# initialize the iterator
# Notice, this function has to "belong to Base", otherwise
# it won't work in a "for a in graph" statement
function Base.iterate(G::MyGraphOut)
    (fromNode,index) = nextArc(G, (1,1))
    if fromNode > length(G.outNodes)
        return nothing
    else
        arc = MyArc(fromNode, G.outNodes[fromNode][index], G.outCost[fromNode][index])
        nextState = nextArc(G, (fromNode,index+1))
        return arc, nextState
    end
end

# iteration step
# Notice, this function has to "belong to Base", otherwise
# it won't work in a "for a in graph" statement
function Base.iterate(G::MyGraphOut, (fromNode,index)::Tuple{Int64,Int64})
    if fromNode > length(G.outNodes)
        return nothing
    else
        arc = MyArc(fromNode, G.outNodes[fromNode][index], G.outCost[fromNode][index])
        nextState = nextArc(G, (fromNode,index+1))
        return arc, nextState
    end
end

function addNode!(G::MyGraphOut)
    push!(G.outNodes, Vector{Int32}())
    push!(G.outCost, Vector{Float64}())
    return convert(Int32, length(G.outNodes))
end

# adds an arc. Does not check for duplication
function addArcNoDup!(G::MyGraphOut,from::Int32,to::Int32, cost::Float64)
    push!(G.outNodes[from], to)
    push!(G.outCost[from], cost)
end

# returns the pair (|V|,|A|)
function size(G::MyGraphOut)
    nArcs = 0
    for i=1:length(G.outNodes)
        nArcs += length(G.outNodes[i])
    end
    return length(G.outNodes), nArcs
end

function getNNodes(G::MyGraphOut)
    return length(G.outNodes)
end

function createMyGraph(g::MyGraphOut)
    g2 = MyGraph()
    n = length(g.outNodes)
    g2.outArcs = [Vector{Int32}() for i=1:n]
    g2.inArcs = [Vector{Int32}() for i=1:n]
    for i=1:n
        for j = 1:length(g.outNodes[i])
            otherNode = g.outNodes[i][j]
            cost = g.outCost[i][j]
            addArcNoDup!(g2,Int32(i),otherNode, cost)
        end
    end
    return g2
end

function testMyGraphOut()
    G = MyGraphOut()
    for i=1:10
        addNode!(G)
    end
    addArcNoDup!(G,Int32(2),Int32(3),3.14)
    addArcNoDup!(G,Int32(2),Int32(5),42.42)
    addArcNoDup!(G,Int32(7),Int32(2),12.34)
    addArcNoDup!(G,Int32(9),Int32(1),12.345)
    #println(iterate(G))
    for arc in G
        println(arc)
    end

    G2 = createMyGraph(G)
    println("number of nodes and arcs in G2: ", size(G2))
    println("Arcs leaving node 2:")
    for arc in outArcs(G2,Int32(2))
        println(arc)
    end
    println("Arcs entering node 1:")
    for arc in inArcs(G2,Int32(1))
        println(arc)
    end

end

function testShortestPath()
    G = MyGraph()
    for i=1:10
        addNode!(G)
    end
    addArcNoDup!(G,Int32(1),Int32(3),3.14)
    addArcNoDup!(G,Int32(1),Int32(10),42.42)
    addArcNoDup!(G,Int32(1),Int32(2),12.34)
    addArcNoDup!(G,Int32(2),Int32(5),12.0)
    addArcNoDup!(G,Int32(3),Int32(5),12.0)
    addArcNoDup!(G,Int32(2),Int32(10),3.0)
    addArcNoDup!(G,Int32(5),Int32(10),1.0)

    nodeBonus = zeros(Float64,10)
    nodeBonus[5] = 2
    #shortestPath(g::MyGraph, nodeOrder::Vector{Int32}, nodeBonus::Vector{Float64}, source::Int32, sink::Int32, maxNPaths::Int64, eps::Float64)
    shortestPath(G, collect(Int32(1):Int32(10)), nodeBonus, Int32(1), Int32(10), 5, 17.0) #-0.00001)

end

#MyGraphModule.testMyGraphOut()
#MyGraphModule.testShortestPath()
end


text:
# This SLR(1) parser was generated by a WIP project
# This project contains a lexer, parser and a compiler in Julia
# Soon on https://github.com/Cc618

# See the syntax file 18.syntax

include("18.yy.jl")

expr = read(stdin, String)
result = parse(expr)

println(result)


text:
using ProgressMeter
using Dates

function main()
    start_t = now()

    out_dir  = "result"
    PARAMDAT = "PARAMDAT.json"
    fwrite   = "write"

    nval = 6    # rho,u,v,p, N2, N
    nch  = 2    # N2, N
    nre  = 2    # N2+N2 -> 2N+N2
                # N2+N  -> 2N+N
    R    = set_gas_const()

    xmax, ymax, nodes, vecAx, vecAy = read_allgrid()
    out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
    time_integ, nt, dt, every_outnum, in_nt, dtau, cfl,
    init_rho, init_u, init_v, init_p, init_N2, init_N, specific_heat_ratio, Rd, bdcon = input_para(PARAMDAT)

    Qbase, cellxmax, cellymax, restartnum, init_step, init_time =
        set_initQbase(xmax, ymax, restart_file, init_rho, init_u, init_v, init_p, init_N2, init_N, dt,
        specific_heat_ratio, out_file_front, out_ext, out_dir, restartnum, R, nch, nval)

    # init Delta_Qcon_hat
    volume = set_volume(nodes, cellxmax, cellymax)
    Minf   = set_Minf(bdcon, specific_heat_ratio)
    dx, dy = set_dx_lts(nodes, cellxmax, cellymax)
    reset_write(fwrite)

    # main loop
    print("threads num : ")
    println(Threads.nthreads())

    # check bd
    check_bd(bdcon)

    #throw(UndefVarError(:x))

    # main loop
    prog = Progress(nt,1)
    for t in 1:nt
        next!(prog)

        evalnum = t+restartnum
        if time_integ == "1"
            Rhat = set_gasconst_hat(Qbase,cellxmax,cellymax,nval,nch,R)

            # exlicit scheme
            Qbase    = set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rhat, nval, nch)
            Qcon     = base_to_conservative(Qbase, cellxmax, cellymax, specific_heat_ratio, nval, nch)
            Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, volume, nval, nch)

            # initial_setup
            # 境界条件分のRhatを定義
            Rhat = set_gasconst_hat(Qbase,cellxmax,cellymax,nval,nch,R)
            chmu, chlambda_tr, chD, chmolef = chemical_value(Qbase, cellxmax, cellymax, Rhat, nval, nch)

            # RHS
            # source term
            wdot = set_wdot(Qbase, cellxmax, cellymax, Rhat, nch, nre, nval)

            # advection_term
            #E_adv_hat, F_adv_hat = AUSM_plusup(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, Minf, volume, nval)
            E_adv_hat, F_adv_hat = AUSM_plus(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)

            # viscos_term
            E_vis_hat, F_vis_hat = central_diff(Qbase, Qcon, cellxmax, cellymax, chmu, chlambda_tr, chD, chmolef,
                                                vecAx, vecAy, specific_heat_ratio, volume, Rhat, nval, nch)

            RHS = setup_RHS(cellxmax, cellymax, E_adv_hat, F_adv_hat, E_vis_hat, F_vis_hat, wdot, nval, volume)

            println(wdot[50,10,:])

            # time integral
            Qcon_hat = time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, nval)

            Qcon  = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, volume, nval, nch)

            println(Qbase[50,10,:])
            Qbase = conservative_to_base(Qcon, cellxmax, cellymax, specific_heat_ratio, nval, nch)
            println(Qbase[50,10,:])

        elseif time_integ == "2"

            Rhat  = set_gasconst_hat(Qbase,cellxmax,cellymax,nval,nch,R)
            Qbase = set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rhat, nval, nch)

            # Qcn
            Qbasen    = copy(Qbase)
            Qconn     = base_to_conservative(Qbasen, cellxmax, cellymax, specific_heat_ratio, nval, nch)
            Qconn_hat = setup_Qcon_hat(Qconn, cellxmax, cellymax, volume, nval, nch)
            Qbasem    = copy(Qbase)

            for tau in 1:in_nt
                # LHS (A_adv_hat=jacobian)
                Rhat = set_gasconst_hat(Qbasem, cellxmax, cellymax, nval, nch, R)

                Qbasem   = set_boundary(Qbasem, cellxmax, cellymax, vecAx, vecAy, bdcon, Rhat, nval, nch)
                Qcon     = base_to_conservative(Qbasem, cellxmax, cellymax, specific_heat_ratio, nval, nch)
                Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, volume, nval, nch)

                # initial_setup
                # 境界条件分のRhatを定義
                Rhat = set_gasconst_hat(Qbasem, cellxmax, cellymax, nval, nch, R)
                chmu, chlambda_tr, chD, chmolef = chemical_value(Qbasem, cellxmax, cellymax, Rhat, nval, nch)
                dtau   = set_lts(Qbasem, cellxmax, cellymax, chmu, dx, dy, vecAx, vecAy, volume, specific_heat_ratio, cfl)

                # RHS
                # sorce term
                wdot, H_hat = set_wdot_for_implicit(Qbasem, cellxmax, cellymax, Rhat, nch, nre, nval, volume)
                println(wdot[50,10,:])

                # advection_term
                #E_adv_hat, F_adv_hat = AUSM_plusup(Qbasem, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, Minf, volume, nval, nch)
                E_adv_hat, F_adv_hat = AUSM_plus(Qbasem, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)

                # viscos_term
                E_vis_hat, F_vis_hat = central_diff(Qbasem, Qcon, cellxmax, cellymax, chmu, chlambda_tr, chD, chmolef,
                                                vecAx, vecAy, specific_heat_ratio, volume, Rhat, nval, nch)

                RHS = setup_RHS(cellxmax, cellymax, E_adv_hat, F_adv_hat, E_vis_hat, F_vis_hat, wdot, nval, volume)

                # lusgs_advection_term
                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_shig, B_beta_shig =
                    one_wave(Qbasem, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval, nch)
                # lusgs_viscos_term
                jalphaP, jbetaP =
                    central_diff_jacobian(Qbasem, Qcon, cellxmax, cellymax, chmu, chlambda_tr, vecAx, vecAy, specific_heat_ratio, volume)

                # LUSGS
                delta_Q = zeros(cellxmax, cellymax, nval)

                ite=0
                norm2 = zeros(nval)
                while true
                    delta_Q_temp = copy(delta_Q)
                    delta_Q      = lusgs_point(dt, dtau, Qcon_hat, Qconn_hat, delta_Q,
                                                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m,
                                                A_beta_shig, B_beta_shig, jalphaP, jbetaP,
                                                H_hat, RHS, cellxmax, cellymax, volume, nval)

                    res   = set_res(delta_Q, delta_Q_temp, cellxmax, cellymax, nval)
                    norm2 = check_converge(res, RHS, cellxmax, cellymax, init_small, nval)

                    ite += 1
                    if norm2[1] < norm_ok && norm2[2] < norm_ok && norm2[3] < norm_ok && norm2[4] < norm_ok
                        break
                    end
                    if ite % 100 ==0
                        println(" now cal norm2 ")
                        println(norm2)
                    end
                end
                output_innertime(fwrite, tau, norm2, nval)

                #println(delta_Q[50,10,:])

                for i in 2:cellxmax-1
                    for j in 2:cellymax-1
                        for l in 1:nval
                            Qcon_hat[i,j,l] = Qcon_hat[i,j,l] + delta_Q[i,j,l]
                        end
                    end
                end

                Qcon   = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, volume, nval, nch)
                Qbasem = conservative_to_base(Qcon, cellxmax, cellymax, specific_heat_ratio, nval, nch)
            end
            println(Qbase[50,10,:])
            Qbase = copy(Qbasem)
            println(Qbase[50,10,:])
        end
        #throw(UndefVarError(:x))

        output_physicaltime(fwrite, t, dt, init_step, init_time)

        check_divrege(Qbase, cellxmax, cellymax, Rhat, fwrite)
        check_small(Qbase, cellxmax, cellymax, nch)

        if round(evalnum) % every_outnum == 0
            println("\n")
            println("nt_______________________________"*string(round(evalnum)))
            Rhat = set_gasconst_hat(Qbase,cellxmax,cellymax,nval,nch,R)
            output_result(evalnum, Qbase, cellxmax, cellymax, specific_heat_ratio, out_file_front, out_ext, out_dir, Rhat, nval)
        end
    end
    end_t = now()
    output_fin(fwrite, start_t, end_t, nt, dt, in_nt, cellxmax, cellymax)
end

# -- main --
main()

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
## linalg.jl: Some generic Linear Algebra definitions

scale(X::AbstractArray, s::Number) = scale!(copy(X), s)
scale(s::Number, X::AbstractArray) = scale!(copy(X), s)

function scale{R<:Real,S<:Complex}(X::AbstractArray{R}, s::S)
    Y = Array(promote_type(R,S), size(X))
    copy!(Y, X)
    scale!(Y, s)
end

scale{R<:Real}(s::Complex, X::AbstractArray{R}) = scale(X, s)

function scale!(X::AbstractArray, s::Number)
    for i in 1:length(X)
        @inbounds X[i] *= s
    end
    X
end
scale!(s::Number, X::AbstractArray) = scale!(X, s)

cross(a::AbstractVector, b::AbstractVector) = [a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]]

triu(M::AbstractMatrix) = triu(M,0)
tril(M::AbstractMatrix) = tril(M,0)
#triu{T}(M::AbstractMatrix{T}, k::Integer)
#tril{T}(M::AbstractMatrix{T}, k::Integer)
triu!(M::AbstractMatrix) = triu!(M,0)
tril!(M::AbstractMatrix) = tril!(M,0)

#diff(a::AbstractVector)
#diff(a::AbstractMatrix, dim::Integer)
diff(a::AbstractMatrix) = diff(a, 1)
diff(a::AbstractVector) = [ a[i+1] - a[i] for i=1:length(a)-1 ]

function diff(A::AbstractMatrix, dim::Integer)
    if dim == 1
        [A[i+1,j] - A[i,j] for i=1:size(A,1)-1, j=1:size(A,2)]
    else
        [A[i,j+1] - A[i,j] for i=1:size(A,1), j=1:size(A,2)-1]
    end
end

gradient(F::AbstractVector) = gradient(F, [1:length(F)])
gradient(F::AbstractVector, h::Real) = gradient(F, [h*(1:length(F))])
#gradient(F::AbstractVector, h::AbstractVector)

diag(A::AbstractVector) = error("use diagm instead of diag to construct a diagonal matrix")
#diag(A::AbstractMatrix)

#diagm{T}(v::AbstractVecOrMat{T})

# special cases of vecnorm; note that they don't need to handle isempty(x)
function vecnormMinusInf(x)
    s = start(x)
    (v, s) = next(x, s)
    minabs = abs(v)
    while !done(x, s)
        (v, s) = next(x, s)
        minabs = Base.scalarmin(minabs, abs(v))
    end
    return float(minabs)
end
function vecnormInf(x)
    s = start(x)
    (v, s) = next(x, s)
    maxabs = abs(v)
    while !done(x, s)
        (v, s) = next(x, s)
        maxabs = Base.scalarmax(maxabs, abs(v))
    end
    return float(maxabs)
end
function vecnorm1(x)
    s = start(x)
    (v, s) = next(x, s)
    av = float(abs(v))
    T = typeof(av)
    sum::promote_type(Float64, T) = av
    while !done(x, s)
        (v, s) = next(x, s)
        sum += abs(v)
    end
    return convert(T, sum)
end
function vecnorm2(x)
    maxabs = vecnormInf(x)
    maxabs == 0 && return maxabs
    s = start(x)
    (v, s) = next(x, s)
    T = typeof(maxabs)
    scale::promote_type(Float64, T) = 1/maxabs
    y = abs(v)*scale
    sum::promote_type(Float64, T) = y*y
    while !done(x, s)
        (v, s) = next(x, s)
        y = abs(v)*scale
        sum += y*y
    end
    return convert(T, maxabs * sqrt(sum))
end
function vecnormp(x, p)
    if p > 1 || p < 0 # need to rescale to avoid overflow/underflow
        maxabs = vecnormInf(x)
        maxabs == 0 && return maxabs
        s = start(x)
        (v, s) = next(x, s)
        T = typeof(maxabs)
        spp::promote_type(Float64, T) = p
        scale::promote_type(Float64, T) = 1/maxabs
        ssum::promote_type(Float64, T) = (abs(v)*scale)^spp
        while !done(x, s)
            (v, s) = next(x, s)
            ssum += (abs(v)*scale)^spp
        end
        return convert(T, maxabs * ssum^inv(spp))
    else # 0 < p < 1, no need for rescaling (but technically not a true norm)
        s = start(x)
        (v, s) = next(x, s)
        av = float(abs(v))
        T = typeof(av)
        pp::promote_type(Float64, T) = p
        sum::promote_type(Float64, T) = av^pp
        while !done(x, s)
            (v, s) = next(x, s)
            sum += abs(v)^pp
        end
        return convert(T, sum^inv(pp))
    end
end
function vecnorm(itr, p::Real=2)
    isempty(itr) && return float(real(zero(eltype(itr))))
    p == 2 && return vecnorm2(itr)
    p == 1 && return vecnorm1(itr)
    p == Inf && return vecnormInf(itr)
    p == 0 && return convert(typeof(float(real(zero(eltype(itr))))),
                             countnz(itr))
    p == -Inf && return vecnormMinusInf(itr)
    vecnormp(itr,p)
end
vecnorm(x::Number, p::Real=2) = p == 0 ? real(x==0 ? zero(x) : one(x)) : abs(x)

norm(x::AbstractVector, p::Real=2) = vecnorm(x, p)

function norm1{T}(A::AbstractMatrix{T})
    m, n = size(A)
    Tnorm = typeof(float(real(zero(T))))
    Tsum = promote_type(Float64,Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:n
            nrmj::Tsum = 0
            for i = 1:m
                nrmj += abs(A[i,j])
            end
            nrm = max(nrm,nrmj)
        end
    end
    return convert(Tnorm, nrm)
end
function norm2{T}(A::AbstractMatrix{T})
    m,n = size(A)
    if m == 1 || n == 1 return vecnorm2(A) end
    Tnorm = typeof(float(real(zero(T))))
    (m == 0 || n == 0) ? zero(Tnorm) : convert(Tnorm, svdvals(A)[1])
end
function normInf{T}(A::AbstractMatrix{T})
    m,n = size(A)
    Tnorm = typeof(float(real(zero(T))))
    Tsum = promote_type(Float64,Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for i = 1:m
            nrmi::Tsum = 0
            for j = 1:n
                nrmi += abs(A[i,j])
            end
            nrm = max(nrm,nrmi)
        end
    end
    return convert(Tnorm, nrm)
end
function norm{T}(A::AbstractMatrix{T}, p::Real=2)
    p == 2 && return norm2(A)
    p == 1 && return norm1(A)
    p == Inf && return normInf(A)
    throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
end

norm(x::Number, p::Real=2) =
    p == 0 ? convert(typeof(real(x)), ifelse(x != 0, 1, 0)) : abs(x)

rank(A::AbstractMatrix, tol::Real) = sum(svdvals(A) .> tol)
function rank(A::AbstractMatrix)
    m,n = size(A)
    (m == 0 || n == 0) && return 0
    sv = svdvals(A)
    return sum(sv .> maximum(size(A))*eps(sv[1]))
end
rank(x::Number) = x==0 ? 0 : 1

function trace(A::AbstractMatrix)
    chksquare(A)
    sum(diag(A))
end
trace(x::Number) = x

#kron(a::AbstractVector, b::AbstractVector)
#kron{T,S}(a::AbstractMatrix{T}, b::AbstractMatrix{S})

#det(a::AbstractMatrix)

inv(a::AbstractVector) = error("argument must be a square matrix")
function inv{T}(A::AbstractMatrix{T})
    S = typeof(one(T)/one(T))
    A_ldiv_B!(convert(AbstractMatrix{S}, A), eye(S, chksquare(A)))
end

function \{TA,TB,N}(A::AbstractMatrix{TA}, B::AbstractArray{TB,N})
    TC = typeof(one(TA)/one(TB))
    A_ldiv_B!(TA == TC ? copy(A) : convert(AbstractMatrix{TC}, A), TB == TC ? copy(B) : convert(AbstractArray{TC,N}, B))
end
\(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b
/(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'
# \(A::StridedMatrix,x::Number) = inv(A)*x Should be added at some point when the old elementwise version has been deprecated long enough
# /(x::Number,A::StridedMatrix) = x*inv(A)

cond(x::Number) = x == 0 ? Inf : 1.0
cond(x::Number, p) = cond(x)

#Skeel condition numbers
condskeel(A::AbstractMatrix, p::Real=Inf) = norm(abs(inv(A))*abs(A), p)
condskeel{T<:Integer}(A::AbstractMatrix{T}, p::Real=Inf) = norm(abs(inv(float(A)))*abs(A), p)
condskeel(A::AbstractMatrix, x::AbstractVector, p::Real=Inf) = norm(abs(inv(A))*abs(A)*abs(x), p)
condskeel{T<:Integer}(A::AbstractMatrix{T}, x::AbstractVector, p::Real=Inf) = norm(abs(inv(float(A)))*abs(A)*abs(x), p)

function issym(A::AbstractMatrix)
    m, n = size(A)
    m==n || return false
    for i = 1:(n-1), j = (i+1):n
        if A[i,j] != A[j,i]
            return false
        end
    end
    return true
end

issym(x::Number) = true

function ishermitian(A::AbstractMatrix)
    m, n = size(A)
    m==n || return false
    for i = 1:n, j = i:n
        if A[i,j] != conj(A[j,i])
            return false
        end
    end
    return true
end

ishermitian(x::Number) = (x == conj(x))

function istriu(A::AbstractMatrix)
    m, n = size(A)
    for j = 1:min(n,m-1), i = j+1:m
        if A[i,j] != 0
            return false
        end
    end
    return true
end

function istril(A::AbstractMatrix)
    m, n = size(A)
    for j = 2:n, i = 1:min(j-1,m)
        if A[i,j] != 0
            return false
        end
    end
    return true
end

istriu(x::Number) = true
istril(x::Number) = true

linreg{T<:Number}(X::StridedVecOrMat{T}, y::Vector{T}) = [ones(T, size(X,1)) X] \ y

# weighted least squares
function linreg(x::AbstractVector, y::AbstractVector, w::AbstractVector)
    sw = sqrt(w)
    [sw sw.*x] \ (sw.*y)
end

# multiply by diagonal matrix as vector
#diagmm!(C::AbstractMatrix, A::AbstractMatrix, b::AbstractVector)

#diagmm!(C::AbstractMatrix, b::AbstractVector, A::AbstractMatrix)

scale!(A::AbstractMatrix, b::AbstractVector) = scale!(A,A,b)
scale!(b::AbstractVector, A::AbstractMatrix) = scale!(A,b,A)

#diagmm(A::AbstractMatrix, b::AbstractVector)
#diagmm(b::AbstractVector, A::AbstractMatrix)

#^(A::AbstractMatrix, p::Number)

#findmax(a::AbstractArray)
#findmin(a::AbstractArray)

#rref{T}(A::AbstractMatrix{T})

function peakflops(n::Integer=2000; parallel::Bool=false)
    a = rand(100,100)
    t = @elapsed a*a
    a = rand(n,n)
    t = @elapsed a*a
    parallel ? sum(pmap(peakflops, [ n for i in 1:nworkers()])) : (2*float64(n)^3/t)
end

# BLAS-like in-place y=alpha*x+y function (see also the version in blas.jl
#                                          for BlasFloat Arrays)
function axpy!(alpha, x::AbstractArray, y::AbstractArray)
    n = length(x)
    n==length(y) || throw(DimensionMismatch(""))
    for i = 1:n
        @inbounds y[i] += alpha * x[i]
    end
    y
end
function axpy!{Ti<:Integer,Tj<:Integer}(alpha, x::AbstractArray, rx::AbstractArray{Ti}, y::AbstractArray, ry::AbstractArray{Tj})
    length(x)==length(y) || throw(DimensionMismatch(""))
    if minimum(rx) < 1 || maximum(rx) > length(x) || minimum(ry) < 1 || maximum(ry) > length(y) || length(rx) != length(ry)
        throw(BoundsError())
    end
    for i = 1:length(rx)
        @inbounds y[ry[i]] += alpha * x[rx[i]]
    end
    y
end

# Elementary reflection similar to LAPACK. The reflector is not Hermitian but ensures that tridiagonalization of Hermitian matrices become real. See lawn72
function elementaryLeft!(A::AbstractMatrix, row::Integer, col::Integer)
    m, n = size(A)
    1 <= row <= m || throw(BoundsError("row cannot be less than one or larger than $(size(A,1))"))
    1 <= col <= n || throw(BoundsError("col cannot be less than one or larger than $(size(A,2))"))
    @inbounds begin
        ξ1 = A[row,col]
        normu = abs2(ξ1)
        for i = row+1:m
            normu += abs2(A[i,col])
        end
        normu = sqrt(normu)
        ν = copysign(normu,real(ξ1))
        A[row,col] += ν
        ξ1 += ν
        A[row,col] = -ν
        for i = row+1:m
            A[i,col] /= ξ1
        end
    end
    ξ1/ν
end
function elementaryRight!(A::AbstractMatrix, row::Integer, col::Integer)
    m, n = size(A)
    1 <= row <= m || throw(BoundsError("row cannot be less than one or larger than $(size(A,1))"))
    1 <= col <= n || throw(BoundsError("col cannot be less than one or larger than $(size(A,2))"))
    row <= col || error("col cannot be larger than row")
    @inbounds begin
        ξ1 = A[row,col]
        normu = abs2(ξ1)
        for i = col+1:n
            normu += abs2(A[row,i])
        end
        normu = sqrt(normu)
        ν = copysign(normu,real(ξ1))
        A[row,col] += ν
        ξ1 += ν
        A[row,col] = -ν
        for i = col+1:n
            A[row,i] /= ξ1
        end
    end
    conj(ξ1/ν)
end
function elementaryRightTrapezoid!(A::AbstractMatrix, row::Integer)
    m, n = size(A)
    1 <= row <= m || throw(BoundsError("row cannot be less than one or larger than $(size(A,1))"))
    @inbounds begin
        ξ1 = A[row,row]
        normu = abs2(A[row,row])
        for i = m+1:n
            normu += abs2(A[row,i])
        end
        normu = sqrt(normu)
        ν = copysign(normu,real(ξ1))
        A[row,row] += ν
        ξ1 += ν
        A[row,row] = -ν
        for i = m+1:n
            A[row,i] /= ξ1
        end
    end
    conj(ξ1/ν)
end

max_stars_repo_name: Maple-Wang/julia

max_issues_repo_name: Maple-Wang/julia

max_forks_repo_name: Maple-Wang/julia


text:
module BenchEquilibria

using BenchmarkTools
using LatticeBoltzmann

import LatticeBoltzmann: equilibrium, equilibrium!

suite = BenchmarkGroup()

function initialize_benchmark(q = D2Q9(), τ = 1.0, scale = 2)
    benchmark_problem = LatticeBoltzmann.TGV(q, τ, scale)

    LatticeBoltzmann.initialize(
        LatticeBoltzmann.ZeroVelocityInitialCondition(),
        q,
        benchmark_problem,
    )
end

for q in LatticeBoltzmann.Quadratures
    suite[string(q)] = BenchmarkGroup([string(q), "maxwell boltzmann equilibrium"])

    suite[string(q)]["initial equilibrium"] = @benchmarkable(
        equilibrium($q, ρ, u, T),
        setup = (ρ = 1.0;
        u = zeros(eltype($q.weights), dimension($q));
        T = 1.0)
    )

    suite[string(q)]["compute equilibrium"] = @benchmarkable(
        equilibrium!($q, ρ, u, T, f),
        setup = (ρ = 1.0;
        u = zeros(eltype($q.weights), dimension($q));
        T = 1.0;
        f = zeros(eltype($q.weights), length($q.weights)))
    )
end

end  # module
BenchEquilibria.suite


text:
ENV["MPLBACKEND"]="agg" # no GUI

using PyPlot, PyCall
using Compat
using Compat.Test

VERSION >= v"0.7.0" && using Base64

plot(1:5, 2:6, "ro-")

line = gca()[:lines][1]
@test line[:get_xdata]() == [1:5;]
@test line[:get_ydata]() == [2:6;]

fig = gcf()
@test isa(fig, PyPlot.Figure)
@test fig[:get_size_inches]() ≈ [6.4, 4.8]

s = stringmime("application/postscript", fig);
m = match(r"%%BoundingBox: *([0-9]+) +([0-9]+) +([0-9]+) +([0-9]+)", s)
@test m !== nothing
boundingbox = map(s -> parse(Int, s), m.captures)
Compat.@info("got plot bounding box ", boundingbox)
@test all([300, 200] .< boundingbox[3:4] - boundingbox[1:2] .< [450,350])

c = get_cmap("viridis")
a = 0.0:0.25:1.0

rgba = pycall(pycall(PyPlot.ScalarMappable, PyObject, cmap=c,
                     norm=PyPlot.Normalize01)["to_rgba"], PyArray, a)
@test rgba ≈ [ 0.267004  0.004874  0.329415  1.0
               0.229739  0.322361  0.545706  1.0
               0.127568  0.566949  0.550556  1.0
               0.369214  0.788888  0.382914  1.0
               0.993248  0.906157  0.143936  1.0 ]


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@info "Running functionality tests..."
@testset "Distributed Stochastic Programs" begin
    for (model, _scenarios, res, name) in problems
        tol = 1e-2
        sp = instantiate(model,
                         _scenarios,
                         optimizer = LShaped.Optimizer)
        @test_throws UnloadableStructure optimize!(sp)
        set_silent(sp)
        set_optimizer_attribute(sp, MasterOptimizer(), GLPK.Optimizer)
        set_optimizer_attribute(sp, SubProblemOptimizer(), GLPK.Optimizer)
        if name == "Infeasible" || name == "Vectorized Infeasible"
            set_optimizer_attribute(sp, FeasibilityStrategy(), FeasibilityCuts())
        end
        @testset "Distributed SP Constructs: $name" begin
            optimize!(sp, cache = true)
            @test termination_status(sp) == MOI.OPTIMAL
            @test isapprox(optimal_decision(sp), res.x̄, rtol = tol)
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp, i), res.ȳ[i], rtol = tol)
            end
            @test isapprox(objective_value(sp), res.VRP, rtol = tol)
            @test isapprox(EWS(sp), res.EWS, rtol = tol)
            @test isapprox(EVPI(sp), res.EVPI, rtol = tol)
            @test isapprox(VSS(sp), res.VSS, rtol = tol)
            @test isapprox(EV(sp), res.EV, rtol = tol)
            @test isapprox(EEV(sp), res.EEV, rtol = tol)
        end
        @testset "Distributed Sanity Check: $name" begin
            sp_nondist = copy(sp, instantiation = StageDecomposition())
            add_scenarios!(sp_nondist, scenarios(sp))
            set_optimizer(sp_nondist, LShaped.Optimizer)
            set_silent(sp_nondist)
            set_optimizer_attribute(sp_nondist, Execution(), Serial())
            set_optimizer_attribute(sp_nondist, MasterOptimizer(), GLPK.Optimizer)
            set_optimizer_attribute(sp_nondist, SubProblemOptimizer(), GLPK.Optimizer)
            if name == "Infeasible" || name == "Vectorized Infeasible"
                set_optimizer_attribute(sp_nondist, FeasibilityStrategy(), FeasibilityCuts())
            end
            optimize!(sp_nondist)
            @test termination_status(sp_nondist) == MOI.OPTIMAL
            @test scenario_type(sp) == scenario_type(sp_nondist)
            @test isapprox(stage_probability(sp), stage_probability(sp_nondist))
            @test num_scenarios(sp) == num_scenarios(sp_nondist)
            @test num_scenarios(sp) == length(scenarios(sp))
            @test num_subproblems(sp) == num_subproblems(sp_nondist)
            @test isapprox(optimal_decision(sp), optimal_decision(sp_nondist))
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp, i), optimal_recourse_decision(sp_nondist, i), rtol = sqrt(tol))
            end
            @test isapprox(objective_value(sp), objective_value(sp_nondist))
        end
        @testset "Distributed Inequalities: $name" begin
            @test EWS(sp) <= VRP(sp)
            @test VRP(sp) <= EEV(sp)
            @test VSS(sp) >= 0
            @test EVPI(sp) >= 0
            @test VSS(sp) <= EEV(sp) - EV(sp)
            @test EVPI(sp) <= EEV(sp) - EV(sp)
        end
        @testset "Distributed Copying: $name" begin
            sp_copy = copy(sp, optimizer = LShaped.Optimizer)
            set_silent(sp_copy)
            add_scenarios!(sp_copy, scenarios(sp))
            @test num_scenarios(sp_copy) == num_scenarios(sp)
            generate!(sp_copy)
            @test num_subproblems(sp_copy) == num_subproblems(sp)
            set_optimizer_attribute(sp_copy, MasterOptimizer(), () -> GLPK.Optimizer())
            set_optimizer_attribute(sp_copy, SubProblemOptimizer(), () -> GLPK.Optimizer())
            if name == "Infeasible" || name == "Vectorized Infeasible"
                set_optimizer_attribute(sp_copy, FeasibilityStrategy(), FeasibilityCuts())
            end
            optimize!(sp)
            optimize!(sp_copy)
            @test termination_status(sp_copy) == MOI.OPTIMAL
            @test isapprox(optimal_decision(sp_copy), optimal_decision(sp), rtol = tol)
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp_copy, i), optimal_recourse_decision(sp, i), rtol = sqrt(tol))
            end
            @test isapprox(objective_value(sp_copy), objective_value(sp), rtol = tol)
            @test isapprox(EWS(sp_copy), EWS(sp), rtol = tol)
            @test isapprox(EVPI(sp_copy), EVPI(sp), rtol = tol)
            @test isapprox(VSS(sp_copy), VSS(sp), rtol = tol)
            @test isapprox(EV(sp_copy), EV(sp), rtol = tol)
            @test isapprox(EEV(sp_copy), EEV(sp), rtol = tol)
        end
    end
end


text:
# Dense vector construction
function frule(
    (_, Δv),
    ::Type{<:GBVector},
    v::Vector{T}
) where {T}
    return GBVector(v), GBVector(unthunk(Δv))
end

function rrule(::Type{<:GBVector}, v::Vector{T}) where {T}
    function vecpullback(ΔΩ)
        return NoTangent(), Vector(unthunk(ΔΩ))
    end
    return GBVector(v), vecpullback
end

# Dense matrix construction
function frule(
    (_, ΔA),
    ::Type{<:GBMatrix},
    A::Matrix{T}
) where {T}
    return GBMatrix(A), GBMatrix(unthunk(ΔA))
end

function rrule(::Type{<:GBMatrix}, A::Matrix{T}) where {T}
    function vecpullback(ΔΩ)
        return NoTangent(), Matrix(unthunk(ΔΩ))
    end
    return GBMatrix(A), vecpullback
end

# Dense matrix from vector (n x 1 matrix)
function frule(
    (_, ΔA),
    ::Type{<:GBMatrix},
    A::Vector{T}
) where {T}
    return GBMatrix(A), GBMatrix(unthunk(ΔA))
end

function rrule(::Type{<:GBMatrix}, A::Vector{T}) where {T}
    sz = size(A)
    function vecpullback(ΔΩ)
        return NoTangent(), reshape(Matrix(unthunk(ΔΩ)), sz)
    end
    return GBMatrix(A), vecpullback
end

# Sparse Vector
function frule(
    (_, _, Δv),
    ::Type{<:GBVector},
    I::AbstractVector{U},
    v::Vector{T}
) where {U<:Integer, T}
    return GBVector(I, v), GBVector(I, unthunk(Δv))
end

function rrule(::Type{<:GBVector}, I::AbstractVector{U}, v::Vector{T}) where {U<:Integer, T}
    function vecpullback(ΔΩ)
        return NoTangent(), NoTangent(), nonzeros(unthunk(ΔΩ))
    end
    return GBVector(I, v), vecpullback
end

# Sparse Matrix
function frule(
    (_,_,_,Δv),
    ::Type{<:GBMatrix},
    I::AbstractVector{U},
    J::AbstractVector{U},
    v::Vector{T}
) where {U<:Integer, T}
    return GBMatrix(I, J, v), GBMatrix(I, J, Δv)
end

function rrule(
    ::Type{<:GBVector},
    I::AbstractVector{U},
    J::AbstractVector{U},
    v::Vector{T}
) where {U<:Integer, T}
    function vecpullback(ΔΩ)
        return NoTangent(), NoTangent(), NoTangent(), nonzeros(unthunk(ΔΩ))
    end
    return GBMatrix(I, J, v), vecpullback
end


text:
function _if1(n)
  if (n>3)
    return n
  else
    return 7
  end
end

function _ifelse(n)
  return ifelse(n>3, n, 7)
end

function _continue1()
  s = 0.0
  for i = 1:9
    if (i == 3)
      continue
    end
    s += 2.1
  end
  return s
end

function _break1()
  s = 0.0
  for i = 1:9
    if (i == 3)
      break
    end
    s += 2.1
  end
  return s
end

function _nestedloop1()
  s = 0.0
  for i = 1:9
      s += 2.1
      for j=1:3
        s += 7.1
      end
  end
  return s
end

function _nestedloopwithinnercontinue1()
  s = 0.0
  for i = 1:9
      s += 2.1
      for j=1:3
        if (j == 2)
          continue
        end
        s += 7.1
      end
  end
  return s
end

function _nestedloopwithinnerbreak1()
  s = 0.0
  for i = 1:9
      s += 2.1
      for j=1:3
        if (j == 2)
          break
        end
        s += 7.1
      end
  end
  return s
end

#function _bitwiseAND()
# |
#end

function _logicalAND(k)
  k == 3 && return 9
  return k
end

function _logicalOR(k)
  k == 3 || return 9
  return k
end

#testing inspiration here: https://docs.julialang.org/en/v1/manual/control-flow/
function _if_elseif_else(x, y)
  if x < y
    return 1
  elseif x > y
    return -1
  else
    return 0
  end
end

#if blocks also return a value
function _ifblock_withreturn(x,y)
  if x < y
      0
  else
      1
  end
end

function _booleanif(b)
  if (b)
    return 2
  else
    return 4
  end
end

_ternary(x,y) = x < y ? 0 : 1
_multiternary(x,y) = x < y ? 1 : x > y ? -1 : 0

function _println(n)
  if n < 3
    println("This is a message visible as console.log().")
  else
    return 1
  end
  return 2
end

function _error(n)
  if n < 3
    error("This is a message visible as console.error(), followed by uncreachable (trap) in wasm execution.")
  else
    return 1.5
  end
end

function fact(n)
  n >= 0 || error("n must be non-negative")
  n == 0 && return 1
  n * fact(n-1) #recursion
end

#=
#wasm cant have Union types
function _unstablereturntype(x)
  x > 0 && return 1
  return 2.3
end
=#

t(x) = x>0 ? true : false
f(x) = x>0 ? false : true

function shortcircuitevaluation(a,b)
  #https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation
  r1 = t(a) && f(b) #false
  r2 = f(a) && t(b) #false
  r3 = f(a) && f(b) #false
  r4 = t(a) || t(b) #true
  r5 = t(a) || f(b) #true
  r6 = f(a) || t(b) #true
  r7 = f(a) || f(b) #false

  if !r1 && !r2 && !r3 && r4 && r5 && r6 && !r7
    #expected.
    return 7
  else
    return 9
  end
end

function shortcircuitevaluation_constant(a,b)
  r1 = t(1) && f(2) #false, always, meaning r1 is actually a constant
  r2 = f(a) && t(b) #false

  if !r1 && !r2
    #expected.
    return 7
  else
    return 9
  end
end

function _unstablevariabletype()
  x = 1
  x = 1.2
  return x
end

function _unstablevariabletype2()
  x = 2 * 2.3
  return x*4.1
end

function exports()

  _if1(2)
  _ifelse(2)
  _continue1()
  _break1()

  _nestedloop1()
  _nestedloopwithinnercontinue1()
  _nestedloopwithinnerbreak1()

  _logicalAND(2)
  _logicalOR(2)

  _if_elseif_else(2,3)

  k = _ifblock_withreturn(2,3)
  k = _ternary(2,3)
  k = _multiternary(2,3)

  k = _println(5)
  k = _error(4)

  fact(5)
  #_unstablereturntype(0) #float
  #_unstablereturntype(1) #int

  #_unstablevariabletype()
  #_unstablevariabletype2()
  shortcircuitevaluation(1,2)
  shortcircuitevaluation_constant(1,2)

  _booleanif(true)

end

exports()


text:
"""
    sliding_window(seq, window_size, step_size, operation)

Perform an operation on sequential windows of a DNA sequence. Return the results of the operation in an array.
"""
function sliding_window(seq::String, window::Int64, step::Int64, op::Function)

    s::Int64 = 1
    w::Int64 = window
    e::Int64 = w
    a::Int64 = step

    results = []

    while e <= length(seq)
        subseq = seq[s:e]
        s = s + (w - (w - a))
        e = e + (w - (w - a))
        r = op(subseq)
        push!(results, r)
    end

    return results

end


text:
module OmegaViz

include("contour.jl")

end


text:
include("Parameter.jl")

export get_index
export get_name
export get_parameter
export get_units
export set_index!
export set_name!
export set_parameter!
export set_units!

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using Documenter
using DemoCards
using ImageMorphology
using OffsetArrays
using ImageBase, ImageShow, TestImages

operators, operators_cb, demo_assets = makedemos("operators")

prettyurls = get(ENV, "CI", nothing) == "true"
assets = Any[demo_assets]

pages = Any[
    "index.md",
    "Concepts" => Any["structuring_element.md"],
    operators,
    "reference.md"
]

makedocs(;
    modules=[ImageMorphology, OffsetArrays],
    format=Documenter.HTML(; prettyurls, assets),
    sitename="ImageMorphology",
    pages,
    doctest=false
)
operators_cb()

deploydocs(;
    forcepush = true,
    push_preview = true,
    repo="github.com/JuliaImages/ImageMorphology.jl.git"
)


text:
function update_r!(s::StateFS, c::ConstantsFS, d::DataFS)
  # update each r_{i, k}
  for i in 1:d.data.I
    ldmix = logdmixture_r(i, s.theta, c.constants, d.data)
    for k in 1:c.constants.K
      update_r_marg_lam!(i, k, ldmix, s, c, d)
      # NOTE: Make sure to always update W immediately after updating r or
      #       W_star!
      update_W!(s, c, d)
    end
  end
end

function compute_w(ws::Vector{Float64}, r::Vector{Bool})
  wr = ws .* r
  return wr ./ sum(wr)
end

function logdmixture_r(i::Integer, s::State, c::Constants, d::Data)
  """
  Precomputes the relevant mixture density, marginalized over gam & lam.
  This is possible because (Z, mus, sig, eta, y_miss) are held constant.
  Note that W is a function of (r, W*), and hence is NOT held constant,
  and so does not appear here.
  """
  # Precompute / reshapes
  mus0 = reshape(mus(false, s, c, d), 1, 1, c.L[0])
  mus1 = reshape(mus(true, s, c, d), 1, 1, c.L[1])
  sig = sqrt.(s.sig2 * c.temper)
  Ni = d.N[i]
  Z = reshape(s.Z, 1, d.J, c.K)
  eta0i = s.eta[false][i:i, :, :]
  eta1i = s.eta[true][i:i, :, :]

  # Ni x J x 1
  yi = reshape(s.y_imputed[i], Ni, d.J, 1)

  # Ni x J x 1
  kernel0 = MCMC.lpdf_gmm(yi, mus0, sig[i], eta0i, dims=3)
  kernel1 = MCMC.lpdf_gmm(yi, mus1, sig[i], eta1i, dims=3)

  # (1 x J x K) .* (Ni x 1 x K) -> Ni x K
  Z_mix = let
    tmp = sum(Z .* kernel1 + (1 .- Z) .* kernel0, dims=2)
    dropdims(tmp, dims=2)
  end

  return Z_mix
end

function update_r_marg_lam!(i::Integer, k::Integer, ldmix::Matrix{Float64},
                            s::StateFS, c::ConstantsFS, d::DataFS)

  # Log likelihood as a function of `r_{i, k}`
  function loglike(r_ik::Bool)::Float64
    # Make a copy of the current (k-dim) vector r_i
    ri = deepcopy(s.r[i, :])

    # Replace r_{i, k} with the provided `r_ik`
    ri[k] = r_ik

    # Compute W given W* and the updated r_i
    wi = compute_w(s.W_star[i, :], ri)

    # Get the indices of the W_i such that W_{i, k} > 0
    non_zero_wri_idx = findall(ri)

    # Initialize log likelihood
    ll = 0.0

    if length(non_zero_wri_idx) > 0  # i.e. if any W_i > 0
      wi_rs = reshape(wi[non_zero_wri_idx], 1, sum(ri))
      lln = MCMC.logsumexp(log.(wi_rs) .+ ldmix[:, non_zero_wri_idx], dims=2)
      ll = sum(lln)
    else
      ll = -Inf
    end

    return ll
  end

  # Log prior as a function of `r_{i, k}`
  function logprior(r_ik::Bool)::Float64
    p_xi = compute_p(d.X[i, :], s.omega)
    return logpdf(Bernoulli(p_xi), r_ik)
  end

  # Log full conditional as a function of `r_{i, k}`
  log_prob(r_ik::Bool)::Float64 = logprior(r_ik) + loglike(r_ik)

  # update r_{i, k} with a metropolis step
  rik_metropolis_update!(i, k, log_prob, s)
end

function rik_metropolis_update!(i::Integer, k::Integer,
                                log_prob::Function, s::StateFS)
  cand = !s.r[i, k]  # Flip bit
  log_acceptance_ratio = log_prob(cand) - log_prob(s.r[i, k])

  # accept with probability `accecptance_ratio`
  if log_acceptance_ratio > log(rand())
    s.r[i, k] = cand
  end
end


text:
using OnlineCCA
using OnlineCCA: readcsv, writecsv
using Test
using Pkg
using DelimitedFiles
using Statistics
using Distributions

tmp = mktempdir()
julia = joinpath(Sys.BINDIR, "julia")
bindir = joinpath(dirname(pathof(OnlineCCA)), "..", "bin")

function testfilesize(remove::Bool, x...)
  for n = 1:length(x)
    @test filesize(x[n]) != 0
    if remove
      rm(x[n])
    end
  end
end

#####################################
X = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 30)))
X[1:50, 1:10] .= 100*X[1:50, 1:10]
X[51:100, 11:20] .= 100*X[51:100, 11:20]
X[101:150, 21:30] .= 100*X[101:150, 21:30]
writecsv(joinpath(tmp, "X.csv"), X)

Y = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 40)))
Y[1:50, 1:15] .= 100*Y[1:50, 1:15]
Y[51:100, 16:24] .= 100*Y[51:100, 16:24]
Y[101:150, 25:40] .= 100*Y[101:150, 25:40]
writecsv(joinpath(tmp, "Y.csv"), Y)

Z = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 50)))
Z[1:50, 1:20] .= 100*Z[1:50, 1:20]
Z[51:100, 21:40] .= 100*Z[51:100, 21:40]
Z[101:150, 41:50] .= 100*Z[101:150, 41:50]
writecsv(joinpath(tmp, "Z.csv"), Z)
#####################################

#####################################
println("####### Binarization (Julia API) #######")
csv2bin(csvfile=joinpath(tmp, "X.csv"),
  binfile=joinpath(tmp, "X.zst"))
csv2bin(csvfile=joinpath(tmp, "Y.csv"),
  binfile=joinpath(tmp, "Y.zst"))
csv2bin(csvfile=joinpath(tmp, "Z.csv"),
  binfile=joinpath(tmp, "Z.zst"))

testfilesize(false, joinpath(tmp, "X.zst"))
testfilesize(false, joinpath(tmp, "Y.zst"))
testfilesize(false, joinpath(tmp, "Z.zst"))
####################################

####################################
println("####### Binarization (Command line) #######")
run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "X.csv")) --binfile $(joinpath(tmp, "X.zst"))`)
run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "Y.csv")) --binfile $(joinpath(tmp, "Y.zst"))`)
run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "Z.csv")) --binfile $(joinpath(tmp, "Z.zst"))`)

testfilesize(false, joinpath(tmp, "X.zst"))
testfilesize(false, joinpath(tmp, "Y.zst"))
testfilesize(false, joinpath(tmp, "Z.zst"))
####################################

####################################
println("####### Summarization (Julia API) #######")
mkdir(joinpath(tmp, "SumrX"))
mkdir(joinpath(tmp, "SumrY"))
mkdir(joinpath(tmp, "SumrZ"))
sumr(binfile=joinpath(tmp, "X.zst"),
  outdir=joinpath(tmp, "SumrX"))
sumr(binfile=joinpath(tmp, "Y.zst"),
  outdir=joinpath(tmp, "SumrY"))
sumr(binfile=joinpath(tmp, "Z.zst"),
  outdir=joinpath(tmp, "SumrZ"))

testfilesize(false,
  joinpath(tmp, "SumrX/Col_Means.csv"),
  joinpath(tmp, "SumrX/Col_LogMeans.csv"),
  joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrX/Col_Vars.csv"),
  joinpath(tmp, "SumrX/Col_LogVars.csv"),
  joinpath(tmp, "SumrX/Col_FTTVars.csv"),
  joinpath(tmp, "SumrX/Col_SDs.csv"),
  joinpath(tmp, "SumrX/Col_LogSDs.csv"),
  joinpath(tmp, "SumrX/Col_FTTSDs.csv"),
  joinpath(tmp, "SumrX/Col_CV2s.csv"),
  joinpath(tmp, "SumrX/Col_NoZeros.csv"),
  joinpath(tmp, "SumrX/Col_NoCounts.csv"),
  joinpath(tmp, "SumrX/Row_Means.csv"),
  joinpath(tmp, "SumrX/Row_LogMeans.csv"),
  joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrX/Row_Vars.csv"),
  joinpath(tmp, "SumrX/Row_LogVars.csv"),
  joinpath(tmp, "SumrX/Row_FTTVars.csv"),
  joinpath(tmp, "SumrX/Row_SDs.csv"),
  joinpath(tmp, "SumrX/Row_LogSDs.csv"),
  joinpath(tmp, "SumrX/Col_FTTVars.csv"),
  joinpath(tmp, "SumrX/Row_CV2s.csv"),
  joinpath(tmp, "SumrX/Row_NoZeros.csv"),
  joinpath(tmp, "SumrX/Row_NoCounts.csv"))

testfilesize(false,
  joinpath(tmp, "SumrY/Col_Means.csv"),
  joinpath(tmp, "SumrY/Col_LogMeans.csv"),
  joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrY/Col_Vars.csv"),
  joinpath(tmp, "SumrY/Col_LogVars.csv"),
  joinpath(tmp, "SumrY/Col_FTTVars.csv"),
  joinpath(tmp, "SumrY/Col_SDs.csv"),
  joinpath(tmp, "SumrY/Col_LogSDs.csv"),
  joinpath(tmp, "SumrY/Col_FTTSDs.csv"),
  joinpath(tmp, "SumrY/Col_CV2s.csv"),
  joinpath(tmp, "SumrY/Col_NoZeros.csv"),
  joinpath(tmp, "SumrY/Col_NoCounts.csv"),
  joinpath(tmp, "SumrY/Row_Means.csv"),
  joinpath(tmp, "SumrY/Row_LogMeans.csv"),
  joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrY/Row_Vars.csv"),
  joinpath(tmp, "SumrY/Row_LogVars.csv"),
  joinpath(tmp, "SumrY/Row_FTTVars.csv"),
  joinpath(tmp, "SumrY/Row_SDs.csv"),
  joinpath(tmp, "SumrY/Row_LogSDs.csv"),
  joinpath(tmp, "SumrY/Col_FTTVars.csv"),
  joinpath(tmp, "SumrY/Row_CV2s.csv"),
  joinpath(tmp, "SumrY/Row_NoZeros.csv"),
  joinpath(tmp, "SumrY/Row_NoCounts.csv"))

testfilesize(false,
  joinpath(tmp, "SumrZ/Col_Means.csv"),
  joinpath(tmp, "SumrZ/Col_LogMeans.csv"),
  joinpath(tmp, "SumrZ/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrZ/Col_Vars.csv"),
  joinpath(tmp, "SumrZ/Col_LogVars.csv"),
  joinpath(tmp, "SumrZ/Col_FTTVars.csv"),
  joinpath(tmp, "SumrZ/Col_SDs.csv"),
  joinpath(tmp, "SumrZ/Col_LogSDs.csv"),
  joinpath(tmp, "SumrZ/Col_FTTSDs.csv"),
  joinpath(tmp, "SumrZ/Col_CV2s.csv"),
  joinpath(tmp, "SumrZ/Col_NoZeros.csv"),
  joinpath(tmp, "SumrZ/Col_NoCounts.csv"),
  joinpath(tmp, "SumrZ/Row_Means.csv"),
  joinpath(tmp, "SumrZ/Row_LogMeans.csv"),
  joinpath(tmp, "SumrZ/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrZ/Row_Vars.csv"),
  joinpath(tmp, "SumrZ/Row_LogVars.csv"),
  joinpath(tmp, "SumrZ/Row_FTTVars.csv"),
  joinpath(tmp, "SumrZ/Row_SDs.csv"),
  joinpath(tmp, "SumrZ/Row_LogSDs.csv"),
  joinpath(tmp, "SumrZ/Col_FTTVars.csv"),
  joinpath(tmp, "SumrZ/Row_CV2s.csv"),
  joinpath(tmp, "SumrZ/Row_NoZeros.csv"),
  joinpath(tmp, "SumrZ/Row_NoCounts.csv"))
####################################

####################################
println("####### Summarization (Command line) #######")
run(`$(julia) $(joinpath(bindir, "sumr")) --binfile $(joinpath(tmp, "X.zst")) --outdir $(joinpath(tmp, "SumrX"))`)
run(`$(julia) $(joinpath(bindir, "sumr")) --binfile $(joinpath(tmp, "Y.zst")) --outdir $(joinpath(tmp, "SumrY"))`)
run(`$(julia) $(joinpath(bindir, "sumr")) --binfile $(joinpath(tmp, "Z.zst")) --outdir $(joinpath(tmp, "SumrZ"))`)

testfilesize(false,
  joinpath(tmp, "SumrX/Col_Means.csv"),
  joinpath(tmp, "SumrX/Col_LogMeans.csv"),
  joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrX/Col_Vars.csv"),
  joinpath(tmp, "SumrX/Col_LogVars.csv"),
  joinpath(tmp, "SumrX/Col_FTTVars.csv"),
  joinpath(tmp, "SumrX/Col_SDs.csv"),
  joinpath(tmp, "SumrX/Col_LogSDs.csv"),
  joinpath(tmp, "SumrX/Col_FTTSDs.csv"),
  joinpath(tmp, "SumrX/Col_CV2s.csv"),
  joinpath(tmp, "SumrX/Col_NoZeros.csv"),
  joinpath(tmp, "SumrX/Col_NoCounts.csv"),
  joinpath(tmp, "SumrX/Row_Means.csv"),
  joinpath(tmp, "SumrX/Row_LogMeans.csv"),
  joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrX/Row_Vars.csv"),
  joinpath(tmp, "SumrX/Row_LogVars.csv"),
  joinpath(tmp, "SumrX/Row_FTTVars.csv"),
  joinpath(tmp, "SumrX/Row_SDs.csv"),
  joinpath(tmp, "SumrX/Row_LogSDs.csv"),
  joinpath(tmp, "SumrX/Col_FTTVars.csv"),
  joinpath(tmp, "SumrX/Row_CV2s.csv"),
  joinpath(tmp, "SumrX/Row_NoZeros.csv"),
  joinpath(tmp, "SumrX/Row_NoCounts.csv"))

testfilesize(false,
  joinpath(tmp, "SumrY/Col_Means.csv"),
  joinpath(tmp, "SumrY/Col_LogMeans.csv"),
  joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrY/Col_Vars.csv"),
  joinpath(tmp, "SumrY/Col_LogVars.csv"),
  joinpath(tmp, "SumrY/Col_FTTVars.csv"),
  joinpath(tmp, "SumrY/Col_SDs.csv"),
  joinpath(tmp, "SumrY/Col_LogSDs.csv"),
  joinpath(tmp, "SumrY/Col_FTTSDs.csv"),
  joinpath(tmp, "SumrY/Col_CV2s.csv"),
  joinpath(tmp, "SumrY/Col_NoZeros.csv"),
  joinpath(tmp, "SumrY/Col_NoCounts.csv"),
  joinpath(tmp, "SumrY/Row_Means.csv"),
  joinpath(tmp, "SumrY/Row_LogMeans.csv"),
  joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrY/Row_Vars.csv"),
  joinpath(tmp, "SumrY/Row_LogVars.csv"),
  joinpath(tmp, "SumrY/Row_FTTVars.csv"),
  joinpath(tmp, "SumrY/Row_SDs.csv"),
  joinpath(tmp, "SumrY/Row_LogSDs.csv"),
  joinpath(tmp, "SumrY/Col_FTTVars.csv"),
  joinpath(tmp, "SumrY/Row_CV2s.csv"),
  joinpath(tmp, "SumrY/Row_NoZeros.csv"),
  joinpath(tmp, "SumrY/Row_NoCounts.csv"))

testfilesize(false,
  joinpath(tmp, "SumrZ/Col_Means.csv"),
  joinpath(tmp, "SumrZ/Col_LogMeans.csv"),
  joinpath(tmp, "SumrZ/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrZ/Col_Vars.csv"),
  joinpath(tmp, "SumrZ/Col_LogVars.csv"),
  joinpath(tmp, "SumrZ/Col_FTTVars.csv"),
  joinpath(tmp, "SumrZ/Col_SDs.csv"),
  joinpath(tmp, "SumrZ/Col_LogSDs.csv"),
  joinpath(tmp, "SumrZ/Col_FTTSDs.csv"),
  joinpath(tmp, "SumrZ/Col_CV2s.csv"),
  joinpath(tmp, "SumrZ/Col_NoZeros.csv"),
  joinpath(tmp, "SumrZ/Col_NoCounts.csv"),
  joinpath(tmp, "SumrZ/Row_Means.csv"),
  joinpath(tmp, "SumrZ/Row_LogMeans.csv"),
  joinpath(tmp, "SumrZ/Col_FTTMeans.csv"),
  joinpath(tmp, "SumrZ/Row_Vars.csv"),
  joinpath(tmp, "SumrZ/Row_LogVars.csv"),
  joinpath(tmp, "SumrZ/Row_FTTVars.csv"),
  joinpath(tmp, "SumrZ/Row_SDs.csv"),
  joinpath(tmp, "SumrZ/Row_LogSDs.csv"),
  joinpath(tmp, "SumrZ/Col_FTTVars.csv"),
  joinpath(tmp, "SumrZ/Row_CV2s.csv"),
  joinpath(tmp, "SumrZ/Row_NoZeros.csv"),
  joinpath(tmp, "SumrZ/Row_NoCounts.csv"))
####################################

####################################
println("####### GD (Julia API) #######")
out_gd1 = gd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="robbins-monro",
  stepsize=1.0e-15,
  numepoch=10,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_gd2 = gd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="momentum",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_gd3 = gd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="nag",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_gd4 = gd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="adagrad",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

@test size(out_gd1[1][1]) == (30, 3)
@test size(out_gd1[1][2]) == (40, 3)
@test size(out_gd1[1][3]) == (50, 3)
@test size(out_gd1[2][1]) == (3, )
@test size(out_gd1[2][2]) == (3, )
@test size(out_gd1[2][3]) == (3, )
@test size(out_gd1[3][1]) == (300, 3)
@test size(out_gd1[3][2]) == (300, 3)
@test size(out_gd1[3][3]) == (300, 3)
@test size(out_gd1[4][1]) == ()
@test size(out_gd1[4][2]) == ()
@test size(out_gd1[4][3]) == ()
@test size(out_gd1[5][1]) == ()
@test size(out_gd1[5][2]) == ()
@test size(out_gd1[5][3]) == ()
@test size(out_gd1[6][1]) == ()
@test size(out_gd1[6][2]) == ()
@test size(out_gd1[6][3]) == ()
@test size(out_gd1[7][1]) == ()
@test size(out_gd1[7][2]) == ()
@test size(out_gd1[7][3]) == ()
@test size(out_gd1[8]) == ()

@test size(out_gd2[1][1]) == (30, 3)
@test size(out_gd2[1][2]) == (40, 3)
@test size(out_gd2[1][3]) == (50, 3)
@test size(out_gd2[2][1]) == (3, )
@test size(out_gd2[2][2]) == (3, )
@test size(out_gd2[2][3]) == (3, )
@test size(out_gd2[3][1]) == (300, 3)
@test size(out_gd2[3][2]) == (300, 3)
@test size(out_gd2[3][3]) == (300, 3)
@test size(out_gd2[4][1]) == ()
@test size(out_gd2[4][2]) == ()
@test size(out_gd2[4][3]) == ()
@test size(out_gd2[5][1]) == ()
@test size(out_gd2[5][2]) == ()
@test size(out_gd2[5][3]) == ()
@test size(out_gd2[6][1]) == ()
@test size(out_gd2[6][2]) == ()
@test size(out_gd2[6][3]) == ()
@test size(out_gd2[7][1]) == ()
@test size(out_gd2[7][2]) == ()
@test size(out_gd2[7][3]) == ()
@test size(out_gd2[8]) == ()

@test size(out_gd3[1][1]) == (30, 3)
@test size(out_gd3[1][2]) == (40, 3)
@test size(out_gd3[1][3]) == (50, 3)
@test size(out_gd3[2][1]) == (3, )
@test size(out_gd3[2][2]) == (3, )
@test size(out_gd3[2][3]) == (3, )
@test size(out_gd3[3][1]) == (300, 3)
@test size(out_gd3[3][2]) == (300, 3)
@test size(out_gd3[3][3]) == (300, 3)
@test size(out_gd3[4][1]) == ()
@test size(out_gd3[4][2]) == ()
@test size(out_gd3[4][3]) == ()
@test size(out_gd3[5][1]) == ()
@test size(out_gd3[5][2]) == ()
@test size(out_gd3[5][3]) == ()
@test size(out_gd3[6][1]) == ()
@test size(out_gd3[6][2]) == ()
@test size(out_gd3[6][3]) == ()
@test size(out_gd3[7][1]) == ()
@test size(out_gd3[7][2]) == ()
@test size(out_gd3[7][3]) == ()
@test size(out_gd3[8]) == ()

@test size(out_gd4[1][1]) == (30, 3)
@test size(out_gd4[1][2]) == (40, 3)
@test size(out_gd4[1][3]) == (50, 3)
@test size(out_gd4[2][1]) == (3, )
@test size(out_gd4[2][2]) == (3, )
@test size(out_gd4[2][3]) == (3, )
@test size(out_gd4[3][1]) == (300, 3)
@test size(out_gd4[3][2]) == (300, 3)
@test size(out_gd4[3][3]) == (300, 3)
@test size(out_gd4[4][1]) == ()
@test size(out_gd4[4][2]) == ()
@test size(out_gd4[4][3]) == ()
@test size(out_gd4[5][1]) == ()
@test size(out_gd4[5][2]) == ()
@test size(out_gd4[5][3]) == ()
@test size(out_gd4[6][1]) == ()
@test size(out_gd4[6][2]) == ()
@test size(out_gd4[6][3]) == ()
@test size(out_gd4[7][1]) == ()
@test size(out_gd4[7][2]) == ()
@test size(out_gd4[7][3]) == ()
@test size(out_gd4[8]) == ()
####################################

####################################
println("####### GD (Command line) #######")
run(`$(julia) $(joinpath(bindir, "gd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "gd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "gd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "gd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))
####################################

####################################
println("####### SGD (Julia API) #######")
out_sgd1 = sgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="robbins-monro",
  stepsize=1.0e-15,
  numepoch=10,
  numbatch=10,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_sgd2 = sgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="momentum",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_sgd3 = sgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="nag",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_sgd4 = sgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="adagrad",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

@test size(out_sgd1[1][1]) == (30, 3)
@test size(out_sgd1[1][2]) == (40, 3)
@test size(out_sgd1[1][3]) == (50, 3)
@test size(out_sgd1[2][1]) == (3, )
@test size(out_sgd1[2][2]) == (3, )
@test size(out_sgd1[2][3]) == (3, )
@test size(out_sgd1[3][1]) == (300, 3)
@test size(out_sgd1[3][2]) == (300, 3)
@test size(out_sgd1[3][3]) == (300, 3)
@test size(out_sgd1[4][1]) == ()
@test size(out_sgd1[4][2]) == ()
@test size(out_sgd1[4][3]) == ()
@test size(out_sgd1[5][1]) == ()
@test size(out_sgd1[5][2]) == ()
@test size(out_sgd1[5][3]) == ()
@test size(out_sgd1[6][1]) == ()
@test size(out_sgd1[6][2]) == ()
@test size(out_sgd1[6][3]) == ()
@test size(out_sgd1[7][1]) == ()
@test size(out_sgd1[7][2]) == ()
@test size(out_sgd1[7][3]) == ()
@test size(out_sgd1[8]) == ()

@test size(out_sgd2[1][1]) == (30, 3)
@test size(out_sgd2[1][2]) == (40, 3)
@test size(out_sgd2[1][3]) == (50, 3)
@test size(out_sgd2[2][1]) == (3, )
@test size(out_sgd2[2][2]) == (3, )
@test size(out_sgd2[2][3]) == (3, )
@test size(out_sgd2[3][1]) == (300, 3)
@test size(out_sgd2[3][2]) == (300, 3)
@test size(out_sgd2[3][3]) == (300, 3)
@test size(out_sgd2[4][1]) == ()
@test size(out_sgd2[4][2]) == ()
@test size(out_sgd2[4][3]) == ()
@test size(out_sgd2[5][1]) == ()
@test size(out_sgd2[5][2]) == ()
@test size(out_sgd2[5][3]) == ()
@test size(out_sgd2[6][1]) == ()
@test size(out_sgd2[6][2]) == ()
@test size(out_sgd2[6][3]) == ()
@test size(out_sgd2[7][1]) == ()
@test size(out_sgd2[7][2]) == ()
@test size(out_sgd2[7][3]) == ()
@test size(out_sgd2[8]) == ()

@test size(out_sgd3[1][1]) == (30, 3)
@test size(out_sgd3[1][2]) == (40, 3)
@test size(out_sgd3[1][3]) == (50, 3)
@test size(out_sgd3[2][1]) == (3, )
@test size(out_sgd3[2][2]) == (3, )
@test size(out_sgd3[2][3]) == (3, )
@test size(out_sgd3[3][1]) == (300, 3)
@test size(out_sgd3[3][2]) == (300, 3)
@test size(out_sgd3[3][3]) == (300, 3)
@test size(out_sgd3[4][1]) == ()
@test size(out_sgd3[4][2]) == ()
@test size(out_sgd3[4][3]) == ()
@test size(out_sgd3[5][1]) == ()
@test size(out_sgd3[5][2]) == ()
@test size(out_sgd3[5][3]) == ()
@test size(out_sgd3[6][1]) == ()
@test size(out_sgd3[6][2]) == ()
@test size(out_sgd3[6][3]) == ()
@test size(out_sgd3[7][1]) == ()
@test size(out_sgd3[7][2]) == ()
@test size(out_sgd3[7][3]) == ()
@test size(out_sgd3[8]) == ()

@test size(out_sgd4[1][1]) == (30, 3)
@test size(out_sgd4[1][2]) == (40, 3)
@test size(out_sgd4[1][3]) == (50, 3)
@test size(out_sgd4[2][1]) == (3, )
@test size(out_sgd4[2][2]) == (3, )
@test size(out_sgd4[2][3]) == (3, )
@test size(out_sgd4[3][1]) == (300, 3)
@test size(out_sgd4[3][2]) == (300, 3)
@test size(out_sgd4[3][3]) == (300, 3)
@test size(out_sgd4[4][1]) == ()
@test size(out_sgd4[4][2]) == ()
@test size(out_sgd4[4][3]) == ()
@test size(out_sgd4[5][1]) == ()
@test size(out_sgd4[5][2]) == ()
@test size(out_sgd4[5][3]) == ()
@test size(out_sgd4[6][1]) == ()
@test size(out_sgd4[6][2]) == ()
@test size(out_sgd4[6][3]) == ()
@test size(out_sgd4[7][1]) == ()
@test size(out_sgd4[7][2]) == ()
@test size(out_sgd4[7][3]) == ()
@test size(out_sgd4[8]) == ()
####################################

####################################
println("####### SGD (Command line) #######")
run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))
####################################

####################################
println("####### RSGD (Julia API) #######")
out_rsgd1 = rsgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="robbins-monro",
  stepsize=1.0e-15,
  numepoch=10,
  numbatch=10,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_rsgd2 = rsgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="momentum",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_rsgd3 = rsgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="nag",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

out_rsgd4 = rsgd(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  scheduling="adagrad",
  stepsize=1.0e-15,
  numepoch=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

@test size(out_rsgd1[1][1]) == (30, 3)
@test size(out_rsgd1[1][2]) == (40, 3)
@test size(out_rsgd1[1][3]) == (50, 3)
@test size(out_rsgd1[2][1]) == (3, )
@test size(out_rsgd1[2][2]) == (3, )
@test size(out_rsgd1[2][3]) == (3, )
@test size(out_rsgd1[3][1]) == (300, 3)
@test size(out_rsgd1[3][2]) == (300, 3)
@test size(out_rsgd1[3][3]) == (300, 3)
@test size(out_rsgd1[4][1]) == ()
@test size(out_rsgd1[4][2]) == ()
@test size(out_rsgd1[4][3]) == ()
@test size(out_rsgd1[5][1]) == ()
@test size(out_rsgd1[5][2]) == ()
@test size(out_rsgd1[5][3]) == ()
@test size(out_rsgd1[6][1]) == ()
@test size(out_rsgd1[6][2]) == ()
@test size(out_rsgd1[6][3]) == ()
@test size(out_rsgd1[7][1]) == ()
@test size(out_rsgd1[7][2]) == ()
@test size(out_rsgd1[7][3]) == ()
@test size(out_rsgd1[8]) == ()

@test size(out_rsgd2[1][1]) == (30, 3)
@test size(out_rsgd2[1][2]) == (40, 3)
@test size(out_rsgd2[1][3]) == (50, 3)
@test size(out_rsgd2[2][1]) == (3, )
@test size(out_rsgd2[2][2]) == (3, )
@test size(out_rsgd2[2][3]) == (3, )
@test size(out_rsgd2[3][1]) == (300, 3)
@test size(out_rsgd2[3][2]) == (300, 3)
@test size(out_rsgd2[3][3]) == (300, 3)
@test size(out_rsgd2[4][1]) == ()
@test size(out_rsgd2[4][2]) == ()
@test size(out_rsgd2[4][3]) == ()
@test size(out_rsgd2[5][1]) == ()
@test size(out_rsgd2[5][2]) == ()
@test size(out_rsgd2[5][3]) == ()
@test size(out_rsgd2[6][1]) == ()
@test size(out_rsgd2[6][2]) == ()
@test size(out_rsgd2[6][3]) == ()
@test size(out_rsgd2[7][1]) == ()
@test size(out_rsgd2[7][2]) == ()
@test size(out_rsgd2[7][3]) == ()
@test size(out_rsgd2[8]) == ()

@test size(out_rsgd3[1][1]) == (30, 3)
@test size(out_rsgd3[1][2]) == (40, 3)
@test size(out_rsgd3[1][3]) == (50, 3)
@test size(out_rsgd3[2][1]) == (3, )
@test size(out_rsgd3[2][2]) == (3, )
@test size(out_rsgd3[2][3]) == (3, )
@test size(out_rsgd3[3][1]) == (300, 3)
@test size(out_rsgd3[3][2]) == (300, 3)
@test size(out_rsgd3[3][3]) == (300, 3)
@test size(out_rsgd3[4][1]) == ()
@test size(out_rsgd3[4][2]) == ()
@test size(out_rsgd3[4][3]) == ()
@test size(out_rsgd3[5][1]) == ()
@test size(out_rsgd3[5][2]) == ()
@test size(out_rsgd3[5][3]) == ()
@test size(out_rsgd3[6][1]) == ()
@test size(out_rsgd3[6][2]) == ()
@test size(out_rsgd3[6][3]) == ()
@test size(out_rsgd3[7][1]) == ()
@test size(out_rsgd3[7][2]) == ()
@test size(out_rsgd3[7][3]) == ()
@test size(out_rsgd3[8]) == ()

@test size(out_rsgd4[1][1]) == (30, 3)
@test size(out_rsgd4[1][2]) == (40, 3)
@test size(out_rsgd4[1][3]) == (50, 3)
@test size(out_rsgd4[2][1]) == (3, )
@test size(out_rsgd4[2][2]) == (3, )
@test size(out_rsgd4[2][3]) == (3, )
@test size(out_rsgd4[3][1]) == (300, 3)
@test size(out_rsgd4[3][2]) == (300, 3)
@test size(out_rsgd4[3][3]) == (300, 3)
@test size(out_rsgd4[4][1]) == ()
@test size(out_rsgd4[4][2]) == ()
@test size(out_rsgd4[4][3]) == ()
@test size(out_rsgd4[5][1]) == ()
@test size(out_rsgd4[5][2]) == ()
@test size(out_rsgd4[5][3]) == ()
@test size(out_rsgd4[6][1]) == ()
@test size(out_rsgd4[6][2]) == ()
@test size(out_rsgd4[6][3]) == ()
@test size(out_rsgd4[7][1]) == ()
@test size(out_rsgd4[7][2]) == ()
@test size(out_rsgd4[7][3]) == ()
@test size(out_rsgd4[8]) == ()
####################################

####################################
println("####### RSGD (Command line) #######")
run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))

run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))
####################################

####################################
println("####### Orthogonal Iteration (Julia API) #######")
out_orthiter = orthiter(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  numepoch=10,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

@test size(out_orthiter[1][1]) == (30, 3)
@test size(out_orthiter[1][2]) == (40, 3)
@test size(out_orthiter[1][3]) == (50, 3)
@test size(out_orthiter[2][1]) == (3, )
@test size(out_orthiter[2][2]) == (3, )
@test size(out_orthiter[2][3]) == (3, )
@test size(out_orthiter[3][1]) == (300, 3)
@test size(out_orthiter[3][2]) == (300, 3)
@test size(out_orthiter[3][3]) == (300, 3)
@test size(out_orthiter[4][1]) == ()
@test size(out_orthiter[4][2]) == ()
@test size(out_orthiter[4][3]) == ()
@test size(out_orthiter[5][1]) == ()
@test size(out_orthiter[5][2]) == ()
@test size(out_orthiter[5][3]) == ()
@test size(out_orthiter[6][1]) == ()
@test size(out_orthiter[6][2]) == ()
@test size(out_orthiter[6][3]) == ()
@test size(out_orthiter[7][1]) == ()
@test size(out_orthiter[7][2]) == ()
@test size(out_orthiter[7][3]) == ()
@test size(out_orthiter[8]) == ()
####################################

####################################
println("####### Orthogonal Iteration (Command line) #######")
run(`$(julia) $(joinpath(bindir, "orthiter")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))
####################################

####################################
println("####### Horst-Jacobi algorithm (Julia API) #######")
out_horst = horst(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  numepoch=10,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")],
  logdir=tmp)

@test size(out_horst[1][1]) == (30, 3)
@test size(out_horst[1][2]) == (40, 3)
@test size(out_horst[1][3]) == (50, 3)
@test size(out_horst[2][1]) == (3, )
@test size(out_horst[2][2]) == (3, )
@test size(out_horst[2][3]) == (3, )
@test size(out_horst[3][1]) == (300, 3)
@test size(out_horst[3][2]) == (300, 3)
@test size(out_horst[3][3]) == (300, 3)
@test size(out_horst[4][1]) == ()
@test size(out_horst[4][2]) == ()
@test size(out_horst[4][3]) == ()
@test size(out_horst[5][1]) == ()
@test size(out_horst[5][2]) == ()
@test size(out_horst[5][3]) == ()
@test size(out_horst[6][1]) == ()
@test size(out_horst[6][2]) == ()
@test size(out_horst[6][3]) == ()
@test size(out_horst[7][1]) == ()
@test size(out_horst[7][2]) == ()
@test size(out_horst[7][3]) == ()
@test size(out_horst[8]) == ()
####################################

####################################
println("####### Horst-Jacobi algorithm (Command line) #######")
run(`$(julia) $(joinpath(bindir, "horst")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv")) --logdir $(tmp)`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))
####################################

####################################
println("####### Out-of-core Multi-set CCA (Julia API) #######")
out_oocmcca = oocmcca(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst"),
    joinpath(tmp, "Z.zst")],
  scale=["ftt", "ftt", "ftt"],
  pseudocount=[1.0, 1.0, 1.0],
  dim=3,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrZ/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv"),
    joinpath(tmp, "SumrZ/Col_FTTVars.csv")])

@test size(out_oocmcca[1][1]) == (30, 3)
@test size(out_oocmcca[1][2]) == (40, 3)
@test size(out_oocmcca[1][3]) == (50, 3)
@test size(out_oocmcca[2][1]) == (3, )
@test size(out_oocmcca[2][2]) == (3, )
@test size(out_oocmcca[2][3]) == (3, )
@test size(out_oocmcca[3][1]) == (300, 3)
@test size(out_oocmcca[3][2]) == (300, 3)
@test size(out_oocmcca[3][3]) == (300, 3)
@test size(out_oocmcca[4][1]) == ()
@test size(out_oocmcca[4][2]) == ()
@test size(out_oocmcca[4][3]) == ()
@test size(out_oocmcca[5][1]) == ()
@test size(out_oocmcca[5][2]) == ()
@test size(out_oocmcca[5][3]) == ()
@test size(out_oocmcca[6][1]) == ()
@test size(out_oocmcca[6][2]) == ()
@test size(out_oocmcca[6][3]) == ()
@test size(out_oocmcca[7][1]) == ()
@test size(out_oocmcca[7][2]) == ()
@test size(out_oocmcca[7][3]) == ()
####################################

####################################
println("####### Out-of-core Multi-set CCA (Command line) #######")
run(`$(julia) $(joinpath(bindir, "oocmcca")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst"))","$(joinpath(tmp, "Z.zst")) --scale ftt,ftt,ftt --pseudocount 1.0,1.0,1.0 --outdir $(tmp) --dim 3 --numepoch 3 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTMeans.csv")) --colvarlist $(joinpath(tmp, "SumrX/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrY/Col_FTTVars.csv"))","$(joinpath(tmp, "SumrZ/Col_FTTVars.csv"))`)

testfilesize(true,
  joinpath(tmp, "Eigen_vectors1.csv"),
  joinpath(tmp, "Eigen_vectors2.csv"),
  joinpath(tmp, "Eigen_vectors3.csv"),
  joinpath(tmp, "Eigen_values1.csv"),
  joinpath(tmp, "Eigen_values2.csv"),
  joinpath(tmp, "Eigen_values3.csv"),
  joinpath(tmp, "ExpCorVar1.csv"),
  joinpath(tmp, "ExpCorVar2.csv"),
  joinpath(tmp, "ExpCorVar3.csv"),
  joinpath(tmp, "LatentVariables1.csv"),
  joinpath(tmp, "LatentVariables2.csv"),
  joinpath(tmp, "LatentVariables3.csv"))
####################################

####################################
println("####### NIPALS (Julia API) #######")
out_nipals = nipals(
  input=[joinpath(tmp, "X.zst"),
    joinpath(tmp, "Y.zst")],
  scale=["ftt", "ftt"],
  pseudocount=[1.0, 1.0],
  dim=3,
  numepoch=10,
  colmeanlist=[joinpath(tmp, "SumrX/Col_FTTMeans.csv"),
    joinpath(tmp, "SumrY/Col_FTTMeans.csv")],
  colvarlist=[joinpath(tmp, "SumrX/Col_FTTVars.csv"),
    joinpath(tmp, "SumrY/Col_FTTVars.csv")],
  logdir=tmp)

@test size(out_nipals[1]) == (300, 3)
@test size(out_nipals[2]) == (30, 3)
@test size(out_nipals[3]) == (300, 3)
@test size(out_nipals[4]) == (40, 3)
@test size(out_nipals[5]) == (3, 3)
####################################

####################################
println("####### NIPALS (Command line) #######")
run(`$(julia) $(joinpath(bindir, "nipals")) --input $(joinpath(tmp, "X.zst"))","$(joinpath(tmp, "Y.zst")) --scale ftt,ftt --pseudocount 1.0,1.0 --outdir $(tmp) --dim 3 --numepoch 10 --colmeanlist $(joinpath(tmp, "SumrX/Col_FTTMeans.csv"))","$(joinpath(tmp, "SumrY/Col_FTTMeans.csv"))`)

testfilesize(true,
  joinpath(tmp, "T.csv"),
  joinpath(tmp, "P.csv"),
  joinpath(tmp, "U.csv"),
  joinpath(tmp, "Q.csv"),
  joinpath(tmp, "B.csv"))
####################################