text:
###################################################################################################################
# GENOME PROPERTIES
###################################################################################################################
"""

    `get_genome_info(FASTA)`

    Get chromosome information.

    # Examples
    ```julia-repl
    julia> CpelNano.get_genome_info(fasta)
    ```
"""
function get_genome_info(fasta::String)::Tuple{Vector{String},Vector{Int64}}

    # Get info from FASTA file
    fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
    chr_sizes = fa_reader.index.lengths
    chr_names = fa_reader.index.names
    close(fa_reader)

    # Return info
    return chr_names,chr_sizes

end
"""

    `get_chr_fa_rec(CHR,FASTA)`

    Get CHR chromosome FASTA record.

    # Examples
    ```julia-repl
    julia> CpelNano.get_chr_fa_rec(chr,fasta)
    ```
"""
function get_chr_fa_rec(chr::String,fasta::String)::FASTA.Record

    # Get info from FASTA file
    fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
    fa_rec = fa_reader[chr]
    close(fa_reader)

    # Return info
    return fa_rec

end
"""
    `get_grps_from_cgs(CPG_POS,MIN_GRP_DIST)`

    Function that returns CpG groups.

    # Examples
    ```julia-repl
    julia> CpelNano.get_grps_from_cgs([10,25,28,50],10)
    3-element Array{CpelNano.CpgGrp,1}:
     CpelNano.CpgGrp(5:15, 10, 1:1, true)
     CpelNano.CpgGrp(20:33, 26, 2:3, false)
     CpelNano.CpgGrp(45:55, 50, 4:4, true)
    ```
"""
function get_grps_from_cgs(cpg_pos::Vector{Int64},min_grp_dist::Int64)::Vector{CpgGrp}

    # Init tmp vars
    cpg_ind = 1:1
    grp_int = cpg_pos[1]:cpg_pos[1]
    length(cpg_pos)==1 && return [CpgGrp(grp_int,cpg_ind)]

    # Find groups
    i = 2
    cpg_grps = Vector{CpgGrp}()
    while i<=length(cpg_pos)

        # Check if grp is complete
        if cpg_pos[i]-maximum(grp_int)>min_grp_dist
            # Store CpG group
            push!(cpg_grps,CpgGrp(grp_int,cpg_ind))
            # Create new tmp vars
            cpg_ind = i:i
            grp_int = cpg_pos[i]:cpg_pos[i]
        else
            # Update tmp vars
            cpg_ind = minimum(cpg_ind):i
            grp_int = minimum(grp_int):cpg_pos[i]
        end

        # If last one, store
        i==length(cpg_pos) && push!(cpg_grps,CpgGrp(grp_int,cpg_ind))

        # Increase counter
        i += 1

    end

    # Adjust intervals for tails
    @inbounds for g=1:length(cpg_grps)

        # Get coordinates of first and last C in CG context
        grp_st = minimum(cpg_grps[g].grp_int)
        grp_end = maximum(cpg_grps[g].grp_int)

        # Update interval field correcting for tails
        cpg_grps[g].grp_int = (grp_st-5):(grp_end+5)

    end

    # Return groups
    return cpg_grps

end
"""
    `get_ρn(DNA_SEQ,CPG_POS_OFFSET)`

    Function that returns CpG density vector ρn.

    # Examples
    ```julia-repl
    julia> CpelNano.get_ρn(dna_seq,cpg_pos_offset)
    ```
"""
function get_ρn(dna_seq::String,cpg_pos_offset::Vector{Int64};width=1000)::Vector{Float64}

    # Loop over each CpG site
    ρn = Vector{Float64}()
    @inbounds for n=1:length(cpg_pos_offset)

        # Get window st and end
        win_st = max(1,cpg_pos_offset[n]-width÷2)
        win_end = min(length(dna_seq),cpg_pos_offset[n]+width÷2)

        # Get num CpG sites
        push!(ρn,length(collect(eachmatch(r"[Cc][Gg]",dna_seq[win_st:win_end]))))

    end

    # Return density
    return ρn./width

end
"""
    `get_ρl(ρn,CPG_GRPS)`

    Function that returns average CpG density per group vector ρl.

    # Examples
    ```julia-repl
    julia> CpelNano.get_ρl(ρl,cpg_grps)
    ```
"""
function get_ρl(ρn::Vector{Float64},cpg_grps::Vector{CpgGrp})::Vector{Float64}

    # Loop over each CpG group
    ρl = Vector{Float64}()
    @inbounds for l=1:length(cpg_grps)

        # Get average
        aux = sum(ρn[cpg_grps[l].cpg_ind])
        aux /= length(cpg_grps[l].cpg_ind)

        # Push
        push!(ρl,aux)

    end

    # Return average group density
    return ρl

end
"""
    `get_dn(CPG_POS)`

    Function that returns distance between contiguous CpG sites.

    # Examples
    ```julia-repl
    julia> CpelNano.get_dn([10,20,30])
    ```
"""
function get_dn(cpg_pos::Vector{Int64})::Vector{Int64}

    # Loop
    dn = Vector{Int64}()
    @inbounds for n=2:length(cpg_pos)
        push!(dn,cpg_pos[n]-cpg_pos[n-1])
    end

    # Return dn
    return dn

end
"""
    `get_dl(CPG_POS,CPG_GRPS)`

    Function that returns distance between contiguous CpG groups.

    # Examples
    ```julia-repl
    julia> CpelNano.get_dl(cpg_pos,cpg_grps)
    ```
"""
function get_dl(cpg_pos::Vector{Int64},cpg_grps::Vector{CpgGrp})::Vector{Int64}

    # Loop
    dl = Vector{Int64}()
    @inbounds for l=2:length(cpg_grps)
        cpg1 = cpg_pos[maximum(cpg_grps[l-1].cpg_ind)]
        cpg2 = cpg_pos[minimum(cpg_grps[l].cpg_ind)]
        push!(dl,cpg2-cpg1)
    end

    # Return dl
    return dl

end
"""
    `get_grp_info!(REG,FA_REC,MIN_GRP_DIST)`

    Function that stores information about each CpG group, separated by a minimum of min_grp_dist, in REG.

    # Examples
    ```julia-repl
    julia> CpelNano.get_grp_info!(reg,fa_rec,min_grp_dist)
    ```
"""
function get_grp_info!(rs::RegStruct,fa_rec::FASTA.Record,min_grp_dist::Int64)::Nothing

    # Get size of chromosome
    chr_size = length(FASTA.sequence(fa_rec))

    # Get position CpG sites in analysis region
    dna_seq = FASTA.sequence(String,fa_rec,rs.chrst:rs.chrend)
    cpg_pos = map(x->getfield(x,:offset),eachmatch(r"[Cc][Gg]",dna_seq)) .+ rs.chrst .- 1
        # print_log("$(cpg_pos)")

    # Set individual CpG site information
    rs.cpg_pos = cpg_pos
    rs.N = length(cpg_pos)
    rs.N>1 || return nothing

    # Set group information
    cpg_grps = get_grps_from_cgs(cpg_pos,min_grp_dist)
    rs.Nl = map(grp->length(grp.cpg_ind),cpg_grps)
    rs.L = length(cpg_grps)
    rs.cpg_grps = cpg_grps
        # print_log("$(cpg_grps)")

    ## Density
    dna_cont_st = max(1,rs.chrst-500)
    dna_cont_end = min(chr_size,rs.chrend+500)
    dna_seq = FASTA.sequence(String,fa_rec,dna_cont_st:dna_cont_end)
        # print_log("$(dna_seq)")

    # Store CpG density
    rs.ρn = get_ρn(dna_seq,cpg_pos.-dna_cont_st)
        # print_log("ρn=$(rs.ρn)")

    # Store average group average density
    rs.ρl = get_ρl(rs.ρn,cpg_grps)
        # print_log("ρl=$(rs.ρl)")

    ## Distances

    # Store CpG distances
    rs.dn = get_dn(cpg_pos)
        # print_log("dn=$(rs.dn)")

    # Store group distances
    rs.dl = get_dl(cpg_pos,cpg_grps)
        # print_log("dl=$(rs.dl)")

        # readline()

    # Return nothing
    return nothing

end
"""
    `get_informME_chr_part(CHR_SIZE,CONFIG)`

    Function that partitions chromosome the same way informME does.

    # Examples
    ```julia-repl
    julia> CpelNano.get_informME_chr_part(chr_size,config)
    ```
"""
function get_informME_chr_part(chr_size::Int64,config::CpelNanoConfig)::Vector{UnitRange{Int64}}

    # Loop until entire chromosome is covered
    ok = false
    ar_st = 1
    ar_end = ar_st + config.size_est_reg - 1
    chr_part = Vector{UnitRange{Int64}}()
    while !ok

        # Check end
        ar_end = ar_end>chr_size ? chr_size : ar_end
        ok = ar_end==chr_size ? true : false

        # Push region
        push!(chr_part,ar_st:ar_end)

        # Update window
        ar_st += config.size_est_reg
        ar_end = ar_st + config.size_est_reg - 1

    end

    # Return partition
    return chr_part

end
"""
    `get_chr_part(CHR,CONFIG,FASTA)`

    Function that partitions chromosome into analysis regions. No CpG group is allowed to be split in the boundary.
    There are two methods for this function: (i) entire chromosome is partitioned, and (ii) only targeted regions
    as defined by the `targ_regs` arguments are partitioned.

    # Examples
    ```julia-repl
    julia> CpelNano.get_chr_part(chr,config,fasta)
    julia> CpelNano.get_chr_part(chr,config,fasta,targ_regs)
    ```
"""
function get_chr_part(chr::String,config::CpelNanoConfig,fasta::String)::Vector{UnitRange{Int64}}

    # Get chromosome sequence
    chr_seq = get_chr_fa_rec(chr,fasta)
    chr_size = length(FASTA.sequence(chr_seq))
    # config.verbose && print_log("chr_seq: $(chr_seq)")

    # Loop over chromosome
    i = 1
    chr_part = Vector{UnitRange{Int64}}()
    while i<chr_size

        # Candidate boundary coordinates
        bndy = i + config.size_est_reg - 1
        # config.verbose && print_log("bndy: $(bndy)")

        # Leave if boundary is beyond size
        if bndy>=chr_size
            push!(chr_part,i:chr_size)
            break
        end

        # Gather DNA sequence from ±1000
        dna_seq = FASTA.sequence(String,chr_seq,(bndy-1000):min(chr_size,bndy+1000))
        # config.verbose && print_log("dna_seq: $(dna_seq[1000:1002])")

        # Obtain absolute position of CpG sites
        cpg_pos = eachmatch(r"[Cc][Gg]",dna_seq)
        cpg_pos = map(x->getfield(x,:offset),cpg_pos) .+ bndy .- 1001
        # config.verbose && print_log("cpg_pos: $(cpg_pos)")

        # Check boundary if there are CpG sites
        if length(cpg_pos)>0

            # Obtain groups in window
            cpg_grps = get_grps_from_cgs(cpg_pos,config.min_grp_dist)
            # config.verbose && print_log("cpg_grps: $(cpg_grps)")

            # Identify problematic group (if any)
            ind_grp = findfirst(grp->bndy in grp.grp_int,cpg_grps)
            # config.verbose && print_log("ind_grp: $(ind_grp)")

            # Define bndy depending on case
            if !isnothing(ind_grp)
                # Set boundary depending on extend included
                grp_st = minimum(cpg_grps[ind_grp].grp_int)
                grp_end = maximum(cpg_grps[ind_grp].grp_int)
                bndy = 2*bndy-grp_st>grp_end ? grp_end+2 : grp_st-1
            end
            # config.verbose && print_log("Modified bndy: $(bndy)")
            # config.verbose && readline()

        end

        # Store and move on to next boundary
        push!(chr_part,i:bndy)
        i = bndy+1

    end

    # Return partition
    return chr_part

end
function get_chr_part(chr::String,config::CpelNanoConfig,fasta::String,targ_regs::Vector{UnitRange{Int64}})::Vector{UnitRange{Int64}}

    # Get chromosome sequence
    chr_seq = get_chr_fa_rec(chr,fasta)
    chr_size = length(FASTA.sequence(chr_seq))

    # Loop over chromosome
    chr_part = Vector{UnitRange{Int64}}()
    for reg in targ_regs

        # Get start and end of target region
        i = minimum(reg)
        reg_end = maximum(reg)

        # Partition target region
        while i<reg_end

            # Candidate boundary coordinates
            bndy = i + config.size_est_reg - 1

            # Leave if boundary is beyond size
            if bndy >= reg_end
                push!(chr_part,i:reg_end)
                break
            end

            # Gather DNA sequence from ±1000
            dna_seq = FASTA.sequence(String,chr_seq,(bndy-1000):min(chr_size,bndy+1000))

            # Obtain absolute position of CpG sites
            cpg_pos = eachmatch(r"[Cc][Gg]",dna_seq)
            cpg_pos = map(x->getfield(x,:offset),cpg_pos) .+ bndy .- 1001
            # config.verbose && print_log("cpg_pos: $(cpg_pos)")

            # Check boundary if there are CpG sites
            if length(cpg_pos)>0

                # Obtain groups in window
                cpg_grps = get_grps_from_cgs(cpg_pos,config.min_grp_dist)
                # config.verbose && print_log("cpg_grps: $(cpg_grps)")

                # Identify problematic group (if any)
                ind_grp = findfirst(grp->bndy in grp.grp_int,cpg_grps)
                # config.verbose && print_log("ind_grp: $(ind_grp)")

                # Define bndy depending on case
                if !isnothing(ind_grp)
                    # Set boundary depending on extend included
                    grp_st = minimum(cpg_grps[ind_grp].grp_int)
                    grp_end = maximum(cpg_grps[ind_grp].grp_int)
                    bndy = 2*bndy-grp_st>grp_end ? grp_end+2 : grp_st-1
                end
                # config.verbose && print_log("Modified bndy: $(bndy)")
                # config.verbose && readline()

            end

            # Store and move on to next boundary
            push!(chr_part,i:bndy)
            i = bndy+1

        end
    end

    # Return partition
    return chr_part

end
###################################################################################################################
# ANALYSIS REGION PROPERTIES
###################################################################################################################
"""
    `get_nls_reg_info!(REG_STRUCT,CONFIG)`

    Function that stores subregion information.

    # Examples
    ```julia-repl
    julia> CpelNano.get_nls_reg_info!(rs,config)
    ```
"""
function get_nls_reg_info!(rs::RegStruct,config::CpelNanoConfig)::Nothing

    # Get size of estimation region
    reg_size = (rs.chrend-rs.chrst)+1

    # Find number of necessary subregions
    k = Int(ceil(reg_size/config.max_size_subreg))

    # Base size of subregions
    base_size = Int(floor(reg_size / k))

    # Set starts
    bnds = fill(base_size,k)

    # Reminder
    rem_size = mod(reg_size,k)

    # Distribute reminder in order
    i = 1
    while rem_size > 0
        bnds[i] += 1
        rem_size -=1
        i = i==k  ? 1 : i+1
    end

    # Add up sizes to obtain offset from start
    bnds = [sum(bnds[1:i]) for i=1:k]

    # Offset all by chromosome start
    bnds .+= rs.chrst .- 1

    # Add chromosome start
    pushfirst!(bnds,rs.chrst)

    # Construct analysis region intervals
    sub_int = [bnds[i]:bnds[i+1] for i=1:k]

    # Divide CpG sites indices
    cpg_inds = fill(0:0,k)
    for i=1:k

        # Find first and last
        fst_ind = findfirst(x->bnds[i]<=x<bnds[i+1],rs.cpg_pos)
        fst_ind = isnothing(fst_ind) ? 0 : fst_ind
        lst_ind = findlast(x->bnds[i]<=x<bnds[i+1],rs.cpg_pos)
        lst_ind = isnothing(lst_ind) ? 0 : lst_ind

        # Store info
        if i<k
            # Excluding the one on right boundary
            cpg_inds[i] = fst_ind:lst_ind
        else
            # Include the CpG on right boundary if exists
            cpg_inds[i] = fst_ind:rs.N
        end

    end
        # print_log("$(rs.cpg_pos)")

    ## Store in struct
    rs.nls_rgs = AnalysisRegions(sub_int,cpg_inds)
        # print_log("$(rs.nls_rgs)")

    # Return nothing
    return nothing

end
"""
    `pmap_subregion_table(REG_INT,FA_REC)`

    Function that returns subregions size and number of CpG sites.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_subregion_table(reg_int,fa_rec)
    ```
"""
function pmap_subregion_table(reg_int::UnitRange{Int64},fa_rec::FASTA.Record,max_size_subreg::Int64)::Array{Int64,2}

    # Get start and end
    chrst = minimum(reg_int)
    chrend = maximum(reg_int)
        # print_log("Range: $(reg_int) ...")

    # Get size of estimation region
    reg_size = chrend-chrst+1

    # Find number of necessary subregions
    k = Int(ceil(reg_size/max_size_subreg))

    # Base size of subregions
    base_size = Int(floor(reg_size / k))

    # Set starts
    bnds = fill(base_size,k)

    # Reminder
    rem_size = mod(reg_size,k)

    # Distribute reminder in order
    i = 1
    while rem_size > 0
        bnds[i] += 1
        rem_size -=1
        i = i==k  ? 1 : i+1
    end

    # Add up sizes to obtain offset from start
    bnds = [sum(bnds[1:i]) for i=1:k]

    # Offset all by chromosome start
    bnds .+= chrst .- 1

    # Add chromosome start
    pushfirst!(bnds,chrst)
        # print_log("Bounds: $(bnds)")

    # Construct subregion intervals
    sub_int = [bnds[i]:bnds[i+1] for i=1:k]
    sub_len = length.(sub_int) .- 1
        # print_log("Subregion intervals: $(sub_int)")
        # print_log("Subregion lenghts: $(sub_len)")

    # Get position CpG sites
    dna_seq = FASTA.sequence(String,fa_rec,reg_int)
    cpg_pos = map(x->getfield(x,:offset),eachmatch(r"[Cc][Gg]",dna_seq)) .+ minimum(reg_int) .- 1

    # If not enough CpG sites for CpelNano return empty
    length(cpg_pos)>9 || return hcat([],[])

    # Divide CpG sites indices
    cpg_inds = fill(0:0,k)
    for i=1:k

        # Find first and last
        fst_ind = findfirst(x->bnds[i]<=x<bnds[i+1],cpg_pos)
        fst_ind = isnothing(fst_ind) ? 0 : fst_ind
        lst_ind = findlast(x->bnds[i]<=x<bnds[i+1],cpg_pos)
        lst_ind = isnothing(lst_ind) ? 0 : lst_ind

        # Store info
        if i<k
            # Excluding the one on right boundary
            cpg_inds[i] = fst_ind:lst_ind
        else
            # Include the CpG on right boundary if exists
            cpg_inds[i] = fst_ind:length(cpg_pos)
        end

    end
    n_cpgs = length.(cpg_inds)
        # print_log("CpG sites: $(cpg_inds)")
        # print_log("Num CpG sites: $(n_cpgs)")

    # Return matrix
    return hcat(sub_len,n_cpgs)

end
"""
    `subregion_table(FASTA,CONFIG)`

    Function that stores a tabulated table with columns:

        | Start of estimation region | End of estimation region | Subregion size | # CpG sites |

    # Examples
    ```julia-repl
    julia> CpelNano.subregion_table(fasta,config)
    ```
"""
function subregion_table(fasta::String,config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names,chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Loop over chromosomes
    for i=1:length(chr_names)

        # Get chromosome name and size
        chr = chr_names[i]
        chr_size = chr_sizes[i]

        # Print info
        print_log("Processing chr $(chr) of length $(chr_size) ...")

        # Get analysis region start positions
        print_log("Partitioning chr $(chr) ...")
        chr_part = get_chr_part(chr,config,fasta)

        # Get FASTA record
        fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
        fa_rec = fa_reader[chr]
        close(fa_reader)

        # Process each analysis region in chr
        print_log("Processing chr $(chr) ...")
        out_pmap = pmap(x->pmap_subregion_table(x,fa_rec,config.max_size_subreg),chr_part)

        # Produce a matrix
        out_pmap = vcat(out_pmap...)

        # Write histogram
        print_log("Writing results from chr $(chr) ...")
        open("$(config.out_dir)/subregion_table.txt","a") do io
            writedlm(io,out_pmap)
        end

    end

    # Return nothing
    return nothing

end
###################################################################################################################
# ESTIMATION REGION ANALYSIS
###################################################################################################################
"""
    `pmap_estimation_region_analysis(REG_INT,CHR,FA_REC,CONFIG)`

    Function that returns a vector with analysis start and end, number of CpG sites, and number of CpG groups.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_estimation_region_analysis(reg_int,chr,fa_rec,config)
    ```
"""
function pmap_estimation_region_analysis(reg_int::UnitRange{Int64},chr::String,fa_rec::FASTA.Record,min_grp_dist::Int64)::Vector{String}

    # Init
    out_vec = fill("",5)
    out_vec[1] = chr
    out_vec[2] = "$(minimum(reg_int))"
    out_vec[3] = "$(maximum(reg_int))"

    # Get position CpG sites
    dna_seq = FASTA.sequence(String,fa_rec,reg_int)
    cpg_pos = map(x->getfield(x,:offset),eachmatch(r"[Cc][Gg]",dna_seq)) .+ minimum(reg_int) .- 1

    # Store CpG sites
    out_vec[4] = "$(length(cpg_pos))"

    # Get groups
    out_vec[5] = length(cpg_pos)>0 ? "$(length(get_grps_from_cgs(cpg_pos,min_grp_dist)))" : "0"

    # Return vector
    return out_vec

end
"""
    `estimation_region_table(FASTA,CONFIG)`

    Function that stores a tabulated table with columns

        | Analysis region size | Number of CpG sites | Number of CpG groups |

    # Examples
    ```julia-repl
    julia> CpelNano.estimation_region_table(fasta,config)
    ```
"""
function estimation_region_table(fasta::String,config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names,chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Loop over chromosomes
    for i=1:length(chr_names)

        # Get chromosome name and size
        chr = chr_names[i]
        chr_size = chr_sizes[i]

        # Print info
        print_log("Processing chr $(chr) of length $(chr_size) ...")

        # Get analysis region start positions
        print_log("Partitioning chr $(chr) ...")
        chr_part = get_chr_part(chr,config,fasta)

        # Get FASTA record
        fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
        fa_rec = fa_reader[chr]
        close(fa_reader)

        # Process each analysis region in chr
        print_log("Analyzing chr $(chr) ...")
        out_pmap = pmap(x->pmap_estimation_region_analysis(x,chr,fa_rec,config.min_grp_dist),chr_part)

        # Produce a matrix
        out_pmap = permutedims(hcat(out_pmap...))

        # Write histogram
        print_log("Writing results from chr $(chr) ...")
        open("$(config.out_dir)/estimation_region_table.txt","a") do io
            writedlm(io,out_pmap)
        end

    end

    # Return nothing
    return nothing

end
"""
    `get_estimation_region_info(FASTA,CHR,CHR_INT,CONFIG)`

    Function that stores a tabulated table with columns

        | Position | Group index | Density | Distance |

    # Examples
    ```julia-repl
    julia> CpelNano.get_estimation_region_info(fasta,chr,chr_int,config)
    ```
"""
function get_estimation_region_info(fasta::String,chr::String,chr_int::UnitRange{Int64},config::CpelNanoConfig)::Nothing

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Get FASTA record
    fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
    fa_rec = fa_reader[chr]
    close(fa_reader)

    # Initizalize empty region structure
    rs = RegStruct()

    # Store chromosomal info (1-based)
    rs.chr = chr
    rs.chrst = minimum(chr_int)
    rs.chrend = maximum(chr_int)

    # Get CpG site information (1-based)
    get_grp_info!(rs, fa_rec, config.min_grp_dist)

    # Generate output matrix
    ρn = rs.ρn
    dn = vcat(rs.dn,0)
    cpg_pos = rs.cpg_pos
    out_mat = [cpg_pos ρn dn]

    # Write histogram
    out_file = "$(config.out_dir)/estimation_region_chr_$(chr)_$(minimum(chr_int))_$(maximum(chr_int)).txt"
    open(out_file,"a") do io
        writedlm(io,out_mat)
    end

    # Return nothing
    return nothing

end
###################################################################################################################
# GROUP LENGTH ANALYSIS
###################################################################################################################
"""
    `pmap_grp_size(CHR_ST,FASTA_REC,CONFIG)`

    Function that returns a vector with CG group sizes between coordinates REG_INT.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_grp_size(chrst,fa_rec,config)
    ```
"""
function pmap_grp_size(chrst::Int64,fa_rec::FASTA.Record,min_grp_dist::Int64)::Vector{Int64}

    # Get position CpG sites
    dna_seq = FASTA.sequence(String,fa_rec,chrst:(chrst+100000))
    cpg_pos = map(x->getfield(x,:offset),eachmatch(r"[Cc][Gg]",dna_seq)) .+ chrst .- 1

    # Set individual CpG site information
    length(cpg_pos)>0 || return Vector{Int64}()

    # Set group information
    cpg_grps = get_grps_from_cgs(cpg_pos,min_grp_dist)

    # Return vector of sizes
    return map(g->length(g.grp_int),cpg_grps)

end
"""
    `hist_group_length(FASTA,CONFIG)`

    Function that produces a plot with size of CpG groups in reference genome FASTA.

    # Examples
    ```julia-repl
    julia> CpelNano.hist_group_length(nano,fasta,config)
    ```
"""
function hist_group_length(fasta::String,config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names,chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Initialize histogram
    histo = fit(Histogram,[],0.5:1.0:50000.5)

    # Loop over chromosomes
    for i=1:length(chr_names)

        # Get chromosome name and size
        chr = chr_names[i]
        chr_size = chr_sizes[i]

        # Print info
        print_log("Processing chr $(chr) of length $(chr_size) ...")

        # Get analysis region start positions
        print_log("Partitioning chr $(chr) ...")
        chrst = 1:100000:(chr_size-100000)

        # Get FASTA record
        fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
        fa_rec = fa_reader[chr]
        close(fa_reader)

        # Process each analysis region in chr
        print_log("Measuring groups in chr $(chr) ...")
        out_pmap = pmap(x->pmap_grp_size(x,fa_rec,config.min_grp_dist),chrst)

        # Produce a single vector
        out_pmap = vcat(out_pmap...)

        # Get histogram
        print_log("Updating histogram with chr $(chr) ...")
        histo_chr = fit(Histogram,out_pmap,0.5:1.0:50000.5)
        histo = merge(histo,histo_chr)

    end

    # Write histogram
    open("$(config.out_dir)/len_group_histogram.txt","w") do io
        writedlm(io,histo.weights)
    end

    # Return nothing
    return nothing

end
###################################################################################################################
# GROUP SIZE (# CpGs) ANALYSIS
###################################################################################################################
"""
    `pmap_grp_cpg_num(CHR_ST,FASTA_REC,CONFIG)`

    Function that returns a vector with CG group sizes in terms of numbers of CpG sites between coordinates REG_INT.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_grp_cpg_num(chrst,fa_rec,config)
    ```
"""
function pmap_grp_cpg_num(chrst::Int64,fa_rec::FASTA.Record,min_grp_dist::Int64)::Vector{Int64}

    # Get position CpG sites
    dna_seq = FASTA.sequence(String,fa_rec,chrst:(chrst+100000))
    cpg_pos = map(x->getfield(x,:offset),eachmatch(r"[Cc][Gg]",dna_seq)) .+ chrst .- 1

    # Set individual CpG site information
    length(cpg_pos)>0 || return Vector{Int64}()

    # Set group information
    cpg_grps = get_grps_from_cgs(cpg_pos,min_grp_dist)

    # Return vector of sizes
    return map(g->length(g.cpg_ind),cpg_grps)

end
"""
    `hist_group_cpg_num(FASTA,CONFIG)`

    Function that produces a histogram of number of CpG sites per group in reference genome FASTA.

    # Examples
    ```julia-repl
    julia> CpelNano.hist_group_cpg_num(nano,fasta,config)
    ```
"""
function hist_group_cpg_num(fasta::String,config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names,chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Initialize histogram
    histo = fit(Histogram,[],0.5:1.0:10000.5)

    # Loop over chromosomes
    for i=1:length(chr_names)

        # Get chromosome name and size
        chr = chr_names[i]
        chr_size = chr_sizes[i]

        # Print info
        print_log("Processing chr $(chr) of length $(chr_size) ...")

        # Get analysis region start positions
        print_log("Partitioning chr $(chr) ...")
        chrst = 1:100000:(chr_size-100000)

        # Get FASTA record
        fa_reader = open(FASTA.Reader,fasta,index=fasta*".fai")
        fa_rec = fa_reader[chr]
        close(fa_reader)

        # Process each analysis region in chr
        print_log("Counting # CpG sites in groups in chr $(chr) ...")
        out_pmap = pmap(x->pmap_grp_cpg_num(x,fa_rec,config.min_grp_dist),chrst)

        # Produce a single vector
        out_pmap = vcat(out_pmap...)

        # Get histogram
        print_log("Updating histogram with chr $(chr) ...")
        histo_chr = fit(Histogram,out_pmap,0.5:1.0:10000.5)
        histo = merge(histo,histo_chr)

    end

    # Write histogram
    open("$(config.out_dir)/num_cpg_per_grp_histogram.txt","w") do io
        writedlm(io,histo.weights)
    end

    # Return nothing
    return nothing

end


text:
using Poirot, Test

@testset "Poirot" begin

@testset "Compiler" begin
  include("compiler.jl")
end

@testset "Inference" begin
  include("inference.jl")
end

end

ext: jl


text:
using Knet, Test
macro timeinclude(x)
    :(print($x); print("\t"); @time include($x))
end

@testset "Knet" begin

#                            commit ee1 c6a 0b5 e3a 9.2 8.3 8.3 6cb 6cb 8.6 6cb 6cb
#                           machine  ci rzr  sc  sc tig ai5 ai4 tr5 tr4 aws osx os4
#@timeinclude("kptr.jl")         #   ?   4   3   8  16   1   1   0   0  20   0   0
#@timeinclude("gpu.jl")          #   3   1   2   2   6   1   1   0   0   2   0   0
@timeinclude("distributions.jl") #   1   1   1   1   2   1   1   2   1   3   3   2
@timeinclude("dropout.jl")       #   6   1   2   8   5                   2
@timeinclude("serialize.jl")     #  16   2  11   1
@timeinclude("jld.jl")           #  11  11   9  26
@timeinclude("gcnode.jl")
@timeinclude("statistics.jl")    #
@timeinclude("cuarray.jl")       #  25
@timeinclude("bmm.jl")           #  10   3   9
@timeinclude("linalg.jl")        #  33  17  22  22  62  24  14  22   7  28  33  19
@timeinclude("batchnorm.jl")     #  32  14  23  22  93
@timeinclude("loss.jl")          #  32  13  19  20  10                   4
@timeinclude("karray.jl")        #  49  20  27  21  55  19  12   -   -  21   -   0
@timeinclude("reduction.jl")     #  67  37  52  55 106  40  21  29  11  57  55  29
@timeinclude("conv.jl")          #  71  33  48  51 107  22  12  62  47  26  44  16
@timeinclude("rnn.jl")           #  95  52  54  41  81                  12
@timeinclude("binary.jl")        # 124  64  75  80  56  34  19 491 119  51  53  25
@timeinclude("unary.jl")         # 172  86 103 103 122  42   6  36   4  56  67  11
@timeinclude("update.jl")        # 316  61  82  60  61  29  26 100  22  72  25  23
#TODO include("data.jl")
#TODO include("hyperopt.jl")
#TODO include("progress.jl")
#TODO include("random.jl")
#TODO include("train.jl")
#TODO include("uva.jl")

end


text:
using brainvisionloader
using Test

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


text:
using Plots
using Unitful
using LaTeXStrings
using DataFrames
using Unitful: N, mm, cm, m, inch, MPa

θ = LinRange(0, 2π, 100)
y = sin.(θ)

p = plot(θ, y,
	title = L"sin(\theta)",
	xlabel = L"\theta",
	ylabel = L"Amplitude",
	xlim = [0, 2π],
	ylim = [-1, 1],
	legend = false
)

display(p)

F = [1000, 3200, -4000, 4300, -2750, 6630]N
r = [300mm, 400mm, 10inch, 15cm, 0.188m, 27cm] .|> mm

df = DataFrame(r=r, F=F)


text:
###############################################################################
#
# Fetch data from the website
#
###############################################################################
"""
    fetch_RAW!(dt::MOD15A2Hv006LAI, year::Int)
    fetch_RAW!(data_url::String, data_loc::String, year::Int, label::String)

Download RAW product data, given
- `dt` [`AbstractUngriddedData`](@ref) type ungridded data type
- `year` Which year data to download
- `data_url` Url to the website where to download the data
- `data_loc` Local folder to store the data
"""
function fetch_RAW!(dt::MOD15A2Hv006LAI, year::Int)
    update_MODIS_password!();

    data_url = "$(MODIS_PORTAL)/MOLT/MOD15A2H.006/";
    data_loc = "$(MODIS_HOME)/MOD15A2H.006/original/";

    fetch_RAW!(data_url, data_loc, year, "MOD15A2H");

    return nothing
end

function fetch_RAW!(dt::MOD09A1v006NIRv, year::Int)
    update_MODIS_password!();

    data_url = "$(MODIS_PORTAL)/MOLT/MOD09A1.006/";
    data_loc = "$(MODIS_HOME)/MOD09A1.006/original/";

    fetch_RAW!(data_url, data_loc, year, "MOD09A1");

    return nothing
end

function fetch_RAW!(
            data_url::String,
            data_loc::String,
            year::Int,
            label::String
)
    # number of days per year
    TDAY = isleapyear(year) ? 366 : 365;

    # fetch file list in each folder
    @info "Fetching file name to download...";
    list_urls = [];
    list_locs = [];
    for doy in 1:8:TDAY
        #folder = parse_date(year, doy, ".") * "/";
        folder = "/";
        try
            @info "Fetching file list from $(data_url * folder)";
            download(data_url * folder, "temp.html");
            for _line in readlines("temp.html")
                if contains(_line, ".hdf\">")
                    _ini = findfirst(label, _line);
                    _end = findfirst(".hdf", _line);
                    _nam = _line[_ini[1]:_end[1]+3];
                    push!(list_urls, data_url * folder * _nam);
                    push!(list_locs, data_loc * string(year) * "/" * _nam);
                end
            end
            rm("temp.html");
        catch err
            @info "Unable to fetch files from $(folder), skip it...";
        end
    end

    # download files if the file does not exist
    @info "Downloading files...";
    @showprogress for i in eachindex(list_locs)
        _url = list_urls[i];
        _loc = list_locs[i];
        if Sys.which("wget") !== nothing
            if !isfile(_loc)
                _lst = `-q $(_url) -O $(_loc)`;
                _psd = `--user $(MODIS_USER_ID) --password $(MODIS_USER_PWD)`;
                run(`wget $(_psd) $(_lst)`);
            end
        else
            @warn "wget not found, exit the loop";
            break
        end
    end

    return nothing
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, Pkg

name = "muparser"
version = v"2.3.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/beltoforion/muparser/archive/refs/tags/v$(version).tar.gz", "b35fc84e3667d432e3414c8667d5764dfa450ed24a99eeef7ee3f6647d44f301")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/muparser-*

CMAKE_FLAGS=(-DCMAKE_INSTALL_PREFIX=$prefix
            -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
            -DCMAKE_BUILD_TYPE=Release
            -DENABLE_SAMPLES=OFF
            -DBUILD_TESTING=OFF
            -DBUILD_SHARED_LIBS=ON)

# Apple's Clang does not support OpenMP? - taken from AMRex build_tarballs.jl
if [[ ${target} == *-apple-* ]]; then

    CMAKE_FLAGS+=(-DENABLE_OPENMP=OFF)

else

    CMAKE_FLAGS+=(-DENABLE_OPENMP=ON)

fi

cmake . ${CMAKE_FLAGS[@]}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
#linux-musl variants hav an same issue as https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/387 on gcc 4, so use higher version
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function constraint(mechanism::Mechanism{T,Nn,Ne,Nb}, body::Body{T}) where {T,Nn,Ne,Nb}
    state = body.state
    timestep= mechanism.timestep

    mass = body.mass
    inertia = body.inertia
    gravity = mechanism.gravity

    x1, q1 = previous_configuration(state)
    x2, q2 = current_configuration(state)
    x3, q3 = next_configuration(state, timestep)

    # dynamics
    D1x = - 1.0 / timestep * mass * (x2 - x1) - 0.5 * timestep * mass * gravity
    D2x =   1.0 / timestep * mass * (x3 - x2) - 0.5 * timestep * mass * gravity
    D1q = -2.0 / timestep * LVᵀmat(q2)' * Lmat(q1) * Vᵀmat() * inertia * Vmat() * Lmat(q1)' * vector(q2)
    D2q = -2.0 / timestep * LVᵀmat(q2)' * Tmat() * Rmat(q3)' * Vᵀmat() * inertia * Vmat() * Lmat(q2)' * vector(q3)

    dynT = D2x + D1x
    dynR = D2q + D1q

    state.d = [dynT; dynR]

    # inputs
    state.d -= [state.JF2; state.Jτ2]

    # impulses
    for id in connections(mechanism.system, body.id)
        Ne < id <= Ne + Nb && continue # body
        impulses!(mechanism, body, get_node(mechanism, id))
    end

    return state.d
end

function constraint_jacobian_configuration(mechanism::Mechanism{T,Nn,Ne,Nb}, body::Body{T}; reg::T=Dojo.REG) where {T,Nn,Ne,Nb}
    state = body.state
    timestep = mechanism.timestep
    mass = body.mass
    inertia = body.inertia

    x2, q2 = current_configuration(state)
    x3, q3 = next_configuration(state, timestep)

    I3 = SMatrix{3,3,T,9}(Diagonal(sones(T,3)))
    Z33 = szeros(T, 3, 3)
    Z34 = szeros(T, 3, 4)

    # dynamics
    dynT = I3 * mass / timestep
    dynR = -2.0 / timestep * LVᵀmat(q2)' * Tmat() * (∂Rᵀmat∂q(Vᵀmat() * inertia * Vmat() * Lmat(q2)' * vector(q3)) + Rmat(q3)' * Vᵀmat() * inertia * Vmat() * Lmat(q2)')

    state.D = [[dynT; Z33] [Z34; dynR]] * integrator_jacobian_velocity(body, timestep)
    state.D += [[reg * I3; Z33] [Z33; reg * I3]]

    # inputs
    nothing

    # impulses
    for id in connections(mechanism.system, body.id)
        Ne < id <= Ne + Nb && continue # body
        impulses_jacobian_velocity!(mechanism, body, get_node(mechanism, id))
    end

    return state.D
end

function integrator_jacobian_velocity(body::Body{T}, timestep) where T
    state = body.state
    x2, v25, q2, ω25 = current_configuration_velocity(state)
    integrator_jacobian_velocity(x2, v25, q2, ω25, timestep)
end

function integrator_jacobian_configuration(body::Body{T},
        timestep; attjac::Bool=true) where T
    state = body.state
    x2, v25, q2, ω25 = current_configuration_velocity(state)
    integrator_jacobian_configuration(x2, v25, q2, ω25, timestep; attjac=attjac)
end

# linear system
function set_matrix_vector_entries!(mechanism, matrix_entry::Entry, vector_entry::Entry, body::Body)
    matrix_entry.value = constraint_jacobian_configuration(mechanism, body)
    vector_entry.value = -constraint(mechanism, body)
end


text:
module Metrics
using SparseArrays

struct CovariantBasis
        # a1
        x_r1::AbstractArray
        y_r1::AbstractArray
        # a2
        x_r2::AbstractArray
        y_r2::AbstractArray
end

struct ContravariantBasis
        # b1
        r1_x::AbstractArray
        r1_y::AbstractArray
        # b2
        r2_x::AbstractArray
        r2_y::AbstractArray
end

struct CovariantMetricTensor
        a11::AbstractArray
        a12::AbstractArray
        a21::AbstractArray
        a22::AbstractArray
end

struct ContravariantMetricTensor
        b11::AbstractArray
        b12::AbstractArray
        b21::AbstractArray
        b22::AbstractArray
end

function build_jacobian(a::CovariantBasis)
        return a.x_r1 .* a.y_r2 - a.x_r2 .* a.y_r1
end

function build_contravariant_basis(J::AbstractArray, a::CovariantBasis)
        return ContravariantBasis( a.y_r2 ./ J,
                                 - a.x_r2 ./ J,
                                 - a.y_r1 ./ J,
                                   a.x_r1 ./ J)

end

function build_covariant_metric_tensor(a::CovariantBasis)
        return CovariantMetricTensor(a.x_r1 .* a.x_r1 + a.y_r1 .* a.y_r1,
                                     a.x_r1 .* a.x_r2 + a.y_r1 .* a.y_r2,
                                     a.x_r2 .* a.x_r1 + a.y_r2 .* a.y_r1,
                                     a.x_r2 .* a.x_r2 + a.y_r2 .* a.y_r2)
end

function build_contravariant_metric_tensor(b::ContravariantBasis)
        return ContravariantMetricTensor(b.r1_x .* b.r1_x + b.r1_y .* b.r1_y,
                                         b.r1_x .* b.r2_x + b.r1_y .* b.r2_y,
                                         b.r2_x .* b.r1_x + b.r2_y .* b.r1_y,
                                         b.r2_x .* b.r2_x + b.r2_y .* b.r2_y)
end

function contravariant_metric_tensor_to_matrix(G::ContravariantMetricTensor)
        return ContravariantMetricTensor(spdiagm(0 => G.b11),
                                         spdiagm(0 => G.b12),
                                         spdiagm(0 => G.b21),
                                         spdiagm(0 => G.b22))
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ 61a5b548-1ddc-11eb-2f7a-0d1b24e3cf48
using Pkg, DrWatson

# ╔═╡ 61d2b2c8-1ddc-11eb-3ae7-a3c92d8a64b3
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Plots
end

# ╔═╡ 0fdee93c-1ddc-11eb-2bde-712854f4263e
md"## Listing 1.11"

# ╔═╡ 61d33860-1ddc-11eb-1c60-eda2aa584f64
function hailLength(x::Int)
    n = 0
    while x != 1
        if x % 2 == 0
            x = Int(x/2)
        else
            x = 3x +1
        end
        n += 1
    end
    return n
end

# ╔═╡ 61dc8cb4-1ddc-11eb-065a-fdfe7782833d
lengths = [hailLength(x0) for x0 in 2:10^7];

# ╔═╡ 61dda098-1ddc-11eb-19ba-c3b600208626
histogram(lengths, bins=1000, normed=:true,
    fill=(:blue, true), la=0, legend=:none,
    xlims=(0, 500), ylims=(0, 0.012),
    xlabel="Length", ylabel="Frequency")

# ╔═╡ 61e6709c-1ddc-11eb-392c-97a365b3a46d
md"## End of listing 1.11"

# ╔═╡ Cell order:
# ╟─0fdee93c-1ddc-11eb-2bde-712854f4263e
# ╠═61a5b548-1ddc-11eb-2f7a-0d1b24e3cf48
# ╠═61d2b2c8-1ddc-11eb-3ae7-a3c92d8a64b3
# ╠═61d33860-1ddc-11eb-1c60-eda2aa584f64
# ╠═61dc8cb4-1ddc-11eb-065a-fdfe7782833d
# ╠═61dda098-1ddc-11eb-19ba-c3b600208626
# ╟─61e6709c-1ddc-11eb-392c-97a365b3a46d


text:
# ------------------------------------
# set volume
# ------------------------------------
function set_volume(nodes, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, volume)
    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                #----------------------------
                #----------------------------
                #----------------------------
                rx = nodes[i+1,j+1,k+1,1] - nodes[i,j,k,1]
                ry = nodes[i+1,j+1,k+1,2] - nodes[i,j,k,2]
                rz = nodes[i+1,j+1,k+1,3] - nodes[i,j,k,3]

                Ax = 0.5 * (vecAx[i,j,k,1] + vecAx[i+1,j,k,1]) +
                     0.5 * (vecAy[i,j,k,1] + vecAy[i,j+1,k,1]) +
                     0.5 * (vecAz[i,j,k,1] + vecAz[i,j,k+1,1])
                Ay = 0.5 * (vecAx[i,j,k,2] + vecAx[i+1,j,k,2]) +
                     0.5 * (vecAy[i,j,k,2] + vecAy[i,j+1,k,2]) +
                     0.5 * (vecAz[i,j,k,2] + vecAz[i,j,k+1,2])
                Az = 0.5 * (vecAx[i,j,k,3] + vecAx[i+1,j,k,3]) +
                     0.5 * (vecAy[i,j,k,3] + vecAy[i,j+1,k,3]) +
                     0.5 * (vecAz[i,j,k,3] + vecAz[i,j,k+1,3])

                volume[i,j,k] = (Ax*rx + Ay*ry + Az*rz) /3
            end
        end
    end
    return volume
end

# ------------------------------------
# set cellcenter
# ------------------------------------
function set_cellcenter(cellcenter, nodes, cellxmax, cellymax, cellzmax)
    for l in 1:3
        for k in 1:cellzmax
            for j in 1:cellymax
                for i in 1:cellxmax
                    xl = 0.125*(nodes[i+1,j+1,  k,l] +
                                nodes[i+1,  j,  k,l] +
                                nodes[  i,j+1,  k,l] +
                                nodes[  i,  j,  k,l] +
                                nodes[i+1,j+1,k+1,l] +
                                nodes[i+1,  j,k+1,l] +
                                nodes[  i,j+1,k+1,l] +
                                nodes[  i,  j,k+1,l])
                    cellcenter[i,j,k,l] = xl
                end
            end
        end
    end
    return cellcenter
end

# ------------------------------------
# set dx and dy for local time stepping
# ------------------------------------
function set_dx_lts(dx, dy, dz, nodes, cellxmax, cellymax, cellzmax, icell)
    # define at cell boundaries
    # dx = zeros(cellxmax+1, cellymax)
    # dy = zeros(cellxmax, cellymax+1)

    point = zeros(3,3)

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax+1 -icell
                #= 点と直線の距離
                   ax+by+cz+d = 0
                   A(x0, y0, z0)
                   => D = abs(ax0 + by0 + cz0 + d) / abs(a^2 + b^2 + c^2)
                =#
                for l in 1:3
                    point[1,l] = nodes[i,  j,  k,l]
                    point[2,l] = nodes[i,j+1,  k,l]
                    point[3,l] = nodes[i,  j,k+1,l]
                end

                a,b,c,d = dis_point_plane(point)

                ccx = 0.125 * (nodes[i,  j,  k,1] + nodes[i+1,  j,  k,1] +
                               nodes[i,j+1,  k,1] + nodes[i+1,j+1,  k,1] +
                               nodes[i,  j,k+1,1] + nodes[i+1,  j,k+1,1] )
                ccy = 0.125 * (nodes[i,  j,  k,2] + nodes[i+1,  j,  k,2] +
                               nodes[i,j+1,  k,2] + nodes[i+1,j+1,  k,2] +
                               nodes[i,  j,k+1,2] + nodes[i+1,  j,k+1,2] )
                ccz = 0.125 * (nodes[i,  j,  k,3] + nodes[i+1,  j,  k,3] +
                               nodes[i,j+1,  k,3] + nodes[i+1,j+1,  k,3] +
                               nodes[i,  j,k+1,3] + nodes[i+1,  j,k+1,3] )

                dx1 = abs(a*ccx + b*ccy + c*ccz + d) / (a^2 + b^2 + c^2)^0.5

                ccx = 0.125 * (nodes[i,  j,  k,1] + nodes[i-1,  j,  k,1] +
                               nodes[i,j+1,  k,1] + nodes[i-1,j+1,  k,1] +
                               nodes[i,  j,k+1,1] + nodes[i-1,  j,k+1,1] )
                ccy = 0.125 * (nodes[i,  j,  k,2] + nodes[i-1,  j,  k,2] +
                               nodes[i,j+1,  k,2] + nodes[i-1,j+1,  k,2] +
                               nodes[i,  j,k+1,2] + nodes[i-1,  j,k+1,2] )
                ccz = 0.125 * (nodes[i,  j,  k,3] + nodes[i-1,  j,  k,3] +
                               nodes[i,j+1,  k,3] + nodes[i-1,j+1,  k,3] +
                               nodes[i,  j,k+1,3] + nodes[i-1,  j,k+1,3] )

                dx2 = abs(a*ccx + b*ccy + c*ccz + d) / (a^2 + b^2 + c^2)^0.5

                dx[i,j,k] = dx1 + dx2
            end
        end
    end

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax+1 -icell
            for i in 1+icell:cellxmax -icell
                for l in 1:3
                    point[1,l] = nodes[  i,  j,  k,l]
                    point[2,l] = nodes[i+1,  j,  k,l]
                    point[3,l] = nodes[  i,  j,k+1,l]
                end

                a,b,c,d = dis_point_plane(point)

                ccx = 0.125 * (nodes[i,  j,  k,1] + nodes[  i,j+1,  k,1] +
                               nodes[i+1,j,  k,1] + nodes[i+1,j+1,  k,1] +
                               nodes[i,  j,k+1,1] + nodes[  i,j+1,k+1,1] )
                ccy = 0.125 * (nodes[i,  j,  k,2] + nodes[  i,j+1,  k,2] +
                               nodes[i+1,j,  k,2] + nodes[i+1,j+1,  k,2] +
                               nodes[i,  j,k+1,2] + nodes[  i,j+1,k+1,2] )
                ccz = 0.125 * (nodes[i,  j,  k,3] + nodes[  i,j+1,  k,3] +
                               nodes[i+1,j,  k,3] + nodes[i+1,j+1,  k,3] +
                               nodes[i,  j,k+1,3] + nodes[  i,j+1,k+1,3] )

                dx1 = abs(a*ccx + b*ccy + c*ccz + d) / (a^2 + b^2 + c^2)^0.5

                ccx = 0.125 * (nodes[i,  j,  k,1] + nodes[  i,j-1,  k,1] +
                               nodes[i+1,j,  k,1] + nodes[i+1,j-1,  k,1] +
                               nodes[i,  j,k+1,1] + nodes[  i,j-1,k+1,1] )
                ccy = 0.125 * (nodes[i,  j,  k,2] + nodes[  i,j-1,  k,2] +
                               nodes[i+1,j,  k,2] + nodes[i+1,j-1,  k,2] +
                               nodes[i,  j,k+1,2] + nodes[  i,j-1,k+1,2] )
                ccz = 0.125 * (nodes[i,  j,  k,3] + nodes[  i,j-1,  k,3] +
                               nodes[i+1,j,  k,3] + nodes[i+1,j-1,  k,3] +
                               nodes[i,  j,k+1,3] + nodes[  i,j-1,k+1,3] )

                dx2 = abs(a*ccx + b*ccy + c*ccz + d) / (a^2 + b^2 + c^2)^0.5

                dy[i,j,k] = dx1 + dx2
            end
        end
    end

    for k in 1+icell:cellzmax+1 -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax -icell
                for l in 1:3
                    point[1,l] = nodes[  i,  j,  k,l]
                    point[2,l] = nodes[i+1,  j,  k,l]
                    point[3,l] = nodes[  i,j+1,  k,l]
                end

                a,b,c,d = dis_point_plane(point)

                ccx = 0.125 * (nodes[i,    j,  k,1] + nodes[  i,  j,k+1,1] +
                               nodes[i+1,  j,  k,1] + nodes[i+1,  j,k+1,1] +
                               nodes[  i,j+1,  k,1] + nodes[  i,j+1,k+1,1] )
                ccy = 0.125 * (nodes[i,    j,  k,2] + nodes[  i,  j,k+1,2] +
                               nodes[i+1,  j,  k,2] + nodes[i+1,  j,k+1,2] +
                               nodes[  i,j+1,  k,2] + nodes[  i,j+1,k+1,2] )
                ccz = 0.125 * (nodes[i,    j,  k,3] + nodes[  i,  j,k+1,3] +
                               nodes[i+1,  j,  k,3] + nodes[i+1,  j,k+1,3] +
                               nodes[  i,j+1,  k,3] + nodes[  i,j+1,k+1,3] )

                dx1 = abs(a*ccx + b*ccy + c*ccz + d) / (a^2 + b^2 + c^2)^0.5

                ccx = 0.125 * (nodes[i,    j,  k,1] + nodes[  i,  j,k-1,1] +
                               nodes[i+1,  j,  k,1] + nodes[i+1,  j,k-1,1] +
                               nodes[  i,j+1,  k,1] + nodes[  i,j+1,k-1,1] )
                ccy = 0.125 * (nodes[i,    j,  k,2] + nodes[  i,  j,k-1,2] +
                               nodes[i+1,  j,  k,2] + nodes[i+1,  j,k-1,2] +
                               nodes[  i,j+1,  k,2] + nodes[  i,j+1,k-1,2] )
                ccz = 0.125 * (nodes[i,    j,  k,3] + nodes[  i,  j,k-1,3] +
                               nodes[i+1,  j,  k,3] + nodes[i+1,  j,k-1,3] +
                               nodes[  i,j+1,  k,3] + nodes[  i,j+1,k-1,3] )

                dx2 = abs(a*ccx + b*ccy + c*ccz + d) / (a^2 + b^2 + c^2)^0.5

                dz[i,j,k] = dx1 + dx2
            end
        end
    end

    return dx, dy, dz
end

# ------------------------------------
# Distance between point and plane
# ------------------------------------
function dis_point_plane(point)
    #= 点と直線の距離
        ax+by+cz+d = 0
        A(x0, y0, z0)
        => D = abs(ax0 + by0 + cz0 + d) / abs(a^2 + b^2 + c^2)

        平面に対する法線ベクトル:vec{n} = (p, q, r)
        p(x−x0) + q(y−y0) + r(z−z0) = 0

        法線ベクトルは外積から求める。
    =#
    a1 = point[2,1] - point[1,1]
    a2 = point[2,2] - point[1,2]
    a3 = point[2,3] - point[1,3]
    b1 = point[3,1] - point[1,1]
    b2 = point[3,2] - point[1,2]
    b3 = point[3,3] - point[1,3]

    p = a2*b3 - a3*b2
    q = a3*b1 - a1*b3
    r = a1*b2 - a2*b1

    a = p
    b = q
    c = r
    d = -p*point[1,1] - q*point[1,2] - r*point[1,3]

    return a, b, c, d
end

# ------------------------------------
# set dtau for local time stepping
# ------------------------------------
function set_lts(dtau, lambda_facex, lambda_facey, lambda_facez, Qbase, cellxmax, cellymax, cellzmax, mu, dx, dy, dz,
                vecAx, vecAy, vecAz, volume, specific_heat_ratio, cfl, icell)
    g = specific_heat_ratio

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax+1 -icell
                rho_av = 0.5 * (Qbase[i,j,k,1] + Qbase[i-1,j,k,1])
                u_av   = 0.5 * (Qbase[i,j,k,2] + Qbase[i-1,j,k,2])
                v_av   = 0.5 * (Qbase[i,j,k,3] + Qbase[i-1,j,k,3])
                w_av   = 0.5 * (Qbase[i,j,k,4] + Qbase[i-1,j,k,4])
                mu_av  = 0.5 * (   mu[i,j,k]   +    mu[i-1,j,k]  )

                ap = (g * Qbase[  i,j,k,5] / Qbase[  i,j,k,1])^0.5
                am = (g * Qbase[i-1,j,k,5] / Qbase[i-1,j,k,1])^0.5
                a_av = 0.5 * (ap + am)

                U   = u_av*vecAx[i,j,k,1] + v_av*vecAx[i,j,k,2] + w_av*vecAx[i,j,k,3]
                lambda_facex[i,j,k] = abs(U) + a_av + 2*mu_av/(rho_av*dx[i,j,k])
            end
        end
    end

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax+1 -icell
            for i in 1+icell:cellxmax -icell
                rho_av = 0.5 * (Qbase[i,j,k,1] + Qbase[i,j-1,k,1])
                u_av   = 0.5 * (Qbase[i,j,k,2] + Qbase[i,j-1,k,2])
                v_av   = 0.5 * (Qbase[i,j,k,3] + Qbase[i,j-1,k,3])
                w_av   = 0.5 * (Qbase[i,j,k,4] + Qbase[i,j-1,k,4])
                mu_av  = 0.5 * (   mu[i,j,k]   +    mu[i,j-1,k]  )

                ap = (g * Qbase[  i,j,k,5] / Qbase[  i,j,k,1])^0.5
                am = (g * Qbase[i,j-1,k,5] / Qbase[i,j-1,k,1])^0.5
                a_av = 0.5 * (ap + am)

                V   = u_av*vecAy[i,j,k,1] + v_av*vecAy[i,j,k,2] + w_av*vecAy[i,j,k,3]
                lambda_facey[i,j,k] = abs(V) + a_av + 2*mu_av/(rho_av*dy[i,j,k])
            end
        end
    end

    for k in 1+icell:cellzmax+1 -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax -icell
                rho_av = 0.5 * (Qbase[i,j,k,1] + Qbase[i,j,k-1,1])
                u_av   = 0.5 * (Qbase[i,j,k,2] + Qbase[i,j,k-1,2])
                v_av   = 0.5 * (Qbase[i,j,k,3] + Qbase[i,j,k-1,3])
                w_av   = 0.5 * (Qbase[i,j,k,4] + Qbase[i,j,k-1,4])
                mu_av  = 0.5 * (   mu[i,j,k]   +    mu[i,j,k-1]  )

                ap = (g * Qbase[  i,j,k,5] / Qbase[  i,j,k,1])^0.5
                am = (g * Qbase[i,j,k-1,5] / Qbase[i,j,k-1,1])^0.5
                a_av = 0.5 * (ap + am)

                W   = u_av*vecAz[i,j,k,1] + v_av*vecAz[i,j,k,2] + w_av*vecAz[i,j,k,3]
                lambda_facez[i,j,k] = abs(W) + a_av + 2*mu_av/(rho_av*dz[i,j,k])
            end
        end
    end

    for k in 1+icell:cellzmax-icell
        for j in 1+icell:cellymax-icell
            for i in 1+icell:cellxmax-icell
                a1 = lambda_facex[  i,  j,  k]
                a2 = lambda_facex[i+1,  j,  k]
                a3 = lambda_facey[  i,  j,  k]
                a4 = lambda_facey[  i,j+1,  k]
                a5 = lambda_facez[  i,  j,  k]
                a6 = lambda_facez[  i,  j,k+1]
                lmax = maximum([a1, a2, a3, a4, a5, a6])

                dtau[i,j,k] = cfl * volume[i,j,k] / lmax
            end
        end
    end
    return dtau
end

# ------------------------------------
# set viscosity by Sutherland's formula
# https://cattech-lab.com/science-tools/sutherland/
# ------------------------------------
function set_mu(mu, Qbase, cellxmax, cellymax, cellzmax,specific_heat_ratio, Rd)
    mu0 = 1.82e-5     # Reference Viscosity, Pa s
    T0  = 293.15      # Reference Temperature, K
    C   = 117         # Sutherland's constant, K

    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                # T = p/(rho Rd )
                T = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd)
                mu[i,j,k] = mu0 * (T/T0)^1.5 * (T0+C)/(T+C)
            end
        end
    end
    return mu
end

# ------------------------------------
# set thermal Conductivity by Sutherland's formula
# https://doi.org/10.11357/jsam1937.37.694
# ------------------------------------
function set_lambda(lambda, Qbase, cellxmax, cellymax, cellzmax, mu, specific_heat_ratio, Rd)
    lam0 = 22.3*10^(-3)  # Reference thermal Conductivity, W/(m K)
    T0   = 273.15        # Reference Temperature, K
    C    = 125           # Sutherland's constant, K

    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                T = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd)
                lambda[i,j,k] = lam0*((T0+C)/(T+C))*(T/T0)^1.5
            end
        end
    end

    return lambda
end

# ------------------------------------
# Conversion from primitive variables to conserved variables
# ------------------------------------
function base_to_conservative(Qbase, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """
    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                Qcon[i,j,k,1] = Qbase[i,j,k,1]
                Qcon[i,j,k,2] = Qbase[i,j,k,1]*Qbase[i,j,k,2]
                Qcon[i,j,k,3] = Qbase[i,j,k,1]*Qbase[i,j,k,3]
                Qcon[i,j,k,4] = Qbase[i,j,k,1]*Qbase[i,j,k,4]
                Qcon[i,j,k,5] = Qbase[i,j,k,5]/(specific_heat_ratio-1) + Qbase[i,j,k,1]*(Qbase[i,j,k,2]^2+Qbase[i,j,k,3]^2+Qbase[i,j,k,3]^2)/2
            end
        end
    end
    return Qcon
end

# ------------------------------------
# Conversion from conserved variables to primitive variables
# ------------------------------------
function conservative_to_base(Qbase, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """
    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                Qbase[i,j,k,1] = Qcon[i,j,k,1]
                Qbase[i,j,k,2] = Qcon[i,j,k,2]/Qcon[i,j,k,1]
                Qbase[i,j,k,3] = Qcon[i,j,k,3]/Qcon[i,j,k,1]
                Qbase[i,j,k,4] = Qcon[i,j,k,4]/Qcon[i,j,k,1]
                Qbase[i,j,k,5] = (Qcon[i,j,k,5]-Qcon[i,j,k,1]*(Qbase[i,j,k,2]^2+Qbase[i,j,k,3]^2+Qbase[i,j,k,4]^2)/2)*(specific_heat_ratio-1)
            end
        end
    end
    return Qbase
end

# ------------------------------------
# Conversion from Cartesian coordinate system to general coordinate system
# ------------------------------------
function setup_Qcon_hat(Qcon, Qcon_hat, cellxmax, cellymax, cellzmax, volume, nval)
    for l in 1:nval
        for k in 1:cellzmax
            for j in 1:cellymax
                for i in 1:cellxmax
                    Qcon_hat[i,j,k,l] = Qcon[i,j,k,l] * volume[i,j,k]
                end
            end
        end
    end
    return Qcon_hat
end

# ------------------------------------
# Conversion from general coordinate system to Cartesian coordinate system
# ------------------------------------
function Qhat_to_Q(Qcon, Qcon_hat, cellxmax, cellymax, cellzmax, volume, nval)
    for l in 1:nval
        for k in 1:cellzmax
            for j in 1:cellymax
                for i in 1:cellxmax
                    Qcon[i,j,k,l] = Qcon_hat[i,j,k,l] / volume[i,j,k]
                end
            end
        end
    end
    return Qcon
end

# ------------------------------------
# cal average convective value
# ------------------------------------
function cal_Qave(Qbase_ave, Qbase, cellxmax, cellymax, cellzmax, nval, stepnum)
    for l in 1:nval
        for k in 1:cellzmax
            for j in 1:cellymax
                for i in 1:cellxmax
                    temp = (Qbase_ave[i,j,k,l]*(stepnum-1) + Qbase[i,j,k,l]) / stepnum
                    Qbase_ave[i,j,k,l] = temp
                end
            end
        end
    end
    return Qbase_ave
end


text:
using ProjectMoana
using Test
import ColorTypes: RGB

include("test_Colors.jl")

include("test_HdrImages.jl")

include("test_Geometry.jl")

include("test_Cameras.jl")

include("test_Shape.jl")

include("test_Surface.jl")

include("test_Renderer.jl")

include("test_PCG.jl")

include("test_Lexer.jl")


text:
@testset "Continuous space" begin

  # Basic model initialization
  space1 = ContinuousSpace(2; periodic = true, extend = (1, 1))
  space2 = ContinuousSpace(2; periodic = false, extend = (1, 1))
  space3 = ContinuousSpace(2; periodic = true, extend = (2, 1))
  space4 = ContinuousSpace(2; periodic = true, extend = (1, 2))
  space5 = ContinuousSpace(2; periodic = true, extend = (2, 2))
  space6 = ContinuousSpace(1; periodic = true, extend = (1,))
  space7 = ContinuousSpace(3; periodic = true, extend = (1,1,1))
  space8 = ContinuousSpace(3; periodic = false, extend = (1,1,1))

  @test space1.D == 2
  @test space2.D == 2
  @test space3.D == 2
  @test space4.D == 2
  @test space5.D == 2
  @test space6.D == 1
  @test space7.D == 3
  @test space8.D == 3

  @test_throws AssertionError ContinuousSpace(2; extend = (-1,1)) # Cannot have negative extent
  @test_throws AssertionError ContinuousSpace(2; extend = [1,1]) # Must be a tuple
  @test_throws AssertionError ContinuousSpace(2; extend = (1,1,1)) # Must be length D
  @test_throws AssertionError ContinuousSpace(2; extend = ("one",1.0)) # Must contain reals

  model1 = ABM(Agent6, space1)
  model2 = ABM(Agent6, space2)
  model3 = ABM(Agent6, space3)
  model4 = ABM(Agent6, space4)
  model5 = ABM(Agent6, space5)
  model6 = ABM(Agent6, space6)

  @test nagents(model1) == 0
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == []

  # add_agent! with no existing agent (the agent is created)
  pos = (0.5, 0.5)
  vel = (0.2, 0.1)
  dia = 0.01
  add_agent!(pos, model1, vel, dia)
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == [1]
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame;
  @test dbrow[1, :a] == 0.5
  @test dbrow[1, :b] == 0.5

  # move_agent! without provided update_vel! function
  move_agent!(model1.agents[1], model1)
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame;
  @test dbrow[1, :a] == 0.7
  @test dbrow[1, :b] == 0.6
  @test dbrow[1, :a] == model1.agents[1].pos[1]
  @test dbrow[1, :b] == model1.agents[1].pos[2]

  kill_agent!(model1.agents[1], model1)
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame
  @test size(dbrow) == (0,0)

  # add_agent! with an existing agent
  agent = Agent6(2, pos, vel, dia)
  add_agent!(agent, model1)
  @test Agents.defvel(agent, model1) == nothing
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == [2]

  # agents within some range are found correctly
  agent2 = model1.agents[2]
  agent3 = Agent6(3, agent2.pos .+ 0.005, vel, dia)
  add_agent_pos!(agent3, model1)
  n_ids = nearby_ids(agent2, model1, agent2.weight)
  @test length(n_ids) == 1
  @test n_ids[1] == 3
  n_ids = nearby_ids(agent2, model1, agent2.weight/10)
  @test length(n_ids) == 0

  # test that it finds both
  n_ids = nearby_ids(agent2.pos, model1, agent2.weight)
  @test sort!(n_ids) == [2, 3]

  # test various metrics
  c = ABM(Agent6, ContinuousSpace(2; extend = (2.0, 2.0), metric = :cityblock))
  e = ABM(Agent6, ContinuousSpace(2; extend = (2.0, 2.0), metric = :euclidean))
  r = sqrt(2) - 0.2
  for (i, φ) in enumerate(range(0; stop = 2π, length = 10))
    a = Agent6(i, (1, 1) .+ r .* sincos(φ), (0.0, 0.0), 0.0)
    add_agent_pos!(a, c)
    add_agent_pos!(a, e)
  end
  @test length(nearby_ids((1, 1), c, r)) == 10
  @test 4 < length(nearby_ids((1, 1), e, r)) < 10
end

mutable struct AgentU1 <: AbstractAgent
  id::Int
  pos::NTuple{2,Float64}
  vel::NTuple{2,Float64}
end

mutable struct AgentU2 <: AbstractAgent
  id::Int
  pos::NTuple{2,Float64}
  vel::NTuple{2,Float64}
end

function ignore_six(model::ABM{A,S,F,P}) where {A,S,F,P}
  [a.id for a in allagents(model) if !(typeof(a) <: Agent6)]
end

@testset "Interacting pairs" begin
  space = ContinuousSpace(2, extend = (10, 10), periodic = false, metric = :euclidean)
  model = ABM(Agent6, space; scheduler = model -> sort!(collect(keys(model.agents));rev=true))
  pos = [
    (7.074386436066224, 4.963014649338054)
    (5.831962448496828, 4.926297135685473)
    (5.122087781793935, 5.300031210394806)
    (3.9715633336430156, 4.8106570045816675)
  ]
  for i in 1:4
    add_agent_pos!(Agent6(i+2, pos[i], (0.0, 0.0), 0), model)
  end
  pairs = interacting_pairs(model, 2.0, :scheduler).pairs
  @test length(pairs) == 2
  fi = [p[1] for p in pairs]
  se = [p[2] for p in pairs]
  @test fi == unique(fi)
  @test se == unique(se)
  for id in fi
    @test id ∉ se
  end
  pairs = interacting_pairs(model, 2.0, :all).pairs
  @test length(pairs) == 5
  @test (3, 6) ∉ pairs

  space2 = ContinuousSpace(2, extend = (10, 10), periodic = false, metric = :euclidean)
  model2 = ABM(Agent6, space2; scheduler = by_id)
  for i in 1:4
    add_agent_pos!(Agent6(i, pos[i], (0.0, 0.0), 0), model2)
  end
  # Note that length here is not the same as the test above with the same function
  # call. This is due to the `scheduler` order of operation
  pairs = interacting_pairs(model2, 2.0, :scheduler).pairs
  @test length(pairs) == 2
  # A more expensive search (in memory, not time), but guarantees true nearest neighbors
  pairs = interacting_pairs(model2, 2.0, :nearest).pairs
  @test length(pairs) == 3
  pairs = interacting_pairs(model2, 2.0, :all).pairs
  @test length(pairs) == 5
  @test (1, 4) ∉ pairs

  space3 = ContinuousSpace(2, extend = (1, 1), periodic = false, metric = :euclidean)
  model3 = ABM(Union{Agent6, AgentU1, AgentU2}, space3; warn = false)
  for i in 1:10
    add_agent_pos!(Agent6(i, (i/10, i/10), (0.0, 0.0), 0), model3)
  end
  for i in 11:20
    add_agent_pos!(AgentU1(i, (i/10-1, 0.5), (0.0, 0.0)), model3)
  end
  for i in 21:30
    add_agent_pos!(AgentU2(i, (0.45, i/10-2), (0.0, 0.0)), model3)
  end
  pairs = interacting_pairs(model3, 0.1, :types).pairs
  @test length(pairs) == 7
  for (a,b) in pairs
      @test typeof(model3[a]) !== typeof(model3[b])
  end
  @test (3, 6) ∉ pairs

  # Test that we have at least some Agent6's in this match
  @test any(typeof(model3[a]) <: Agent6 || typeof(model3[b]) <: Agent6 for (a,b) in pairs)
  pairs = interacting_pairs(model3, 0.2, :types; scheduler = ignore_six).pairs
  @test length(pairs) == 12
  # No Agent6's when using the ignore_six scheduler
  @test all(!(typeof(model3[a]) <: Agent6) && !(typeof(model3[b]) <: Agent6) for (a,b) in pairs)

  # Fix #288
  space = ContinuousSpace(2; periodic = true, extend = (1,1), metric=:euclidean)
  model = ABM(Agent6, space)
  pos = [(0.0, 0.0),(0.2, 0.2),(0.5, 0.5)]
  for i in pos
    add_agent!(i,model,(0.0,0.0),1.0)
  end
  pairs = interacting_pairs(model, .29, :all)
  @test length(pairs) == 1
  (a,b) = first(pairs)
  @test (a.id, b.id) == (1,2)
  # Before the #288 fix, this would return (2,3) as a pair
  # which has a euclidean distance of 0.42
  pairs = interacting_pairs(model, .3, :all)
  @test length(pairs) == 1
  (a,b) = first(pairs)
  @test (a.id, b.id) == (1,2)
end

@testset "nearest neighbor" begin
  space = ContinuousSpace(2; periodic = true, extend = (1,1), metric = :euclidean)
  model = ABM(Agent9, space)
  pos = [(0.0, 0.0),(0.2, 0.0),(0.2, 0.2),(0.5, 0.5)]
  for i in pos
    add_agent!(i,model,(0.0,0.0),nothing)
  end

  for agent in allagents(model)
    agent.f1 = nearest_neighbor(agent, model, sqrt(2)).id
  end

  @test model.agents[1].f1 == 2
  @test model.agents[2].f1 == 1
  @test model.agents[3].f1 == 2
  @test model.agents[4].f1 == 3
end


text:
passports = readlines("day4/input.txt")

passports_dicts = []
current_passport = Dict()
for line in passports
    global current_passport, passports_dicts
    if line == ""
        push!(passports_dicts, current_passport)
        current_passport = Dict()
        continue
    end
    push!(current_passport, [s[1] => s[2] for s in split.(split(line, " "), ":")]...)
end
push!(passports_dicts, current_passport)

println(count(p -> issubset(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"], keys(p)), passports_dicts))

valid_passports = 0
for passport in passports_dicts
    if !issubset(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"], keys(passport))
        continue
    end
    for (k, v) in pairs(passport)
        passport[k] = convert(String, v)
    end
    valid = length(passport["byr"]) == 4 && 1920 <= parse(Int, passport["byr"]) <= 2002
    valid = valid && length(passport["iyr"]) == 4 && 2010 <= parse(Int, passport["iyr"]) <= 2020
    valid = valid && length(passport["eyr"]) == 4 && 2020 <= parse(Int, passport["eyr"]) <= 2030
    valid = valid && begin
        result = false
        match_string = convert(String, passport["hgt"])
        if length(match_string) >= 4 && match_string[end-1:end] == "cm"
            result = 150 <= parse(Int, match_string[1:end-2]) <= 193
        elseif length(match_string) >= 4 && match_string[end-1:end] == "in"
            result = 59 <= parse(Int, match_string[1:end-2]) <= 76
        end
        result
    end
    valid = valid && match(r"^#[0-9a-f]{6}$", passport["hcl"]) !== nothing
    valid = valid && passport["ecl"] in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
    valid = valid && match(r"^[0-9]{9}$", passport["pid"]) !== nothing
    global valid_passports += valid
end

println(valid_passports)

# oh my god that was horrible lol


text:
import Base: dotview, view, unsafe_view, copyto!, eachindex, _maybe_reshape_parent, reshape, to_shape, SubArray, compute_stride1, axes, check_parent_index_match, IndexStyle, LinearIndices
using Base: unalias, index_ndims, @_inline_meta, @boundscheck, ViewIndex, OneTo, rdims, viewindexing, ensure_indexable, index_dimsum, fill_to_length
using Base.Broadcast: Broadcasted

### k[1:2,3:4] .= 0 => materialize!(dotview(k,1:2,3:4),broadcasted(identity,0))

# The following adapted from base: views.jl broadcast.jl subarray.jl
# Much of this will be unnecessary if we can inherit from AbstractArray

dotview(A::KnetArray,I...) = view(A,I...)

function view(A::KnetArray, I::Vararg{Any,N}) where {N}
    @_inline_meta
    J = map(i->unalias(A,i), to_indices(A, I))
    #TODO: @boundscheck checkbounds(A, J...)
    unsafe_view(_maybe_reshape_parent(A, index_ndims(J...)), J...)
end

function unsafe_view(A::KnetArray, I::Vararg{ViewIndex,N}) where {N}
    @_inline_meta
    SubArray(A, I)
end

eachindex(::IndexLinear, A::KnetArray) = (@_inline_meta; OneTo(length(A)))
LinearIndices(A::KnetArray) = LinearIndices(axes(A))

_maybe_reshape_parent(A::KnetArray, ::NTuple{1, Bool}) = reshape(A, Val(1))
_maybe_reshape_parent(A::KnetArray{<:Any,1}, ::NTuple{1, Bool}) = reshape(A, Val(1))
_maybe_reshape_parent(A::KnetArray{<:Any,N}, ::NTuple{N, Bool}) where {N} = A
_maybe_reshape_parent(A::KnetArray, ::NTuple{N, Bool}) where {N} = reshape(A, Val(N))

reshape(parent::KnetArray{T,N}, ndims::Val{N}) where {T,N} = parent
reshape(parent::KnetArray, ndims::Val{N}) where N = reshape(parent, rdims(Val(N), axes(parent)))
reshape(parent::KnetArray, shp::Tuple{Union{Integer,OneTo}, Vararg{Union{Integer,OneTo}}}) = reshape(parent, to_shape(shp))

function SubArray(parent::KnetArray, indices::Tuple)
    @_inline_meta
    SubArray(IndexStyle(viewindexing(indices), IndexStyle(parent)), parent, ensure_indexable(indices), index_dimsum(indices...))
end

IndexStyle(::KnetArray)=IndexLinear()

compute_stride1(parent::KnetArray, I::NTuple{N,Any}) where {N} =
    (@_inline_meta; compute_stride1(1, fill_to_length(axes(parent), OneTo(1), Val(N)), I))

function axes(A::KnetArray{T,N}, d) where {T,N}
    @_inline_meta
    d <= N ? axes(A)[d] : OneTo(1)
end

check_parent_index_match(parent::KnetArray{T,N}, ::NTuple{N, Bool}) where {T,N} = nothing

# dotview(P::KnetArray{T,N},I...) where {T,N} =SubArray{T,N,typeof(P),typeof(I),false}(P,I,0,0)
# check_parent_index_match(parent::KnetArray{T,N}, ::NTuple{N, Bool}) where {T,N} = nothing


text:
using CUDA: DenseCuArray
using CUDA.CUFFT: CuFFTPlan

"""
Function overloads for `batched_mul` provided by `NNlib` so that FFTW Plans can be handled.

For `mul!`, `FFTW.jl` already provides an implementation. However, we need to loop over the batch dimension.
"""

# Extension for CPU Arrays
for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
    # Note: use $FORWARD and $BACKWARD below because of issue #9775
    @eval begin
        function NNlib.batched_mul!(y::StridedArray{$Tc}, p::rFFTWPlan{$Tr,$FORWARD}, x::StridedArray{$Tr})
            assert_applicable(p, x[:,:,1], y[:,:,1]) # no need to check every batch dim
            @inbounds for k ∈ 1:size(y,3)
                @views unsafe_execute!(p, x[:,:,k], y[:,:,k])
            end
            return y
        end
        function NNlib.batched_mul!(y::StridedArray{$Tr}, p::rFFTWPlan{$Tc,$BACKWARD}, x::StridedArray{$Tc})
            assert_applicable(p, x[:,:,1], y[:,:,1]) # no need to check every batch dim
            @inbounds for k ∈ 1:size(y,3)
                @views unsafe_execute!(p, x[:,:,k], y[:,:,k])
            end
            return y
        end
    end
end

# Methods for GPU Arrays, borrowed from CUDA.jl -> fft.jl #490
function NNlib.batched_mul!(y::DenseCuArray{Ty}, p::CuFFTPlan{T,K,false}, x::DenseCuArray{T}
                           ) where {Ty,T,K}
    CUFFT.assert_applicable(p, x[:,:,1], y[:,:,1])
    @inbounds for k ∈ 1:size(y,3)
        @views CUFFT.unsafe_execute!(p, x[:,:,k] ,y[:,:,k])
    end
    return y
end


text:
mutable struct apriltag_t <: LCMType
  id::Int32
  hammingDistance::Int32
  p::SMatrix{4,2,Float64}
  cxy::SVector{2, Float64}
  homography::SMatrix{3,3,Float64}
  pose::SMatrix{4,4,Float64}
  familyName::String
end

@lcmtypesetup(apriltag_t)


text:
function FieldIntegratedEphemeris(arg0::FieldAbsoluteDate, arg1::FieldAbsoluteDate, arg2::FieldAbsoluteDate, arg3::FieldStateMapper, arg4::PropagationType, arg5::FieldDenseOutputModel, arg6::Map, arg7::List, arg8::Vector{JString})
    return FieldIntegratedEphemeris((FieldAbsoluteDate, FieldAbsoluteDate, FieldAbsoluteDate, FieldStateMapper, PropagationType, FieldDenseOutputModel, Map, List, Vector{JString}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function add_event_detector(obj::FieldAbstractAnalyticalPropagator, arg0::FieldEventDetector)
    return jcall(obj, "addEventDetector", void, (FieldEventDetector,), arg0)
end

function clear_events_detectors(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function get_events_detectors(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_frame(obj::FieldIntegratedEphemeris)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "getGeneratedEphemeris", FieldBoundedPropagator, ())
end

function get_initial_state(obj::FieldIntegratedEphemeris)
    return jcall(obj, "getInitialState", FieldSpacecraftState, ())
end

function get_managed_additional_states(obj::FieldAbstractPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_max_date(obj::FieldIntegratedEphemeris)
    return jcall(obj, "getMaxDate", FieldAbsoluteDate, ())
end

function get_min_date(obj::FieldIntegratedEphemeris)
    return jcall(obj, "getMinDate", FieldAbsoluteDate, ())
end

function get_pv_coordinates(obj::FieldIntegratedEphemeris, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function get_pv_provider(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "getPvProvider", FieldPVCoordinatesProvider, ())
end

function is_additional_state_managed(obj::FieldAbstractPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::FieldAbstractAnalyticalPropagator, arg0::FieldAbsoluteDate, arg1::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate, FieldAbsoluteDate), arg0, arg1)
end

function propagate(obj::FieldAbstractPropagator, arg0::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate,), arg0)
end

function reset_initial_state(obj::FieldIntegratedEphemeris, arg0::FieldSpacecraftState)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState,), arg0)
end

function set_attitude_provider(obj::FieldAbstractPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::FieldAbstractPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_master_mode(obj::FieldAbstractPropagator, arg0::FieldOrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (FieldOrekitStepHandler,), arg0)
end

function set_slave_mode(obj::FieldAbstractPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end


text:
using SACFiles
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
include("test_utilities.jl")
import SACTestUtilities

include("header.jl")
include("data_types.jl")
include("io.jl")


text:
###############################################################################
#
# Reflectance related
#
###############################################################################
"""
    BLUE(can_rad::CanopyRads{FT},
         wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the BLUE @ 469 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function BLUE(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return REF_WL(can_rad, wls, FT(469))
end

"""
    EVI(can_rad::CanopyRads{FT},
        wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the EVI, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function EVI(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    _BLUE = BLUE(can_rad, wls);
    _NIR  = NIR(can_rad, wls);
    _RED  = RED(can_rad, wls);

    return FT(2.5) * (_NIR - _RED) / (_NIR + 6 * _RED - FT(7.5) * _BLUE + 1)
end

"""
    EVI2(can_rad::CanopyRads{FT},
        wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the EVI2, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function EVI2(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    _BLUE = BLUE(can_rad, wls);
    _NIR  = NIR(can_rad, wls);
    _RED  = RED(can_rad, wls);

    return FT(2.5) * (_NIR - _RED) / (_NIR + FT(2.4) * _RED + 1)
end

"""
    LSWI(can_rad::CanopyRads{FT},
         wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the LSWI, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function LSWI(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    _BLUE = BLUE(can_rad, wls);
    _NIR  = NIR(can_rad, wls);
    _RED  = RED(can_rad, wls);
    _SWIR = SWIR(can_rad, wls);

    return (_NIR - _SWIR) / (_NIR + _SWIR)
end

"""
    NDVI(can_rad::CanopyRads{FT},
         wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the NDVI, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function NDVI(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    _NIR = NIR(can_rad, wls);
    _RED = RED(can_rad, wls);

    return (_NIR - _RED) / (_NIR + _RED)
end

"""
    NIR(can_rad::CanopyRads{FT},
        wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the NIR @ 858.5 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function NIR(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return REF_WL(can_rad, wls, FT(858.5))
end

"""
    NIRv(can_rad::CanopyRads{FT},
         wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the NIRv, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function NIRv(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    _NIR = NIR(can_rad, wls);
    _RED = RED(can_rad, wls);

    return (_NIR - _RED) / (_NIR + _RED) * _NIR
end

"""
    RED(can_rad::CanopyRads{FT},
        wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the RED @ 645 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function RED(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return REF_WL(can_rad, wls, FT(645))
end

"""
    REF_WL(wls::WaveLengths{FT},
           can_rad::CanopyRads{FT}
           wls::WaveLengths{FT},
           twl::FT
    ) where {FT<:AbstractFloat}

Return the Reflectance, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
- `twl` Target wave length in nm
"""
function REF_WL(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT},
            twl::FT
) where {FT<:AbstractFloat}
    @unpack WL,nWL = wls;
    @unpack alb_obs = can_rad;

    # find the index where twl nm is
    ind = 0;
    for i in 1:(nWL-1)
        if WL[i] <= twl <= WL[i+1]
            ind = i;
            break;
        end
    end

    # warning if ind == 0
    if ind==0
        @warn "target wave length out of bounds, please check the set up!";
        return FT(NaN)
    else
        return (WL[ind+1] - twl) / (WL[ind+1] - WL[ind]) * alb_obs[ind] +
               (twl - WL[ind]) / (WL[ind+1] - WL[ind]) * alb_obs[ind+1]
    end
end

"""
    SWIR(can_rad::CanopyRads{FT},
        wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the SWIR @ 2130 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function SWIR(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return REF_WL(can_rad, wls, FT(2130))
end

###############################################################################
#
# SIF related
#
###############################################################################
"""
    SIF_WL(wls::WaveLengths{FT},
           can_rad::CanopyRads{FT}
           wls::WaveLengths{FT},
           twl::FT
    ) where {FT<:AbstractFloat}

Return the SIF, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
- `twl` Target SIF wave length in nm
"""
function SIF_WL(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT},
            twl::FT
) where {FT<:AbstractFloat}
    @unpack WLF,nWLF = wls;
    @unpack SIF_obs = can_rad;

    # find the index where twl nm is
    ind = 0;
    for i in 1:(nWLF-1)
        if WLF[i] <= twl <= WLF[i+1]
            ind = i;
            break;
        end
    end

    # warning if ind == 0
    if ind==0
        @warn "target wave length out of bounds, please check the set up!";
        return FT(NaN)
    else
        return (WLF[ind+1] - twl) / (WLF[ind+1] - WLF[ind]) * SIF_obs[ind] +
               (twl - WLF[ind]) / (WLF[ind+1] - WLF[ind]) * SIF_obs[ind+1]
    end
end

"""
    SIF_740(can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the SIF @ 740 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function SIF_740(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return SIF_WL(can_rad, wls, FT(740))
end

"""
    SIF_757(can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the SIF @ 757 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function SIF_757(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return SIF_WL(can_rad, wls, FT(757))
end

"""
    SIF_771(can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Return the SIF @ 771 nm, given
- `can_rad` [`CanopyRads`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function SIF_771(
            can_rad::CanopyRads{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    return SIF_WL(can_rad, wls, FT(771))
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
export doImport, doExport
export copyAndPermute, packAndPrepare, unpackAndCombine, checkSize
export releaseViews, createViews, createViewsNonConst
export getMap
export CombineMode, ADD, INSERT, REPLACE, ABSMAX, ZERO

# Note that all packet size information was removed due to the use of julia's
# built in serialization/objects

"""
Tells JuliaPetra how to combine data received from other processes with existing data on the calling process for specific import or export options.

Here is the list of combine modes:
  * ADD: Sum new values into existing values
  * INSERT: Insert new values that don't currently exist
  * REPLACE: REplace existing values with new values
  * ABSMAX: If ``x_{old}`` is the old value and ``x_{new}`` the incoming new value, replace ``x_{old}`` with ``\\max(x_{old}, x_{new})``
  * ZERO: Replace old values with zero
"""
@enum CombineMode ADD=1 INSERT=2 REPLACE=3 ABSMAX=4 ZERO=5

"""
An interface for providing a source when constructing and using multi-vectors and matrices in parallel.

    getMap(::SrcDistObject)
Gets the map of the indices of the object

See [`DistObject`](@ref)
"""
const SrcDistObject = Any

"""
An interface for providing a target when constructing and using multi-vectors, and matrices in parallel.

To support transfers the following methods must be implemented for the combination of source type and the target type

    getMap(::DistObject)
Gets the map of the indices of the object

    checkSizes(source::<:SrcDistObject{GID, PID, LID}, target::<:DistObject{GID, PID, LID})::Bool
Whether the source and target are compatible for a transfer

    copyAndPermute(source::<:SrcDistObject{GID, PID, LID}, target::<:DistObject{GID, PID, LID}, numSameIDs::LID, permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1})
Perform copies and permutations that are local to this process.

    packAndPrepare(source::<:SrcDistObject{GID, PID, LID}, target::<:DistObjectGID, PID, LID}, exportLIDs::AbstractArray{LID, 1}, distor::Distributor{GID, PID, LID})::AbstractArray
Perform any packing or preparation required for communications.  The
method returns the array of objects to export

    unpackAndCombine(target::<:DistObject{GID, PID, LID}, importLIDs::AbstractArray{LID, 1}, imports::AAbstractrray, distor::Distributor{GID, PID, LID}, cm::CombineMode)
Perform any unpacking and combining after communication

See [`SrcDistObject`](@ref)
"""
const DistObject = Any

"""
Returns true if this object is a distributed global
"""
distributedGlobal(obj) = distributedGlobal(getMap(obj))

getComm(obj) = getComm(getMap(obj))

"""
    getMap(::SrcDistObject)
    getMap(::DistObject)
Gets the map of the indices of the object
"""
function getMap end

## import/export interface ##

"""
    doImport(target, source, importer::Import{GID, PID, LID}, cm::CombineMode)
    doImport(target, source, importer::Union{Import{GID, PID, LID}, Nothing}, cm::CombineMode)

Import data into this object using an Import object ("forward mode").
A null `importer` indicates a trivial import
"""
function doImport(source, target, importer::Import{GID, PID, LID},
        cm::CombineMode) where {GID <:Integer, PID <: Integer, LID <: Integer}
    doTransfer(source, target, cm, numSameIDs(importer), permuteToLIDs(importer),
        permuteFromLIDs(importer), remoteLIDs(importer), exportLIDs(importer),
        distributor(importer), false)
end

function doImport(source, target, importer::Union{Import{GID, PID, LID}, Nothing},
        cm::CombineMode) where {GID <:Integer, PID <: Integer, LID <: Integer}
    if importer != nothing
        doImport(source, target, importer, cm)
    end
end

"""
    doExport(target, source, exporter::Export{GID, PID, LID}, cm::CombineMode)
    doExport(target, source, exporter::Union{Export{GID, PID, LID}, Nothing}, cm::CombineMode)

Export data into this object using an Export object ("forward mode")
A null `exporter` indicates a trivial import
"""
function doExport(source, target,
        exporter::Export{GID, PID, LID}, cm::CombineMode) where {
        GID <:Integer, PID <: Integer, LID <: Integer}
    doTransfer(source, target, cm, numSameIDs(exporter), permuteToLIDs(exporter),
        permuteFromLIDs(exporter), remoteLIDs(exporter), exportLIDs(exporter),
        distributor(exporter), false)
end

function doExport(source, target, exporter::Union{Export{GID, PID, LID}, Nothing},
        cm::CombineMode) where {GID <:Integer, PID <: Integer, LID <: Integer}
    if exporter != nothing
        doImport(source, target, exporter, cm)
    end
end

"""
    doImport(source, target, exporter::Export{GID, PID, LID}, cm::CombineMode)
    doImport(source, target, exporter::Union{Export{GID, PID, LID}, Nothing}, cm::CombineMode)

Import data into this object using an Export object ("reverse mode").
A null `exporter` indicates a trivial import.
"""
function doImport(source, target,
        exporter::Export{GID, PID, LID}, cm::CombineMode) where {
            GID <:Integer, PID <: Integer, LID <: Integer}
    doTransfer(source, target, cm, numSameIDs(exporter), permuteToLIDs(exporter),
        permuteFromLIDs(exporter), remoteLIDs(exporter), exportLIDs(exporter),
        distributor(exporter), true)
end

function doImport(source, target, exporter::Union{Export{GID, PID, LID}, Nothing},
        cm::CombineMode) where {GID <:Integer, PID <: Integer, LID <: Integer}
    if exporter != nothing
        doImport(source, target, exporter, cm)
    end
end

"""
    doExport(source, target, importer::Import{GID, PID, LID}, cm::CombineMode)
    doExport(source, target, importer::Union{Import{GID, PID, LID}, Nothing}, cm::CombineMode)

Export data into this object using an Import object ("reverse mode")
A null `importer` indicates a trivial export.
"""
function doExport(source, target,
        importer::Import{GID, PID, LID}, cm::CombineMode) where {
            GID <:Integer, PID <: Integer, LID <: Integer}
    doTransfer(source, target, cm, numSameIDs(importer), permuteToLIDs(importer),
        permuteFromLIDs(importer), remoteLIDs(importer), exportLIDs(importer),
        distributor(importer), true)
end

function doExport(source, target, importer::Union{Import{GID, PID, LID}, Nothing},
        cm::CombineMode) where {GID <:Integer, PID <: Integer, LID <: Integer}
    if importer != nothing
        doImport(source, target, importer, cm)
    end
end

## import/export functionality ##

"""
    checkSizes(source, target)::Bool

Compare the source and target objects for compatiblity.  By default, returns false.  Override this to allow transfering to/from subtypes
"""
checkSizes(source, target) = false

"""
    copyAndPermute(source::<:SrcDistObject{GID, PID, LID}, target::<:DistObject{GID, PID, LID}, numSameIDs::LID, permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1})
Perform copies and permutations that are local to this process.
"""
function copyAndPermute end

"""
    packAndPrepare(source::<:SrcDistObject{GID, PID, LID}, target::<:DistObjectGID, PID, LID}, exportLIDs::AbstractArray{LID, 1}, distor::Distributor{GID, PID, LID})::AbstractArray
Perform any packing or preparation required for communications.  The
method returns the array of objects to export
"""
function packAndPrepare end

"""
    unpackAndCombine(target::<:DistObject{GID, PID, LID}, importLIDs::AbstractArray{LID, 1}, imports::AAbstractrray, distor::Distributor{GID, PID, LID}, cm::CombineMode)
Perform any unpacking and combining after communication
"""
function unpackAndCombine end

"""
    doTransfer(src, target, cm::CombineMode, numSameIDs::LID, permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1}, remoteLIDs::AbstractArray{LID, 1}, exportLIDs::AbstractArray{LID, 1}, distor::Distributor{GID, PID, LID}, reversed::Bool)

Perform actual redistribution of data across memory images
"""
function doTransfer(source, target, cm::CombineMode,
        numSameIDs::LID, permuteToLIDs::AbstractArray{LID, 1},
        permuteFromLIDs::AbstractArray{LID, 1}, remoteLIDs::AbstractArray{LID, 1},
        exportLIDs::AbstractArray{LID, 1}, distor::Distributor{GID, PID, LID},
        reversed::Bool) where {GID <: Integer, PID <: Integer, LID <: Integer}

    if !checkSizes(source, target)
        throw(InvalidArgumentError("checkSize() indicates that the destination " *
                "object is not a legal target for redistribution from the " *
                "source object.  This probably means that they do not have " *
                "the same dimensions.  For example, MultiVectors must have " *
                "the same number of rows and columns."))
    end

    readAlso = true #from TPetras rwo
    if cm == INSERT || cm == REPLACE
        numIDsToWrite = numSameIDs + length(permuteToLIDs) + length(remoteLIDs)
        if numIDsToWrite == numMyElements(getMap(target))
            # overwriting all local data in the destination, so write-only suffices

            #TODO look at FIXME on line 503
            readAlso = false
        end
    end

    #TODO look at FIXME on line 514
    createViews(source)

    #tell target to create a view of its data
    #TODO look at FIXME on line 531
    createViewsNonConst(target, readAlso)

    if numSameIDs + length(permuteToLIDs) != 0
        copyAndPermute(source, target, numSameIDs, permuteToLIDs, permuteFromLIDs)
    end

    # only need to pack & send comm buffers if combine mode is not ZERO
    # ZERO combine mode indicates results are the same as if all zeros were recieved
    if cm != ZERO
        exports = packAndPrepare(source, target, exportLIDs, distor)

        if ((reversed && distributedGlobal(target))
                || (!reversed && distributedGlobal(source)))
            if reversed
                #do exchange of remote data
                imports = resolveReverse(distor, exports)
            else
                imports = resolve(distor, exports)
            end

            unpackAndCombine(target, remoteLIDs, imports, distor, cm)
        end
    end

    releaseViews(source)
    releaseViews(target)

    nothing
end

"""
    createViews(obj)

doTransfer calls this on the source object.  By default it does nothing, but the source object can use this as a hint to fetch data from a compute buffer on an off-CPU decice (such as GPU) into host memory
"""
function createViews(obj)
end

"""
    createViewsNonConst(obj, readAlso::Bool)

doTransfer calls this on the target object.  By default it does nothing, but the target object can use this as a hint to fetch data from a compute buffer on an off-CPU decice (such as GPU) into host memory
readAlso indicates whether the doTransfer might read from the original buffer
"""
function createViewsNonConst(obj, readAlso::Bool)
end

"""
    releaseViews(obj)

doTransfer calls this on the target and source as it completes to allow any releasing of buffers or views.  By default it does nothing
"""
function releaseViews(obj)
end


text:
export
    TrajdataSegment,

    nsteps,
    pull_record,
    pull_continuous_segments,

    sample_random_subinterval

struct TrajdataSegment
    trajdata_index :: Int # index within EvaluationData.trajdatas of the relevant trajdata
    egoid          :: Int # the active vehicle
    frame_lo       :: Int # the starting frame in trajdata (does not count any sort of history)
    frame_hi       :: Int # the ending frame in trajdata, simulation is from frame_lo and propagates until frame_hi
end
function Base.:(==)(a::TrajdataSegment, b::TrajdataSegment)
    a.trajdata_index   == b.trajdata_index &&
    a.egoid            == b.egoid &&
    a.frame_lo         == b.frame_lo &&
    a.frame_hi         == b.frame_hi
end
function Base.show(io::IO, seg::TrajdataSegment)
    println(io, "TrajdataSegment")
    @printf(io, "\ttrajdata_index:  %d\n", seg.trajdata_index)
    @printf(io, "\tegoid:           %d\n", seg.egoid)
    @printf(io, "\tframe_lo:        %d\n", seg.frame_lo)
    @printf(io, "\tframe_hi:        %d\n", seg.frame_hi)
end

nsteps(seg::TrajdataSegment) = seg.frame_hi - seg.frame_lo # total number of sim steps
Records.nframes(seg::TrajdataSegment) = nsteps(seg) + 1 # total number of frames spanned by trajdata segment

function pull_record(seg::TrajdataSegment, trajdata::Trajdata, prime_history::Int=0)
    rec = SceneRecord(nframes(seg)+prime_history, get_timestep(trajdata))
    scene = Scene()

    # prime_history
    for i in prime_history : -1 : 1
        frame = seg.frame_lo - i
        get!(scene, trajdata, frame)
        update!(rec, scene)
    end

    for frame in seg.frame_lo : seg.frame_hi
        get!(scene, trajdata, frame)
        update!(rec, scene)
    end
    rec
end

"""
    sample_random_subinterval(seg::TrajdataSegment, nframes::Int)
Returns a new TrajdataSegment containing a uniformly random subinterval of length nframes
from the given segment.
Throws an error if the TrajdataSegment is too short
"""
function sample_random_subinterval(seg::TrajdataSegment, nframes::Int)
    len = AutoCore.nframes(seg)

    dom_hi = len - nframes + 1
    domain = 1 : dom_hi

    frame_lo = rand(domain) + seg.frame_lo - 1
    frame_hi = frame_lo + nframes - 1

    TrajdataSegment(seg.trajdata_index, seg.egoid, frame_lo, frame_hi)
end

"""
    pull_continuous_segments(trajdata::Trajdata, trajdata_index::Int)
Run through the Trajdata and pull TrajdataSegments for each id for
continuous frames in which they appear
"""
mutable struct SegmentInProgress
    id::Int
    frame_lo::Int
    updated::Bool
end
function pull_continuous_segments(trajdata::Trajdata, trajdata_index::Int)

    segments = TrajdataSegment[]
    segments_in_progress = Dict{Int, SegmentInProgress}()

    for (frameindex, tdframe) in enumerate(trajdata.frames)

        for iip in values(segments_in_progress)
            iip.updated = false
        end

        for frame in tdframe.lo : tdframe.hi
            state = trajdata.states[frame]

            if haskey(segments_in_progress, state.id)
                segments_in_progress[state.id].updated = true
            else
                segments_in_progress[state.id] = SegmentInProgress(state.id, frameindex, true)
            end
        end

        for iip in values(segments_in_progress)
            if !iip.updated
                push!(segments, TrajdataSegment(trajdata_index, iip.id, iip.frame_lo, frameindex-1))
            end
        end
        filter!((id,iip)->iip.updated, segments_in_progress)
    end

    for iip in values(segments_in_progress)
        push!(segments, TrajdataSegment(trajdata_index, iip.id, iip.frame_lo, length(trajdata.frames)))
    end

    segments
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
using Pkg
Pkg.add("JSON2")

max_stars_repo_name: kubetrail/julia

max_issues_repo_name: kubetrail/julia

max_forks_repo_name: kubetrail/julia