text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

const ContactElements3D = Union{Tri3,Tri6,Quad4,Quad8,Quad9}

function create_orthogonal_basis(n)
    I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    k = argmax([norm(cross(n,I[:,k])) for k in 1:3])
    t1 = cross(n, I[:,k])/norm(cross(n, I[:,k]))
    t2 = cross(n, t1)
    return t1, t2
end

""" Create rotation matrix Q for element nodes rotating quantities to nt coordinaet system. """
function create_rotation_matrix(element::Element{Tri3}, time::Float64)
    n = element("normal", time)
    t11, t21 = create_orthogonal_basis(n[1])
    t12, t22 = create_orthogonal_basis(n[2])
    t13, t23 = create_orthogonal_basis(n[3])
    Q1_ = [n[1] t11 t21]
    Q2_ = [n[2] t12 t22]
    Q3_ = [n[3] t13 t23]
    Z = zeros(3, 3)
    Q = [
        Q1_ Z Z
        Z Q2_ Z
        Z Z Q3_]
    return Q
end

function create_rotation_matrix(element::Element{Quad4}, time::Float64)
    n = element("normal", time)
    t11, t21 = create_orthogonal_basis(n[1])
    t12, t22 = create_orthogonal_basis(n[2])
    t13, t23 = create_orthogonal_basis(n[3])
    t14, t24 = create_orthogonal_basis(n[4])
    Q1_ = [n[1] t11 t21]
    Q2_ = [n[2] t12 t22]
    Q3_ = [n[3] t13 t23]
    Q4_ = [n[4] t14 t24]
    Z = zeros(3, 3)
    Q = [
        Q1_ Z Z Z
        Z Q2_ Z Z
        Z Z Q3_ Z
        Z Z Z Q4_]
    return Q
end

function create_rotation_matrix(element::Element{Tri6}, time::Float64)
    n = element("normal", time)
    t11, t21 = create_orthogonal_basis(n[1])
    t12, t22 = create_orthogonal_basis(n[2])
    t13, t23 = create_orthogonal_basis(n[3])
    t14, t24 = create_orthogonal_basis(n[4])
    t15, t25 = create_orthogonal_basis(n[5])
    t16, t26 = create_orthogonal_basis(n[6])
    Q1_ = [n[1] t11 t21]
    Q2_ = [n[2] t12 t22]
    Q3_ = [n[3] t13 t23]
    Q4_ = [n[4] t14 t24]
    Q5_ = [n[5] t15 t25]
    Q6_ = [n[6] t16 t26]
    Z = zeros(3, 3)
    Q = [
        Q1_ Z Z Z Z Z
        Z Q2_ Z Z Z Z
        Z Z Q3_ Z Z Z
        Z Z Z Q4_ Z Z
        Z Z Z Z Q5_ Z
        Z Z Z Z Z Q6_]
    return Q
end

""" Create a contact segmentation between one slave element and list of master elements.

Returns
-------

Vector with tuples: (master_element, polygon_clip_vertices, polygon_clip_centroid, polygon_clip_area)
"""
function create_contact_segmentation(slave_element, master_elements, x0, n0, time::Float64; deformed=false)
    result = []
    x1 = slave_element("geometry", time)
    if deformed
        x1 = map(+, x1, slave_element("displacement", time))
    end
    S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x1]
    for master_element in master_elements
        x2 = master_element("geometry", time)
        if deformed
            x2 = map(+, x2, master_element("displacement", time))
        end
        M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2]
        P = get_polygon_clip(S, M, n0)
        length(P) < 3 && continue # no clipping or shared edge (no volume)
        check_orientation!(P, n0)
        N_P = length(P)
        P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
        if isapprox(P_area, 0.0)
            error("Polygon P has zero area")
        end
        C0 = calculate_centroid(P)
        push!(result, (master_element, P, C0, P_area))
    end
    return result
end

function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time::Float64)

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)

    nsl = length(slave_element)
    X1 = slave_element("geometry", time)
    u1 = slave_element("displacement", time)
    x1 = map(+, X1, u1)
    n1 = slave_element("normal", time)
    la = slave_element("lambda", time)

    Q3 = create_rotation_matrix(slave_element, time)

    # project slave nodes to auxiliary plane (x0, Q)
    xi = get_mean_xi(slave_element)
    N = vec(get_basis(slave_element, xi, time))
    x0 = interpolate(N, X1)
    n0 = interpolate(N, n1)

    # create contact segmentation
    segmentation = create_contact_segmentation(slave_element, slave_element("master elements", time), x0, n0, time)

    if length(segmentation) == 0 # no overlapping surface in slave and maters
        return
    end

    Ae = Matrix{Float64}(I, nsl, nsl)

    if problem.properties.dual_basis # construct dual basis

        De = zeros(nsl, nsl)
        Me = zeros(nsl, nsl)

        # loop all polygons
        for (master_element, P, C0, P_area) in segmentation

            # loop integration cells
            for cell in get_cells(P, C0)
                virtual_element = Element(Tri3, Int[])
                update!(virtual_element, "geometry", tuple(cell...))
                for ip in get_integration_points(virtual_element, 3)
                    detJ = virtual_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ
                    x_gauss = virtual_element("geometry", ip, time)
                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
                    N1 = slave_element(xi_s, time)
                    De += w*Matrix(Diagonal(vec(N1)))
                    Me += w*N1'*N1
                end # integration points done

            end # integration cells done

        end # master elements done

        Ae = De*inv(Me)

    end

    # loop all polygons
    for (master_element, P, C0, P_area) in segmentation

        nm = length(master_element)
        X2 = master_element("geometry", time)
        u2 = master_element("displacement", time)
        x2 = map(+, X2, u2)

        De = zeros(nsl, nsl)
        Me = zeros(nsl, nm)
        ce = zeros(field_dim*nsl)
        ge = zeros(field_dim*nsl)

        # loop integration cells
        for cell in get_cells(P, C0)
            virtual_element = Element(Tri3, Int[])
            update!(virtual_element, "geometry", tuple(cell...))
            # loop integration point of integration cell
            for ip in get_integration_points(virtual_element, 3)

                # project gauss point from auxiliary plane to master and slave element
                x_gauss = virtual_element("geometry", ip, time)
                xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
                xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)

                detJ = virtual_element(ip, time, Val{:detJ})
                w = ip.weight*detJ

                # add contributions
                N1 = vec(get_basis(slave_element, xi_s, time))
                N2 = vec(get_basis(master_element, xi_m, time))
                Phi = Ae*N1
                De += w*Phi*N1'
                Me += w*Phi*N2'

                x_s = interpolate(N1, map(+,X1,u1))
                x_m = interpolate(N2, map(+,X2,u2))
                ge += w*vec((x_m-x_s)*Phi')

            end # integration points done

        end # integration cells done

        # add contribution to contact virtual work
        sdofs = get_gdofs(problem, slave_element)
        mdofs = get_gdofs(problem, master_element)
        nsldofs = length(sdofs)
        nmdofs = length(mdofs)
        D3 = zeros(nsldofs, nsldofs)
        M3 = zeros(nsldofs, nmdofs)
        for i=1:field_dim
            D3[i:field_dim:end, i:field_dim:end] += De
            M3[i:field_dim:end, i:field_dim:end] += Me
        end

        add!(problem.assembly.C1, sdofs, sdofs, D3)
        add!(problem.assembly.C1, sdofs, mdofs, -M3)
        add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
        add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
        add!(problem.assembly.g, sdofs, Q3'*ge)

    end # master elements done

end

""" Assemble quadratic surface element to contact problem. """
function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time::Float64)

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)

    alp = props.alpha

    if alp != 0.0
        T = [
                1.0 0.0 0.0 0.0 0.0 0.0
                0.0 1.0 0.0 0.0 0.0 0.0
                0.0 0.0 1.0 0.0 0.0 0.0
                alp alp 0.0 1.0-2*alp 0.0 0.0
                0.0 alp alp 0.0 1.0-2*alp 0.0
                alp 0.0 alp 0.0 0.0 1.0-2*alp
            ]
    else
        T = Matrix(1.0*I, 6, 6)
    end

    nsl = length(slave_element)
    Xs = slave_element("geometry", time)
    n1 = slave_element("normal", time)

    Q3 = create_rotation_matrix(slave_element, time)

    Ae = Matrix(1.0*I, nsl, nsl)

    if problem.properties.dual_basis # construct dual basis

        nsl = length(slave_element)
        De = zeros(nsl, nsl)
        Me = zeros(nsl, nsl)

        for sub_slave_element in split_quadratic_element(slave_element, time)

            slave_element_nodes = get_connectivity(sub_slave_element)
            nsl = length(sub_slave_element)

            X1 = sub_slave_element("geometry", time)
            #u1 = sub_slave_element("displacement", time)
            #x1 = X1 + u1
            n1 = sub_slave_element("normal", time)
            #la = sub_slave_element("lambda", time)

            # create auxiliary plane
            xi = get_mean_xi(sub_slave_element)
            N = vec(get_basis(sub_slave_element, xi, time))
            x0 = interpolate(N, X1)
            n0 = interpolate(N, n1)

            # project slave nodes to auxiliary plane
            S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]

            # 3. loop all master elements
            for master_element in slave_element("master elements", time)

                Xm = master_element("geometry", time)

                if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
                    continue
                end

                # split master element to linear sub-elements and loop
                for sub_master_element in split_quadratic_element(master_element, time)

                    master_element_nodes = get_connectivity(sub_master_element)
                    nm = length(sub_master_element)
                    X2 = sub_master_element("geometry", time)
                    #u2 = sub_master_element("displacement", time)
                    #x2 = X2 + u2

                    # 3.1 project master nodes to auxiliary plane and create polygon clipping
                    M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
                    P = get_polygon_clip(S, M, n0)
                    length(P) < 3 && continue # no clipping or shared edge (no volume)
                    check_orientation!(P, n0)

                    N_P = length(P)
                    P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
                    if isapprox(P_area, 0.0)
                        error("Polygon P has zero area")
                    end

                    C0 = calculate_centroid(P)

                    # 4. loop integration cells
                    for cell in get_cells(P, C0)
                        virtual_element = Element(Tri3, Int[])
                        update!(virtual_element, "geometry", tuple(cell...))
                        for ip in get_integration_points(virtual_element, 3)
                            detJ = virtual_element(ip, time, Val{:detJ})
                            w = ip.weight*detJ
                            x_gauss = virtual_element("geometry", ip, time)
                            xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
                            N1 = vec(slave_element(xi_s, time)*T)
                            De += w*Matrix(Diagonal(N1))
                            Me += w*N1*N1'
                        end # integration points done

                    end # integration cells done

                end # sub master elements done

            end # master elements done

        end # sub slave elements done

        Ae = De*inv(Me)

    end

    # split slave element to linear sub-elements and loop
    for sub_slave_element in split_quadratic_element(slave_element, time)

        slave_element_nodes = get_connectivity(sub_slave_element)
        nsl = length(sub_slave_element)
        X1 = sub_slave_element("geometry", time)
        n1 = sub_slave_element("normal", time)

        # create auxiliary plane
        xi = get_mean_xi(sub_slave_element)
        N = vec(get_basis(sub_slave_element, xi, time))
        x0 = interpolate(N, X1)
        n0 = interpolate(N, n1)

        # project slave nodes to auxiliary plane
        S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]

        # 3. loop all master elements
        for master_element in slave_element("master elements", time)

            Xm = master_element("geometry", time)

            if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
                continue
            end

            # split master element to linear sub-elements and loop
            for sub_master_element in split_quadratic_element(master_element, time)

                master_element_nodes = get_connectivity(sub_master_element)
                nm = length(master_element)
                X2 = sub_master_element("geometry", time)
                #u2 = master_element("displacement", time)
                #x2 = X2 + u2

                # 3.1 project master nodes to auxiliary plane and create polygon clipping
                M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
                P = get_polygon_clip(S, M, n0)
                length(P) < 3 && continue # no clipping or shared edge (no volume)
                check_orientation!(P, n0)

                N_P = length(P)
                P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
                if isapprox(P_area, 0.0)
                    error("Polygon P has zero area")
                end

                C0 = calculate_centroid(P)

                # integration is done in quadratic elements
                nsl = length(slave_element)
                nm = length(master_element)
                De = zeros(nsl, nsl)
                Me = zeros(nsl, nm)
                ge = zeros(field_dim*nsl)

                # 4. loop integration cells
                for cell in get_cells(P, C0)
                    virtual_element = Element(Tri3, Int[])
                    update!(virtual_element, "geometry", tuple(cell...))

                    # 5. loop integration point of integration cell
                    for ip in get_integration_points(virtual_element, 3)

                        # project gauss point from auxiliary plane to master and slave element
                        x_gauss = virtual_element("geometry", ip, time)
                        xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
                        xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time)

                        detJ = virtual_element(ip, time, Val{:detJ})
                        w = ip.weight*detJ

                        # add contributions
                        N1 = vec(get_basis(slave_element, xi_s, time)*T)
                        N2 = vec(get_basis(master_element, xi_m, time))
                        Phi = Ae*N1

                        De += w*Phi*N1'
                        Me += w*Phi*N2'

                        us = slave_element("displacement", time)
                        um = master_element("displacement", time)
                        xs = interpolate(N1, map(+,Xs,us))
                        xm = interpolate(N2, map(+,Xs,um))
                        ge += w*vec((xm-xs)*Phi')

                    end # integration points done

                end # integration cells done

                # 6. add contribution to contact virtual work
                sdofs = get_gdofs(problem, slave_element)
                mdofs = get_gdofs(problem, master_element)
                nsldofs = length(sdofs)
                nmdofs = length(mdofs)
                D3 = zeros(nsldofs, nsldofs)
                M3 = zeros(nsldofs, nmdofs)
                for i=1:field_dim
                    D3[i:field_dim:end, i:field_dim:end] += De
                    M3[i:field_dim:end, i:field_dim:end] += Me
                end

                add!(problem.assembly.C1, sdofs, sdofs, D3)
                add!(problem.assembly.C1, sdofs, mdofs, -M3)
                add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
                add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
                add!(problem.assembly.g, sdofs, Q3'*ge)

            end # sub master elements done

        end # master elements done

    end # sub slave elements done

end

"""
Frictionless 3d small sliding contact.

problem
time
dimension
finite_sliding
friction
use_forwarddiff
"""
function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::Type{Val{false}}, ::Type{Val{false}}, ::Type{Val{false}})

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)

    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
    normals = calculate_normals(slave_elements, time, Val{2};
                                rotate_normals=props.rotate_normals)
    update!(slave_elements, "normal", time => normals)

    # 2. loop all slave elements
    for slave_element in slave_elements
        assemble!(problem, slave_element, time)
    end # slave elements done, contact virtual work ready

    S = sort(collect(keys(normals))) # slave element nodes
    weighted_gap = Dict{Int64, Vector{Float64}}()
    contact_pressure = Dict{Int64, Vector{Float64}}()
    complementarity_condition = Dict{Int64, Vector{Float64}}()
    is_active = Dict{Int64, Int}()
    is_inactive = Dict{Int64, Int}()
    is_slip = Dict{Int64, Int}()
    is_stick = Dict{Int64, Int}()

    la = problem.assembly.la

    # FIXME: for matrix operations, we need to know the dimensions of the
    # final matrices
    ndofs = 0
    ndofs = max(ndofs, size(problem.assembly.K, 2))
    ndofs = max(ndofs, size(problem.assembly.C1, 2))
    ndofs = max(ndofs, size(problem.assembly.C2, 2))
    ndofs = max(ndofs, size(problem.assembly.D, 2))
    ndofs = max(ndofs, size(problem.assembly.g, 2))
    ndofs = max(ndofs, size(problem.assembly.c, 2))

    C1 = sparse(problem.assembly.C1, ndofs, ndofs)
    C2 = sparse(problem.assembly.C2, ndofs, ndofs)
    D = sparse(problem.assembly.D, ndofs, ndofs)
    g = Vector(problem.assembly.g, ndofs)
    c = Vector(problem.assembly.c, ndofs)

    maxdim = maximum(size(C1))
    if problem.properties.alpha != 0.0
        alp = problem.properties.alpha
        Te = [
                1.0 0.0 0.0 0.0 0.0 0.0
                0.0 1.0 0.0 0.0 0.0 0.0
                0.0 0.0 1.0 0.0 0.0 0.0
                alp alp 0.0 1.0-2*alp 0.0 0.0
                0.0 alp alp 0.0 1.0-2*alp 0.0
                alp 0.0 alp 0.0 0.0 1.0-2*alp
            ]
        invTe = [
            1.0 0.0 0.0 0.0 0.0 0.0
            0.0 1.0 0.0 0.0 0.0 0.0
            0.0 0.0 1.0 0.0 0.0 0.0
            -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
            0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
            -alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
        ]
        # construct global transformation matrices T and invT
        T = SparseMatrixCOO()
        invT = SparseMatrixCOO()
        for element in slave_elements
            dofs = get_gdofs(problem, element)
            for i=1:field_dim
                ldofs = dofs[i:field_dim:end]
                add!(T, ldofs, ldofs, Te)
                add!(invT, ldofs, ldofs, invTe)
            end
        end
        T = sparse(T, maxdim, maxdim, (a, b) -> b)
        invT = sparse(invT, maxdim, maxdim, (a, b) -> b)
        # fill diagonal
        d = ones(size(T, 1))
        d[get_nonzero_rows(T)] .= 0.0
        T += sparse(Diagonal(d))
        invT += sparse(Diagonal(d))
        #invT2 = sparse(inv(full(T)))
        #@info("invT == invT2? ", invT == invT2)
        #maxabsdiff = maximum(abs(invT - invT2))
        #@info("max diff = $maxabsdiff")
        C1 = C1*invT
        C2 = C2*invT
    end

    tol = problem.properties.drop_tolerance
    SparseArrays.droptol!(C1, tol)
    SparseArrays.droptol!(C2, tol)

    for j in S
        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
        weighted_gap[j] = g[dofs]
    end

    state = problem.properties.contact_state_in_first_iteration
    if problem.properties.iteration == 1
        @info("First contact iteration, initial contact state = $state")

        if state == :AUTO
            avg_gap = mean([weighted_gap[j][1] for j in S])
            std_gap = std([weighted_gap[j][1] for j in S])
            if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
                state = :ACTIVE
            else
                state = :UNKNOWN
            end
            @info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
        end

    end

    # active / inactive node detection
    for j in S
        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
        weighted_gap[j] = g[dofs]
        if length(la) != 0
            normal = normals[j]
            tangent1, tangent2 = create_orthogonal_basis(normal)
            p = dot(normal, la[dofs])
            t1 = dot(tangent1, la[dofs])
            t2 = dot(tangent2, la[dofs])
            contact_pressure[j] = [p, t1, t2]
        else
            contact_pressure[j] = [0.0, 0.0, 0.0]
        end
        complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]

        if complementarity_condition[j][1] > 0.0
            is_inactive[j] = 0
            is_active[j] = 1
            is_slip[j] = 1
            is_stick[j] = 0
        else
            is_inactive[j] = 1
            is_active[j] = 0
            is_slip[j] = 0
            is_stick[j] = 0
        end
    end

    if (problem.properties.iteration == 1) && (state == :ACTIVE)
        for j in S
            is_inactive[j] = 0
            is_active[j] = 1
            is_slip[j] = 1
            is_stick[j] = 0
        end
    end

    if (problem.properties.iteration == 1) && (state == :INACTIVE)
        for j in S
            is_inactive[j] = 1
            is_active[j] = 0
            is_slip[j] = 0
            is_stick[j] = 0
        end
    end

    @info("# | active | stick | slip | gap | pres | comp")
    for j in S
        str1 = "$j | $(is_active[j]) | $(is_stick[j]) | $(is_slip[j]) | "
        str2 = "$(round(weighted_gap[j][1]; digits=3)) | $(round(contact_pressure[j][1]; digits=3)) | $(round(complementarity_condition[j][1]; digits=3))"
        @info(str1 * str2)
    end

    for j in S
        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
        tdofs = [3*(j-1)+2, 3*(j-1)+3]
        if is_inactive[j] == 1
            # remove inactive nodes from assembly
            C1[dofs,:] .= 0.0
            C2[dofs,:] .= 0.0
            D[dofs,:] .= 0.0
            g[dofs,:] .= 0.0
        elseif (is_active[j] == 1) && (is_slip[j] == 1)
            # constitutive modelling in tangent direction, frictionless contact
            C2[tdofs,:] .= 0.0
            g[tdofs] .= 0.0
            normal = normals[j]
            tangent1, tangent2 = create_orthogonal_basis(normal)
            D[tdofs[1], dofs] .= tangent1
            D[tdofs[2], dofs] .= tangent2
        end
    end

    problem.assembly.C1 = C1
    problem.assembly.C2 = C2
    problem.assembly.D = D
    problem.assembly.g = g

end

function postprocess!(problem::Problem{Contact}, time::Float64, ::Type{Val{Symbol("contact pressure")}})
    n = problem("normal", time)
    la = problem("lambda", time)
    node_ids = keys(n)
    cp = Dict(nid => dot(n[nid], la[nid]) for nid in node_ids)
    # FIXME: have to define zero contact pressure & lambda to master elements
    # elements because interface.elements = [slave_elements; master_elements]
    for nid in keys(la)
        if !haskey(cp, nid)
            cp[nid] = 0.0
        end
    end
    update!(problem, "contact pressure", time => cp)
end


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

mutable struct LoadBalancerBackendAddressPoolListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{BackendAddressPool} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function LoadBalancerBackendAddressPoolListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(LoadBalancerBackendAddressPoolListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(LoadBalancerBackendAddressPoolListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type LoadBalancerBackendAddressPoolListResult

const _property_map_LoadBalancerBackendAddressPoolListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_LoadBalancerBackendAddressPoolListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{BackendAddressPool}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ LoadBalancerBackendAddressPoolListResult }) = collect(keys(_property_map_LoadBalancerBackendAddressPoolListResult))
Swagger.property_type(::Type{ LoadBalancerBackendAddressPoolListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancerBackendAddressPoolListResult[name]))}
Swagger.field_name(::Type{ LoadBalancerBackendAddressPoolListResult }, property_name::Symbol) =  _property_map_LoadBalancerBackendAddressPoolListResult[property_name]

function check_required(o::LoadBalancerBackendAddressPoolListResult)
    true
end

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

ext: jl


text:
module test_stdlib_libdl

using Test
using Libdl

libllvm_paths = filter(Libdl.dllist()) do lib
    occursin("LLVM", basename(lib))
end

filename = basename(first(libllvm_paths))
if Sys.iswindows()
    @test filename == "LLVM.dll"
else
    (name, ext) = splitext(filename)
    @test startswith(name, "libLLVM")
end

end # module test_stdlib_libdl


text:
struct Id
    id::Int64
end

using AutoHashEquals
@auto_hash_equals mutable struct Term
    head::Symbol
    args::Array{Id}
end

mutable struct EClass
    nodes::Array{Term}
    parents::Array{Tuple{Term,Id}}
end

using DataStructures

mutable struct EGraph
    unionfind::IntDisjointSets
    memo::Dict{Term,Id} # int32 UInt32?
    classes::Dict{Id,EClass} # Use array?
    dirty_unions::Array{Id}
end

# Build an empty EGraph
function EGraph()
    EGraph(IntDisjointSets(0), Dict(), Dict(), [])
end

struct PatVar
    id::Symbol
end

struct PatTerm
    head::Symbol
    args::Array{Union{PatTerm, PatVar}}
end

Pattern = Union{PatTerm,PatVar}


text:
using Statistics, StatsBase
using StructArrays
using Recombinase
using Recombinase: compute_summary, aroundindex, discrete,
    prediction, density
using Test
using IndexedTables
using OnlineStatsBase

@testset "discrete" begin
    x = [1, 2, 3, 1, 2, 3]
    y = [0.3, 0.1, 0.3, 0.4, 0.2, 0.1]
    across = [1, 1, 1, 2, 2, 2]
    res = compute_summary(
        discrete(prediction)(min_nobs = 1),
        across,
        (x, y),
        stats = Mean()
    )
    xcol, ycol = fieldarrays(res)
    @test xcol == [1, 2, 3]
    @test map(Recombinase._first, ycol) ≈ [0.35, 0.15, 0.2]
    res = compute_summary(
        discrete(density),
        across,
        (x,),
        stats = Mean()
    )
    xcol, ycol = fieldarrays(res)
    @test xcol == [1, 2, 3]
    @test map(Recombinase._first, ycol) ≈ [1, 1, 1]./3
end

# example function to test initstats and fitvecmany!
function fitvec(stats, iter, ranges)
    start = iterate(iter)
    start === nothing && error("Nothing to fit!")
    val, state = start
    init = Recombinase.initstats(stats, Recombinase.merge_tups(axes(val), ranges))
    Recombinase.fitvecmany!(init, Iterators.rest(iter, state))
    StructArray(((nobs = nobs(el), value = value(el)) for el in init);
        unwrap = t -> t <: Union{Tuple, NamedTuple})
end

@testset "timeseries" begin
    v1 = rand(1000) # day 1
    v2 = rand(50) # day 2
    traces = [v1, v1, v1, v2, v2, v2]
    ts = [10, 501, 733, 1, 20, 30]
    trims = [7:13, 498:504, 730:736, 1:4, 17:23, 27:33]
    stats = Series(mean = Mean(), variance = Variance())
    s = fitvec(stats, (aroundindex(trace, t) for (trace, t) in zip(traces, ts)), -5:5);
    @test axes(s) == (-5:5,)
    @test s[-3].nobs == 4
    @test s[-3].value isa NamedTuple{(:mean, :variance)}
    s = fitvec(stats, (aroundindex(trace, t, trim) for (trace, t, trim) in zip(traces, ts, trims)), -5:5);
    @test axes(s) == (-5:5,)
    @test s[-3].nobs == 4
    @test s[-3].value isa NamedTuple{(:mean, :variance)}
    @test s[-4].nobs == 0

    stats = Mean()
    s = fitvec(stats, (aroundindex(trace, t) for (trace, t) in zip(traces, ts)), -5:5);
    @test axes(s) == (-5:5,)
    @test s[-3].nobs == 4
    @test s[-3].value isa Float64

    g() = aroundindex(rand(4, 4, 4), (2, 2), (1:3, 1:3))
    @inferred g()
    @test axes(g()) == (-1:1, -1:1, 1:4)
end


text:
"""

trueG2_test() は，分割表のマス目のうちに 0 があるような場合でも，正確に lambda=0.0 にした場合の，正しい答えを返す。

"""

using Distributions
function trueG2_test(x; correct=false)
    # セルに 0 があっても正しい答えを出す
    # correct=true で，連続性の補正も行うことができる
    ln(n) = sum(n .== 0 ? 0 : n .* log(n) for n in vcat(n...))
    nrows, ncols = size(x)
    n = sum(x) # 全サンプルサイズ
    n1 = sum(x, dims=2) # 行和
    n2 = sum(x, dims=1) # 列和
    G2 = 2*(ln(x) - ln(n1) - ln(n2) + ln(n)) # G 統計量
    correct && (G2 /= 1 + (n * sum(1 ./ n1) - 1) * (n * sum(1 ./ n2) - 1) / (6n * nrows * ncols)) # 連続性の補正
    df = (nrows - 1) * (ncols - 1) # G の自由度
    p = ccdf(Chisq(df), G2)
    name = correct ? "corrected G-sq." : "G-sq."
    println("$name = $G2,  df = $df,  p value = $p")
end

"""
使用法

julia> z = [4 5 2 0
            0 7 6 1
            1 0 3 1]
3×4 Matrix{Int64}:
 4  5  2  0
 0  7  6  1
 1  0  3  1

julia> trueG2_test(z)
G-sq. = 15.364591286599591,  df = 6,  p value = 0.017602886503051453

また，ピアソンの χ2
検定では 2×2 分割表の場合にしか連続性の補正は行えないが，G2 検定は分割表の大きさに関わらず，連続性の補正を行うことができる。

上に示した関数 trueG2_test() で correct=true を指定すればよい。

julia> trueG2_test(z, correct=true) # 連続性の補正
corrected G-sq. = 13.776490649307341,  df = 6,  p value = 0.032235015595581236

"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
using Weilianqi
using Test

include("lovelidalu.jl")


text:
"""
Set of functions for generating signals.
"""

function gendata_RW(process_noise,
                    measurement_noise,
                    mean_state0,
                    precision_state0;
                    time_horizon=100)
    "Generate data according to a random walk"

    # Dimensionality
    d = size(process_noise, 1)

    # Preallocate
    y = zeros(time_horizon, d)
    x = zeros(time_horizon+1, d)

    # Initialize state
    x[1, :] = randn(d)./sqrt(precision_state0) .+ mean_state0

    for t = 1:time_horizon

        # Evolve state
        x[t+1, :] = randn(d).*sqrt(process_noise) + x[t, :]

        # Observe
        y[t, :] = randn(d).*sqrt(measurement_noise) + x[t+1, :]

    end
    return y, x
end

function gendata_LGDS(transition_matrix,
                      emission_matrix,
                      process_noise,
                      measurement_noise,
                      mean_state0,
                      precision_state0;
                      time_horizon=100)
    "Generate data according to a linear Gaussian dynamical system"

    # Dimensionality
    d = size(transition_matrix, 1)

    # Preallocate
    y = zeros(time_horizon, d)
    x = zeros(time_horizon+1, d)

    # Initialize state
    x[1, :] = randn(d)./sqrt(precision_state0) .+ mean_state0

    for t = 1:time_horizon

        # Evolve state
        x[t+1, :] = randn(d).*sqrt(process_noise) + transition_matrix*x[t, :]

        # Observe
        y[t, :] = randn(d).*sqrt(measurement_noise) + emission_matrix*x[t+1, :]

    end
    return y, x
end

function gendata_LGDSBO(transition_matrix,
                       emission_matrix,
                       process_noise,
                       measurement_noise,
                       mean_state0,
                       precision_state0;
                       time_horizon=100)
    "Generate data according to a linear Gaussian dynamical system with binary observations"

    # Dimensionality
    d = size(transition_matrix, 1)

    # Preallocate
    y = zeros(time_horizon, d)
    x = zeros(time_horizon+1, d)

    # Initialize state
    x[1, :] = randn(d)./sqrt(precision_state0) .+ mean_state0

    for t = 1:time_horizon

        # Evolve state
        x[t+1, :] = randn(d).*sqrt(process_noise) + transition_matrix*x[t, :]

        # Observe
        y[t, :] = round.(1 ./(1 .+ exp.(-x[t,:]))).*2 .-1

    end

    return y, x
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module Loops

import ..DebugMsg
DebugMsg.init()

using CompilerTools
import Base.show
export DomLoops, Loop

PRETEST = 1
POSTTEST = 2
INFINITE = 3

function classify(head, back_edge, bbs)
    if CompilerTools.CFGs.classifyBlock(bbs[head]) == CompilerTools.CFGs.BLOCK_GOTOIFNOT
        return PRETEST
    elseif CompilerTools.CFGs.classifyBlock(bbs[back_edge]) == CompilerTools.CFGs.BLOCK_GOTOIFNOT
        return POSTTEST
    else
        return INFINITE
    end
end

"""
A type to hold information about a loop.
A loop has a "head" that dominates all the other blocks in the loop.
A loop has a back_edge which is a block that has "head" as one of its successors.
It also contains "members" which is a set of basic block labels of all the basic blocks that are a part of this loop.
"""
type Loop
    head :: Int           # The head basic block ID of the loop.
    back_edge :: Int      # The end of the loop that branches back to the head.
    members :: Set{Int}   # All the basic blocks that are part of the loop.
    exits :: Set{Int}     # Blocks outside the loop that are successors to some blocks inside the loop.
    blocks_that_exit :: Set{Int} # The blocks inside the loop that have succesors outside the loop.
    ltype :: Int

    function Loop(h :: Int, b :: Int, m :: Set{Int}, ex :: Set{Int}, be :: Set{Int}, lt :: Int)
        new(h, b, m, ex, be, lt)
    end
end

"""
A type that holds information about which basic blocks dominate which other blocks.
It also contains an array "loops" of all the loops discovered within the function.
The same basic block may occur as a member in multiple loop entries if those loops are nested.
This module doesn't currently help identify these nesting levels but by looking at loop members it is easy enough to figure out.
"""
type DomLoops
    dom_dict :: Dict{Int,Set}
    loops    :: Array{Loop,1}
end

function getPreLoopBlock(l :: Loop, bl :: CompilerTools.CFGs.CFG)
    head_bb = bl.basic_blocks[l.head]
    return setdiff(map(x -> x.label, head_bb.preds), l.members)
end

function getPostLoopBlock(l :: Loop, bl :: CompilerTools.CFGs.CFG)
    back_edge_bb = bl.basic_blocks[l.back_edge]
    return setdiff(map(x -> x.label, back_edge_bb.succs), l.members)
end

function insertNewBlockBeforeLoop(l :: Loop, bl :: CompilerTools.CFGs.CFG, stmts :: Array{Any,1})
    @dprintln(3, "insertNewBlockBeforeLoop l = ", l, " bl = ", bl, " stmts = ", stmts)
    new_bb, new_goto_stmt = CompilerTools.CFGs.insertBefore(bl, l.head, true, l.back_edge)
    @dprintln(3, "new_bb = ", new_bb, " new_goto_stmt = ", new_goto_stmt.expr)
    for stmt in stmts
        CompilerTools.CFGs.addStatementToEndOfBlock(bl, new_bb, stmt)
    end
    if new_goto_stmt != nothing
        push!(new_bb.statements, new_goto_stmt)
    end
    return new_bb
end

"""
Takes a DomLoops object containing loop information about the function.
Returns true if the given basic block label "bb" is in some loop in the function.
"""
function isInLoop(dl :: DomLoops, bb :: Int)
    for i in dl.loops
        if in(bb, i.members)
            return true
        end
    end

    return false
end

"""
Finds those computations within a loop that are iteration invariant.
Takes as input:
   l - the Loop to find invariants in
   udinfo - UDChains (use-definition chains) for the basic blocks of the function.
   bl - LivenessAnalysis.BlockLiveness with liveness information for variables within the function.
"""
function findLoopInvariants(l :: Loop,
                            udinfo :: Dict{CompilerTools.LivenessAnalysis.BasicBlock,CompilerTools.UDChains.UDInfo},
                            bl :: CompilerTools.LivenessAnalysis.BlockLiveness)
    all_uses = Set{Symbol}()

    @dprintln(3,l.members)
    @dprintln(3,collect(keys(bl.basic_blocks)))
    @dprintln(3,bl)

    # Accumulate a set of all variables used in the loop.
    for i in l.members
        @dprintln(3,"i = ", i)
        bb = bl.basic_blocks[i]  # get the basic block from its index
        all_uses = union(all_uses, bb.use)
    end

    @dprintln(3, "all_uses = ", all_uses)

    # Eliminate all variable usages in "all_uses" that have a definition in the loop.
    for i in l.members
        @dprintln(3,"i = ", i)
        bb = bl.basic_blocks[i]  # get the basic block from its index
        bb_ud = udinfo[bb]
        @dprintln(3, bb_ud)
        # Get all the inputs to this block.
        for j in bb_ud.live_in
            @dprintln(3, "symbol in block live_in = ", j[1])
            # If one of those inputs is used in this basic block of the loop.
            if in(j[1], all_uses)
                @dprintln(3, "symbol in all_uses")
                for k in j[2]
                    if k == nothing
                        continue
                    end
                    @dprintln(3, "this block = ", k.label)
                    if in(k.label, l.members)
                        temp_set = Set{Symbol}()
                        push!(temp_set, j[1])
                        @dprintln(3, "this block belongs to members all_uses = ", all_uses, " type = ", typeof(all_uses), " type j[1] = ", typeof(j[1]), " temp set = ", temp_set)
                        all_uses = setdiff(all_uses, temp_set)
                        break
                    end
                end
            end
        end
    end

    # At this point, we have an initial seed for loop-invariant computation in a set of variables that are used in the loop but never defined anywhere in the loop.

    changed = true
    while(changed)
        # While we keep finding more loop invariant computations, keep iterating.
        changed = false
        @dprintln(0,"starting changed loop")
        non_invariant = Set()
        invariant     = Set()

        # Try to find calculations that depend only on constants or other loop invariant calculatings and put those in "invariant".
        for i in l.members
            bb = bl.basic_blocks[i]  # get the basic block from its index
            for stmt in bb.statements
                if !isempty(stmt.def)
                    # Compute symbols used in the statement that aren't invariant.
                    non_invariants = setdiff(stmt.use, all_uses)
                    @dprintln(0,"Def statement use = ", stmt.use, " all_uses = ", all_uses, " non_invariants = ", non_invariants, " stmt = ", stmt)
                    if isempty(non_invariants)
                        # def could be invariants if they aren't non_invariant elsewhere.
                        for def in stmt.def
                            if !in(def, non_invariant) && !in(def, all_uses)
                                push!(invariant, def)
                            end
                        end
                    else
                        # def isn't invariant so remove from invariant if a prior statement that was invariant was added.
                        for def in stmt.def
                            if in(def, invariant)
                                delete!(invariant, def)
                            end
                            push!(non_invariant, def)
                        end
                    end
                    @dprintln(0,"non_invariant set = ", non_invariant, " invariant set = ", invariant)
                end
            end
        end

        # If we found additional loop invariant calculations then add them to the set of invariant calculation in all_uses and indicate we found a change
        # so that the loop will continue searching for more.
        if !isempty(invariant)
            @dprintln(0,"Found some new invariants = ", invariant)
            changed = true
            all_uses = union(all_uses, invariant)
        end
    end

    return all_uses
end

"""
Add to the "members" of the loop being accumulated given "cur_bb" which is known to be a member of the loop.
"""
function flm_internal(cur_bb, members, bbs)
    # If we haven't previously added the current basic block to the member set.
    if !in(cur_bb, members)
        bb = bbs[cur_bb]
        # Add the BasicBlock object for cur_bb to the members set.
        push!(members, cur_bb)

        # Add all predecessors of the current block as loop members.
        for pred in bb.preds
            flm_internal(pred.label, members, bbs)
        end
    end
    members
end

"""
Find all the members of the loop as specified by the "head" basic block and the "back_edge" basic block.
Also takes the dictionary of labels to basic blocks.
Start with just the head of the loop as a member and then starts recursing with the back_edge using flm_internal.
"""
function findLoopMembers(head, back_edge, bbs)
    members = Set(head)
    flm_internal(back_edge, members, bbs)
end

"""
Scans all loop members to find their successors which themselves are not loop members.
This computes the set of exit blocks.
"""
function findExitBlocks(members, bl)
    blocks_that_exit = Set{Int}()
    exits = Set{Int}()

    for lm in members
        mbb = bl.basic_blocks[lm]
        for s in mbb.succs
            if !in(s.label, members)
                push!(exits, s.label)
                push!(blocks_that_exit, lm)
                assert(CompilerTools.CFGs.classifyBlock(mbb) == CompilerTools.CFGs.BLOCK_GOTOIFNOT)
            end
        end
    end

    return (exits, blocks_that_exit)
end

"""
Find the loops in a CFGs.CFG in "bl".
"""
function compute_dom_loops(bl :: CompilerTools.CFGs.CFG)
    # Get the depth-first numering for the CFG.
    bbs_df_order = bl.depth_first_numbering
    # Get the number of basic blocks.
    num_bb = length(bl.basic_blocks)
    assert(num_bb == length(bbs_df_order))

    dom_dict = CompilerTools.CFGs.compute_dominators(bl)

    loops = Loop[]

    # Find the loop in the function.
    for i = 1:num_bb
        # For each basic block ...
        bb_index = bbs_df_order[i]
        bb = bl.basic_blocks[bb_index]
        succ_array = collect(bb.succs)
        for j in succ_array
            local succ_id = j.label
            # If one of the block's successors dominates the current basic block then we've found a loop.
            if in(succ_id, dom_dict[bb_index])
                # Identify the members of the loop.
                members = findLoopMembers(succ_id, bb_index, bl.basic_blocks)
                @dprintln(3,"loop members = ", members, " type = ", typeof(members))
                new_loop = Loop(succ_id, bb_index, members, findExitBlocks(members, bl)..., classify(succ_id, bb_index, bl.basic_blocks))
                @dprintln(3,"new_loop = ", new_loop, " type = ", typeof(new_loop))
                # Store the Loop information in the array of loops to be returned.
                push!(loops, new_loop)
            end
        end
    end

    DomLoops(dom_dict, loops)  # Return the dominator information and the array of loops found.
end

end


text:
#create a path ending at source taking each vertex only one time.
function makecirclepath(G::DictGraph, S = Set{Int}(), T = DictGraph(directed = true); source = 0, ending = 0)
    if length(G.relation) - length(S) <= 1
        l = getweight(G, ending, source)
        addedge!(T, ending, source; weight = l)
        return T, source, ending
    end

    ℒ = minweight(G, ending; restrict = S)
    addedge!(T, ending, ℒ.u; weight = ℒ.w)
    S = S ∪ ending
    makecirclepath(G, S, T; ending = ℒ.u)
end

function makepath(D::AbstractDict, source::Int, target::Int, G::DictGraph = DictGraph(); pred = 2)
    if s == v
        nothing
    elseif d[v][pred] == -1
        error("no path from $s to $v exists.")
    else
        makepath(d, s, d[v][2], G)
        addedge!(G, d[v][2], s; weight = d[v][3])
    end
end

function alg2(G::DictGraph, D::AbstractDict, source::Int, target::Int)

end


text:
""" Undirected wiring diagrams as symmetric monoidal categories.
"""
module MonoidalUndirectedWiringDiagrams
export HypergraphDiagram, HypergraphDiagramOb, HypergraphDiagramHom,
  ObUWD, HomUWD, cospan_action, dom_mask, codom_mask

using AutoHashEquals
using Compat: isnothing

using ...GAT, ...Present, ...CategoricalAlgebra.CSets
using ...Theories: HypergraphCategory
import ...Theories: dom, codom, compose, id, ⋅, ∘, otimes, ⊗, munit, braid, σ,
  mcopy, Δ, mmerge, ∇, delete, ◊, create, □, dunit, dcounit, dagger
using ...CategoricalAlgebra.CSets: disjoint_union
using ..UndirectedWiringDiagrams
using ..UndirectedWiringDiagrams: AbstractUWD, TheoryUWD, TheoryTypedUWD
import ..UndirectedWiringDiagrams: singleton_diagram, junction_diagram

# Data types
############

@present TheoryUntypedHypergraphDiagram <: TheoryUWD begin
  Name::Data
  name::Attr(Box,Name)
end

const UntypedHypergraphDiagram = ACSetType(TheoryUntypedHypergraphDiagram,
  index=[:box, :junction, :outer_junction])

""" Diagram suitable for representing morphisms in hypergraph categories.

An undirected wiring diagram where the ports and junctions are typed and the
boxes have names identifying the morphisms.
"""
@present TheoryHypergraphDiagram <: TheoryTypedUWD begin
  Name::Data
  name::Attr(Box,Name)
end

const HypergraphDiagram = ACSetType(TheoryHypergraphDiagram,
  index=[:box, :junction, :outer_junction])

# Cospan algebra
################

""" Act on undirected wiring diagram with a cospan, as in a cospan algebra.

A *cospan algebra* is a lax monoidal functor from a category of (typed) cospans
(Cospan,⊕) to (Set,×) (Fong & Spivak, 2019, Hypergraph categories, §2.1: Cospans
and cospan-algebras). Undirected wiring diagrams are a cospan algebra in
essentially the same way that any hypergraph category is (Fong & Spivak, 2019,
§4.1: Cospan-algebras and hypergraph categories are equivalent). In fact, we use
this action to implement the hypergraph category structure, particularly
composition, on undirected wiring diagrams.
"""
function cospan_action(f::UWD, args...) where UWD <: AbstractUWD
  ocompose(cospan_diagram(UWD, args...), 1, f)
end

function compose(f::UWD, g::UWD,
                 f_codom::BitVector, g_dom::BitVector) where UWD <: AbstractUWD
  f_dom, g_codom = .!f_codom, .!g_dom
  k, ℓ, ℓ′, m = sum(f_dom), sum(f_codom), sum(g_dom), sum(g_codom)
  @assert ℓ == ℓ′
  f_inner, g_inner = zeros(Int, k+ℓ), zeros(Int, ℓ+m)
  f_inner[f_dom] = 1:k; g_inner[g_codom] = k+ℓ .+ (1:m)
  f_inner[f_codom] = g_inner[g_dom] = k .+ (1:ℓ)
  inner, outer = [f_inner; g_inner], [1:k; k+ℓ .+ (1:m)]
  junction_types = [ port_type(f, f_dom, outer=true);
                     port_type(f, f_codom, outer=true);
                     port_type(g, g_codom, outer=true) ]
  cospan_action(f⊗g, inner, outer, junction_types)
end

otimes(f::UWD, g::UWD) where UWD <: AbstractUWD = disjoint_union(f, g)
⊗(f::UWD, g::UWD) where UWD <: AbstractUWD = otimes(f, g)

# Categorical interface
#######################

""" List of port types representing outer boundary of undirected wiring diagram.
"""
@auto_hash_equals struct ObUWD{UWD <: AbstractUWD, T}
  types::Vector{T}
end
ObUWD{UWD}(types::Vector{T}) where {UWD<:AbstractUWD,T} = ObUWD{UWD,T}(types)

Base.length(A::ObUWD) = length(A.types)
Base.cat(A::ObUWD{UWD}, B::ObUWD{UWD}) where UWD =
  ObUWD{UWD}([A.types; B.types])

""" Undirected wiring diagram with specified domain and codomain.

The outer ports of the undirected wiring diagram are partitioned into domain and
codomain by masks (bit vectors).
"""
@auto_hash_equals struct HomUWD{UWD <: AbstractUWD}
  diagram::UWD
  dom::BitVector

  function HomUWD{UWD}(diagram::UWD;
      dom::Union{BitVector,Nothing}=nothing,
      codom::Union{BitVector,Nothing}=nothing) where UWD <: AbstractUWD
    if !(isnothing(dom) || isnothing(codom))
      @assert codom == .!dom "Domain and codomain masks must form partition"
    elseif !isnothing(codom)
      dom = .!codom
    elseif isnothing(dom)
      error("Must supply at least one of domain and codomain masks")
    end
    length(dom) == length(ports(diagram, outer=true)) ||
      error("Length of (co)domain masks must equal number of outer ports")
    new{UWD}(diagram, dom)
  end
end
HomUWD(diagram::UWD; kw...) where UWD <: AbstractUWD =
  HomUWD{UWD}(diagram; kw...)

dom_mask(f::HomUWD) = f.dom
codom_mask(f::HomUWD) = .!f.dom

const HypergraphDiagramOb{T,Name} = ObUWD{HypergraphDiagram{T,Name},T}
const HypergraphDiagramHom{T,Name} = HomUWD{HypergraphDiagram{T,Name}}

dom_ports(f::HomUWD) = ports(f.diagram, outer=true)[dom_mask(f)]
codom_ports(f::HomUWD) = ports(f.diagram, outer=true)[codom_mask(f)]
dom_port_types(f::HomUWD) = port_type(f.diagram, dom_ports(f), outer=true)
codom_port_types(f::HomUWD) = port_type(f.diagram, codom_ports(f), outer=true)

dom(f::HomUWD{UWD}) where UWD = ObUWD{UWD}(dom_port_types(f))
codom(f::HomUWD{UWD}) where UWD = ObUWD{UWD}(codom_port_types(f))

""" Undirected wiring diagrams as a hypergraph category.

The objects are lists of port types and morphisms are undirected wiring diagrams
whose outer ports are partitioned into domain and codomain.
"""
@instance HypergraphCategory{ObUWD, HomUWD} begin
  @import dom, codom

  function compose(f::HomUWD, g::HomUWD)
    k, ℓ = sum(dom_mask(f)), sum(codom_mask(g))
    HomUWD(compose(f.diagram, g.diagram, codom_mask(f), dom_mask(g)),
           dom=(1:k+ℓ .<= k))
  end

  otimes(A::ObUWD, B::ObUWD) = cat(A, B)
  munit(::Type{ObUWD}) = ObUWD(Symbol[])

  function otimes(f::HomUWD, g::HomUWD)
    HomUWD(otimes(f.diagram, g.diagram), dom=[dom_mask(f); dom_mask(g)])
  end

  id(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, repeat(1:n, 2)), dom=[trues(n); falses(n)])
  end
  braid(A::ObUWD, B::ObUWD) = let m = length(A), n = length(B)
    HomUWD(junction_diagram(A⊗B, [1:m+n; m+1:m+n; 1:m]),
           dom=[trues(m+n); falses(m+n)])
  end
  mcopy(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, repeat(1:n, 3)), dom=[trues(n); falses(2n)])
  end
  mmerge(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, repeat(1:n, 3)), dom=[trues(2n); falses(n)])
  end
  delete(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, 1:n), dom=trues(n))
  end
  create(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, 1:n), dom=falses(n))
  end
  dunit(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, repeat(1:n, 2)), dom=falses(2n))
  end
  dcounit(A::ObUWD) = let n = length(A)
    HomUWD(junction_diagram(A, repeat(1:n, 2)), dom=trues(2n))
  end
  dagger(f::HomUWD) = HomUWD(f.diagram, codom=dom_mask(f))
end

munit(::Type{ObUWD{UWD}}) where UWD = ObUWD{UWD}(Symbol[])
munit(::Type{ObUWD{UWD,T}}) where {UWD,T} = ObUWD{UWD}(T[])

singleton_diagram(A::ObUWD{UWD}; data...) where UWD =
  singleton_diagram(UWD, A.types; data...)
junction_diagram(A::ObUWD{UWD}, outer) where UWD =
  junction_diagram(UWD, outer, A.types)

end


text:
@testset "version" begin

@test isa(CUDAdrv.version(), VersionNumber)

@test isa(CUDAdrv.release(), VersionNumber)
@test CUDAdrv.release().patch == 0

end


text:
using MacroTools

export render, setdisplay, unsetdisplay, getdisplay, current_input, Media, @media, media,
        @render

VERSION < v"0.5-" && const supertype = super

# Some type system utils

distance(S, T) =
  !(S <: T) ? Inf :
  S == T ? 0. :
  1 + distance(supertype(S), T)

nearest(T, U, V) =
  distance(T, U) < distance(T, V) ? U : V

nearest(T, Ts) =
  reduce((U, V) -> nearest(T, U, V), Ts)

function compare2desc(ex)
  (isexpr(ex, :comparison) && ex.args[2] == :(<:)) || return ex
  return Expr(:(<:), ex.args[1], ex.args[3])
end

"""
Similar to `abstract`:

    @media Foo

defines `Foo`, as well as `FooT`, the type representing `Foo`
and its descendants (which is useful for dispatch).

    @media Bar <: Foo
    Bar::FooT
"""
macro media(def)
  T = namify(def)
  def = compare2desc(def)
  quote
    abstract $def
    typealias $(Symbol(string(T, "T"))){T<:$T} Type{T}
    nothing
  end |> esc
end

# The media heirarchy defines an extensible set of possible
# output types. Displayable types are associated with a media
# type as a trait.

@media Graphical
@media Plot <: Graphical
@media Image <: Graphical

@media Textual
@media Numeric <: Textual
@media RichText <: Textual

@media Tabular
@media Matrix <: Tabular
@media List <: Tabular
@media Dataset <: Tabular

"""
`media(T)` gives the media type of the type `T`.
The default is `Textual`.

    media(Gadfly.Plot) == Media.Plot
"""
media(x) = media(typeof(x))

media(T, M) =
  eval(current_module(), :(Media.media{T<:$T}(::Type{T}) = $(Any[M])[1]))

media(Any, Media.Textual)

media(AbstractMatrix, Media.Matrix)
media(AbstractVector, Media.List)

# A "pool" simply associates types with output devices. Obviously
# the idea is to use media types for genericity, but object types
# (e.g. `Float64`, `AbstractMatrix`) can also be used (which will
# override the media trait of the relevant objects).

const _pool = Dict()

defaultpool() = _pool

"""
    setdisplay([input], T, output)

Display `T` objects using `output` when produced by `input`.

`T` is an object type or media type, e.g. `Gadfly.Plot` or `Media.Graphical`.

    display(Editor(), Image, Console())
"""
setdisplay(T, output) =
  defaultpool()[T] = output

unsetdisplay(T) =
  haskey(defaultpool(), T) && delete!(defaultpool(), T)

"""
    getdisplay(T)

Find out what output device `T` will display on.
"""
function getdisplay(T, pool; default = nothing)
  K = nearest(T, [Any, keys(pool)...])
  K == Any && (K = nearest(media(T), keys(pool)))
  K == Any && default ≠ nothing && return default
  return pool[K]
end

# There should be a pool associated with each input device. Normally,
# it should take into account the global pool. The device should
# also override `setdisplay(input, T, output)`

# This design allows e.g. terminal and IJulia displays to be linked
# to their respective inputs, so that both can be used simultaneously
# with the same kernel. At the same time you can link a global display
# to both (e.g. a popup window for plots).

macro defpool(D)
  :(let pool = Dict()
      Media.pool(::$D) = merge(Media.defaultpool(), pool)
      Media.setdisplay(::$D, T, input) = pool[T] = input
      Media.unsetdisplay(::$D, T) = delete!(pool, T)
    end) |> esc
end

# In order to actually display things, we need to know what the current
# input device is. This is stored as a dynamically-bound global variable,
# with the intention that an input device will rebind the current input
# to itself whenever it evaluates code.

# This will also be useful for device-specific functionality like reading
# input and producing warnings.

@dynamic input::Any = nothing

current_input() = input[]

# e.g.

# @dynamic let Media.input = REPL
#   eval(:(render(x)))
# end

# `render` is a stand-in for `display` here.
# Displays should override `render` to display the given object appropriately.

pool() = pool(current_input())
pool(::Void) = defaultpool()

getdisplay(x; default = nothing) =
  getdisplay(x, pool(), default = default)

primarytype(x) = typeof(x).name.primary

render(x) =
  render(getdisplay(primarytype(x)), x)

# Most of the time calls to `render` will defer to rendering some lower-level
# type. `@render` reduces the boilerplate here by automatically calling
# `render` again on the return type.

macro render(T, x, f)
  @capture(T, d_::T′_ | T′_)
  d == nothing && (d = gensym())
  @gensym result
  quote
    function Media.render($d::$T′, $x)
      Media.render($d, $f)
    end
  end |> esc
end


text:
using DiffEqBase
using DiffEqJump
sr = [1.0,2.0,50.]
maj = MassActionJump(sr,[[1 => 1],[1 => 1],[0 => 1]], [[1 => 1], [1 => -1], [1 => 1]])
params = (1.0,2.0,50.)
tspan = (0.,4.)
u0 = [5]
dprob = DiscreteProblem(u0, tspan, params)
jprob = JumpProblem(dprob, Direct(), maj)
solve(EnsembleProblem(jprob), SSAStepper(), EnsembleThreads(); trajectories=10)
solve(EnsembleProblem(jprob;safetycopy=true), SSAStepper(), EnsembleThreads(); trajectories=10)


text:
function import_bibtex(file::AbstractString)
    return BibParser.BibTeX.parse_file(file)[1]
end

function int_to_spaces(n::Int)
    str = ""
    for i in 1:n
        str *= " "
    end
    return str
end

# Dictionnary to handle spaces while exporting BibTeX
const spaces = Dict{AbstractString,AbstractString}(map(
    s -> (string(s) => int_to_spaces(BibInternal.space(s))),
    BibInternal.fields)
)

# Function to write required fields
function field_to_bibtex(
    key::AbstractString,
    value::AbstractString
    )
    space = get(spaces, key, int_to_spaces(BibInternal.space(Symbol(key))))
    return value == "" ? "" : " $key$space = {$value},\n"
end

function export_bibtex(entry::Article)
    str  = "@article{" * entry.id * ",\n"
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("pages", entry.pages)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("volume", entry.volume)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::Book)
    str  = "@book{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("edition", entry.edition)
    str *= field_to_bibtex("editor", entry.editor)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("publisher", entry.publisher)
    str *= field_to_bibtex("series", entry.series)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("volume", entry.volume)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::Booklet)
    str  = "@booklet{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("howpublished", entry.howpublished)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::InBook)
    str  = "@inbook{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("chapter", entry.chapter)
    str *= field_to_bibtex("edition", entry.edition)
    str *= field_to_bibtex("editor", entry.editor)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("pages", entry.pages)
    str *= field_to_bibtex("publisher", entry.publisher)
    str *= field_to_bibtex("series", entry.series)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("type", entry.type)
    str *= field_to_bibtex("volume", entry.volume)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::InCollection)
    str  = "@incollection{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("booktitle", entry.booktitle)
    str *= field_to_bibtex("chapter", entry.chapter)
    str *= field_to_bibtex("edition", entry.edition)
    str *= field_to_bibtex("editor", entry.editor)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("pages", entry.pages)
    str *= field_to_bibtex("publisher", entry.publisher)
    str *= field_to_bibtex("series", entry.series)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("type", entry.type)
    str *= field_to_bibtex("volume", entry.volume)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::InProceedings)
    str  = "@inproceedings{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("booktitle", entry.booktitle)
    str *= field_to_bibtex("editor", entry.editor)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("pages", entry.pages)
    str *= field_to_bibtex("publisher", entry.publisher)
    str *= field_to_bibtex("series", entry.series)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("volume", entry.volume)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::Manual)
    str  = "@manual{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("edition", entry.edition)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("organization", entry.organization)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::MasterThesis)
    str  = "@masterthesis{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("school", entry.school)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("type", entry.type)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::Misc)
    str  = "@misc{" * entry.id * ",\n"
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("howpublished", entry.howpublished)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::PhDThesis)
    str  = "@phdthesis{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("school", entry.school)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("type", entry.type)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::Proceedings)
    str  = "@proceedings{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("editor", entry.editor)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("organization", entry.organization)
    str *= field_to_bibtex("publisher", entry.publisher)
    str *= field_to_bibtex("series", entry.series)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("volume", entry.volume)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::TechReport)
    str  = "@techreport{" * entry.id * ",\n"
    str *= field_to_bibtex("address", entry.address)
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("institution", entry.institution)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("number", entry.number)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("type", entry.type)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(entry::Unpublished)
    str  = "@unpublished{" * entry.id * ",\n"
    str *= field_to_bibtex("author", entry.author)
    str *= field_to_bibtex("key", entry.key)
    str *= field_to_bibtex("month", entry.month)
    str *= field_to_bibtex("note", entry.note)
    str *= field_to_bibtex("title", entry.title)
    str *= field_to_bibtex("year", entry.year)
    for (key, value) in pairs(entry.other)
        str *= field_to_bibtex(key, value)
    end
    return str[1:end - 2] * "\n}"
end

function export_bibtex(bibliography::DataStructures.OrderedSet{AbstractEntry})
    str = ""
    for e in bibliography
        str *= export_bibtex(e) * "\n"
    end
    return str[1:end - 1]
end


text:
getdayfields(days::Int) = (d->d<0 ? Symbol("DM",abs(d)) : Symbol("D",d)).(collect((-1*days):days))

function formeventwindows(event::DataFrame, days::Int;
  refresheventcrsp::Bool = true,
  Fret::Symbol = error("Fret is required"),
  validpermnos::Vector{Int} = unique(event[!, :permno]),
  Fdate::Symbol = error("Fdate is required"),
  retbound::NFloat64 = nothing,
  prefix = error("prefix is required")
  )

  days>0 || error("days must be greater than 0")
  Nevent::Int = size(event,1)
  event.eid = collect(1:Nevent)

  #load the crsp df
  local crsp::DataFrame = prepcrsp(
    crspcolumns = [:date, :permno, :ret, :net, :price, :Lprice],
    crsptype = REPLICATION_TYPE[],
    refreshcrsp=refresheventcrsp,
    validpermnos = validpermnos,
    prefix = "$(prefix)event") do df::DataFrame
    computenetreturns!(df)
    sort!(df,[:permno,:date])

    #lag the price in the simplest fashion
    df.Lprice = Vector{MFloat64}(missing, size(df,1))
    for sdf ∈ groupby(df, :permno)
      Nsdf::Int = size(sdf,1)
      (Nsdf ≤ 1) && continue

      sdf.Lprice[2:end] .= sdf.price[1:(end-1)]
    end
  end

  #make the appropriate sorts
  crspdates::Vector{Date} = sort(unique(crsp.date))

  #winsorize to an fixed vlaue if desired
  deleterows!(crsp, (ismissing).(crsp[!,Fret]))
  if !isnothing(retbound)
    (retbound > 0.0) || error("retbound must be > 0")
    crsp.net .= (f::Float64->(max(min(f,retbound), -retbound))).(crsp.net)
    crsp.ret .= (f::Float64->(max(min(f,retbound), -retbound))).(crsp.ret)
  end

  sort!(event, [:permno, Fdate])

  #assume a fixed set of date for each event
  dateranges::Dict =Dict{Date,NTuple{2, Date}}(
    crspdates[i]=> (crspdates[i-days], crspdates[i+days])
    for i ∈ (1+days):(length(crspdates)-days))

  #get the begin and end date of each range
  filter!(r->haskey(dateranges, r[Fdate]), event)
  event.begindate = (d->dateranges[d][1]).(event[!,Fdate])
  event.enddate = (d->dateranges[d][2]).(event[!,Fdate])
  event.tokeep=trues(size(event,1))

  sevents::GroupedDataFrame = groupby(event, :permno)

  #now check for overlapping dates
  for sevent ∈ groupby(event, :permno)
    Nsevent::Int = size(sevent,1)
    for (i,r) ∈ enumerate(eachrow(sevent))
      t₀::Date = r[Fdate]
      (i > 1) &&  (r[:tokeep] *= t₀>sevent[i-1,:enddate])
      (i < Nsevent) && (r[:tokeep] *= t₀<sevent[i+1,:begindate])
    end
  end
  event = event[event.tokeep, :]

  #now build the panel
  eventXday::Vector{NTuple{2,Int}} = vec(collect(Iterators.product(event.eid, -days:days)))
  newevent::DataFrame = DataFrame(eid= (t->t[1]).(eventXday), day = (t->t[2]).(eventXday))
  event = join(event, newevent, on=:eid)

  dateidx::Dict = Dict(d=>i for (i,d) ∈ enumerate(crspdates))
  event.date = (r::DataFrameRow->crspdates[dateidx[r[Fdate]] + r.day]).(eachrow(event))

  Nevent = size(event,1)
  event = join(event, crsp, on=[:permno, :date], kind=:left)
  @assert Nevent == size(event,1) #make sure we picked the correct keys

  #qc on the returns and net returns
  (s->finiteormissing!(event, s)).([:ret, :net, :price, :Lprice])

  return event
end

function formeventwindowsold(event::DataFrame, days::Int;
  refresheventcrsp::Bool = true,
  Fret::Symbol = error("Fret is required"),
  validpermnos::Vector{Int} = unique(event[!, :permno]),
  Fdate::Symbol = error("Fdate is required"),
  retbound::NFloat64 = nothing,
  prefix = error("prefix is required")
  )

  days>0 || error("days must be greater than 0")
  Nevent::Int = size(event,1)

  #create the fields that will hold the returns
  local dayfields::Vector{Symbol} = getdayfields(days)
  for f::Symbol ∈ dayfields
    event[!, f] = Vector{MFloat64}(missing, Nevent)
  end
  event.Lprice = Vector{MFloat64}(missing, Nevent)
  event.price = Vector{MFloat64}(missing, Nevent)

  #load the crsp df
  local crsp::DataFrame = prepcrsp(
    crspcolumns = [:date, :permno, :ret, :net, :price],
    crsptype = REPLICATION_TYPE[],
    refreshcrsp=refresheventcrsp,
    validpermnos = validpermnos,
    prefix = "$(prefix)event") do df::DataFrame
    computenetreturns!(df)
    sort!(df,[:permno,:date])
  end

  #make the appropriate sorts
  crspdates::Vector{Date} = sort(unique(crsp.date))
  sort!(event, [:permno, Fdate])

  #winsorize to an fixed vlaue if desired
  deleterows!(crsp, (ismissing).(crsp[!,Fret]))
  if !isnothing(retbound)
    (retbound > 0.0) || error("retbound must be > 0")
    crsp[!,Fret] .= (f::Float64->(max(min(f,retbound), -retbound))).(crsp[!,Fret])
  end

  validdateranges::Dict =Dict{Date,NTuple{2, Date}}(
    crspdates[i]=> (crspdates[i-days], crspdates[i+days])
    for i ∈ (1+days):(length(crspdates)-days))

  #allocate space for the information we will get for the study
  event.begindate = Vector{Union{Missing,Date}}(missing, Nevent)
  event.enddate = Vector{Union{Missing,Date}}(missing, Nevent)
  event.tokeep=trues(Nevent)

  #index the dataframe rows
  eventidx::Dict{Int, SubDataFrame} = Dict{Int, SubDataFrame}(
    sevent.permno[1] => sevent for sevent ∈ groupby(event, :permno))

  scrsps::GroupedDataFrame = groupby(crsp, :permno)
  dateidx::Dict = Dict(d=>i for (i,d) ∈ enumerate(crspdates))

  #for each valid permno we have and each event for those permnos
  @mpar for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]
    (!haskey(eventidx, scrsp.permno[1])) && continue
    for r::DataFrameRow ∈ eachrow(eventidx[scrsp.permno[1]])

      t₀::Date = r[Fdate]
      Nscrsp::Int = size(scrsp,1)

      #locate the appropriate entry, and beform first set of validity checks
      i₀::Int = searchsortedfirst(scrsp.date, t₀)
      #println("i₀: $i₀ Nscrsp: $Nscrsp")
      #println("i₀: $i₀\ndays: $days\nlength(scrsp.date): $(length(scrsp.date))
      #  scrsp.date: $(scrsp.date[i₀])\nt₀: $t₀")
      if ((i₀ - days < 1) || #also takes care of certain "not found" cases
        (i₀ + days > Nscrsp) ||
        (scrsp.date[i₀] ≠ t₀) ||
        (!haskey(validdateranges, t₀)))

        r.tokeep = false
        #error("break")
        continue
      end

      #check that the range is complete
      (tₗ::Date,tₕ::Date) = validdateranges[t₀]
      if (tₗ ≠ scrsp.date[i₀-days]) || (tₕ ≠ scrsp.date[i₀+days])
        r.tokeep = false
        continue
      end

      #now copy the return info
      r.price = scrsp[i₀, Fret]
      r.Lprice = scrsp[i₀-1, Fret]
      for (i::Int, f::Symbol) ∈ enumerate(dayfields)
        iₜ::Int = i₀ + (i - days - 1)
        ret::MFloat64 = scrsp[iₜ, Fret]
        r[f] = isfinite(ret) ? ret : missing
      end

    end
  end

  #drop the rows that missed the prior check
  event = event[event.tokeep, :]

  #drop row with incomplete returns
  event = event[completecases(event[:, dayfields]), :]

  return event
end

function formevent(;refreshevent::Bool = true,
    eventpath::String = EVENT_PATH,
    refreshseocrsp::Bool=true,
    refreshdistcrsp::Bool=true,
    eventname::String = EVENT_NAME,
    injlsstream::Function = IN_JLS_STREAM,
    outjlsstream::Function = OUT_JLS_STREAM,
    seodays::Int = SEO_DAYS,
    distdays::Int = DIST_DAYS
    )::NTuple{2,DataFrame}

  local bb::DataFrame, dist::DataFrame
  eventname = "$eventname-$(REPLICATION_TYPE[])"

  if   refreshevent
    bb = readbbseo()
    bb |> CSV.write("output\\bbevent-$(REPLICATION_TYPE[]).csv")
    dist = readcrspevent()

    #keep track of permnos so that we only refresheventcrsp once
    #validpermnos::Vector{Int} = unique([bb.permno; dist.permno])
    bb = formeventwindows(bb,seodays,
      refresheventcrsp=refreshseocrsp,
      validpermnos=unique(bb.permno),
      Fret = SEO_RET,
      Fdate=:date,
      prefix="seo"
      )
    dist = formeventwindows(dist,distdays,
      refresheventcrsp=refreshdistcrsp,
      validpermnos=unique(dist.permno),
      Fret = DIST_RET,
      Fdate=:exdate,
      retbound = DIST_RET_BOUND,
      prefix="dist"
      )
    print("process: ")
    @time processcrspevent!(dist)

    outjlsstream("$eventpath\\$eventname.jls.lz4") do s
      serialize(s, (bb,dist))
    end
  else
    injlsstream("$eventpath\\$eventname.jls.lz4") do s
      (bb,dist) = deserialize(s) #NOTE: DEBUG ONLY
    end
  end

  bb |> CSV.write("output\\cleanbb-$(REPLICATION_TYPE[]).csv")
  #dist |> CSV.write("output\\cleandist-$(REPLICATION_TYPE[]).csv")
  return (bb,dist)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


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

mutable struct RetrieveBootDiagnosticsDataResult <: SwaggerModel
    consoleScreenshotBlobUri::Any # spec type: Union{ Nothing, String } # spec name: consoleScreenshotBlobUri
    serialConsoleLogBlobUri::Any # spec type: Union{ Nothing, String } # spec name: serialConsoleLogBlobUri

    function RetrieveBootDiagnosticsDataResult(;consoleScreenshotBlobUri=nothing, serialConsoleLogBlobUri=nothing)
        o = new()
        validate_property(RetrieveBootDiagnosticsDataResult, Symbol("consoleScreenshotBlobUri"), consoleScreenshotBlobUri)
        setfield!(o, Symbol("consoleScreenshotBlobUri"), consoleScreenshotBlobUri)
        validate_property(RetrieveBootDiagnosticsDataResult, Symbol("serialConsoleLogBlobUri"), serialConsoleLogBlobUri)
        setfield!(o, Symbol("serialConsoleLogBlobUri"), serialConsoleLogBlobUri)
        o
    end
end # type RetrieveBootDiagnosticsDataResult

const _property_map_RetrieveBootDiagnosticsDataResult = Dict{Symbol,Symbol}(Symbol("consoleScreenshotBlobUri")=>Symbol("consoleScreenshotBlobUri"), Symbol("serialConsoleLogBlobUri")=>Symbol("serialConsoleLogBlobUri"))
const _property_types_RetrieveBootDiagnosticsDataResult = Dict{Symbol,String}(Symbol("consoleScreenshotBlobUri")=>"String", Symbol("serialConsoleLogBlobUri")=>"String")
Base.propertynames(::Type{ RetrieveBootDiagnosticsDataResult }) = collect(keys(_property_map_RetrieveBootDiagnosticsDataResult))
Swagger.property_type(::Type{ RetrieveBootDiagnosticsDataResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RetrieveBootDiagnosticsDataResult[name]))}
Swagger.field_name(::Type{ RetrieveBootDiagnosticsDataResult }, property_name::Symbol) =  _property_map_RetrieveBootDiagnosticsDataResult[property_name]

function check_required(o::RetrieveBootDiagnosticsDataResult)
    true
end

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


text:
immutable Multinomial <: DiscreteMultivariateDistribution
    n::Int
    prob::Vector{Float64}

    function Multinomial(n::Integer, p::Vector{Float64})
        if n <= 0
            throw(ArgumentError("n must be a positive integer."))
        end
        if !isprobvec(p)
            throw(ArgumentError("p = $p is not a probability vector."))
        end
        new(int(n), p)
    end

    Multinomial(n::Integer, p::Vector{Float64}, ::NoArgCheck) = new(int(n), p)
    Multinomial(n::Integer, k::Integer) = new(int(n), fill(1.0 / k, k))
end

# Properties

dim(d::Multinomial) = length(d.prob)

mean(d::Multinomial) = d.n .* d.prob

function var(d::Multinomial)
    p = d.prob
    k = length(p)
    v = Array(Float64, k)
    n = d.n
    for i = 1:k
        @inbounds pi = p[i]
        v[i] = n * pi * (1.0 - pi)
    end
    v
end

function cov(d::Multinomial)
    p = d.prob
    k = length(p)
    C = Array(Float64, k, k)
    n = d.n

    for j = 1:k
        pj = p[j]
        for i = 1:j-1
            @inbounds C[i,j] = - n * p[i] * pj
        end

        @inbounds C[j,j] = n * pj * (1.0-pj)
    end

    for j = 1:k-1
        for i = j+1:k
            @inbounds C[i,j] = C[j,i]
        end
    end
    C
end

function mgf(d::Multinomial, t::AbstractVector)
    p = d.prob
    n = d.n
    k = length(p)
    s = 0.0
    for i in 1:k
        s += p[i] * exp(t[i])
    end
    return s^n
end

function cf(d::Multinomial, t::AbstractVector)
    p = d.prob
    n = d.n
    k = length(p)
    s = 0.0 + 0.0im
    for i in 1:k
        s += p[i] * exp(im * t[i])
    end
    return s^n
end

# Evaluation

function insupport{T<:Real}(d::Multinomial, x::Vector{T})
    n = length(x)
    if length(d.prob) != n
        return false
    end
    s = 0.0
    for i in 1:n
        @inbounds xi = x[i]
        if !(isinteger(xi) && xi >= 0)
            return false
        end
        s += xi
    end
    return s == d.n  # integer computation would not yield truncation errors
end

function logpdf{T<:Real}(d::Multinomial, x::Vector{T})
    n = d.n
    p = d.prob
    k = length(p)
    length(x) == k || throw(ArgumentError("Invalid dimension of x."))

    s = lgamma(n + 1.0)
    t = zero(T)
    for i = 1 : k
        @inbounds xi = x[i]
        @inbounds pi = p[i]
        t += xi
        s -= lgamma(xi + 1.0)
        @inbounds s += xi * log(pi)
    end
    return (t == n ? s : -Inf)::Float64
end

pdf{T <: Real}(d::Multinomial, x::Vector{T}) = exp(logpdf(d, x))

# Sampling

function multinom_rand!{T<:Real}(n::Int, p::Vector{Float64}, x::AbstractVector{T})
    k = length(p)
    length(x) == k || throw(ArgumentError("Invalid argument dimension."))

    rp = 1.0  # remaining total probability
    i = 0
    km1 = k - 1

    while i < km1 && n > 0
        i += 1
        @inbounds pi = p[i]
        if pi < rp
            xi = rand(Binomial(n, pi / rp))
            @inbounds x[i] = xi
            n -= xi
            rp -= pi
        else
            # In this case, we don't even have to sample
            # from Binomial. Just assign remaining counts
            # to xi.

            @inbounds x[i] = n
            n = 0
            # rp = 0.0 (no need for this, as rp is no longer needed)
        end
    end

    if i == km1
        @inbounds x[k] = n
    else  # n must have been zero
        z = zero(T)
        for j = i+1 : k
            @inbounds x[j] = z
        end
    end

    return x
end

rand!{T<:Real}(d::Multinomial, x::Vector{T}) = multinom_rand!(d.n, d.prob, x)

function rand!{T<:Real}(d::Multinomial, x::Matrix{T})
    k = dim(d)
    size(x,1) == k || throw(ArgumentError("Invalid argument dimension."))

    for i = 1 : size(x, 2)
        multinom_rand!(d.n, d.prob, view(x, :, i))
    end
    return x
end

immutable MultinomialSampler <: DiscreteMultivariateDistribution
    d::Multinomial
    alias::AliasTable
    function MultinomialSampler(d::Multinomial)
        new(d, AliasTable(d.prob))
    end
end

function rand!{T <: Real}(s::MultinomialSampler, x::Vector{T})
    d::Multinomial = s.d
    n::Int = d.n
    k = dim(s)

    if n^2 > k
        d = s.d
        multinom_rand!(n, d.prob, x)
    else
        # Use an alias table
        fill!(x, convert(T, 0))
        a = s.alias
        for i = 1:n
            x[rand(a)] += 1
        end
    end
    return x
end

dim(s::MultinomialSampler) = length(s.d.prob)
sampler(d::Multinomial) = MultinomialSampler(d)

rand(s::MultinomialSampler) = rand!(s, zeros(Int, dim(s)))

## Fit model

immutable MultinomialStats <: SufficientStats
    n::Int  # number of trials in each experiment
    scnts::Vector{Float64}  # sum of counts
    tw::Float64  # total sample weight

    MultinomialStats(n::Int, scnts::Vector{Float64}, tw::Real) = new(n, scnts, float64(tw))
end

function suffstats{T<:Real}(::Type{Multinomial}, x::Matrix{T})
    K = size(x, 1)
    n::T = zero(T)
    scnts = zeros(K)

    for j = 1:size(x,2)
        nj = zero(T)
        for i = 1:K
            @inbounds xi = x[i,j]
            @inbounds scnts[i] += xi
            nj += xi
        end

        if j == 1
            n = nj
        elseif nj != n
            error("Each sample in X should sum to the same value.")
        end
    end
    MultinomialStats(n, scnts, size(x,2))
end

function suffstats{T<:Real}(::Type{Multinomial}, x::Matrix{T}, w::Array{Float64})
    length(w) == size(x, 2) || throw(ArgumentError("Inconsistent argument dimensions."))

    K = size(x, 1)
    n::T = zero(T)
    scnts = zeros(K)
    tw = 0.

    for j = 1:size(x,2)
        nj = zero(T)
        @inbounds wj = w[j]
        tw += wj
        for i = 1:K
            @inbounds xi = x[i,j]
            @inbounds scnts[i] += xi * wj
            nj += xi
        end

        if j == 1
            n = nj
        elseif nj != n
            error("Each sample in X should sum to the same value.")
        end
    end
    MultinomialStats(n, scnts, tw)
end

fit_mle(::Type{Multinomial}, ss::MultinomialStats) = Multinomial(ss.n, ss.scnts * inv(ss.tw * ss.n))

function fit_mle{T<:Real}(::Type{Multinomial}, x::Matrix{T})
    ss = suffstats(Multinomial, x)
    Multinomial(ss.n, multiply!(ss.scnts, inv(ss.tw * ss.n)))
end

function fit_mle{T<:Real}(::Type{Multinomial}, x::Matrix{T}, w::Array{Float64})
    ss = suffstats(Multinomial, x, w)
    Multinomial(ss.n, multiply!(ss.scnts, inv(ss.tw * ss.n)))
end


text:
function optimize(je_mmle_model::JointEstimationMMLEModel)
    local parameters = je_mmle_model.parameters
    local latents = je_mmle_model.latents
    local dist = je_mmle_model.dist
    local n_index = je_mmle_model.n_index
    local i_index = je_mmle_model.i_index
    local responses_per_item = je_mmle_model.responses_per_item
    local responses_per_examinee = je_mmle_model.responses_per_examinee

    I = size(parameters, 1)
    N = size(latents, 1)

    batch_den = max(1,(nprocs()-1))
    batch_size_N = Int(ceil(N/batch_den))
    batch_size_I = Int(ceil(I/batch_den))

    iter = 1

    # starting values
    before_time = time()
    stop = false
    old_likelihood = 0.0
    old_pars = hcat(_get_parameters_vals.(parameters)...)
    start_time = time()

    latents = pmap(
        (l, idx, r) ->
            l.assessed || _update_posterior(l, parameters[idx], r),
        latents,
        i_index,
        responses_per_examinee,
        batch_size = batch_size_N,
        distributed = true
    )
    # Threads.@threads for n in 1:N
    #     latents[n] = latents[n].assessed || _update_posterior(latents[n], parameters[i_index[n]], responses_per_examinee[n])
    # end
    # latents = @sync @distributed (vcat) for n in 1:N
    #     if latents[n].assessed
    #         latents[n]
    #     else
    #         _update_posterior(latents[n], parameters[i_index[n]], responses_per_examinee[n])
    #     end
    # end
    while !stop
        # calibrate items

        # parameters = @sync @distributed (vcat) for i in 1:I
        #     if parameters[i].calibrated
        #         parameters[i]
        #     else
        #         m_step(parameters[i], latents[n_index[i]], responses_per_item[i], je_mmle_model.int_opt_settings)
        #     end
        # end
        # Threads.@threads for i in 1:I
        #     parameters[i] = parameters[i].calibrated || _calibrate_item_mmle(parameters[i], latents[n_index[i]], responses_per_item[i], je_mmle_model.int_opt_settings)
        # end
        parameters = pmap(
             (p, idx, r) ->
                 p.calibrated || m_step(p, latents[idx], r, je_mmle_model.int_opt_settings),
             parameters,
             n_index,
             responses_per_item,
             batch_size = batch_size_I,
             distributed = true
        )

        # #update posteriors

        # latents = @sync @distributed (vcat) for n in 1:N
        #     if latents[n].assessed
        #         latents[n]
        #     else
        #         _update_posterior(latents[n], parameters[i_index[n]], responses_per_examinee[n])
        #     end
        # end
        # Threads.@threads for n in 1:N
        #     latents[n] = latents[n].assessed || _update_posterior(latents[n], parameters[i_index[n]], responses_per_examinee[n])
        # end

        latents = pmap(
            (l, idx, r) ->
                l.assessed || _update_posterior(l, parameters[idx], r),
            latents,
            i_index,
            responses_per_examinee,
            batch_size = batch_size_N,
            distributed = true
        )
        #rescale dist
        if je_mmle_model.rescale_latent
            dist = _rescale(
                dist,
                latents;
                metric = je_mmle_model.metric
            )
            for n = 1:N
                latents[n].prior = dist
            end
        end

        if any([
            check_iter(iter; max_iter = Int64(je_mmle_model.ext_opt_settings[1]), verbosity = Int(je_mmle_model.ext_opt_settings[5])),
            check_time(start_time; max_time = Int64(je_mmle_model.ext_opt_settings[2]), verbosity = Int(je_mmle_model.ext_opt_settings[5])),
            check_f_tol_rel!(
                latents,
                old_likelihood;
                f_tol_rel = je_mmle_model.ext_opt_settings[3],
                verbosity = Int(je_mmle_model.ext_opt_settings[5])
                ),
            check_x_tol_rel!(
                parameters,
                old_pars;
                x_tol_rel = je_mmle_model.ext_opt_settings[4],
                verbosity = Int(je_mmle_model.ext_opt_settings[5])
            )]
            )
            stop = true
        end
        iter += 1
    end

    return parameters::Vector{<:AbstractParameters}, latents::Vector{<:AbstractLatent}, dist::Distributions.DiscreteUnivariateDistribution
end


text:
UI = UncertainIndexDataset(example_uidxs)

n = 3
@test resample(UI) isa Vector{<:Real}
@test resample(UI, n) isa Vector{Vector{T}} where T <: Real


text:
#!/usr/bin/env julia
#=
Using 17k map and the reduced set of loci, randomly select 1000 ID as reference
and randomly select 50 ID to mask. Them impute them back and check error rates.

This used 4,796 17k data.
=#

using DelimitedFiles, Plots, Statistics
dat = readdlm("num.txt")
dat = reshape(dat, 2, :)'
n = length(dat[:, 1])
ave = mean(dat[:, 1])

plot(dat[:, 1], label="Error rates", dpi=300)
plot!([1, n], [ave, ave], label="Mean error rates", ylabel="Error rate", dpi=300)

savefig("reduced-l2m-v3.png")

max_stars_repo_name: xijiang/Refactered-NSG

max_issues_repo_name: xijiang/Refactered-NSG

max_forks_repo_name: xijiang/Refactered-NSG


text:
# Defines different storages for zarr arrays. Currently only regular files (DirectoryStore)
# and Dictionaries are supported
abstract type AbstractStore end

#Define the interface
"""
  storagesize(d::AbstractStore)

This function shall return the size of all data files in a store.
"""
function storagesize end

"""
    Base.getindex(d::AbstractStore,i::String)

Returns the data stored in the given key as a Vector{UInt8}
"""
Base.getindex(d::AbstractStore,i::AbstractString) = error("getindex not implemented for store $(typeof(d))")

"""
    Base.setindex!(d::AbstractStore,v,i::String)

Writes the values in v to the given store and key.
"""
Base.setindex!(d::AbstractStore,v,i::AbstractString) = error("setindex not implemented for store $(typeof(d))")

"""
    subdirs(d::AbstractStore, p)

Returns a list of keys for children stores in the given store at path p.
"""
function subdirs end

"""
    subkeys(d::AbstractStore, p)

Returns the keys of files in the given store.
"""
function subkeys end

"""
    Base.delete!(d::AbstractStore, k::String)

Deletes the given key from the store.
"""

citostring(i::CartesianIndex) = join(reverse((i - oneunit(i)).I), '.')
_concatpath(p,s) = isempty(p) ? s : rstrip(p,'/') * '/' * s

Base.getindex(s::AbstractStore, p, i::CartesianIndex) = s[p, citostring(i)]
Base.getindex(s::AbstractStore, p, i) = s[_concatpath(p,i)]
Base.delete!(s::AbstractStore, p, i::CartesianIndex) = delete!(s, p, citostring(i))
Base.delete!(s::AbstractStore, p, i) = delete!(s, _concatpath(p,i))
Base.haskey(s::AbstractStore, k) = isinitialized(s,k)
Base.setindex!(s::AbstractStore,v,p,i) = setindex!(s,v,_concatpath(p,i))
Base.setindex!(s::AbstractStore,v,p,i::CartesianIndex) = s[p, citostring(i)]=v

maybecopy(x) = copy(x)
maybecopy(x::String) = x

function getattrs(s::AbstractStore, p)
  atts = s[p,".zattrs"]
  if atts === nothing
    Dict()
  else
    JSON.parse(replace(String(maybecopy(atts)),": NaN,"=>": \"NaN\","))
  end
end
function writeattrs(s::AbstractStore, p, att::Dict)
  b = IOBuffer()
  JSON.print(b,att)
  s[p,".zattrs"] = take!(b)
  att
end

is_zgroup(s::AbstractStore, p) = isinitialized(s,_concatpath(p,".zgroup"))
is_zarray(s::AbstractStore, p) = isinitialized(s,_concatpath(p,".zarray"))

isinitialized(s::AbstractStore, p, i::CartesianIndex)=isinitialized(s,p,citostring(i))
isinitialized(s::AbstractStore, p, i) = isinitialized(s,_concatpath(p,i))
isinitialized(s::AbstractStore, i) = s[i] !== nothing

getmetadata(s::AbstractStore, p,fill_as_missing) = Metadata(String(maybecopy(s[p,".zarray"])),fill_as_missing)
function writemetadata(s::AbstractStore, p, m::Metadata)
  met = IOBuffer()
  JSON.print(met,m)
  s[p,".zarray"] = take!(met)
  m
end

isemptysub(s::AbstractStore, p) = isempty(subkeys(s,p)) && isempty(subdirs(s,p))

#Here different storage backends can register regexes that are checked against
#during auto-check of storage format when doing zopen
storageregexlist = Pair[]

include("directorystore.jl")
include("dictstore.jl")
include("s3store.jl")
include("gcstore.jl")
include("consolidated.jl")
include("http.jl")


text:
Combinatorics.combinations(a, c, k::Integer) = MultisetCombinations(a, c, k)

struct MultisetCombinations{T,U,V}
    a::T
    c::U
    k::V
end

Base.start(c::MultisetCombinations) = begin
    s = [fill(i, c.c[i]) for i in eachindex(c.a)] |> flatten
    c.k <= length(s) ? s[1:c.k] : [length(c.a) + 1]
end

Base.next(c::MultisetCombinations, s) = begin
    combination = [c.a[si] for si in s]
    c.k > 0 || return combination, [length(c.a) + 1]

    s = copy(s)
    for i in length(s):-1:1
        s[i] += 1
        while s[i] <= length(c.a) && countnz(s .== s[i]) > c.c[s[i]]
            s[i] += 1
        end

        if s[i] <= length(c.a)
            t = [fill(i, c.c[i] - countnz(s .== i)) for i in s[i]:length(c.a)] |> flatten
            if length(t) >= length(s)-i
                s[i+1:end] = t[1:(length(s)-i)]
            else
                s = [length(c.a) + 1]
            end
            break
        end
    end

    combination, s
end

Base.done(c::MultisetCombinations, s) = !isempty(s) && s[1] > length(c.a)

Base.eltype(c::MultisetCombinations) = eltype(typeof(c))
Base.eltype(::Type{MultisetCombinations{T,U,V}}) where {T,U,V} = Array{eltype(T),1}

Base.length(c::MultisetCombinations) = begin
    c.k == 0 && return big(1)

    csum = sum(c.c)
    c.k > csum && return big(0)
    c.k == csum && return big(1)

    # http://www.m-hikari.com/ams/ams-2011/ams-17-20-2011/siljakAMS17-20-2011.pdf
    R, z = PowerSeriesRing(ZZ, c.k + 1, "z")
    gf = prod(sum(z^j for j in 0:m) for m in c.c)
    l = coeff(gf, c.k)

    convert(BigInt, l)
end


text:
module TestModelStructure_VariableTimestep

#tests the framework of components and connections

using Test
using Mimi

import Mimi:
    connect_param!, unconnected_params, set_dimension!,
    reset_compdefs, numcomponents, get_connections, internal_param_conns, dim_count,
    compdef, getproperty, setproperty!, dimension, dimensions, compdefs

reset_compdefs()

@defcomp A begin
    varA::Int = Variable(index=[time])
    parA::Int = Parameter()

    function run_timestep(p, v, d, t)
        v.varA[t] = p.parA
    end
end

@defcomp B begin
    varB::Int = Variable()

    function run_timestep(p, v, d, t)
        if t.t < 10
            v.varB = 1
        else
            v.varB = 10
        end
    end
end

@defcomp C begin
    varC::Int = Variable()
    parC::Int = Parameter()

    function run_timestep(p, v, d, t)
        v.varC = p.parC
    end
end

years = [2015:5:2100; 2110:10:2200]
first_A = 2050
last_A = 2150

m = Model()
set_dimension!(m, :time, years)

@test_throws ErrorException add_comp!(m, A, last = 2210)
@test_throws ErrorException add_comp!(m, A, first = 2010)
@test_throws ArgumentError add_comp!(m, A, after=:B)
# @test_throws ErrorException add_comp!(m, A, after=:B)
add_comp!(m, A, first = first_A, last = last_A) #test specific last and first

add_comp!(m, B, before=:A)

add_comp!(m, C, after=:B) # test a later first than model
# Component order is B -> C -> A.

connect_param!(m, :A, :parA, :C, :varC)

unconn = unconnected_params(m)
@test length(unconn) == 1
@test unconn[1] == (:C, :parC)

connect_param!(m, :C => :parC, :B => :varB)

@test_throws ErrorException add_comp!(m, C, after=:A, before=:B)

@test numcomponents(m.md) == 3

@test length(internal_param_conns(m)) == 2
@test get_connections(m, :A, :incoming)[1].src_comp_name == :C

@test length(get_connections(m, :B, :incoming)) == 0
@test get_connections(m, :B, :outgoing)[1].dst_comp_name == :C

@test length(get_connections(m, :A, :all)) == 1

@test length(unconnected_params(m)) == 0

#############################################
#  Tests for connecting scalar parameters   #
#############################################

run(m)

# @test all([m[:A, :varA][t] == 1 for t in 1:2])
# @test all([m[:A, :varA][t] == 10 for t in 3:16])

dim = Mimi.Dimension(years)
varA = m[:A, :varA][dim[first_A]:dim[last_A]]
@test all([varA[i] == 1 for i in 1:2])
@test all([varA[i] == 10 for i in 3:16])

##########################
#   tests for indexing   #
##########################

@test dim_count(m.md, :time) == 28

@test m[:A, :parA] == 10
@test_throws ErrorException m[:A, :xx]

time = dimension(m, :time)
a = collect(keys(time))
@test all([a[i] == years[i] for i in 1:28])

@test dimensions(m, :A, :varA)[1] == :time
@test length(dimensions(m, :A, :parA)) == 0

################################
#  tests for delete! function  #
################################

@test_throws ErrorException delete!(m, :D)
@test length(m.md.internal_param_conns) == 2
delete!(m, :A)
@test length(m.md.internal_param_conns) == 1
@test !(:A in compdefs(m))
@test length(compdefs(m)) == 2

#######################################
#   Test check for unset parameters   #
#######################################

@defcomp D begin
  varD = Variable(index=[time])
  parD = Parameter()
end

add_comp!(m, D)
@test_throws ErrorException Mimi.build(m)

end # module


text:
using Distributed
addprocs(48; exeflags=["--project"])
@everywhere begin
    using Serialization
    using RenyiEntropyOfMagic
end

N = parse(Int, ARGS[1])
results = scan_tfim_renyi_entropy(N; nsamples=10^6)
serialize(pkgdir(RenyiEntropyOfMagic, "data", "$N-sites-result.dat"), results)