text:
#!/usr/bin/env julia

function readbam(bam::String, strandness::String, k::Int64, cluster::Bool; mappingquality=0, verbose=false)

    reader = open(BAM.Reader, bam)
    record = BAM.Record()
    treads = Int64(0)
    passreads = Int64(0)
    pareads = Int64(0)
    pclus = Dict{String,Array{Int16,1}}()
    pclus2 = Dict()
    mal_reads = Int64(0)

    if !(strandness in ["+" "-"])
        println(STDERR, "ERROR: provided srandness is not correct: $strandness, [+/-].")
        exit(1)
    end

    while !eof(reader)
        read!(reader, record)

        if BAM.ismapped(record)
            rfname::String = BAM.refname(record)
            flag = BAM.flag(record)
            mapq = BAM.mappingquality(record)
            treads += 1

            iftrue::Bool = false

            # if pair end
            if flag&1 == 1
                # if strandness +, take first read
                if strandness == "+" && (flag&64==64)
                    iftrue = true
                # if strandness -, take second read
                elseif strandness == "-" && (flag&128==128)
                    iftrue = true
                else
                    iftrue = false
                end
            else
                iftrue = true
            end

            # skip not primary aligned, supplementary alignment.
            if mapq > mappingquality && (flag&256 == 0) && (flag&2048 == 0) && iftrue
                spl = split(BAM.tempname(record),":")
                passreads += 1

                if flag&16 == 0
                    strness = "+"
                    pos = BAM.rightposition(record)
                elseif flag&16 == 16
                    strness = "-"
                    pos = BAM.position(record)
                else
                    strness = "."
                end

                if length(spl) > 1
                    if spl[2] == "A" && strness != "."
                        pareads += 1
                        cps = rfname * "::" * string(pos) * "::" * strness
                        pclus = clust!(pclus, cps, parse(Int16, spl[1]))

                        if cluster
                            pclus2 = cluster!(pclus2, rfname, pos, [parse(Int64, spl[1])],
                                              strness, k; verbose=verbose)
                        end
                    end

                else
                    mal_reads += 1
                    if verbose
                        println(STDERR, "WARNING: READ NAME is NOT FORMATED PROPERLY: ", BAM.tempname(record))
                        println(STDERR, "TOTAL PRIMARY, POLYA, MALFORMATED READS: $passreads $pareads $mal_reads")
                    end
                end
            end
        end
    end

    return treads, passreads, pareads, pclus, pclus2
end

"""
    cluster!(d::Dict, chr::String, pos::Int64, l::Array{Int64,1},
             str::String, k::Int64; verbose::Bool=false)::Dict

    Function adds new termination site (TS) to a given dataset by searching most
    adjacent clusters and if needed correcting clusters to meat definition of
    +/-k from the center.

    # Arguments
    - `d::Dict`: a collection of clusters where key is defined by chromosome,
       mapping position and strand.
    - `chr::String`: Chromosome.
    - `pos::Int64`: Mapping position of TS.
    - `l::Array{Int64,1}`: Lengths of detected PolyA. First values is count.
    - `str::String`: Strand
    - `k::Int64`: k - distance from center allowed.
    - `verbose::Bool=false`: print clustering proceeding to STDOUT.
"""
function cluster!(d::Dict, chr::String, pos::Int64, l::Array{Int64,1},
                  str::String, k::Int64; verbose::Bool=false)::Dict

    if verbose
        println(STDOUT, "SEARCHING FOR: ", chr*"::"*string(pos)*"::"*str, " TS in +/- $k")
    end

    # search if TS fall into any clusters in range of +/- k, rember keys if so.
    a = Array{String,1}()
    v = Array{Int,1}()
    found::Bool = false

    for i in pos-k:pos+k
        if i >= 0
            key = chr*"::"*string(i)*"::"str
            if key in keys(d)
                push!(a, key)
                push!(v, abs(pos-i))
                found = true
                if verbose
                    println(STDOUT, "FOUND EXISTING IN +/- $k RANGE CLUSTER: ", key)
                end
            end
        end
    end

    if !found
        if verbose
            println(STDOUT, "NO EXISTING CLUSTERS FOUND: ", chr*"::"*string(pos)*"::"*str)
            println(STDOUT, "CREATE NEW. CHECK ADJ.")
        end

        # in some cases with recursion l might be already with frequencies.
        if length(l) > 1
            d[chr*"::"*string(pos)*"::"*str] = Dict(pos => l)
            search_adj_clust!(d, pos, chr*"::"*string(pos)*"::"*str,
                              chr, str, k; verbose=verbose)
            return d
        else
            d[chr*"::"*string(pos)*"::"*str] = Dict(pos => vcat([1], l))
            search_adj_clust!(d, pos, chr*"::"*string(pos)*"::"*str,
                              chr, str, k; verbose=verbose)
            return d
        end
    end

    minv = k + 1
    near_center = ""
    ct = 0

    # find which is the nearest cluster if it is and remember its key.
    for i in v
        ct += 1
        if i < minv
            minv = i
            near_center = a[ct]
            if verbose
                println(STDOUT, "NEAREST CLUSTER FOUND TO BE: ", a[ct])
            end

        elseif i == minv
            # means that this TS is in the midle between two clusters
            if verbose
                println(STDOUT, "POSITION TO BE IN THE MIDLE OF 2 CLUSTERS: ",
                        chr*"::"*string(pos - i)*"::"*"str",
                        " --> ",
                        chr*"::"*string(pos)*"::"*"str",
                        " <-- ",
                        chr*"::"*string(i + pos)*"::"*"str"
                        )
            end

            # try to decide by weight of cluster.
            weight1 = 0
            for x in keys(d[chr*"::"*string(pos-i)*"::"*str])
                weight1 += d[chr*"::"*string(pos-i)*"::"*str][x][1]
            end

            weight2 = 0
            for x in keys(d[chr*"::"*string(pos+i)*"::"*str])
                weight2 += d[chr*"::"*string(pos+i)*"::"*str][x][1]
            end

            if weight1 > weight2
                near_center = chr*"::"*string(pos-i)*"::"*str
            elseif weight2 > weight1
                near_center = chr*"::"*string(pos+i)*"::"*str

            else
                # in some cases with recursion l might be already with frequencies.
                if length(l) > 1
                    d[chr*"::"*string(pos)*"::"*str] = Dict(pos => l)
                    search_adj_clust!(d, pos, chr*"::"*string(pos)*"::"*str,
                                      chr, str, k; verbose=verbose)
                    return d

                else
                    d[chr*"::"*string(pos)*"::"*str] = Dict(pos => vcat([1], l))
                    search_adj_clust!(d, pos, chr*"::"*string(pos)*"::"*str,
                                      chr, str, k; verbose=verbose)
                    return d
                end
            end
        end
    end

    # search in a given cluster
    # check if TS is in this cluster
    # else add new
    # check center and recenter if needed.

    if pos in keys(d[near_center])
        if verbose
            println(STDOUT, "FOUND SAME TS IN THE CLUSTER: ", pos, " in ", near_center)
        end

        d[near_center][pos][1] += 1
        d[near_center][pos] = vcat(d[near_center][pos], l)
        recenter!(d, d[near_center], near_center, k, chr, str; verbose=verbose)

    else
        if verbose
            println(STDOUT, "NO TS FOUND IN THE CLUSTER, ADDING NEW: ", pos, " in ", near_center)
        end

        # in some cases with recursion l might be already with frequencies.
        if length(l) > 1
            d[near_center][pos] = l
            recenter!(d, d[near_center], near_center, k, chr, str; verbose=verbose)
        else
            d[near_center][pos] = vcat([1], l)
            recenter!(d, d[near_center], near_center, k, chr, str; verbose=verbose)
        end
    end

    return d
end

function recenter!(d1::Dict, d2::Dict, near_center::String, k::Int64,
                   chr::String, str::String; verbose::Bool=false)::Dict

    # find new most frequent element
    center_pos = parse(Int64, split(near_center, "::")[2])

    if center_pos in keys(d2)
        center_freq = d2[center_pos][1]
    else
        if verbose
            println(STDOUT, "CENTER NOT FOUND, MUST HAVE CHANGED IN RECURSION: $center_pos")
        end

        return d1
    end

    new_center = 0
    new_center_freq = 0
    found_new_center::Bool = false
    equals = Int64[]
    all_positions = Int64[]

    if verbose
        println(STDOUT, "CHECK CENTER IF PERSISTS.")
    end

    for position in keys(d2)
        push!(all_positions, position)

        if position != center_pos
            if !found_new_center && d2[position][1] > center_freq
                if verbose
                    println(STDOUT, "PREVIOUS POSITION: $center_pos with FREQUENCY: $center_freq")
                    println(STDOUT, "NEW POSITION: $position with FREQUENCY: ", d2[position][1])
                end
                new_center = position
                new_center_freq = d2[position][1]
                found_new_center = true
                equals = []

            elseif !found_new_center && d2[position][1] == center_freq
                push!(equals, position)
                if verbose
                    println(STDOUT, "FOUND EQUAL FREQUENCIES: ", equals)
                end

            elseif found_new_center && d2[position][1] > new_center_freq
                new_center = position
                new_center_freq = d2[position][1]
                equals = []
                if verbose
                    println(STDOUT, "PREVIOUS POSITION: $center_pos with FREQUENCY: $center_freq")
                    println(STDOUT, "NEW POSITION: $position with FREQUENCY: ", d2[position][1])
                end

            elseif found_new_center && d2[position][1] == new_center_freq
                push!(equals, position)
            end
        end
    end

    # find witch one is nearest to median if found > 1 most frequent
    # and select it as center.
    if length(equals) > 0
        med = round(median(all_positions))

        if found_new_center
            diff = abs(new_center - med)
        else
            diff = abs(center_pos - med)
        end

        for i in equals
            if abs(i - med) < diff
                new_center = i
                found_new_center = true
                if verbose
                    println(STDOUT, "FOUND NEW CENTER MOST ADJ TO MEDIAN: ", i)
                end
            end
        end
    end

    if found_new_center
        new_cluster_id = chr*"::"*string(new_center)*"::"*str
    else
        # for checking adjacent clusters even if new center not found.
        new_cluster_id = near_center
        new_center = center_pos
    end

    # != arrises from recursive call when two adj clusters are merged.
    # in this case we want to cat all values of TS.
    if d1[near_center] != d2
        if verbose
            println(STDOUT, "CAT VALUES of $d2 to ", d1[near_center])
        end

        for position in keys(d2)
            if position in keys(d1[near_center])
                d1[near_center][position] = vcat(d2[position], d1[near_center][position][2:end])
                d1[near_center][position][1] += d2[position][1]
            else
                d1[near_center][position] = d2[position]
            end
        end

        if verbose
            println(STDOUT, "CATTED ARRAY: ", d1[near_center])
        end

        if found_new_center
            d1[new_cluster_id] = d1[near_center]
            delete!(d1, near_center)
            recenter!(d1, d1[new_cluster_id], new_cluster_id, k, chr, str; verbose=verbose)

            return d1

        else

            return d1
        end
    end

    # check if values of a cluster is still holds +/- k
    # if not - reject those TS to recursion.
    # evaluate near stand alone TS and join if holds.
    rejected_TS1 = Dict{Int64,Array{Int64,1}}()
    if found_new_center
        for i in all_positions
            if i > new_center + k || i < new_center - k

                if verbose
                    println(STDOUT, "REJECTING POSITION FOR RECURSION: ", i)
                end

                rejected_TS1[i] = d2[i]
                if verbose
                    println(STDOUT, "!!! DELETING: $i --> ", d2[i])
                end

                delete!(d2, i)
            end
        end

        # correct main collection of cluster.
        if verbose
            println(STDOUT, "!!! DELETING FROM MAIN DATASET: $near_center --> ", d1[near_center])
            println(STDOUT, "RECREATE IN MAIN DATASET: $new_cluster_id --> ", d2)
        end

        delete!(d1, near_center)
        d1[new_cluster_id] = d2
    end

    search_adj_clust!(d1, new_center, new_cluster_id, chr, str, k; verbose=verbose)

    # recursively correct rejected
    if found_new_center
        for position in keys(rejected_TS1)

            if verbose
                println(STDOUT, "CLUSTER RECURSION OF REJECTED VALUES:")
                println(STDOUT, "NEW CLUSTER: $new_cluster_id, POSITION FOR RECURSION: $position")
            end

            cluster!(d1, chr, position, rejected_TS1[position], str, k; verbose=verbose)
            if verbose
                println(STDOUT, "RETURNED FROM CLUSTER RECURSION2.")
            end
        end
    end

    if verbose && !found_new_center
        println(STDOUT, "NEW CENTER WAS NOT FOUND.")
    end

    return d1
end

function search_adj_clust!(d1, new_center::Int64, new_cluster_id::String,
                           chr::String, str::String, k::Int64; verbose=false)

    # new cluster id might change from recursion.
    if !(new_cluster_id in keys(d1))
        return d1
    end

    found_adj_clusters = String[]
    distancies = Int64[]
    adj_centers = Int64[]
    found::Bool = false
    adj_cluster_id = ""
    adj_center = -1

    # search for adjacent clusters
    for adj_center in new_center-k*2:new_center+k*2

        adj_cluster_id = chr*"::"*string(adj_center)*"::"str

        if adj_cluster_id != new_cluster_id && adj_cluster_id in keys(d1)
            if verbose
                println(STDOUT, "FOUND ADJ CLUSTER IN +/- 2*$k RANGE: $adj_cluster_id near $new_cluster_id")
                println(STDOUT, "$adj_cluster_id:")
            end

            push!(found_adj_clusters, adj_cluster_id)
            push!(distancies, abs(new_center - adj_center))
            push!(adj_centers, adj_center)
            found = true
        end
    end

    if !found
        return d1
    end

    # search wich is more adjacent
    minv = 2*k+1
    ct = 0

    for i in distancies
        ct += 1
        if i < minv
            minv = i
            adj_cluster_id = found_adj_clusters[ct]
            adj_center = adj_centers[ct]

            if verbose
                println(STDOUT, "ADJ CLUSTER FOUND TO BE: ", found_adj_clusters[ct])
            end

        elseif i == minv
            if verbose
                println(STDOUT, "POSITION TO BE IN THE MIDLE OF 2 ADJ CLUSTERS: ",
                        chr*"::"*string(new_center - i)*"::"*"str",
                        " --> ",
                        chr*"::"*string(new_center)*"::"*"str",
                        " <-- ",
                        chr*"::"*string(i + new_center)*"::"*"str"
                        )
            end

            # try to decide by weight of cluster.
            weight1 = 0
            for x in keys(d1[chr*"::"*string(new_center-i)*"::"*str])
                weight1 += d1[chr*"::"*string(new_center-i)*"::"*str][x][1]
            end

            weight2 = 0
            for x in keys(d1[chr*"::"*string(new_center+i)*"::"*str])
                weight2 += d1[chr*"::"*string(new_center+i)*"::"*str][x][1]
            end

            if weight1 > weight2
                adj_cluster_id  = chr*"::"*string(new_center-i)*"::"*str
            elseif weight2 > weight1
                adj_cluster_id  = chr*"::"*string(new_center+i)*"::"*str
            end
        end
    end

    ck = collect(keys(d1[adj_cluster_id]))
    ck2 = collect(keys(d1[new_cluster_id]))
    min = minimum(ck)
    max = maximum(ck)

    if verbose
        println(STDOUT, d1[adj_cluster_id])
        println(STDOUT, ck)
    end

    # if found cluster falls in new cluster - add and check center
    # with recursion.
    if min in new_center-k:new_center+k && max in new_center-k:new_center+k
        if verbose
            println(STDOUT, "FOUND CLUSTER FALLS IN: $adj_cluster_id --> $new_cluster_id")
            println(STDOUT, "RECENTERING WITH NEW CLUSTER ID: $new_cluster_id")
        end

        recenter!(d1, d1[adj_cluster_id], new_cluster_id, k, chr, str; verbose=verbose)

        if verbose
            println(STDOUT, "RETURNED FROM RECENTER RECURSION: $new_cluster_id")
        end

    # if whole cluster does not fall, still some
    # values might be more adjacent to either cluster center
    else
        if verbose
            println(STDOUT, "FOUND CLUSTER DID NOT FALL IN: $adj_cluster_id --> $new_cluster_id")
            println(STDOUT, "CHECKING IF ANY VALUES ARE MORE ADJACENT.")
        end

        check_adj_values!(d1, ck, new_center, adj_center,
                          adj_cluster_id, chr, str,
                          k; verbose=verbose)

        check_adj_values!(d1, ck2, adj_center, new_center,
                          new_cluster_id, chr, str,
                          k; verbose=verbose)

        if verbose
            println(STDOUT, "FINISHED CHECKING ADJ VALUES.")
        end
    end

    return d1
end

function check_adj_values!(d1::Dict, ck::Array{Int64,1}, new_center::Int64,
                           adj_center::Int64, adj_cluster_id::String,
                           chr::String, str::String, k::Int64; verbose=false)

    for position in ck
        if position != new_center && abs(adj_center - position) > abs(new_center - position)
            if verbose
                println(STDOUT, "FOUND VALUE TO BE MORE ADJ TO $new_center RATHER THEN $adj_cluster_id: $position")
                println(STDOUT, ck)

                try
                    println(STDOUT, d1[adj_cluster_id])
                catch
                    KeyError(adj_cluster_id)
                    println(STDERR, "$adj_cluster_id WAS ALREADY REMOVED.")
                end
            end

            # sometimes keys are deleted in recursion
            if adj_cluster_id in keys(d1)
                if position in keys(d1[adj_cluster_id])
                    val = values(d1[adj_cluster_id][position])
                    if verbose
                        println(STDOUT, "!!! DELETE $position in $adj_cluster_id")
                        println(STDOUT, "AND PUSH IT FOR RECURSION.")
                    end

                    delete!(d1[adj_cluster_id], position)
                    cluster!(d1, chr, position, val, str, k; verbose=verbose)

                    if verbose
                        println(STDOUT, "RETURNED FROM CLUSTER RECURSION1")
                    end

                elseif verbose
                    println(STDOUT, "POSITION $position WAS ALREADY REMOVED.")
                end

            elseif verbose
                println(STDOUT, "CLUSTER CENTER $adj_cluster_id WAS ALREADY REMOVED.")
            end
        end
    end

    return d1
end

"""
    merge_adj_clusters!(d::Dict, m::Int64; verbose=false)::Dict
    Searches for adj clusters with <= m distance and merge.
    # Arguments
    - `d::Dict`: a collection of clusters where key is defined by chromosome,
       mapping position and strand.
    - `m::Int64`: a min distance allowed between adjacent cluster values.
    - `verbose::Bool=false`: print clustering proceeding to STDOUT.
"""
function merge_adj_clusters!(d::Dict, m=Int64; verbose::Bool=false)::Dict

    chroms = Dict()

    # collect and sort centers
    for center in collect(keys(d))
        csplit = split(center,"::")
        chr = csplit[1]
        center_pos = parse(Int64,csplit[2])
        str = csplit[3]
        if chr*"::"*str in keys(chroms)
            chroms[chr*"::"*str] = push!(chroms[chr*"::"*str], center_pos)
            sort!(chroms[chr*"::"*str])
        else
            chroms[chr*"::"*str] =  [center_pos]
        end
    end

    # check distancies between clusters and merge if <= m
    for cs in keys(chroms)
        csplit = split(cs,"::")
        chr = csplit[1]
        str = csplit[2]
        ct = 1

        for idx in 1:length(chroms[cs])

            if length(chroms[cs]) == ct
                break
            end

            center_pos = chroms[cs][ct]
            positions = collect(keys(d[chr*"::"*string(center_pos)*"::"str]))
            right = maximum(positions)

            # check next adj center min value
            positions2 = collect(keys(d[chr*"::"*string(chroms[cs][ct+1])*"::"str]))
            left = minimum(positions2)

            if left - right <= m
                if verbose
                    println(STDOUT, "")
                    println(STDOUT, "FOUND TWO ADJ CLUSTERS .... WITH DISTANCE <= $m .... ")
                    println(STDOUT, "MERGING $chr::",chroms[cs][ct+1],"::$str WITH $chr::$center_pos::$str .... ")
                end

                temp_cluster::Dict = merge(d[chr*"::"*string(center_pos)*"::"str],
                                           d[chr*"::"*string(chroms[cs][ct+1])*"::"str])

                new_center_pos = 0
                most_fp = Int64[]
                most_fv = 0
                ct2 = 0
                all_c = Int64[]

                # get most frequent positions
                for k in keys(temp_cluster)
                    push!(all_c, k)
                    if temp_cluster[k][1] > most_fv
                        most_fp = Int64[]
                        most_fv = temp_cluster[k][1]
                        push!(most_fp, k)

                    elseif temp_cluster[k][1] == most_fv
                        push!(most_fp, k)
                    end
                end

                if verbose
                    println(STDOUT, "MOST FREQUENT POSITIONS IN A MERGED CUSTER:")
                    println(STDOUT, "$most_fp")
                end

                # search for new center
                if length(most_fp) > 1
                    med = median(all_c)
                    most_fp_m = most_fp .- med
                    unqs = unique(most_fp_m)
                    l = length(most_fp_m)

                    if unqs == l
                        imin = indmin(most_fp_m)
                        new_center_pos = most_fp[imin]

                        if verbose
                            println(STDOUT, "UNIQUE: ", most_fp[imin])
                        end

                    else
                        if verbose
                            println(STDOUT, "NO UNIQUE CENTERS FOUND .... ")
                            println(STDOUT, "SELECTING MEDIAN POSITION ....")
                        end

                        new_center_pos = trunc(Int64, round(med, 0))

                        # if median is non existing position
                        if length(find(all_c .== new_center_pos)) == 0
                            if verbose
                                println(STDOUT, "MEDIAN POSITION IS EMPTY ....")
                                println(STDOUT, "SELECTING MOST WEIGHTED CENTER ....")
                            end

                            # get sizes of original clusters
                            # center will be of most weighted cluster
                            size1 = cluster_size(d[chr*"::"*string(center_pos)*"::"str])
                            size2 = cluster_size(d[chr*"::"*string(chroms[cs][ct+1])*"::"str])

                            if size1 > size2
                                if verbose
                                    println(STDOUT, "FOUND WEIGHTED NEW CENTER: $center_pos ....")
                                end
                                new_center_pos = center_pos

                            elseif size2 > size1
                                if verbose
                                    println(STDOUT, "FOUND WEIGHTED NEW CENTER: ", chroms[cs][ct+1], " ....")
                                end
                                new_center_pos = chroms[cs][ct+1]

                            else
                                if verbose
                                    println(STDOUT, "NO WEIGHTED CENTER FOUND ....")
                                end

                                if new_center_pos%2 == 0
                                    if verbose
                                        println(STDOUT, "SELECTED $center_pos ....")
                                    end
                                    new_center_pos = center_pos

                                else
                                    if verbose
                                        println(STDOUT, "SELECTED ", chroms[cs][ct+1], " ....")
                                    end
                                    new_center_pos = chroms[cs][ct+1]
                                end
                            end
                        end
                    end

                else
                    new_center_pos = most_fp[1]
                end

                new_center = chr*"::"*string(new_center_pos)*"::"*str
                if verbose
                    println(STDOUT, "MERGED NEW CLUSTER ID: $new_center ....")
                end

                delete!(d, chr*"::"*string(center_pos)*"::"str)
                if verbose
                    println(STDOUT, "REMOVED FIRST PREVIOUS CLUSTER: $chr::$center_pos::str ....")
                end

                delete!(d, chr*"::"*string(chroms[cs][ct+1])*"::"str)
                if verbose
                    println(STDOUT, "REMOVED SECOND PREVIOUS CLUSTER: ", chr*"::"*string(chroms[cs][ct+1])*"::"str, " ....")
                end

                d[new_center] = temp_cluster
                if verbose
                    println(STDOUT, "REMOVING CENTER FROM CHROMS STRUCTURE: ", chroms[cs][ct], " ....")
                end
                deleteat!(chroms[cs], ct)

                if verbose
                    println(STDOUT, "REMOVING CENTER FROM CHROMS STRUCTURE: ", chroms[cs][ct], " ....")
                end
                deleteat!(chroms[cs], ct)

                if verbose
                    println(STDOUT, "ADDING NEW CENTER TO CHROMS STRUCTURE: $new_center_pos ....")
                end
                push!(chroms[cs], new_center_pos)

                if verbose
                    println(STDOUT, "SORTING CHROMS STRUCTURE: ....")
                end
                sort!(chroms[cs])

            else
                if verbose
                    println(STDOUT, "NO ADJ CLUSTER FOUND TO $chr::$center_pos::$str")
                end
                ct += 1
            end
        end
    end

    return d
end

function cluster_size(d::Dict)::Int64

    size = 0

    for k in keys(d)
        size += d[k][1]
    end

    return size
end

function clust!(d::Dict, p::Any, l::Any)::Dict

    if haskey(d, p)
        v = d[p]
        d[p] = push!(v, l)
    else
        d[p] = [l]
    end

    return d
end

"""
    stats_poly_a(d::Dict{String,Array{Int16,1}})

    Collects simple statistics on polyA sites.
    Returns a dataframe.
"""
function stats_poly_a(d::Dict{String,Array{Int16,1}})::DataFrame
    dframe = DataFrame(Chrmosome=String[], Position=Int32[], Strand=String[], Median=Float16[], Minimum=Int16[], Maximum=Int16[], Counts=Int32[])

    for k in keys(d)
        md = median(d[k])
        mn = minimum(d[k])
        mx = maximum(d[k])
        l = length(d[k])
        push!(dframe, [split(k,"::")[1] parse(Int32,split(k,"::")[2]) split(k,"::")[3] md mn mx l])
    end

    return dframe
end

function stats_clusters(d::Dict)::DataFrame
    dframe = DataFrame(Chrmosome=String[], Start=Int64[], End=Int64[],
                       Strand=String[], ClusterCenter=Int64[], ClusterSize=Int32[],
                       ClusterMedian=Float64[], ClusterMean=Float64[],
                       ClusterMin=Int64[], ClusterMax=Int64[],
                       Cluster1stQuartile=Float64[], Cluster3rdQuartile=Float64[],
                       TSMedian=Float64[], TSMean=Float64[],
                       TSMin=Int64[], TSMax=Int64[],
                       TS1stQuartile=Float64[], TS3rdQuartile=Float64[])

    for k in keys(d)
        id = split(k,"::")
        chr = id[1]
        center = parse(Int64, id[2])
        str = id[3]

        all_positions = Int64[]
        frequencies = Int64[]
        ts_lengths = Int64[]
        for position in keys(d[k])
            push!(all_positions, position)
            push!(frequencies, d[k][position][1])
            ts_lengths = vcat(ts_lengths, d[k][position][2:end])
         end

        start = minimum(all_positions)
        end_  = maximum(all_positions)
        md = median(frequencies)
        avg = mean(frequencies)
        mn = minimum(frequencies)
        mx = maximum(frequencies)
        fq = nquantile(frequencies,4)[2]
        tq = nquantile(frequencies,4)[4]
        md2 = median(ts_lengths)
        avg2 = mean(ts_lengths)
        mn2 = minimum(ts_lengths)
        mx2 = maximum(ts_lengths)
        fq2 = nquantile(ts_lengths,4)[2]
        tq2 = nquantile(ts_lengths,4)[4]
        push!(dframe, [chr start end_ str center sum(frequencies) md avg mn mx fq tq md2 avg2 mn2 mx2 fq2 tq2])
    end

    return dframe
end

function parseGFF3(gff3file::String)::Dict{String, IntervalCollection{String}}
    intcol = Dict{String, IntervalCollection{String}}()

    for record in open(GFF3.Reader, gff3file)
        chr::String = GFF3.seqid(record)
        start::Int = GFF3.seqstart(record)
        seqend::Int = GFF3.seqend(record)
        strand::Strand = GFF3.strand(record)
        split_k = get_split_key(chr, start, seqend)
        mdstr = parserecord(record)

        for i in split_k
            if haskey(intcol, i)
                push!(intcol[i], Interval(chr, start, seqend, strand, mdstr))
            else
                intcol[i] = IntervalCollection{String}()
                push!(intcol[i], Interval(chr, start, seqend, strand, mdstr))
            end
        end
    end

    return intcol
end

"""
    get_split_key(chr::String, x::Int64, y::Int64; step::Int64=10000)

    Calculates ranges depending on the step and gff3 coordinates.
"""
function get_split_key(chr::String, x::Int64, y::Int64; step::Int64=10000)::Array{String,1}

    a = Array{String,1}()

    function get_interval(x::Int64,y::Int64)
        if x >= y
            if x%y == 0
                ix1 = x - y + 1
                ix2 = x
            else
                ix1 = x - x%y + 1
                ix2 = x - x%y + y
            end
        else
            ix1 = 1
            ix2 = y
        end

        return ix1:ix2
    end

    interv = get_interval(x, step)

    if y in interv
        push!(a, chr*":"*string(interv))
    else
        push!(a, chr*":"*string(interv))
        interv2 = get_interval(y, step)

        if interv[end] == interv2[1]
            push!(a, chr*":"*string(interv2))
        else
            for i in interv[end]:step:interv2[1]-1
                ix1=i+1
                ix2=i+step
                push!(a, chr*":"*string(ix1:ix2))
            end
        end
    end

    return a
end

function parserecord(r::GFF3.Record)::String

    md::String = ""

    try s = GFF3.attributes(r, "ID")
        md = string(md,"ID=")
        for x in s
            if x == s[end]
                md = string(md,x)
            else
                md = string(md,x,",")
            end
        end
        md = string(md,";")
    catch md = string(md,"ID=;")
    end

    try s = GFF3.attributes(r, "Parent")
        md = string(md,"Parent=")
        for x in s
            if x == s[end]
                md = string(md,x)
            else
                md = string(md,x,",")
            end
        end
        md = string(md,";")
    catch md = string(md,"Parent=;")
    end

    try s = GFF3.attributes(r, "gene_id")
        md = string(md,"gene_id=")
        for x in s
            if x == s[end]
                md = string(md,x)
            else
                md = string(md,x,",")
            end
        end
        md = string(md,";")
    catch md = string(md,"gene_id=;")
    end

    try s = GFF3.attributes(r, "transcript_id")
        md = string(md,"transcript_id=")
        for x in s
            if x == s[end]
                md = string(md,x)
            else
                md = string(md,x,",")
            end
        end
        md = string(md,";")
    catch md = string(md,"transcript_id=;")
    end

    try s = GFF3.attributes(r, "gene_type")
        md = string(md,"gene_type=")
        for x in s
            if x == s[end]
                md = string(md,x)
            else
                md = string(md,x,",")
            end
        end
        md = string(md,";")
    catch md = string(md,"gene_type=;")
    end

    try s = GFF3.attributes(r, "gene_name")
        md = string(md,"gene_name=")
        for x in s
            if x == s[end]
                md = string(md,x)
            else
                md = string(md,x,",")
            end
        end
        md = string(md,";")
    catch md = string(md,"gene_name=;")
    end

    md = string(md,"ftype=",GFF3.featuretype(r))

    return md
end

"""
    getintervals(dframe::DataFrame)

    Returns a collection of intervals from a polyA sites dataframe.
"""
function getintervals(dframe::DataFrame)::IntervalCollection{String}

    intcol = IntervalCollection{String}()

    for i in eachrow(dframe)

        if i[3] == "+"
            s = Strand('+')
        elseif i[3] == "-"
            s = Strand('-')
        elseif i[3] == "."
            s = Strand('.')
        else
            s = Strand('?')
        end

        push!(intcol, Interval(i[1], i[2], i[2], s,
              "Median="*string(i[4])*";"*"Min="*string(i[5])*
              ";"*"Max="*string(i[6])*";"*"Counts="*string(i[7])))
    end

    return intcol
end

"""
    annotate_polya_sites(a::IntervalCollection, b::Dict{String, IntervalCollection{String}})

    Function searches for overlaping intervals. If overlapped - adds annotation
    to the dataframe. If not - add NA.

    # Arguments
    - `a::IntervalCollection`: a collection of intervals of PolyA sites.
    - `b::Dict{String, IntervalCollection{String}}`: a dictionary with a sorted collections of intervals from gff3.
"""
function annotate_polya_sites(a::IntervalCollection,
                              b::Dict{String, IntervalCollection{String}};
                              verbose::Bool=false)::DataFrame

    # Two dataframes for '+' and '-' strands.
    dfp = DataFrame(Chr=String[], Start=Int64[], End=Int64[], Name=String[],
                    Counts=Int32[], Strand=String[], Feature=String[],
                    Median=Float32[], Min=Int16[], Max=Int16[], Biotype=String[])
    dfn = DataFrame(Chr=String[], Start=Int64[], End=Int64[], Name=String[],
                    Counts=Int32[], Strand=String[], Feature=String[],
                    Median=Float32[], Min=Int16[], Max=Int16[], Biotype=String[])
    # Dataframe for not annotated sites.
    dfna = DataFrame(Chr=String[], Start=Int64[], End=Int64[], Name=String[],
                    Counts=Int32[], Strand=String[], Feature=String[],
                    Median=Float32[], Min=Int16[], Max=Int16[], Biotype=String[])

    for i1 in a
        ct = Int64(0)
        # Parse metadata from polyA site.
        splt1::Array{SubString{String},1} = split(metadata(i1),";")
        c = parse(Int32,split(splt1[4],"=")[2])
        str = string(strand(i1))
        me = parse(Float32,split(splt1[1],"=")[2])
        mi = parse(Int16,split(splt1[2],"=")[2])
        mx = parse(Int16,split(splt1[3],"=")[2])
        chr = seqname(i1)
        it1 = first(i1)
        it2 = last(i1)

        # Assuming that start and end is the same of the polyA site.
        split_k = get_split_key(chr, it1, it2)[1]

        try
            for i2 in b[split_k]
                if chechoverlap(i1, i2)
                    # parse metadata from GFF3 anotation
                    splt::Array{SubString{String},1} = split(metadata(i2),";")
                    ft = split(splt[7],"=")[2]
                    gn = split(splt[6],"=")[2]
                    bt = split(splt[5],"=")[2]
                    ct += 1

                    if strand(i1) == Strand('-')
                        push!(dfn, [seqname(i1) it1-1 it2 gn c str ft me mi mx bt])
                    else
                        push!(dfp, [seqname(i1) it1-1 it2 gn c str ft me mi mx bt])
                    end
                end
            end
        catch
            KeyError(split_k)

            if verbose
                println(STDERR, "KeyError: $split_k")
            end
        end

        # if no overlaping of polyA coordinates with GFF3 annotated coordinates - NA is added.
        if ct == 0
            ft = "NA"
            gn = "NA"
            bt = "NA"

            push!(dfna, [seqname(i1) first(i1)-1 last(i1) gn c str ft me mi mx bt])
        end
    end

    sort!(dfn, [:Chr, :Start, :End, :Name], rev=true)
    sort!(dfp, [:Chr, :Start, :End, :Name])
    sort!(dfna, [:Chr, :Start, :End, :Name])

    # remove duplicated gene features
    dfp = rmdups(dfp)
    dfn = rmdups(dfn)

    enumeratenames!(dfn)
    enumeratenames!(dfp)
    enumeratenames!(dfna)

    df = vcat(dfp, dfn, dfna)

    return df
end

function annotate_clusters(a::DataFrame, b::Dict{String, IntervalCollection{String}};
                           verbose::Bool=false)::DataFrame

    # Two dataframes for '+' and '-' strands.
    dfp = DataFrame(Chr=String[], Start=Int64[], End=Int64[], GeneName=String[],
                    ClusterSize=Int32[], Strand=String[], Feature=String[],
                    ClusterCenter=Int64[], Biotype=String[],
                    ClusterMedian=Float64[], ClusterMean=Float64[],
                    ClusterMin=Int64[], ClusterMax=Int64[],
                    Cluster1stQuartile=Float64[], Cluster3rdQuartile=Float64[],
                    TSMedian=Float64[], TSMean=Float64[],
                    TSMin=Int64[], TSMax=Int64[],
                    TS1stQuartile=Float64[], TS3rdQuartile=Float64[])

    dfn = DataFrame(Chr=String[], Start=Int64[], End=Int64[], GeneName=String[],
                    ClusterSize=Int32[], Strand=String[], Feature=String[],
                    ClusterCenter=Int64[], Biotype=String[],
                    ClusterMedian=Float64[], ClusterMean=Float64[],
                    ClusterMin=Int64[], ClusterMax=Int64[],
                    Cluster1stQuartile=Float64[], Cluster3rdQuartile=Float64[],
                    TSMedian=Float64[], TSMean=Float64[],
                    TSMin=Int64[], TSMax=Int64[],
                    TS1stQuartile=Float64[], TS3rdQuartile=Float64[])

    # Dataframe for not annotated sites.
    dfna = DataFrame(Chr=String[], Start=Int64[], End=Int64[], GeneName=String[],
                    ClusterSize=Int32[], Strand=String[], Feature=String[],
                    ClusterCenter=Int64[], Biotype=String[],
                    ClusterMedian=Float64[], ClusterMean=Float64[],
                    ClusterMin=Int64[], ClusterMax=Int64[],
                    Cluster1stQuartile=Float64[], Cluster3rdQuartile=Float64[],
                    TSMedian=Float64[], TSMean=Float64[],
                    TSMin=Int64[], TSMax=Int64[],
                    TS1stQuartile=Float64[], TS3rdQuartile=Float64[])
    multiple_features = Int64(0)

    for i1 in eachrow(a)
        ct = Int64(0)
        # Parse metadata from polyA site.
        chr = i1[1]
        it1 = i1[2]
        it2 = i1[3]
        str = i1[4]
        c = i1[6]
        center = i1[5]
        cluster_multiple_features = Int64(0)

        split_keys = get_split_key(chr, it1, it2)

        for split_k in split_keys
            try
                for i2 in b[split_k]
                    if chechoverlap(Interval(chr, it1, it2, str[1]), i2)
                        cluster_multiple_features += 1
                        # parse metadata from GFF3 anotation
                        splt::Array{SubString{String},1} = split(metadata(i2),";")
                        ft = split(splt[7],"=")[2]
                        gn = split(splt[6],"=")[2]
                        bt = split(splt[5],"=")[2]
                        ct += 1

                        if str == "-"
                            push!(dfn, hcat([chr it1-1 it2 gn c str ft center bt],
                                            [i1[7] i1[8] i1[9] i1[10] i1[11]],
                                            [i1[12] i1[13] i1[14] i1[15] i1[16]],
                                            [i1[17] i1[18]]))

                        else
                            push!(dfp, hcat([chr it1-1 it2 gn c str ft center bt],
                                            [i1[7] i1[8] i1[9] i1[10] i1[11]],
                                            [i1[12] i1[13] i1[14] i1[15] i1[16]],
                                            [i1[17] i1[18]]))
                        end
                    end
                end

            catch
                KeyError(split_k)
                if verbose
                    println(STDERR, "KeyError: $split_k")
                end
            end
        end

        if cluster_multiple_features > 1
            multiple_features += 1
        end

        # if no overlaping of polyA coordinates with GFF3 annotated coordinates - NA is added.
        if ct == 0
            ft = "NA"
            gn = "NA"
            bt = "NA"

            push!(dfna, hcat([chr it1-1 it2 gn c str ft center bt],
                            [i1[7] i1[8] i1[9] i1[10] i1[11]],
                            [i1[12] i1[13] i1[14] i1[15] i1[16]],
                            [i1[17] i1[18]]))
        end
    end

    sort!(dfn, [:Chr, :Start, :End, :GeneName], rev=true)
    sort!(dfp, [:Chr, :Start, :End, :GeneName])
    sort!(dfna, [:Chr, :Start, :End, :GeneName])

    # remove duplicated gene features
    dfp = rmdups(dfp; cluster=true, verbose=verbose)
    dfn = rmdups(dfn; cluster=true, verbose=verbose)

    enumeratenames!(dfn)
    enumeratenames!(dfp)
    enumeratenames!(dfna)

    df = vcat(dfp, dfn, dfna)

    return df
end

"""
    chechoverlap(i1::Interval, i2::Interval)

    Checks if intervals overlap.
"""
function chechoverlap(i1::Interval, i2::Interval)::Bool

    # only if chr match and strand match
    if cmp(seqname(i1), seqname(i2)) == 0 && strand(i1) == strand(i2)
        if i2.first <= i1.last && i2.first >= i1.first
            return true
        elseif i2.last <= i1.last && i2.last >= i1.first
            return true
        elseif i1.first <= i2.last && i1.first >= i2.first
            return true
        elseif i1.last >= i2.first && i1.last <= i2.last
            return true
        else
            return false
        end
    else
        return false
    end

end

"""
    enumeratenames!(df::DataFrame)

    Enumerates gene names in the data frame by adding incrementing number.
"""
function enumeratenames!(df::DataFrame)::DataFrame

    ct = Int64(0)
    l = size(df)[1]
    na_ct = Int64(0)

    for (i, x) in enumerate(eachrow(df))
        ct += 1

        if x[4] == "NA" && i+1 <= l
            na_ct += 1
            x[4] = x[4]*".$na_ct"

        elseif i+1 <= l
            if x[4] == eachrow(df)[i+1][4]
                x[4] = x[4]*".$ct"
            else
                x[4] = x[4]*".$ct"
                ct = 0
            end

        else
            x[4] = x[4]*".$ct"
            ct = 0
        end
    end

    return df
end

"""
    rmdups(dframe::DataFrame)

    Removes duplicated entries differentiating by feature.
    DataFrame should be sorted by Chr, Start, End and gene name.
    Featues are rated by their annotation hierarchy.
    Keeps most accurate feature (gene has exons and introns, etc.).

    # Arguments
    - `dframe::DataFrame`: a sorted dataframe.
"""
function rmdups(dframe::DataFrame; cluster::Bool=false, verbose::Bool=false)::DataFrame

    if cluster
        df = DataFrame(Chr=String[], Start=Int64[], End=Int64[], GeneName=String[],
                        ClusterSize=Int32[], Strand=String[], Feature=String[],
                        ClusterCenter=Int64[], Biotype=String[],
                        ClusterMedian=Float64[], ClusterMean=Float64[],
                        ClusterMin=Int64[], ClusterMax=Int64[],
                        Cluster1stQuartile=Float64[], Cluster3rdQuartile=Float64[],
                        TSMedian=Float64[], TSMean=Float64[],
                        TSMin=Int64[], TSMax=Int64[],
                        TS1stQuartile=Float64[], TS3rdQuartile=Float64[])

    else
        df = DataFrame(Chr=String[], Start=Int64[], End=Int64[], Name=String[],
                       Counts=Int32[], Strand=String[], Feature=String[],
                       Median=Float32[], Min=Int16[], Max=Int16[], Biotype=String[])
    end

    d = Dict("gene"=>13, "transcript"=>12, "mRNA"=>12, "exon"=>11, "intron"=>10, "CDS"=>9,
             "UTR"=>8, "three_prime_UTR"=>7, "five_prime_UTR"=>6,
             "polyA_sequence"=>5, "polyA_site"=>4, "start_codon"=>3,
             "stop_codon"=>2, "Selenocysteine"=>1)

    ct = Int64(0)
    ct2 = Int64(0)

    for y in eachrow(dframe[:,1:6])
        ct2 += 1

        if ct == 0
            if cluster
                push!(df, hcat([dframe[1][ct2] dframe[2][ct2] dframe[3][ct2]],
                               [dframe[4][ct2] dframe[5][ct2] dframe[6][ct2]],
                               [dframe[7][ct2] dframe[8][ct2] dframe[9][ct2]],
                               [dframe[10][ct2] dframe[11][ct2]],
                               [dframe[12][ct2] dframe[13][ct2]],
                               [dframe[14][ct2] dframe[15][ct2]],
                               [dframe[16][ct2] dframe[17][ct2]],
                               [dframe[18][ct2] dframe[19][ct2]],
                               [dframe[20][ct2] dframe[21][ct2]]))

            else
                push!(df, hcat([dframe[1][ct2] dframe[2][ct2] dframe[3][ct2]],
                               [dframe[4][ct2] dframe[5][ct2] dframe[6][ct2]],
                               [dframe[7][ct2] dframe[8][ct2] dframe[9][ct2]],
                               [dframe[10][ct2] dframe[11][ct2]]))
            end
            ct += 1
        end

        if ct > 0
            if (df[end,:][1] == [y[1]] && df[end,:][2] == [y[2]]
                && df[end,:][3] == [y[3]] && df[end,:][4] == [y[4]]
                && df[end,:][5] == [y[5]] && df[end,:][6] == [y[6]])

                try
                    if d[dframe[7][ct2]] < d[df[end,:][7][1]]
                        df[7][ct] = dframe[7][ct2]
                    end

                catch
                    KeyError(dframe[7][ct2])
                    if verbose
                        println(STDERR, "KeyError: $ct2 ")
                    end
                end

            else
                ct += 1

                if cluster
                    push!(df, hcat([dframe[1][ct2] dframe[2][ct2] dframe[3][ct2]],
                                   [dframe[4][ct2] dframe[5][ct2] dframe[6][ct2]],
                                   [dframe[7][ct2] dframe[8][ct2] dframe[9][ct2]],
                                   [dframe[10][ct2] dframe[11][ct2]],
                                   [dframe[12][ct2] dframe[13][ct2]],
                                   [dframe[14][ct2] dframe[15][ct2]],
                                   [dframe[16][ct2] dframe[17][ct2]],
                                   [dframe[18][ct2] dframe[19][ct2]],
                                   [dframe[20][ct2] dframe[21][ct2]]))

                else
                    push!(df, hcat([dframe[1][ct2] dframe[2][ct2] dframe[3][ct2]],
                                   [dframe[4][ct2] dframe[5][ct2] dframe[6][ct2]],
                                   [dframe[7][ct2] dframe[8][ct2] dframe[9][ct2]],
                                   [dframe[10][ct2] dframe[11][ct2]]))

                end

            end
        end
    end

    return df
end

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
pprint(matrix) = for i = 1:size(matrix, 1) println(join(matrix[i, :])) end
function printmaze(maze)
    walls = split("╹ ╸ ┛ ╺ ┗ ━ ┻ ╻ ┃ ┓ ┫ ┏ ┣ ┳ ╋")
    h, w = size(maze)
    f = cell -> 2 ^ ((3cell[1] + cell[2] + 3) / 2)
    wall(i, j) = if maze[i,j] == 0 " " else
        walls[Int(sum(f, filter(x -> maze[x...] != 0, neighbors([i, j], [h, w], 1)) .- [[i, j]]))]
    end
    mazewalls = collect(wall(i, j) for i in 1:2:h, j in 1:w)
    pprint(mazewalls)
end

printmaze(maze(10, 10))

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module CheqFilters

using ..GPUUtils, ..Utilities, JuAFEM
using CuArrays, ..FEA, Statistics
using GPUArrays: GPUVector
import ..GPUUtils: whichdevice
using Parameters: @unpack
import CUDAdrv

export  AbstractCheqFilter,
        CheqFilter

const dev = CUDAdrv.device()
const ctx = CUDAdrv.CuContext(dev)

abstract type AbstractCheqFilter end

@params struct FilterMetadata
    cell_neighbouring_nodes
    cell_node_weights
end
@define_cu(FilterMetadata, :cell_neighbouring_nodes, :cell_node_weights)
GPUUtils.whichdevice(m::FilterMetadata) = whichdevice(m.cell_neighbouring_nodes)

@params struct CheqFilter{_filtering, T, TV <: AbstractVector{T}} <: AbstractCheqFilter
    filtering::Val{_filtering}
    metadata::FilterMetadata
    rmin::T
    nodal_grad::TV
    last_grad::TV
    cell_weights::TV
end
@define_cu(CheqFilter, :metadata, :nodal_grad, :last_grad, :cell_weights)
GPUUtils.whichdevice(c::CheqFilter) = whichdevice(c.nodal_grad)

function FilterMetadata(::Type{T}, ::Type{TI}) where {T, TI}
    cell_neighbouring_nodes = Vector{TI}[]
    cell_node_weights = Vector{T}[]

    return FilterMetadata(cell_neighbouring_nodes, cell_node_weights)
end

function FilterMetadata(solver, rmin::T, ::Type{TI}) where {T, TI}
    problem = solver.problem
    cell_neighbouring_nodes, cell_node_weights = get_neighbour_info(problem, rmin)
    return FilterMetadata(RaggedArray(cell_neighbouring_nodes), RaggedArray(cell_node_weights))
end

function get_neighbour_info(problem, rmin::T) where {T}
    @unpack black, white, varind = problem
    dh = problem.ch.dh
    grid = dh.grid
    @unpack node_cells = problem.metadata
    TI = eltype(node_cells.offsets)

    all_neighbouring_nodes = Vector{TI}[]
    all_node_weights = Vector{T}[]
    sizehint!(all_neighbouring_nodes, getncells(grid))
    sizehint!(all_node_weights, getncells(grid))

    cells_to_traverse = TI[]
    neighbouring_nodes = TI[]
    node_weights = T[]
    visited_cells = Set{Int}()
    for cell in CellIterator(dh)
        current_cell_id = cell.current_cellid.x
        if black[current_cell_id] || white[current_cell_id]
            continue
        end
        empty!(cells_to_traverse)
        empty!(visited_cells)
        empty!(neighbouring_nodes)
        empty!(node_weights)

        center = Node(mean(cell.coords).data)

        push!(cells_to_traverse, current_cell_id)
        push!(visited_cells, current_cell_id)
        while !isempty(cells_to_traverse)
            # Takes first element and removes it -> breadth first traversal
            cell_id = popfirst!(cells_to_traverse)
            for n in grid.cells[cell_id].nodes
                node = getnodes(grid, n)
                dist = norm(node.x - center.x)
                if dist < rmin
                    push!(neighbouring_nodes, n)
                    push!(node_weights, max(rmin-dist, zero(T)))
                    n_cells = node_cells[n]
                    for c in n_cells
                        next_cell_id = c[1]
                        if !(next_cell_id in visited_cells)
                            push!(cells_to_traverse, next_cell_id)
                            push!(visited_cells, next_cell_id)
                        end
                    end
                end
            end
        end
        push!(all_neighbouring_nodes, copy(neighbouring_nodes))
        push!(all_node_weights, copy(node_weights))
    end

    return all_neighbouring_nodes, all_node_weights
end

CheqFilter{true}(args...) = CheqFilter(Val(true), args...)
CheqFilter{false}(args...) = CheqFilter(Val(false), args...)

function CheqFilter(::Val{filtering}, solver::AbstractFEASolver, args...) where {filtering}
    CheqFilter(Val(filtering), whichdevice(solver), solver, args...)
end

function CheqFilter(::Val{true}, ::CPU, solver::TS, rmin::T, ::Type{TI}=Int) where {T, TI<:Integer, TS<:AbstractFEASolver}
    metadata = FilterMetadata(solver, rmin, TI)
    TM = typeof(metadata)
    problem = solver.problem
    grid = problem.ch.dh.grid
    nnodes = getnnodes(grid)
    nodal_grad = zeros(T, nnodes)
    TV = typeof(nodal_grad)

    black = problem.black
    white = problem.white
    nel = length(black)
    nfc = sum(black) + sum(white)
    last_grad = zeros(T, nel-nfc)

    cell_weights = zeros(T, nnodes)

    return CheqFilter(Val(true), metadata, rmin, nodal_grad, last_grad, cell_weights)
end

function CheqFilter(::Val{false}, ::CPU, solver::TS, rmin::T, ::Type{TI}=Int) where {T, TS<:AbstractFEASolver, TI<:Integer}
    metadata = FilterMetadata(T, TI)
    nodal_grad = T[]
    last_grad = T[]
    cell_weights = T[]
    return CheqFilter(Val(false), metadata, rmin, nodal_grad, last_grad, cell_weights)
end

function (cf::CheqFilter{true, T})(grad, vars, elementinfo) where {T}
    cf.rmin <= 0 && return
    @unpack nodal_grad, cell_weights, metadata = cf
    @unpack black, white, varind, cellvolumes, cells = elementinfo
    @unpack cell_neighbouring_nodes, cell_node_weights = metadata
    node_cells = elementinfo.metadata.node_cells

    update_nodal_grad!(nodal_grad, node_cells, cell_weights, cells, cellvolumes, black, white, varind, grad)
    normalize_grad!(nodal_grad, cell_weights)
    update_grad!(grad, black, white, varind, cell_neighbouring_nodes, cell_node_weights, nodal_grad)
end

(cf::CheqFilter{false})(args...) = nothing

function update_nodal_grad!(nodal_grad::AbstractVector, node_cells, cell_weights, cells, cellvolumes, black, white, varind, grad)
    T = eltype(nodal_grad)
    for n in 1:length(nodal_grad)
        nodal_grad[n] = zero(T)
        cell_weights[n] = zero(T)
        r = node_cells.offsets[n]:node_cells.offsets[n+1]-1
        for i in r
            c = node_cells.values[i][1]
            if black[c] || white[c]
                continue
            end
            ind = varind[c]
            w = cellvolumes[c]
            cell_weights[n] += w
            nodal_grad[n] += w * grad[ind]
        end
    end
    return
    #=
function update_nodal_grad!(nodal_grad::AbstractVector, cell_weights, cells, cellvolumes, black, white, varind, grad)
    T = eltype(nodal_grad)
    nodal_grad .= zero(T)
    cell_weights .= 0
    @inbounds for i in 1:length(cells)
        if black[i] || white[i]
            continue
        end
        nodes = cells[i].nodes
        for n in nodes
            ind = varind[i]
            w = cellvolumes[i]
            cell_weights[n] += w
            nodal_grad[n] += w*grad[ind]
        end
    end
    =#
end

function update_nodal_grad!(nodal_grad::GPUVector, node_cells, args...)
    T = eltype(nodal_grad)
    allargs = (nodal_grad, node_cells.offsets, node_cells.values, args...)
    callkernel(dev, cheq_kernel1, allargs)
    CUDAdrv.synchronize(ctx)
    return
end

function cheq_kernel1(nodal_grad, node_cells_offsets, node_cells_values, cell_weights,
        cells, cellvolumes, black, white, varind, grad)
    T = eltype(nodal_grad)
    n = @thread_global_index()
    offset = @total_threads()
    while n <= length(nodal_grad)
        nodal_grad[n] = zero(T)
        cell_weights[n] = zero(T)
        r = node_cells_offsets[n]:node_cells_offsets[n+1]-1
        for i in r
            c = node_cells_values[i][1]
            if black[c] || white[c]
                continue
            end
            ind = varind[c]
            w = cellvolumes[c]
            cell_weights[n] += w
            nodal_grad[n] += w * grad[ind]
        end
        n += offset
    end
end

function normalize_grad!(nodal_grad::AbstractVector, cell_weights)
    for n in 1:length(nodal_grad)
        if cell_weights[n] > 0
            nodal_grad[n] /= cell_weights[n]
        end
    end
end
function normalize_grad!(nodal_grad::GPUVector, cell_weights)
    T = eltype(nodal_grad)
    args = (nodal_grad, cell_weights)
    callkernel(dev, cheq_kernel2, args)
    CUDAdrv.synchronize(ctx)
    return
end
function cheq_kernel2(nodal_grad, cell_weights)
    T = eltype(nodal_grad)
    n = @thread_global_index()
    offset = @total_threads()
    while n <= length(nodal_grad)
        w = cell_weights[n]
        w = ifelse(w > 0, w, one(T))
        nodal_grad[n] /= w
        n += offset
    end
    return
end

function update_grad!(grad::AbstractVector, black, white, varind, cell_neighbouring_nodes, cell_node_weights, nodal_grad)
    @inbounds for i in 1:length(black)
        if black[i] || white[i]
            continue
        end
        ind = varind[i]
        nodes = cell_neighbouring_nodes[ind]
        if length(nodes) == 0
            continue
        end
        weights = cell_node_weights[ind]
        grad[ind] = dot(view(nodal_grad, nodes), weights) / sum(weights)
    end

    return
end

function update_grad!(grad::GPUVector, black, white, varind, cell_neighbouring_nodes, cell_node_weights, nodal_grad)
    T = eltype(grad)
    allargs = (grad, black, white, varind, cell_neighbouring_nodes.offsets, cell_neighbouring_nodes.values, cell_node_weights.values, nodal_grad)
    callkernel(dev, cheq_kernel3, allargs)
    CUDAdrv.synchronize(ctx)
    return
end
function cheq_kernel3(grad, black, white, varind, cell_neighbouring_nodes_offsets, cell_neighbouring_nodes_values, cell_node_weights_values, nodal_grad)
    T = eltype(nodal_grad)
    i = @thread_global_index()
    offset = @total_threads()
    while i <= length(black)
        if black[i] || white[i]
            continue
        end
        ind = varind[i]
        r = cell_neighbouring_nodes_offsets[ind]:cell_neighbouring_nodes_offsets[ind+1]-1
        length(r) == 0 && continue
        grad[ind] = zero(T)
        sum_weights = zero(T)
        for linear_ind in r
            w = cell_node_weights_values[linear_ind]
            sum_weights += w
            node_ind = cell_neighbouring_nodes_values[linear_ind]
            grad[ind] += nodal_grad[node_ind] * w
        end
        grad[ind] /= sum_weights
        i += offset
    end

    return
end

end


text:
import MPI
MPI.Init()

comm = MPI.COMM_WORLD
N = 5
root = 0

if MPI.Comm_rank(comm) == root
    print(" Running on $(MPI.Comm_size(comm)) processes\n")
end
MPI.Barrier(comm)

if MPI.Comm_rank(comm) == root
    A = [i*(1.0 + im*2.0) for i = 1:N]
else
    A = Array{ComplexF64}(undef, N)
end

MPI.Bcast!(A, root, comm)

print("rank = $(MPI.Comm_rank(comm)), A = $A\n")

if MPI.Comm_rank(comm) == root
    B = Dict("foo" => "bar")
else
    B = nothing
end

B = MPI.bcast(B, root, comm)
print("rank = $(MPI.Comm_rank(comm)), B = $B\n")

if MPI.Comm_rank(comm) == root
    f = x -> x^2 + 2x - 1
else
    f = nothing
end

f = MPI.bcast(f, root, comm)
print("rank = $(MPI.Comm_rank(comm)), f(3) = $(f(3))\n")


text:
using QuakeML, Test
using Dates: DateTime

@testset "Field access" begin
    @testset "CreationInfo" begin
        let ci = QuakeML.CreationInfo()
            @test ci.agency_id === missing
            ci.agency_id = "a"^64
            @test ci.agency_id == "a"^64
            @test_throws ArgumentError ci.agency_id = "a"^65
            @test_throws ArgumentError ci.version = "a"^65
            @test_throws ArgumentError ci.author = "a"^129
        end
    end

    @testset "WaveformStreamID" begin
        let ws = QuakeML.WaveformStreamID(network_code="AN", station_code="QML")
            @test ws.network_code == "AN"
            @test ws.station_code == "QML"
            @test ws.uri === missing
            @test_throws ArgumentError ws.uri = "not a URI"
            ws.uri = "smi:local/uri"
            @test ws.uri == QuakeML.ResourceReference("smi:local/uri")
            @test_throws ArgumentError ws.network_code = "123456789"
            @test_throws ArgumentError ws.station_code = "123456789"
            @test_throws ArgumentError ws.location_code = "123456789"
            @test_throws ArgumentError ws.channel_code = "123456789"
        end
    end

    @testset "Amplitude" begin
        let a = QuakeML.Amplitude(generic_amplitude=1)
            a.snr = 1
            @test a.snr === 1.0
            a.snr = 1.f0
            @test a.snr === 1.0
            @test_throws ArgumentError a.type = "a"^33
            @test_throws ArgumentError a.magnitude_hint = "a"^33
        end
    end

    @testset "Magnitude" begin
        let a = QuakeML.Magnitude(mag=1)
            a.mag == QuakeML.RealQuantity(value=1.0)
            @test_throws ArgumentError a.type = "a"^33
        end
    end

    @testset "StationMagnitude" begin
        let a = QuakeML.StationMagnitude(mag=1)
            @test a.mag == QuakeML.RealQuantity(value=1)
            @test_throws ArgumentError a.type = "a"^33
        end
    end

    @testset "Origin" begin
        let a = QuakeML.Origin(time=DateTime(2000), longitude=0, latitude=1)
            @test a.time == QuakeML.TimeQuantity(value=DateTime(2000))
            @test a.longitude == QuakeML.RealQuantity(value=0.0)
            @test a.latitude == QuakeML.RealQuantity(value=1.0)
            @test_throws ArgumentError a.region = "a"^129
        end
    end

    @testset "OriginQuality" begin
        let oq = QuakeML.OriginQuality()
            @test_throws ArgumentError oq.ground_truth_level = "a"^33
            oq.ground_truth_level = "a"^32
            @test oq.ground_truth_level == "a"^32
        end
    end
end

max_stars_repo_name: anowacki/QuakeML


text:
using AutoAligns
using Base.Test

import AutoAligns: align_string, get_alignment

## test alignment building blocks
@test align_string(left, "foo", 8) == "foo     "
@test align_string(right, "foo", 8) == "     foo"
@test align_string(center, "foo", 8) == "   foo  "

## test default alignment determination
@test get_alignment(left, 0) == left
@test get_alignment(right, 0) == right
@test get_alignment(center, 0) == center
@test_throws MethodError get_alignment(119, 0)
let v = [left, right, center]
    @test get_alignment(v, 1) == left
    @test get_alignment(v, 2) == right
    @test get_alignment(v, 3) == center
    @test get_alignment(v, 4) == center # default
end
let d = Dict(1 => left, 2 => right, 3 => center, :default => center)
    @test get_alignment(d, 1) == left
    @test get_alignment(d, 2) == right
    @test get_alignment(d, 3) == center
    @test get_alignment(d, 4) == center # default
end

## example 1
let s = """
1   2     3
foo bar baz
 5   55 555
"""
    aa = AutoAlign(align = Dict(:default => left, 5 => right))
    println(aa, 1, " ", 2, " ", 3)
    println(aa, "foo", " ", "bar", " ", "baz")
    aa.align = center
    println(aa, 5, " ", 55, " ", 555)
    @test string(aa) == s
end

## example 2 (ragged, no final newline, default from Dict)
let s = """
     1 2   3
ragged
0      44 55"""
    aa = AutoAlign(align = Dict(:default => right, 3 => left))
    print(aa, 1)
    print(aa, " "; align=left)
    println(aa, 2, " ", 3)
    println(aa, "ragged")
    print(aa, 0; align=left)
    print(aa, " ", 44, " ", 55)
    @test string(aa) == s
end


text:
#!/usr/bin/env julia

"""
gff_feature.jl - Houses functions that process GFF features.

By: Shaun Adkins (sadkins@som.umaryland.edu)
"""

isstranded(strand_type::String) = return (strand_type == "no" ? false : true)

function create_uniq_coords_dict(features::Array{GFF3.Record,1}, stranded_type::String)
    uniq_coords = Dict{String, Dict}()
    seqids = Set(map(x -> GFF3.seqid(x), features))
    for seqid in seqids
        uniq_coords[seqid] = get_nonoverlapping_coords_by_seqid(features, seqid)
        # If alignment data is unstranded, then symdiff both strands to + strand
        if !isstranded(stranded_type)
            # Ensure a default empty set is used if that strand has no coordinates
            symdiff!(get!(uniq_coords[seqid],'+', BitSet()), pop!(uniq_coords[seqid], '-', BitSet()))
        end
    end
    return uniq_coords
end

function get_feature_coords_set(feature::GFF3.Record)
    """Get the range of coordinates for this feature, returned as a Set."""
    return BitSet(GFF3.leftposition(feature) : GFF3.rightposition(feature))
end

function get_featurename_from_attrs(feature::GFF3.Record, attr_type::String)
    """Get attribute ID to use as the feature name."""
    gene_vector = GFF3.attributes(feature, attr_type)    # col 9
    validate_feature_attribute(gene_vector)
    return gene_vector[1]
end

function get_feature_nonoverlapping_coords(feature::GFF3.Record, uniq_coords::Dict{String, Dict}, stranded::Bool)
    """Get set of nonoverlapping coordinates for given feature."""
    seqid = GFF3.seqid(feature)
    strand = getstrand(feature, stranded)
    feat_coords = get_feature_coords_set(feature)
    return intersect(feat_coords, uniq_coords[seqid][strand])
end

function get_nonoverlapping_coords_by_seqid(features::Array{GFF3.Record,1}, seqid::String)
    """ Get dictionary nonoverlapping coordinates based on the features on this sequence ID."""
    uniq_seqid_coords = Dict{Char,BitSet}()
    # Get only features with this reference id
    seqid_feats = filter(x -> GFF3.seqid(x) == seqid, features)
    strand = Set(map(x -> GFF3.strand(x), seqid_feats))
    for s in strand
        strandchar = convert(Char, s)
        uniq_seqid_coords[strandchar] = get_nonoverlapping_coords_by_seqid_and_strand(seqid_feats, s)
    end
    return uniq_seqid_coords
end

function get_nonoverlapping_coords_by_seqid_and_strand(features::Array{GFF3.Record,1}, strand::Strand)
    """ Get set of nonoverlapping coordinates based on the features for this strand of this sequence ID."""
    # Get features one strand at a time
    seqid_feats_by_strand = filter(x-> GFF3.strand(x) == strand, features)
    return BitSet(mapreduce(x -> get_feature_coords_set(x), symdiff, seqid_feats_by_strand))
end

function getstrand(feature::GFF3.Record, stranded::Bool)
    """Get strand of GFF3 feature with respect to strandedness arguments."""
    return getstrand(convert(Interval, feature), stranded)
end

function validate_feature_attribute(gene_vector::Vector{String})
    """Ensure only one attribute of this type exists for this feature record."""
    length(gene_vector) > 0 || error("ERROR - Attribute field 'ID' found to have no entries.\n")
    length(gene_vector) == 1 || error("ERROR - Attribute field 'ID' found to have multiple entries.\n")
end

max_stars_repo_name: IGS/FADU

max_issues_repo_name: adkinsrs/FADU

max_forks_repo_name: adkinsrs/FADU


text:
module Clapeyron
using StaticArrays
using LinearAlgebra
using SparseArrays
#for the assoc solver and the sparse packed VofV
import PackedVectorsOfVectors
const PackedVofV = PackedVectorsOfVectors.PackedVectorOfVectors

#for non allocating vectors of zeros and ones
using FillArrays: FillArrays

using Roots: Roots
using NLSolvers
import Metaheuristics
using DiffResults, ForwardDiff
using Scratch

include("solvers/Solvers.jl")
using .Solvers
using .Solvers: log, sqrt
∂Tag = Solvers.∂Tag
using Unitful
import LogExpFunctions

#The Base of Clapeyron: EoSModel and eos(model,V,T,z)
include("base/constants.jl")
include("base/EoSModel.jl")
include("base/eosshow.jl")

#EoSParam, ClapeyronParam
include("params/paramvectors.jl")
include("params/ClapeyronParam.jl")
include("params/combiningrules.jl")

using CSV, Tables
#getparams machinery
include("database/ParamOptions.jl")
include("database/database.jl")
include("database/UserReader.jl")

#macros, used for defining models
include("utils/macros.jl")

#splitting models, useful for methods.
include("utils/split_model.jl")

#Clapeyron methods (AD, property solvers, etc)
include("methods/methods.jl")

#=
the dependency chain is the following:

base --> params -|-> database ----------------|
                 |-> split_model --> methods -|-> models
                 |-> macros ------------------|

=#

#Clapeyron EoS collection
include("models/types.jl") #type hierarchy
include("models/ideal/ideal.jl")
include("models/ideal/BasicIdeal.jl")
include("models/ideal/MonomerIdeal.jl")
include("models/ideal/ReidIdeal.jl")
include("models/ideal/WalkerIdeal.jl")
include("models/ideal/JobackIdeal.jl")

include("models/SAFT/PCSAFT/PCSAFT.jl")
include("models/SAFT/PCSAFT/variants/sPCSAFT.jl")
include("models/SAFT/ogSAFT/ogSAFT.jl")
include("models/SAFT/CPA/CPA.jl")
include("models/SAFT/CPA/variants/sCPA.jl")
include("models/SAFT/SAFTVRSW/SAFTVRSW.jl")
include("models/SAFT/LJSAFT/LJSAFT.jl")
include("models/SAFT/softSAFT/softSAFT.jl")
include("models/SAFT/SAFTVRMie/SAFTVRMie.jl")
include("models/SAFT/SAFTVRMie/variants/SAFTVRQMie.jl")
include("models/SAFT/SAFTgammaMie/SAFTgammaMie.jl")
include("models/SAFT/CKSAFT/CKSAFT.jl")
include("models/SAFT/CKSAFT/variants/sCKSAFT.jl")
include("models/SAFT/BACKSAFT/BACKSAFT.jl")
include("models/SAFT/equations.jl")

include("models/cubic/equations.jl")
include("models/cubic/vdW.jl")
include("models/cubic/RK/RK.jl")
include("models/cubic/PR/PR.jl")

include("models/Activity/Wilson/Wilson.jl")
include("models/Activity/NRTL/NRTL.jl")
include("models/Activity/UNIQUAC/UNIQUAC.jl")
include("models/Activity/UNIFAC/UNIFAC.jl")
include("models/Activity/UNIFAC/variants/ogUNIFAC.jl")
include("models/Activity/UNIFAC/variants/PSRK.jl")
include("models/Activity/UNIFAC/variants/VTPR.jl")
include("models/Activity/equations.jl")

include("models/Activity/COSMOSAC/utils.jl")
include("models/Activity/COSMOSAC/COSMOSAC02.jl")
include("models/Activity/COSMOSAC/COSMOSAC10.jl")
include("models/Activity/COSMOSAC/COSMOSACdsp.jl")

include("models/cubic/alphas/alphas.jl")
include("models/cubic/mixing/mixing.jl")
include("models/cubic/translation/translation.jl")

include("models/cubic/RK/variants/SRK.jl")
include("models/cubic/RK/variants/PSRK.jl")
include("models/cubic/PR/variants/PR78.jl")
include("models/cubic/PR/variants/VTPR.jl")
include("models/cubic/PR/variants/UMRPR.jl")

include("models/EmpiricHelmholtz/IAPWS95.jl")
include("models/EmpiricHelmholtz/PropaneRef.jl")
include("models/EmpiricHelmholtz/LJRef.jl")
include("models/EmpiricHelmholtz/MultiFluid/multifluid.jl")

include("models/LatticeFluid/SanchezLacombe/SanchezLacombe.jl")

include("models/SPUNG/SPUNG.jl")
include("models/UFTheory/UFTheory.jl")

include("models/cached/CachedEoS.jl")

end # module


text:
using RunTests
exit(run_tests())


text:
module InfinitesimalGenerators

using Arpack
using Distributions
using FillArrays
using FiniteDiff
using KrylovKit
using LinearAlgebra
using Roots

include("MarkovProcess.jl")
include("AdditiveFunctional.jl")
include("feynman_kac.jl")
include("principal_eigenvalue.jl")

export
MarkovProcess,
generator,
state_space,
stationary_distribution,
feynman_kac,
DiffusionProcess,
OrnsteinUhlenbeck,
CoxIngersollRoss,
AdditiveFunctional,
cgf,
tail_index,
AdditiveFunctionalDiffusion
end


text:
using WAV

struct Sound
  samples
  sr
end

"""
  Constructor for a `Sound` object
    Sound(fname)

  # Args
  * **fname** The filename of the wav file to load

  # Returns
  A `Sound` object containing the samples and the the sampling rate

  Examples
  =========

  s = Sound("sound.wav")
"""
function Sound(fname)
  s, sr = wavread(fname)
  s = vec(s)
  return Sound(s, sr)
end


text:
using UtilitiesForMRI, LinearAlgebra, CUDA, Test
CUDA.allowscalar(false)

# Cartesian domain
n = (256, 256, 256)
h = (abs(randn()), abs(randn()), abs(randn()))
X = spatial_sampling(Float64, n; h=h)

# Cartesian sampling in k-space
phase_encoding = (1,2)
K = kspace_Cartesian_sampling(X; phase_encoding=phase_encoding)

# Adjoint test (linear operator)
R = rotation()
nt, nk = size(K)
φ = randn(Float64, nt, 3)
Rφ = R(φ)
k = randn(Float64, nt, nk, 3)
l = randn(Float64, nt, nk, 3)
@test dot(Rφ*k, l) ≈ dot(k, Rφ'*l) rtol=1e-6

# Jacobian
RφK, ∂RφK = ∂(R()*K, φ)
@test RφK ≈ Rφ*coord(K) rtol=1e-6

# Adjoint test (Jacobian)
Δφ = randn(Float64, nt, 3); Δφ *= norm(φ)/norm(Δφ)
ΔR_ = ∂RφK*Δφ; ΔR = randn(ComplexF64, size(ΔR_)); ΔR *= norm(ΔR_)/norm(ΔR)
@test dot(∂RφK*Δφ, ΔR) ≈ dot(Δφ, ∂RφK'*ΔR) rtol=1e-6

# Gradient test
Δφ = randn(Float64, nt, 3); Δφ *= norm(φ)/norm(Δφ)
t = 1e-6
RφK_p1 = R(φ+0.5*t*Δφ)*coord(K)
RφK_m1 = R(φ-0.5*t*Δφ)*coord(K)
@test (RφK_p1-RφK_m1)/t ≈ ∂RφK*Δφ rtol=1e-5

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
export cwt, icwt, wts_dvv

# """
#
#   wtdtw(ref,cur,t,window,freqmin,freqmax)
#
# dv/v with dynamic time warping method from continuous wavelet transformation.
#
# This function uses the dynamic time warping method at each frequency from the
# continuous wavelet transform to compare the Reference waveform to the current
# waveform to get the relative seismic velocity variation (and associated error).
#
# # Arguments
# - `ref::AbstractArray`: Reference correlation.
# - `cur::AbstractArray`: Current correlation.
# - `t::AbstractArray`: time vector, common to both `ref` and `cur`.
# - `window::AbstractArray`: vector of the indices of the `cur` and `ref` windows
#                           on which you want to do the measurements
# - `freqmin::Float64`: minimum frequency in the correlation [Hz]
# - `freqmax::Float64`: maximum frequency in the correlation [Hz]
# - `maxlag::Int64`: number of maxlag id to search the distance.
# - `b::Int64`: value to control in distance calculation algorithm (see Mikesell et al. 2015).
# - `direction::Int64`: direction to accumulate errors (1=forward, -1=backward, 0=double to smooth)
# - `f0::Real`: Nondimensional frequency from Torrence & Campo, 1998 eq. 1 [Hz].
# - `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.
# - `standardize::Bool`: Remove mean and std from wavelet spectrum or not.
#
# # Returns
# - `dvv::AbstactArray`: Relative Velocity Change dv/v (in %)
# - `err::AbstractArray`: Errors in the dv/v measurements
#
# Originally written in python by Congcong Yuan (30 Jun, 2019)
# """
# function wtdtw(ref,cur,t,window,freqmin,freqmax;maxlag=80,b=1,direction=1,f0=6,dj=1/12,
#                standardize=true)
#
#     T = eltype(ref)
#     dt = mean(diff(t))
#     # apply cwt on two traces
#     cwt1,sj,freqs,coi = cwt(ref,dt,freqmin,freqmax, f0=f0,dj=dj)
#     cwt2,sj,freqs,coi = cwt(cur,dt,freqmin,freqmax, f0=f0,dj=dj)
#
#     # extract real values of cwt
#     rcwt1, rcwt2 = real.(cwt1), real.(cwt2)
#
#     # Use DTW method to extract dvv
#     Nfreq = length(freqs)
#     dvv = zeros(T,Nfreq)
#     err = zeros(T,Nfreq)
#
#     if standardize
#         standardize!(rcwt1)
#         standardize!(rcwt2)
#     end
#
#     for ii = 1:Nfreq
#
#         # find cone of influence
#         indcoi = findall(freqs[ii] .>= 1. ./coi)
#         indt = findall(x->x in indcoi, window)
#
#         if length(indt) == 0
#             continue
#         end
#
#         dvv[ii], err[ii] = dtw_dvv(rcwt1[window[indt],ii], rcwt2[window[indt],ii],t[window[indt]],maxlag, b, direction)
#     end
#     return dvv, err, freqs
# end

"""

  wts_dvv(ref,cur,t,window,freqmin,freqmax)

dv/v with stretching method from continuous wavelet transformation.

This function uses the stretching method at each frequency from the
continuous wavelet transform to compare the Reference waveform to the current
waveform to get the relative seismic velocity variation (and associated error).

# Arguments
- `ref::AbstractArray`: Reference correlation.
- `cur::AbstractArray`: Current correlation.
- `t::AbstractArray`: time vector, common to both `ref` and `cur`.
- `window::AbstractArray`: vector of the indices of the `cur` and `ref` windows
                          on which you want to do the measurements
- `freqmin::Float64`: minimum frequency in the correlation [Hz]
- `freqmax::Float64`: maximum frequency in the correlation [Hz]
- `f0::Real`: Nondimensional frequency from Torrence & Campo, 1998 eq. 1 [Hz].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.
- `standardize::Bool`: Remove mean and std from wavelet spectrum or not.

# Returns
- `dvv::AbstactArray`: Relative Velocity Change dv/v (in %)
- `err::AbstractArray`: Errors in the dv/v measurements

Originally written in python by Congcong Yuan (30 Jun, 2019)
"""
function wts_dvv(ref,cur,t,window,freqmin,freqmax;f0=6,dj=1/12,
               standardize=true)

    T = eltype(ref)
    dt = mean(diff(t))
    # apply cwt on two traces
    cwt1,sj,freqs,coi = cwt(ref,dt,freqmin,freqmax, f0=f0,dj=dj)
    cwt2,sj,freqs,coi = cwt(cur,dt,freqmin,freqmax, f0=f0,dj=dj)

    # Use DTW method to extract dvv
    Nfreq = length(freqs)
    dvv = zeros(T,Nfreq)
    err = zeros(T,Nfreq)

    for ii = 1:Nfreq
        icwt1 = icwt(cwt1[:,ii],sj[ii],dt)
        icwt2 = icwt(cwt2[:,ii],sj[ii],dt)

        if standardize
            standardize!(icwt1)
            standardize!(icwt2)
        end

        # find cone of influence
        indcoi = findall(freqs[ii] .>= 1. ./coi)
        indt = findall(x->x in indcoi, window)

        if length(indt) == 0
            continue
        end

        # calculate dvv with stretching
        dvv[ii], cc,cdp, ep, err[ii],allC = stretching(icwt1, icwt2,t,window[indt],
               freqmin,freqmax,ntrial=30,dvmin=-0.01,dvmax=0.01)
    end
    return dvv, err, freqs
end

# """
#
#   dtw_dvv(ref,cur,t,maxlag,d,direction)
#
# dv/v with dynamic time warping method.
#
# This function uses the dynamic time warping method to compare the optimal stretching
# between Reference waveform to the current waveform to get the relative seismic
# velocity variation (and associated error).
#
# # Arguments
# - `ref::AbstractArray`: Reference correlation.
# - `cur::AbstractArray`: Current correlation.
# - `t::AbstractArray`: time vector, common to both `ref` and `cur`.
# - `maxlag::Int64`: number of maxlag id to search the distance.
# - `b::Int64`: value to control in distance calculation algorithm (see Mikesell et al. 2015).
# - `direction::Int64`: direction to accumulate errors (1=forward, -1=backward, 0=double to smooth)
#
# # Returns
# - `dvv::Float64`: Relative Velocity Change dv/v (in %)
# - `err::Float64`: Errors in the dv/v measurements
#
# """
# function dtw_dvv(ref,cur,t,maxlag, b, direction)
#     dt = mean(diff(t))
#     stbarTime, stbar, dist, error = dtwdt(ref, cur, dt, maxLag=maxlag, b=b, direction=direction)
#
#     # perform linear regression
#     model = glm(@formula(Y ~0 + X),DataFrame(X=t,Y=stbarTime),Normal(),
#                 IdentityLink(),wts=ones(length(t)))
#
#     return coef(model)[1], stderror(model)[1]
# end

"""

  cwt(signal,dt,freqmin,freqmax)

Continuous wavelet transform of the signal at specified scales.
Note: uses the Morlet wavelet ONLY.

# Arguments
- `signal::AbstractArray`: Input signal array.
- `dt::AbstractFloat`: Sampling interval [s].
- `freqmin::AbstractFloat`: Minimum frequency for cwt [Hz].
- `freqmax::AbstractFloat`: Maximum frequency for cwt [Hz].
- `f0::Real`: Nondimensional frequency from Torrence & Campo, 1998 eq. 1 [Hz].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.

# Returns
- `W::AbstractArray`: Wavelet transform from Morlet wavelet.
- `freqs::AbstractArray`: Fourier frequencies for wavelet scales [Hz].
- `coi::AbstractArray`: Cone of influence - maximum period (in s) of useful
    information at that particular time. Periods greater than
    those are subject to edge effects.

This is a Julia translation of the cwt in pycwt https://github.com/regeirk/pycwt
"""
function cwt(signal::AbstractArray{T,1},dt::AbstractFloat,freqmin::AbstractFloat,
             freqmax::AbstractFloat;f0=6.,dj=1/12) where T <: AbstractFloat
    n0 = length(signal)
    flambda = T.((4 .* π) ./ (f0 .+ sqrt(2 .+ f0^2)))
    s0 = 2 * dt / flambda
    J = convert(Int,round(log2(n0 .* dt ./ s0) ./ dj))

    # The scales as of Mallat 1999
    sj = T.(s0 .* 2 .^ ((0:J) .* dj))
    freqs = 1 ./ (flambda .* sj)

    # subset by freqmin & freqmax
    ind = findall((freqs .> freqmin) .& (freqs .< freqmax))
    sj = sj[ind]
    freqs = freqs[ind]

    # signal fft
    signal_ft = fft(signal,1)
    N = length(signal_ft)
    # Fourier angular frequencies
    ftfreqs = T.(2 .* π * FFTW.fftfreq(n0,1/dt))

    # Creates wavelet transform matrix as outer product of scaled transformed
    # wavelets and transformed signal according to the convolution theorem.
    # (i)   Transform scales to column vector for outer product;
    # (ii)  Calculate 2D matrix [s, f] for each scale s and Fourier angular
    #       frequency f;
    # (iii) Calculate wavelet transform;
    psi_ft_bar = ((sj' .* ftfreqs[2] .* T(N)) .^ T(0.5)) .* conj.(psi_ft(ftfreqs * sj',f0))
    W = ifft(signal_ft .* psi_ft_bar,1)

    # Checks for NaN in transform results and removes them from the scales if
    # needed, frequencies and wavelet transform. Trims wavelet transform at
    # length `n0`.
    sel = findall(.!all(isnan.(W),dims=1)[:])
    if length(sel) > 0
        sj = sj[sel]
        freqs = freqs[sel]
        W = W[:,sel]
    end

    # Determines the cone-of-influence. Note that it is returned as a function
    # of time in Fourier periods. Uses triangualr Bartlett window with
    # non-zero end-points.
    coi = n0 ./ 2 .- abs.((0:n0-1) .- (n0 - 1.) ./2.)
    coi = T.(flambda ./ sqrt(2.) .* dt .* coi)
    return W, sj, freqs, coi
end

"""

  icwt(W,sj,dt)

Inverse continuous wavelet transform at specified scales.
Note: uses the Morlet wavelet ONLY.

# Arguments
- `W::AbstractArray`: Wavelet transform, the result of the `cwt` function.
- `sj::AbstractArray`: Scale indices as returned by the `cwt` function.
- `dt::AbstractFloat`: Sampling interval [s].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.

# Returns
- `iW::AbstractArray`: Inverse wavelet transform from Morlet wavelet.

This is a Julia translation of the icwt in pycwt https://github.com/regeirk/pycwt
Note that the pycwt version has incorrect scaling.
"""
function icwt(W::AbstractArray, sj::AbstractArray, dt::AbstractFloat;dj=1/12)
    T = real(eltype(W))
    # As of Torrence and Compo (1998), eq. (11)
    iW = T(dj .* sqrt(dt) ./ 0.776 .* (pi ^ 0.25)) .* sum(real.(W) ./ (sj .^ T(0.5))',dims=2)
    return iW
end

"""

  icwt(W,sj,dt)

Inverse continuous wavelet transform at specified scale.
Note: uses the Morlet wavelet ONLY.

# Arguments
- `W::AbstractArray`: Wavelet transform, the result of the `cwt` function.
- `sj::AbstractFloat`: Scale index as returned by the `cwt` function.
- `dt::AbstractFloat`: Sampling interval [s].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.

# Returns
- `iW::AbstractArray`: Inverse wavelet transform from Morlet wavelet.

This is a Julia translation of the icwt in pycwt https://github.com/regeirk/pycwt
Note that the pycwt version has incorrect scaling.
"""
function icwt(W::AbstractArray, sj::AbstractFloat, dt::AbstractFloat;dj=1/12)
    T = real(eltype(W))
    # As of Torrence and Compo (1998), eq. (11)
    iW = T(dj .* sqrt(dt) ./ 0.776 .* (pi ^ 0.25)) .* real.(W) ./ sj ^ T(0.5)
    return iW
end

function psi_ft(A::AbstractArray{T},f0::Real) where T <: AbstractFloat
    return exp.(T(-0.5) .* (A .- T(f0)) .^2) .* T(π ^ -0.25)
end


text:
macro proto( name )
    ex = quote
        struct $name{NT<:NamedTuple}
            properties::NT
        end # struct

        $name(;kwargs...) = $name(kwargs.data)

        function Base.getproperty( o::$name, s::Symbol )
            return getproperty( getfield(o, :properties), s )
        end # function

        function Base.propertynames( o::$name )
            return propertynames( getfield(o, :properties) )
        end # function
    end # quote
    esc(ex)
end # macro


text:
struct RasterFlags
    is_raster::Bool
    is_pairwise::Bool
    is_advanced::Bool
    is_onetoall::Bool
    is_alltoone::Bool
    grnd_file_is_res::Bool
    policy::Symbol
    four_neighbors::Bool
    avg_res::Bool
    solver::String
    outputflags::OutputFlags
end

function raster_pairwise(T, V, cfg)::Matrix{T}

    # Get input
    rasterdata = load_raster_data(T, V, cfg)

    # Get compute flags
    flags = get_raster_flags(cfg)

    pt_file_contains_polygons = length(rasterdata.points_rc[1]) !=
                                length(unique(rasterdata.points_rc[3]))

    if pt_file_contains_polygons
        _pt_file_polygons_path(rasterdata, flags, cfg)
    else
        _pt_file_no_polygons_path(rasterdata, flags, cfg)
    end
end

function get_raster_flags(cfg)

    # Computation flags
    is_raster = true
    is_pairwise = cfg["scenario"] in PAIRWISE
    is_advanced = cfg["scenario"] in ADVANCED
    is_onetoall = cfg["scenario"] in ONETOALL
    is_alltoone = cfg["scenario"] in ALLTOONE
    four_neighbors = cfg["connect_four_neighbors_only"] in TRUELIST
    avg_res = cfg["connect_using_avg_resistances"] in TRUELIST
    solver = cfg["solver"]
    ground_file_is_resistances =
        cfg["ground_file_is_resistances"] in TRUELIST
    policy = Symbol(cfg["remove_src_or_gnd"])

    # Output Flags
    o = get_output_flags(cfg)

    RasterFlags(is_raster, is_pairwise, is_advanced,
                is_onetoall, is_alltoone,
                ground_file_is_resistances, policy,
                four_neighbors, avg_res, solver, o)
end

function _pt_file_no_polygons_path(rasterdata::RasData{T,V},
                    flags, cfg)::Matrix{T} where {T,V}

    graphdata = compute_graph_data_no_polygons(rasterdata, flags)
    r = single_ground_all_pairs(graphdata, flags, cfg)

    write_cum_maps(graphdata.cum, rasterdata.cellmap, cfg, rasterdata.hbmeta,
                    flags.outputflags.write_max_cur_maps,
                    flags.outputflags.write_cum_cur_map_only)

    r
end

function _pt_file_polygons_path(rasterdata::RasData{T,V},
                        flags, cfg)::Matrix{T} where {T,V}

    # get unique list of points
    # for every point pair do
        # construct new polymap
        # construct new nodemap
        # construct new graph
        # solve for two points
    # end

    # Data
    gmap = rasterdata.cellmap
    polymap = rasterdata.polymap
    points_rc = rasterdata.points_rc
    avg_res = flags.avg_res
    four_neighbors = flags.four_neighbors

    # Cumulative maps
    cum = initialize_cum_maps(gmap, flags.outputflags.write_max_cur_maps)

    pts = unique(points_rc[3])
    resistances = -1 * ones(length(pts), length(pts))

    n = calc_num_pairs(pts)
    csinfo("Total number of pair solves = $n", cfg["suppress_messages"] in TRUELIST)

    k = 1
    for i = 1:size(pts, 1)
        pt1 = pts[i]
        for j = i+1:size(pts, 1)
            pt2 = pts[j]
            csinfo("Solving pair $k of $n", cfg["suppress_messages"] in TRUELIST)
            k += 1
            graphdata = compute_graph_data_polygons(rasterdata, flags, pt1, pt2, cum)
            pairwise_resistance = single_ground_all_pairs(graphdata, flags, cfg, false)
            resistances[i,j] = resistances[j,i] = pairwise_resistance[2,3]
        end
    end
    for i = 1:size(pts, 1)
        resistances[i,i] = 0
    end
    P = [0, pts...]
    r = hcat(P, vcat(pts', resistances))

    write_cum_maps(cum, gmap, cfg, rasterdata.hbmeta,
                   flags.outputflags.write_max_cur_maps,
                   flags.outputflags.write_cum_cur_map_only)

    # save resistances
    save_resistances(r, cfg)

    r
end

function calc_num_pairs(pts)
    n = 0
    for i = 1:size(pts, 1)
        for j = i+1:size(pts, 1)
            n += 1
        end
    end
    n
end

function compute_graph_data_polygons(rasterdata::RasData{T,V},
                            flags, pt1, pt2, cum)::GraphData{T,V} where {T,V}

    # Data
    gmap = rasterdata.cellmap
    polymap = rasterdata.polymap
    points_rc = rasterdata.points_rc
    hbmeta = rasterdata.hbmeta

    # Flags
    avg_res = flags.avg_res
    four_neighbors = flags.four_neighbors

    # Construct new polymap
    newpoly = create_new_polymap(gmap, polymap, points_rc, pt1, pt2)
    nodemap = construct_node_map(gmap, newpoly)

    # Construct graph
    a = construct_graph(gmap, nodemap, avg_res, four_neighbors)
    G = laplacian(a)

    # Find connected components
    cc = connected_components(SimpleWeightedGraph(a))

    # Construct points vector
    x,y = 0,0
    # x = find(x -> x == pt1, points_rc[3])[1]
    # y = find(x -> x == pt2, points_rc[3])[1]
    x = something(findfirst(isequal(pt1), points_rc[3]), 0)
    y = something(findfirst(isequal(pt2), points_rc[3]), 0)
    c1 = nodemap[points_rc[1][x], points_rc[2][x]]
    c2 = nodemap[points_rc[1][y], points_rc[2][y]]
    points = V[c1, c2]

    # Exclude pairs array
    exclude_pairs = Tuple{V,V}[]

    GraphData(G, cc, points, [pt1, pt2],
            exclude_pairs, nodemap, newpoly, hbmeta, gmap, cum)
end

#=function compute_graph_data(rasterdata, flags)

    points_rc = rasterdata.points_rc

    pt_file_contains_polygons =
            length(points_rc[1]) != length(unique(points_rc[3]))

    if !pt_file_contains_polygons
        graphdata = _pairwise_no_polygons(rasterdata, flags)
    end

    graphdata
end=#

function compute_graph_data_no_polygons(data::RasData{T,V},
                                        flags)::GraphData{T,V} where {T,V}

    # Data
    cellmap = data.cellmap
    polymap = data.polymap
    points_rc = data.points_rc
    included_pairs = data.included_pairs
    hbmeta = data.hbmeta

    # Flags
    avg_res = flags.avg_res
    four_neighbors = flags.four_neighbors
    write_max_cur_maps = flags.outputflags.write_max_cur_maps

    # Nodemap and graph construction
    nodemap = construct_node_map(cellmap, polymap)
    G = construct_graph(cellmap, nodemap, avg_res, four_neighbors)
    G = laplacian(G)

    # Connected Components
    cc = connected_components(SimpleWeightedGraph(G))

    # Generate exclude pairs array
    if !isempty(included_pairs)
        exclude_pairs = generate_exclude_pairs(points_rc, included_pairs)
    else
        exclude_pairs = Tuple{V,V}[]
    end

    points = zeros(V, length(points_rc[3]))
    for (i,v) in enumerate(zip(points_rc[1], points_rc[2]))
        points[i] = nodemap[v...]
    end

    # Cumulative current maps
    cum = initialize_cum_maps(cellmap, write_max_cur_maps)

    GraphData(G, cc, points, points_rc[3],
                exclude_pairs, nodemap, polymap,
                hbmeta, cellmap, cum)
end
Base.isempty(t::IncludeExcludePairs) = t.mode == :undef

function generate_exclude_pairs(points_rc, included_pairs::IncludeExcludePairs{V}) where V

    exclude_pairs_array = Tuple{V,V}[]
    mat = included_pairs.include_pairs
    mode = included_pairs.mode == :include ? 0 : 1

    prune_points!(points_rc, included_pairs.point_ids)
    for j = 1:size(mat, 2)
        for i = 1:size(mat, 1)
            if mat[i,j] == mode && mat[j,i] == mode
                push!(exclude_pairs_array, (i,j))
            end
        end
    end

    exclude_pairs_array
end

function construct_node_map(gmap, polymap::Matrix{V}) where V

    nodemap = zeros(V, size(gmap))
    ind = gmap .> 0
    nodemap[ind] = 1:sum(ind)

    if isempty(polymap)
        return nodemap
    end

    #=d = Dict{Int, Vector{Int}}()
    #=for i in unique(polymap)
        d[i] = find(x -> x == i, polymap)
    end=#
    m, n = size(polymap)
    I, J, V = findnz(polymap)
    for (i,v) in enumerate(V)
        if v != -9999
            if haskey(d, v)
                push!(d[v], sub2ind((m, n), I[i], J[i]))
            else
                d[v] = [sub2ind((m, n), I[i], J[i])]
            end
        end
    end
    d[0] = find(x -> x == 0, polymap)
    k = 1
    for i in find(gmap)
        if i in d[0]
            nodemap[i] = k
            k += 1
        else
            for key in keys(d)
                if i in d[key]
                    if i == first(d[key])
                        nodemap[i] = k
                        k += 1
                    else
                        f = first(d[key])
                        if polymap[f] != 0 && nodemap[f] == 0
                            nodemap[i] = k
                            nodemap[f] = k
                            k += 1
                        else
                            nodemap[i] = nodemap[f]
                        end
                    end
                end
            end
        end
    end
    for i in eachindex(polymap)
        if polymap[i] != 0 && nodemap[i] == 0
            val::Float64 = polymap[i]
            index::INT = findfirst(x -> x == val, polymap)
            nodemap[i] = nodemap[index]
        end
    end=#

    idx = gmap .> 0
    polymap_pruned = zeros(V, size(gmap))
    polymap_pruned[idx] = polymap[idx]

    polynums = unique(polymap)
    for i = 1:size(polynums, 1)
        polynum = polynums[i]
        if polynums[i] != 0
            idx1 = findall(x -> x == polynum, polymap_pruned)
            idx2 = findall(x -> x == polynum, polymap)
            if length(idx1) > 0
                nodemap[idx2] .= nodemap[idx1[1]]
            end
        end
    end
    relabel!(nodemap, V(1))

    nodemap
end

function relabel!(nodemap::Matrix{V}, offset = V(0)) where V
    oldlabels = nodemap[findall(x->x!=0,nodemap)]
    newlabels = zeros(V, size(oldlabels))
    s = sort(oldlabels)
    perm = sortperm(oldlabels)
    prepend!(s, s[1] - 1)
    f = findall(x->x!=0,diff(s))
    newlabels[f] .= 1
    newlabels = cumsum(newlabels)
    newlabels[perm] = copy(newlabels)
    nodemap[findall(x->x!=0,nodemap)] = newlabels .- V(1) .+ offset
end

function construct_graph(gmap, nodemap::Matrix{S}, avg_res, four_neighbors) where S
    f1 = avg_res ? res_avg : cond_avg
    f2 = avg_res ? weirder_avg : weird_avg
    I = Vector{S}()
    J = Vector{S}()
    V = Vector{eltype(gmap)}()
    for j = 1:size(gmap, 2)
        for i = 1:size(gmap, 1)
            if nodemap[i,j] == 0
                continue
            else
                # Horizontal neighbour
                if j != size(gmap, 2) && nodemap[i,j+1] != 0
                    push!(I, nodemap[i,j])
                    push!(J, nodemap[i,j+1])
                    push!(V, f1(gmap[i,j], gmap[i,j+1]))
                end

                # Vertical neighbour
                if i != size(gmap, 1) && nodemap[i+1, j] != 0
                    push!(I, nodemap[i,j])
                    push!(J, nodemap[i+1,j])
                    push!(V, f1(gmap[i,j], gmap[i+1,j]))
                end

                if !four_neighbors
                    # Diagonal neighbour
                    if i != size(gmap, 1) && j != size(gmap, 2) && nodemap[i+1, j+1] != 0
                        push!(I, nodemap[i,j])
                        push!(J, nodemap[i+1,j+1])
                        push!(V, f2(gmap[i,j], gmap[i+1,j+1]))
                    end

                    if i != 1 && j != size(gmap, 2) && nodemap[i-1, j+1] != 0
                        push!(I, nodemap[i,j])
                        push!(J, nodemap[i-1,j+1])
                        push!(V, f2(gmap[i,j], gmap[i-1,j+1]))
                    end
                end
            end
        end
    end
    m = maximum(nodemap)
    a = sparse(I,J,V, m, m)
    a = a + a'
    a
end

res_avg(x, y) = 1 / ((1/x + 1/y) / 2)
cond_avg(x, y) = (x + y) / 2
weird_avg(x,y) = (x + y) / (2*√2)
weirder_avg(x, y) = 1 / (√2 * (1/x + 1/y) / 2)

function create_new_polymap(gmap, polymap::Matrix{V}, points_rc,
                pt1 = 0, pt2 = 0, point_map = Matrix{V}(undef,0,0)) where V

    f(x) = (points_rc[1][x], points_rc[2][x])

    if !isempty(point_map)
        # Combine polymap and pointmap
        newpoly = deepcopy(polymap)
        point_file_no_polygons = length(points_rc[3]) ==
                        length(unique(points_rc[3]))
        if isempty(polymap)
            newpoly = point_map
        elseif point_file_no_polygons
            k = maximum(polymap)
            for i in findall(x->x!=0,point_map)
                if polymap[i] == 0
                    newpoly[i] = point_map[i] + k
                end
            end
        else
            k = max(maximum(polymap), maximum(point_map))
            for i in findall(x->x!=0,point_map)
                v1 = point_map[i]
                v2 = newpoly[i]
                if v2 == 0
                    newpoly[i] = k + v1
                    continue
                end
                if v1 != v2
                    ind = findall(x -> x == v2, newpoly)
                    newpoly[ind] .= v1
                end
            end
        end
        return newpoly
    end

    if isempty(polymap)
        newpoly = zeros(V, size(gmap)...)
        id1 = findall(x -> x == pt1, points_rc[3])
        id2 = findall(x -> x == pt2, points_rc[3])
        map(x -> newpoly[f(x)...] = pt1, id1)
        map(x -> newpoly[f(x)...] = pt2, id2)
        return newpoly
    else
        newpoly = deepcopy(polymap)
        k = maximum(polymap)
        for p in (pt1, pt2)
            # find the locations of the point
            idx = findall(x -> x == p, points_rc[3])

            if length(idx) == 1
                continue
            end
            allzero = mapreduce(x -> polymap[f(x)...] == 0, &, idx)
            if allzero
                map(x -> newpoly[f(x)...] = k + 1, idx)
                k += 1
            else
                nz = filter(x -> polymap[f(x)...]!= 0, idx)
                if length(nz) == 1
                    map(x -> newpoly[f(x)...] = polymap[overlap[1]], idx)
                else
                    coords = map(x -> f(x), nz)
                    vals = map(x -> polymap[x...], coords)
                    overlap = findall(in(vals), polymap)
                    newpoly[overlap] .= k + 1
                    k += 1
                end
            end
        end
        return newpoly
    end
end