text:
# include("utils.jl")
# include("problem.jl")

function solve!(prb::inner_problem)
    # Set up BLAS
    BLAS.set_num_threads(4)

    # prb.X = prb.X
    # prb.y = prb.y
    # prb.S = prb.S

    #=
    Compute the corresponding linear indices for each block.
    =#
    nblocks = length(prb.block_dims)
    blk_sizes = [(elem > 0) ? elem*(elem+1)÷2 : abs(elem) for elem in [0; prb.block_dims]]  # |n| if diagonal, n(n+1)÷2 else
    blk_indices = cumsum(blk_sizes)
    blocks = [blk_indices[j]+1:blk_indices[j+1] for j = 1:nblocks]                          # linear indices for each block
    n = sum(blk_sizes)
    N = sum(abs.(prb.block_dims))
    ns = [0; prb.block_dims]

    #=
    Allocate the residuals and compute initial values.
    =#
    rp = Vector{Float64}(undef, prb.m)
    Rd = Vector{Float64}(undef, n)
    Rc = Vector{Float64}(undef, n)

    BLAS.blascopy!(prb.m, prb.b, 1, rp, 1)  # b ← rp
    BLAS.blascopy!(n, prb.C, 1, Rd, 1)      # svec(C) ← Rd
    BLAS.axpy!(-1.0, prb.S, Rd)                # Rd - prb.S ← Rd

    MKLSparse.BLAS.cscmv!('T', -1.0, "GU2F", prb.A, prb.X, 1.0, rp)  # b - 𝓐prb.X ← rp
    MKLSparse.BLAS.cscmv!('N', -1.0, "GU2F", prb.A, prb.y, 1.0, Rd)  # Rd - prb.S - 𝓐ᵀprb.y ← Rd

    #=
    Compute μ and other statistics.
    =#
    μ = BLAS.dot(prb.X, prb.S) / N
    relgap = BLAS.dot(prb.X, prb.S) / (1.0 + max(abs(BLAS.dot(prb.C, prb.X)), abs(BLAS.dot(prb.b, prb.y))))
    pinfeas = norm(rp) / (1.0 + norm(prb.b))
    dinfeas = norm(Rd) / (1.0 + norm(prb.C))
    φ = max(pinfeas, dinfeas)
    solution_relgap = prb.ε_relgap  # default value: 1e-6
    solution_φ = prb.ε_φ            # default value: 1e-6
    g = 0.9                         # if adaptative
    # g = 0.98                        # if constant

    sdp_blocks_indices = zeros(Int64, nblocks)
    n_sdp_blocks = 0

    for j = 1:nblocks
        if ns[j+1] > 0
            sdp_blocks_indices[j] = n_sdp_blocks+1
            n_sdp_blocks += 1
        end
    end

    #=
    Preallocate additional memory to perform in-place matrix operations for semidefinite blocks.
    Using Vectors of Matrices is more julianic than using slices.
    =#
    L = [Matrix{Float64}(undef, abs(ns[j+1]), abs(ns[j+1])) for j in 1:nblocks if sdp_blocks_indices[j] > 0]
    R = [Matrix{Float64}(undef, abs(ns[j+1]), abs(ns[j+1])) for j in 1:nblocks if sdp_blocks_indices[j] > 0]
    U = [Matrix{Float64}(undef, abs(ns[j+1]), abs(ns[j+1])) for j in 1:nblocks if sdp_blocks_indices[j] > 0]
    U2 = [Matrix{Float64}(undef, abs(ns[j+1]), abs(ns[j+1])) for j in 1:nblocks if sdp_blocks_indices[j] > 0]
    U3 = [Matrix{Float64}(undef, abs(ns[j+1]), abs(ns[j+1])) for j in 1:nblocks if sdp_blocks_indices[j] > 0]

    # MAX_SDP_BLOCK_SIZE = maximum(abs(ns[j+1]) for j in 1:nblocks if sdp_blocks_indices[j] > 0)
    # SDP_SIZE = [ns[j+1] for j in 1:nblocks if sdp_blocks_indices[j] > 0]
    # U_ = Matrix{Float64}(undef, MAX_SDP_BLOCK_SIZE, MAX_SDP_BLOCK_SIZE)
    # U2_ = Matrix{Float64}(undef, MAX_SDP_BLOCK_SIZE, MAX_SDP_BLOCK_SIZE)
    # U3_ = Matrix{Float64}(undef, MAX_SDP_BLOCK_SIZE, MAX_SDP_BLOCK_SIZE)

    m_indices = [Int64[] for j = 1:nblocks if sdp_blocks_indices[j] > 0]
    for j = 1:nblocks
        if sdp_blocks_indices[j] > 0
            index = sdp_blocks_indices[j]
            for k = 1:prb.m
                if nnz(prb.A[blocks[j], k]) > 0
                    push!(m_indices[index], k)
                end
            end
        end
    end

    #=
    Additional memory to solve for Δy.
    =#
    BB = zeros(Float64, size(prb.A))
    B = zeros(Float64, size(prb.A, 2), size(BB, 2))
    h = Vector{Float64}(undef, prb.m)
    hh = Matrix{Float64}(undef, prb.m, 1)

    H1 = zeros(Float64, size(Rd))
    H2 = zeros(Float64, size(prb.X))

    #=
    Steps
    =#
    δX = Vector{Float64}(undef, size(prb.X))
    δS = Vector{Float64}(undef, size(prb.S))
    δy = zeros(Float64, prb.m)

    αs = zeros(Float64, length(ns)-1)  # Computation of max. step-length
    βs = zeros(Float64, length(ns)-1)

    τ = Vector{Float64}(undef, prb.m)  # Stores the Householder vectors
    iter = 1
    iterMax = prb.maxIter

    #=
    Extra-memory for sparse devectorized matrices Aᵢ.
    =#

    # We take the matrix Aᵢ with most non-zero elements.
    # NNZ_MAX = maximum(nnz(prb.A[:, col]) for col = 1:size(prb.A,2))
    NNZ_MAX = maximum(nnz(prb.A[blocks[j], col]) for col = 1:size(prb.A, 2) for j = 1:nblocks if sdp_blocks_indices[j] > 0)
    # println("new, old: ", NNZ_MAX, ", ", NNZ_MAX_prev)
    # We take the biggest semidefinite block.
    N_MAX = (n_sdp_blocks > 0) ? maximum(abs(ns[j+1]) for j in 1:nblocks if sdp_blocks_indices[j] > 0) : 1
    # The vectors needed for a CSC sparse matrix.
    col_cache = Vector{BlasInt}(undef, N_MAX+1)
    row_cache = Vector{BlasInt}(undef, NNZ_MAX)
    nnz_cache = Vector{Float64}(undef, NNZ_MAX)

    # Print the header if verbose
    if prb.verbose == true
        @printf("iter   |   pobj (scal.)  |   dobj (scal.)\n")
    end

    #=
    Main loop
    =#
    @inbounds while (relgap > solution_relgap) || (φ > solution_φ)
        @views @inbounds for j = 1:nblocks
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]

                fast_smat!(L[index], prb.X[blocks[j]])  # Devectorize the matrices.
                fast_smat!(R[index], prb.S[blocks[j]])

                # BLAS.symm!('L','U', -1.0, L[index], R[index], 0.0, U2[index])
                LAPACK.potrf!('U', L[index])         # Compute Cholesky factorizations.
                LAPACK.potrf!('U', R[index])

                # copyto!(U5[index], 1.0I)             # R\I ← U5
                # LAPACK.trtrs!('U', 'N', 'N', R[index], U5[index])

                # triu!(L[index])
                triu!(R[index])

                # BLAS.gemm!('N','T', 1.0, R[index], L[index], 0.0, U[index])
                # copyto!(U[index], R[index])
                BLAS.blascopy!(length(U[index]), R[index], 1, U[index], 1)
                BLAS.trmm!('R','U','T','N', 1.0, L[index], U[index])
                BLAS.syrk!('U','N', -1.0, U[index], 0.0, U2[index])

                svec!(Rc[blocks[j]], U2[index])
            # @time begin
                # BLAS.blascopy!(SDP_SIZE[index]^2, R[index], 1, U_[1:SDP_SIZE[index],1:SDP_SIZE[index]], 1)
                # BLAS.trmm!('R','U','T','N', 1.0, L[index], U_[1:SDP_SIZE[index],1:SDP_SIZE[index]])
                # BLAS.syrk!('U','N', -1.0, U_[1:SDP_SIZE[index],1:SDP_SIZE[index]], 0.0, U2_[1:SDP_SIZE[index],1:SDP_SIZE[index]])
                #
                # svec!(Rc[blocks[j]], U2_[1:SDP_SIZE[index],1:SDP_SIZE[index]])
            # end
            else
                for i in blocks[j]
                    Rc[i] = -prb.X[i]*prb.S[i]
                end
            end
        end

        #=
        Update the residuals.
        =#
        BLAS.blascopy!(prb.m, prb.b, 1, rp, 1)
        BLAS.blascopy!(n, prb.C, 1, Rd, 1)
        BLAS.axpy!(-1.0, prb.S, Rd)

        MKLSparse.BLAS.cscmv!('T', -1.0, "GU2F", prb.A, prb.X, 1.0, rp)
        MKLSparse.BLAS.cscmv!('N', -1.0, "GU2F", prb.A, prb.y, 1.0, Rd)

        μ = BLAS.dot(prb.X, prb.S) / N
        relgap = BLAS.dot(prb.X, prb.S) / (1.0 + max(abs(BLAS.dot(prb.C, prb.X)), abs(BLAS.dot(prb.b, prb.y))))
        pinfeas = norm(rp) / (1.0 + norm(prb.b))
        dinfeas = norm(Rd) / (1.0 + norm(prb.C))
        φ = max(pinfeas, dinfeas)

        for j = 1:nblocks
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]

                @views fast_smat!(U[index], prb.X[blocks[j]])

                # @views copyto!(U2[index], 1.0I)
                # @views BLAS.trsm!('L','U','T','N', 1.0, R[index], U2[index])
                # @views BLAS.trsm!('L','U','N','N', 1.0, R[index], U2[index])

                # for k = 1:prb.m
                for k in m_indices[index]
                    # if iter== 3 && all(prb.A[blocks[j], k] .== 0.0)
                    #     count_zero += 1
                    # end
                    sparse_A_mul(U2[index], prb.A[blocks[j], k], U[index], col_cache, row_cache, nnz_cache)
                    # @views BLAS.symm!('L','U', 1.0, U2[index], U4[index], 0.0, U3[index])
                    @views BLAS.trsm!('L','U','T','N', 1.0, R[index], U2[index])
                    @views BLAS.trsm!('L','U','N','N', 1.0, R[index], U2[index])
                    @views add_transpose!(U2[index])
                    @views svec!(BB[blocks[j], k], U2[index])
                end

                @views fast_smat!(U2[index], Rd[blocks[j]])
                @views BLAS.symm!('R','U', 1.0, U[index], U2[index], 0.0, U3[index])
                # @views BLAS.symm!('L','U', 1.0, U2[index], U4[index], 0.0, U3[index])
                @views BLAS.trsm!('L','U','T','N', 1.0, R[index], U3[index])
                @views BLAS.trsm!('L','U','N','N', 1.0, R[index], U3[index])
                @views add_transpose!(U3[index])
                @views svec!(H1[blocks[j]], U3[index])

                @views fast_smat!(U2[index], Rc[blocks[j]])
                # @views BLAS.trmm!('R','U','T','N', 1.0, U5[index], U3[index])
                # @views BLAS.trmm!('L','U','N','N', 1.0, U5[index], U3[index])
                @views BLAS.trsm!('L','U','N','N', 1.0, R[index], U2[index])
                @views BLAS.trsm!('R','U','T','N', 1.0, R[index], U2[index])
                @views add_transpose!(U2[index])
                @views svec!(H2[blocks[j]], U2[index])
            else
                @views for i in blocks[j]
                    for k = 1:prb.m
                        BB[i, k] = prb.A[i,k]*prb.X[i]/prb.S[i]
                    end
                    H1[i] = Rd[i]*prb.X[i]/prb.S[i]
                    H2[i] = Rc[i]/prb.S[i]
                end
            end
        end

        MKLSparse.BLAS.cscmm!('T', 1.0, "GU2F", prb.A, BB, 0.0, B)      # 𝓐BB ← B
        BLAS.blascopy!(prb.m, rp, 1, h, 1)
        MKLSparse.BLAS.cscmv!('T', 1.0, "GU2F", prb.A, H1-H2, 1.0, h)   # h

        #=
        Solve the linear system for δy (Predictor step).
        =#
        BLAS.blascopy!(prb.m, h, 1, hh, 1)
        LAPACK.geqrf!(B, τ)                         # QR factorization.
        LAPACK.ormqr!('L','T', B, τ, hh)            # Apply factors to the RHS.
        BLAS.trsm!('L','U','N','N', 1.0, B, hh)     # Triangular solve.
        BLAS.blascopy!(prb.m, hh, 1, δy, 1)         # Copy solution.

        # Compute δS
        BLAS.blascopy!(n, Rd, 1, δS, 1)
        MKLSparse.BLAS.cscmv!('N', -1.0, "GU2F", prb.A, δy, 1.0, δS)

        # Compute δX
        @views for j = 1:nblocks
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]
                fast_smat!(U2[index], δS[blocks[j]]) # S\I * δS * X
                BLAS.symm!('L','U', 1.0, U2[index], U[index], 0.0, U3[index])
                # BLAS.symm!('L','U', 1.0, U2[index], U4[index], 0.0, U3[index])
                BLAS.trsm!('L','U','T','N', 1.0, R[index], U3[index])
                BLAS.trsm!('L','U','N','N', 1.0, R[index], U3[index])
                add_transpose!(U3[index])
                svec!(δX[blocks[j]], U3[index])
            else
                for i in blocks[j]
                    δX[i] = δS[i]*prb.X[i]/prb.S[i]
                end
            end
        end

        BLAS.axpby!(-1.0, prb.X, -1.0, δX)

        #=
        Computation of the minimum eigenvalues.
        =#
        @views for j = 1:length(ns)-1
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]
                fast_smat!(U[index], δX[blocks[j]])

                BLAS.trsm!('R','U','N','N', 1.0, L[index], U[index])
                BLAS.trsm!('L','U','T','N', 1.0, L[index], U[index])
                λ_a = LAPACK.syevr!('N','I','U', U[index], 0.0, 0.0, 1, 1, 1e-6)[1][1]

                fast_smat!(U[index], δS[blocks[j]])
                BLAS.trsm!('R','U','N','N', 1.0, R[index], U[index])
                BLAS.trsm!('L','U','T','N', 1.0, R[index], U[index])
                λ_b = LAPACK.syevr!('N','I','U', U[index], 0.0, 0.0, 1, 1, 1e-6)[1][1]
            else
                λ_a = minimum(δX[blocks[j]] ./ prb.X[blocks[j]])
                λ_b = minimum(δS[blocks[j]] ./ prb.S[blocks[j]])
            end
            αs[j] = (λ_a < 0.0) ? -1.0/λ_a : Inf
            βs[j] = (λ_b < 0.0) ? -1.0/λ_b : Inf
        end

        α = min(1.0, g*minimum(αs))     # max. step lengths
        β = min(1.0, g*minimum(βs))

        #=
        Update of σ.
        =#
        if μ > 1e-6
            if min(α,β) < 1.0/sqrt(3.0)
                exp = 1
            else
                exp = max(1.0, 3*min(α,β)^2)
            end
        else
            exp = 1
        end

        if dot(prb.X+α*δX, prb.S+β*δS) < 0.0
            σ = 0.8
        else
            frac = dot(prb.X+α*δX, prb.S+β*δS)/dot(prb.X, prb.S)
            σ = min(1.0, frac^exp)
        end

        #=
        Update of the residual Rc for the corrector step.
        =#
        @views for j = 1:nblocks
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]
                fast_smat!(U[index], δX[blocks[j]])
                fast_smat!(U2[index], δS[blocks[j]])
                BLAS.symm!('L','U', 1.0, U[index], U2[index], 0.0, U3[index])
                BLAS.trmm!('L','U','N','N', 1.0, R[index], U3[index])
                BLAS.trsm!('R','U','N','N', 1.0, R[index], U3[index])
                add_transpose!(U3[index])
                svec_add!(Rc[blocks[j]], μ*σ*I - U3[index])
            else
                for i in blocks[j]
                    Rc[i] += μ*σ - δX[i]*δS[i]
                end
            end
        end

        @views for j = 1:nblocks
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]
                fast_smat!(U[index], prb.X[blocks[j]])
                fast_smat!(U2[index], Rd[blocks[j]])

                BLAS.symm!('L','U', 1.0, U2[index], U[index], 0.0, U3[index])
                # BLAS.symm!('L','U', 1.0, U3[index], U4[index], 0.0, U2[index])
                BLAS.trsm!('L','U','T','N', 1.0, R[index], U3[index])
                BLAS.trsm!('L','U','N','N', 1.0, R[index], U3[index])

                add_transpose!(U3[index])
                svec!(H1[blocks[j]], U3[index])

                fast_smat!(U3[index], Rc[blocks[j]])
                # BLAS.trmm!('R','U','T','N', 1.0, U5[index], U4[index])
                # BLAS.trmm!('L','U','N','N', 1.0, U5[index], U4[index])
                BLAS.trsm!('L','U','N','N', 1.0, R[index], U3[index])
                BLAS.trsm!('R','U','T','N', 1.0, R[index], U3[index])

                add_transpose!(U3[index])
                svec!(H2[blocks[j]], U3[index])
            else
                for i in blocks[j]
                    H1[i] = Rd[i]*prb.X[i]/prb.S[i]
                    H2[i] = Rc[i]/prb.S[i]
                end
            end
        end

        BLAS.blascopy!(prb.m, rp, 1, h, 1)
        MKLSparse.BLAS.cscmv!('T', 1.0, "GU2F", prb.A, H1-H2, 1.0, h)

        #=
        Solving the linear system for Δy (Corrector step).
        =#
        BLAS.blascopy!(prb.m, h, 1, hh, 1)
        LAPACK.ormqr!('L','T', B, τ, hh)
        BLAS.trsm!('L','U','N','N', 1.0, B, hh)
        BLAS.blascopy!(prb.m, hh, 1, δy, 1)

        # ΔS
        BLAS.blascopy!(n, Rd, 1, δS, 1)
        MKLSparse.BLAS.cscmv!('N', -1.0, "GU2F", prb.A, δy, 1.0, δS)

        # ΔX
        @views @inbounds for j = 1:nblocks
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]
                fast_smat!(U2[index], δS[blocks[j]])

                # copyto!(U3[index], 1.0I)
                # BLAS.trsm!('L','U','T','N', 1.0, R[index], U3[index])
                # BLAS.trsm!('L','U','N','N', 1.0, R[index], U3[index])
                #
                # BLAS.symm!('L','U', 1.0, U2[index], U[index], 0.0, U4[index])
                # BLAS.symm!('L','U', 1.0, U3[index], U4[index], 0.0, U2[index])
                BLAS.trsm!('L','U','T','N', 1.0, R[index], U2[index])
                BLAS.trsm!('L','U','N','N', 1.0, R[index], U2[index])
                BLAS.symm!('R','U', 1.0, U[index], U2[index], 0.0, U3[index])

                add_transpose!(U3[index])
                svec!(H1[blocks[j]], U3[index])
            else
                for i in blocks[j]
                    H1[i] = δS[i]*prb.X[i]/prb.S[i]
                end
            end
        end

        δX .= H2 - H1

        #=
        Minimum eigenvalues for the max. step-length.
        =#
        @views for j = 1:length(ns)-1
            if sdp_blocks_indices[j] > 0
                index = sdp_blocks_indices[j]

                fast_smat!(U[index], δX[blocks[j]])
                BLAS.trsm!('R','U','N','N', 1.0, L[index], U[index])
                BLAS.trsm!('L','U','T','N', 1.0, L[index], U[index])
                λ_a = LAPACK.syev!('N','U', U[index])[1]

                fast_smat!(U[index], δS[blocks[j]])
                BLAS.trsm!('R','U','N','N', 1.0, R[index], U[index])
                BLAS.trsm!('L','U','T','N', 1.0, R[index], U[index])
                λ_b = LAPACK.syev!('N','U', U[index])[1]

            else
                λ_a = minimum(δX[blocks[j]] ./ prb.X[blocks[j]])
                λ_b = minimum(δS[blocks[j]] ./ prb.S[blocks[j]])
            end
            αs[j] = (λ_a < 0.0) ? -1.0/λ_a : Inf
            βs[j] = (λ_b < 0.0) ? -1.0/λ_b : Inf
        end

        α = min(1.0, g*minimum(αs))
        β = min(1.0, g*minimum(βs))

        g = 0.9 + 0.09*min(α, β)  # adaptative update

        BLAS.axpy!(α, δX, prb.X)
        BLAS.axpy!(β, δS, prb.S)
        BLAS.axpy!(β, δy, prb.y)

        pinfeas = norm(rp) / (1.0 + norm(prb.b))
        dinfeas = norm(Rd) / (1.0 + norm(prb.C))

        # Print iteration info if verbose
        if prb.verbose == true
            @printf("%2d     |  %.3e     |  %.3e\n", iter, BLAS.dot(prb.C, prb.X), BLAS.dot(prb.b, prb.y))
        end

        iter += 1
        if iter == iterMax
            break
        end

    end

    # Print the last iteration info
    if prb.verbose == true
        @printf("%2d     |  %.3e     |  %.3e\n", iter, BLAS.dot(prb.C, prb.X), BLAS.dot(prb.b, prb.y))
    end

    # Infeasibility tests
    inftol = 1e-8
    if BLAS.dot(prb.b, prb.y) / norm(prb.A*prb.y + prb.S) > 1/inftol
        println("Primal likely infeasible")
    end
    if -BLAS.dot(prb.C, prb.X) / norm(prb.A'*prb.X) > 1/inftol
        println("Dual likely infeasible")
    end

    prb.num_iters = iter

    @views for j = 1:length(ns)-1
        prb.X[blocks[j]] *= (prb.normb / prb.normsA[j])
        prb.S[blocks[j]] *= (prb.normC * prb.normsA[j])
    end

    prb.y *= prb.normC

    # @. prb.X = prb.X
    # @. prb.y = prb.y
    # @. prb.S = prb.S
    lmul!(prb.normb, prb.b)

    # Unscaling
    @views for j = 1:length(ns)-1
        prb.A[blocks[j],:] *= prb.normsA[j]
        prb.C[blocks[j]] *= prb.normC*prb.normsA[j]
    end

    if prb.verbose == true
        println("")
        println("Primal value: ", BLAS.dot(prb.C, prb.X))
        println("Dual value: ", BLAS.dot(prb.b, prb.y))
        println("N iters: ", iter)
    end
end

max_stars_repo_name: TomaszWasTaken/MySDPSolver

max_issues_repo_name: TomaszWasTaken/MySDPSolver

max_forks_repo_name: TomaszWasTaken/MySDPSolver


text:
using Compat

const LIBMBHFMM2D = string(Base.Filesystem.dirname(Base.source_path()), "/../bin/libmbhfmm2d")

function modbh_qn(x::Array{Float64,1},
                   lambda::Float64,rscale::Float64,n::Int)

    ifders = zero(Int32)
    ders = zeros(Float64,1)
    n = convert(Int32,n)

    n2 = length(x)

    qtemp = zeros(Float64,n+2)
    qn = zeros(Float64,n+1,n2)

    kvec = zeros(Float64,n+5)

    for i = 1:n2
        ccall( (:diffslogbk_fast_,LIBMBHFMM2D), Cvoid,
               (Ref{Float64},Ref{Float64},Ref{Float64},
                Ref{Float64},Ref{Int32},Ref{Float64},
                Ref{Float64},Ref{Int32}),
               x[i],lambda,rscale,qtemp,ifders,ders,
               kvec,n)
        qn[:,i] = qtemp[1:n+1]
    end

    return qn
end

function besselq(n::Int,x::Float64,lambda::Float64;
                 rscale::Float64=1.0)
    xvec = zeros(Float64,1)
    xvec[1] = x
    qn = modbh_qn(xvec,lambda,rscale,n)

    return qn[n+1,1]
end

function modbh_qn_test()
    xtest1 = [0.1]
    lamtest1 = 1.0
    q0test1 = 0.1244839317079709285005145657

    rscale = 1.0
    n = 1

    qn = modbh_qn(xtest1,lamtest1,rscale,n)

    err = abs(q0test1-qn[1,1])
    relerr = err/abs(q0test1)

    return err, relerr
end

function modbh_pn(x::Array{Float64,1},
                   lambda::Float64,rscale::Float64,n::Int)

    ifders = zero(Int32)
    ders = zeros(Float64,1)
    n = convert(Int32,n)

    n2 = length(x)

    ptemp = zeros(Float64,n+2)
    pn = zeros(Float64,n+1,n2)

    ivec = zeros(Float64,n+5)

    for i = 1:n2
        ccall( (:diffszkik_fast_,LIBMBHFMM2D), Cvoid,
               (Ref{Float64},Ref{Float64},Ref{Float64},
                Ref{Float64},Ref{Int32},Ref{Float64},
                Ref{Float64},Ref{Int32}),
               x[i],lambda,rscale,ptemp,ifders,ders,
               ivec,n)
        pn[:,i] = ptemp[1:n+1]
    end

    return pn
end

function modbh_pn_test()
    xtest1 = [0.1]
    lamtest1 = 1.0
    p0test1 = 0.00250156293409560140021055

    rscale = 1.0
    n = 1

    pn = modbh_pn(xtest1,lamtest1,rscale,n)

    err = abs(p0test1-pn[1,1])
    relerr = err/abs(p0test1)

    return err, relerr
end

function besselp(n::Int,x::Float64,lambda::Float64;
                 rscale::Float64=1.0)
    xvec = zeros(Float64,1)
    xvec[1] = x
    pn = modbh_pn(xvec,lambda,rscale,n)

    return pn[n+1,1]
end

function besselj0m1(x::Float64)
    # evaluate J_0(x) - 1 stably

    xh2 = (0.5*x)^2
    if (xh2 <= 1)
        tot = 0.0
        xterm = 1.0
        for i = 1:12
            xterm = -xterm*xh2/(i*i)
            tot = tot + xterm
        end
    else
        tot = besselj(0,x)-1
    end

    return tot
end

function besselj0m1_test()
    testin = 0.1
    testval = -0.002498437933959967718713101525
    funval = besselj0m1(testin)
    err = abs(testval-funval)
    relerr = err/abs(testval)

    return err, relerr
end

ext: jl


text:
two_point_forward(f, x::Real, h::Real=1.0) = (f(x+h) - f(x)) / h
three_point_endpoint(f, x::Real, h::Real=1.0) = (-3f(x) + 4f(x+h)-f(x+2h)) / 2h
three_point_midpoint(f, x::Real, h::Real=1.0) = (f(x+h) - f(x-h)) / 2h
five_point_endpoint(f, x::Real, h::Real=1.0) = (-25f(x) + 48f(x+h) - 36f(x+2h) + 16f(x+3h) - 3f(x+4h)) / 12h
five_point_midpoint(f, x::Real, h::Real=1.0) = (f(x-2h) - 8f(x-h) + 8f(x+h)-f(x+2h)) / 12h

_beta_inc(z1::Real, z2::Real, a::Real, b::Real) = inv(a) * (z2^a * _₂F₁(a,1-b,a+1,z2) - z1^a * _₂F₁(a,1-b,a+1,z1))
_beta_inc(x::Real, a::Real, b::Real) = _beta_inc(0.0, x, a, b)
_beta_inc(a::Real, b::Real) = _beta_inc(0.0, prevfloat(1.0), a, b)

_gsd_quantile(p::Real, F₀::Real, x₀::Real, α::Real, g::Real, k::Real, γ::Real) = x₀ + _beta_inc(F₀^k, p^k, (1-g)/k, 1-γ) / (α*k)

function _gsd_mean(F₀::Real, x₀::Real, α::Real, g::Real, k::Real, γ::Real)
	z1 = F₀^k
	a1 = (1 - g) / k
	b1 = 1 - γ
	a2 = (2 - g) / k
	b2 = 1 - γ
	x₀ + (_beta_inc(z1, prevfloat(1.0), a1, b1) - _beta_inc(a2, b2)) * inv(α * k)
end

function _gsd_moment(F₀::Real, x₀::Real, α::Real, g::Real, k::Real, γ::Real; moment::Int=1)
    z1 = F₀^k
    a = (1 - g) / k
    b = 1 - γ
    c = inv(α * k)

    f = q -> (x₀ + c*_beta_inc(z1, q^k, a, b))^moment
    quadgk(f, 0, 1, atol=1e-4)[1]
end


text:
module EzFFT
using FFTW, Plots

export ezfft, ezplot
export make_log_ln

function make_log_ln(st, ed)
    mxx, zxx = [], []

    stj = floor(Int64, log10(st))
    enj = floor(Int64, log10(ed))

    sti = floor(Int64, st / 10^Float64(stj))
    eni = floor(Int64, ed / 10^Float64(enj))

    for j in stj:enj, i in 1:9
        if j == enj && i > eni
        else
            push!(mxx, i*10^Float64(j))
            x = 10^Float64(j)
            if x >= 1.0
                x = convert(Int64, x)
            end
            push!(zxx, ifelse(i==1, string(x), ""))
        end
    end

    return mxx, zxx
end

function rectangular(dl::Array{Float64,1})
    return dl
end

function ezfft(ds::Array{Float64, 1}, ts::Float64; window=:rectangular)
    maxhz = 1/ts
    minhz = length(ds)/ts

    res = abs.(fft(eval(Expr(:call, window, ds)))./ length(ds) .*2.0)
    res = res[1:round(Int64, length(res)/2.0)]
    popfirst!(res)
    return res
end

function ezplot(ds::Array{Float64, 1}, ts::Float64; ylabel="", xlabel="MHz", fontsize=16)
    if xlabel == "MHz"
        wxx = 1e6
    elseif xlabel == "kHz"
        wxx = 1e3
    else
        wxx = 1.0
    end
    frqmax = round(Int64, 1/ts)
    frqmax = round(Int64, 1/ts)/wxx
    frqmin = frqmax/length(ds)/2
    frqset = (frqmin:frqmin:frqmax/2)

    mxx, zxx = make_log_ln(frqmin, frqmax)
    plot(frqset, ds, xticks=(mxx, zxx), xscale=:log10, legend=false)
    plot!(xlabel="[$(xlabel)]", ylabel="[$(ylabel)]", tickfont=font(fontsize), legendfont=font(fontsize), guidefont=font(fontsize))
end

end # module

max_forks_repo_name: tatssuki/EzFFT


text:
################################
######    Polymorphism    ######
################################
# Expected number of polymorphism above frequency x from the standard diffusion theory
# f(x) = ∫(s) θs* (1/x(1-x)) * ( (ℯ^(Ns)*(1-ℯ^(-4Ns(1-x))) ) / (ℯ^(4Ns-1)))
# Convolote with the binomial to obtain the downsampled sfs
# E[P(x)] = ∑_x=x*→x=1 fB(x*)
############Neutral#############

"""

	sfs_neut()

Expected rate of neutral allele frequency reduce by backgrou	nd selection. The spectrum depends on the number of individual []

```math
\\mathbb{E}[Ps_{(x)}] = \\sum{x^{*}=x}{x^{*}=1}f_{B}(x)
```
# Input
 - `param::parameters`: mutable structure containing the model
 - `SparseMatrixCSC{Float64, Int64}`: convoluted SFS given the B value.
# Return:
 - `Array{Float64}`: expected rate of neutral alleles frequencies.
"""
function sfs_neut(param::parameters,binom::SparseMatrixCSC{Float64,Int64})

	NN2 = convert(Int64,ceil(param.NN*param.B))

	# Allocating variables
	neutral_sfs(i::Int64) = 1.0/(i)

	x = collect(0:NN2)
	solved_neutral_sfs = neutral_sfs.(x)
	replace!(solved_neutral_sfs, Inf => 0.0)

	out::Array{Float64,1} = param.B*(param.θ_coding)*0.25*(binom*solved_neutral_sfs)

	return out
end

############Positive############
# Variable gamma in function changed to s to avoid problem with exported SpecialFunctions.gamma
"""

	sfs_pos(s,p)

Expected rate of positive selected allele frequency reduce by background selection. The spectrum depends on the number of individuals.

# Arguments
 - `param::parameters`: mutable structure containing the model.
 - `s::Int64`: selection strength.
 - `p::Float64`: positive selected alleles probabilty.
 - `binom::SparseMatrixCSC{Float64, Int64}`: convoluted SFS given the B value.
# Return:
 - `Array{Float64}`: expected positive selected alleles frequencies.
"""
function sfs_pos(param::parameters,s::Int64,p::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if p == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else

		red_plus = Φ(param,s)

		# Solving sfs
		NN2 = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		s_corrected = s*param.B

		s_exp1 = exp(s_corrected*2)
		s_exp2 = exp(s_corrected*-2)

		positive_sfs(i::Float64,g1::Float64=s_exp1,g2::Float64=s_exp2,p::Float64=p) = Float64(p*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))

		# Original
		# p*0.5*(ℯ^(2*s_corrected)*(1-ℯ^(-2.0*s_corrected*(1.0-i)))/((ℯ^(2*s_corrected)-1.0)*i*(1.0-i)))

		# Allocating outputs
		solved_positive_sfs::Array{Float64,1} = (1.0/(NN2)) * (positive_sfs.(xa2))
		replace!(solved_positive_sfs, NaN => 0.0)

		out::Array{Float64,1} = (param.θ_coding)*red_plus*0.75*(binom*solved_positive_sfs)

	end

	return out
end

function sfs_pos_float(param::parameters,s::Int64,p::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if p == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else
		red_plus = Φ(param,s)

		# Solving sfs
		NN2 = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		s_corrected = s*param.B
		s_exp1::Quadmath.Float128 = exp(Quadmath.Float128(s_corrected*2))
		s_exp2::Quadmath.Float128 = exp(Quadmath.Float128(s_corrected*-2))

		positive_sfs(i::Float64,g1::Quadmath.Float128=s_exp1,g2::Quadmath.Float128=s_exp2,p::Float64=p) = Float64(p*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))
		# Allocating outputs
		solved_positive_sfs::Array{Float64,1} = (1.0/(NN2)) * (positive_sfs.(xa2))
		replace!(solved_positive_sfs, NaN => 0.0)
		out::Array{Float64,1} = (param.θ_coding)*red_plus*0.75*(binom*solved_positive_sfs)
		# out = out[2:end-1]

	end

	return out
end

######Slightly deleterious######
"""

	sfs_neg(param,p)

Expected rate of positive selected allele frequency reduce by background selection. Spectrum drawn on a gamma DFE. It depends on the number of individuals.

# Arguments
 - `param::parameters`: mutable structure containing the model
 - `p::Float64`: positive selected alleles probabilty.
 - `binom::SparseMatrixCSC{Float64, Int64}`: convoluted SFS given the B value.
# Return:
 - `Array{Float64}`: expected negative selected alleles frequencies.
"""
function sfs_neg(param::parameters,p::Float64,binom::SparseMatrixCSC{Float64,Int64})

	beta     = param.be/(1.0*param.B)
	NN2      = convert(Int64, ceil(param.NN*param.B))
	xa       = collect(0:NN2)/NN2

	solve_z   = similar(xa);

	z(x::Float64,p::Float64=p) = (1.0-p)*(2.0^-param.al)*(beta^param.al)*(-zeta(param.al,x+beta/2.0) + zeta(param.al,(2+beta)/2.0))/((-1.0+x)*x)

	solve_z   = z.(xa)

	if (solve_z[1] == Inf || isnan(solve_z[1]))
		solve_z[1] = 0.0
	end
	if (solve_z[lastindex(solve_z)] == Inf || isnan(solve_z[lastindex(solve_z)]))
		solve_z[lastindex(solve_z)] = 0.0
	end

	solved_negative =  1.0/(NN2+0.0).*solve_z

	out = param.B*(param.θ_coding)*0.75*(binom*solved_negative)

	return out
end

"""
	cumulative_sfs(sfs_tmp)

Changing SFS considering all values above a frequency *x*. The original asymptotic-MK approach takes Pn(x) and Ps(x) as the number of polymorphic sites at frequency *x* rather than above *x*, but this approach scales poorly as sample size increases. We define the polymorphic spectrum as stated above since these quantities trivially have the same asymptote but are less affected by changing sample size.

# Arguments
 - `sfs_tmp::Vector{Float64}`: SFS array
 - `freqs::Bool`: true/false wheter the input SFS array containing a first columns with the frequencies.
# Output
 - `sfs_tmp::Vector{Float64}`: SFS array
"""
function cumulative_sfs(sfs_tmp::Union{Array,SubArray},freqs::Bool=true)

	out      = Array{Float64}(undef, size(sfs_tmp,1),size(sfs_tmp,2))

	if freqs
		idx = 2
	else
		idx = 1
	end

	out[1,idx:end] = sum(sfs_tmp[:,idx:end],dims=1)

	@simd for i in 2:(size(sfs_tmp)[1])

		app = out[i-1,idx:end] .- sfs_tmp[i-1,idx:end]

		if sum(app) > 0.0
			out[i,idx:end] = app
		else
			out[i,idx:end] = zeros(length(app))
		end
	end

	if freqs
		out[:,1] = sfs_tmp[:,1]
	end

	return out
end

"""
	reduce_sfs(sfs_tmp,bins)

Function to reduce the SFS into N bins.

# Arguments
 - `sfs_tmp::Vector{Float64}`: SFS array
 - `bins::Int64`: bin size to collapse the SFS
# Output
 - `sfs_tmp::Vector{Float64}`: Binned SFS array
"""
function reduce_sfs(sfs_tmp::Array,bins::Int64)

	freq  = collect(0:(size(sfs_tmp,1)-1))/size(sfs_tmp,1)
	h1    = StatsBase.fit(StatsBase.Histogram,freq,0:(1/(bins-1)):1)
	xmap1 = StatsBase.binindex.(Ref(h1), freq)

	tmp = hcat(sfs_tmp,xmap1)
	out = Array{Float64}(undef,bins-1 ,size(sfs_tmp,2))
	vIter =  convert(Array,unique(xmap1)')
	@simd for i = eachindex(vIter)
		@inbounds out[i,2:end] = sum(tmp[tmp[:,end] .== i,2:end-1],dims=1)
	end

	out[:,1] = collect(1:(bins-1)) ./ bins

	return (out)
end


text:
using Kalman
using Base.Test
using StaticArrays

#include(joinpath("..", "docs", "make.jl"))

# write your own tests here
include("testgenerate.jl")
include("testfilter.jl")
include("testsmoother.jl")
include("testiterator.jl")
include("testtrack.jl")


text:
function linsolve(operator, b, x₀, alg::BiCGStab, a₀::Number=0, a₁::Number=1)
    # Initial function operation and division defines number type
    y₀ = apply(operator, x₀)
    T = typeof(dot(b, y₀) / norm(b) * one(a₀) * one(a₁))
    α₀ = convert(T, a₀)
    α₁ = convert(T, a₁)
    # Continue computing r = b - a₀ * x₀ - a₁ * operator(x₀)
    r = one(T)*b # r = mul!(similar(b, T), b, 1)
    r = iszero(α₀) ? r : axpy!(-α₀, x₀, r)
    r = axpy!(-α₁, y₀, r)
    x = mul!(similar(r), x₀, 1)
    normr = norm(r)
    S = typeof(normr)

    # Algorithm parameters
    maxiter = alg.maxiter
    tol::S = alg.tol
    numops = 1 # operator has been applied once to determine r
    numiter = 0

    # Check for early return
    if normr < tol
        if alg.verbosity > 0
            @info """BiCGStab linsolve converged without any iterations:
             *  norm of residual = $normr
             *  number of operations = 1"""
        end
        return (x, ConvergenceInfo(1, r, normr, numiter, numops))
    end

    # First iteration
    numiter += 1
    r_shadow = mul!(similar(r), r, 1)     # shadow residual
    ρ = dot(r_shadow, r)

    # Method fails if ρ is zero.
    if ρ ≈ 0.0
        @warn """BiCGStab linsolve errored after $numiter iterations:
        *   norm of residual = $normr
        *   number of operations = $numops"""
        return (x, ConvergenceInfo(0, r, normr, numiter, numops))
    end

    ## BiCG part of the algorithm.
    p = mul!(similar(r), r, 1)
    v = do_apply!(operator, α₀, α₁, p)
    numops += 1

    σ = dot(r_shadow, v)
    α = ρ / σ

    s = mul!(similar(r), r, 1)
    axpy!(-α, v, s)  # half step residual

    xhalf = mul!(similar(x), x, 1)
    axpy!(+α, p, xhalf) # half step iteration

    normr = norm(s)

    # Check for early return at half step.
    if normr < tol
        # Replace approximate residual with the actual residual.
        s = mul!(similar(b), b, 1)
        axpy!(-1, do_apply!(operator, α₀, α₁, xhalf), s)
        numops += 1

        normr_act = norm(s)
        if normr_act < tol
            if alg.verbosity > 0
                @info """BiCGStab linsolve converged at iteration $(numiter-1/2):
                 *  norm of residual = $normr_act
                 *  number of operations = $numops"""
            end
            return (xhalf, ConvergenceInfo(1, s, normr_act, numiter, numops))
        end
    end

    ## GMRES part of the algorithm.
    t = do_apply!(operator, α₀, α₁, s)
    numops += 1

    ω = dot(t, s) / dot(t, t)

    mul!(x, xhalf, 1)
    axpy!(+ω, s, x)    # full step iteration

    mul!(r, s, 1)
    axpy!(-ω, t, r)      # full step residual

    # Check for early return at full step.
    normr = norm(r)
    if normr < tol
        # Replace approximate residual with the actual residual.
        mul!(r, b, 1)
        axpy!(-1, do_apply!(operator, α₀, α₁, x), r)
        numops += 1

        normr_act = norm(r)
        if normr_act < tol
            if alg.verbosity > 0
                @info """BiCGStab linsolve converged at iteration $(numiter):
                *  norm of residual = $normr_act
                *  number of operations = $numops"""
            end
            return (x, ConvergenceInfo(1, r, normr_act, numiter, numops))
        end
    end

    while numiter < maxiter
        if alg.verbosity > 0
            msg = "BiCGStab linsolve in iter $numiter: "
            msg *= "normres = "
            msg *= @sprintf("%12e", normr)
            @info msg
        end

        numiter += 1
        ρold = ρ
        ρ = dot(r_shadow, r)
        β = (ρ / ρold) * (α / ω)

        axpy!(-ω, v, p)
        axpby!(1, r, β, p)

        v = do_apply!(operator, α₀, α₁, p)
        numops += 1

        σ = dot(r_shadow, v)
        α = ρ / σ

        s = mul!(s, r, 1)
        s = axpy!(-α, v, s)     # half step residual

        xhalf = mul!(xhalf, x, 1)
        xhalf = axpy!(+α, p, xhalf)     # half step iteration

        normr = norm(s)

        if alg.verbosity > 0
            msg = "BiCGStab linsolve in iter $(numiter-1/2): "
            msg *= "normres = "
            msg *= @sprintf("%12e", normr)
            @info msg
        end

        # Check for return at half step.
        if normr < tol
            # Compute non-approximate residual.
            #= s = mul!(similar(b), b, 1)
            s = iszero(α₀) ? s : axpy!(-α₀, xhalf, s)
            axpy!(-α₁, apply(operator, xhalf), s)
            numops += 1 =#
            s = mul!(similar(b), b, 1)
            axpy!(-1, do_apply!(operator, α₀, α₁, xhalf), s)
            numops += 1

            normr_act = norm(s)
            if normr_act < tol
                if alg.verbosity > 0
                    @info """BiCGStab linsolve converged at iteration $(numiter-1/2):
                    *  norm of residual = $normr_act
                    *  number of operations = $numops"""
                end
                return (xhalf, ConvergenceInfo(1, s, normr_act, numiter, numops))
            end
        end

        ## GMRES part of the algorithm.
        t = do_apply!(operator, α₀, α₁, s)
        numops += 1

        ω = dot(t, s) / dot(t, t)

        mul!(x, xhalf, 1)
        axpy!(+ω, s, x)    # full step iteration

        mul!(r, s, 1)
        axpy!(-ω, t, r)    # full step residual

        # Check for return at full step.
        normr = norm(r)
        if normr < tol
            # Replace approximate residual with the actual residual.
            mul!(r, b, 1)
            axpy!(-1, do_apply!(operator, α₀, α₁, x), r)
            numops += 1

            normr_act = norm(r)
            if normr_act < tol
                if alg.verbosity > 0
                    @info """BiCGStab linsolve converged at iteration $(numiter):
                    *  norm of residual = $normr_act
                    *  number of operations = $numops"""
                end
                return (x, ConvergenceInfo(1, r, normr_act, numiter, numops))
            end
        end
    end

    if alg.verbosity > 0
        @warn """BiCGStab linsolve finished without converging after $numiter iterations:
        *   norm of residual = $normr
        *   number of operations = $numops"""
    end
    return (x, ConvergenceInfo(0, r, normr, numiter, numops))
end

function do_apply!(operator, α₀, α₁, x)
    y = apply(operator, x)
    if α₀ != zero(α₀) || α₁ != one(α₁)
        axpby!(α₀, x, α₁, y)
    end

    return y
end


text:
"""
    GibbsConditional(sym, conditional)

A "pseudo-sampler" to manually provide analytical Gibbs conditionals to `Gibbs`.
`GibbsConditional(:x, cond)` will sample the variable `x` according to the conditional `cond`, which
must therefore be a function from a `NamedTuple` of the conditioned variables to a `Distribution`.

The `NamedTuple` that is passed in contains all random variables from the model in an unspecified
order, taken from the [`VarInfo`](@ref) object over which the model is run. Scalars and vectors are
stored in their respective shapes. The tuple also contains the value of the conditioned variable
itself, which can be useful, but using it creates something that is not a Gibbs sampler anymore (see
[here](https://github.com/TuringLang/Turing.jl/pull/1275#discussion_r434240387)).

# Examples

```julia
α_0 = 2.0
θ_0 = inv(3.0)
x = [1.5, 2.0]
N = length(x)

@model function inverse_gdemo(x)
    λ ~ Gamma(α_0, θ_0)
    σ = sqrt(1 / λ)
    m ~ Normal(0, σ)
    @. x ~ \$(Normal(m, σ))
end

# The conditionals can be formulated in terms of the following statistics:
x_bar = mean(x) # sample mean
s2 = var(x; mean=x_bar, corrected=false) # sample variance
m_n = N * x_bar / (N + 1)

function cond_m(c)
    λ_n = c.λ * (N + 1)
    σ_n = sqrt(1 / λ_n)
    return Normal(m_n, σ_n)
end

function cond_λ(c)
    α_n = α_0 + (N - 1) / 2 + 1
    β_n = s2 * N / 2 + c.m^2 / 2 + inv(θ_0)
    return Gamma(α_n, inv(β_n))
end

m = inverse_gdemo(x)

sample(m, Gibbs(GibbsConditional(:λ, cond_λ), GibbsConditional(:m, cond_m)), 10)
```
"""
struct GibbsConditional{S, C}
    conditional::C

    function GibbsConditional(sym::Symbol, conditional::C) where {C}
        return new{sym, C}(conditional)
    end
end

DynamicPPL.getspace(::GibbsConditional{S}) where {S} = (S,)

function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:GibbsConditional},
    vi::AbstractVarInfo;
    kwargs...
)
    return nothing, vi
end

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:GibbsConditional},
    vi::AbstractVarInfo;
    kwargs...
)
    condvals = conditioned(tonamedtuple(vi))
    conddist = spl.alg.conditional(condvals)
    updated = rand(rng, conddist)
    # Setindex allows only vectors in this case.
    vi = setindex!!(vi, [updated;], spl)
    # Update log joint probability.
    vi = last(DynamicPPL.evaluate!!(model, rng, vi, SampleFromPrior()))

    return nothing, vi
end

"""
    conditioned(θ::NamedTuple)

Extract a `NamedTuple` of the values in `θ`; i.e., all names of `θ`, mapping to their respective
values.

`θ` is assumed to come from `tonamedtuple(vi)`, which returns a `NamedTuple` of the form

```julia
t = (m = ([0.234, -1.23], ["m[1]", "m[2]"]), λ = ([1.233], ["λ"])
```

and this function implements the cleanup of indexing. `conditioned(t)` will therefore return

```julia
(λ = 1.233, m = [0.234, -1.23])
```
"""
@generated function conditioned(θ::NamedTuple{names}) where {names}
    condvals = [:($n = extractparam(θ.$n)) for n in names]
    return Expr(:tuple, condvals...)
end

"""Takes care of removing the `tonamedtuple` indexing form."""
extractparam(p::Tuple{Vector{<:Array{<:Real}}, Vector{String}}) = foldl(vcat, p[1])
function extractparam(p::Tuple{Vector{<:Real}, Vector{String}})
    values, strings = p
    if length(values) == length(strings) == 1 && !occursin(r".\[.+\]$", strings[1])
        # if m ~ MVNormal([1], I), we could have have ([1], ["m[1]"])!
        return values[1]
    else
        return values
    end
end


text:
includet("OpCodes.jl")
using .OpCodes

function testsettings(data, settings, input)
    ics = [IntComputer(data) for i in 1:5]
    addinput!.(ics, settings)
    for i in eachindex(settings)
        compute!(ics[i], input)
        input = getoutput(ics[i])[1]
    end
    input
end

# do it
using Combinatorics
data = readline("data/Dec7_data.txt")
findmax(map(x->testsettings(data, x, 0), permutations(0:4)))

# Question 2
function circular(data, settings, input)
    channels = [Channel{Int}(Inf) for i in 1:5]
    ics = [IntComputer(data, channels[i], channels[i == 5 ? 1 : i + 1]) for i in 1:5]
    addinput!.(ics, settings)
    addinput!(ics[1], input)
    @sync for i in 1:5
        @async compute!(ics[i])
    end
    getoutput(ics[end])
end

findmax(map(x->circular(data, x, 0), permutations(5:9)))

# Test cases

# for 1
test1 = ps("3,15,3,16,1002,16,10,16,1,16,15,15,4,15,99,0,0")
settings = [4,3,2,1,0]
testsettings(test1, settings, 0)

test2 = ps("3,23,3,24,1002,24,10,24,1002,23,-1,23,101,5,23,23,1,24,23,23,4,23,99,0,0")
settings = 0:4
testsettings(test2, settings, 0)

# for 2
settings = 9:-1:5
data = "3,26,1001,26,-4,26,3,27,1002,27,2,27,1,27,26,27,4,27,1001,28,-1,28,1005,28,6,99,0,0,5"
circular(data, settings, 0)[1]

settings = [9,7,8,5,6]
data = "3,52,1001,52,-5,52,3,53,1,52,56,54,1007,54,5,55,1005,55,26,1001,54,-5,54,1105,1,12,1,53,54,53,1008,54,0,55,1001,55,1,55,2,53,55,53,4,53,1001,56,-1,56,1005,56,6,99,0,0,0,0,10"
circular(data, settings, 0)[1]


text:
#!/usr/bin/env julia

using NarrativeTest

# Handle filesystem paths under Windows.
subs = NarrativeTest.common_subs()
if Sys.iswindows()
    drive = splitdrive(pwd())[1]
    pushfirst!(subs, r" /?…/" => SubstitutionString(" $drive\\\\…\\\\"))
end

NarrativeTest.runtests(subs=subs)


text:
include("/Users/maximilianbittens/Documents/Unity/Unity.jl/src/MeshReaders.jl")
using GeometryTypes
using ColorTypes
import Unity: UnityMesh, PyramidMesh, UnityVector3, UnityText
import MeshReaders: GmshReader, load

file1 = "/Users/maximilianbittens/Documents/IBNM/Dissertation/Modelle/Knochen\ heile.msh"
file2 = "/Users/maximilianbittens/Documents/IBNM/Dissertation/Modelle/Mayo.msh"

bone_mesh = MeshReaders.load(GmshReader, file1)
mayo_mesh = MeshReaders.load(GmshReader, file2)
bone_triangles = vcat(map(x->x[4:end],filter(x->x[1:3]==[2,3,1],bone_mesh.elements))-1...)
mayo_triangles = vcat(map(x->x[4:end],filter(x->x[1:3]==[2,3,1],mayo_mesh.elements))-1...)
bone_lines = vcat(map(x->x[4:end],filter(x->x[1:3]==[1,3,1],bone_mesh.elements))-1...)
mayo_lines = vcat(map(x->x[4:end],filter(x->x[1:3]==[1,3,1],mayo_mesh.elements))-1...)
bone_vertices = UnityVector3[UnityVector3(bone_mesh.nodes[1,i],bone_mesh.nodes[2,i],bone_mesh.nodes[3,i]) for i = 1:size(bone_mesh.nodes,2)]
#mayo_vertices = UnityVector3[UnityVector3(mayo_mesh.nodes[1,i],mayo_mesh.nodes[2,i],mayo_mesh.nodes[3,i]) for i = 1:size(mayo_mesh.nodes,2)]
mayo_vertices = [Point3f0(mayo_mesh.nodes[1,i],mayo_mesh.nodes[2,i],mayo_mesh.nodes[3,i]) for i = 1:size(mayo_mesh.nodes,2)]
#mayo_vertices = map(x->pcamat*x,mayo_vertices)
bone_color = [RGBA{Float32}(.8,.8,.8,.3) for vert in bone_vertices]
mayo_color = [RGBA{Float32}(.8,0,.3,.3) for vert in mayo_vertices]
#mayo_color = [rand(RGBA{Float32}) for vert in mayo_vertices]
#mayo_color = [rand(RGBA{Float32}) for i = 1:3:length(mayo_triangles)]
bone_options = ["surface_shader = wireframe"]
mayo_options = ["surface_shader = transparent"]

unity_text = UnityText("",UnityVector3(-30.,400.,-1418.),UnityVector3(5,5,1.),UnityVector3(90.,180.,1.))

unity_bone_mesh = UnityMesh("Bone",bone_vertices,Int32[],bone_lines,bone_triangles,bone_color,bone_options,UnityText[unity_text])
unity_mayo_mesh = UnityMesh("Mayo",mayo_vertices,Int32[],mayo_lines,mayo_triangles,mayo_color,mayo_options)
unity_mayo_mesh = UnityMesh("Mayo",mayo_vertices,Int32[],Int32[],mayo_triangles,mayo_color,mayo_options)

function plot(um::UnityMesh)
	socket = connect(8052)
	write(socket,um)
	close(socket)
end

function rand!(um::UnityMesh)

	for vert in um.vertices
		vert.x += randn()/2
		vert.y += randn()/2
		vert.z += randn()/2
	end
	um.colors = rand(RGBA{Float32},length(um.vertices))
end

#sleep(1)
#socket = connect(8052)
#for i = 1:10
#	unity_bone_mesh.text[1].text = "$i"
#	rand!(unity_bone_mesh)
#	#plot(unity_bone_mesh)
#	write(socket,unity_bone_mesh)
#	sleep(.75)
#end
#close(socket)


text:
using JsonGrinder, JSON, Test, SparseArrays, Flux, Random, HierarchicalUtils
using JsonGrinder: ExtractScalar, ExtractCategorical, ExtractArray, ExtractDict, ExtractVector
using JsonGrinder: extractempty
using Mill
using Mill: catobs, nobs, MaybeHotMatrix
using Flux: OneHotMatrix
using LinearAlgebra

function less_categorical_scalar_extractor()
	[
	(e -> (keys_len = length(keys(e)); keys_len / e.updated < 0.5 && length(keys(e)) <= 10 && JsonGrinder.is_numeric_or_numeric_string(e)),
		(e, uniontypes) -> ExtractCategorical(keys(e), uniontypes)),
	(e -> JsonGrinder.is_intable(e),
		(e, uniontypes) -> JsonGrinder.extractscalar(Int32, e, uniontypes)),
	(e -> JsonGrinder.is_floatable(e),
	 	(e, uniontypes) -> JsonGrinder.extractscalar(JsonGrinder.FloatType, e, uniontypes)),
	# it's important that condition here would be lower than maxkeys
	(e -> (keys_len = length(keys(e)); keys_len / e.updated < 0.1 && keys_len < 10000 && !JsonGrinder.is_numeric_or_numeric_string(e)),
		(e, uniontypes) -> ExtractCategorical(keys(e), uniontypes)),
	(e -> true,
		(e, uniontypes) -> JsonGrinder.extractscalar(JsonGrinder.unify_types(e), e, uniontypes)),]
end

testing_settings = (; scalar_extractors = less_categorical_scalar_extractor())

function with_emptyismissing(f::Function, a)
    orig_val = Mill.emptyismissing()
    Mill.emptyismissing!(a)
    f()
    Mill.emptyismissing!(orig_val)
end

@testset "ExtractScalar" begin
    @testset "with uniontypes" begin
	    sc = ExtractScalar(Float64,2,3,true)
	    @test sc.uniontypes == true
	    @test all(sc("5").data .== [9])
	    @test all(sc(5).data .== [9])
	    @test all(sc(nothing).data .=== [missing])
	    @test all(sc(missing).data .=== [missing])
	    @test nobs(sc(missing)) == 1
	    @test nobs(sc(nothing)) == 1
	    @test sc(extractempty).data isa Matrix{Union{Missing, Float64}}
	    @test nobs(sc(extractempty)) == 0
	    @test nobs(sc(5)) == 1

	    sc = ExtractScalar(Float32, 0.5, 4.0, true)
	    @test sc.uniontypes == true
	    @test sc(1).data isa Matrix{Union{Missing, Float32}}
	    @test sc(extractempty).data isa Matrix{Union{Missing, Float32}}

	    sc = JsonGrinder.extractscalar(Float32, true)
	    @test sc.uniontypes == true
	    @test sc(1).data isa Matrix{Union{Missing, Float32}}
	    @test sc(Dict(1=>1)).data isa Matrix{Union{Missing, Float32}}
	    @test length(sc) == 1

		@testset "store_input" begin
    		@test sc("5", store_input=true).data == sc("5", store_input=false).data
			@test sc("5", store_input=true).metadata == fill("5",1,1)
			@test isnothing(sc("5", store_input=false).metadata)
			@test sc(5, store_input=true).data == sc(5, store_input=false).data
			@test sc(5, store_input=true).metadata == fill(5,1,1)
			@test isnothing(sc(5, store_input=false).metadata)
			@test sc(nothing, store_input=true).data ≃ sc(nothing, store_input=false).data
			@test sc(nothing, store_input=true).metadata == fill(nothing,1,1)
			@test isnothing(sc(nothing, store_input=false).metadata)
		end
    end

    @testset "without uniontypes" begin
    	sc1 = ExtractScalar(Float64,2,3,false)
    	@test sc1.uniontypes == false
    	@test all(sc1("5").data .== [9])
    	@test all(sc1(5).data .== [9])
    	@test_throws ErrorException sc1(nothing)
    	@test_throws ErrorException sc1(missing)
		@test nobs(sc1(5)) == 1

    	sc2 = ExtractScalar(Float32, 0.5, 4.0, false)
    	@test sc2.uniontypes == false
    	@test sc2(1).data isa Matrix{Float32}
    	sc3 = JsonGrinder.extractscalar(Float32,false)
    	@test sc3.uniontypes == false
    	@test sc3(1).data isa Matrix{Float32}
    	@test_throws ErrorException sc3(Dict(1=>1))
    	@test length(sc3) == 1

		@testset "extractempty" begin
    		@test sc1(extractempty).data isa Matrix{Float64}
    		@test nobs(sc1(extractempty)) == 0
			@test sc2(extractempty).data isa Matrix{Float32}
	    	@test isnothing(sc2(extractempty, store_input=false).metadata)
	    	@test sc2(extractempty, store_input=true).metadata isa Matrix{UndefInitializer}
	    	@test size(sc2(extractempty, store_input=true).metadata) == (1,0)
		end

		@testset "store_input" begin
    		@test sc3("5", store_input=true).data == sc3("5", store_input=false).data
       		@test sc3("5", store_input=true).metadata == fill("5",1,1)
       		@test isnothing(sc3("5", store_input=false).metadata)
       		@test sc3(5, store_input=true).data == sc3(5, store_input=false).data
       		@test sc3(5, store_input=true).metadata == fill(5,1,1)
       		@test isnothing(sc3(5, store_input=false).metadata)
       		@test_throws ErrorException sc3(nothing, store_input=true)
       		@test_throws ErrorException sc3(nothing, store_input=true)
		end
    end
end
# todo: add for each extractor tests with and without uniontypes, and with and without store_input
# todo: categorize tests w.r.t. individual extractors, make a checklist of what each extractor needs to test
@testset "ExtractArray" begin
	sc = ExtractArray(ExtractCategorical(2:4))
	with_emptyismissing(false) do
		e234 = sc([2,3,4], store_input=false)
		e234 = sc([2,3,4], store_input=false)
		en = sc(nothing, store_input=false)
		e234s = sc([2,3,4], store_input=true)
		ens = sc(nothing, store_input=true)
		@test all(e234.data.data .== Matrix(1.0I, 4, 3))
		@test nobs(en.data) == 0
		@test en.data.data isa MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
		@test nobs(en.data.data) == 0
		@test all(en.bags.bags .== [0:-1])

		@testset "store_input" begin
			@test e234s.data.data == e234.data.data
			@test e234s.data.metadata == [2,3,4]
			@test e234s.data[1].metadata == [2]
			@test e234s.data[2].metadata == [3]
			@test e234s.data[3].metadata == [4]
			@test ens.data.data == en.data.data
			@test ens.metadata == [nothing]
			@test isnothing(en.metadata)
			@test isnothing(ens.data.metadata)
			@test isnothing(e234.data.metadata)
			@test isnothing(en.data.metadata)
		end
	end
	with_emptyismissing(true) do
		en = sc(nothing, store_input=false)
		ens = sc(nothing, store_input=true)
		em = sc(missing, store_input=false)
		ems = sc(missing, store_input=true)
		@test en.data isa Missing
		@test en.data ≃ em.data
		@test ens.data isa Missing
		@test ens.data ≃ ems.data
		@test isnothing(en.metadata)
		@test isnothing(em.metadata)
		@test ens.metadata == [nothing]
		@test ems.metadata ≃ [missing]
		@test all(en.bags.bags .== [0:-1])
		@test all(em.bags.bags .== [0:-1])
	end

	@testset "extractempty" begin
		@test nobs(sc(extractempty).data.data) == 0
		@test nobs(sc(extractempty).data) == 0
		@test isempty(sc(extractempty).bags.bags)
		@test sc(extractempty).data.data isa MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
		with_emptyismissing(true) do
			@test nobs(sc(extractempty)) == 0
		end
		with_emptyismissing(false) do
			@test nobs(sc(extractempty)) == 0
		end
	end
	sc = ExtractArray(ExtractScalar(Float32))
	e234 = sc([2,3,4], store_input=false)
	en = sc(nothing, store_input=false)
	e234s = sc([2,3,4], store_input=true)
	ens = sc(nothing, store_input=true)
	@test all(e234.data.data .== [2 3 4])
	@test nobs(en.data) == 0
	@test all(en.bags.bags .== [0:-1])
	@test nobs(sc(Dict(1=>1)).data) == 0

	@test nobs(sc(extractempty).data.data) == 0
	@test nobs(sc(extractempty).data) == 0
	@test isempty(sc(extractempty).bags.bags)
	@test sc(extractempty).data.data isa Matrix{Union{Missing, Float32}}

	@test en.data.data isa Matrix{Union{Missing, Float32}}
	@test e234s.data.data == e234.data.data
	@test e234s.data.metadata == [2 3 4]
	@test e234s.data[1].metadata == fill(2,1,1)
	@test e234s.data[2].metadata == fill(3,1,1)
	@test e234s.data[3].metadata == fill(4,1,1)
	@test ens.data.data == en.data.data
	@test isnothing(en.data.metadata)
	@test ens.metadata == [nothing]
	@test isnothing(e234.data.metadata)
	@test isnothing(en.data.metadata)

	# testing https://github.com/CTUAvastLab/JsonGrinder.jl/issues/76
	e = ExtractArray(ExtractScalar(Float32))
	@test e([1,2,3]) == BagNode(
		ArrayNode([1. 2. 3.]),
		[1:3]
	)
	@test_throws MethodError e((1,2,3))
	@test e(nothing) == BagNode(
		ArrayNode(Matrix(fill(zero(Float32),1,0))),
		[0:-1]
	)
	@test_throws MethodError e(e)
end

@testset "ExtractVector" begin
	@testset "with uniontypes" begin
		sc1 = ExtractVector(5, true)
		e1 = sc1([1, 2, 2, 3, 4], store_input=false)
		e1s = sc1([1, 2, 2, 3, 4], store_input=true)
		e2s = sc1([1, 2, 2, 3], store_input=true)
		n1 = sc1(missing, store_input=false)
		n1s = sc1(missing, store_input=true)
		@test sc1.uniontypes == true
		@test e1.data == e1s.data
		@test e1.data isa Matrix{Union{Missing, Float32}}

		sc2 = ExtractVector{Int64}(5, true)
		e12234 = sc2([1, 2, 2, 3, 4], store_input=false)
		e12234s = sc2([1, 2, 2, 3, 4], store_input=true)
		en = sc2(nothing, store_input=false)
		ens = sc2(nothing, store_input=true)
		@test sc2.uniontypes == true
		@test e12234.data ≈ [1, 2, 2, 3, 4]
		@test e12234.data isa Matrix{Union{Missing, Int64}}
		@test e12234s.data ≈ e12234.data
		@test ens.data ≃ en.data
		@test en.data isa Matrix{Union{Missing, Int64}}
		@test en.data ≃ [missing missing missing missing missing]'

		sc3 = ExtractVector{Float32}(5, true)
		@test sc3.uniontypes == true
		# feature vector longer than expected
		sc122345 = sc3([1, 2, 2, 3, 4, 5], store_input=false)
		sc122345s = sc3([1, 2, 2, 3, 4, 5], store_input=true)
		sc12345 = sc3([1, 2, 3, 4, 5])
		@test sc122345.data == [1f0 2f0 2f0 3f0 4f0]'
		@test sc12345.data isa Matrix{Union{Missing, Float32}}
		e56 = sc3([5, 6], store_input=false)
		e56s = sc3([5, 6], store_input=true)
		@test e56.data ≃ [5f0 6f0 missing missing missing]'
		@test e56s.data ≃ e56.data
		@test sc122345s.data ≈ sc122345.data
		@test all(sc3(Dict(1=>2)).data .=== missing)

		@testset "extractempty" begin
			@test sc1(extractempty).data isa Matrix{Union{Missing, Float32}}
			@test nobs(sc1(extractempty).data) == 0
			@test nobs(sc1(extractempty)) == 0
			@test sc2(extractempty).data isa Matrix{Union{Missing, Int64}}
			@test nobs(sc2(extractempty).data) == 0
			@test nobs(sc2(extractempty)) == 0
		end

		@testset "store_input" begin
			@test e1s.metadata == [[1, 2, 2, 3, 4]]
			@test isnothing(e1.metadata)
			@test catobs(e1s, e1s).data == [1 1; 2 2; 2 2; 3 3; 4 4]
			@test catobs(e1s, e1s).metadata == [[1, 2, 2, 3, 4], [1, 2, 2, 3, 4]]
			@test catobs(e1s, e1s)[1].data == [1 2 2 3 4]'
			@test catobs(e1s, e1s)[1].metadata == [[1, 2, 2, 3, 4]]
			@test n1s.metadata ≃ [missing]
			@test catobs(e1s, n1s).metadata ≃ [[1, 2, 2, 3, 4], missing]
			@test e2s.data ≃ [1 2 2 3 missing]'
			@test e2s.metadata == [[1, 2, 2, 3]]
			@test e12234s.metadata == [[1, 2, 2, 3, 4]]
			@test ens.metadata == [nothing]
			@test sc122345s.metadata == [[1, 2, 2, 3, 4, 5]]
			@test isnothing(e56.metadata)
			@test e56s.metadata == [[5, 6]]
		end
	end
	@testset "without uniontypes" begin
		sc1 = ExtractVector(5, false)
		@test sc1.uniontypes == false
		e1 = sc1([1, 2, 2, 3, 4], store_input=false)
		e1s = sc1([1, 2, 2, 3, 4], store_input=true)
		@test e1.data == [1f0 2f0 2f0 3f0 4f0]'
		@test e1.data isa Matrix{Float32}
		@test e1.data == e1s.data

		sc2 = ExtractVector{Int64}(5, false)
		@test sc2.uniontypes == false
		@test all(sc2([1, 2, 2, 3, 4]).data .== [1, 2, 2, 3, 4])
		@test sc2([1, 2, 2, 3, 4]).data isa Array{Int64, 2}

		sc3 = ExtractVector{Float32}(5, false)
		@test sc3.uniontypes == false

		# feature vector longer than expected
		@test sc3([1, 2, 2, 3, 4, 5]).data == [1f0 2f0 2f0 3f0 4f0]'
		@test sc3([1, 2, 3, 4, 5]).data isa Matrix{Float32}

		@testset "extractempty" begin
			@test sc1(extractempty).data isa Matrix{Float32}
			@test nobs(sc1(extractempty).data) == 0
			@test sc2(extractempty).data isa Matrix{Int64}
			@test nobs(sc2(extractempty).data) == 0
			@test sc3(extractempty).data isa Matrix{Float32}
			@test nobs(sc3(extractempty).data) == 0
		end

		@testset "store_input" begin
			@test e1s.metadata == [[1, 2, 2, 3, 4]]
			@test catobs(e1s, e1s).data == [1 1; 2 2; 2 2; 3 3; 4 4]
			@test catobs(e1s, e1s).metadata == [[1, 2, 2, 3, 4], [1, 2, 2, 3, 4]]
			@test catobs(e1s, e1s)[1].data == [1 2 2 3 4]'
			@test catobs(e1s, e1s)[1].metadata == [[1, 2, 2, 3, 4]]
			@test isnothing(e1.metadata)
		end

		@testset "missing" begin
			@test_throws ErrorException sc1(nothing, store_input=false)
			@test_throws ErrorException sc1(nothing, store_input=true)
			@test_throws ErrorException sc2(nothing, store_input=false)
			@test_throws ErrorException sc2(nothing, store_input=true)
			@test_throws ErrorException sc1(missing)
			@test_throws ErrorException sc2(missing)
			@test_throws ErrorException sc1([5, 6])
			@test_throws ErrorException sc3([5, 6])
			@test_throws ErrorException sc3(Dict(1=>2))
		end
	end
end

@testset "ExtractDict" begin
	dict = Dict("a" => ExtractScalar(Float64,2,3),
				"b" => ExtractScalar(Float64),
				"c" => ExtractArray(ExtractScalar(Float64,2,3)))
	br = ExtractDict(dict)
	a1 = br(Dict("a" => 5, "b" => 7, "c" => [1,2,3,4]))
	a2 = br(Dict("a" => 5, "b" => 7))
	a3 = br(Dict("a" => 5, "c" => [1,2,3,4]))

	@test catobs(a1,a1)[:a].data ≈ [9 9]
	@test catobs(a1,a1)[:b].data ≈ [7 7]
	@test catobs(a1,a1)[:c].data.data ≈  [-3 0 3 6 -3 0 3 6]
	@test all(catobs(a1,a1)[:c].bags .== [1:4,5:8])

	@test catobs(a1,a2)[:a].data ≈ [9 9]
	@test catobs(a1,a2)[:b].data ≈ [7 7]
	@test catobs(a1,a2)[:c].data.data ≈ [-3 0 3 6]
	@test all(catobs(a1,a2)[:c].bags .== [1:4,0:-1])

	@test catobs(a2,a3)[:a].data ≈ [9 9]
	@test catobs(a2,a3)[:b].data ≃ [7.0 missing]
	@test catobs(a2,a3)[:c].data.data ≈ [-3 0 3 6]
	@test all(catobs(a2,a3)[:c].bags .== [0:-1,1:4])

	@test catobs(a1,a3)[:a].data ≈ [9 9]
	@test catobs(a1,a3)[:b].data ≃ [7.0 missing]
	@test catobs(a1,a3)[:c].data.data ≈ [-3 0 3 6 -3 0 3 6]
	@test all(catobs(a1,a3)[:c].bags .== [1:4,5:8])

	@test a1[:a].data ≈ [9]
	@test a1[:b].data ≈ [7]
	@test a2[:a].data ≈ [9]
	@test a2[:b].data ≈ [7]
	@test a3[:a].data ≈ [9]
	@test a3[:b].data ≃ hcat(missing)
	@test a1[:c].data.data ≈ [-3 0 3 6]
	@test all(a1[:c].bags .== [1:4])

	@test a3[:c].data.data ≈ [-3 0 3 6]
	@test all(a3[:c].bags .== [1:4])
	@test catobs(a3,a3)[:c].data.data ≈ [-3 0 3 6 -3 0 3 6]
	@test all(catobs(a3,a3)[:c].bags .== [1:4,5:8])

	@testset "extractempty" begin
		a4 = br(extractempty)
		@test nobs(a4) == 0
		@test nobs(a4[:a]) == 0
		@test a4[:a].data isa Matrix{Union{Missing, Float64}}
		@test nobs(a4[:b]) == 0
		@test a4[:b].data isa Matrix{Union{Missing, Float64}}
		@test nobs(a4[:c]) == 0
		@test nobs(a4[:c].data) == 0
		@test a4[:c].data.data isa Matrix{Union{Missing, Float64}}
	end
	@testset "Nested Missing Arrays" begin
		dict = Dict("a" => ExtractArray(ExtractScalar(Float32,2,3)),
			"b" => ExtractArray(ExtractScalar(Float32,2,3)))
		br = ExtractDict(dict)
		a1 = br(Dict("a" => [1,2,3], "b" => [1,2,3,4]), store_input=false)
		a2 = br(Dict("b" => [2,3,4]), store_input=false)
		a3 = br(Dict("a" => [2,3,4]), store_input=false)
		a4 = br(Dict{String, Any}(), store_input=false)
		a1s = br(Dict("a" => [1,2,3], "b" => [1,2,3,4]), store_input=true)
		a2s = br(Dict("b" => [2,3,4]), store_input=true)
		a3s = br(Dict("a" => [2,3,4]), store_input=true)
		a4s = br(Dict{String, Any}(), store_input=true)

		@test all(catobs(a1,a2).data[1].data.data .== [-3.0  0.0  3.0  6.0  0.0  3.0  6.0])
		@test all(catobs(a1,a2).data[1].bags .== [1:4, 5:7])
		@test all(catobs(a1,a2).data[2].data.data .== [-3.0  0.0  3.0])
		@test all(catobs(a1,a2).data[2].bags .== [1:3, 0:-1])

		@test all(catobs(a2,a3).data[1].data.data .== [0.0  3.0  6.0])
		@test all(catobs(a2,a3).data[1].bags .== [1:3, 0:-1])
		@test all(catobs(a2,a3).data[2].data.data .== [0 3 6])
		@test all(catobs(a2,a3).data[2].bags .== [0:-1, 1:3])

		@test all(catobs(a1,a4).data[1].data.data .== [-3.0  0.0  3.0  6.0])
		@test all(catobs(a1,a4).data[1].bags .== [1:4, 0:-1])
		@test all(catobs(a1,a4).data[2].data.data .== [-3.0  0.0  3.0])
		@test all(catobs(a1,a4).data[2].bags .== [1:3, 0:-1])

		@test all(a4.data[2].data.data isa Matrix{Union{Missing, Float32}})
		# todo: ověřit catobsování empty bagů s plnými bagy s metadaty jestli to funguje u všech možných extraktorů uvnitř bagu
		@test catobs(a1,a4).data[1].data.data == catobs(a1s,a4s).data[1].data.data
		@test catobs(a1,a4).data[1].bags == catobs(a1s,a4s).data[1].bags
		@test catobs(a1,a4).data[2].data.data == catobs(a1s,a4s).data[2].data.data
		@test catobs(a1,a4).data[2].bags == catobs(a1,a4).data[2].bags

		@test catobs(a1s,a4s).data[1].data.metadata == [1 2 3 4]
		@test catobs(a1s,a4s).data[2].data.metadata == [1 2 3]
		@test a1s.data[1].data.metadata == [1 2 3 4]
		@test a1s.data[2].data.metadata == [1 2 3]
		@test a4s.data[1].data.metadata == zeros(1,0)
		@test a4s.data[2].data.metadata == zeros(1,0)

		a4 = br(extractempty)
		@test nobs(a4) == 0
		@test nobs(a4[:a]) == 0
		@test nobs(a4[:a].data) == 0
		@test a4[:a].data.data isa Matrix{Union{Missing, Float32}}
		@test nobs(a4[:b]) == 0
		@test nobs(a4[:b].data) == 0
		@test a4[:b].data.data isa Matrix{Union{Missing, Float32}}
	end
end

@testset "ExtractCategorical" begin
	e = ExtractCategorical(["a","b"], true)
	ea = e("a", store_input=false)
	eb = e("b", store_input=false)
	ez = e("z", store_input=false)
	en = e(nothing, store_input=false)
	em = e(missing, store_input=false)
	eas = e("a", store_input=true)
	ebs = e("b", store_input=true)
	ezs = e("z", store_input=true)
	ens = e(nothing, store_input=true)
	@test e.uniontypes == true
	@test ea.data ≈ [1, 0, 0]
	@test eb.data ≈ [0, 1, 0]
	@test ez.data ≈ [0, 0, 1]
	@test en.data ≃ [missing missing missing]'
	@test em.data ≃ [missing missing missing]'
	@test typeof(ea.data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test typeof(en.data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test typeof(em.data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test e(extractempty).data isa MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test nobs(e(extractempty)) == 0

	@test e(["a", "b"]).data ≃ [missing missing missing]'
	@test mapreduce(e, catobs, ["a", "b"]).data ≈ [1 0; 0 1; 0 0]
	@test e(["a", missing]).data ≃ [missing missing missing]'
	@test mapreduce(e, catobs, ["a", missing]).data ≃ [true missing; false missing; false missing]
	@test e(["a", missing, "x"]).data ≃ [missing missing missing]'
	@test mapreduce(e, catobs, ["a", missing, "x"]).data ≃ [true missing false; false missing false; false missing true]
	@test typeof(e(["a", "b"]).data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test typeof(mapreduce(e, catobs, ["a", "b"]).data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test typeof(e(["a", "b", nothing]).data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
	@test typeof(mapreduce(e, catobs, ["a", "b", nothing]).data) == MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}

	@test catobs(ea, eb).data ≈ [1 0; 0 1; 0 0]
	@test reduce(catobs, [ea.data, eb.data]) ≈ [1 0; 0 1; 0 0]
	@test hcat(ea.data, eb.data) ≈ [1 0; 0 1; 0 0]
	@test e(Dict(1=>2)).data |> collect ≃ [missing missing missing]'

	@test catobs(eas, ebs).data == catobs(ea, eb).data
	@test reduce(catobs, [eas.data, ebs.data]) == reduce(catobs, [ea.data, eb.data])
	@test e(Dict(1=>2)).data ≃ [missing missing missing]'

	@test catobs(eas, ebs).metadata == ["a", "b"]
	@test e(Dict(1=>2), store_input=true).metadata == [Dict(1=>2)]

	@test nobs(ea) == 1
	@test nobs(eb) == 1
	@test nobs(ez) == 1
	@test nobs(en) == 1
	@test nobs(em) == 1
	@test nobs(em.data) == 1
	@test nobs(e([missing, nothing])) == 1
	@test nobs(mapreduce(e, catobs, [missing, nothing])) == 2
	@test nobs(e([missing, nothing, "a"])) == 1
	@test nobs(mapreduce(e, catobs, [missing, nothing, "a"])) == 3

	@test isnothing(ExtractCategorical([], true))
	e2 = ExtractCategorical(JsonGrinder.Entry(Dict("a"=>1,"c"=>1), 2), true)
	@test e2.uniontypes == true
	@test e2("a").data ≈ [1, 0, 0]
	@test e2("c").data ≈ [0, 1, 0]
	@test e2("b").data ≈ [0, 0, 1]
	@test e2(nothing).data ≃ [missing missing missing]'
	@test e2(missing).data ≃ [missing missing missing]'

	e3 = ExtractCategorical(JsonGrinder.Entry(Dict(1=>1,2=>1), 2), true)
	@test e3.uniontypes == true
	@test e3(1).data ≈ [1, 0, 0]
	@test e3(2).data ≈ [0, 1, 0]
	@test e3(4).data ≈ [0, 0, 1]
	@test e3(1.).data ≈ [1, 0, 0]
	@test e3(2.).data ≈ [0, 1, 0]
	@test e3(4.).data ≈ [0, 0, 1]
	@test e3([]).data ≃ [missing missing missing]'

	e4 = ExtractCategorical(JsonGrinder.Entry(Dict(1.0=>1,2.0=>1), 2), true)
	@test e4.uniontypes == true
	@test e4(1).data ≈ [1, 0, 0]
	@test e4(2).data ≈ [0, 1, 0]
	@test e4(4).data ≈ [0, 0, 1]
	@test e4(1.).data ≈ [1, 0, 0]
	@test e4(2.).data ≈ [0, 1, 0]
	@test e4(4.).data ≈ [0, 0, 1]
	@test e4([]).data ≃ [missing missing missing]'

	e = ExtractCategorical(["a","b"], false)
	@test e.uniontypes == false
	@test e("a").data ≈ [1, 0, 0]
	@test e("b").data ≈ [0, 1, 0]
	@test e("z").data ≈ [0, 0, 1]
	@test_throws ErrorException e(nothing)
	@test_throws ErrorException e(missing)
	@test typeof(e("a").data) == OneHotMatrix{UInt32, UInt32(3), Vector{UInt32}}
	@test e(extractempty).data isa OneHotMatrix{UInt32, UInt32(3), Vector{UInt32}}
	@test nobs(e(extractempty)) == 0

	@test mapreduce(e, catobs, ["a", "b"]).data ≈ [1 0; 0 1; 0 0]
	@test_throws ErrorException e(["a", "b"]).data
	@test_throws ErrorException e(["a", missing])
	@test_throws ErrorException e(["a", missing, "x"])
	@test_throws ErrorException typeof(e(["a", "b"]).data)
	@test typeof(mapreduce(e, catobs, ["a", "b"]).data) == OneHotMatrix{UInt32, UInt32(3), Vector{UInt32}}
	@test_throws ErrorException e(["a", "b", nothing])

	@test isnothing(ExtractCategorical([], false))
	e2 = ExtractCategorical(JsonGrinder.Entry(Dict("a"=>1,"c"=>1), 2), false)
	@test e2.uniontypes == false
	@test e2("a").data ≈ [1, 0, 0]
	@test e2("c").data ≈ [0, 1, 0]
	@test e2("b").data ≈ [0, 0, 1]
	@test_throws ErrorException e2(nothing)
	@test_throws ErrorException e2(missing)

	@test catobs(e("a"), e("b")).data ≈ [1 0; 0 1; 0 0]
	@test reduce(catobs, [e("a").data, e("b").data]) ≈ [1 0; 0 1; 0 0]
	@test hcat(e("a").data, e("b").data) ≈ [1 0; 0 1; 0 0]
	@test_throws ErrorException e(Dict(1=>2))

	@test nobs(e("a")) == 1
	@test nobs(e("b")) == 1
	@test nobs(e("z")) == 1
	@test_throws ErrorException nobs(e(nothing))
	@test_throws ErrorException nobs(e(missing))
	@test_throws ErrorException nobs(e([missing, nothing]))
	@test_throws ErrorException nobs(e([missing, nothing, "a"]))

	e3 = ExtractCategorical(JsonGrinder.Entry(Dict(1=>1,2=>1), 2), false)
	@test e3.uniontypes == false
	@test e3(1).data ≈ [1, 0, 0]
	@test e3(2).data ≈ [0, 1, 0]
	@test e3(4).data ≈ [0, 0, 1]
	@test e3(1.).data ≈ [1, 0, 0]
	@test e3(2.).data ≈ [0, 1, 0]
	@test e3(4.).data ≈ [0, 0, 1]
	@test_throws ErrorException isequal(e3([]).data, [missing missing missing]')

	e4 = ExtractCategorical(JsonGrinder.Entry(Dict(1.0=>1,2.0=>1), 2), false)
	@test e4.uniontypes == false
	@test e4(1).data ≈ [1, 0, 0]
	@test e4(2).data ≈ [0, 1, 0]
	@test e4(4).data ≈ [0, 0, 1]
	@test e4(1.).data ≈ [1, 0, 0]
	@test e4(2.).data ≈ [0, 1, 0]
	@test e4(4.).data ≈ [0, 0, 1]
	@test_throws ErrorException isequal(e4([]).data, [missing missing missing]')

	@testset "type conversions" begin
		j1 = JSON.parse("""{"a": "4"}""")
		j2 = JSON.parse("""{"a": "11.5"}""")
		j3 = JSON.parse("""{"a": 7}""")
		j4 = JSON.parse("""{"a": 4.5}""")

		sch = JsonGrinder.schema([j1,j2,j3,j4])

		ext = suggestextractor(sch)

		@test ext(Dict("a"=>4)) == ext(Dict("a"=>4.0))
		@test ext(Dict("a"=>4)) == ext(Dict("a"=>4f0))
		@test ext(Dict("a"=>4)) == ext(Dict("a"=>"4"))
		@test ext(Dict("a"=>4)) == ext(Dict("a"=>"4.0"))
		# todo: add here metadata test
	end
end

@testset "ExtractString" begin
	e = ExtractString(true)
	ehello = e("Hello", store_input=false)
	ehellos = e("Hello", store_input=true)
	@test ehello.data.S == ["Hello"]
	@test ehello.data == ehellos.data
	@test e(Symbol("Hello")).data.S == ["Hello"]
	@test e(["Hello", "world"]).data.S ≃ [missing]
	@test mapreduce(e, catobs, ["Hello", "world"]).data.S == ["Hello", "world"]

	@test e(missing).data.S ≃ [missing]
	@test e(nothing).data.S ≃ [missing]
	@test isequal(e(Dict(1=>2)), e(missing))
	@test ehello isa ArrayNode{NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}},Nothing}
	@test ehello.data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}
	@test e(missing).data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}
	@test e(nothing).data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}
	@test ehellos.metadata == ["Hello"]
	@test nobs(e(extractempty)) == 0

	e = ExtractString(false)
	@test e("Hello").data.S == ["Hello"]
	@test e(Symbol("Hello")).data.S == ["Hello"]
	@test_throws ErrorException e(["Hello", "world"]).data.S
	@test mapreduce(e, catobs, ["Hello", "world"]).data.S == ["Hello", "world"]
	@test_throws ErrorException e(missing)
	@test_throws ErrorException e(nothing)
	@test_throws ErrorException e(Dict(1=>2))
	@test e("Hello") isa ArrayNode{NGramMatrix{String,Vector{String},Int64},Nothing}
	@test e("Hello").data isa NGramMatrix{String,Vector{String},Int64}
end

@testset "ExtractKeyAsField" begin
	JsonGrinder.updatemaxkeys!(1000)
	js = [Dict(randstring(5) => rand()) for _ in 1:1000]
	sch = JsonGrinder.schema(js)
	ext = JsonGrinder.suggestextractor(sch, (;key_as_field = 500))
	@test ext isa JsonGrinder.ExtractKeyAsField

	b = ext(js[1])
	k = only(keys(js[1]))
	@test b.data[:item].data[1] ≈ (js[1][k] - ext.item.c) * ext.item.s
	@test b.data[:key].data.S[1] == k

	with_emptyismissing(true) do
		b = ext(nothing)
		@test nobs(b) == 1
		@test ismissing(b.data)
		b = ext(Dict())
		@test nobs(b) == 1
		@test ismissing(b.data)
	end
	with_emptyismissing(false) do
		b = ext(nothing)
		@test nobs(b) == 1
		@test nobs(b.data) == 0
		@test b.data[:key].data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}
		b = ext(Dict())
		@test nobs(b) == 1
		@test nobs(b.data) == 0
		@test b.data[:key].data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}
	end

	b = ext(extractempty)
	@test nobs(b) == 0
	@test nobs(b.data) == 0
	@test nobs(b.data[:item]) == 0
	@test b.data[:item].data isa Matrix{Float32}
	@test nobs(b.data[:key]) == 0
	@test b.data[:key].data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}

	js = [Dict(randstring(5) => Dict("a" => rand(), "b" => randstring(1))) for _ in 1:1000]
	sch = JsonGrinder.schema(js)
	ext = JsonGrinder.suggestextractor(sch, (;key_as_field = 500))

	b = ext(js[1])
	k = only(keys(js[1]))
	i = ext.item(js[1][k])
	@test b.data[:item][:a].data == i[:a].data
	@test b.data[:item][:b].data == i[:b].data
	@test b.data[:key].data.S[1] == k
	@test isnothing(b.data[:key].metadata)
	@test isnothing(b.data[:item][:a].metadata)
	@test isnothing(b.data[:item][:b].metadata)

	b = ext(nothing)
	@test nobs(b) == 1
	@test nobs(b.data) == 0
	b = ext(Dict())
	@test nobs(b) == 1
	@test nobs(b.data) == 0

	b = ext(js[1], store_input=true)
	k = only(keys(js[1]))
	@test b.data[:key].metadata == [k]
	@test b.data[:item][:a].metadata == fill(js[1][k]["a"],1,1)
	@test b.data[:item][:b].metadata == [js[1][k]["b"]]

	b = ext(Dict(), store_input=true)
	@test b.metadata == [Dict()]
	@test isnothing(b.data[:key].metadata)
	@test isnothing(b.data[:item].metadata)

	b = ext(nothing, store_input=true)
	@test b.metadata == [nothing]
	@test isnothing(b.data.data.key.metadata)
	@test isnothing(b.data.data.item.metadata)

	j1 = JSON.parse("""{"a": 1}""")
	j2 = JSON.parse("""{"b": 2.5}""")
	j3 = JSON.parse("""{"c": 3.1}""")
	j4 = JSON.parse("""{"d": 5}""")
	j5 = JSON.parse("""{"e": 4.5}""")
	j6 = JSON.parse("""{"f": 5}""")
	sch = JsonGrinder.schema([j1, j2, j3, j4, j5, j6])
	e = JsonGrinder.key_as_field(sch, NamedTuple(), path = "")
	@test e isa JsonGrinder.ExtractKeyAsField

	e2 = JsonGrinder.key_as_field(sch, NamedTuple(), path = "")
	@test hash(e) === hash(e2)
	@test e == e2
end

@testset "equals and hash test" begin
	other1 = Dict(
		"a" => ExtractArray(ExtractScalar(Float64,2,3)),
		"b" => ExtractArray(ExtractScalar(Float64,2,3)),
		"c" => ExtractCategorical(["a","b"]),
		"d" => ExtractVector(4),
	)
	br1 = ExtractDict(other1)
	other11 = Dict(
		"a" => ExtractArray(br1),
		"b" => ExtractScalar(Float64,2,3),
	)
	br11 = ExtractDict(other11)

	other2 = Dict(
		"a" => ExtractArray(ExtractScalar(Float64,2,3)),
		"b" => ExtractArray(ExtractScalar(Float64,2,3)),
		"c" => ExtractCategorical(["a","b"]),
		"d" => ExtractVector(4),
	)
	br2 = ExtractDict(other2)
	other22 = Dict("a" => ExtractArray(br2), "b" => ExtractScalar(Float64,2,3))
	br22 = ExtractDict(other22)

	@test hash(br11) === hash(br22)
	@test hash(br11) !== hash(br1)
	@test hash(br11) !== hash(br2)
	@test hash(br1) === hash(br2)

	@test br11 == br22
	@test br11 != br1
	@test br11 != br2
	@test br1 == br2
end
# todo: add tests for settings of suggestextractor!!!
# todo: add tests for single keyed dict, that model generated from such jsons behaves as we want
@testset "Extractor skip empty lists" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}], "b": [], "c": [[], []], "d": 1, "e": {"a": 1}, "f": {"a": []}}""")
	j2 = JSON.parse("""{"a": [{"a":3},{"b":4}], "b": [], "c": [[]], "d": 1, "e": {"a": 2}, "f": {"a": []}}""")
	j3 = JSON.parse("""{"a": [{"a":1},{"b":3}], "b": [], "c": [[], [], [], []], "d": 2, "e": {"a": 3}, "f": {"a": []}}""")
	j4 = JSON.parse("""{"a": [{"a":2},{"b":4}], "b": [], "c": [], "d": 4, "e": {"a": 3}, "f": {"a": []}}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch)
	@test ext[:a] isa ExtractArray
	@test isnothing(ext[:b])
	@test isnothing(ext[:c])
	@test ext[:d] isa ExtractCategorical
	@test ext[:e] isa ExtractDict
	@test isnothing(ext[:f])
end

@testset "Extractor of numbers as strings" begin
	j1 = JSON.parse("""{"a": "1", "b": "a", "c": "1.1", "d": 1.1, "e": "1.2", "f": 1, "g": "1"}""")
	j2 = JSON.parse("""{"a": "2", "b": "b", "c": "2",   "d": 2, "e": "1.3", "f": 2, "g": "2"}""")
	j3 = JSON.parse("""{"a": "3", "b": "c", "c": "2.3", "d": 2.3, "e": "1.4", "f": 3, "g": "3"}""")
	j4 = JSON.parse("""{"a": "4", "b": "c", "c": "5",   "d": 5, "e": "1.4", "f": 3, "g": "4"}""")
	j5_200 = [JSON.parse("""{"a": "5", "b": "$("c"^i)", "c": "5", "d": $(5+i/10), "e": "1.4", "f": 3, "g": "$i"}""") for i in 5:200]
	sch = JsonGrinder.schema([j1, j2, j3, j4, j5_200...])
	ext = suggestextractor(sch)

	@test ext[:a] isa ExtractCategorical
	@test ext[:b] isa ExtractString
	@test ext[:c] isa ExtractCategorical
	@test ext[:d] isa ExtractScalar{Float32}
	@test ext[:e] isa ExtractCategorical
	@test ext[:f] isa ExtractCategorical
	@test ext[:g] isa ExtractScalar{Float32}

	ext_j1 = ext(j1)
	ext_j2 = ext(j2)
	ext_j3 = ext(j3)
	ext_j4 = ext(j4)

	@test eltype(ext_j1[:a].data) <: Bool
	@test eltype(ext_j1[:b].data) <: Int64
	@test eltype(ext_j1[:c].data) <: Bool
	@test eltype(ext_j1[:d].data) <: Float32
	@test eltype(ext_j1[:e].data) <: Bool
	@test eltype(ext_j1[:f].data) <: Bool
	@test eltype(ext_j1[:g].data) <: Float32

	@test eltype(ext_j2[:a].data) <: Bool
	@test eltype(ext_j2[:b].data) <: Int64
	@test eltype(ext_j2[:c].data) <: Bool
	@test eltype(ext_j2[:d].data) <: Float32
	@test eltype(ext_j2[:e].data) <: Bool
	@test eltype(ext_j2[:f].data) <: Bool
	@test eltype(ext_j2[:g].data) <: Float32

	@test eltype(ext_j3[:a].data) <: Bool
	@test eltype(ext_j3[:b].data) <: Int64
	@test eltype(ext_j3[:c].data) <: Bool
	@test eltype(ext_j3[:d].data) <: Float32
	@test eltype(ext_j3[:e].data) <: Bool
	@test eltype(ext_j3[:f].data) <: Bool
	@test eltype(ext_j3[:g].data) <: Float32

	@test eltype(ext_j4[:a].data) <: Bool
	@test eltype(ext_j4[:b].data) <: Int64
	@test eltype(ext_j4[:c].data) <: Bool
	@test eltype(ext_j4[:d].data) <: Float32
	@test eltype(ext_j4[:e].data) <: Bool
	@test eltype(ext_j4[:f].data) <: Bool
	@test eltype(ext_j4[:g].data) <: Float32

	@test ext_j1["U"].data ≈ [0]
	@test ext_j2["U"].data ≈ [(2-1.1)/23.9]
	@test ext_j3["U"].data ≈ [(2.3-1.1)/23.9]
	@test ext_j4["U"].data ≈ [(5-1.1)/23.9]

	ext_j1s = ext(j1, store_input=true)
	ext_j2s = ext(j2, store_input=true)
	ext_j3s = ext(j3, store_input=true)
	ext_j4s = ext(j4, store_input=true)

	@test ext_j1["U"].data == ext_j1s["U"].data
	@test ext_j2["U"].data == ext_j2s["U"].data
	@test ext_j3["U"].data == ext_j3s["U"].data
	@test ext_j4["U"].data == ext_j4s["U"].data

	@test Mill.metadata.(values(ext_j1s.data)) == ([1], ["a"], ["1"], fill(1.1,1,1), ["1.2"], ["1.1"], fill("1",1,1))
	@test Mill.metadata.(values(ext_j2s.data)) == ([2], ["b"], ["2"], fill(2,1,1), ["1.3"], ["2"], fill("2",1,1))
	@test Mill.metadata.(values(ext_j3s.data)) == ([3], ["c"], ["3"], fill(2.3,1,1), ["1.4"], ["2.3"], fill("3",1,1))
	@test Mill.metadata.(values(ext_j4s.data)) == ([3], ["c"], ["4"], fill(5,1,1), ["1.4"], ["5"], fill("4",1,1))

	m = reflectinmodel(sch, ext)
	@test buf_printtree(m) == """
	ProductModel ↦ ArrayModel(Dense(52, 10)) \t# 2 arrays, 530 params, 2.148 KiB
	  ├── f: ArrayModel(Dense(4, 10)) \t# 2 arrays, 50 params, 280 bytes
	  ├── b: ArrayModel(Dense(2053, 10)) \t# 2 arrays, 20_540 params, 80.312 KiB
	  ├── a: ArrayModel(Dense(6, 10)) \t# 2 arrays, 70 params, 360 bytes
	  ├── d: ArrayModel(identity)
	  ├── e: ArrayModel(Dense(4, 10)) \t# 2 arrays, 50 params, 280 bytes
	  ├── c: ArrayModel(Dense(5, 10)) \t# 2 arrays, 60 params, 320 bytes
	  └── g: ArrayModel(identity)
	"""
end

# todo: add separate tests for reflectinmodel(sch, ext) to test behavior with various missing and non-missing stuff in schema
# e.g. empty string, missing keys, irregular schema and MultiRepresentation

@testset "Suggest feature vector extraction" begin
	j1 = JSON.parse("""{"a": "1", "b": [1, 2, 3], "c": [1, 2, 3]}""")
	j2 = JSON.parse("""{"a": "2", "b": [2, 2, 3], "c": [1, 2, 3, 4]}""")
	j3 = JSON.parse("""{"a": "3", "b": [3, 2, 3], "c": [1, 2, 3, 4, 5]}""")
	j4 = JSON.parse("""{"a": "4", "b": [2, 3, 4], "c": [1, 2, 3]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch, testing_settings)

	@test ext[:a] isa ExtractScalar{Float32}
	@test ext[:b] isa ExtractVector
	@test ext[:b].n == 3
	@test ext[:c] isa ExtractArray{ExtractCategorical{Number,UInt32}}

	@test buf_printtree(ext, trav=true) ==
	"""
	Dict [""]
	  ├── a: Float32 ["E"]
	  ├── b: FeatureVector with 3 items ["U"]
	  └── c: Array of ["k"]
	           └── Categorical d = 6 ["s"]
	"""

	ext_j1 = ext(j1, store_input=false)
	ext_j2 = ext(j2, store_input=false)
	ext_j3 = ext(j3, store_input=false)
	ext_j4 = ext(j4, store_input=false)

	ext_j1s = ext(j1, store_input=true)
	ext_j2s = ext(j2, store_input=true)
	ext_j3s = ext(j3, store_input=true)
	ext_j4s = ext(j4, store_input=true)

	@test ext_j1["E"].data ≈ [0]
	@test ext_j2["E"].data ≈ [1/3]
	@test ext_j3["E"].data ≈ [2/3]
	@test ext_j4["E"].data ≈ [1]

	@test ext_j1["U"].data ≈ [1, 2, 3]
	@test ext_j2["U"].data ≈ [2, 2, 3]
	@test ext_j3["U"].data ≈ [3, 2, 3]
	@test ext_j4["U"].data ≈ [2, 3, 4]

	@test ext_j1["s"].data ≈ [
	1  0  0
 	0  1  0
 	0  0  1
 	0  0  0
 	0  0  0
 	0  0  0]
	@test ext_j2["s"].data ≈ [
	1  0  0  0
 	0  1  0  0
 	0  0  1  0
 	0  0  0  1
 	0  0  0  0
 	0  0  0  0
	]
	@test ext_j3["s"].data ≈ [
	1  0  0  0  0
	0  1  0  0  0
	0  0  1  0  0
	0  0  0  1  0
	0  0  0  0  1
	0  0  0  0  0
	]
	@test ext_j4["s"].data ≈ [
	1  0  0
 	0  1  0
 	0  0  1
 	0  0  0
 	0  0  0
 	0  0  0
	]

	for (a,b) in zip([ext_j1, ext_j2, ext_j3, ext_j4], [ext_j1s, ext_j2s, ext_j3s, ext_j4s])
		@test a["E"].data == b["E"].data
		@test a["U"].data == b["U"].data
		@test a["s"].data == b["s"].data
	end

	@test ext_j1s["U"].metadata == [[1, 2, 3]]
	@test ext_j2s["U"].metadata == [[2, 2, 3]]
	@test ext_j3s["U"].metadata == [[3, 2, 3]]
	@test ext_j4s["U"].metadata == [[2, 3, 4]]

	@test ext_j1s["s"].metadata == [1, 2, 3]
	@test ext_j2s["s"].metadata == [1, 2, 3, 4]
	@test ext_j3s["s"].metadata == [1, 2, 3, 4, 5]
	@test ext_j4s["s"].metadata == [1, 2, 3]
end

@testset "Suggest complex" begin
	JsonGrinder.updatemaxkeys!(1000)
	js = [Dict("a" => rand(), "b" => Dict(randstring(5) => rand()), "c"=>[rand(), rand()], "d"=>[rand() for i in 1:rand(1:10)]) for _ in 1:1000]
	sch = JsonGrinder.schema(js)
	ext = JsonGrinder.suggestextractor(sch, (;key_as_field = 500))

	@test ext[:a] isa ExtractScalar{Float32}
	@test ext[:b] isa JsonGrinder.ExtractKeyAsField
	@test ext[:b].key isa ExtractString
	@test ext[:b].item isa ExtractScalar{Float32}
	@test ext[:c] isa ExtractVector
	@test ext[:c].n == 2
	@test ext[:d] isa ExtractArray
end

@testset "Extractor typed missing bags" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [1, 2]}""")
	j3 = JSON.parse("""{"a": [1, 2, 3]}""")
	j4 = JSON.parse("""{"a": [1, 2, 4, 5]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch, testing_settings)

	ext_j1 = ext(j1)
	ext_j2 = ext(j2)
	ext_j3 = ext(j3)
	ext_j4 = ext(j4)

	@test ext_j1[:a].data.data isa Matrix{Float32}
	@test ext_j2[:a].data.data isa Matrix{Float32}
	@test ext_j3[:a].data.data isa Matrix{Float32}
	@test ext_j4[:a].data.data isa Matrix{Float32}
end

@testset "testing irregular extractor" begin
	j1 = JSON.parse("""{"a": 4}""")
	j2 = JSON.parse("""{"a": { "a": "hello", "b":[5,6]}}""")
	j3 = JSON.parse("""{"a": [1, 2, 3, 4]}""")
	j4 = JSON.parse("""{"a": 5}""")
	j5 = JSON.parse("""{"a": 3}""")

	sch = schema([j1,j2,j3,j4,j5])
	ext = suggestextractor(sch, testing_settings)
	a = ext(j1)
	@test a[:a][:e1].data[1] == 0.5
	@test ismissing(a[:a][:e2].data[:a].data.S[1])
	@test nobs(a[:a][:e3]) == 1
	@test nobs(a[:a][:e3].data) == 1
end

@testset "Mixed scalar extraction" begin
    # todo: add test for extraction of numeric scalars from string and also missings, if it's all matrix and catobsing properly
	j1 = JSON.parse("""{"a": "1"}""")
	j2 = JSON.parse("""{"a": 4}""")
	j3 = JSON.parse("""{"a": "3.1"}""")
	j4 = JSON.parse("""{"a": 2.5}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	sch_hash = hash(sch)
	ext = suggestextractor(sch, testing_settings)
	@test sch_hash == hash(sch)	# testing that my entry parkour does not modify schema
	@test nchildren(ext[:a]) == 1

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	@test e1["k"].data ≈ [0]
	@test e2["k"].data ≈ [1]
	@test e3["k"].data ≈ [0.7]
	@test e4["k"].data ≈ [0.5]
end

@testset "merging numeric scalars" begin
	j1 = JSON.parse("""{"a": "1", "b":1}""")
	j2 = JSON.parse("""{"a": 4, "b":2}""")
	j3 = JSON.parse("""{"a": "3.1", "b":3}""")
	j4 = JSON.parse("""{"b":4}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch, testing_settings)
	@test nchildren(ext[:a]) == 1

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	@test e1["M"].data ≈ [0]
	@test e2["M"].data ≈ [1]
	@test e3["M"].data ≈ [0.7]
	@test e4["M"].data ≃ fill(missing,1,1)
	@test e1["U"].data ≈ [0]
	@test e2["U"].data ≈ [1/3]
	@test e3["U"].data ≈ [2/3]
	@test e4["U"].data ≈ [1]
end

@testset "Mixed scalar extraction with other types" begin
	j1 = JSON.parse("""{"a": "1"}""")
	j2 = JSON.parse("""{"a": 2.5}""")
	j3 = JSON.parse("""{"a": "3.1"}""")
	j4 = JSON.parse("""{"a": 5}""")
	j5 = JSON.parse("""{"a": "4.5"}""")
	j6 = JSON.parse("""{"a": [1, 2, 3, 4, 5]}""")
	j7 = JSON.parse("""{"a": {"Sylvanas is the worst warchief ever": "yes"}}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4, j5, j6, j7])
	ext = suggestextractor(sch, testing_settings)

	@test buf_printtree(sch) ==
	"""
	[Dict] \t# updated = 7
	  └── a: [MultiEntry] \t# updated = 7
	           ├── 1: [Scalar - String], 3 unique values \t# updated = 3
	           ├── 2: [Scalar - Float64,Int64], 2 unique values \t# updated = 2
	           ├── 3: [List] \t# updated = 1
	           │        └── [Scalar - Int64], 5 unique values \t# updated = 5
	           └── 4: [Dict] \t# updated = 1
	                    └── Sylvanas is the worst warchief ever: [Scalar - String], 1 unique values \t# updated = 1
	"""

	@test buf_printtree(ext, trav=true) ==
    """
	Dict [""]
	  └── a: MultiRepresentation ["U"]
	           ├── e1: FeatureVector with 5 items ["c"]
	           ├── e2: Dict ["k"]
	           │         └── Sylvanas is the worst warchief ever: String ["o"]
	           └── e3: Float32 ["s"]
	"""

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	e5 = ext(j5)
	@test e1["s"].data ≈ [0]
	@test e2["s"].data ≈ [0.375]
	@test e3["s"].data ≈ [0.525]
	@test e4["s"].data ≈ [1.0]
	@test e5["s"].data ≈ [0.875]
	@test buf_printtree(e1) ==
	"""
	ProductNode \t# 1 obs, 48 bytes
	  └── a: ProductNode \t# 1 obs, 48 bytes
	           ├── e1: ArrayNode(5×1 Array with Union{Missing, Float32} elements) \t# 1 obs, 73 bytes
	           ├── e2: ProductNode \t# 1 obs, 32 bytes
	           │         └── Sylvanas is the worst warchief ever: ArrayNode(2053×1 NGramMatrix with Union{Missing, Int64} elements) \t# 1 obs, 112 bytes
	           └── e3: ArrayNode(1×1 Array with Union{Missing, Float32} elements) \t# 1 obs, 53 bytes
	"""
	@test ext[:a][1] == ext["c"]
end

@testset "mixing numeric and non-numeric strings" begin
	j1 = JSON.parse("""{"a": "hello"}""")
	j2 = JSON.parse("""{"a": "4"}""")
	j3 = JSON.parse("""{"a": 5}""")
	j4 = JSON.parse("""{"a": 3.0}""")
	j5 = JSON.parse("""{"a": [1, 2, 3, 4, 5]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4, j5])
	ext = suggestextractor(sch, testing_settings)

	@test buf_printtree(sch) ==
	"""
	[Dict] \t# updated = 5
	  └── a: [MultiEntry] \t# updated = 5
	           ├── 1: [Scalar - String], 2 unique values \t# updated = 2
	           ├── 2: [Scalar - Float64,Int64], 2 unique values \t# updated = 2
	           └── 3: [List] \t# updated = 1
	                    └── [Scalar - Int64], 5 unique values \t# updated = 5
	"""

	@test buf_printtree(ext) ==
	"""
	Dict
	  └── a: MultiRepresentation
	           ├── e1: String
	           ├── e2: Float32
	           └── e3: FeatureVector with 5 items
	"""

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	e5 = ext(j5)
	@test all(e1["k"].data .=== [missing])
	@test e2["k"].data ≈ [0.5]
	@test e3["k"].data ≈ [1.0]
	@test e4["k"].data ≈ [0.0]
	@test all(e5["k"].data .=== [missing])

	@test hash(ext) !== hash(suggestextractor(JsonGrinder.schema([j1, j2, j4, j5])))
end

@testset "empty string and substring" begin
	a, b, c = split("a b ", " ")
	d = "d"
	e = ""
	f = 5.2
	@test a isa SubString{String}
	@test b isa SubString{String}
	@test c isa SubString{String}
	@test d isa String
	@test e isa String
	@test f isa AbstractFloat
	@test c == ""
	ext = JsonGrinder.extractscalar(AbstractString)
	@test SparseMatrixCSC(ext(c).data) == SparseMatrixCSC(ext(e).data)
	@test all(ext(f).data.S .=== ext(nothing).data.S)

	@test ext(f, store_input=true) != ext(nothing, store_input=true)
	@test ext(f, store_input=true).data ≃ ext(nothing, store_input=true).data
	@test ext(a, store_input=true).metadata == ["a"]
	@test ext(b, store_input=true).metadata == ["b"]
	@test ext(c, store_input=true).metadata == [""]
	@test ext(d, store_input=true).metadata == ["d"]
	@test ext(e, store_input=true).metadata == [""]
	@test ext(f, store_input=true).metadata == [5.2]
	@test ext(nothing, store_input=true).metadata == [nothing]
	#todo: add tests for all extractors for store_input=true
end

@testset "AuxiliaryExtractor" begin
	e2 = ExtractCategorical(["a","b"])
	e = AuxiliaryExtractor(e2, (ext, sample; store_input=false)->ext(String(sample); store_input))

	@test e("b") == e("bbb"[1:1])
	@test e("b").data ≈ [0, 1, 0]
	@test e("bbb"[1:1]).data ≈ [0, 1, 0]

    @test buf_printtree(e) ==
	"""
	Auxiliary extractor with
	  └── Categorical d = 3
	"""
end

@testset "Not skipping single dict key" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [1, 2]}""")
	j3 = JSON.parse("""{"a": [1, 2, 3]}""")
	j4 = JSON.parse("""{"a": [1, 2, 4, 5]}""")
	j5 = JSON.parse("""{"a": [1, 2, 3, 6, 7]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4, j5])
	ext = suggestextractor(sch, testing_settings)

	@test buf_printtree(ext) ==
	"""
	Dict
	  └── a: Array of
	           └── Float32
	"""

	ext_j2 = ext(j2)
	@test buf_printtree(ext_j2) ==
    """
	ProductNode \t# 1 obs, 16 bytes
	  └── a: BagNode \t# 1 obs, 80 bytes
	           └── ArrayNode(1×2 Array with Float32 elements) \t# 2 obs, 56 bytes
	"""

	m = reflectinmodel(sch, ext)
	@test buf_printtree(m) ==
    """
	ProductModel ↦ ArrayModel(identity)
	  └── a: BagModel ↦ BagCount([SegmentedMean(1); SegmentedMax(1)]) ↦ ArrayModel(Dense(3, 10)) \t# 4 arrays, 42 params, 328 bytes
	           └── ArrayModel(identity)
	"""
end

@testset "default_scalar_extractor" begin

end

@testset "type stability of extraction missings" begin
	# todo: make tests for categorical, scalar, string, vector
	j1 = JSON.parse("""{"a": 1}""")
	j2 = JSON.parse("""{"b": "a"}""")
	j3 = JSON.parse("""{"a": 3, "b": "b"}""")
	j4 = JSON.parse("""{"a": 4, "b": "c"}""")
	j5 = JSON.parse("""{"a": 5, "b": "d"}""")

	sch = schema([j1,j2,j3,j4,j5])
	ext = suggestextractor(sch, testing_settings)
	m = reflectinmodel(sch, ext)

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e12 = catobs(e1, e2)
	e23 = catobs(e2, e3)
	e13 = catobs(e1, e3)
	typeof(e1)
	typeof(e2)
	typeof(e3)
	typeof(e12)
	typeof(e23)
	typeof(e13)
	@test typeof(e1) == typeof(e2)
	@test typeof(e1) == typeof(e3)
	@test typeof(e1) == typeof(e12)
	@test typeof(e1) == typeof(e23)
	@test typeof(e1) == typeof(e13)
end


text:
# Box first
@test begin
    p = Point(0.5, 0.5)
    s = Box(Point(0, 0), Point(1, 1))

    p ∈ s
end

@test begin
    p = Point(0.0, 0.5)
    s = Box(Point(0, 0), Point(1, 1))

    p ∈ s
end

@test begin
    p = Point(-1E-10, 0.5)
    s = Box(Point(0, 0), Point(1, 1))

    p ∉ s
end

@test begin
    s = Box(Point(0, 0, 0), Point(1, 1, 1))
    p = Point(0.5, 0.5, 0.5)

    p ∈ s
end

@test begin
    s = Box(Point(0, 0, 0), Point(1, 1, 1))
    p = Point(1, 1, 1)

    p ∈ s
end

@test begin
    s = Box(Point(0, 0, 0), Point(1, 1, 1))
    p = Point(1, 0.5, 1)

    p ∈ s
end

@test begin
    s = Box(Point(0, 0, 0), Point(1, 1, 1))
    p = Point(3, 0.5, 0.5)

    p ∉ s
end

max_stars_repo_name: MiroK/fractional-dim