text:
function segway_rk4(xt::Vector{Float64}, ut::AbstractVector, dt::Number, params::Vector)
    xt1 = copy(xt)
    segway_rk4!(xt1, ut, dt, params)
    return xt1
end

function segway_rk4!(xt::AbstractVector, ut::AbstractVector, dt::Number, params::Vector)
    temp = Array{Float64}(undef, 7)

    #get partial results dy1-dy4
    dy1 = dxdt_segway(xt, ut, params);
    temp .= dy1.*(dt/2.0) .+ xt;

    dy2 = dxdt_segway(temp, ut, params);
    temp .= dy2.*(dt/2.0) .+ xt;

    dy3 = dxdt_segway(temp, ut, params);
    temp .= dy3.*dt .+ xt;

    dy4 = dxdt_segway(temp, ut, params);

    # update xt
    @. temp = dt*(dy1+2.0*(dy2+dy3)+dy4)/6.0

    #Check if segway has fallen #TODO limits anpassen
    if -pi/2>=xt[4]+temp[4]
        xt[4] = -pi/2
		xt[5:7] = 0.0
    elseif pi/2<=xt[4]+temp[4]
        xt[4] = pi/2
		xt[5:7] = 0.0
    else
        # Limit angle to lie between 0 and 2*pi
        @. xt += temp
        #xt[3] = mod(xt[3], 2*pi)
    end
    return nothing
end

"""
Segway dynamics, taken from Pathak2005.
State vector is [x0, y0, phi, alpha, dalpha, v, dphi]
"""
function dxdt_segway(x::Vector{Float64}, u::AbstractVector, params::Vector)
    g = 9.81 #N/kg

    # Rename variables
    Mb = params[1]
    Mw = params[2]
    R = params[3]
    cz = params[4]
    b = params[5]
    Ixx = params[6]
    Iyy = params[7]
    Izz = params[8]
    Iwa = params[9]
    Iwd = params[10]

    phi = x[3]
    alpha = x[4]
    dalpha = x[5]
    v = x[6]
    dphi = x[7]

    # Rigid Body Dynamics
    Dalpha = Mb^2*cos(alpha)^2*cz^2*R^2+((-Mb^2-2*Mw*Mb)*cz^2-2*Iyy*Mw-Iyy*Mb)*R^2-2*Mb*cz^2*Iwa-2*Iyy*Iwa
    Galpha = (-Mb*cz^2+Izz-Ixx)*R^2*cos(alpha)^2 + (Mb*cz^2+Ixx+2*Iwd+2*b^2*Mw)*R^2+2*b^2*Iwa
    H = 1/2*Mb*R^2*Izz+Iwa*Izz-Mw*R^2*Ixx-Iwa*Ixx-Mb*cz^2*Mw*R^2 - Mb*cz^2*Iwa -1/2*Mb*R^2*Ixx+Mw*R^2*Izz
    Kalpha = (-4*Iyy*Mb*R^2*cz-3*R^2*Mb^2*cz^3+Mb*R^2*cz*(Ixx-Izz))*sin(alpha) + (Mb*R^2*cz*(Ixx-Izz)+R^2*Mb^2*cz^3)*sin(3*alpha)
    f21 = sin(2*alpha)*dphi^2*H/Dalpha + Mb^2*cz^2*R^2*sin(2*alpha)*dalpha^2/(2*Dalpha) + (-2*Mb^2*R^2*cz-4*Iwa*Mb*cz-4*Mw*R^2*Mb*cz)*g*sin(alpha)/(2*Dalpha)
    f22 = Kalpha*dphi^2+(Mb^2*cz^2*R^2*g*sin(2*alpha))/(2*Dalpha)+(-4*Iyy*Mb*R^2*cz-4*R^2*Mb^2*cz^3)*sin(alpha)*dalpha^2/(4*Dalpha)
    f23 = (-(Ixx-Izz)*R^2-Mb*cz^2*R^2)*sin(2*alpha)*dalpha*dphi/Galpha-sin(alpha)*R^2*Mb*cz*v*dphi/Galpha
    g21 = (u[1]+u[2])*(Mb*R^2+2*Mw*R^2+2*Iwa+Mb*cos(alpha)*cz*R)/Dalpha
    g22 = -(u[1]+u[2])*R*(Mb*cos(alpha)*cz*R+Iyy+Mb*cz^2)/Dalpha
    g23 = (u[1]-u[2])*R*b/Galpha

    xdot = Array{Float64}(undef, 7)
    xdot[1] = cos(phi)*v
    xdot[2] = sin(phi)*v
    xdot[3] = x[7]
    xdot[4] = x[5]
    xdot[5] = f21 + g21
    xdot[6] = f22 + g22
    xdot[7] = f23 + g23

    return xdot
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
function boundaries_to_sparse(boundary)
    n = length(boundary)
    I = Vector{Int}(undef, n)
    J = Vector{Int}(undef, n)
    V = Vector{Bool}(undef, n)
    for (idx, el) in enumerate(boundary)
        cell, face = el
        I[idx] = face
        J[idx] = cell
        V[idx] = true
    end
    return sparse(I, J, V)
end

"""
    generate_grid(celltype::Cell, nel::NTuple, [left::Vec, right::Vec)

Return a `Grid` for a rectangle in 1, 2 or 3 dimensions. `celltype` defined the type of cells,
e.g. `Triangle` or `Hexahedron`. `nel` is a tuple of the number of elements in each direction.
`left` and `right` are optional endpoints of the domain. Defaults to -1 and 1 in all directions.
"""
generate_grid

# Line
function generate_grid(::Type{Line}, nel::NTuple{1,Int}, left::Vec{1,T}=Vec{1}((-1.0,)), right::Vec{1,T}=Vec{1}((1.0,))) where {T}
    nel_x = nel[1]
    n_nodes = nel_x + 1

    # Generate nodes
    coords_x = collect(range(left[1], stop=right[1], length=n_nodes))
    nodes = Node{1,T}[]
    for i in 1:n_nodes
        push!(nodes, Node((coords_x[i],)))
    end

    # Generate cells

    cells = Line[]
    for i in 1:nel_x
        push!(cells, Line((i, i+1)))
    end

    # Cell faces
    boundary = Vector([(1, 1),
                       (nel_x, 2)])

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    facesets = Dict("left"  => Set{Tuple{Int,Int}}([boundary[1]]),
                    "right" => Set{Tuple{Int,Int}}([boundary[2]]))
    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# QuadraticLine
function generate_grid(::Type{QuadraticLine}, nel::NTuple{1,Int}, left::Vec{1,T}=Vec{1}((-1.0,)), right::Vec{1,T}=Vec{1}((1.0,))) where {T}
    nel_x = nel[1]
    n_nodes = 2*nel_x + 1

    # Generate nodes
    coords_x = collect(range(left[1], stop=right[1], length=n_nodes))
    nodes = Node{1,T}[]
    for i in 1:n_nodes
        push!(nodes, Node((coords_x[i],)))
    end

    # Generate cells
    cells = QuadraticLine[]
    for i in 1:nel_x
        push!(cells, QuadraticLine((2*i-1, 2*i+1, 2*i)))
    end

    # Cell faces
    boundary = Tuple{Int,Int}[(1, 1),
                         (nel_x, 2)]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    facesets = Dict("left"  => Set{Tuple{Int,Int}}([boundary[1]]),
                    "right" => Set{Tuple{Int,Int}}([boundary[2]]))
    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

function _generate_2d_nodes!(nodes, nx, ny, LL, LR, UR, UL)
      for i in 0:ny-1
        ratio_bounds = i / (ny-1)

        x0 = LL[1] * (1 - ratio_bounds) + ratio_bounds * UL[1]
        x1 = LR[1] * (1 - ratio_bounds) + ratio_bounds * UR[1]

        y0 = LL[2] * (1 - ratio_bounds) + ratio_bounds * UL[2]
        y1 = LR[2] * (1 - ratio_bounds) + ratio_bounds * UR[2]

        for j in 0:nx-1
            ratio = j / (nx-1)
            x = x0 * (1 - ratio) + ratio * x1
            y = y0 * (1 - ratio) + ratio * y1
            push!(nodes, Node((x, y)))
        end
    end
end

function generate_grid(C::Type{Cell{2,M,N}}, nel::NTuple{2,Int}, X::Vector{Vec{2,T}}) where {M,N,T}
    @assert length(X) == 4
    generate_grid(C, nel, X[1], X[2], X[3], X[4])
end

function generate_grid(C::Type{Cell{2,M,N}}, nel::NTuple{2,Int}, left::Vec{2,T}=Vec{2}((-1.0,-1.0)), right::Vec{2,T}=Vec{2}((1.0,1.0))) where {M,N,T}
    LL = left
    UR = right
    LR = Vec{2}((UR[1], LL[2]))
    UL = Vec{2}((LL[1], UR[2]))
    generate_grid(C, nel, LL, LR, UR, UL)
end

# Quadrilateral
function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = nel_x*nel_y
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = Quadrilateral[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, Quadrilateral((node_array[i,j], node_array[i+1,j], node_array[i+1,j+1], node_array[i,j+1])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
                              [(cl, 2) for cl in cell_array[end,:]];
                              [(cl, 3) for cl in cell_array[:,end]];
                              [(cl, 4) for cl in cell_array[1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# QuadraticQuadrilateral
function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = nel_x*nel_y
    n_nodes_x = 2*nel_x + 1; n_nodes_y = 2*nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = QuadraticQuadrilateral[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, QuadraticQuadrilateral((node_array[2*i-1,2*j-1],node_array[2*i+1,2*j-1],node_array[2*i+1,2*j+1],node_array[2*i-1,2*j+1],
                                             node_array[2*i,2*j-1],node_array[2*i+1,2*j],node_array[2*i,2*j+1],node_array[2*i-1,2*j],
                                             node_array[2*i,2*j])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
                              [(cl, 2) for cl in cell_array[end,:]];
                              [(cl, 3) for cl in cell_array[:,end]];
                              [(cl, 4) for cl in cell_array[1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# Hexahedron
function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1; n_nodes_z = nel_z + 1
    n_nodes = n_nodes_x * n_nodes_y * n_nodes_z

    # Generate nodes
    coords_x = range(left[1], stop=right[1], length=n_nodes_x)
    coords_y = range(left[2], stop=right[2], length=n_nodes_y)
    coords_z = range(left[3], stop=right[3], length=n_nodes_z)
    nodes = Node{3,T}[]
    for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
        push!(nodes, Node((coords_x[i], coords_y[j], coords_z[k])))
    end

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z))
    cells = Hexahedron[]
    for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x
        push!(cells, Hexahedron((node_array[i,j,k], node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k],
                                 node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y, nel_z))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,:,1][:]];
                              [(cl, 2) for cl in cell_array[:,1,:][:]];
                              [(cl, 3) for cl in cell_array[end,:,:][:]];
                              [(cl, 4) for cl in cell_array[:,end,:][:]];
                              [(cl, 5) for cl in cell_array[1,:,:][:]];
                              [(cl, 6) for cl in cell_array[:,:,end][:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,1][:]))   .+ offset]); offset += length(cell_array[:,:,1][:])
    facesets["front"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1,:][:]))   .+ offset]); offset += length(cell_array[:,1,:][:])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:,:][:])) .+ offset]); offset += length(cell_array[end,:,:][:])
    facesets["back"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end,:][:])) .+ offset]); offset += length(cell_array[:,end,:][:])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# Triangle
function generate_grid(::Type{Triangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = 2*nel_x*nel_y
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = Triangle[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, Triangle((node_array[i,j], node_array[i+1,j], node_array[i,j+1]))) # ◺
        push!(cells, Triangle((node_array[i+1,j], node_array[i+1,j+1], node_array[i,j+1]))) # ◹
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
                               [(cl, 1) for cl in cell_array[2,end,:]];
                               [(cl, 2) for cl in cell_array[2,:,end]];
                               [(cl, 3) for cl in cell_array[1,1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# QuadraticTriangle
function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = 2*nel_x*nel_y
    n_nodes_x = 2*nel_x + 1; n_nodes_y = 2*nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = QuadraticTriangle[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, QuadraticTriangle((node_array[2*i-1,2*j-1], node_array[2*i+1,2*j-1], node_array[2*i-1,2*j+1],
                                        node_array[2*i,2*j-1], node_array[2*i,2*j], node_array[2*i-1,2*j]))) # ◺
        push!(cells, QuadraticTriangle((node_array[2*i+1,2*j-1], node_array[2*i+1,2*j+1], node_array[2*i-1,2*j+1],
                                        node_array[2*i+1,2*j], node_array[2*i,2*j+1], node_array[2*i,2*j]))) # ◹
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
                              [(cl, 1) for cl in cell_array[2,end,:]];
                              [(cl, 2) for cl in cell_array[2,:,end]];
                              [(cl, 3) for cl in cell_array[1,1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# Tetrahedron
function generate_grid(::Type{Tetrahedron}, cells_per_dim::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nodes_per_dim = cells_per_dim .+ 1

    cells_per_cube = 6
    total_nodes = prod(nodes_per_dim)
    total_elements = cells_per_cube * prod(cells_per_dim)

    n_nodes_x, n_nodes_y, n_nodes_z = nodes_per_dim
    n_cells_x, n_cells_y, n_cells_z = cells_per_dim

    # Generate nodes
    coords_x = range(left[1], stop=right[1], length=n_nodes_x)
    coords_y = range(left[2], stop=right[2], length=n_nodes_y)
    coords_z = range(left[3], stop=right[3], length=n_nodes_z)
    numbering = reshape(1:total_nodes, nodes_per_dim)

    # Pre-allocate the nodes & cells
    nodes = Vector{Node{3,T}}(undef, total_nodes)
    cells = Vector{Tetrahedron}(undef, total_elements)

    # Generate nodes
    node_idx = 1
    @inbounds for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
        nodes[node_idx] = Node((coords_x[i], coords_y[j], coords_z[k]))
        node_idx += 1
    end

    # Generate cells, case 1 from: http://www.baumanneduard.ch/Splitting%20a%20cube%20in%20tetrahedras2.htm
    # cube = (1, 2, 3, 4, 5, 6, 7, 8)
    # left = (1, 4, 5, 8), right = (2, 3, 6, 7)
    # front = (1, 2, 5, 6), back = (3, 4, 7, 8)
    # bottom = (1, 2, 3, 4), top = (5, 6, 7, 8)
    cell_idx = 0
    @inbounds for k in 1:n_cells_z, j in 1:n_cells_y, i in 1:n_cells_x
        cell = (
            numbering[i  , j  , k],
            numbering[i+1, j  , k],
            numbering[i+1, j+1, k],
            numbering[i  , j+1, k],
            numbering[i  , j  , k+1],
            numbering[i+1, j  , k+1],
            numbering[i+1, j+1, k+1],
            numbering[i  , j+1, k+1]
        )

        cells[cell_idx + 1] = Tetrahedron((cell[1], cell[2], cell[4], cell[8]))
        cells[cell_idx + 2] = Tetrahedron((cell[1], cell[5], cell[2], cell[8]))
        cells[cell_idx + 3] = Tetrahedron((cell[2], cell[3], cell[4], cell[8]))
        cells[cell_idx + 4] = Tetrahedron((cell[2], cell[7], cell[3], cell[8]))
        cells[cell_idx + 5] = Tetrahedron((cell[2], cell[5], cell[6], cell[8]))
        cells[cell_idx + 6] = Tetrahedron((cell[2], cell[6], cell[7], cell[8]))

        cell_idx += cells_per_cube
    end

    # Order the cells as c_nxyz[n, x, y, z] such that we can look up boundary cells
    c_nxyz = reshape(1:total_elements, (cells_per_cube, cells_per_dim...))

    @views le = [map(x -> (x,4), c_nxyz[1, 1, :, :][:])   ; map(x -> (x,2), c_nxyz[2, 1, :, :][:])]
    @views ri = [map(x -> (x,1), c_nxyz[4, end, :, :][:]) ; map(x -> (x,1), c_nxyz[6, end, :, :][:])]
    @views fr = [map(x -> (x,1), c_nxyz[2, :, 1, :][:])   ; map(x -> (x,1), c_nxyz[5, :, 1, :][:])]
    @views ba = [map(x -> (x,3), c_nxyz[3, :, end, :][:]) ; map(x -> (x,3), c_nxyz[4, :, end, :][:])]
    @views bo = [map(x -> (x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> (x,1), c_nxyz[3, :, :, 1][:])]
    @views to = [map(x -> (x,3), c_nxyz[5, :, :, end][:]) ; map(x -> (x,3), c_nxyz[6, :, :, end][:])]

    boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba])

    facesets = Dict(
        "left" => Set(le),
        "right" => Set(ri),
        "front" => Set(fr),
        "back" => Set(ba),
        "bottom" => Set(bo),
        "top" => Set(to),
    )

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

ext: jl


text:
# ---
# title: 451. Sort Characters By Frequency
# id: problem451
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Heap
# link: <https://leetcode.com/problems/sort-characters-by-frequency/description/>
# hidden: true
# ---
#
# Given a string, sort it in decreasing order based on the frequency of
# characters.
#
# **Example 1:**
#
#
#
#     Input:
#     "tree"
#
#     Output:
#     "eert"
#
#     Explanation:
#     'e' appears twice while 'r' and 't' both appear once.
#     So 'e' must appear before both 'r' and 't'. Therefore "eetr" is also a valid answer.
#
#
# **Example 2:**
#
#
#
#     Input:
#     "cccaaa"
#
#     Output:
#     "cccaaa"
#
#     Explanation:
#     Both 'c' and 'a' appear three times, so "aaaccc" is also a valid answer.
#     Note that "cacaca" is incorrect, as the same characters must be together.
#
#
# **Example 3:**
#
#
#
#     Input:
#     "Aabb"
#
#     Output:
#     "bbAa"
#
#     Explanation:
#     "bbaA" is also a valid answer, but "Aabb" is incorrect.
#     Note that 'A' and 'a' are treated as two different characters.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
__precompile__()

module Includes

export @include_once

mtimes = Dict{String,Float64}()

macro include_once(file::String)
    mt = mtime(file);
    if mt > get!(mtimes, file, 0)
        mtimes[file] = mt
        return :( include($file) )
    elseif mt == 0
        error("File $file does not exist")
    else
        return
    end
end

end

max_stars_repo_name: Ashymad/JuliaModules

max_issues_repo_name: Ashymad/JuliaModules

max_forks_repo_name: Ashymad/JuliaModules


text:
using LibCEED.C, Printf, UnsafeArrays
using LibCEED.C: CeedInt, CeedScalar

# A structure used to pass additional data to f_build_mass
mutable struct BuildContextC
    dim::CeedInt
    space_dim::CeedInt
end

# libCEED Q-function for building quadrature data for a mass operator
function f_build_mass_c(
    ctx_ptr::Ptr{Cvoid},
    Q::CeedInt,
    in_ptr::Ptr{Ptr{CeedScalar}},
    out_ptr::Ptr{Ptr{CeedScalar}},
)
    # in[0] is Jacobians with shape [dim, nc=dim, Q]
    # in[1] is quadrature weights, size (Q)
    ctx = unsafe_load(Ptr{BuildContextC}(ctx_ptr))
    J = UnsafeArray(unsafe_load(in_ptr, 1), (Int(Q), Int(ctx.dim^2)))
    w = UnsafeArray(unsafe_load(in_ptr, 2), (Int(Q),))
    qdata = UnsafeArray(unsafe_load(out_ptr, 1), (Int(Q),))
    if ctx.dim == 1
        @inbounds @simd for i = 1:Q
            qdata[i] = J[i]*w[i]
        end
    elseif ctx.dim == 2
        @inbounds @simd for i = 1:Q
            qdata[i] = (J[i, 1]*J[i, 4] - J[i, 2]*J[i, 3])*w[i]
        end
    elseif ctx.dim == 3
        @inbounds @simd for i = 1:Q
            qdata[i] =
                (
                    J[i, 1]*(J[i, 5]*J[i, 9] - J[i, 6]*J[i, 8]) -
                    J[i, 2]*(J[i, 4]*J[i, 9] - J[i, 6]*J[i, 7]) +
                    J[i, 3]*(J[i, 4]*J[i, 8] - J[i, 5]*J[i, 7])
                )*w[i]
        end
    else
        error("Bad dimension")
    end
    return CeedInt(0)
end

# libCEED Q-function for applying a mass operator
function f_apply_mass_c(
    ctx,
    Q::CeedInt,
    in_ptr::Ptr{Ptr{CeedScalar}},
    out_ptr::Ptr{Ptr{CeedScalar}},
)
    u = UnsafeArray(unsafe_load(in_ptr, 1), (Int(Q),))
    qdata = UnsafeArray(unsafe_load(in_ptr, 2), (Int(Q),))
    v = UnsafeArray(unsafe_load(out_ptr, 1), (Int(Q),))
    @inbounds @simd for i = 1:Q
        v[i] = qdata[i]*u[i]
    end
    return CeedInt(0)
end

function get_cartesian_mesh_size_c(dim, order, prob_size)
    dims = zeros(Int, dim)
    # Use the approximate formula:
    #    prob_size ~ num_elem * order^dim
    num_elem = div(prob_size, order^dim)
    s = 0 # find s: num_elem/2 < 2^s <= num_elem
    while num_elem > 1
        num_elem = div(num_elem, 2)
        s += 1
    end
    r = s%dim
    for d = 1:dim
        sd = div(s, dim)
        if r > 0
            sd += 1
            r -= 1
        end
        dims[d] = 2^sd
    end
    dims
end

function build_cartesian_restriction_c(
    ceed,
    dim,
    nxyz,
    order,
    ncomp,
    num_qpts;
    form_strided=false,
)
    p = order
    pp1 = p + 1
    nnodes = pp1^dim # number of scal. nodes per element
    elem_qpts = num_qpts^dim # number of qpts per element
    num_elem = 1
    scalar_size = 1

    nd = p*nxyz .+ 1
    num_elem = prod(nxyz)
    scalar_size = prod(nd)
    size = scalar_size*ncomp

    # elem:         0             1                 n-1
    #        |---*-...-*---|---*-...-*---|- ... -|--...--|
    # nnodes:   0   1    p-1  p  p+1       2*p             n*p

    el_nodes = zeros(C.CeedInt, num_elem*nnodes)
    exyz = zeros(Int, dim)
    @inbounds for e = 0:(num_elem-1)
        re = e
        for d = 1:dim
            exyz[d] = re%nxyz[d]
            re = div(re, nxyz[d])
        end
        for lnodes = 0:(nnodes-1)
            gnodes = 0
            gnodes_stride = 1
            rnodes = lnodes
            for d = 1:dim
                gnodes += (exyz[d]*p + rnodes%pp1)*gnodes_stride
                gnodes_stride *= nd[d]
                rnodes = div(rnodes, pp1)
            end
            el_nodes[e*nnodes+lnodes+1] = gnodes
        end
    end

    restr = Ref{C.CeedElemRestriction}()
    C.CeedElemRestrictionCreate(
        ceed[],
        num_elem,
        nnodes,
        ncomp,
        scalar_size,
        ncomp*scalar_size,
        C.CEED_MEM_HOST,
        C.CEED_COPY_VALUES,
        el_nodes,
        restr,
    )
    if form_strided
        restr_i = Ref{C.CeedElemRestriction}()
        err = C.CeedElemRestrictionCreateStrided(
            ceed[],
            num_elem,
            elem_qpts,
            ncomp,
            ncomp*elem_qpts*num_elem,
            C.CEED_STRIDES_BACKEND[],
            restr_i,
        )
        return size, restr, restr_i
    else
        return size, restr
    end
end

function set_cartesian_mesh_coords_c(dim, nxyz, mesh_order, mesh_coords)
    p = mesh_order
    nd = p*nxyz .+ 1
    num_elem = prod(nxyz)
    scalar_size = prod(nd)

    coords_ref = Ref{Ptr{C.CeedScalar}}()
    C.CeedVectorGetArray(mesh_coords[], C.CEED_MEM_HOST, coords_ref)
    coords = unsafe_wrap(Array, coords_ref[], scalar_size*dim)

    nodes = zeros(C.CeedScalar, p + 1)
    # The H1 basis uses Lobatto quadrature points as nodes.
    C.CeedLobattoQuadrature(p + 1, nodes, C_NULL) # nodes are in [-1,1]
    nodes = 0.5 .+ 0.5*nodes
    for gsnodes = 0:(scalar_size-1)
        rnodes = gsnodes
        for d = 1:dim
            d1d = rnodes%nd[d]
            coords[gsnodes+scalar_size*(d-1)+1] = (div(d1d, p) + nodes[d1d%p+1])/nxyz[d]
            rnodes = div(rnodes, nd[d])
        end
    end
    C.CeedVectorRestoreArray(mesh_coords[], coords_ref)
end

function transform_mesh_coords_c(dim, mesh_size, mesh_coords)
    coords_ref = Ref{Ptr{C.CeedScalar}}()
    C.CeedVectorGetArray(mesh_coords[], C.CEED_MEM_HOST, coords_ref)
    coords = unsafe_wrap(Array, coords_ref[], mesh_size)

    if dim == 1
        for i = 1:mesh_size
            # map [0,1] to [0,1] varying the mesh density
            coords[i] = 0.5 + 1.0/sqrt(3.0)*sin((2.0/3.0)*pi*(coords[i] - 0.5))
        end
        exact_volume = 1
    else
        num_nodes = div(mesh_size, dim)
        for i = 1:num_nodes
            # map (x,y) from [0,1]x[0,1] to the quarter annulus with polar
            # coordinates, (r,phi) in [1,2]x[0,pi/2] with area = 3/4*pi
            u = coords[i]
            v = coords[i+num_nodes]
            u = 1.0 + u
            v = pi/2*v
            coords[i] = u*cos(v)
            coords[i+num_nodes] = u*sin(v)
        end
        exact_volume = 3.0/4.0*pi
    end

    C.CeedVectorRestoreArray(mesh_coords[], coords_ref)
    return exact_volume
end

function run_ex1_c(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
    ncompx = dim
    prob_size < 0 && (prob_size = 256*1024)

    gallery = false

    ceed = Ref{C.Ceed}()
    C.CeedInit(ceed_spec, ceed)

    mesh_basis = Ref{C.CeedBasis}()
    sol_basis = Ref{C.CeedBasis}()
    C.CeedBasisCreateTensorH1Lagrange(
        ceed[],
        dim,
        ncompx,
        mesh_order + 1,
        num_qpts,
        C.CEED_GAUSS,
        mesh_basis,
    )
    C.CeedBasisCreateTensorH1Lagrange(
        ceed[],
        dim,
        1,
        sol_order + 1,
        num_qpts,
        C.CEED_GAUSS,
        sol_basis,
    )

    # Determine the mesh size based on the given approximate problem size.
    nxyz = get_cartesian_mesh_size_c(dim, sol_order, prob_size)
    println("Mesh size: ", nxyz)

    # Build CeedElemRestriction objects describing the mesh and solution discrete
    # representations.
    mesh_size, mesh_restr =
        build_cartesian_restriction_c(ceed, dim, nxyz, mesh_order, ncompx, num_qpts)
    sol_size, sol_restr, sol_restr_i = build_cartesian_restriction_c(
        ceed,
        dim,
        nxyz,
        sol_order,
        1,
        num_qpts,
        form_strided=true,
    )
    println("Number of mesh nodes     : ", div(mesh_size, dim))
    println("Number of solution nodes : ", sol_size)

    # Create a C.CeedVector with the mesh coordinates.
    mesh_coords = Ref{C.CeedVector}()
    C.CeedVectorCreate(ceed[], mesh_size, mesh_coords)
    set_cartesian_mesh_coords_c(dim, nxyz, mesh_order, mesh_coords)
    # Apply a transformation to the mesh.
    exact_vol = transform_mesh_coords_c(dim, mesh_size, mesh_coords)

    # Create the Q-function that builds the mass operator (i.e. computes its
    # quadrature data) and set its context data.
    build_qfunc = Ref{C.CeedQFunction}()

    build_ctx = BuildContextC(dim, dim)
    qf_ctx = Ref{C.CeedQFunctionContext}()
    C.CeedQFunctionContextCreate(ceed[], qf_ctx)
    C.CeedQFunctionContextSetData(
        qf_ctx[],
        C.CEED_MEM_HOST,
        C.CEED_USE_POINTER,
        sizeof(build_ctx),
        pointer_from_objref(build_ctx),
    )

    if !gallery
        qf_build_mass = @cfunction(
            f_build_mass_c,
            C.CeedInt,
            (Ptr{Cvoid}, C.CeedInt, Ptr{Ptr{C.CeedScalar}}, Ptr{Ptr{C.CeedScalar}})
        )
        # This creates the QFunction directly.
        C.CeedQFunctionCreateInterior(ceed[], 1, qf_build_mass, "julia", build_qfunc)
        C.CeedQFunctionAddInput(build_qfunc[], "dx", ncompx*dim, C.CEED_EVAL_GRAD)
        C.CeedQFunctionAddInput(build_qfunc[], "weights", 1, C.CEED_EVAL_WEIGHT)
        C.CeedQFunctionAddOutput(build_qfunc[], "qdata", 1, C.CEED_EVAL_NONE)
        C.CeedQFunctionSetContext(build_qfunc[], qf_ctx[])
    else
        # This creates the QFunction via the gallery.
        name = "Mass$(dim)DBuild"
        C.CeedQFunctionCreateInteriorByName(ceed[], name, build_qfunc)
    end

    # Create the operator that builds the quadrature data for the mass operator.
    build_oper = Ref{C.CeedOperator}()
    C.CeedOperatorCreate(
        ceed[],
        build_qfunc[],
        C.CEED_QFUNCTION_NONE[],
        C.CEED_QFUNCTION_NONE[],
        build_oper,
    )
    C.CeedOperatorSetField(
        build_oper[],
        "dx",
        mesh_restr[],
        mesh_basis[],
        C.CEED_VECTOR_ACTIVE[],
    )
    C.CeedOperatorSetField(
        build_oper[],
        "weights",
        C.CEED_ELEMRESTRICTION_NONE[],
        mesh_basis[],
        C.CEED_VECTOR_NONE[],
    )
    C.CeedOperatorSetField(
        build_oper[],
        "qdata",
        sol_restr_i[],
        C.CEED_BASIS_COLLOCATED[],
        C.CEED_VECTOR_ACTIVE[],
    )

    # Compute the quadrature data for the mass operator.
    qdata = Ref{C.CeedVector}()
    elem_qpts = num_qpts^dim
    num_elem = prod(nxyz)
    C.CeedVectorCreate(ceed[], num_elem*elem_qpts, qdata)

    print("Computing the quadrature data for the mass operator ...")
    flush(stdout)
    GC.@preserve build_ctx C.CeedOperatorApply(
        build_oper[],
        mesh_coords[],
        qdata[],
        C.CEED_REQUEST_IMMEDIATE[],
    )
    println(" done.")

    # Create the Q-function that defines the action of the mass operator.
    apply_qfunc = Ref{C.CeedQFunction}()
    if !gallery
        qf_apply_mass = @cfunction(
            f_apply_mass_c,
            C.CeedInt,
            (Ptr{Cvoid}, C.CeedInt, Ptr{Ptr{C.CeedScalar}}, Ptr{Ptr{C.CeedScalar}})
        )
        # This creates the QFunction directly.
        C.CeedQFunctionCreateInterior(ceed[], 1, qf_apply_mass, "julia", apply_qfunc)
        C.CeedQFunctionAddInput(apply_qfunc[], "u", 1, C.CEED_EVAL_INTERP)
        C.CeedQFunctionAddInput(apply_qfunc[], "qdata", 1, C.CEED_EVAL_NONE)
        C.CeedQFunctionAddOutput(apply_qfunc[], "v", 1, C.CEED_EVAL_INTERP)
    else
        # This creates the QFunction via the gallery.
        C.CeedQFunctionCreateInteriorByName(ceed[], "MassApply", apply_qfunc)
    end

    # Create the mass operator.
    oper = Ref{C.CeedOperator}()
    C.CeedOperatorCreate(
        ceed[],
        apply_qfunc[],
        C.CEED_QFUNCTION_NONE[],
        C.CEED_QFUNCTION_NONE[],
        oper,
    )
    C.CeedOperatorSetField(oper[], "u", sol_restr[], sol_basis[], C.CEED_VECTOR_ACTIVE[])
    C.CeedOperatorSetField(
        oper[],
        "qdata",
        sol_restr_i[],
        C.CEED_BASIS_COLLOCATED[],
        qdata[],
    )
    C.CeedOperatorSetField(oper[], "v", sol_restr[], sol_basis[], C.CEED_VECTOR_ACTIVE[])

    # Compute the mesh volume using the mass operator: vol = 1^T \cdot M \cdot 1
    print("Computing the mesh volume using the formula: vol = 1^T.M.1 ...")
    flush(stdout)
    # Create auxiliary solution-size vectors.
    u = Ref{C.CeedVector}()
    v = Ref{C.CeedVector}()
    C.CeedVectorCreate(ceed[], sol_size, u)
    C.CeedVectorCreate(ceed[], sol_size, v)

    # Initialize 'u' with ones.
    C.CeedVectorSetValue(u[], 1.0)

    # Apply the mass operator: 'u' -> 'v'.
    C.CeedOperatorApply(oper[], u[], v[], C.CEED_REQUEST_IMMEDIATE[])

    # Compute and print the sum of the entries of 'v' giving the mesh volume.
    v_host_ref = Ref{Ptr{C.CeedScalar}}()
    C.CeedVectorGetArrayRead(v[], C.CEED_MEM_HOST, v_host_ref)
    v_host = unsafe_wrap(Array, v_host_ref[], sol_size)
    vol = sum(v_host)
    C.CeedVectorRestoreArrayRead(v[], v_host_ref)

    println(" done.")
    @printf("Exact mesh volume    : % .14g\n", exact_vol)
    @printf("Computed mesh volume : % .14g\n", vol)
    @printf("Volume error         : % .14g\n", vol - exact_vol)

    # Free dynamically allocated memory.
    C.CeedVectorDestroy(u)
    C.CeedVectorDestroy(v)
    C.CeedVectorDestroy(qdata)
    C.CeedVectorDestroy(mesh_coords)
    C.CeedOperatorDestroy(oper)
    C.CeedQFunctionDestroy(apply_qfunc)
    C.CeedOperatorDestroy(build_oper)
    C.CeedQFunctionDestroy(build_qfunc)
    C.CeedElemRestrictionDestroy(sol_restr)
    C.CeedElemRestrictionDestroy(mesh_restr)
    C.CeedElemRestrictionDestroy(sol_restr_i)
    C.CeedBasisDestroy(sol_basis)
    C.CeedBasisDestroy(mesh_basis)
    C.CeedDestroy(ceed)
end

run_ex1_c(
    ceed_spec="/cpu/self",
    dim=3,
    mesh_order=4,
    sol_order=4,
    num_qpts=4 + 2,
    prob_size=-1,
)

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
module Example

using ForwardDiff

greet() = print("Hello World!")

include("extra_file.jl")

export my_f, derivative

end # module


text:
abstract type SystemCategory <: PowerSystemCaseBuilderType end
struct PSYTestSystems <: SystemCategory end
struct PSITestSystems <: SystemCategory end
struct SIIPExampleSystems <: SystemCategory end
struct PSIDTestSystems <: SystemCategory end
struct PSSETestSystems <: SystemCategory end
struct MatPowerTestSystems <: SystemCategory end

const PACKAGE_DIR = joinpath(dirname(dirname(pathof(PowerSystemCaseBuilder))))

const SYSTEM_DESCRIPTORS_FILE = joinpath(PACKAGE_DIR, "src", "system_descriptor.jl")

const SERIALIZED_DIR = joinpath(PACKAGE_DIR, "data", "serialized_system")
const SERIALIZE_NOARGS_DIR = joinpath(PACKAGE_DIR, "data", "serialized_system", "NoArgs")
const SERIALIZE_FORECASTONLY_DIR =
    joinpath(PACKAGE_DIR, "data", "serialized_system", "ForecastOnly")
const SERIALIZE_RESERVEONLY_DIR =
    joinpath(PACKAGE_DIR, "data", "serialized_system", "ReserveOnly")
const SERIALIZE_FORECASTRESERVE_DIR =
    joinpath(PACKAGE_DIR, "data", "serialized_system", "ForecastReserve")

const SEARCH_DIRS = [
    SERIALIZE_NOARGS_DIR,
    SERIALIZE_FORECASTONLY_DIR,
    SERIALIZE_RESERVEONLY_DIR,
    SERIALIZE_FORECASTRESERVE_DIR,
]

const SERIALIZE_FILE_EXTENSIONS =
    [".json", "_validation_descriptors.json", "_time_series_storage.h5"]

function download_RTS(; kwargs...)
    PowerSystems.download(
        "https://github.com/GridMod/RTS-GMLC",
        "master",
        joinpath(PACKAGE_DIR, "data"),
    )
end

function download_modified_tamu_ercot_da(; kwargs...)
    directory = abspath(normpath(joinpath(PACKAGE_DIR, "data")))
    data = joinpath(directory, "tamu_ercot")
    # This is temporary place for hosting the dataset.
    data_urls = Dict(
        "DA_sys.json" => "https://www.dropbox.com/sh/uzohjqzoyinyyas/AAC40qKEowAbGax-yYiB_4wna/DA_sys.json?dl=1",
        "DA_sys_validation_descriptors.json" => "https://www.dropbox.com/sh/uzohjqzoyinyyas/AADWU21wuWW62Fl5SP4ubo8Va/DA_sys_validation_descriptors.json?dl=1",
        "DA_sys_time_series_storage.h5" => "https://www.dropbox.com/sh/uzohjqzoyinyyas/AADURazsNKxO5l4_1wBiW8qsa/DA_sys_time_series_storage.h5?dl=1",
    )
    if !isdir(data)
        @info "Downloading TAMU ERCOT dataset."
        mkpath(data)
        for (file, urls) in data_urls
            tempfilename = Base.download(urls)
            mv(tempfilename, joinpath(data, file), force = true)
        end
    end
    return data
end

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
using NeuralVerification: compute_output
using LinearAlgebra
function sample_based_bounds(network, cell, coefficients, num_samples)
    xs = sample(cell, num_samples)
    min_obj = Inf
    max_obj = -Inf
    for x in xs
        output = compute_output(network, x)
        obj = dot(output, coefficients)
        min_obj = min(min_obj, obj)
        max_obj = max(max_obj, obj)
    end
    return [min_obj, max_obj]
end


text:
struct EnumerableGather{T,S,F,I,A} <: Enumerable
    source::S
    fields::F
    indexFields::I
    savedFields::A
    key::Symbol
    value::Symbol
end

struct Not{T}
    val::T
end

function gather(source::Enumerable, args...; key::Symbol = :key, value::Symbol = :value)
    fields = fieldnames(eltype(source))
    if length(args) > 0
        indexFields_vector = Vector{Symbol}(undef, 0)
        firstArg = true
        for arg in args
            if typeof(arg) == Symbol
                push!(indexFields_vector, arg)
            else typeof(arg) == Not{Symbol}
                if firstArg
                    indexFields_vector = [a for a in fields if a != arg.val]
                else
                    indexFields_vector = [a for a in indexFields_vector if a != arg.val]
                end
            end
            firstArg = false
        end
        indexFields = tuple(indexFields_vector...)
    else
        indexFields = fields
    end

    savedFields = (n for n in fields if !(n in indexFields)) # fields that are not in `indexFields`
    savedFieldsType = (fieldtype(eltype(source), savedField) for savedField in savedFields)

    valueTypes = (fieldtype(eltype(source), indexField) for indexField in indexFields)
    valueType = reduce(promote_type, valueTypes)

    T = NamedTuple{(savedFields..., key, value), Tuple{savedFieldsType..., Symbol, valueType}}
    return EnumerableGather{T, typeof(source), typeof(fields), typeof(indexFields), typeof(savedFields)}(source,
        fields, indexFields, savedFields, key, value)
end

function Base.iterate(iter::EnumerableGather{T, S, F, I, A}) where {T, S, F, I, A}
    source_iterate = iterate(iter.source)
    if source_iterate == nothing || length(iter.indexFields) == 0
        return nothing
    end
    key = iter.indexFields[1]
    current_source_row = source_iterate[1]
    value = current_source_row[key]
    return (T((Base.map(n->current_source_row[n], iter.savedFields)..., key, value)),
        (current_source_row=current_source_row, source_state=source_iterate[2], current_index_field_index=1))
end

function Base.iterate(iter::EnumerableGather{T, S, F, I, A}, state) where {T, S, F, I, A}
    current_index_field_index = state.current_index_field_index + 1
    if current_index_field_index > length(iter.indexFields)
        source_iterate = iterate(iter.source, state.source_state)
        if source_iterate == nothing || length(iter.indexFields) == 0
            return nothing
        end
        current_index_field_index = 1
        source_state = source_iterate[2]
        current_source_row = source_iterate[1]
    else
        source_state = state.source_state
        current_source_row = state.current_source_row
    end
    key = iter.indexFields[current_index_field_index]
    value = current_source_row[key]
    return (T((Base.map(n->current_source_row[n], iter.savedFields)..., key, value)),
        (current_source_row=current_source_row, source_state=source_state, current_index_field_index=current_index_field_index))
end

function Base.eltype(iter::EnumerableGather{T, S, F, I, A}) where {T, S, F, I, A}
    return T
end


text:
# =====================================
# MNA 5
#
# system type: affine continuous system
# state dimension: 10913
# =====================================
using ReachabilityBenchmarks, MathematicalSystems, MAT, SparseArrays

function mna5_model()
    file = matopen(@relpath "mna5.mat")

    # system matrix
    A = sparse(read(file, "A"))

    # affine term
    b = sparsevec(19:27, [fill(-0.1, 5); fill(-0.2, 4)], 10913)

    # affine continuous system
    S = @system(x' = Ax + b)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
blockvalid(x, v) = count(isequal(v), x) ≤ 1

function backtrack!(x)
    pos = findfirst(isequal(0), x)
    isa(pos, Nothing) && return true
    iloc = 3div(pos[1]-1, 3) .+ (1:3)
    jloc = 3div(pos[2]-1, 3) .+ (1:3)
    for k in 1:9
        x[pos] = k
        blockvalid(view(x, pos[1], :), k) || continue
        blockvalid(view(x, :, pos[2]), k) || continue
        blockvalid(view(x, iloc, jloc), k) || continue
        backtrack!(x) && return true
    end
    x[pos] = 0
    return false
end

function ssolve(lines, i)
    t = [lines[10i-j][k] - '0' for j in 8:-1:0, k in 1:9]
    backtrack!(t)
    sum([100, 10, 1] .* t[1, 1:3])
end

lines = readlines("p096_sudoku.txt")
@time sum(ssolve(lines, i) for i in 1:50)
@time sum(ssolve(lines, i) for i in 1:50)

max_stars_repo_path: Chapter04/13. Using array views to avoid memory allocation/sudoku.jl


text:
src = [46.8523, -121.7603]
src2 = [48.7767, -121.8144]
rec = [45.5135 -122.6801; 44.0442 -123.0925; 42.3265 -122.8756]

G0 = gcdist(src, rec)                           # vec, arr
G1 = gcdist(src[1], src[2], rec)                # lat, lon, arr
G2 = gcdist(src[1], src[2], rec[1,1], rec[1,2]) # s_lat, s_lon, r_lat, r_lon
G3 = gcdist(vcat(src',src2'), rec[1,:])             # arr, arr
G4 = gcdist(vcat(src',src2'), rec[1,:])             # arr, rec[1,:]

@test G0 == G1
@test G0[1,:] == G1[1,:] == G2[1,:]
@test G3 == G4

printstyled("  phase_utils\n", color=:light_green)
GC.gc()
spad = 10.0
epad = 10.0
to = 30
src = "IRIS"
sta = "PB.B004..EH?,PB.B004..BS?,PB.B001..BS?,PB.B001..EH?"

# First, a well-formatted string
(H,R) = FDSNevq("2018-11-30T17:29:29.00", nev=1, src="IRIS")
H = H[1]

# Now, not so hot
(H,R) = FDSNevq("201103110547", mag=[3.0, 9.9], nev=1, src="IRIS")
H = H[1]

# Create channel data
s = H.ot                                      # Start time for FDSNsta is event origin time
t = u2d(d2u(s) + 3600.0)                      # End time is 60 minutes later; we'll truncate
S = FDSNsta(sta, s=s, t=t, to=to)

# Check that nothing is initially in the phase catalog
Ev = SeisEvent(hdr=H, data=S[1:1])
@test length(Ev.data.pha[1]) == 0

printstyled("    request with invalid parameter\n", color=:light_green)
redirect_stdout(out) do
  get_pha!(Ev, pha="", model="do.my.little.dance.on.the.catwalk", to=to, v=2)
end
@test length(Ev.data.pha[1]) == 0

printstyled("    request with user-specified phase list\n", color=:light_green)
Ev = SeisEvent(hdr=H, data=S[1:1])
get_pha!(Ev, pha="P,S", to=to)
@test length(Ev.data.pha[1]) == 2
println("")
show_phases(Ev.data.pha[1])
println("")

printstyled("    request styles that return all phases\n", color=:light_green)
Ev = SeisEvent(hdr=H, data=S[1:1])
get_pha!(Ev, pha="", to=to)
pcat1 = Ev.data.pha[1]

Ev = SeisEvent(hdr=H, data=S[1:1])
get_pha!(Ev, pha="all", to=to)
pcat2 = Ev.data.pha[1]

Ev = SeisEvent(hdr=H, data=S[1:1])
get_pha!(Ev, pha="ttall", to=to)
pcat3 = Ev.data.pha[1]

@test pcat1 == pcat2 == pcat3

# This should work for all stations, yielding only the S time
printstyled("    multi-channel request\n", color=:light_green)
J = findall([endswith(id, "EHZ") for id in S.id])

printstyled("      default phase\n", color=:light_green)
Ev = SeisEvent(hdr=H, data=S[J])
get_pha!(Ev, pha="S", to=to)
for i = 1:Ev.data.n
  @test length(Ev.data.pha[i]) == 1
  @test haskey(Ev.data.pha[i], "S")
end

printstyled("      user-specified phase list\n", color=:light_green)
Ev = SeisEvent(hdr=H, data=S[J])
get_pha!(Ev, pha="pP,PP", to=to)
for i = 1:Ev.data.n
  @test length(Ev.data.pha[i]) == 2
  @test haskey(Ev.data.pha[i], "PP")
  @test haskey(Ev.data.pha[i], "pP")
end

# This should work for all stations, yielding all times
printstyled("      all phases\n", color=:light_green)
Ev = SeisEvent(hdr=H, data=S[J])
get_pha!(Ev, pha="all", to=to)
for i = 1:Ev.data.n
  @test length(Ev.data.pha[i]) > 1
  @test haskey(Ev.data.pha[i], "P")
  @test haskey(Ev.data.pha[i], "S")
  @test haskey(Ev.data.pha[i], "pP")
  @test haskey(Ev.data.pha[i], "PP")
end
println("")
show_phases(Ev.data.pha[1])
println("")


text:
mutable struct DataMMM
  n::Int
  hc::Vector{Float64} # nominal cost
  d::Vector{Float64} # deviation

  C::Int
  w::Vector{Int}

  from::Vector{Int}
  to::Vector{Int}
  δ⁺::Array{Set,1}
  δ⁻::Array{Set,1}
  nV::Int
  A::Matrix{Int}
  s::Int
  t::Int

  DataMMM() = new()
end

mutable struct solution
  TTime::Int64
  Tsol::Int64
  LP::Float64
  opt::Float64
  n_sol::Int64
  n_sol_red::Int64

  solution() = new()
end

function read_dataMMM(instance::AbstractString, number::Int)
  data = DataMMM()
  f = open(string(dirname(@__FILE__),"/data/",instance))
    ln = readline(f)
    this_instance = false
    while !this_instance
      ln = readline(f)
      parse = split(ln)
      this_number = convert(Int64,float(parse[2]))
      this_instance = (this_number == number)
      if !this_instance
        for i in 1:10 ln = readline(f) end
      end
    end
    ln = readline(f)
    parse = split(ln)
    data.n = convert(Int64,float(parse[2]))
    ln = readline(f)
    #parse = split(ln)
    #data.Γ = convert(Int64,float(parse[6]))
    ln = readline(f)
    ln = readline(f)
    data.w = extract_vector(ln)
    ln = readline(f)
    parse = split(ln)
    data.C = convert(Int64,float(parse[2]))
    ln = readline(f)
    ln = readline(f)
    data.hc = extract_vector(ln)
    ln = readline(f)
    ln = readline(f)
    data.d = extract_vector(ln)
  close(f)
  return data
end

function extract_vector(line::String)
  line = replace(line, "(", "")
  line = replace(line, ")", "")
  line = replace(line, ",", " ")
  vector = Float64[]
  parse = split(line)
  for element in parse push!(vector,convert(Float64,float(element))) end
  return vector
end

function read_dataMMM_SP(instance::AbstractString)
  data = DataMMM()
  file = readdlm(string(dirname(@__FILE__),instance))
  println("********** Read instance $(file[2,1]) ****************")
  data.nV = file[4,1]
  data.n = file[6,1]
  coord = file[8:(7+data.nV),1:2]
  data.A = file[(9+data.nV):(8+data.nV+data.n),1:data.nV]
  data.s = file[10+data.nV+data.n,1]
  data.t = file[12+data.nV+data.n,1]

  println(data.n)
  data.from = Vector{Int}(undef,data.n)
  data.to = Vector{Int}(undef,data.n)

  for i in 1:data.n
    data.from[i] = findall(data.A[i,1:data.nV] .== 1)[1]
    data.to[i] = findall(data.A[i,1:data.nV] .== -1)[1]
  end

  data.δ⁺ = Array{Set,1}(undef,data.nV)
  data.δ⁻ = Array{Set,1}(undef,data.nV)
  for i in 1:data.nV
    data.δ⁺[i] = Set()
    data.δ⁻[i] = Set()
  end
  for i in 1:data.n
    push!(data.δ⁺[data.from[i]],i)
    push!(data.δ⁻[data.to[i]],i)
  end

  #println(data.δ⁺)
  #println("from $from")
  #println("to $to")
  #println(coord)

  data.d = ones(data.n)
  data.hc = ones(data.n)
  for i in 1:data.n
    data.hc[i] = norm(coord[data.from[i],1:2]-coord[data.to[i],1:2])
    data.d[i] = data.hc[i]/2
  end
  return data
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
function Simulation(initial_time::DateTime)
  Simulation(Dates.datetime2epochms(initial_time))
end

function run(env::Environment, until::DateTime)
  run(env, Dates.datetime2epochms(until))
end

function timeout(env::Environment, delay::Period; priority::Int=0, value::Any=nothing)
  time = now(env)
  del = Dates.datetime2epochms(Dates.epochms2datetime(time)+delay)-time
  timeout(env, del; priority=priority, value=value)
end

function nowDatetime(env::Environment)
  Dates.epochms2datetime(now(env))
end


text:
include("tt_tools.jl")

"""
sparse implementation of one_body and two_body matrices
"""

function undigits(x;base=2)
    return sum([x[k]*base^(k-1) for k in 1:length(x)])
end

function sp_one_body_to_matrix(k,l,L,N)
    Is = spzeros(2^L)
    Js = spzeros(2^L)
    Vs = spzeros(2^L)
    j = 0 #counter
    for occ_list in combinations(1:L,N)
        if (l in occ_list) && ((k==l) || !(k in occ_list))
            j+=1
            x = [m in occ_list for m in 1:L]
            Js[j] = undigits(x)+1
            Vs[j] = (-1)^sum(x[1:l-1])
            x[l] = 0
            x[k] = 1
            Is[j] = undigits(x)+1
            Vs[j] *= (-1)^sum(x[1:k-1])
        end
    end
    return sparse(Is[1:j],Js[1:j],Vs[1:j],2^L,2^L)
end

#assume k<l, m<n
function sp_two_body_to_matrix(k,l,m,n,L,N)
    Is = spzeros(2^L)
    Js = spzeros(2^L)
    Vs = spzeros(2^L)
    j = 0 #counter
    for occ_list in combinations(1:L,N)
        if (m in occ_list) && (n in occ_list) && ((k==m) || (k==n) || !(k in occ_list)) && ((l==m) || (l==n) || !(l in occ_list))
            j+=1
            x = [r in occ_list for r in 1:L]
            Js[j] = undigits(x)+1
            Vs[j] = (-1)^sum(x[1:n-1])
            x[n] = 0
            Vs[j] *= (-1)^sum(x[1:m-1])
            x[m] = 0
            Vs[j] *= (-1)^sum(x[1:l-1])
            x[l] = 1
            Vs[j] *= (-1)^sum(x[1:k-1])
            x[k] = 1
            Is[j] = undigits(x)+1
        end
    end
    return sparse(Is[1:j],Js[1:j],Vs[1:j],2^L,2^L)
end

#h and V sparse
function sp_sc_to_mat(h,V,N)
    L = size(h,1)
    A = spzeros(2^L,2^L)
    for i in findall(!iszero,h)
        H = sp_one_body_to_matrix(i[1],i[2],L,N)
        A += h[i]*(H+H')
    end
    for i in findall(!iszero,V)
        H = sp_two_body_to_matrix(i[1],i[2],i[3],i[4],L,N)
        A += V[i]*(H+H')
    end
    return dropzeros!(A)
end

function ham_to_mpo(h,V,N)
    L = size(h,1)
    H = sp_sc_to_mat(h,V,N)
    dim = 2*ones(Int,L)
    return tto_spdecomp(sparsetensor_mat(dim,dim,H),L)
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
#### Function to define parameters

##### FUNCTION to create tuple of parameters needed for dBdt

# # assign a value for Z, I_K and T as it would be in a loop, just to check if ModelParameters works to generate parameters
# Z = 10
# I_K = 1
# T = 293.15
function ModelParameters(param, T, I_K, Z)
    ### Calculate body masses
    M = BodyMasses(Z=Z)

    ### Calculate carrying capacity
    K = carryingcapacity_BA(I_K, param, M, T)
    #@assert K[2] == 0 # defensive programming
    #@assert K[3] == 0

    ### Calculate producer growth rate
    r = growth_BA(param, M, T)
    @assert r[2] == 0.0 # defensive programming
    @assert r[3] == 0.0

    ### Calculate metabolism
    x = metabolism_BA(param, M, T)
    @assert x[1] == 0.0

    ### Calculate maximum ingestion
    y = max_ingestion_BA(param, M, Z, T)

    ### Calculate half saturation density
    B0 = half_saturation_BA(param, M, Z, T)

    # Combine parameters into Tuple
    ModelParameters = (M = M, K = K, r = r, x = x, y = y, B0 = B0)
    return ModelParameters
end

#### NEW ####
# Add new argument C for chemical concentration

function ModelParameters_CG(param, T, I_K, Z, dose, target_rate::String ; slope = 5.0, lower = 0.0, xmid = 50.0)
    ### Calculate body masses
    M = BodyMasses(Z=Z)

    ### Calculate effect proportion from dose
    response = log_logistic(dose = dose, slope = slope, lower = lower, xmid = xmid)

    ### Calculate carrying capacity
    K = carryingcapacity_BA(I_K, param, M, T)
    #@assert K[2] == 0 # defensive programming
    #@assert K[3] == 0

    ### Calculate producer growth rate
    r = growth_BA(param, M, T)
    @assert r[2] == 0.0 # defensive programming
    @assert r[3] == 0.0

    ### Calculate metabolism
    x = metabolism_BA(param, M, T)
    @assert x[1] == 0.0

    ### Calculate maximum ingestion
    y = max_ingestion_BA(param, M, Z, T)

    ### Calculate half saturation density
    B0 = half_saturation_BA(param, M, Z, T)

    ### Alter target rate in proportion to dose response
    if target_rate == "r" # growth
        r = response .*  r
    end

    if target_rate == "x" # metabolism
        x = response .* x
    end

    if target_rate == "y" # max consumption
        y = response .* y
    end

    if target_rate == "B0" # half saturation density
        B0 = response .* B0
    end

    if target_rate == "FR" # functional response term
        y = y
        B0 = B0
    end

    # Combine parameters into Tuple
    ModelParameters = (M = M, K = K, r = r, x = x, y = y, B0 = B0, dose = dose, slope = slope, lower = lower, xmid = xmid, response = response)
    return ModelParameters
end

#dose = 50.0
# ModelParameters_CG(param, T, I_K, Z, 30.0, "y", slope = 5.0, lower = 0.0, xmid = 50.0)
# ModelParameters_CG(param, T, I_K, Z, 30.0, "B0", slope = 5.0, lower = 0.0, xmid = 50.0)
p = ModelParameters_CG(param, T, I_K, Z, 30.0, "x", slope = 5.0, lower = 0.0, xmid = 50.0)

#log_logistic(dose = dose, slope = 5.0, lower = 0.0, xmid = 50.0)
# p = ModelParameters(param, T, I_K, Z)
# typeof(p)
# p.B0

max_stars_repo_name: hanamayall/simple-feeding-DR

max_issues_repo_name: hanamayall/simple-feeding-DR

max_forks_repo_name: hanamayall/simple-feeding-DR


text:
"""
    EmptyDict() :: AbstractDict{Union{},Union{}}
    EmptyDict(itr) :: AbstractDict
    EmptyDict{K,V}() :: AbstractDict{K,V}
    EmptyDict{K,V}(itr) :: AbstractDict{K,V}
"""
EmptyDict

"""
    SingletonDict(k => v) :: AbstractDict
    SingletonDict(itr) :: AbstractDict
    SingletonDict{K,V}(k => v) :: AbstractDict{K,V}
    SingletonDict{K,V}(itr) :: AbstractDict{K,V}
"""
SingletonDict

abstract type AbstractMicroDict{K,V,L} <: AbstractDict{K,V} end
abstract type AbstractEmptyDict{K,V,L} <: AbstractMicroDict{K,V,L} end
abstract type AbstractSingletonDict{K,V,L} <: AbstractMicroDict{K,V,L} end

struct EmptyDict{K,V,L} <: AbstractEmptyDict{K,V,L} end
EmptyDict() = EmptyDict{Union{},Union{},Dict}()
EmptyDict{K,V}() where {K,V} = EmptyDict{K,V,Dict}()
EmptyDict(::Type{Pair{K,V}}) where {K,V} = EmptyDict{K,V,Dict}()

struct SingletonDict{K,V,L} <: AbstractSingletonDict{K,V,L}
    k::K
    v::V
end
SingletonDict(kv::Pair{K,V}) where {K,V} = SingletonDict{K,V,Dict}(first(kv), last(kv))
SingletonDict{K,V}((k, v)::Pair) where {K,V} =
    SingletonDict{K,V,Dict}(convert(K, k)::K, convert(V, v)::V)

@inline getvalue(d::SingletonDict) = d.k => d.v
@inline upcast(d::AbstractMicroDict{K,V,L}) where {K,V,L} = L(d)::AbstractDict{K,V}

Base.length(::AbstractEmptyDict) = 0
@inline Base.get(::AbstractEmptyDict, _, default) = default
@inline Base.getindex(::AbstractEmptyDict, k) = throw(KeyError(k))
@inline Base.iterate(::AbstractEmptyDict) = nothing

Base.length(::AbstractSingletonDict) = 1
@inline Base.get(d::AbstractSingletonDict, k, default) =
    isequal(first(getvalue(d)), k) ? last(getvalue(d)) : default
@inline function Base.getindex(d::AbstractSingletonDict, k)
    @boundscheck isequal(k, first(getvalue(d))) || throw(KeyError(k))
    return last(getvalue(d))
end
@inline Base.iterate(d::AbstractSingletonDict) = (getvalue(d), nothing)
@inline Base.iterate(d::AbstractSingletonDict, ::Nothing) = nothing

emptyshim(::Type{L}) where {L<:AbstractDict} = EmptyDict{Union{},Union{},L}()
emptyshim(::Type{L}, ::Type{Pair{K,V}}) where {K,V,L<:AbstractDict} = EmptyDict{K,V,L}()
singletonshim(::Type{L}, kv::Pair{K,V}) where {K,V,L<:AbstractDict} =
    SingletonDict{K,V,L}(first(kv), last(kv))

function EmptyDict(itr)
    validate_empty_iterator(itr, :EmptyDict)
    if Base.IteratorEltype(itr) isa Base.HasEltype
        T = eltype(itr)
        if T === Union{}
            return EmptyDict()
        elseif T <: Pair
            return EmptyDict(T)
        end
        throw(ArgumentError("Unexpected `eltype`: `$T`"))
    else
        return EmptyDict()
    end
end

function EmptyDict{K,V}(itr) where {K,V}
    validate_empty_iterator(itr, :EmptyDict)
    return EmptyDict{K,V}()
end

function SingletonDict(itr)
    kv = validate_singleton_iterator(itr, :SingletonDict)
    safe_length(kv) == 2 || throw(ArgumentError(string(
        "Given iterator does not contain a pair. ",
        "`SingletonDict` requires an iterator of exactly one pair.",
    )))
    k, v = kv
    return SingletonDict(k => v)
end

function SingletonDict{K,V}(itr) where {K,V}
    kv = validate_singleton_iterator(itr, :SingletonDict)
    safe_length(kv) == 2 || throw(ArgumentError(string(
        "Given iterator does not contain a pair. ",
        "`SingletonDict` requires an iterator of exactly one pair.",
    )))
    k, v = kv
    return SingletonDict{K,V}(k => v)
end

function Base.showarg(io::IO, d::EmptyDict{<:Any,<:Any,Dict}, toplevel::Bool)
    @nospecialize
    print(io, "EmptyDict")
    toplevel && print(io, '{', keytype(d), ',', valtype(d), '}')
end

function Base.showarg(io::IO, d::SingletonDict{<:Any,<:Any,Dict}, toplevel::Bool)
    @nospecialize
    print(io, "SingletonDict")
    toplevel && print(io, '{', keytype(d), ',', valtype(d), '}')
end

function Base.show(io::IO, ::MIME"text/plain", d::EmptyDict{<:Any,<:Any,Dict})
    @nospecialize
    Base.showarg(io, d, true)
    print(io, "()")
end


text:
#
# Autogenerated by Thrift Compiler (0.11.0)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

struct _enum__Type
  BOOLEAN::Int32
  INT32::Int32
  INT64::Int32
  INT96::Int32
  FLOAT::Int32
  DOUBLE::Int32
  BYTE_ARRAY::Int32
  FIXED_LEN_BYTE_ARRAY::Int32
end
const _Type = _enum__Type(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

struct _enum_ConvertedType
  UTF8::Int32
  MAP::Int32
  MAP_KEY_VALUE::Int32
  LIST::Int32
  ENUM::Int32
  DECIMAL::Int32
  DATE::Int32
  TIME_MILLIS::Int32
  TIME_MICROS::Int32
  TIMESTAMP_MILLIS::Int32
  TIMESTAMP_MICROS::Int32
  UINT_8::Int32
  UINT_16::Int32
  UINT_32::Int32
  UINT_64::Int32
  INT_8::Int32
  INT_16::Int32
  INT_32::Int32
  INT_64::Int32
  JSON::Int32
  BSON::Int32
  INTERVAL::Int32
end
const ConvertedType = _enum_ConvertedType(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7), Int32(8), Int32(9), Int32(10), Int32(11), Int32(12), Int32(13), Int32(14), Int32(15), Int32(16), Int32(17), Int32(18), Int32(19), Int32(20), Int32(21))

struct _enum_FieldRepetitionType
  REQUIRED::Int32
  OPTIONAL::Int32
  REPEATED::Int32
end
const FieldRepetitionType = _enum_FieldRepetitionType(Int32(0), Int32(1), Int32(2))

struct _enum_Encoding
  PLAIN::Int32
  PLAIN_DICTIONARY::Int32
  RLE::Int32
  BIT_PACKED::Int32
  DELTA_BINARY_PACKED::Int32
  DELTA_LENGTH_BYTE_ARRAY::Int32
  DELTA_BYTE_ARRAY::Int32
  RLE_DICTIONARY::Int32
end
const Encoding = _enum_Encoding(Int32(0), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7), Int32(8))

struct _enum_CompressionCodec
  UNCOMPRESSED::Int32
  SNAPPY::Int32
  GZIP::Int32
  LZO::Int32
  BROTLI::Int32
  LZ4::Int32
  ZSTD::Int32
end
const CompressionCodec = _enum_CompressionCodec(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6))

struct _enum_PageType
  DATA_PAGE::Int32
  INDEX_PAGE::Int32
  DICTIONARY_PAGE::Int32
  DATA_PAGE_V2::Int32
end
const PageType = _enum_PageType(Int32(0), Int32(1), Int32(2), Int32(3))

struct _enum_BoundaryOrder
  UNORDERED::Int32
  ASCENDING::Int32
  DESCENDING::Int32
end
const BoundaryOrder = _enum_BoundaryOrder(Int32(0), Int32(1), Int32(2))

mutable struct Statistics <: Thrift.TMsg
  max::Vector{UInt8}
  min::Vector{UInt8}
  null_count::Int64
  distinct_count::Int64
  max_value::Vector{UInt8}
  min_value::Vector{UInt8}
  Statistics() = (o=new(); fillunset(o); o)
end # mutable struct Statistics
meta(t::Type{Statistics}) = meta(t, Symbol[:max,:min,:null_count,:distinct_count,:max_value,:min_value], Int[], Dict{Symbol,Any}())

mutable struct StringType <: Thrift.TMsg
end # mutable struct StringType

mutable struct UUIDType <: Thrift.TMsg
end # mutable struct UUIDType

mutable struct MapType <: Thrift.TMsg
end # mutable struct MapType

mutable struct ListType <: Thrift.TMsg
end # mutable struct ListType

mutable struct EnumType <: Thrift.TMsg
end # mutable struct EnumType

mutable struct DateType <: Thrift.TMsg
end # mutable struct DateType

mutable struct NullType <: Thrift.TMsg
end # mutable struct NullType

mutable struct DecimalType <: Thrift.TMsg
  scale::Int32
  precision::Int32
  DecimalType() = (o=new(); fillunset(o); o)
end # mutable struct DecimalType

mutable struct MilliSeconds <: Thrift.TMsg
end # mutable struct MilliSeconds

mutable struct MicroSeconds <: Thrift.TMsg
end # mutable struct MicroSeconds

mutable struct TimeUnit <: Thrift.TMsg
  MILLIS::MilliSeconds
  MICROS::MicroSeconds
  TimeUnit() = (o=new(); fillunset(o); o)
end # mutable struct TimeUnit
meta(t::Type{TimeUnit}) = meta(t, Symbol[:MILLIS,:MICROS], Int[], Dict{Symbol,Any}())

mutable struct TimestampType <: Thrift.TMsg
  isAdjustedToUTC::Bool
  unit::TimeUnit
  TimestampType() = (o=new(); fillunset(o); o)
end # mutable struct TimestampType

mutable struct TimeType <: Thrift.TMsg
  isAdjustedToUTC::Bool
  unit::TimeUnit
  TimeType() = (o=new(); fillunset(o); o)
end # mutable struct TimeType

mutable struct IntType <: Thrift.TMsg
  bitWidth::UInt8
  isSigned::Bool
  IntType() = (o=new(); fillunset(o); o)
end # mutable struct IntType

mutable struct JsonType <: Thrift.TMsg
end # mutable struct JsonType

mutable struct BsonType <: Thrift.TMsg
end # mutable struct BsonType

mutable struct LogicalType <: Thrift.TMsg
  STRING::StringType
  MAP::MapType
  LIST::ListType
  ENUM::EnumType
  DECIMAL::DecimalType
  DATE::DateType
  TIME::TimeType
  TIMESTAMP::TimestampType
  INTEGER::IntType
  UNKNOWN::NullType
  JSON::JsonType
  BSON::BsonType
  LogicalType() = (o=new(); fillunset(o); o)
end # mutable struct LogicalType
meta(t::Type{LogicalType}) = meta(t, Symbol[:STRING,:MAP,:LIST,:ENUM,:DECIMAL,:DATE,:TIME,:TIMESTAMP,:INTEGER,:UNKNOWN,:JSON,:BSON], Int[1,2,3,4,5,6,7,8,10,11,12,13], Dict{Symbol,Any}())

mutable struct SchemaElement <: Thrift.TMsg
  _type::Int32
  type_length::Int32
  repetition_type::Int32
  name::String
  num_children::Int32
  converted_type::Int32
  scale::Int32
  precision::Int32
  field_id::Int32
  logicalType::LogicalType
  SchemaElement() = (o=new(); fillunset(o); o)
end # mutable struct SchemaElement
meta(t::Type{SchemaElement}) = meta(t, Symbol[:_type,:type_length,:repetition_type,:num_children,:converted_type,:scale,:precision,:field_id,:logicalType], Int[], Dict{Symbol,Any}())

mutable struct DataPageHeader <: Thrift.TMsg
  num_values::Int32
  encoding::Int32
  definition_level_encoding::Int32
  repetition_level_encoding::Int32
  statistics::Statistics
  DataPageHeader() = (o=new(); fillunset(o); o)
end # mutable struct DataPageHeader
meta(t::Type{DataPageHeader}) = meta(t, Symbol[:statistics], Int[], Dict{Symbol,Any}())

mutable struct IndexPageHeader <: Thrift.TMsg
end # mutable struct IndexPageHeader

mutable struct DictionaryPageHeader <: Thrift.TMsg
  num_values::Int32
  encoding::Int32
  is_sorted::Bool
  DictionaryPageHeader() = (o=new(); fillunset(o); o)
end # mutable struct DictionaryPageHeader
meta(t::Type{DictionaryPageHeader}) = meta(t, Symbol[:is_sorted], Int[], Dict{Symbol,Any}())

mutable struct DataPageHeaderV2 <: Thrift.TMsg
  num_values::Int32
  num_nulls::Int32
  num_rows::Int32
  encoding::Int32
  definition_levels_byte_length::Int32
  repetition_levels_byte_length::Int32
  is_compressed::Bool
  statistics::Statistics
  DataPageHeaderV2() = (o=new(); fillunset(o); o)
end # mutable struct DataPageHeaderV2
meta(t::Type{DataPageHeaderV2}) = meta(t, Symbol[:is_compressed,:statistics], Int[], Dict{Symbol,Any}(:is_compressed => true))

mutable struct PageHeader <: Thrift.TMsg
  _type::Int32
  uncompressed_page_size::Int32
  compressed_page_size::Int32
  crc::Int32
  data_page_header::DataPageHeader
  index_page_header::IndexPageHeader
  dictionary_page_header::DictionaryPageHeader
  data_page_header_v2::DataPageHeaderV2
  PageHeader() = (o=new(); fillunset(o); o)
end # mutable struct PageHeader
meta(t::Type{PageHeader}) = meta(t, Symbol[:crc,:data_page_header,:index_page_header,:dictionary_page_header,:data_page_header_v2], Int[], Dict{Symbol,Any}())

mutable struct KeyValue <: Thrift.TMsg
  key::String
  value::String
  KeyValue() = (o=new(); fillunset(o); o)
end # mutable struct KeyValue
meta(t::Type{KeyValue}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

mutable struct SortingColumn <: Thrift.TMsg
  column_idx::Int32
  descending::Bool
  nulls_first::Bool
  SortingColumn() = (o=new(); fillunset(o); o)
end # mutable struct SortingColumn

mutable struct PageEncodingStats <: Thrift.TMsg
  page_type::Int32
  encoding::Int32
  count::Int32
  PageEncodingStats() = (o=new(); fillunset(o); o)
end # mutable struct PageEncodingStats

mutable struct ColumnMetaData <: Thrift.TMsg
  _type::Int32
  encodings::Vector{Int32}
  path_in_schema::Vector{String}
  codec::Int32
  num_values::Int64
  total_uncompressed_size::Int64
  total_compressed_size::Int64
  key_value_metadata::Vector{KeyValue}
  data_page_offset::Int64
  index_page_offset::Int64
  dictionary_page_offset::Int64
  statistics::Statistics
  encoding_stats::Vector{PageEncodingStats}
  ColumnMetaData() = (o=new(); fillunset(o); o)
end # mutable struct ColumnMetaData
meta(t::Type{ColumnMetaData}) = meta(t, Symbol[:key_value_metadata,:index_page_offset,:dictionary_page_offset,:statistics,:encoding_stats], Int[], Dict{Symbol,Any}())

mutable struct EncryptionWithFooterKey <: Thrift.TMsg
end # mutable struct EncryptionWithFooterKey

mutable struct EncryptionWithColumnKey <: Thrift.TMsg
  path_in_schema::Vector{String}
  column_key_metadata::Vector{UInt8}
  EncryptionWithColumnKey() = (o=new(); fillunset(o); o)
end # mutable struct EncryptionWithColumnKey
meta(t::Type{EncryptionWithColumnKey}) = meta(t, Symbol[:column_key_metadata], Int[], Dict{Symbol,Any}())

mutable struct ColumnCryptoMetaData <: Thrift.TMsg
  ENCRYPTION_WITH_FOOTER_KEY::EncryptionWithFooterKey
  ENCRYPTION_WITH_COLUMN_KEY::EncryptionWithColumnKey
  ColumnCryptoMetaData() = (o=new(); fillunset(o); o)
end # mutable struct ColumnCryptoMetaData
meta(t::Type{ColumnCryptoMetaData}) = meta(t, Symbol[:ENCRYPTION_WITH_FOOTER_KEY,:ENCRYPTION_WITH_COLUMN_KEY], Int[], Dict{Symbol,Any}())

mutable struct ColumnChunk <: Thrift.TMsg
  file_path::String
  file_offset::Int64
  meta_data::ColumnMetaData
  offset_index_offset::Int64
  offset_index_length::Int32
  column_index_offset::Int64
  column_index_length::Int32
  crypto_meta_data::ColumnCryptoMetaData
  ColumnChunk() = (o=new(); fillunset(o); o)
end # mutable struct ColumnChunk
meta(t::Type{ColumnChunk}) = meta(t, Symbol[:file_path,:meta_data,:offset_index_offset,:offset_index_length,:column_index_offset,:column_index_length,:crypto_meta_data], Int[], Dict{Symbol,Any}())

mutable struct RowGroup <: Thrift.TMsg
  columns::Vector{ColumnChunk}
  total_byte_size::Int64
  num_rows::Int64
  sorting_columns::Vector{SortingColumn}
  RowGroup() = (o=new(); fillunset(o); o)
end # mutable struct RowGroup
meta(t::Type{RowGroup}) = meta(t, Symbol[:sorting_columns], Int[], Dict{Symbol,Any}())

mutable struct TypeDefinedOrder <: Thrift.TMsg
end # mutable struct TypeDefinedOrder

mutable struct ColumnOrder <: Thrift.TMsg
  TYPE_ORDER::TypeDefinedOrder
  ColumnOrder() = (o=new(); fillunset(o); o)
end # mutable struct ColumnOrder
meta(t::Type{ColumnOrder}) = meta(t, Symbol[:TYPE_ORDER], Int[], Dict{Symbol,Any}())

mutable struct PageLocation <: Thrift.TMsg
  offset::Int64
  compressed_page_size::Int32
  first_row_index::Int64
  PageLocation() = (o=new(); fillunset(o); o)
end # mutable struct PageLocation

mutable struct OffsetIndex <: Thrift.TMsg
  page_locations::Vector{PageLocation}
  OffsetIndex() = (o=new(); fillunset(o); o)
end # mutable struct OffsetIndex

mutable struct ColumnIndex <: Thrift.TMsg
  null_pages::Vector{Bool}
  min_values::Vector{Vector{UInt8}}
  max_values::Vector{Vector{UInt8}}
  boundary_order::Int32
  null_counts::Vector{Int64}
  ColumnIndex() = (o=new(); fillunset(o); o)
end # mutable struct ColumnIndex
meta(t::Type{ColumnIndex}) = meta(t, Symbol[:null_counts], Int[], Dict{Symbol,Any}())

mutable struct FileMetaData <: Thrift.TMsg
  version::Int32
  schema::Vector{SchemaElement}
  num_rows::Int64
  row_groups::Vector{RowGroup}
  key_value_metadata::Vector{KeyValue}
  created_by::String
  column_orders::Vector{ColumnOrder}
  FileMetaData() = (o=new(); fillunset(o); o)
end # mutable struct FileMetaData
meta(t::Type{FileMetaData}) = meta(t, Symbol[:key_value_metadata,:created_by,:column_orders], Int[], Dict{Symbol,Any}())

mutable struct AesGcmV1 <: Thrift.TMsg
  aad_metadata::Vector{UInt8}
  AesGcmV1() = (o=new(); fillunset(o); o)
end # mutable struct AesGcmV1
meta(t::Type{AesGcmV1}) = meta(t, Symbol[:aad_metadata], Int[], Dict{Symbol,Any}())

mutable struct AesGcmCtrV1 <: Thrift.TMsg
  aad_metadata::Vector{UInt8}
  AesGcmCtrV1() = (o=new(); fillunset(o); o)
end # mutable struct AesGcmCtrV1
meta(t::Type{AesGcmCtrV1}) = meta(t, Symbol[:aad_metadata], Int[], Dict{Symbol,Any}())

mutable struct EncryptionAlgorithm <: Thrift.TMsg
  AES_GCM_V1::AesGcmV1
  AES_GCM_CTR_V1::AesGcmCtrV1
  EncryptionAlgorithm() = (o=new(); fillunset(o); o)
end # mutable struct EncryptionAlgorithm
meta(t::Type{EncryptionAlgorithm}) = meta(t, Symbol[:AES_GCM_V1,:AES_GCM_CTR_V1], Int[], Dict{Symbol,Any}())

mutable struct FileCryptoMetaData <: Thrift.TMsg
  encryption_algorithm::EncryptionAlgorithm
  encrypted_footer::Bool
  footer_key_metadata::Vector{UInt8}
  footer_offset::Int64
  iv_prefix::Vector{UInt8}
  FileCryptoMetaData() = (o=new(); fillunset(o); o)
end # mutable struct FileCryptoMetaData
meta(t::Type{FileCryptoMetaData}) = meta(t, Symbol[:footer_key_metadata,:iv_prefix], Int[], Dict{Symbol,Any}())


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

mutable struct DateAfterModification <: SwaggerModel
    daysAfterModificationGreaterThan::Any # spec type: Union{ Nothing, Float32 } # spec name: daysAfterModificationGreaterThan

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

const _property_map_DateAfterModification = Dict{Symbol,Symbol}(Symbol("daysAfterModificationGreaterThan")=>Symbol("daysAfterModificationGreaterThan"))
const _property_types_DateAfterModification = Dict{Symbol,String}(Symbol("daysAfterModificationGreaterThan")=>"Float32")
Base.propertynames(::Type{ DateAfterModification }) = collect(keys(_property_map_DateAfterModification))
Swagger.property_type(::Type{ DateAfterModification }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DateAfterModification[name]))}
Swagger.field_name(::Type{ DateAfterModification }, property_name::Symbol) =  _property_map_DateAfterModification[property_name]

function check_required(o::DateAfterModification)
    (getproperty(o, Symbol("daysAfterModificationGreaterThan")) === nothing) && (return false)
    true
end

function validate_property(::Type{ DateAfterModification }, name::Symbol, val)
    if name === Symbol("daysAfterModificationGreaterThan")
        Swagger.validate_param(name, "DateAfterModification", :minimum, val, 0, false)
    end
end


text:
abstract Localization

include("dim_select.jl")
export iter_dim_seletct

function iter_select (ii :: Tuple{Int64, Int64}, loc :: Localization,
                      raw :: Array{Float64,2},
                      obs :: Array{Float64,1}, error_obs :: Array{Float64,1})
  (iobs, ixx ) = ii
  llrr = iter_dim_seletct(ii, loc, size(raw))
  return (obs[iobs], error_obs[iobs], raw[llrr,:], squeeze(raw[ixx,:],1))
end

function iter_refill!(ii :: Tuple{Int64, Int64}, m :: Localization,
                     raw :: Array{Float64,2}, upyy :: Array{Float64,2} )
  (iobs, ixx ) = ii
  llrr = iter_dim_seletct(ii, loc, size(raw))
  raw[llrr,:] = upyy
end

export iter_select, iter_refill!