text:
module GmshParser

export extract_mesh

function getNodes(meshString::Array{String})
    """
    takes meshString and returns all Nodes with their corresponding tags
    as 4th element in the Array element
    """

    nodeStart = findall(meshString .== "\$Nodes")[1] #\$ escape string key $
    nodeEnd = findall(meshString .== "\$EndNodes")[1]
    _, nnodes,min_node_tag,max_node_tag = parse.(Int, split(meshString[nodeStart+1]))
    @assert max_node_tag == nnodes "nodes tags are not contiguous"
    nodes = []
    node = []
    tag = []
    nodeBlock::Int = 0
    n_entity_nodes = 0
    entity_dim = 0
    entity_tag = 0
    tmp = 0
    first = true

    for nodeLine in meshString[nodeStart+2:nodeEnd-1]
        if nodeBlock == 0 && parse(Int, split(nodeLine)[4]) == 0
            continue
        elseif nodeBlock == 0 && parse(Int, split(nodeLine)[4]) != 0
            entity_dim, entity_tag, _, n_entity_nodes = parse.(Int, split(nodeLine))
            #nodeBlock = parse(Int, split(nodeLine)[4])
            nodeBlock = 2 * n_entity_nodes
            tmp = nodeBlock
            if first == false
                #block = push!.(node, tag)
                block = vcat.(tag, node)
                push!(nodes, block)
                tag = []
                node = []
            end
            first = false
        else
            if nodeBlock > tmp // 2
                push!(tag, [entity_dim, entity_tag, parse(Int, nodeLine)])
            else
                push!(node, parse.(Float64, split(nodeLine)))
            end
            nodeBlock -= 1.0
        end
    end
    # The arithmetic above does not account the last node Block
    # hence we need to explicitely append it
    #block = push!.(node, tag)
    block = vcat.(tag, node)
    push!(nodes, block)
    # flatten the array
    return collect(Iterators.flatten(nodes)) # [entDim entTag nodeTag X Y Z]
end

function map_entities_to_physical_names(meshstring)
    linestart = findall(meshstring .== "\$Entities")[1] #\$ escape string key $
    linestop = findall(meshstring .== "\$EndEntities")[1]
    phys_tag2name = map_physical_tags_to_names(meshstring)
    npoints, ncurves, nsurfs, nvols = parse.(Int, split(meshstring[linestart+1]))
    limits = cumsum((npoints, ncurves, nsurfs, nvols))
    entities_tag2name = (Dict{Int,String}(), Dict{Int,String}(), Dict{Int,String}(), Dict{Int,String}())
    dim_id = 1
    for (i,line) in enumerate(meshstring[linestart+2 : linestop-1])
        dim_id = get_dim_id_recursively(i,dim_id,limits)
        if dim_id == 1
            pointline = split(line)
            if length(pointline) > 5
                n_phys_tags = parse(Int,pointline[5])
                for itag in 1:n_phys_tags
                    entity_tag = parse(Int,pointline[1])
                    phys_tag = parse(Int,pointline[5+itag])
                    push!(entities_tag2name[dim_id], entity_tag => phys_tag2name[phys_tag])
                end
            end
        elseif dim_id >= 2
            line_vec = split(line)
            if length(line_vec) > 8
                n_phys_tags = parse(Int,line_vec[8])
                for itag in 1:n_phys_tags
                    entity_tag = parse(Int,line_vec[1])
                    phys_tag = parse(Int,line_vec[8+itag])
                    push!(entities_tag2name[dim_id], entity_tag => phys_tag2name[phys_tag])
                end
            end
        end
    end
    return entities_tag2name
end

function map_physical_tags_to_names(meshstring)
    linestart = findall(meshstring .== "\$PhysicalNames")[1] #\$ escape string key $
    linestop = findall(meshstring .== "\$EndPhysicalNames")[1]
    phys_tag2name = Dict{Int,String}()
    for line in meshstring[linestart+2 : linestop-1]
        entity_line = split(line)
        phys_tag2name[parse(Int,entity_line[2])] = entity_line[3][2:end-1]
    end
    return phys_tag2name
end

function get_dim_id_recursively(i,dim_id,limits)
    if i > limits[dim_id]
        return get_dim_id_recursively(i,dim_id + 1,limits)
    else
        return dim_id
    end
end

function getElements(meshString::Array{String})
    dim = get_dim(meshString)
    elementStart = findall(meshString .== "\$Elements")[1] #\$ escape string key $
    elementStop = findall(meshString .== "\$EndElements")[1]
    elementTotal = parse(Int, split(meshString[elementStart+1])[2]) #no. of Nodes
    nodeElements = []
    boundaryElements = []
    domainElements = []
    ele = []
    eleBlock::Int = 0
    tmp = 0
    eleDim = 0
    first = true
    eletype = 0
    entity_tag = 1

    for eleLine in meshString[elementStart+2:elementStop-1]
        if eleBlock == 0 && parse(Int, split(eleLine)[4]) == 0 # Entity definition and containing no elements
            continue
        elseif eleBlock == 0 && parse(Int, split(eleLine)[4]) != 0 # entity def and contains elements
            eleBlock = parse(Int, split(eleLine)[4]) # number of element for this entity
            tmp = eleBlock
            if first == false
                if eleDim == dim-1
                    push!(boundaryElements, pushfirst!.(ele,entity_tag))
                elseif eleDim == dim
                    push!(domainElements, ele)
                end
            end
            eleDim, entity_tag, eletype, param = parse.(Int, split(eleLine)) # dimension of the entity
            ele = []
            first = false # TODO move that before loop ends
        else
            push!(ele, parse.(Int, split(eleLine)))  # elem id + nodes global indices
            eleBlock -= 1.0
        end
    end
    if eleDim == dim - 2
        push!(boundaryElements, pushfirst!.(ele,entity_tag))
    elseif eleDim == dim - 1
        push!(boundaryElements, pushfirst!.(ele,entity_tag)) # [ entity_tag  elem id  nodes_gidx... ]
    elseif eleDim == dim
        push!(domainElements, ele) # [ elem id  nodes_idx... ]
    end

    # flatten the array
    nodeElements = collect(Iterators.flatten(nodeElements))
    boundaryElements = collect(Iterators.flatten(boundaryElements))
    domainElements = collect(Iterators.flatten(domainElements))
    return nodeElements, boundaryElements, domainElements
end

# function nodeToDataFrame(nodes)
#     df = DataFrame(tag = [], x = [], y = [], z = [])
#     for node in nodes
#         push!(df, [node[1], node[2], node[3], node[4]])
#     end
#     return df
# end

function get_facesets(meshstring, elements, el_type)
    nodeElements, boundaryElements, domainElements = elements
    n_nodeElements = size(nodeElements,1)
    n_boundaryElements = size(boundaryElements,1)
    n_domainElements, nnodes_domainElements = size(domainElements)
    nnodes_domainElements -= 1
    boundary_faces = Tuple{String,Int,Int}[]
    n_boundaries_entities = length(unique(boundaryElements[:,1]))
    #boundaries_names = get_boundaries_entities_names(meshstring, dim)
    curves_tag2name = map_entities_to_physical_names(meshstring)[2] # get this for curves entities
    for i = 1:n_boundaryElements
        boundary_entity_tag = boundaryElements[i,1]
        boundary_name = curves_tag2name[boundary_entity_tag]
        boundary_nodes = boundaryElements[i,3:end]
        for j = 1:n_domainElements
            element_nodes = domainElements[j,2:end]
            #println("element type is : $(domainElements[j,1]) $el_type")
            if issubset(boundary_nodes, element_nodes)
                faces_nodes_indices = faces_indices(el_type)
                face_id = findall([issubset(element_nodes[face_idxs], boundary_nodes) for face_idxs in faces_nodes_indices])
                @assert length(face_id) == 1 "length(face_id) = $(length(face_id))"
                push!(boundary_faces, (boundary_name, j, face_id[1])) # pushes (physical_entity_tag, elem_id, face_number)
            end # domainElements[j,1] - (n_nodeElements + n_boundaryElements)
        end
    end

    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    for name in unique([face[1] for face in boundary_faces])#unique(two_c[:,1])
        slice = filter(x->x[1]==name, boundary_faces)
        facesets[name] = Set{Tuple{Int,Int}}(getindex.(slice,[2:3]))
    end
    return facesets
end

function get_nodesets(meshstring, nodes) #[entDim entTag nodeTag X Y Z]
    nodes_in_nodalentity = [node for node in nodes if node[1] == 0]#nodes[nodes[:,1] .== 0 ,:]
    nodesets = Dict{String,Set{Int}}()
    points_tag2name = map_entities_to_physical_names(meshstring)[1]
    nnodessets = length(points_tag2name)
    for (entity_tag,entity_name) in pairs(points_tag2name)
        nodes_w_entity_tag = [node for node in nodes_in_nodalentity if node[2] == entity_tag]#[nodes_in_nodalentity[:,2] .== entity_tag,:]
        nodes_tags = Set([node[3] for node in nodes_w_entity_tag])
        push!(nodesets, entity_name => nodes_tags)
    end
    return nodesets
end

function faces_indices(c::String)
    if c in ("Line","QuadraticLine")
        return [1, 2]
    elseif c in ("Triangle","QuadraticTriangle")
        return ([1,2], [2,3], [3,1])
    elseif c in ("Quadrilateral","QuadraticQuadrilateral")
        return ([1,2], [2,3], [3,4], [4,1])
    elseif c in ("Tetrahedron","QuadraticTetrahedron")
        return ([1,2,3], [1,2,4], [2,3,4], [1,4,3])
    elseif c in ("Hexahedron","QuadraticHexahedron")
        return ([1,4,3,2], [1,2,6,5], [2,3,7,6], [3,4,8,7], [1,5,8,4], [5,6,7,8])
    else

    end
end

get_el_type_str(meshstring,higher_order_elem) = _get_el_type_str(get_dim(meshstring),length(higher_order_elem[end,2:end]))

function _get_el_type_str(dim, nnodes)
    if dim == 2
        if nnodes == 3
            el_type = "Triangle" #TRI3
        elseif nnodes == 6
            el_type = "QuadraticTriangle" #TRI3
        elseif nnodes == 4
            el_type = "Quadrilateral" #QUAD4
        elseif nnodes == 9
            el_type = "QuadraticQuadrilateral" #QUAD4
        else
            @error "element type identifier $el_type_number from .msh doesn't have any JuAFEM correspondance ."
        end
    else
        @error "dimension != 2 are not implemented yet"
    end
    return el_type
end

function get_dim(meshstring)
    linestart = findall(meshstring .== "\$Entities")[1] #\$ escape string key $
    ndim_entities = parse.(Int, split(meshstring[linestart+1]))
    domain_dim = findlast(ndim_entities .!= 0) - 1
    return domain_dim
end

function extract_mesh(meshfile)
    io = open(meshfile)
    lines = readlines(io)
    close(io)

    meshnodes = getNodes(lines)

    el_0D, el_1D, el_2D = getElements(lines)

    el_0D_mat = hcat(el_0D...)'
    el_1D_mat = hcat(el_1D...)'
    el_2D_mat = hcat(el_2D...)'
    elements_mat = (el_0D_mat, el_1D_mat, el_2D_mat)

    higher_order_el_mat = elements_mat[findlast(.!isempty.(elements_mat))]
    #println(higher_order_el_mat)
    el_type_str = get_el_type_str(lines,higher_order_el_mat)

    facesets = get_facesets(lines,elements_mat,el_type_str)
    nodesets = get_nodesets(lines, meshnodes)

    return el_type_str, higher_order_el_mat, meshnodes, facesets, nodesets

end

end # module


text:
using DEBTrait
using CSV, DataFrames, Statistics

## this is all still very clumsy

df_monomers = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/data/monomers_all.csv", DataFrame, missingstring="N/A")
chemFormBiom = [1, 1.8, 0.2, 0.5, 0, 0, 0]

n_monomers = size(df_monomers,1)
dGcox = zeros(n_monomers)
for i in 1:n_monomers
    elementstring = df_monomers.Formula[i]
    out = DEBTrait.ThermoStoichWizard.get_lambda(elementstring, chemFormBiom)
    dGcox[i] = out[2][3]
end
df_monomers.dGcox = dGcox

dGcat = zeros(n_monomers)
for i in 1:n_monomers
    elementstring = df_monomers.Formula[i]
    out = DEBTrait.ThermoStoichWizard.get_lambda(elementstring, chemFormBiom)
    dGcat[i] = out[2][5]
end
df_monomers.dGcat = dGcat

dGAn = zeros(n_monomers)
for i in 1:n_monomers
    elementstring = df_monomers.Formula[i]
    out = DEBTrait.ThermoStoichWizard.get_lambda(elementstring, chemFormBiom)
    dGAn[i] = out[2][8]
end
df_monomers.dGAn = dGAn

λ = zeros(n_monomers)
for i in 1:n_monomers
    elementstring = df_monomers.Formula[i]
    out = DEBTrait.ThermoStoichWizard.get_lambda(elementstring, chemFormBiom)
    λ[i] = out[1][1]
end
df_monomers.lambda = λ

η    = 0.43
y_DE = @. (λ*η*dGcat)/dGcox
df_monomers.y_DE = y_DE

N_C = zeros(n_monomers)
for i in 1:n_monomers
    elementstring = df_monomers.Formula[i]
    N_C[i]        = DEBTrait.ThermoStoichWizard.extract_composition(elementstring)[1]
end
df_monomers.N_C = N_C

N_N = zeros(n_monomers)
for i in 1:n_monomers
    elementstring = df_monomers.Formula[i]
    N_N[i]        = DEBTrait.ThermoStoichWizard.extract_composition(elementstring)[3]
end
df_monomers.N_N = N_N

D_S           = DEBTrait.aqueous_diffusivity(convert(Array{Float64,1}, df_monomers.Molecular_weight))

df            = DataFrame()
df.y_DE       = y_DE
df.D_S        = D_S
df.N_C        = N_C
df.N_N        = N_N
df.Ontology   = df_monomers.Ontology
CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/thermodynamic_props_all.csv", df)

############ calculate median values for monomer classes
tmp           = convert(Array{Int64,1},1:n_monomers)
sugars        = tmp[df_monomers.Ontology .== "sugar"]
organics      = tmp[df_monomers.Ontology .== "organic acid"]
aminos        = tmp[df_monomers.Ontology .== "amino acid"]
fattys        = tmp[df_monomers.Ontology .== "fatty acid"]
nucleos       = tmp[df_monomers.Ontology .== "nucleotide"]
auxins        = tmp[df_monomers.Ontology .== "auxin"]

y_DE          = vcat(median(y_DE[sugars]),median(y_DE[organics]), median(y_DE[aminos]), median(y_DE[fattys]), median(y_DE[nucleos]), median(y_DE[auxins]))
D_S           = vcat(median(D_S[sugars]),median(D_S[organics]), median(D_S[aminos]), median(D_S[fattys]), median(D_S[nucleos]), median(D_S[auxins]))
N_C           = vcat(median(N_C[sugars]),median(N_C[organics]), median(N_C[aminos]), median(N_C[fattys]), median(N_C[nucleos]), median(N_C[auxins]))
N_N           = vcat(median(N_N[sugars]),median(N_N[organics]), median(N_N[aminos]), median(N_N[fattys]), median(N_N[nucleos]), median(N_N[auxins]))

df_classes    = DataFrame()
df_classes.y_DE = y_DE
df_classes.D_S = D_S
df_classes.N_C = N_C
df_classes.N_N = N_N
df_classes.Ontology = ["sugar", "organic acid", "amino acid", "fatty acid", "nucletotide", "auxin"]

CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/thermodynamic_props_class.csv", df_classes)

ext: jl


text:
using Compat

if VERSION < v"0.7"
    Base.read(x, ::Type{String}) = readstring(x)
end

function openssl_ripemd160(x::AbstractString)
    read(pipeline(`printf $x`, `openssl ripemd160`), String)[10:end - 1]
end

function openssl_ripemd160(x::Array{UInt8, 1})
    read(pipeline(`printf $(String(x))`,
                  `openssl rmd160`),
         String)[10:end - 1]
end

function openssl_ripemd160(x::NTuple{N, UInt8}) where N
    t = String([x...])
    read(pipeline(`printf $t`, `openssl rmd160`), String)[10:end - 1]
end

function vs_openssl(x)
    bytes2hex(Ripemd.ripemd160(x)) == openssl_ripemd160(x)
end

@testset "Ripemd160 vs openssl a's" begin
    for i in 1:1000
        x = [0x61 for j in 1:i]
        @test vs_openssl(x)
    end
    for i in 1:1000
        x = *(["a" for j in 1:i]...,)
        @test vs_openssl(x)
    end
    # This takes a long time, because Julia will compile a different function
    # for each length/loop iterations
    for i in 1:10
        x = ntuple(x -> 0x61, i)
        @test vs_openssl(x)
    end
end

@testset "Ripemd160 vs openssl abc" begin
    d = Ripemd.codeunits("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789")
    for i in 1:1000
        x = [d[j % length(d) + 1] for j in 1:i]
        @test vs_openssl(x)
    end
    for i in 1:1000
        x = *([Char(d[j % length(d) + 1]) for j in 1:i]...,)
        @test vs_openssl(x)
    end
    # This takes a long time, because Julia will compile a different function
    # for each length/loop iterations
    for i in 1:10
        x = ntuple(j -> d[j % length(d) + 1], i)
        @test vs_openssl(x)
    end
end

@testset "Ripemd160 1M a's" begin

    d1 = [0x61 for i in 1:1_000_000];
    d2 = ntuple((i) -> 0x61, 1_000_000);
    r = "52783243c1697bdbe16d37f97f68f08325dc1528"

    @test bytes2hex(Ripemd.ripemd160(d1)) == r
    @test bytes2hex(Ripemd.ripemd160(d2)) == r
end

@testset "Ripemd160" begin
    @test bytes2hex(Ripemd.ripemd160("asdf")) ==
        "0ef2aed6346def670a8019e4ea42cf4c76018139"
    @test bytes2hex(Ripemd.ripemd160("")) ==
        "9c1185a5c5e9fc54612808977ee8f548b2258d31"
    @test bytes2hex(Ripemd.ripemd160("a")) ==
        "0bdc9d2d256b3ee9daae347be6f4dc835a467ffe"
    @test bytes2hex(Ripemd.ripemd160("abc")) ==
        "8eb208f7e05d987a9b044a8e98c6b087f15a0bfc"
    @test bytes2hex(Ripemd.ripemd160("message digest")) ==
        "5d0689ef49d2fae572b881b123a85ffa21595f36"
    @test bytes2hex(Ripemd.ripemd160("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789")) ==
        "b0e20b6e3116640286ed3a87a5713079b21f5189"
    @test bytes2hex(Ripemd.ripemd160("I'd hold you up to say to your mother, 'this kid's gonna be the best kid in the world. This kid's gonna be somebody better than anybody I ever knew.' And you grew up good and wonderful. It was great just watching you, every day was like a privilege. Then the time come for you to be your own man and take on the world, and you did. But somewhere along the line, you changed. You stopped being you. You let people stick a finger in your face and tell you you're no good. And when things got hard, you started looking for something to blame, like a big shadow. Let me tell you something you already know. The world ain't all sunshine and rainbows. It's a very mean and nasty place and I don't care how tough you are it will beat you to your knees and keep you there permanently if you let it. You, me, or nobody is gonna hit as hard as life. But it ain't about how hard ya hit. It's about how hard you can get hit and keep moving forward. How much you can take and keep moving forward. That's how winning is done! Now if you know what you're worth then go out and get what you're worth. But ya gotta be willing to take the hits, and not pointing fingers saying you ain't where you wanna be because of him, or her, or anybody! Cowards do that and that ain't you! You're better than that! I'm always gonna love you no matter what. No matter what happens. You're my son and you're my blood. You're the best thing in my life. But until you start believing in yourself, ya ain't gonna have a life. Don't forget to visit your mother.")) ==
        "fff55c23c197b4fded67e09424e5aef9dafad1c6"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits("asdf"))) ==
        "0ef2aed6346def670a8019e4ea42cf4c76018139"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits(""))) ==
        "9c1185a5c5e9fc54612808977ee8f548b2258d31"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits("a"))) ==
        "0bdc9d2d256b3ee9daae347be6f4dc835a467ffe"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits("abc"))) ==
        "8eb208f7e05d987a9b044a8e98c6b087f15a0bfc"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits("message digest"))) ==
        "5d0689ef49d2fae572b881b123a85ffa21595f36"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789"))) ==
        "b0e20b6e3116640286ed3a87a5713079b21f5189"
    @test bytes2hex(Ripemd.ripemd160(Ripemd.codeunits("I'd hold you up to say to your mother, 'this kid's gonna be the best kid in the world. This kid's gonna be somebody better than anybody I ever knew.' And you grew up good and wonderful. It was great just watching you, every day was like a privilege. Then the time come for you to be your own man and take on the world, and you did. But somewhere along the line, you changed. You stopped being you. You let people stick a finger in your face and tell you you're no good. And when things got hard, you started looking for something to blame, like a big shadow. Let me tell you something you already know. The world ain't all sunshine and rainbows. It's a very mean and nasty place and I don't care how tough you are it will beat you to your knees and keep you there permanently if you let it. You, me, or nobody is gonna hit as hard as life. But it ain't about how hard ya hit. It's about how hard you can get hit and keep moving forward. How much you can take and keep moving forward. That's how winning is done! Now if you know what you're worth then go out and get what you're worth. But ya gotta be willing to take the hits, and not pointing fingers saying you ain't where you wanna be because of him, or her, or anybody! Cowards do that and that ain't you! You're better than that! I'm always gonna love you no matter what. No matter what happens. You're my son and you're my blood. You're the best thing in my life. But until you start believing in yourself, ya ain't gonna have a life. Don't forget to visit your mother."))) ==
        "fff55c23c197b4fded67e09424e5aef9dafad1c6"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("asdf")...,))) ==
        "0ef2aed6346def670a8019e4ea42cf4c76018139"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("")...,))) ==
        "9c1185a5c5e9fc54612808977ee8f548b2258d31"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("a")...,))) ==
        "0bdc9d2d256b3ee9daae347be6f4dc835a467ffe"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("abc")...,))) ==
        "8eb208f7e05d987a9b044a8e98c6b087f15a0bfc"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("message digest")...,))) ==
        "5d0689ef49d2fae572b881b123a85ffa21595f36"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789")...,))) ==
        "b0e20b6e3116640286ed3a87a5713079b21f5189"
    @test bytes2hex(Ripemd.ripemd160((Ripemd.codeunits("I'd hold you up to say to your mother, 'this kid's gonna be the best kid in the world. This kid's gonna be somebody better than anybody I ever knew.' And you grew up good and wonderful. It was great just watching you, every day was like a privilege. Then the time come for you to be your own man and take on the world, and you did. But somewhere along the line, you changed. You stopped being you. You let people stick a finger in your face and tell you you're no good. And when things got hard, you started looking for something to blame, like a big shadow. Let me tell you something you already know. The world ain't all sunshine and rainbows. It's a very mean and nasty place and I don't care how tough you are it will beat you to your knees and keep you there permanently if you let it. You, me, or nobody is gonna hit as hard as life. But it ain't about how hard ya hit. It's about how hard you can get hit and keep moving forward. How much you can take and keep moving forward. That's how winning is done! Now if you know what you're worth then go out and get what you're worth. But ya gotta be willing to take the hits, and not pointing fingers saying you ain't where you wanna be because of him, or her, or anybody! Cowards do that and that ain't you! You're better than that! I'm always gonna love you no matter what. No matter what happens. You're my son and you're my blood. You're the best thing in my life. But until you start believing in yourself, ya ain't gonna have a life. Don't forget to visit your mother.")...,))) ==
        "fff55c23c197b4fded67e09424e5aef9dafad1c6"
end


text:
include("yeast_dyn.jl")
include("coop.jl")
include("shuffle.jl")
#srand(1)

## --- a (wrong?) analytical model for mean cooperativity
## --- gamma = 2.5 is the best fit to the data for n=100
## --- with basin weighted averaging - coopbasinweighted()
function coop_anal(p,a)
return p^a + (1-p)^a
end

# --- calculate average cooperativity over random networks
n = 11
n = length(selfact)
ne = 29                                # number of edges
ntotal = 50000;                        # ensemble size
listc=Float64[]                    # list of cooperativity values
lista=Float64[]                    # list of active node fractions
notconn=0;
nofp=0;
samenw=0;
for (np = 15:1:15)                     # scan np = 0 .. ne
#    interc = randomnkp(selfact,ne,np) # generate one sample with given (nkp)
    netdict = Dict();
    i=0
    while (i<ntotal)
        (intersh,selfsh) = (switchnet(interc),selfact) # preserve in/out degrees
#        (intersh,selfsh) = shufflenet(interc,selfact) # shuffle
        if (isconnected(symadj(intersh))) # check if connected
            netkey=join(sort(real(eigvals(intersh+diagm(selfsh)))))
            if (haskey(netdict,netkey)) samenw += 1  # generated before
            else
                i += 1                  # increase network count
                netdict[netkey] = 1      # record this network
                (ncop,actratio) = cooppointattr(intersh,selfsh) # get cooperativity
                if (ncop >= 0)    # will return -1 if no nontrivial fp
                    push!(listc,ncop)
                    push!(lista,actratio)
                else nofp += 1 end
            end
        else notconn += 1 end
        if i%(ntotal/100)==0 print(".") end
    end
    println("\n same network count = ",samenw)
end
(bins,vals) = hist(listc)
#display(plot(bins,vals))

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories, Catlab.Present
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq

@present ThSIR(FreeSchema) begin
    (S,I,R)::Ob
end

"""This code will work for any choice of 'type graph'"""
function update_state(state::StructACSet{S}) where S
    rules = []
    while !isempty(rules)
        state = apply_rule(state, pop!(rules))
    end
end

@acset_type SIR(ThSIR)

# infect = Rule("name", L, R)

I = @acset SIR begin I=1 end
I2 = @acset SIR begin I=2 end
SI = @acset SIR begin S=1; I=1 end
# L_infect =  homomorphism(I, SI) # equivalent to below b/c only one morphism
L_infect =  ACSetTransformation(I, SI; I=[1])
R_infect = ACSetTransformation(I, I2; I=[1])

ex_state = @acset SIR begin S=10; I=3; R=1 end

# rewrite(L_infect, R_infect, ex_state) # execute immediately

match = homomorphism(SI, ex_state)

matches = [match_1, match_2] # pointing to ex_state

"""
L    I    R
|
G <- K -> H (updated G)
   kg  kh

"""

_, kg, _, kh = rewrite_match_maps(L_infect, R_infect, match)

# pretend state_2 is a homomorphism from state_1 -> state_2

matches = [match for match in matches]

state_3 = rewrite_match(L_infect, R_infect, match_2)

# # type 2 --------------------------------------------------
# @present ThSIR2(FreeSchema) begin
#     State::Ob
#     Agents::Ob
#     state::Hom(Agents, State)
# end

# @acset_type SIR2(ThSIR2)

# state0 = vcat(fill.(1:3, [5,4,1])...)

# SIR2_state = @acset SIR2 begin State=3; Agents=10; state=state0 end

# I = @acset SIR2 begin Agents=1; State=3; state=2 end
# I2 = @acset SIR2 begin Agents=2; State=3; state=[2,2] end
# SI = @acset SIR2 begin Agents=2; State=3; state=[1,2] end

# # L_infect =  homomorphism(I, SI) # equivalent to below b/c only one morphism
# homomorphism(I, SI)
# homomorphism(I, I2)

# L_infect =  ACSetTransformation(I, SI; Agents=[1])
# R_infect = ACSetTransformation(I, I2)


text:
#=
using .Soss
using .Soss: ConditionalModel, getmoduletypencoding, argvals, observations,
             Assign, Return, LineNumber, Sample, testvalue, @q, type2model,
             loadvals, getntkeys, buildSource, isleaf, _unwrap_type
using NestedTuples
using MacroTools
ascube(m::ConditionalModel{A, B}, _data::NamedTuple) where {A,B} = ascube(m | _data)

function ascube(m::ConditionalModel{A, B}) where {A,B}
    return _ascube(getmoduletypencoding(m), Model(m), argvals(m), observations(m))
end

# function ascube(m::Model{EmptyNTtype, B}) where {B}
#     return ascube(m,NamedTuple())
# end

ascube(d, _data) = nothing

ascube(μ::MT.AbstractMeasure,  _data::NamedTuple=NamedTuple()) = ascube(μ)
ascube(μ::Dists.Distribution,  _data::NamedTuple=NamedTuple()) = ascube(μ)

ascube(d::Dists.AbstractMvNormal, _data::NamedTuple=NamedTuple()) = ascube(d)
#ascube(::NamedTuple{(), Tuple{}}) = nothing

export sourceascube

sourceascube(m::Model) = sourceascube()(m)

function sourceascube(_data=NamedTuple())
    function(_m::Model)

        _datakeys = getntkeys(_data)
        proc(_m, st::Assign)        = :($(st.x) = $(st.rhs))
        proc(_m, st::Return)     = nothing
        proc(_m, st::LineNumber) = nothing

        function proc(_m, st::Sample)
            x = st.x
            xname = QuoteNode(x)
            rhs = st.rhs

            thecode = @q begin
                _t = ascube($rhs, get(_data, $xname, NamedTuple()))
                if !isnothing(_t)
                    _result = merge(_result, ($x=_t,))
                end
            end

            # Non-leaves might be referenced later, so we need to be sure they
            # have a value
            isleaf(_m, st.x) || pushfirst!(thecode.args, :($x = Soss.testvalue($rhs)))

            return thecode
        end

        wrap(kernel) = @q begin
            _result = NamedTuple()
            $kernel
            $ascube(_result)
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten

    end
end

@gg function _ascube(M::Type{<:GeneralizedGenerated.TypeLevel}, _m::Model{Asub,B}, _args::A, _data) where {Asub,A,B}
    body = type2model(_m) |> sourceascube(_data) |> loadvals(_args, _data)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end
=#


text:
using Test

@testset "characters construction" begin
    @test Characters{3, UInt16}("abα") == "abα"
    @test Characters{2}("abc") == "ab"
    @test length(Characters{12}("12     ")) == 2
    @test length(Characters{3}("helloworld")) == 3
    @test String(Characters{12}(" abc  ")) == " abc"
    @test isequal(Characters(""), missing)
    @test isequal(Characters{3, UInt16}.(["a", "b", "", missing]), ["a","b", missing, missing])
    @test isequal(Characters{3, UInt8}.(["a", "b", "", missing]), ["a","b", missing, missing])
end

@testset "characters comparison" begin
    @test c"12" == "12"
    @test c"12  " == "12"
    @test "ab cd e" == c"ab cd e"
    @test c"gh   " == c"gh"
    @test "12" < c"13"
    @test !isless(c"12  ", "12")
    @test !isless(c"12  ", c"12")
    @test isless("abc", c"xy    ")
    @test isless(c"abc", "x y z")
    @test isless(c"1", missing)
end


text:
#==============================================================================

  Lifcycle Model based on code from Ben Moll's website:
      www.princeton.edu/~moll/HACTproject.htm

==============================================================================#

using LinearAlgebra, SparseArrays, Plots

# First set up model parameters
γ = 2 #CRRA utility parameter
σ_sq = (0.8)^2 # for the processes for productivity
corr = exp(-.9) # persistence for productivity
ρ = 0.05 # discount rate
r = 0.035 #interest rate
w =1 #wages

zmean = exp(σ_sq/2) # mean productivity
θ = -log(corr) # parameter for OU process for productivity

J =15 # number of points in gridspace for productivity
zmin = 0.75
zmax = 2.5

H = 300 #number of points in gridspace for assets
amin = 0
amax=100

T = 75 # maximum age of agents
N = 300 # number of age steps
dt = T/N

maxit =1000 #max number of iterations for HJB
crit = 10^(-10) # convergence criterion for HJB
convergence_criterion=10^(-5) #criterion for the KF and HJB?

# Set up some key variables
a = LinRange(amin,amax, H)
da = (amax-amin)/(H-1)
aa = a*ones(1,J)

z = LinRange(zmin,zmax,J)
dz = (zmax-zmin)/(J-1)
dz_sq = dz^2
zz=z*ones(1,H)

μ = -θ.*z.*log.(z)+σ_sq/2*z # drift for z
Σ_sq = σ_sq.*z.^2 #variance from Ito's lemma

# Create the matrix that summarizes the evolution of our system in the z dimension
 yy = (-Σ_sq/dz_sq - μ/dz)
 χ = Σ_sq/(2*dz_sq)
 ζ = μ/dz + Σ_sq/(2*dz_sq)

 # Define the Diagonals of this matrix
 updiag_z = zeros(H,1)
 	for j = 1:J
		global updiag_z =[updiag_z; repeat([ζ[j]], H, 1)]
	end

 centerdiag_z=repeat([χ[1]+yy[1]],H,1)
	for j = 2:J-1
		global centerdiag_z = [centerdiag_z; repeat([yy[j]], H, 1)]
	end
 centerdiag_z=[centerdiag_z; repeat([yy[J]+ζ[J]], H, 1)]

lowdiag_z = repeat([χ[2]], H, 1)
	for j=3:J
		global lowdiag_z = [lowdiag_z; repeat([χ[j]],H,1)]
	end
# spdiags in Matlab allows for automatic trimming/adding of zeros
    # spdiagm does not do this

Aswitch =  spdiagm(0=>centerdiag_z[:], -H=>lowdiag_z[:], H=>updiag_z[H+1:end-H])

# Set terminal condition, i.e. value of death ≈ 0
small_number=10^(-8)
v_terminal = small_number*(small_number.+aa).^(1-γ)/(1-γ)

global V=v_terminal

# Set up spaces for variables from the finite difference loop
Vaf, Vab, c = [zeros(H,J) for i in 1:3]

v = zeros(H,J,N)
gg =[]
A_t=[]
c_t=[]
ss_t=[]

#================ Main loop =================================================#

for n in N:-1:1
	global V=V
    v[:,:,n]=V
    #forward difference
    Vaf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
    Vaf[H,:] = (w*z .+ r.*amax).^(-γ) #will never be used, but impose state constraint a<=amax just in case
    #backward difference
    Vab[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
    Vab[1,:] = (w*z .+ r.*amin).^(-γ) #state constraint boundary condition

    #I_concave = Vab .> Vaf #indicator whether value function is concave (problems arise if this is not the case)

    # consumption and savings with forward difference
    cf = Vaf.^(-1/γ)
    ssf = zz' + r.*aa - cf
    # consumption and savings with backward difference
    cb = Vab.^(-1/γ)
    ssb = zz' + r.*aa - cb
    # consumption and derivative of value function at steady state
    c0 = zz' + r.*aa
    V0 = c0.^(-γ)

    #= dV_upwind makes a choice of forward or backward differences based on
     the sign of the drift    =#
    If = ssf .> 0 #positive drift → forward difference
    Ib = ssb .< 0 #negative drift → backward difference
    I0 = (1.0 .-If-Ib) #at steady state

    dV_Upwind = Vaf.*If + Vab.*Ib + V0.*I0  #important to include third term
    c = dV_Upwind.^(-1/γ)
    u = c.^(1-γ)/(1-γ)

    # CONSTRUCT MATRIX
    X = - min.(ssb,0)/da
    Y = - max.(ssf,0)/da + min.(ssb,0)/da
    Z = max.(ssf,0)/da

	# Define the Diagonals of this transition matrix
	updiag = 0
       for j = 1:J
           updiag =[updiag; Z[1:H-1,j]; 0]
       end
    updiag =(updiag[:])

    centerdiag=reshape(Y, H*J, 1)
    centerdiag = (centerdiag[:]) # for tuples

   lowdiag = X[2:H, 1]
       for j = 2:J
           lowdiag = [lowdiag; 0; X[2:H,j]]
       end
   lowdiag=(lowdiag)

   # spdiags in Matlab allows for automatic trimming/adding of zeros
       # spdiagm does not do this
  	AA = sparse(Diagonal(centerdiag))+ [zeros(1, H*J); sparse(Diagonal(lowdiag)) zeros(H*J-1,1)] + sparse(Diagonal(updiag))[2:end, 1:(H*J)] # trim first element
	A = AA + Aswitch

    push!(A_t,A)
    B = (1/dt + ρ)*sparse(I, H*J, H*J) - A

    u_stacked = reshape(u, H*J,1)
    V_stacked = reshape(V, H*J,1)

    b = u_stacked + V_stacked/dt
    V_stacked = B\b #SOLVE SYSTEM OF EQUATIONS

    global V = reshape(V_stacked, H,J)

    push!(ss_t,zz'+ r.*aa -c)
	push!(c_t, c)
end

# Plots

#inital consumption
plot(a, c_t[end], legend=false)
#inital savings
plot(a, ss_t[end], legend=false)
plot!(a, zeros(1,H))

# Create plots from lifecycle.pdf
dt_inv = convert(Int64, 1/dt)
dt_inv40 = convert(Int64, 40/dt)
dt_inv70 = convert(Int64, 70/dt)

# Policy function for consumption

p1=plot(a, c_t[300-dt_inv][:,1], label="Age 1, Lowest Income", ylims=(0,5), xlims=(0,80))
plot!(a, c_t[300-dt_inv][:,J], label="Age 1, Highest Income", xlabel="Wealth")
plot!(a, c_t[300-dt_inv40][:,1], label="Age 40 Lowest Income", ylabel="Consumption" )
plot!(a, c_t[300-dt_inv40][:,J], label="Age 40 Highest Income", legend=:bottomright)
plot!(a, c_t[300-dt_inv70][:,1], label="Age 70 Lowest Income")
plot!(a, c_t[300-dt_inv70][:,J], label="Age 70 Highest Income")

# Policy function for savings
p2= plot(a, ss_t[300-dt_inv][:,1], label="Age 1, Lowest Income", ylims=(-5,2), xlims=(0,100))
plot!(a, ss_t[300-dt_inv][:,J], label="Age 1, Highest Income", xlabel="Wealth")
plot!(a, ss_t[300-dt_inv40][:,1], label="Age 40 Lowest Income", ylabel="Savings" )
plot!(a, ss_t[300-dt_inv40][:,J], label="Age 40 Highest Income", legend=:bottomright)
plot!(a, ss_t[300-dt_inv70][:,1], label="Age 70 Lowest Income" )
plot!(a, ss_t[300-dt_inv70][:,J], label="Age 70 Highest Income")
plot!(a, zeros(H,1), line=:dash, color=:black, label="")

plot(p1,p2,legend=:bottomright,layout = (1,2))


text:
export twri_objective, TWRIOptions

# TWRI options
mutable struct TWRIOptions
    grad_corr::Bool
    comp_alpha::Bool
    weight_fun
    eps
    params
    Invq::String
end

"""
    TWRIOptions
        grad_corr::Bool
        comp_alpha::Bool
        weight_fun
        eps
        params::Symbol
        Invq::String

Options structure for TWRI.

`grad_corr`: Whether to add the gradient correction J'(m0, q)*∇_y

`comp_alpha`: Whether to compute optimal alpha (alpha=1 if not)

`weight_fun`: Whether to apply focusing/weighting function to F(m0)'*y and its norm

`eps`: Epsilon (noise level) value (default=0)

`Invq`: How to compute F'Y, either as full field or as a rank 1 approximation `w(t)*q(x)` using the source wavelet for w

`param`: Which gradient to compute. Choices are `nothing` (objective only), `:m`, `:y` or `:all`

Constructor
==========

All arguments are optional keyword arguments with the following default values:

TWRIOptions(;grad_corr=false, comp_alpha=true, weight_fun=nothing, eps=0, params=:m)
"""

TWRIOptions(;grad_corr=false, comp_alpha=true,
            weight_fun=nothing, eps=0, params=:m, Invq="standard")=
            TWRIOptions(grad_corr, comp_alpha, weight_fun, eps, params, Invq)

function subsample(opt::TWRIOptions, srcnum::Int)
    eloc = length(opt.eps) == 1 ? opt.eps : opt.eps[srcnum]
    return TWRIOptions(opt.grad_corr, opt.comp_alpha, opt.weight_fun, eloc, opt.params, opt.Invq)
end

twri_objective(model_full::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing}) =
    twri_objective(model_full, source, dObs, y, Options(), TWRIOptions())
twri_objective(model_full::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing}, opt::Options) =
    twri_objective(model_full, source, dObs, y, opt, TWRIOptions())
twri_objective(model_full::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing}, opt::TWRIOptions) =
    twri_objective(model_full, source, dObs, y, Options(), opt)
twri_objective(model_full::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing}, twri_opt::TWRIOptions, opt::Options) =
    twri_objective(model_full, source, dObs, y, opt, twri_opt)

function twri_objective(model_full::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing},
                        options::Options, optionswri::TWRIOptions)
    # Load full geometry for out-of-core geometry containers
    dObs.geometry = Geometry(dObs.geometry)
    source.geometry = Geometry(source.geometry)

    # Limit model to area with sources/receivers
    if options.limit_m == true
        model = deepcopy(model_full)
        model = limit_model_to_receiver_area(source.geometry, dObs.geometry, model, options.buffer_size)
    else
        model = model_full
    end

    # Set up Python model structure
    modelPy = devito_model(model, options)
    dtComp = get_dt(model; dt=options.dt_comp)

    # Extrapolate input data to computational grid
    qIn = time_resample(source.data[1],source.geometry,dtComp)[1]
    obsd = typeof(dObs.data[1]) == SegyIO.SeisCon ? convert(Array{Float32,2}, dObs.data[1][1].data) : dObs.data[1]
    dObserved = time_resample(obsd, dObs.geometry, dtComp)[1]

    isnothing(y) ? Y = nothing : Y = time_resample(y.data[1], y.geometry, dtComp)[1]

    # Set up coordinates
    src_coords = setup_grid(source.geometry, model.n)  # shifts source coordinates by origin
    rec_coords = setup_grid(dObs.geometry, model.n)    # shifts rec coordinates by origin

    ~isempty(options.frequencies) ? freqs = options.frequencies : freqs = nothing
    ~isempty(options.frequencies) ? (wfilt, freqs) =  filter_w(qIn, dtComp, freqs) : wfilt = nothing
    obj, gradm, grady = pycall(ac."wri_func", PyObject,
                               modelPy, src_coords, qIn, rec_coords, dObserved, Y,
                               t_sub=options.subsampling_factor, space_order=options.space_order,
                               grad=optionswri.params, grad_corr=optionswri.grad_corr, eps=optionswri.eps,
                               alpha_op=optionswri.comp_alpha, w_fun=optionswri.weight_fun,
                               freq_list=freqs, wfilt=wfilt)

    if (optionswri.params in [:m, :all])
        gradm = remove_padding(gradm, modelPy.padsizes; true_adjoint=options.sum_padding)
        options.limit_m==true && (gradm = extend_gradient(model_full, model, gradm))
    end
    if ~isnothing(grady)
        grady = time_resample(grady, dtComp, dObs.geometry)
        grady = judiVector(dObs.geometry, grady)
    end

    return filter_out(obj, gradm, grady)
end

function filter_w(qIn, dt, freqs)
    ff = FFTW.fftfreq(length(qIn), 1/dt)
    df = ff[2] - ff[1]
    inds = [findmin(abs.(ff.-f))[2] for f in freqs]
    DFT = joDFT(length(qIn); DDT=Float32)
    R = joRestriction(size(DFT,1), inds; RDT=Complex{Float32}, DDT=Complex{Float32})
    qfilt = DFT'*R'*R*DFT*qIn
    return qfilt, freqs
end

filter_out(obj, ::Nothing, ::Nothing) = obj
filter_out(obj, m, ::Nothing) = obj, m
filter_out(obj, ::Nothing, y) = obj, y
filter_out(obj, m, y) = obj, m, y

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
mutable struct Matrix{T} <: ElementalMatrix{T}
    obj::Ptr{Cvoid}
end

for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:ComplexF32, :c),
                    (:ComplexF64, :z))
    @eval begin

        # destructor to be used in finalizer. Don't call explicitly
        function destroy(A::Matrix{$elty})
            ElError(ccall(($(string("ElMatrixDestroy_", ext)), libEl), Cuint,
                (Ptr{Cvoid},), A.obj))
            return nothing
        end

        function Matrix(::Type{$elty})
            obj = Ref{Ptr{Cvoid}}(0)
            ElError(ccall(($(string("ElMatrixCreate_", ext)), libEl), Cuint,
                (Ref{Ptr{Cvoid}},), obj))
            A = Matrix{$elty}(obj[])
            finalizer(destroy, A)
            return A
        end

        function getindex(A::Matrix{$elty}, i::Integer, j::Integer)
            x = Ref{$elty}(0)
            ElError(ccall(($(string("ElMatrixGet_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt, Ref{$elty}),
                A.obj, i - 1, j - 1, x))
            return x[]
        end

        function resize!(A::Matrix{$elty}, i::Integer, j::Integer = 1) # to mimic vector behavior
            ElError(ccall(($(string("ElMatrixResize_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt),
                A.obj, i, j))
            return A
        end

        function height(A::Matrix{$elty})
            rs = Ref{ElInt}(0)
            ElError(ccall(($(string("ElMatrixHeight_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
            A.obj, rs))
            return rs[]
        end

        function width(A::Matrix{$elty})
            rs = Ref{ElInt}(0)
            ElError(ccall(($(string("ElMatrixWidth_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
            A.obj, rs))
            return rs[]
        end

        function lockedBuffer(A::Matrix{$elty})
            rp = Ref{Ptr{$elty}}(0)
            ElError(ccall(($(string("ElMatrixLockedBuffer_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{Ptr{$elty}}),
                A.obj, rp))
            return rp[]
        end

        function setindex!(A::Matrix{$elty}, x::Number, i::Integer, j::Integer)
            ElError(ccall(($(string("ElMatrixSet_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt, $elty),
                A.obj, i - 1, j - 1, x))
            return A
        end
    end
end

Matrix() = Matrix(Float64)

# Julia convenience

pointer(A::Matrix) = lockedBuffer(A)

function similar(::Matrix, ::Type{T}, sz::Dims) where {T}
    A = Matrix(T)
    resize!(A, sz...)
    return A
end


text:
using GmshParser
using Test

@testset "GmshParser.jl" begin
    # Write your tests here.
end


text:
function costheta_uv(u,varr)
    #=Evaluates the cosine of theta between a fixe vector u and and array of vector v
    Returns the dot product of the u with each vector in v
    =#
    a = fill(zero(eltype(varr)),size(varr,1) )
    for k in 1:size(varr, 2)
        for j in 1:size(varr, 1)
            a[j] += u[k] * varr[j, k]
        end
    end
    return a
end

function norm(varr)
    a = fill(zero(eltype(varr)),size(varr,1) )
    for k in 1:size(varr, 2)
        for j in 1:size(varr, 1)
            a[j] += varr[j,k] * varr[j, k]
        end
    end
    return a
end

function assign_sky_zones(ra,dec,dRAzone,dDECzone)
    #assigns the zone for each object
    RA_ed=collect(0:dRAzone:360)
    DEC_ed=collect(-90:dDECzone:90)
    nRA=size(RA_ed,1)
    nDEC=size(DEC_ed,1)

    zone_num=zeros(Int32,size(ra,1))

    for ii in 1:nRA-1
        ibin=findall( (ra.> RA_ed[ii]) .& (ra.<=RA_ed[ii+1]))
        zone_num[ibin].=ii
    end

    for jj in 1:nDEC-1
        ibin=findall((dec.> DEC_ed[jj]) .& (dec.<=DEC_ed[jj+1]))
        zone_num[ibin].= zone_num[ibin] .+ (nRA*jj)
    end

    return zone_num
end

function assign_priority(nobj,priority,priority_frac)
    #= Generates uniform random numbers
    split from priority fraction
    assign priority
    =#

    @assert sum(priority_frac)==1

    obj_priority=zeros(Int32,nobj)
    rand_num=rand(nobj)
    #println(minimum(rand_num),' ',maximum(rand_num))
    n_class=size(priority,1)

    max_lim=0;
    #count_ass=zeros(Int32,nobj)
    for ii in 1:n_class
        min_lim=max_lim
        max_lim=min_lim+priority_frac[ii]
        ind_this=findall( (rand_num.>= min_lim) .& (rand_num.<=max_lim))
        obj_priority[ind_this] .= priority[ii]
        #count_ass[ind_this] .= count_ass[ind_this].+1
    end
    return obj_priority
end

function sort_by_zone_priority(indsel,zone_num,obj_priority)
    #= first sort by zone and then sort by priority for each zone
    This construct a new variable to sort finally
     first sort by main variable, find minimum difference
    take second variable and set in the minimum difference
    add the scale second variable to first one and sort by this new quantity
    =#
    #sort by zone
    izone_sort=sortperm(zone_num)
    mindif=zone_num[izone_sort[2]]-zone_num[izone_sort[1]]

    #scale the second variable
    pr_min=minimum(obj_priority)
    pr_max=maximum(obj_priority)

    #new quantity for sorting
    new_q=(obj_priority .- pr_min)./(pr_max-pr_min)
    new_q=new_q+zone_num
    izone_sort=sortperm(new_q,rev=:true)
    return izone_sort
end

function approx_convert_ra_dec_to_gal_b(ra,dec)
   #=Convert ra,dec in degree to galactic latitude only b
    This is only approximate and can have error of about 0.01 degree
    =#
    deg2rad=pi/180
    alpha_NGP=192.85*deg2rad #ra in deg
    delta_NGP=27.14*deg2rad  #dec in deg
    l_NCP=6.28308866829052 # in radian

    sinb=sin(delta_NGP)*sin(dec*deg2rad)+ cos(delta_NGP)*cos(dec*deg2rad)*cos(ra*deg2rad-alpha_NGP)
    #cosb=sqrt(1-sinb*sinb)
    #cos_sin_lb= cos(dec)*sin(ra*deg2rad-alpha_NGP)

    #println(sinb,' ',cosb,' ',cos_sin_lb,' ',cos_sin_lb/cosb)

    #get l longitude and b latitude
    b=asin(sinb)/deg2rad
    #l=l_NCP-asin(cos_sin_lb/cosb)
    return b #in degree
end

function split_Galactic_cap(ra,dec,rettype)
    #return the indices of objects in NGC or SGC
    b= @. approx_convert_ra_dec_to_gal_b(ra,dec)

    if(rettype=="boolean")
        ind_NGC=b.>0
        ind_SGC=b.<=0
    else
        ind_NGC=findall(b.>0)
        ind_SGC=findall(b.<=0)
    end
    #ind_SGC=deleteat!([1:size(ra);], ind_NGC)

    return ind_NGC,ind_SGC
end

function write_JLD2_file(outfile,zone_num,quant_dic)
    #=Write the structured data in jld2 format with subgroup feature
     For any tracer write following structure
      -> NGC/SGC -> groupd in ra-dec zones
    =#
    us_zone=sort(unique(zone_num))
    #append!(size(zone_num,1),fel_zone)
    #@show(fel_zone,size(fel_zone),size(us_zone))

    quant_write=["RA","DEC","PRIORITY","FITS_index_julia"]

    ind_NS=Dict()
    ind_NS["NGC"],ind_NS["SGC"]= split_Galactic_cap(quant_dic["RA"],quant_dic["DEC"],"boolean")

    groups_order=Dict()
    groups_order["NGC"]=Dict()
    groups_order["SGC"]=Dict()

    jldopen(outfile, "a+") do file
        for zz in 1:size(us_zone,1)
            zone=us_zone[zz]
            for sky in ["NGC","SGC"]
                isel=findall((zone_num .==zone) .&  (ind_NS[sky]))
                if(size(isel,1)==0)
                    continue
                end

                tgroup="$(sky)/zone_$(zone)"
                if(sky in keys(file))
                    if("zone_$(zone)" in keys(file[sky]))
                        gtmp=file[sky]["zone_$(zone)"]
                    else
                        gtmp=JLD2.Group(file, tgroup)
                    end
                else
                    gtmp=JLD2.Group(file, tgroup)
                end

                if(!("zone_order" in keys(groups_order[sky])) )
                        groups_order[sky]["zone_order"]=[]
                        groups_order[sky]["nobj"]=[]
                end

                append!(groups_order[sky]["zone_order"],["$(tgroup)/"])
                append!(groups_order[sky]["nobj"],size(isel,1))

                for quant in quant_write
                    gtmp["$(quant)"] = quant_dic[quant][isel]
                end #ends the quant loop
            end#ends the sky loop
        end #ends the zone loop

        #Now write the zone order
        groups_order["NS"]=Dict()
        groups_order["NS"]["zone_order"]=[]
        groups_order["NS"]["nobj"]=[]
        groups_order["NS"]["beg_index"]=[]

        count_NS=1
        for sky in ["NGC","SGC"]
            groups_order[sky]["beg_index"]=[]
            count=1
            #println("$(sky) $(size(groups_order[sky]["zone_order"])) $(size(groups_order[sky]["nobj"])))")
            for tt in 1:size(groups_order[sky]["zone_order"],1)
                append!(groups_order[sky]["beg_index"],count)
                count=count+groups_order[sky]["nobj"][tt]

                #Also append for count_NS
                append!(groups_order["NS"]["zone_order"],[groups_order[sky]["zone_order"][tt]])
                append!(groups_order["NS"]["nobj"],groups_order[sky]["nobj"][tt])
                append!(groups_order["NS"]["beg_index"],count_NS)
                count_NS=count_NS+groups_order[sky]["nobj"][tt]
            end

            file[sky]["order"]=groups_order[sky]
        end
        file["order"]=groups_order["NS"]

    end #closes the JLD2 file

end

function embed_sphere(ra, dec)
    #= Adapted from desimode.io """Embed `ra`, `dec` to a uniform sphere in three dimensions.
    Note that intentionallly x,y,z are kept at three arrays for efficiency as Julia uses column-major
    and three arrays will be easier to manipulate compared to one
    """
    =#

    deg2rad=pi/180.0

    phi = deg2rad.*ra
    theta = deg2rad.*(90.0 .- dec)
    r = @. sin(theta)
    nobj=size(ra,1)
    x_unit =@.  r * cos(phi)
    y_unit =@.  r * sin(phi)
    z_unit =@. cos(theta)

    return x_unit,y_unit,z_unit
end

function embed_sphere!(ra,dec,x_unit,y_unit,z_unit)
    #= Adapted from desimode.io """Embed `ra`, `dec` to a uniform sphere in three dimensions.
    Note that intentionallly x,y,z are kept at three arrays for efficiency as Julia uses column-major
    and three arrays will be easier to manipulate compared to one
    This version uses pre-allocated arrrays to store data and hence modifies its argument
    """
    =#

    deg2rad=pi/180.0

    phi = deg2rad.*ra
    theta = deg2rad.*(90.0 .- dec)
    r = @. sin(theta)
    nobj=size(ra,1)
    @. x_unit = r * cos(phi)
    @. y_unit = r * sin(phi)
    @. z_unit = cos(theta)
end

function embed_sphere_xyzmatrix(ra,dec)
    #= calles embed sphere but returns a matrix which is sometimes convenient=#
    nra=size(ra,1)
    txyz=zeros(nra,3)
    embed_sphere!(ra, dec,view(txyz,1:nra,1),view(txyz,1:nra,2),view(txyz,:,3))

    return txyz
end

function load_DESI_tiles(tile_file,program,pass,sky)
    #=load the tiles from a tile file and select for program and pass if given=#

    #open fits file
    fin=FITS(tile_file)
    #println("Loading tiles from: $(tile_file)")

    nrow=size(read(fin[2],"RA"),1)

    tiles_dic=Dict()
    quant_need=["TILEID","RA","DEC","PASS"]

    sel_dic=Dict()
    if(program!=nothing)
        sel_dic["PROGRAM"]=program
    end
    if(pass!=nothing)
        sel_dic["PASS"]=pass
    end

    if(program==nothing && pass==nothing)
        indsel=1:nrow
    else
        indbool=trues(nrow)
        for tkey in keys(sel_dic)
            tval=read(fin[2],tkey)
            indbool .= indbool .& (tval .==sel_dic[tkey])
        end
        indsel=findall(indbool)
    end

    for qq in quant_need
        tiles_dic[qq]=keepat!(read(fin[2],qq),indsel)
    end

    #Now apply sky selection if needed
    if(sky!=nothing)
        ind_NS=Dict()
        ind_NS["NGC"],ind_NS["SGC"]=split_Galactic_cap(tiles_dic["RA"],tiles_dic["DEC"],"index")
        for qq in quant_need
            tiles_dic[qq]=tiles_dic[qq][ind_NS[sky]]
        end
    end

    close(fin)

    return tiles_dic
end

function deg2rad(rdeg)
    return pi*rdeg/180
end

function Assign_tile_id(tiles_dic,tile_radius_indeg,ra_tracer, dec_tracer)
    #= assign the object to the tile id, assumes tiles are non-overlapping=
    1) get xyz_unit for tile centers and ra, dec_tracer
    2) define the scales , assign zones, select tiles
    3) scan the object determine dot product and assign tiles
    =#

    nobj=size(ra_tracer,1)
    ntile=size(tiles_dic["RA"],1)

    #min_RA=minimum(ra_tracer)-5.0*tile_radius_indeg
    #determine the radius in dot product unit
    radius_in_costheta=cos(deg2rad(tile_radius_indeg))
    #@show radius_in_costheta

    #determine the tiles centers in unit sphere
    tiles_xyz= embed_sphere_xyzmatrix(tiles_dic["RA"],tiles_dic["DEC"])

    #determine the tracers on unit sphere
    tracer_xyz=embed_sphere_xyzmatrix(ra_tracer,dec_tracer)

    count_assign=zeros(Int32,nobj)
    tileid_assign=zeros(Int32,nobj)

    #scan the two
    #rdist=zeros(Float64,nobj)
    for tt in 1:ntile
        rdist=costheta_uv(tiles_xyz[tt,:],tracer_xyz)
        indsel=findall(rdist .> radius_in_costheta)
        #@show(tiles_xyz[tt,:],size(indsel),nobj,minimum(rdist),maximum(rdist))
        if(size(indsel,1)==0)
            continue
        end
        tileid_assign[indsel] .= tiles_dic["TILEID"][tt]
        count_assign[indsel] .= count_assign[indsel] .+ 1
    end
    @assert maximum(count_assign)<=1
    return tileid_assign
end

function traverse_group(jldfile,group)
    #= Traverse a JLD2 group for all the file within by loading order
    =#
    #jldopen("mocks/JLD2_data/LRG.jld2", "r") do file
    #    tgroup=file["$(group)/order"]
    #end

    return load(jldfile,"$(group)/order")
end

function check_variable_exists(jldfile,var_path)
    #= checks if a given variable path exists in the file=#
    var_arr=split(var_path,"/")
    file=jldopen(jldfile,"r")
    res=true
    tfile=file
    for tcol in var_arr
        if(!(tcol in keys(tfile)))
            res=false
            break
        end
        tfile=tfile[tcol]
    end

    close(file)
    return res
end

function load_groupzone_column(jldfile,group_zone,column)
    #=loads a column for a particular group and zone
    column should be a single column name=#

    this_col=load(jldfile,"$(group_zone)/$(column)")

    return this_col
end

function Append_groupzone_xyz_unit(jldfile,group_zone)
    #=Add xyz in unit sphere for a particular group and zone
    First check for if the column name exists=#
    columns=["X_UNIT","Y_UNIT","Z_UNIT"]

    for col_name in columns
        if(check_variable_exists(jldfile,"$(group_zone)$(col_name)"))
            println("exists: $(group_zone)$(col_name)")
            return
        end
    end

    #loads the ra,dec
    ra=load_groupzone_column(jldfile,group_zone,"RA")
    dec=load_groupzone_column(jldfile,group_zone,"DEC")
    #ra,dec=load(jldfile,"$(group)/RA","$(group)/DEC")

    #embeds the sphere
    unit_xyz=embed_sphere_xyzmatrix(ra, dec)

    #Now write this to the file
    jldopen(jldfile, "a+") do file
        for cc in 1:size(columns,1)
            col_name=columns[cc]
            write(file,"$(group_zone)/$(col_name)",unit_xyz[:,cc])
        end
    end
end

function load_group_column(jldfile,group,columns)
    #=loads a column for a particular group scanning all the zones
    columns should be Array of available properties=#

    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    #initialize the ditionary
    col_dic=Dict()
    for col in columns
        col_dic[col]=[]
    end

    for iz in 1:size(jld_order["zone_order"],1)
        tg=jld_order["zone_order"][iz]

        for col in columns
            append!(col_dic[col],load_groupzone_column(jldfile,tg,col))
        end

    end
    return col_dic
end

function Indexing_column(jldfile,group,columns,mask)
    #=Generate index file for groups
    writes a ditionary for each minimum and maximum values=#

    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    #initialize the ditionary
    minmax_dic=Dict()
    for col in columns
        minmax_dic[col]=zeros(size(jld_order["zone_order"],1),2)
    end

    for iz in 1:size(jld_order["zone_order"],1)
        tg=jld_order["zone_order"][iz]

        for col in columns
            this_col=load_groupzone_column(jldfile,tg,col)
            if(mask!=nothing)
                ind_sub=falses(size(this_col,1))
                for mm in mask
                    ind_sub .= ind_sub .|| (this_col .== mm)
                end
                indsel=findall(@. !ind_sub)
                keepat!(this_col,indsel)
            end
            if(size(this_col,1)==0)#This will happen only when mask is available
                minmax_dic[col][iz,:]=[mask[1],mask[1]]
            else
                minmax_dic[col][iz,:]=[minimum(this_col),maximum(this_col)]
            end
        end
    end

    #Now check if indexing exists
    update_dictionary_jld2file(jldfile,group,minmax_dic,"indexing")

    return minmax_dic

end

function LISTING_column(jldfile,group,outname,columns_use)
    #=For each zone generate the unique list of values in a particular colums
    This should be only used for columns which has small number of unique values
    for example TILEID=#

    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    #initialize the ditionary
    unique_dic=Dict()
    unique_dic[outname]=Dict()
    #for col in columns
    #    minmax_dic[col]=zeros(size(jld_order["zone_order"],1),2)
    #end

    for iz in 1:size(jld_order["zone_order"],1)
        tg=jld_order["zone_order"][iz]

        comb_col=[]
        for col in columns_use
            this_col=load_groupzone_column(jldfile,tg,col)
            append!(comb_col,unique(this_col))
        end
        unique_dic[outname][iz]=unique(comb_col)
    end

    #Now check if indexing exists
    update_dictionary_jld2file(jldfile,group,unique_dic,"listing")

    return unique_dic

end

function update_dictionary_jld2file(jldfile,group,new_index,key)
    #= check if dictionary exists then load and merge with new dictionary and write the new one
    first delete the old one before writing=#

    dic_exists=check_variable_exists(jldfile,"$(group)/$(key)")
    old_dic=Dict()

    if( dic_exists)
        old_dic=load(jldfile,"$(group)/$(key)")
    end

    #Include old indexing if not updated in new_index
    for tcol in keys(old_dic)
        if(tcol in keys(new_index))
            continue
        end
        new_index[tcol]=old_dic[tcol]
    end

    jldopen(jldfile, "a+") do file
        #delete if needed or just write
        if(dic_exists)
            Base.delete!(file,"$(group)/$(key)")
        end
        write(file,"$(group)/$(key)",new_index)
    end
end

function Append_column(jldfile,group,columns)
    #=Add the given propety if doesn't already exists
    A reference group needs to be given which can empty string or NGC/SGC=#

    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    for iz in 1:size(jld_order["zone_order"],1)
        tg=jld_order["zone_order"][iz]

        for col in columns
            if(col=="XYZ_UNIT")
                Append_groupzone_xyz_unit(jldfile,tg)
            else
                println("Not Implemented Append for $(col)")
            end
        end
    end
end

function Append_TileID(jldfile,tile_file,group,program,npass,tile_radius_indeg,overwrite)
    #= Assign tile id for each pass to a given tracer
    =#

    if("SGC" in split(group,"/"))
        sky="SGC"
    elseif("NGC" in split(group,"/"))
        sky="NGC"
    else
        sky=nothing
    end

    #load tiles for all the passes
    tiles_dic=Dict()
    for pass in 0:npass-1
        tiles_dic[pass]=load_DESI_tiles(tile_file,program,pass,sky)
    end

    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    for iz in 1:size(jld_order["zone_order"],1)
        tg=jld_order["zone_order"][iz]

        ra_tracer=load_groupzone_column(jldfile,tg,"RA")
        dec_tracer=load_groupzone_column(jldfile,tg,"DEC")

        for pass in 0:npass-1
            tile_pass_exists=check_variable_exists(jldfile,"$(tg)TILEID_PASS$(pass)")

            if(tile_pass_exists & !overwrite)
                println("$(tg)/TILEID_PASS$(pass) exists with no-overwrite")
                continue
            end
            #@show pass,tg,tile_pass_exists & !overwrite, tile_pass_exists

            tileid_assign=Assign_tile_id(tiles_dic[pass],tile_radius_indeg,ra_tracer, dec_tracer)
            jldopen(jldfile, "a+") do file
                #delete if needed or just write
                if(tile_pass_exists)
                    Base.delete!(file,"$(tg)TILEID_PASS$(pass)")
                end
                write(file,"$(tg)TILEID_PASS$(pass)",tileid_assign)
            end
        end
    end
end

""" load all objects in a particular tile of given tracer
    also returns indices with reference to the group which can be used to load
    or write back additional properties later
    if the tile_pass>0 then also look for object assignemnt in lower pass tiles and return a collision
    array
"""
function Load_tracers_intile(config,tracer,tile_id,tile_pass,group,columns;numobs_need=1,num_fba=1)
    #= load all objects in a particular tile
    also returns indices with reference to the group which can be used to load
    or write back additional properties later
    =#

    jldfile=config["target"][tracer]["JLDfile" ]

    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    #load the listing to identify the zones which has tile
    listing=load(jldfile,"$(group)/listing")

    indsel=[]
    #Dictionary to store the tracer properties
    tracer_tile=Dict()
    tracer_tile["ref_group"]=group
    tracer_tile["TILEID"]=tile_id
    tracer_tile["PASS"]=tile_pass
    tracer_tile["index"]=[]
    for col in columns
        tracer_tile[col]=[]
    end

    #list of tile which will have over lap in this region
    collided_tiles=zeros(Int32,0)

    for iz in 1:size(jld_order["zone_order"],1)
        tg=jld_order["zone_order"][iz]

        #if tileid not in listing then continue
        if( !(tile_id in listing["TILEID"][iz]))
            continue
        end

        #Now load the tile_id is in this zone
        this_id_arr=load(jldfile,"$(tg)TILEID_PASS$(tile_pass)")
        #It is important to make sure the reference groups are used consistently for indices to work
        indsel_this=findall(this_id_arr .== tile_id)

	append!(tracer_tile["index"], indsel_this.+ (jld_order["beg_index"][iz]-1) )

        for col in columns
            tcol_val=keepat!(load(jldfile,"$(tg)$(col)"),indsel_this)
            append!(tracer_tile[col],tcol_val)
        end

        #find all the tiles colliding with lower passes passes starts with 0
        if(tile_pass>0)
            for pre_pass in 0:tile_pass-1
                pre_id_arr=keepat!(load(jldfile,"$(tg)TILEID_PASS$(pre_pass)"),indsel_this)
                utiles=unique(pre_id_arr)
                for ut in utiles
                    if((!(ut in collided_tiles)) & (ut>0))
                        append!(collided_tiles,ut)
                    end
                end
                #println(pre_pass,' ',utiles)
            end
            #println("colide: ",collided_tiles)
        end
    end #end of iz loop

    #Now load the information whether object is already assigned
    if(size(collided_tiles,1)>0)
        count_assigned=count_assignment_tracer(config,tracer,tracer_tile,collided_tiles,num_fba)
        #Now populate the collided boolean array:
        #one set to true will not be considered for assignements
        tracer_tile["collided"] = count_assigned .>= numobs_need
    else
        tracer_tile["collided"] = zeros(Bool,size(tracer_tile["index"],1),num_fba)
    end

    return tracer_tile

end

function debug_minmax(collided,index,tdic;tmsg="")
   if(sum(collided[:,1])==0)
      return
   end
   indices=index[collided[:,1]]
   index_minmax=[minimum(indices),maximum(indices)]
   ind12=indexin(index_minmax,index)
   for tkey in keys(tdic)
      if(tdic[tkey] isa Array)
	 if(size(tdic[tkey],1) == size(index,1))
             string(tmsg,tkey,tdic[tkey][ind12[1]],tdic[tkey][ind12[2]],",")
	 end
      end
   end
   println("Debug:",tmsg)
   return
end

"""Counts the number of time any object is assigned in the given list of tiles
Only counts for the objects passed in the dictionary
fba_jldfile : JLD2 file storing assignmnet for the individual tiles
tracer_tile: A dictionary loaded for a tile along with its index
tile_list: list of tileID to count the assignmnet in
num_fba: number of fba realizations was run
"""
function count_assignment_tracer(config,tracer,tracer_tile,tile_list,num_fba)
    #total number od objects
    nobj_tile=size(tracer_tile["index"],1)
    #index to sort the objects in tile by index field
    tile_isort=sortperm(tracer_tile["index"])

    #allocate memory
    count_assigned=zeros(nobj_tile,num_fba)

    #scan the colliding tile and count number of assignment
    for col_tile in tile_list
        #FBA jldfile for the tile
	fba_jldfile=TileFBA_FileName(config,col_tile)

	index_pre_tile=load(fba_jldfile,"$(tracer)_index")
	bool_ass=load(fba_jldfile,"$(tracer)_Bool_ass")
        #make sure the tile was assigned with consistent num_fba
	@assert num_fba==size(bool_ass,2) "inconsistent num_fba \nfor tile: $(col_tile) in $(fba_jldfile)"

        #find the sorting permutation
        isort_perm=sortperm(index_pre_tile)

        #rearrange with sorted index for this tile
        index_pre_tile=index_pre_tile[isort_perm]
        bool_ass=bool_ass[isort_perm,:]

        #append the counts #follow the logic of sorted indices
        ctile=1
        for ti in 1:size(index_pre_tile,1)
            while(index_pre_tile[ti]>tracer_tile["index"][tile_isort[ctile]])
                ctile +=1
                if(ctile>nobj_tile)
                    break
                end
            end

            if(ctile>nobj_tile)
                break
            end

            if(index_pre_tile[ti]==tracer_tile["index"][tile_isort[ctile]])
                count_assigned[tile_isort[ctile],:] .+= bool_ass[ti,:]
            end
        end
    end #end of the col_tile loop
    return count_assigned
end# end of function

function pre_process_FITS2JLD_tracer(fname,priority,priority_frac,outfile)
    #= Pre-process a given tracer in following step
    priority is array containing list of priority
    priority_frac is simply the fraction for each priority should sum to 1
      1) read the fits file
      2) select y5 (2^1) foot-print and n(z) (2^0)
      3)  create zone ins ra-dec 20binx20bin
      3) split in NGC/SGC and next steps for each zone
        a) In each zone sort them according to priority (if multiple priority exists)
        b) Create a group and write each zone in seperate group
    =#

    fin=FITS(fname)
    println("$(now(UTC)) PREP Loading file: $(fname)")

    #apply footprint and nz selection
    status=read(fin[2],"STATUS")
    indsel=findall((status .& 2^1) .* (status .& 2^0).>1)
    #Only for testing to reduce size
    #indsel=indsel[1:1000:size(indsel,1)]
    nobj=size(indsel,1)
    println("$(now(UTC)) PREP nobj selection (footprint and nz): $(nobj)")

    #read ra,dec of selection
    ra=keepat!(read(fin[2],"RA"),indsel)
    dec=keepat!(read(fin[2],"DEC"),indsel)
    #redshift=keepat!(read(fin[2],"Z"),indsel)

    println("$(now(UTC)) PREP Assigning zones and priority and sorting")
    #Assign zone number for each object
    zone_num=assign_sky_zones(ra,dec,20,20)
    #asign_priority
    obj_priority=assign_priority(nobj,priority,priority_frac)

    #sort by zone and then priority within zone in reverse order (higher priority first)
    indsel_pr_zone=indsel_pr_zone=sort_by_zone_priority(indsel,zone_num,obj_priority)
    #for each priority further sorting can be done here if needed

    #re-ordering everything
    zone_num=zone_num[indsel_pr_zone]

    quant_dic=Dict()
    quant_dic["RA"]=ra[indsel_pr_zone]
    quant_dic["DEC"]=dec[indsel_pr_zone]
    #quant_dic["Z"]=redshift[indsel_pr_zone]
    quant_dic["PRIORITY"]=obj_priority[indsel_pr_zone]
    quant_dic["FITS_index_julia"]=indsel[indsel_pr_zone]

    println("$(now(UTC)) PREP Writing to JLD2 file: $(outfile)")
    write_JLD2_file(outfile,zone_num,quant_dic)

    close(fin)

end

"""
find the true index same as inbuilt findall but does this inplace avoiding new memory allocation
returns the number of elements in the output index_arr to be used
"""
function find_index!(bool_array,index_arr)
    nel=0
    @inbounds for ii in 1:size(bool_array,1)
        if(bool_array[ii])
            nel +=1
            index_arr[nel]=ii
        end
    end
    @assert nel<=size(index_arr,1) "The input index_array $(size(index_arr)) is smaller than needed for output $(nel)"
    return nel
end

"""
find the true index same as inbuilt findall but does this inplace avoiding new memory allocation
returns the number of elements in the output index_arr to be used
It performs an AND operation between the two inputs
"""
function find_index_AND!(bool_array1,bool_array2,index_arr)
    nel=0
    @inbounds for ii in 1:size(bool_array1,1)
        if(bool_array1[ii] & bool_array2[ii])
            nel +=1
            index_arr[nel]=ii
        end
    end
    @assert nel<=size(index_arr,1) "The input index_array $(size(index_arr)) is smaller than needed for output $(nel)"
    return nel
end

"""
find the true index same as inbuilt findall but does this inplace avoiding new memory allocation
returns the number of elements in the output index_arr to be used
It performs OR operation between the two inputs
"""
function find_index_OR!(bool_array1,bool_array2,index_arr)
    nel=0
    @inbounds for ii in 1:size(bool_array1,1)
        if(bool_array1[ii] || bool_array2[ii])
            nel +=1
            index_arr[nel]=ii
        end
    end
    @assert nel<=size(index_arr,1) "The input index_array $(size(index_arr)) is smaller than needed for output $(nel)"
    return nel
end

using Dates

"""
Converts a dictionary of events to the msg string desribing the time use
msg: any string to be prepended
events_dic: dictionar with array of events and times
depends in the Dates Package
example:
edic=Dict()
edic["events"]=["full","prepare","Assing"]
edic["times"] = [now(UTC)]
sleep(0.5)
append!(edic["times"],[now(UTC)])
sleep(0.2)
append!(edic["times"],[now(UTC)])

convert_events_to_msg(edic,"hello ")
wil give output:
"hello  full:0.709 sec(prepare:0.504 sec,Assing:0.205 sec,)"
"""
function convert_events_to_msg(events_dic,msg)

    nevent=length(events_dic["events"])
    if(nevent==0)
        return ""
    end

    #tdiff=string(Dates.canonicalize(Dates.CompoundPeriod(events_dic["times"][end] - events_dic["times"][1])))
    tdiff=Dates.DateTime(events_dic["times"][end]) - Dates.DateTime(events_dic["times"][1])
    tdiff=string(tdiff/Millisecond(1000), " sec")
    msg=string(msg," ",events_dic["events"][1],":",tdiff)

    msg_local=""
    for ii in 2:nevent
        #tdiff=string(Dates.canonicalize(Dates.CompoundPeriod(events_dic["times"][ii] - events_dic["times"][ii-1])))
        tdiff=Dates.DateTime(events_dic["times"][ii]) - Dates.DateTime(events_dic["times"][ii-1])
        tdiff=string(tdiff/Millisecond(1000), " sec,")
        msg_local=string(msg_local,events_dic["events"][ii],":",tdiff)
    end

    if(msg_local=="")
        return msg
    else
        msg_local=string(msg,"(",msg_local,")")
        return msg_local
    end
end

"""applly partitions and return indices
nobj::Integer Numbe of objects needs to be divided
npartition::Integer number of partitions needed
mypart: The partition needs to be selected should be between 0 and npartition-1
"""
function partition_indices(nobj,npartition,mypart)
    @assert ((mypart<npartition) & (mypart>=0)) "Selected partion $(mypart) is not in the range 0-$(npartition-1)"
    #How many should I execute
    nobj_per_part= Int(floor(nobj/npartition))

    if(nobj_per_part==0)
        println("Too many partition asked for, setting one object per partition,>$(nobj) partiitons will be empty ")
        nobj_per_part=1
        if(mypart<nobj)
            beg_index=mypart+nobj_per_part
            end_index=minimum([beg_index,nobj])
        else
            beg_index=1
            end_index=0
        end
    else
        #The indices of tiles I will execute
        beg_index=(mypart*nobj_per_part)+1
        end_index=minimum([beg_index+nobj_per_part-1,nobj])
        if(mypart==npartition-1)
            end_index=nobj
        end
    end

    return beg_index,end_index
end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
# Displays

type Inline end
type Clipboard end

type Editor end
type Console end
type PlotPane end

render(::Clipboard, x) = stringmime(MIME"text/plain"(), x)

render(::Editor, x) =
  render(Inline(), Copyable(x))

render(::Console, x) =
  Atom.msg("result", render(Inline(), Copyable(x)))

render(::PlotPane, x) =
  Atom.msg("plot", render(Inline(), x))

# View data structures

type Model
  data
end

immutable Tree
  head
  children::Vector{Any}
end

immutable LazyTree
  head
  children::Function
end

immutable SubTree
  label
  child
end

type Copyable
  view
  text::String
  Copyable(view, text::AbstractString) = new(view, text)
end

Copyable(view, text) = Copyable(view, render(Clipboard(), text))
Copyable(view) = Copyable(view, view)

immutable Link
  file::String
  line::Int
  contents::Vector{Any}
  Link(file::AbstractString, line::Integer, contents...) =
    new(file, line, [contents...])
end

Link(file::AbstractString, contents...) = Link(file, 0, contents...)

link(a...) = Link(a...)

type Row
  xs::Vector{Any}
  Row(xs...) = new(collect(xs))
end

@render Inline l::Row begin
  span([render(Inline(), x) for x in l.xs])
end


text:
# This file was generated, do not modify it. # hide
using MLJ
using UrlDownload
using PyPlot
ioff() # hide
import DataFrames: DataFrame, describe, names, select!
using Statistics

LinearRegressor = @load LinearRegressor pkg=MLJLinearModels

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module VectorDataUtils

using LoopVectorization
using LinearAlgebra
using Distances

export l2inner, l2norm, sql2norm, l2dist, sql2dist
export l1norm, l1dist
export xAx, mahaldist, sqmahaldist
export running_mean, running_mean!, update_running_mean!
export running_meanvar, running_meanvar!, update_running_meanvar!
export running_meancov, running_meancov!, update_running_meancov!
export sql2, l2, l1
export generate_4_clusters, generate_circle_cluster, generate_4_circle_clusters
export kmeans, kmeans!, kmpp, kmrand, kmeans_init, closest_centre, closest_centre!, mindist2cntrs

include("Lp.jl")
include("quadratic_forms.jl")
include("running_stats.jl")
include("distances.jl")
include("toy_data.jl")
include("kmeans.jl")
include("helper_functions.jl")

end


text:
# this is not testing hijack(), but load(revise=true), but is in this
# Hijack package because it's there that Revise gets loaded
module HijackTestsFail # ReTest.load will search for HijackTests, so this fails
using ReTest

@testset "Hijack: load(revise=true)" begin
    @test true
end
end # HijackTests

nothing # including this file must not return the module


text:
function _transform(x::Array{nf_elem}, y::fmpz_mat)
  n = length(x)
  @assert n == rows(y)
  m = cols(y)
  z = Array{nf_elem}(m)
  for i in 1:m
    z[i] = x[1]^y[1, i]
    for j in 2:n
      z[i] = z[i]*x[j]^y[j, i]
    end
  end
  return z
end

function _make_row_primitive(x::fmpz_mat, j::Int)
  y = x[j, 1]
  for i in 1:cols(x)
    y = gcd(y, x[j, i])
  end
  if y > 1
    for i in 1:cols(x)
      x[j, i] = div(x[j, i], y)
    end
  end
end

function _make_row_primitive!(x::Array{fmpz, 1})
  y = x[1]
  for i in 2:length(x)
    y = gcd(y, x[i])
    if y == 1
      return x
    end
  end
  if y > 1
    for i in 1:cols(x)
      x[i] = div(x[i], y)
    end
  end
end

function _primitive_element(F::FqNmodFiniteField)
  #println("Computing primitive element of $F")
  #println("Have to factor $(order(F) - 1)")
  fac = factor(order(F) - 1)
  while true
    a = rand(F)
    if iszero(a)
      continue
    end
    is_primitive = true
    for (l, ll) in fac
      if isone(a^(div(order(F) - 1, l)))
        is_primitive = false
        break
      end
    end
    if is_primitive
      return a
    end
  end
end


text:
for f in (:sampling_grid, :interpolation_grid, :platform_grid)
    @eval begin
        $f(ss::ProductSamplingStyle, ap::ApproximationProblem; options...) =
            ProductGrid(elements($f, ss, ap; options...)...)
    end
end

sampling_grid(ss::SamplingStyle, ap::ApproximationProblem; options...) =
    sampling_grid(ss, platform(ap), parameter(ap), samplingparameter(ss, ap; options...); options...)

sampling_grid(ss::InterpolationStyle, ap::ApproximationProblem; options...) =
    interpolation_grid(ss, platform(ap), parameter(ap); dict = ap.dict, options...)

# - generic grid: we invoke platform_grid on the platform
sampling_grid(samplingstyle::GridStyle, ap::ApproximationProblem; options...) =
    haskey(options,:grid) ? options[:grid] : platform_grid(samplingstyle, platform(ap), parameter(ap); options...)
platform_grid(samplingstyle::GridStyle, ap::ApproximationProblem, args...; options...) =
    haskey(options,:grid) ? options[:grid] : platform_grid(samplingstyle, platform(ap), parameter(ap); options...)
sampling_grid(ss::DiscreteGramStyle, ap::ApproximationProblem; options...) =
    oversampling_grid(ss, platform(ap), parameter(ap), samplingparameter(ss,ap;options...); options...)

sampling_grid(ss::OversamplingStyle, ap::ApproximationProblem; options...) =
    oversampling_grid(ss, platform(ap) ,parameter(ap), samplingparameter(ss,ap;options...);options...)

sampling_grid(samplingstyle::GridStyle, platform::Platform, param, args...; options...) =
    haskey(options,:grid) ?  options[:grid] : platform_grid(samplingstyle, platform(ap), parameter(ap); options...)
sampling_grid(samplingstyle::OversamplingStyle, platform::Platform, param, L, args...; options...) =
    oversampling_grid(samplingstyle, platform, param, L; options...)


text:
# LSP speech synthesis digital filter
mutable struct LSPDF <: LinearPredictionVariantsSynthesisFilter
    delay::Vector{Float64}

    function LSPDF(order)
        new(SPTK.lspdf_delay(order))
    end
end

function filt!(f::LSPDF, x::Float64, coef::Vector{Float64})
    SPTK.lspdf(x, coef, delay(f))
end

function to_filtcoef(f::LSPDF, c::SpectralParamState{LineSpectralPair})
    b = copy(c)
    rawdata(loggain!(b))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CoolProp"
version = v"6.4.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://sourceforge.net/projects/coolprop/files/CoolProp/$version/source/CoolProp_sources.zip", "b10b2be2f88675b7e46cae653880be93558009c8970d23e50ea917ce095791f6"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir

sed -i 's/Windows/windows/' CoolProp.sources/dev/Tickets/60.cpp
sed -i 's/Windows/windows/' CoolProp.sources/src/CPfilepaths.cpp
sed -i 's/.*-m.*BITNESS.*//' CoolProp.sources/CMakeLists.txt

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DCOOLPROP_SHARED_LIBRARY=ON ../CoolProp*/
VERBOSE=ON cmake --build . --config Release --target CoolProp -- -j${nproc}
mkdir -p ${libdir}
cp -a *CoolProp* ${libdir}
install_license $WORKSPACE/srcdir/CoolProp*/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libCoolProp", :libcoolprop)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
@enum HealthStatus begin
    hsMinValue
    hsUninitialized  # /< Health status has not been set.
    hsUnavailable    # /< Orbit information was not available, PVT invalid.
    hsUnused        # /< Sat health is not used in computing this PVT.
    hsUnknown       # /< Health state is unknown.
    hsUnhealthy     # /< Sat is marked unhealthy, do not use PVT.
    hsDegraded      # /< Sat is in a degraded state, recommend do not use.
    hsHealthy       # /< Satellite is healthy, PVT valid.
    hsMaxValue
end

@enum Frames begin
    # add new frames BEFORE count, then add to Strings[] in ReferenceFrame.cpp
    # and make parallel to this enum.

    # Unknown MUST BE FIRST, and MUST = 0
    fmUnknown    # /< unknown frame
    fmWGS84      # /< WGS84, assumed to be the latest version
    fmWGS84G730  # /< WGS84, GPS week 730 version
    fmWGS84G873  # /< WGS84, GPS week 873 version
    fmWGS84G1150 # /< WGS84, GPS week 1150 version
    fmITRF       # /< ITRF, assumed to be the latest version
    fmPZ90       # /< PZ90 (GLONASS)
    fmPZ90KGS    # /< PZ90 the "original"
    fmCGCS2000   # /< CGCS200 (BDS)
    # count MUST BE LAST
    fmcount        # /< the number of frames - not a frame
end

struct Xvt
    x::ECEF
    v::Tuple{Float64,Float64,Float64}
    clkbias::Float64
    clkdrift::Float64
    relcorr::Float64
    frame::Frames
    health::HealthStatus

end

@with_kw struct GNSSEllipsoid
    # / Defined in TR8350.2, Appendix A.1
    # / @return semi-major axis of Earth in meters.
    # a= 6378137.0,
    # / Derived from TR8350.2, Appendix A.1
    # / @return semi-major axis of Earth in km.
    # a_km = a / 1000.0,
    #= /**
    * Derived from TR8350.2, Appendix A.1
    * @return flattening (ellipsoid parameter).
    */ =#
    # flattening = 0.335281066475e-2,
    # / Defined in TR8350.2, Table 3.3
    # / @return eccentricity (ellipsoid parameter).
    # eccentricity= 8.1819190842622e-2,
    # / Defined in TR8350.2, Table 3.3
    # / @return eccentricity squared (ellipsoid parameter).
    # eccSquared = 6.69437999014e-3,
    ell::Ellipsoid = wgs84_ellipsoid
    # / Defined in TR8350.2, 3.2.4 line 3-6, or Table 3.1
    # / @return angular velocity of Earth in radians/sec.
    angVelocity::Float64 = 7.292115e-5
    # / Defined in TR8350.2, Table 3.1
    # / @return geocentric gravitational constant in m**3 / s**2
    gm::Float64 = 3986004.418e8
    # / Defined in TR8350.2, 3.3.2 line 3-11
    # / @return Speed of light in m/s.
    c::Float64 =  299792458.0
end
const WGS84Ellipsoid = GNSSEllipsoid()
const GPSEllipsoid = GNSSEllipsoid(wgs84_ellipsoid,
                    7.2921151467e-5, 3.986005e14, C_MPS)
const PZ90Ellipsoid = (ge = GNSSEllipsoid(Ellipsoid(a="6378136.0",
                        f_inv="298.25784",),
                    7.292115e-5, 398600.4418e9, 299792458.0,),
                    j20 = -1.08262575e-3,
                    gm_km = 398600.4418,
                    a_km = 6378136.0 / 1000.0,)

include("RinexNavCalcOrb.jl")
include("RinexNavCalcFit.jl")

function calculate_pos(ees::EphEpochStore, dt::DateTime, )
    if ees.si.ss == sstGlonass
        return calculate_Fit_pos(ees, dt)
    elseif ees.si.ss == sstGeosync
        for (k, v) in ees.eph
            vd = v.d
            return ECEF(vd[4], vd[8], vd[12], )
        end
    else
        return calculate_Orbit_pos(ees, dt)
    end
    return ECEF{Float64}(NaN, NaN, NaN)
end
# function calculate_pos(nav::RinexNav, sat::SatID, jd::Float64)

# end


text:
struct QuantizationTable{T <: Union{UInt8, UInt16}} <: AbstractBlock{T}
    id::UInt8
    table::SVector{64,T}
    function QuantizationTable{T}(table::AbstractArray; id::Integer=0) where T
        id < 4 || error("`id` must be less than 4.")
        if size(table) == (8, 8)
            new{T}(id, [table[i] for i in ZIGZAG_INDICES])
        else
            new{T}(id, table)
        end
    end
end
QuantizationTable{T}(table::QuantizationTable{T}; id::Integer=table.id) where T =
    QuantizationTable{T}(table.table, id=id)

const QuantizationTable8 = QuantizationTable{UInt8}
const QuantizationTable16 = QuantizationTable{UInt16}

qt_precision(::QuantizationTable{T}) where T = 8sizeof(T)

Base.sizeof(::Type{QuantizationTable{T}}) where T = 1 + 64sizeof(T)
Base.sizeof(qt::QuantizationTable) = sizeof(typeof(qt))

Base.getindex(qt::QuantizationTable, i::Integer) = qt.table[i]

function Base.show(io::IO, qt::QuantizationTable{T}) where T
    id = Int(qt.id)
    bits = qt_precision(qt)
    print(io, "QuantizationTable", bits, "([")
    for y = 1:7
        print(io, join(string.(qt[y,:]), " "), "; ")
    end
    print(io, join(string.(qt[8,:]), " "), "], id=", qt.id, ")")
end

function Base.show(io::IO, ::MIME"text/plain", qt::QuantizationTable{T}) where T
    id = Int(qt.id)
    bits = qt_precision(qt)
    print(io, "8×8 ", bits, "-bit QuantizationTable{", T, "}(id=", id, "):")
    len = ndigits(maximum(qt)) + 2
    for y = 1:8
        println(io)
        for x = 1:8
            print(io, lpad(string(qt[y,x]), len))
        end
    end
end

function Base.read(io::IO, t::Type{<:QuantizationTable})
    pq, tq = u4u4(io)
    isconcretetype(t) && sizeof(eltype(t)) != pq - 1 && error("precision mismatch")
    if pq == 0x0
        return QuantizationTable8(read(io, 64), id=tq)
    elseif pq == 0x1
        return QuantizationTable16(hton.(read(io, UInt16, 64)), id=tq)
    else
        error()
    end
end

function Base.write(io::IO, qt::QuantizationTable)
    p = qt_precision(qt) == 8 ? 0x00 : 0x10
    write(io, qt.id | p)
    write(io, hton.(qt.table))
end

quantize(s::S, q::T) where {S, T} = q == zero(T) ? zero(S) : convert(S, round(s / q))
dequantize(s::S, q::T) where {S, T} = convert(S, s * q)

quantize(b::DCTBlock, qt::QuantizationTable) = DCTBlock(quantize.(b, qt))
dequantize(b::DCTBlock, qt::QuantizationTable) = DCTBlock(dequantize.(b, qt))