text:
## TODO micro-optimization, reuse buffers, etc.
## TODO write a version that doesn't assume that B is well-conditionned, and doesn't reuse B applications
## TODO it seems there is a lack of orthogonalization immediately after locking, maybe investigate this to save on some choleskys

# vprintln(args...) = println(args...)  # Uncomment for output
vprintln(args...) = nothing

using LinearAlgebra

# In the LOBPCG we deal with a subspace (X, R, P), where X, R and P are views
# into bigger arrays. We don't store these subspace vectors as a dense matrix
# or a special struct, but just as a tuple.
# For efficiently matrix-matrix products of this subspace with itself or with
# dense matrices, we introduce the functions `bmul` and `boverlap` below.

# Given A as a list of blocks [A1, A2, A3] this forms the matrix-matrix product
# A * B avoiding a concatenation of the blocks to a dense array
@views @timing "block multiplication" function bmul(blocks::Tuple, B)
    res = blocks[1] * B[1:size(blocks[1], 2), :]  # First multiplication
    offset = size(blocks[1], 2)
    for block in blocks[2:end]
        mul!(res, block, B[offset .+ (1:size(block, 2)), :], 1, 1)
        offset += size(block, 2)
    end
    res
end
bmul(A, B::Tuple) = error("not implemented")
bmul(A::Tuple, B::Tuple) = error("not implemented")
bmul(A, B) = A * B

# Given A and B as two lists of blocks [A1, A2, A3], [B1, B2, B3] form the matrix
# A' B optionally assuming that the output is a Hermitian matrix and will be used
# with the Hermitian view, such that only the upper triangle needs to be computed.
@views function boverlap(blocksA::Tuple, blocksB::Tuple; assume_hermitian=false)
    rows = sum(size(blA, 2) for blA in blocksA)
    cols = sum(size(blB, 2) for blB in blocksB)
    ret = similar(blocksA[1], rows, cols)

    orow = 0  # row offset
    for (iA, blA) in enumerate(blocksA)
        ocol = 0  # column offset
        for (iB, blB) in enumerate(blocksB)
            if !assume_hermitian || iA ≤ iB
                mul!(ret[orow .+ (1:size(blA, 2)), ocol .+ (1:size(blB, 2))], blA', blB)
            end
            ocol += size(blB, 2)
        end
        orow += size(blA, 2)
    end
    ret
end
boverlap(blocksA::Tuple, B) = boverlap(blocksA, (B, ))
boverlap(A, B) = A' * B

# Perform a Rayleigh-Ritz for the N first eigenvectors.
@timing function rayleigh_ritz(X, AX, BX, N)
    F = eigen(Hermitian(boverlap(X, AX, assume_hermitian=true)))  # == Hermitian(X' AX)
    F.vectors[:,1:N], F.values[1:N]
end

import LinearAlgebra: cholesky
function cholesky(X::Union{Matrix{ComplexF16}, Hermitian{ComplexF16,Matrix{ComplexF16}}})
    # Cholesky factorization above may promote the type
    # (e.g. Float16 is promoted to Float32. This undoes it)
    # See https://github.com/JuliaLang/julia/issues/16446
    U = cholesky(ComplexF32.(X)).U
    (U=convert.(ComplexF16, U), )
end

# Orthogonalizes X to tol
# Returns the new X, the number of Cholesky factorizations algorithm, and the
# growth factor by which small perturbations of X can have been
# magnified
@timing function ortho(X; tol=2eps(real(eltype(X))))
    local R

    # # Uncomment for "gold standard"
    # U,S,V = svd(X)
    # return U*V', 1, 1

    growth_factor = 1

    success = false
    nchol = 0
    while true
        O = Hermitian(X'X)
        try
            R = cholesky(O).U
            nchol += 1
            success = true
        catch err
            @assert isa(err, PosDefException)
            vprintln("fail")
            # see https://arxiv.org/pdf/1809.11085.pdf for a nice analysis
            # We are not being very clever here; but this should very rarely happen
            # so it should be OK
            α = 100
            nbad = 0
            while true
                O += α*eps(real(eltype(X)))*norm(X)^2*I
                α *= 10
                try
                    R = cholesky(O).U
                    nchol += 1
                    break
                catch err
                    @assert isa(err, PosDefException)
                end
                nbad += 1
                if nbad > 10
                    error("Cholesky shifting is failing badly, this should never happen")
                end
            end
            success = false
        end
        invR = inv(R)
        rmul!(X, invR) # we do not use X/R because we use invR next

        normest(M) = maximum(abs.(diag(M))) + norm(M - Diagonal(diag(M)))

        # We would like growth_factor *= opnorm(inv(R)) but it's too
        # expensive, so we use an upper bound which is sharp enough to
        # be close to 1 when R is close to I, which we need for the
        # outer loop to terminate
        # ||invR|| = ||D + E|| ≤ ||D|| + ||E|| ≤ ||D|| + ||E||_F,
        # where ||.|| is the 2-norm and ||.||_F the Frobenius

        norminvR = normest(invR)
        growth_factor *= norminvR

        # condR = 1/LAPACK.trcon!('I', 'U', 'N', Array(R))
        condR = normest(R)*norminvR # in practice this seems to be an OK estimate

        vprintln("Ortho(X) success? $success ", eps(real(eltype(X)))*condR^2, " < $tol")

        # a good a posteriori error is that X'X - I is eps()*κ(R)^2;
        # in practice this seems to be sometimes very overconservative
        success && eps(real(eltype(X)))*condR^2 < tol && break

        nchol > 10 && error("Ortho(X) is failing badly, this should never happen")
    end

    # @assert norm(X'X - I) < tol

    X, nchol, growth_factor
end

# Randomize the columns of X if the norm is below tol
function drop!(X, tol=2eps(real(eltype(X))))
    dropped = Int[]
    for i=1:size(X,2)
        n = norm(@views X[:,i])
        if n <= tol
            X[:,i] = randn(eltype(X), size(X,1))
            push!(dropped, i)
        end
    end
    dropped
end

# Find X that is orthogonal, and B-orthogonal to Y, up to a tolerance tol.
function ortho(X, Y, BY; tol=2eps(real(eltype(X))))
    Threads.@threads for i=1:size(X,2)
        n = norm(@views X[:,i])
        X[:,i] ./= n
    end

    niter = 1
    ninners = zeros(Int,0)
    while true
        BYX = boverlap(BY, X)   # = BY' X
        X .-= bmul(Y, BYX)  # = Y * BYX
        # If the orthogonalization has produced results below 2eps, we drop them
        # This is to be able to orthogonalize eg [1;0] against [e^iθ;0],
        # as can happen in extreme cases in the ortho(cP, cX)
        dropped = drop!(X)
        @views if dropped != []
            # X[:, dropped] = X[:, dropped] - Y * BY' * X[:, dropped]
            X[:, dropped] .= X[:, dropped] .- bmul(Y, boverlap(BY, X[:, dropped]))
        end
        if norm(BYX) < tol && niter > 1
            push!(ninners, 0)
            break
        end
        X, ninner, growth_factor = ortho(X, tol=tol)
        push!(ninners, ninner)

        # norm(BY'X) < tol && break should be the proper check, but
        # usually unnecessarily costly. Instead, we predict the error
        # according to the growth factor of the orthogonalization.
        # Assuming BY'Y = 0 to machine precision, after the
        # Y-orthogonalization, BY'X is O(eps), so BY'X will be bounded
        # by O(eps * growth_factor).

        # If we're at a fixed point, growth_factor is 1 and if tol >
        # eps(), the loop will terminate, even if BY'Y != 0
        growth_factor*eps(real(eltype(X))) < tol && break

        niter > 10 && error("Ortho(X,Y) is failing badly, this should never happen")
        niter += 1
    end
    vprintln("ortho choleskys: ", ninners) # get how many Choleskys are performed

    # @assert (norm(BY'X)) < tol
    # @assert (norm(X'X-I)) < tol

    X
end

function final_retval(X, AX, resid_history, niter, n_matvec)
    λ = real(diag(X' * AX))
    residuals = AX .- X*Diagonal(λ)
    (λ=λ, X=X,
     residual_norms=[norm(residuals[:, i]) for i in 1:size(residuals, 2)],
     residual_history=resid_history[:, 1:niter+1],
     n_matvec=n_matvec)
end

### The algorithm is Xn+1 = rayleigh_ritz(hcat(Xn, A*Xn, Xn-Xn-1))
### We follow the strategy of Hetmaniuk and Lehoucq, and maintain a B-orthonormal basis Y = (X,R,P)
### After each rayleigh_ritz step, the B-orthonormal X and P are deduced by an orthogonal rotation from Y
### R is then recomputed, and orthonormalized explicitly wrt BX and BP
### We reuse applications of A/B when it is safe to do so, ie only with orthogonal transformations

@timing function LOBPCG(A, X, B=I, precon=((Y, X, R)->R), tol=1e-10, maxiter=100;
                        miniter=1, ortho_tol=2eps(real(eltype(X))),
                        n_conv_check=nothing, display_progress=false)
    N, M = size(X)

    # If N is too small, we will likely get in trouble
    N >= 3M || @warn "Your problem is too small, and LOBPCG might
        fail; use a full diagonalization instead"

    n_conv_check === nothing && (n_conv_check = M)
    resid_history = zeros(real(eltype(X)), M, maxiter+1)
    buf_X = zero(X)
    buf_P = zero(X)

    X = ortho(X, tol=ortho_tol)[1]
    n_matvec = M   # Count number of matrix-vector products
    AX = A*X
    # full_X/AX/BX will always store the full (including locked) X.
    # X/AX/BX only point to the active part
    P = zero(X)
    AP = zero(X)
    R = zero(X)
    AR = zero(X)
    if B != I
        BR = zero(X)
        BX = B*X
        BP = zero(X)
    else
        # The B arrays share the data
        BR = R
        BX = X
        BP = P
    end
    nlocked = 0
    niter = 0  # the first iteration is fake
    λs = @views [(X[:,n]'*AX[:,n]) / (X[:,n]'BX[:,n]) for n=1:M]
    new_X = X
    new_AX = AX
    new_BX = BX
    full_X = X
    full_AX = AX
    full_BX = BX

    while true
        if niter > 0 # first iteration is just to compute the residuals: no X update
            ###  Perform the Rayleigh-Ritz
            mul!(AR, A, R)
            n_matvec += size(R, 2)

            # Form Rayleigh-Ritz subspace
            if niter > 1
                Y = (X, R, P)
                AY = (AX, AR, AP)
                BY = (BX, BR, BP)  # data shared with (X, R, P) in non-general case
            else
                Y = (X, R)
                AY = (AX, AR)
                BY = (BX, BR)  # data shared with (X, R) in non-general case
            end
            cX, λs = rayleigh_ritz(Y, AY, BY, M-nlocked)

            # Update X. By contrast to some other implementations, we
            # wait on updating P because we have to know which vectors
            # to lock (and therefore the residuals) before computing P
            # only for the unlocked vectors. This results in better convergence.
            new_X  = bmul(Y , cX)  # = Y * cX
            new_AX = bmul(AY, cX)  # = AY * cX,   no accuracy loss, since cX orthogonal
            new_BX = (B == I) ? new_X : bmul(BY, cX)
        end

        ### Compute new residuals
        new_R = new_AX .- new_BX .* λs'
        # it is actually a good question of knowing when to
        # precondition. Here seems sensible, but it could plausibly be
        # done before or after
        @views for i=1:size(X,2)
            resid_history[i + nlocked, niter+1] = norm(new_R[:, i])
        end
        vprintln(niter, "   ", resid_history[:, niter+1])
        precondprep!(precon, X)
        ldiv!(precon, new_R)

        ### Compute number of locked vectors
        prev_nlocked = nlocked
        if niter ≥ miniter  # No locking if below miniter
            for i=nlocked+1:M
                if resid_history[i, niter+1] < tol
                    nlocked += 1
                    vprintln("locked $nlocked")
                else
                    # We lock in order, assuming that the lowest
                    # eigenvectors converge first; might be tricky otherwise
                    break
                end
            end
        end

        if display_progress
            println("Iter $niter, converged $(nlocked)/$(n_conv_check), resid ",
                    norm(resid_history[1:n_conv_check, niter+1]))
        end

        if nlocked >= n_conv_check  # Converged!
            X .= new_X  # Update the part of X which is still active
            AX .= new_AX
            return final_retval(full_X, full_AX, resid_history, niter, n_matvec)
        end
        newly_locked = nlocked - prev_nlocked
        active = newly_locked+1:size(X,2) # newly active vectors

        if niter > 0
            ### compute P = Y*cP only for the newly active vectors
            Xn_indices = newly_locked+1:M-prev_nlocked
            # TODO understand this for a potential save of an
            # orthogonalization, see Hetmaniuk & Lehoucq, and Duersch et. al.
            # cP = copy(cX)
            # cP[Xn_indices,:] .= 0
            e = zeros(eltype(X), size(cX, 1), M - prev_nlocked)
            for i in 1:length(Xn_indices)
                e[Xn_indices[i], i] = 1
            end
            cP = cX .- e
            cP = cP[:, Xn_indices]
            # orthogonalize against all Xn (including newly locked)
            cP = ortho(cP, cX, cX, tol=ortho_tol)

            # Get new P
            new_P = bmul(Y, cP)    # = Y * cP
            new_AP = bmul(AY, cP)  # = AY * cP
            if B != I
                new_BP = bmul(BY, cP)  # = BY * cP
            else
                new_BP = new_P
            end
        end

        # Update all X, even newly locked
        X .= new_X
        AX .= new_AX
        if B != I
            BX .= new_BX
        end

        # Sanity check
        for i = 1:size(X, 2)
            if abs(norm(@view X[:, i]) - 1) >= sqrt(eps(real(eltype(X))))
                error("LOBPCG is badly failing to keep the vectors normalized; this should never happen")
            end
        end

        # Restrict all views to active
        @views begin
            X = X[:, active]
            AX = AX[:, active]
            BX = BX[:, active]
            R = new_R[:, active]
            AR = AR[:, active]
            BR = BR[:, active]
            P = P[:, active]
            AP = AP[:, active]
            BP = BP[:, active]
            λs = λs[active]
        end

        # Update newly active P
        if niter > 0
            P .= new_P
            AP .= new_AP
            if B != I
                BP .= new_BP
            end
        end

        # Orthogonalize R wrt all X, newly active P
        if niter > 0
            Z  = (full_X, P)
            BZ = (full_BX, BP)  # data shared with (full_X, P) in non-general case
        else
            Z  = full_X
            BZ = full_BX
        end
        R .= ortho(R, Z, BZ, tol=ortho_tol)

        if B != I
            # At this point R is orthogonalized but not B-orthogonalized.
            # We assume that B is relatively well-conditioned so that R is
            # close to be B-orthogonal. Therefore one step is OK, and B R
            # can be re-used
            mul!(BR, B, R)
            O = Hermitian(R'*BR)
            U = cholesky(O).U
            rdiv!(R, U)
            rdiv!(BR, U)
        end

        niter < maxiter || break
        niter = niter + 1
    end

    final_retval(full_X, full_AX, resid_history, maxiter, n_matvec)
end


text:
function DescriptiveStatistics()
    return DescriptiveStatistics(())
end

function DescriptiveStatistics(arg0::Vector{jdouble})
    return DescriptiveStatistics((Vector{jdouble},), arg0)
end

function DescriptiveStatistics(arg0::jint)
    return DescriptiveStatistics((jint,), arg0)
end

function accept(obj::DescriptiveStatistics, arg0::jdouble)
    return jcall(obj, "accept", void, (jdouble,), arg0)
end

function add_value(obj::DescriptiveStatistics, arg0::jdouble)
    return jcall(obj, "addValue", void, (jdouble,), arg0)
end

function apply(obj::DescriptiveStatistics, arg0::UnivariateStatistic)
    return jcall(obj, "apply", jdouble, (UnivariateStatistic,), arg0)
end

function clear(obj::DescriptiveStatistics)
    return jcall(obj, "clear", void, ())
end

function copy(obj::DescriptiveStatistics)
    return jcall(obj, "copy", DescriptiveStatistics, ())
end

function get_element(obj::DescriptiveStatistics, arg0::jint)
    return jcall(obj, "getElement", jdouble, (jint,), arg0)
end

function get_geometric_mean(obj::DescriptiveStatistics)
    return jcall(obj, "getGeometricMean", jdouble, ())
end

function get_kurtosis(obj::DescriptiveStatistics)
    return jcall(obj, "getKurtosis", jdouble, ())
end

function get_max(obj::DescriptiveStatistics)
    return jcall(obj, "getMax", jdouble, ())
end

function get_mean(obj::DescriptiveStatistics)
    return jcall(obj, "getMean", jdouble, ())
end

function get_min(obj::DescriptiveStatistics)
    return jcall(obj, "getMin", jdouble, ())
end

function get_n(obj::DescriptiveStatistics)
    return jcall(obj, "getN", jlong, ())
end

function get_percentile(obj::DescriptiveStatistics, arg0::jdouble)
    return jcall(obj, "getPercentile", jdouble, (jdouble,), arg0)
end

function get_population_variance(obj::DescriptiveStatistics)
    return jcall(obj, "getPopulationVariance", jdouble, ())
end

function get_quadratic_mean(obj::DescriptiveStatistics)
    return jcall(obj, "getQuadraticMean", jdouble, ())
end

function get_skewness(obj::DescriptiveStatistics)
    return jcall(obj, "getSkewness", jdouble, ())
end

function get_sorted_values(obj::DescriptiveStatistics)
    return jcall(obj, "getSortedValues", Vector{jdouble}, ())
end

function get_standard_deviation(obj::DescriptiveStatistics)
    return jcall(obj, "getStandardDeviation", jdouble, ())
end

function get_sum(obj::DescriptiveStatistics)
    return jcall(obj, "getSum", jdouble, ())
end

function get_sum_of_squares(obj::DescriptiveStatistics)
    return jcall(obj, "getSumOfSquares", jdouble, ())
end

function get_values(obj::DescriptiveStatistics)
    return jcall(obj, "getValues", Vector{jdouble}, ())
end

function get_variance(obj::DescriptiveStatistics)
    return jcall(obj, "getVariance", jdouble, ())
end

function get_window_size(obj::DescriptiveStatistics)
    return jcall(obj, "getWindowSize", jint, ())
end

function remove_most_recent_value(obj::DescriptiveStatistics)
    return jcall(obj, "removeMostRecentValue", void, ())
end

function replace_most_recent_value(obj::DescriptiveStatistics, arg0::jdouble)
    return jcall(obj, "replaceMostRecentValue", jdouble, (jdouble,), arg0)
end

function set_window_size(obj::DescriptiveStatistics, arg0::jint)
    return jcall(obj, "setWindowSize", void, (jint,), arg0)
end

function to_string(obj::DescriptiveStatistics)
    return jcall(obj, "toString", JString, ())
end

ext: jl


text:
function plotbenchmarks(; write_results = OPTIONS["write_results"],
                          test = OPTIONS["test"],
                          blas_num_threads = OPTIONS["blas_num_threads"],
                          blocksparse_num_threads = OPTIONS["blocksparse_num_threads"],
                          maxdims = OPTIONS["maxdims"],
                          cpp_or_julia = OPTIONS["cpp_or_julia"],
                          descriptions = OPTIONS["descriptions"],
                          benchmarks = OPTIONS["benchmarks"],
                          splitblocks = OPTIONS["splitblocks"],
                          cpp_itensor_version = OPTIONS["cpp_itensor_version"],
                          cpp_itensor_dir = OPTIONS["cpp_itensor_dir"],
                          julia_itensor_version = OPTIONS["julia_itensor_version"])
  println(stdout, "Plotting benchmarks with the following options:")
  @show write_results
  @show test
  @show blas_num_threads
  @show blocksparse_num_threads
  print(stdout, "maxdims = ")
  show(stdout, MIME("text/plain"), maxdims)
  println(stdout)
  @show cpp_or_julia
  @show benchmarks
  @show cpp_itensor_version
  @show julia_itensor_version
  println(stdout)

  data_dir = joinpath(pkgdir(@__MODULE__), "data")

  for benchmark in benchmarks
    println("Plotting benchmark $benchmark for bond dimensions $(maxdims[benchmark])")
    println()

    maxdims_benchmark = maxdims[benchmark]
    if test ≠ false
      maxdims_benchmark = maxdims_benchmark[test]
    end

    p = plot(; title = descriptions[benchmark],
             legend = :topleft,
             legendfontsize = 9,
             xlabel = "Maximum bond dimension",
             ylabel = "Computation time (seconds)",
             xtickfontsize = 12,
             ytickfontsize = 12,
             xguidefontsize = 14,
             yguidefontsize = 14)

    markers = [:circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline, :+, :x]
    marker_ind = 1
    for blas_num_thread in blas_num_threads,
        blocksparse_num_thread in blocksparse_num_threads
      marker, marker_ind = iterate(markers, marker_ind)

      times_julia = Float64[]
      times_cpp = Float64[]

      for maxdim in maxdims_benchmark

        # Read the data
        datapath = "bench_$benchmark"
        datapath = joinpath(datapath, "blas_num_threads_$(blas_num_thread)")
        datapath = joinpath(datapath, "blocksparse_num_threads_$(blocksparse_num_thread)")
        if splitblocks
          datapath = joinpath(datapath, "splitblocks_$(splitblocks)")
        end
        datapath = joinpath(datapath, "maxdim_$(maxdim)")
        datapath *= ".txt"
        if isnothing(cpp_or_julia) || cpp_or_julia == "julia"
          filepath = joinpath(data_dir, "julia")
          filepath = joinpath(filepath, "v$julia_itensor_version")
          filepath = joinpath(filepath, datapath)
          time = readdlm(filepath)[]
          push!(times_julia, time)
        end
        if isnothing(cpp_or_julia) || cpp_or_julia == "c++"
          filepath = joinpath(data_dir, "c++")
          filepath = joinpath(filepath, "v$cpp_itensor_version")
          if splitblocks
            # C++ doesn't currently support splitblocks
            datapath = replace(datapath, "splitblocks_$(splitblocks)" => "")
          end
          filepath = joinpath(filepath, datapath)
          time = readdlm(filepath)[]
          push!(times_cpp, time)
        end
      end

      label_details = ""
      if length(blas_num_threads) > 1
        label_details *= ", $blas_num_thread BLAS thread(s)"
      end
      if length(blocksparse_num_threads) > 1
        label_details *= ", $blocksparse_num_thread block sparse thread(s)"
      end
      if splitblocks
        label_details *= ", splitblocks"
      end

      if isnothing(cpp_or_julia) || cpp_or_julia == "julia"
        plot!(p, maxdims_benchmark, times_julia;
              line = (:solid, 4),
              marker = (marker, 8),
              color = :blue,
              label = "Julia" * label_details)
      end
      if isnothing(cpp_or_julia) || cpp_or_julia == "c++"
        if splitblocks
          # C++ doesn't currently support splitblocks
          label_details = replace(label_details, ", splitblocks" => "")
        end

        plot!(p, maxdims_benchmark, times_cpp;
              line = (:dash, 4),
              marker = (marker, 8),
              color = :red,
              label = "C++" * label_details)
      end
    end
    filename = "plot_benchmark_$(benchmark)_maxdims_$(maxdims_benchmark)_blas_num_threads_$(blas_num_threads)_blocksparse_num_threads_$(blocksparse_num_threads)"
    if splitblocks
      filename *= "_splitblocks_$(splitblocks)"
    end
    filename *= ".png"
    filepath = joinpath(pkgdir(@__MODULE__), "plots", filename)
    println(stdout, "Saving plot to $filepath")
    savefig(p, joinpath(pkgdir(@__MODULE__), "plots", filename))
  end
  return nothing
end


text:
module Word2Vec

using Reexport

include("./doctools/tokenize.jl")
export tokenize_document

include("inputprocessing.jl")
export process_input

include("huffman/Huffman.jl")
@reexport using .Huffman

include("contextdata/ContextDataTypes.jl")
@reexport using .ContextDataTypes

include("heirarchical_softmax.jl")
export HeirarchicalSoftmax

include("wordvectors.jl")
export
    # types
    WordVectors,

    # functions
    wordvector
    index, size,
    vocabulary, in_vocabulary,
    get_vector, cosine, cosine_similar_words,
    analogy, analogy_words, similarity

end


text:
const N = 10000000

function list_comprehension(test_list)
    new_list = [t*t for t in test_list]
    return new_list
end

function for_loop(test_list)
    new_list = similar(test_list)
    length = size(test_list)[1]
    for i in 1:length
        new_list[i] = (test_list[i]*test_list[i])
    end
    return new_list
end

function for_loop_append(test_list)
    new_list = []
    for t in test_list
        push!(new_list, t*t)
    end
    return new_list
end

function test()
    test_list = randn(N)

    println("--------------")
    println("リスト内包表記")
    @time array1 = list_comprehension(test_list)
    println("--------------")
    println("for文 (先に確保)")
    @time array2 = for_loop(test_list)
    println("--------------")
    println("for文（push）")
    @time array3 = for_loop_append(test_list)
    println("--------------")
    println("読み取り")
    @time for i in 1:N
        @assert array1[i] == array2[i]
        @assert array1[i] == array3[i]
    end

end

test()


text:
# Controls whether tuples and non-pointerfree immutables, which Julia
# stores as references, are stored inline in compound types when
# possible. Currently this is problematic because Julia fields of these
# types may be undefined.
const INLINE_TUPLE = false
const INLINE_POINTER_IMMUTABLE = false

const JLD_REF_TYPE = JldDatatype(HDF5Datatype(HDF5.H5T_STD_REF_OBJ, false), 0)
const BUILTIN_TYPES = Set([Symbol, Type, UTF16String, BigFloat, BigInt])
const H5CONVERT_DEFINED = ObjectIdDict()
const JLCONVERT_DEFINED = ObjectIdDict()

## Helper functions

# Holds information about the mapping between a Julia and HDF5 type
immutable JldTypeInfo
    dtypes::Vector{JldDatatype}
    offsets::Vector{Int}
    size::Int
end

# Get information about the HDF5 types corresponding to Julia types
function JldTypeInfo(parent::JldFile, types::(Type...), commit::Bool)
    dtypes = Array(JldDatatype, length(types))
    offsets = Array(Int, length(types))
    offset = 0
    for i = 1:length(types)
        dtype = dtypes[i] = h5fieldtype(parent, types[i], commit)
        offsets[i] = offset
        offset += HDF5.h5t_get_size(dtype)
    end
    JldTypeInfo(dtypes, offsets, offset)
end
JldTypeInfo(parent::JldFile, T::ANY, commit::Bool) =
    JldTypeInfo(parent, T.types, commit)

# Write an HDF5 datatype to the file
function commit_datatype(parent::JldFile, dtype::HDF5Datatype, T::ANY)
    pparent = parent.plain
    if !exists(pparent, pathtypes)
        gtypes = g_create(pparent, pathtypes)
    else
        gtypes = pparent[pathtypes]
    end

    id = length(gtypes)+1
    try
        HDF5.t_commit(gtypes, @sprintf("%08d", id), dtype)
    finally
        close(gtypes)
    end
    a_write(dtype, name_type_attr, full_typename(parent, T))

    # Store in map
    parent.jlh5type[T] = JldDatatype(dtype, id)
end

# If parent is nothing, we are creating the datatype in memory for
# validation, so don't commit it
commit_datatype(parent::Nothing, dtype::HDF5Datatype, T::ANY) =
    JldDatatype(dtype, -1)

# The HDF5 library loses track of relationships among committed types
# after the file is saved. We mangle the names by appending a
# sequential identifier so that we can recover these relationships
# later.
mangle_name(jtype::JldDatatype, jlname) =
    jtype.index <= 0 ? string(jlname, "_") : string(jlname, "_", jtype.index)
Base.convert(::Type{HDF5.Hid}, x::JldDatatype) = x.dtype.id

## Serialization of datatypes to JLD
##
## h5fieldtype - gets the JldDatatype corresponding to a given
## Julia type, when the Julia type is stored as an element of an HDF5
## compound type or array. This is the only function that can operate
## on non-leaf types.
##
## h5type - gets the JldDatatype corresponding to an object of the
## given Julia type. For pointerfree types, this is usually the same as
## the h5fieldtype.
##
## h5convert! - converts data from Julia to HDF5 in a buffer. Most
## methods are dynamically generated by gen_h5convert, but methods for
## special built-in types are predefined.
##
## jlconvert - converts data from HDF5 to a Julia object.
##
## jlconvert! - converts data from HDF5 to Julia in a buffer. This is
## only applicable in cases where fields of that type may not be stored
## as references (e.g., not plain types).

## Special types
##
## To create a special serialization of a datatype, one should:
##
## - Define a method of h5fieldtype that dispatches to h5type
## - Define a method of h5type that constructs the type
## - Define a no-op method for gen_h5convert
## - Define h5convert! and jlconvert
## - If the type is an immutable, define jlconvert!
## - Add the type to BUILTIN_TYPES

## HDF5 bits kinds

# This construction prevents these methods from getting called on type unions
@eval typealias BitsKindTypes Union($(map(x->Type{x}, HDF5.HDF5BitsKind.types)...))

h5fieldtype(parent::JldFile, T::BitsKindTypes, ::Bool) =
    h5type(parent, T, false)

h5type(::JldFile, T::BitsKindTypes, ::Bool) =
    JldDatatype(HDF5Datatype(HDF5.hdf5_type_id(T), false), 0)

gen_h5convert(::JldFile, ::BitsKindTypes) = nothing
h5convert!{T<:HDF5.HDF5BitsKind}(out::Ptr, ::JldFile, x::T, ::JldWriteSession) =
    unsafe_store!(convert(Ptr{T}, out), x)

_jlconvert_bits{T}(::Type{T}, ptr::Ptr) = unsafe_load(convert(Ptr{T}, ptr))
_jlconvert_bits!{T}(out::Ptr, ::Type{T}, ptr::Ptr) =
    (unsafe_store!(convert(Ptr{T}, out), unsafe_load(convert(Ptr{T}, ptr))); nothing)

jlconvert(T::BitsKindTypes, ::JldFile, ptr::Ptr) = _jlconvert_bits(T, ptr)
jlconvert!(out::Ptr, T::BitsKindTypes, ::JldFile, ptr::Ptr) = _jlconvert_bits!(out, T, ptr)

## ByteStrings

h5fieldtype{T<:ByteString}(parent::JldFile, ::Type{T}, ::Bool) =
    h5type(parent, T, false)

# Stored as variable-length strings
function h5type{T<:ByteString}(::JldFile, ::Type{T}, ::Bool)
    type_id = HDF5.h5t_copy(HDF5.hdf5_type_id(T))
    HDF5.h5t_set_size(type_id, HDF5.H5T_VARIABLE)
    HDF5.h5t_set_cset(type_id, HDF5.cset(T))
    JldDatatype(HDF5Datatype(type_id, false), 0)
end

gen_h5convert{T<:ByteString}(::JldFile, ::Type{T}) = nothing
h5convert!(out::Ptr, ::JldFile, x::ByteString, ::JldWriteSession) =
    unsafe_store!(convert(Ptr{Ptr{Uint8}}, out), pointer(x))

function jlconvert(T::Union(Type{ASCIIString}, Type{UTF8String}), ::JldFile, ptr::Ptr)
    strptr = unsafe_load(convert(Ptr{Ptr{Uint8}}, ptr))
    n = int(ccall(:strlen, Csize_t, (Ptr{Uint8},), strptr))
    T(pointer_to_array(strptr, n, true))
end

function jlconvert(T::Union(Type{ByteString}), ::JldFile, ptr::Ptr)
    strptr = unsafe_load(convert(Ptr{Ptr{Uint8}}, ptr))
    str = bytestring(strptr)
    c_free(strptr)
    str
end

## UTF16Strings

h5fieldtype(parent::JldFile, ::Type{UTF16String}, commit::Bool) =
    h5type(parent, UTF16String, commit)

# Stored as compound types that contain a vlen
function h5type(parent::JldFile, ::Type{UTF16String}, commit::Bool)
    haskey(parent.jlh5type, UTF16String) && return parent.jlh5type[UTF16String]
    vlen = HDF5.h5t_vlen_create(HDF5.H5T_NATIVE_UINT16)
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, HDF5.h5t_get_size(vlen))
    HDF5.h5t_insert(id, "data_", 0, vlen)
    HDF5.h5t_close(vlen)
    dtype = HDF5Datatype(id, parent.plain)
    commit ? commit_datatype(parent, dtype, UTF16String) : JldDatatype(dtype, -1)
end

gen_h5convert(::JldFile, ::Type{UTF16String}) = nothing
h5convert!(out::Ptr, ::JldFile, x::UTF16String, ::JldWriteSession) =
    unsafe_store!(convert(Ptr{HDF5.Hvl_t}, out), HDF5.Hvl_t(length(x.data), pointer(x.data)))

function jlconvert(::Type{UTF16String}, ::JldFile, ptr::Ptr)
    hvl = unsafe_load(convert(Ptr{HDF5.Hvl_t}, ptr))
    UTF16String(pointer_to_array(convert(Ptr{Uint16}, hvl.p), hvl.len, true))
end

## Symbols

h5fieldtype(parent::JldFile, ::Type{Symbol}, commit::Bool) =
    h5type(parent, Symbol, commit)

# Stored as a compound type that contains a variable length string
function h5type(parent::JldFile, ::Type{Symbol}, commit::Bool)
    haskey(parent.jlh5type, Symbol) && return parent.jlh5type[Symbol]
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 8)
    HDF5.h5t_insert(id, "symbol_", 0, h5fieldtype(parent, UTF8String, commit))
    dtype = HDF5Datatype(id, parent.plain)
    commit ? commit_datatype(parent, dtype, Symbol) : JldDatatype(dtype, -1)
end

gen_h5convert(::JldFile, ::Type{Symbol}) = nothing
function h5convert!(out::Ptr, file::JldFile, x::Symbol, wsession::JldWriteSession)
    str = string(x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end

jlconvert(::Type{Symbol}, file::JldFile, ptr::Ptr) = symbol(jlconvert(UTF8String, file, ptr))

## BigInts and BigFloats

h5fieldtype(parent::JldFile, T::Union(Type{BigInt}, Type{BigFloat}), commit::Bool) =
    h5type(parent, T, commit)

# Stored as a compound type that contains a variable length string
function h5type(parent::JldFile, T::Union(Type{BigInt}, Type{BigFloat}), commit::Bool)
    haskey(parent.jlh5type, T) && return parent.jlh5type[T]
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 8)
    HDF5.h5t_insert(id, "data_", 0, h5fieldtype(parent, ASCIIString, commit))
    dtype = HDF5Datatype(id, parent.plain)
    commit ? commit_datatype(parent, dtype, T) : JldDatatype(dtype, -1)
end

gen_h5convert(::JldFile, ::Union(Type{BigInt}, Type{BigFloat})) = nothing
function h5convert!(out::Ptr, file::JldFile, x::BigInt, wsession::JldWriteSession)
    str = base(62, x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end
function h5convert!(out::Ptr, file::JldFile, x::BigFloat, wsession::JldWriteSession)
    str = string(x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end

jlconvert(::Type{BigInt}, file::JldFile, ptr::Ptr) =
    Base.parseint_nocheck(BigInt, jlconvert(ASCIIString, file, ptr), 62)
jlconvert(::Type{BigFloat}, file::JldFile, ptr::Ptr) =
    BigFloat(jlconvert(ASCIIString, file, ptr))

## Types

h5fieldtype{T<:Type}(parent::JldFile, ::Type{T}, commit::Bool) =
    h5type(parent, Type, commit)

# Stored as a compound type that contains a variable length string
function h5type{T<:Type}(parent::JldFile, ::Type{T}, commit::Bool)
    haskey(parent.jlh5type, Type) && return parent.jlh5type[Type]
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 8)
    HDF5.h5t_insert(id, "typename_", 0, h5fieldtype(parent, UTF8String, commit))
    dtype = HDF5Datatype(id, parent.plain)
    commit ? commit_datatype(parent, dtype, Type) : JldDatatype(dtype, -1)
end

gen_h5convert{T<:Type}(::JldFile, ::Type{T}) = nothing

function h5convert!(out::Ptr, file::JldFile, x::Type, wsession::JldWriteSession)
    str = full_typename(file, x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end

jlconvert{T<:Type}(::Type{T}, file::JldFile, ptr::Ptr) = julia_type(jlconvert(UTF8String, file, ptr))

## Pointers

h5type{T<:Ptr}(parent::JldFile, ::Type{T}, ::Bool) = throw(PointerException())

## Union()

h5fieldtype(parent::JldFile, ::Type{Union()}, ::Bool) = JLD_REF_TYPE

## Arrays

# These show up as having T.size == 0, hence the need for
# specialization.
h5fieldtype{T,N}(parent::JldFile, ::Type{Array{T,N}}, ::Bool) = JLD_REF_TYPE

## User-defined types
##
## Similar to special types, but h5convert!/jl_convert are dynamically
## generated.

## Tuples

if INLINE_TUPLE
    h5fieldtype(parent::JldFile, T::(Type...), commit::Bool) =
        isleaftype(T) ? h5type(parent, T, commit) : JLD_REF_TYPE
else
    h5fieldtype(parent::JldFile, T::(Type...), ::Bool) = JLD_REF_TYPE
end

function h5type(parent::JldFile, T::(Type...), commit::Bool)
    !isa(T, (DataType...)) && unknown_type_err()
    T = T::(DataType...)

    haskey(parent.jlh5type, T) && return parent.jlh5type[T]
    isleaftype(T) || error("unexpected non-leaf type $T")

    typeinfo = JldTypeInfo(parent, T, commit)
    if isopaque(T)
        id = HDF5.h5t_create(HDF5.H5T_OPAQUE, opaquesize(T))
    else
        id = HDF5.h5t_create(HDF5.H5T_COMPOUND, typeinfo.size)
    end
    for i = 1:length(typeinfo.offsets)
        fielddtype = typeinfo.dtypes[i]
        HDF5.h5t_insert(id, mangle_name(fielddtype, i), typeinfo.offsets[i], fielddtype)
    end

    dtype = HDF5Datatype(id, parent.plain)
    if commit
        jlddtype = commit_datatype(parent, dtype, T)
        if isempty(T)
            # to allow recovery of empty tuples, which HDF5 does not allow
            a_write(dtype, "empty", uint8(1))
        end
        jlddtype
    else
        JldDatatype(dtype, -1)
    end
end

function gen_jlconvert(typeinfo::JldTypeInfo, T::(Type...))
    haskey(JLCONVERT_DEFINED, T) && return

    ex = Expr(:block)
    args = ex.args
    tup = Expr(:tuple)
    tupargs = tup.args
    for i = 1:length(typeinfo.dtypes)
        h5offset = typeinfo.offsets[i]
        field = symbol(string("field", i))

        if HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
            push!(args, :($field = read_ref(file, unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset))))
        else
            push!(args, :($field = jlconvert($(T[i]), file, ptr+$h5offset)))
        end
        push!(tupargs, field)
    end
    @eval jlconvert(::Type{$T}, file::JldFile, ptr::Ptr) = ($ex; $tup)
    JLCONVERT_DEFINED[T] = true
    nothing
end

## All other objects

# For cases not defined above: If the type is mutable and non-empty,
# this is a reference. If the type is immutable, this is a type itself.
if INLINE_POINTER_IMMUTABLE
    h5fieldtype(parent::JldFile, T::ANY, commit::Bool) =
        isleaftype(T) && (!T.mutable || T.size == 0) ? h5type(parent, T, commit) : JLD_REF_TYPE
else
    h5fieldtype(parent::JldFile, T::ANY, commit::Bool) =
        isleaftype(T) && (!T.mutable || T.size == 0) && T.pointerfree ? h5type(parent, T, commit) : JLD_REF_TYPE
end

function h5type(parent::JldFile, T::ANY, commit::Bool)
    !isa(T, DataType) && unknown_type_err()
    T = T::DataType

    haskey(parent.jlh5type, T) && return parent.jlh5type[T]
    isleaftype(T) || error("unexpected non-leaf type ", T)

    if isopaque(T)
        # Empty type or non-basic bitstype
        id = HDF5.h5t_create(HDF5.H5T_OPAQUE, opaquesize(T))
    else
        # Compound type
        typeinfo = JldTypeInfo(parent, T.types, commit)
        id = HDF5.h5t_create(HDF5.H5T_COMPOUND, typeinfo.size)
        for i = 1:length(typeinfo.offsets)
            fielddtype = typeinfo.dtypes[i]
            HDF5.h5t_insert(id, mangle_name(fielddtype, T.names[i]), typeinfo.offsets[i], fielddtype)
        end
    end

    dtype = HDF5Datatype(id, parent.plain)
    if commit
        jlddtype = commit_datatype(parent, dtype, T)
        if T.size == 0
            # to allow recovery of empty types, which HDF5 does not allow
            a_write(dtype, "empty", uint8(1))
        end
        jlddtype
    else
        JldDatatype(dtype, -1)
    end
end

# Normal objects
function _gen_jlconvert_type(typeinfo::JldTypeInfo, T::ANY)
    ex = Expr(:block)
    args = ex.args
    for i = 1:length(typeinfo.dtypes)
        h5offset = typeinfo.offsets[i]

        if HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
            push!(args, quote
                ref = unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset)
                if ref != HDF5.HDF5ReferenceObj_NULL
                    out.$(T.names[i]) = read_ref(file, ref)::$(T.types[i])
                end
            end)
        else
            push!(args, :(out.$(T.names[i]) = jlconvert($(T.types[i]), file, ptr+$h5offset)))
        end
    end
    @eval function jlconvert(::Type{$T}, file::JldFile, ptr::Ptr)
        out = ccall(:jl_new_struct_uninit, Any, (Any,), $T)::$T
        $ex
        out
    end
    nothing
end

# Immutables
function _gen_jlconvert_immutable(typeinfo::JldTypeInfo, T::ANY)
    ex = Expr(:block)
    args = ex.args
    jloffsets = fieldoffsets(T)
    for i = 1:length(typeinfo.dtypes)
        h5offset = typeinfo.offsets[i]
        jloffset = jloffsets[i]

        if HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
            obj = gensym("obj")
            push!(args, quote
                ref = unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset)
                local $obj # must keep alive to prevent collection
                if ref == HDF5.HDF5ReferenceObj_NULL
                    unsafe_store!(convert(Ptr{Int}, out)+$jloffset, 0)
                else
                    # The typeassert ensures that the reference type is
                    # valid for this type
                    $obj = read_ref(file, ref)::$(T.types[i])
                    unsafe_store!(convert(Ptr{Ptr{Void}}, out)+$jloffset, pointer_from_objref($obj))
                end
            end)
        elseif uses_reference(T.types[i])
            # Tuple fields and non-pointerfree immutables are stored
            # inline by JLD if INLINE_TUPLE/INLINE_POINTER_IMMUTABLE is
            # true, but not by Julia
            obj = gensym("obj")
            push!(args, quote
                $obj = jlconvert($(T.types[i]), file, ptr+$h5offset)
                unsafe_store!(convert(Ptr{Ptr{Void}}, out)+$jloffset, pointer_from_objref($obj))
            end)
        else
            push!(args, :(jlconvert!(out+$jloffset, $(T.types[i]), file, ptr+$h5offset)))
        end
    end
    @eval begin
        jlconvert!(out::Ptr, ::Type{$T}, file::JldFile, ptr::Ptr) = ($ex; nothing)
        $(
        if T.pointerfree
            quote
                function jlconvert(::Type{$T}, file::JldFile, ptr::Ptr)
                    out = Array($T, 1)
                    jlconvert!(pointer(out), $T, file, ptr)
                    out[1]
                end
            end
        else
            # XXX can this be improved?
            quote
                function jlconvert(::Type{$T}, file::JldFile, ptr::Ptr)
                    out = ccall(:jl_new_struct_uninit, Any, (Any,), $T)::$T
                    jlconvert!(pointer_from_objref(out)+sizeof(Int), $T, file, ptr)
                    out
                end
            end
        end
        )
    end
    nothing
end

function gen_jlconvert(typeinfo::JldTypeInfo, T::ANY)
    haskey(JLCONVERT_DEFINED, T) && return

    if isempty(T.names)
        if T.size == 0
            @eval begin
                jlconvert(::Type{$T}, ::JldFile, ::Ptr) = $T()
                jlconvert!(out::Ptr, ::Type{$T}, ::JldFile, ::Ptr) =
                    unsafe_store!(convert(Ptr{Ptr{Void}}, out), pointer_from_objref($T()))
            end
        else
            @eval begin
               jlconvert(::Type{$T}, ::JldFile, ptr::Ptr) =  _jlconvert_bits($T, ptr)
               jlconvert!(out::Ptr, ::Type{$T}, ::JldFile, ptr::Ptr) =  _jlconvert_bits!(out, $T, ptr)
            end
        end
        nothing
    elseif T.size == 0
        @eval begin
            jlconvert(::Type{$T}, ::JldFile, ::Ptr) = ccall(:jl_new_struct_uninit, Any, (Any,), $T)::$T
            jlconvert!(out::Ptr, ::Type{$T}, ::JldFile, ::Ptr) = nothing
        end
    elseif T.mutable
        _gen_jlconvert_type(typeinfo, T)
    else
        _gen_jlconvert_immutable(typeinfo, T)
    end
    JLCONVERT_DEFINED[T] = true
    nothing
end

## Common functions for all non-special types (including gen_h5convert)

# Whether this datatype should be stored as opaque
isopaque(t::(DataType...)) = isa(t, ())
isopaque(t::DataType) = isempty(t.names)

# The size of this datatype in the HDF5 file (if opaque)
opaquesize(t::(DataType...)) = 1
opaquesize(t::DataType) = max(1, t.size)

# Whether a type that is stored inline in HDF5 should be stored as a
# reference in Julia. This will only be called such that it returns
# true for some unions of special types defined above, unless either
# INLINE_TUPLE or INLINE_POINTER_IMMUTABLE is true.
uses_reference(T::DataType) = !T.pointerfree
uses_reference(::Tuple) = true
uses_reference(::UnionType) = true

unknown_type_err() =
    error("""$T is not of a type supported by JLD
             Please report this error at https://github.com/timholy/HDF5.jl""")

gen_h5convert(parent::JldFile, T) =
    haskey(H5CONVERT_DEFINED, T) || _gen_h5convert(parent, T)

# There is no point in specializing this
function _gen_h5convert(parent::JldFile, T::ANY)
    dtype = parent.jlh5type[T].dtype
    istuple = isa(T, Tuple)
    if istuple
        types = T
    else
        if isopaque(T::DataType)
            if (T::DataType).size == 0
                @eval h5convert!(out::Ptr, ::JldFile, x::$T, ::JldWriteSession) = nothing
            else
                @eval h5convert!(out::Ptr, ::JldFile, x::$T, ::JldWriteSession) =
                    unsafe_store!(convert(Ptr{$T}, out), x)
            end
            return
        end
        types = (T::DataType).types
    end

    getindex_fn = istuple ? (:getindex) : (:getfield)
    ex = Expr(:block)
    args = ex.args
    n = HDF5.h5t_get_nmembers(dtype.id)
    for i = 1:n
        offset = HDF5.h5t_get_member_offset(dtype.id, i-1)
        if HDF5.h5t_get_member_class(dtype.id, i-1) == HDF5.H5T_REFERENCE
            if istuple
                push!(args, :(unsafe_store!(convert(Ptr{HDF5ReferenceObj}, out)+$offset,
                                            write_ref(file, $getindex_fn(x, $i), wsession))))
            else
                push!(args, quote
                    if isdefined(x, $i)
                        ref = write_ref(file, $getindex_fn(x, $i), wsession)
                    else
                        ref = HDF5.HDF5ReferenceObj_NULL
                    end
                    unsafe_store!(convert(Ptr{HDF5ReferenceObj}, out)+$offset, ref)
                end)
            end
        else
            gen_h5convert(parent, types[i])
            push!(args, :(h5convert!(out+$offset, file, $getindex_fn(x, $i), wsession)))
        end
    end
    @eval h5convert!(out::Ptr, file::JldFile, x::$T, wsession::JldWriteSession) = ($ex; nothing)
    H5CONVERT_DEFINED[T] = true
    nothing
end

## Find the corresponding Julia type for a given HDF5 type

# Type mapping function. Given an HDF5Datatype, find (or construct) the
# corresponding Julia type.
function jldatatype(parent::JldFile, dtype::HDF5Datatype)
    class_id = HDF5.h5t_get_class(dtype.id)
    if class_id == HDF5.H5T_STRING
        cset = HDF5.h5t_get_cset(dtype.id)
        if cset == HDF5.H5T_CSET_ASCII
            return ASCIIString
        elseif cset == HDF5.H5T_CSET_UTF8
            return UTF8String
        else
            error("character set ", cset, " not recognized")
        end
    elseif class_id == HDF5.H5T_INTEGER || class_id == HDF5.H5T_FLOAT
        # This can be a performance hotspot
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_DOUBLE) > 0 && return Float64
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT64) > 0 && return Int64
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_FLOAT) > 0 && return Float32
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT32) > 0 && return Int32
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT8) > 0 && return Uint8
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT64) > 0 && return Uint64
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT32) > 0 && return Uint32
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT8) > 0 && return Int8
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT16) > 0 && return Int16
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT16) > 0 && return Uint16
        error("unrecognized integer or float type")
    elseif class_id == HDF5.H5T_COMPOUND || class_id == HDF5.H5T_OPAQUE
        addr = HDF5.objinfo(dtype).addr
        haskey(parent.h5jltype, addr) && return parent.h5jltype[addr]

        typename = a_read(dtype, name_type_attr)
        T = julia_type(typename)
        if T == UnsupportedType
            warn("type $typename not present in workspace; reconstructing")
            T = reconstruct_type(parent, dtype, gensym(typename))
        end

        if !(T in BUILTIN_TYPES)
            # Call jldatatype on dependent types to validate them and
            # define jlconvert
            if class_id == HDF5.H5T_COMPOUND
                for i = 0:HDF5.h5t_get_nmembers(dtype.id)-1
                    member_name = HDF5.h5t_get_member_name(dtype.id, i)
                    idx = rsearchindex(member_name, "_")
                    if idx != sizeof(member_name)
                        member_dtype = HDF5.t_open(parent.plain, string(pathtypes, '/', lpad(member_name[idx+1:end], 8, '0')))
                        jldatatype(parent, member_dtype)
                    end
                end
            end

            gen_jlconvert(JldTypeInfo(parent, T, false), T)
        end

        # Verify that types match
        newtype = h5type(parent, T, false).dtype
        dtype == newtype || throw(TypeMismatchException(typename))

        # Store type in type index
        index = typeindex(parent, addr)
        parent.jlh5type[T] = JldDatatype(dtype, index)
        parent.h5jltype[addr] = T
        T
    else
        error("unrecognized HDF5 datatype class ", class_id)
    end
end

# Create a Julia type based on the HDF5Datatype from the file. Used
# when the type is no longer available.
function reconstruct_type(parent::JldFile, dtype::HDF5Datatype, name::Symbol)
    class_id = HDF5.h5t_get_class(dtype.id)
    if class_id == HDF5.H5T_OPAQUE
        if exists(dtype, "empty")
            @eval (immutable $name; end; $name)
        else
            sz = int(HDF5.h5t_get_size(dtype.id))*8
            @eval (bitstype $sz $name; $name)
        end
    else
        fields = Expr(:block)
        for i = 0:HDF5.h5t_get_nmembers(dtype.id)-1
            member_name = HDF5.h5t_get_member_name(dtype.id, i)
            idx = rsearchindex(member_name, "_")
            field_name = symbol(member_name[1:idx-1])
            if idx != sizeof(member_name)
                member_dtype = HDF5.t_open(parent.plain, string(pathtypes, '/', lpad(member_name[idx+1:end], 8, '0')))
                push!(fields.args, :($(field_name)::$(jldatatype(parent, member_dtype))))
            else
                member_class = HDF5.h5t_get_member_class(dtype.id, i)
                if member_class == HDF5.H5T_REFERENCE
                    push!(fields.args, field_name)
                else
                    member_dtype = HDF5Datatype(HDF5.h5t_get_member_type(dtype.id, i), parent.plain)
                    push!(fields.args, :($(field_name)::$(jldatatype(parent, member_dtype))))
                end
            end
        end

        @eval (immutable $name; $fields; end; $name)
    end
end

# Get the index of a type in the types group. This could be cached, but
# it's already many times faster than calling H5Iget_name with a lot of
# data in the file, and it only needs to be called once per type.
# Revisit if this ever turns out to be a bottleneck.
function typeindex(parent::JldFile, addr::HDF5.Haddr)
    gtypes = parent.plain[pathtypes]
    i = 1
    for x in gtypes
        if HDF5.objinfo(x).addr == addr
            return i
        end
        i += 1
    end
end


text:
using SparseIR
using Test
using Random

@testset "spr.transform" begin
    for stat in [fermion, boson]
        beta = 1e+4
        wmax = 1.0
        eps = 2e-8
        basis = FiniteTempBasis(stat, beta, wmax, eps)
        spr = SparsePoleRepresentation(basis)

        Random.seed!(4711)

        num_poles = 10
        poles = wmax * (2 * rand(num_poles) .- 1)
        coeffs = 2 * rand(num_poles) .- 1
        @assert maximum(abs.(poles)) <= wmax

        Gl = to_IR(SparsePoleRepresentation(basis, poles), coeffs)

        g_spr = from_IR(spr, Gl)

        # Comparison on Matsubara frequencies
        smpl = MatsubaraSampling(basis)
        smpl_for_spr = MatsubaraSampling(spr, smpl.sampling_points)
        giv = evaluate(smpl_for_spr, g_spr)
        giv_ref = evaluate(smpl, Gl; axis=0)
        @test isapprox(giv, giv_ref; atol=300 * eps, rtol=0)

        # Comparison on tau
        smpl_tau = TauSampling(basis)
        gtau = evaluate(smpl_tau, Gl)
        smpl_tau_for_spr = TauSampling(spr)
        gtau2 = evaluate(smpl_tau_for_spr, g_spr)

        @test isapprox(gtau, gtau2; atol=300 * eps, rtol=0)
    end
end


text:
# scan and accumulate

## COV_EXCL_START

# partial scan of individual thread blocks within a grid
# work-efficient implementation after Blelloch (1990)
#
# number of threads needs to be a power-of-2
#
# performance TODOs:
# - shuffle
# - warp-aggregate atomics
# - the ND case is quite a bit slower than the 1D case (not using Cartesian indices,
#   before 35fcbde1f2987023229034370b0c9091e18c4137). optimize or special-case?
function partial_scan(op::Function, output::CuDeviceArray{T}, input::CuDeviceArray,
                      Rdim, Rpre, Rpost, Rother, neutral, init,
                      ::Val{inclusive}=Val(true)) where {T, inclusive}
    threads = blockDim().x
    thread = threadIdx().x
    block = blockIdx().x

    temp = @cuDynamicSharedMem(T, (2*threads,))

    # iterate the main dimension using threads and the first block dimension
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    # iterate the other dimensions using the remaining block dimensions
    j = (blockIdx().z-1) * gridDim().y + blockIdx().y

    if j > length(Rother)
        return
    end

    @inbounds begin
        I = Rother[j]
        Ipre = Rpre[I[1]]
        Ipost = Rpost[I[2]]
    end

    # load input into shared memory (apply `op` to have the correct type)
    @inbounds temp[thread] = if i <= length(Rdim)
        op(neutral, input[Ipre, i, Ipost])
    else
        op(neutral, neutral)
    end

    # build sum in place up the tree
    offset = 1
    d = threads>>1
    while d > 0
        sync_threads()
        @inbounds if thread <= d
            ai = offset * (2*thread-1)
            bi = offset * (2*thread)
            temp[bi] = op(temp[ai], temp[bi])
        end
        offset *= 2
        d >>= 1
    end

    # clear the last element
    @inbounds if thread == 1
        temp[threads] = neutral
    end

    # traverse down tree & build scan
    d = 1
    while d < threads
        offset >>= 1
        sync_threads()
        @inbounds if thread <= d
            ai = offset * (2*thread-1)
            bi = offset * (2*thread)

            t = temp[ai]
            temp[ai] = temp[bi]
            temp[bi] = op(t, temp[bi])
        end
        d *= 2
    end

    sync_threads()

    # write results to device memory
    @inbounds if i <= length(Rdim)
        val = if inclusive
            op(temp[thread], input[Ipre, i, Ipost])
        else
            temp[thread]
        end
        if init !== nothing
            val = op(something(init), val)
        end
        output[Ipre, i, Ipost] = val
    end

    return
end

# aggregate the result of a partial scan by applying preceding block aggregates
function aggregate_partial_scan(op::Function, output::CuDeviceArray,
                                aggregates::CuDeviceArray, Rdim, Rpre, Rpost, Rother,
                                init)
    threads = blockDim().x
    thread = threadIdx().x
    block = blockIdx().x

    # iterate the main dimension using threads and the first block dimension
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    # iterate the other dimensions using the remaining block dimensions
    j = (blockIdx().z-1) * gridDim().y + blockIdx().y

    @inbounds if block > 1 && i <= length(Rdim) && j <= length(Rother)
        I = Rother[j]
        Ipre = Rpre[I[1]]
        Ipost = Rpost[I[2]]

        val = op(aggregates[Ipre, block-1, Ipost], output[Ipre, i, Ipost])
        if init !== nothing
            val = op(something(init), val)
        end
        output[Ipre, i, Ipost] = val
    end

    return
end

## COV_EXCL_STOP

function scan!(f::Function, output::CuArray{T}, input::CuArray;
               dims::Integer, init=nothing, neutral=GPUArrays.neutral_element(f, T)) where {T}
    dims > 0 || throw(ArgumentError("dims must be a positive integer"))
    inds_t = axes(input)
    axes(output) == inds_t || throw(DimensionMismatch("shape of B must match A"))
    dims > ndims(input) && return copyto!(output, input)
    isempty(inds_t[dims]) && return output

    f = cufunc(f)

    # iteration domain across the main dimension
    Rdim = CartesianIndices((size(input, dims),))

    # iteration domain for the other dimensions
    Rpre = CartesianIndices(size(input)[1:dims-1])
    Rpost = CartesianIndices(size(input)[dims+1:end])
    Rother = CartesianIndices((length(Rpre), length(Rpost)))

    # determine how many threads we can launch for the scan kernel
    args = (f, output, input, Rdim, Rpre, Rpost, Rother, neutral, init, Val(true))
    kernel_args = cudaconvert.(args)
    kernel_tt = Tuple{Core.Typeof.(kernel_args)...}
    kernel = cufunction(partial_scan, kernel_tt)
    kernel_config = launch_configuration(kernel.fun; shmem=(threads)->2*threads*sizeof(T))

    # determine the grid layout to cover the other dimensions
    if length(Rother) > 1
        dev = CUDAdrv.device(kernel.fun.mod.ctx)
        max_other_blocks = attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_GRID_DIM_Y)
        blocks_other = (min(length(Rother), max_other_blocks),
                        cld(length(Rother), max_other_blocks))
    else
        blocks_other = (1, 1)
    end

    # does that suffice to scan the array in one go?
    full = nextpow(2, length(Rdim))
    if full <= kernel_config.threads
        @cuda(threads=full, blocks=(1, blocks_other...), shmem=2*full*sizeof(T), name="scan",
              partial_scan(f, output, input, Rdim, Rpre, Rpost, Rother, neutral, init, Val(true)))
    else
        # perform partial scans across the scanning dimension
        partial = prevpow(2, kernel_config.threads)
        blocks_dim = cld(length(Rdim), partial)
        @cuda(threads=partial, blocks=(blocks_dim, blocks_other...), shmem=2*partial*sizeof(T),
              partial_scan(f, output, input, Rdim, Rpre, Rpost, Rother, neutral, init, Val(true)))

        # get the total of each thread block (except the first) of the partial scans
        aggregates = fill(neutral, Base.setindex(size(input), blocks_dim, dims))
        copyto!(aggregates, selectdim(output, dims, partial:partial:length(Rdim)))

        # scan these totals to get totals for the entire partial scan
        accumulate!(f, aggregates, aggregates; dims=dims)

        # add those totals to the partial scan result
        # NOTE: we assume that this kernel requires fewer resources than the scan kernel.
        #       if that does not hold, launch with fewer threads and calculate
        #       the aggregate block index within the kernel itself.
        @cuda(threads=partial, blocks=(blocks_dim, blocks_other...),
              aggregate_partial_scan(f, output, aggregates, Rdim, Rpre, Rpost, Rother, init))

        CuArrays.unsafe_free!(aggregates)
    end

    return output
end

## Base interface

Base._accumulate!(op, output::CuArray, input::CuVector, dims::Nothing, init::Nothing) =
    scan!(op, output, input; dims=1)

Base._accumulate!(op, output::CuArray, input::CuArray, dims::Integer, init::Nothing) =
    scan!(op, output, input; dims=dims)

Base._accumulate!(op, output::CuArray, input::CuVector, dims::Nothing, init::Some) =
    scan!(op, output, input; dims=1, init=init)

Base._accumulate!(op, output::CuArray, input::CuArray, dims::Integer, init::Some) =
    scan!(op, output, input; dims=dims, init=init)

Base.accumulate_pairwise!(op, result::CuVector, v::CuVector) = accumulate!(op, result, v)


text:
function day10()
    part = [0, 0]
    day10lines = filter(!isempty, strip.(readlines("AoCdata/AoC_2021_day10.txt")))
    rightop, vals, scores = [')', ']', '}', '>'], [3, 57, 1197, 25137], Int[]
    for line in day10lines
        corrupted = false
        rights = Char[]
        for c in line
            if c == '('
                push!(rights, ')')
                continue
            elseif c == '['
                push!(rights, ']')
                continue
            elseif c == '{'
                push!(rights, '}')
                continue
            elseif c == '<'
                push!(rights, '>')
                continue
            else
                if isempty(rights) || c != rights[end]
                    part[1] += vals[findfirst(==(c), rightop)]
                    corrupted = true
                else
                    pop!(rights)
                end
            end
        end
        !corrupted && push!(
            scores,
            reduce((x, y) -> 5x + findfirst(==(y), rightop), reverse(rights); init=0),
        )
    end
    part[2] = sort(scores)[(length(scores) + 1) ÷ 2]
    return part
end

part = day10()
println("Part 1:", part[1])
println("Part 2:", part[2])


text:
#==============================================================================#
# AWSSNS.jl
#
# SNS API. See http://aws.amazon.com/documentation/sns/
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

__precompile__()

module AWSSNS

export sns_list_topics, sns_delete_topic, sns_create_topic, sns_subscribe_sqs,
       sns_subscribe_email, sns_subscribe_lambda, sns_publish,
       sns_list_subscriptsion, sns_unsubscribe

using AWSCore
using SymDict
using Retry

sns_arn(aws::AWSConfig, topic_name) = arn(aws, "sns", topic_name)

function sns(aws::AWSConfig, query)
    query["ContentType"] = "JSON"
    do_request(post_request(aws, "sns", "2010-03-31", query))
end

function sns(aws::AWSConfig, action, topic; args...)

    sns(aws, merge(StringDict(args),
                   "Action" => action,
                   "Name" => topic,
                   "TopicArn" => sns_arn(aws, topic)))
end

function sns_list_topics(aws::AWSConfig)
    r = sns(aws, Dict("Action" => "ListTopics"))
    [split(t["TopicArn"],":")[6] for t in r["Topics"]]
end

function sns_create_topic(aws::AWSConfig, topic_name)
    sns(aws, "CreateTopic", topic_name)
end

function sns_delete_topic(aws::AWSConfig, topic_name)
    sns(aws, "DeleteTopic", topic_name)
end

function sns_publish(aws::AWSConfig, topic_name, message, subject="No Subject")

    if length(subject) > 100
        subject = subject[1:100]
    end
    sns(aws, "Publish", topic_name, Message = message, Subject = subject)
end

function sns_subscribe_sqs(aws::AWSConfig, topic_name, queue; raw=false)

    r = sns(aws, Dict("Action" => "Subscribe",
                      "Name" => topic_name,
                      "TopicArn" => sns_arn(aws, topic_name),
                      "Endpoint" => arn(aws, "sqs", queue),
                      "Protocol" => "sqs"))

    if raw
        sns(aws, "SetSubscriptionAttributes", topic_name,
                  SubscriptionArn = r["SubscriptionArn"],
                  AttributeName = "RawMessageDelivery",
                  AttributeValue = "true")
    end

#=
    import AWSSQS: sqs, sqs_get_queue

    This is probably a bad idea.
    It might wipe out a policy already attached to the queue

    q = sqs_get_queue(aws, queue)

    sqs(q, Dict(
        "Action" => "SetQueueAttributes",
        "Attribute.Name" => "Policy",
        "Attribute.Value" => """{
          "Version": "2008-10-17",
          "Statement": [
            {
              "Effect": "Allow",
              "Principal": {
                "AWS": "*"
              },
              "Action": "SQS:SendMessage",
              "Resource": "$(arn(aws, "sqs", queue))",
              "Condition": {
                "ArnEquals": {
                  "aws:SourceArn": "$(sns_arn(aws, topic_name))"
                }
              }
            }
          ]
        }"""
    ))
=#
end

function sns_subscribe_email(aws::AWSConfig, topic_name, email)

    sns(aws, "Subscribe", topic_name, Endpoint = email, Protocol = "email")
end

import AWSLambda: lambda

function sns_subscribe_lambda(aws::AWSConfig, topic_name, lambda_name)

    if ismatch(r"^arn", lambda_name)
        lambda_arn = lambda_name
        lambda_name = split(lambda_arn, ":")[7]
        laws = copy(aws)
        laws[:region] = arn_region(lambda_arn)
    else
        lambda_arn = arn(aws, "lambda", "function:$lambda_name")
        laws = aws
    end

    sns(aws, "Subscribe", topic_name, Endpoint = lambda_arn, Protocol = "lambda")

# FIXME does not scale. End up with:
#   ERROR: LoadError: 400 -- The final policy size is bigger than the limit.
#=
    @protected try

        lambda(laws, "POST"; path="$lambda_name/policy", query=Dict(
               "Action" => "lambda:InvokeFunction",
               "Principal" => "sns.amazonaws.com",
               "SourceArn" => sns_arn(aws, topic_name),
               "StatementId" => "sns_$topic_name"))
    catch e
        @ignore if e.code == "409" end
    end
=#
end

function sns_list_subscriptsion(aws::AWSConfig, topic_name)
    r = sns(aws, "ListSubscriptionsByTopic", topic_name)
    r["Subscriptions"]
end

function sns_unsubscribe(aws::AWSConfig, topic_name, SubscriptionArn)

    sns(aws, "Unsubscribe", topic_name, SubscriptionArn = SubscriptionArn)
end

function sns_unsubscribe(aws::AWSConfig, topic_name, pattern::Regex)

    for s in sns_list_subscriptsion(aws, topic_name)
        if ismatch(pattern, s["Endpoint"])
            sns_unsubscribe(aws, topic_name, s["SubscriptionArn"])
        end
    end
end

end # module AWSSNS

#==============================================================================#
# End of file.
#==============================================================================#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ServiceSpecification <: SwaggerModel
    metricSpecifications::Any # spec type: Union{ Nothing, Vector{MetricSpecification} } # spec name: metricSpecifications

    function ServiceSpecification(;metricSpecifications=nothing)
        o = new()
        validate_property(ServiceSpecification, Symbol("metricSpecifications"), metricSpecifications)
        setfield!(o, Symbol("metricSpecifications"), metricSpecifications)
        o
    end
end # type ServiceSpecification

const _property_map_ServiceSpecification = Dict{Symbol,Symbol}(Symbol("metricSpecifications")=>Symbol("metricSpecifications"))
const _property_types_ServiceSpecification = Dict{Symbol,String}(Symbol("metricSpecifications")=>"Vector{MetricSpecification}")
Base.propertynames(::Type{ ServiceSpecification }) = collect(keys(_property_map_ServiceSpecification))
Swagger.property_type(::Type{ ServiceSpecification }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ServiceSpecification[name]))}
Swagger.field_name(::Type{ ServiceSpecification }, property_name::Symbol) =  _property_map_ServiceSpecification[property_name]

function check_required(o::ServiceSpecification)
    true
end

function validate_property(::Type{ ServiceSpecification }, name::Symbol, val)
end


text:
#=
ndgrid.jl

I know that using ndgrid is not the "julian way" but when porting code
from matlab sometimes it is easier just to use it instead of refactoring
to do some other way.

The better way would be to use CartesianIndices
https://julialang.org/blog/2016/02/iteration

2019-06-18, Jeff Fessler, University of Michigan
=#

export ndgrid

"""
`(xx,yy) = ndgrid(x::AbstractVector{<:Any}, y::AbstractVector{<:Any})`
todo - improve?
"""
function ndgrid(
		x::AbstractVector{<:Any},
		y::AbstractVector{<:Any})

	tmp = Iterators.product(x, y)
	return [p[1] for p in tmp], [p[2] for p in tmp]
end

"""
`(xx,yy,zz) = ndgrid(x::AbstractVector{<:Any}, y::..., z::...)`
todo - improve?
"""
function ndgrid(
		x::AbstractVector{<:Any},
		y::AbstractVector{<:Any},
		z::AbstractVector{<:Any})

	tmp = Iterators.product(x, y, z)
	return [p[1] for p in tmp], [p[2] for p in tmp], [p[3] for p in tmp]
end