text:
function ugly_print(i::UInt64)
    bit_array = cvt_to_bitarray(i)
    uglyPrint(bit_array)
end

function ugly_print(pieces_array::Array{UInt64,1})
    pieces_uint = cvt_to_uint(pieces_array)
    ugly_print(pieces_uint)
end

function ugly_print(valid_moves::Array{Any,1})
    pieces = EMPTY
    for move in valid_moves
        pieces |= move[2]
    end
    ugly_print(pieces)
end

function uglyPrint(b::BitArray)

    r_b = Int.(transpose(reshape(b, 8, :)))
    ranks = ["8", "7", "6", "5", "4", "3", "2", "1"]

    @printf("\n  o-----------------o\n")
    for i = 1:8
        @printf("%s | ", ranks[i])
        for j = 1:8
            @printf("%d ", r_b[i,j])
        end
        @printf("|\n")
    end
    @printf("  o-----------------o\n")
    @printf("    a b c d e f g h\n")
end

function pretty_print(board::Bitboard, player_color::String="white")
    ranks = ["8", "7", "6", "5", "4", "3", "2", "1"]

    free = transpose(reshape(cvt_to_bitarray(board.free), 8, :))
    taken = transpose(reshape(cvt_to_bitarray(board.taken), 8, :))

    p = transpose(reshape(cvt_to_bitarray(board.p), 8, :))
    r = transpose(reshape(cvt_to_bitarray(board.r), 8, :))
    n = transpose(reshape(cvt_to_bitarray(board.n), 8, :))
    b = transpose(reshape(cvt_to_bitarray(board.b), 8, :))
    q = transpose(reshape(cvt_to_bitarray(board.q), 8, :))
    k = transpose(reshape(cvt_to_bitarray(board.k), 8, :))
    P = transpose(reshape(cvt_to_bitarray(board.P), 8, :))
    R = transpose(reshape(cvt_to_bitarray(board.R), 8, :))
    N = transpose(reshape(cvt_to_bitarray(board.N), 8, :))
    B = transpose(reshape(cvt_to_bitarray(board.B), 8, :))
    Q = transpose(reshape(cvt_to_bitarray(board.Q), 8, :))
    K = transpose(reshape(cvt_to_bitarray(board.K), 8, :))
    white = transpose(reshape(cvt_to_bitarray(board.white), 8, :))
    black = transpose(reshape(cvt_to_bitarray(board.black), 8, :))

    pieces = Dict("pawn"=>" o",
        "rook"=>" Π",
        "knight"=>" ζ",
        "bishop"=>" Δ",
        "queen"=>" Ψ",
        "king"=>" +")
    # pieces = Dict("pawn"=>" o",
    #   "rook"=>" R",
    #   "knight"=>" N",
    #   "bishop"=>" B",
    #   "queen"=>" Q",
    #   "king"=>" K")
    labels = ["pawn", "knight", "bishop", "rook", "queen", "king"]

    @printf("\n  o-------------------------o\n")

    if player_color == "white"
        idxs = 1:8
    else
        idxs = 8:-1:1
    end

    bgc = "w"
    for i in idxs
        @printf(Crayon(reset=true), "%s | ", ranks[i])
        for j in idxs
            if free[i,j]
                if bgc == "w"
                    @printf(Crayon(reset=true, background=:dark_gray), "   ")
                    bgc = "b"
                else
                    @printf(Crayon(reset=true, background=:default), "   ")
                    bgc = "w"
                end
            else
                if black[i,j]
                    if p[i,j]
                        c = pieces["pawn"]
                    elseif r[i,j]
                        c = pieces["rook"]
                    elseif n[i,j]
                        c = pieces["knight"]
                    elseif b[i,j]
                        c = pieces["bishop"]
                    elseif q[i,j]
                        c = pieces["queen"]
                    elseif k[i,j]
                        c = pieces["king"]
                    else
                        error("Black piece not found")
                    end
                    color = :light_magenta
                else
                    if P[i,j]
                        c = pieces["pawn"]
                    elseif R[i,j]
                        c = pieces["rook"]
                    elseif N[i,j]
                        c = pieces["knight"]
                    elseif B[i,j]
                        c = pieces["bishop"]
                    elseif Q[i,j]
                        c = pieces["queen"]
                    elseif K[i,j]
                        c = pieces["king"]
                    else
                        error("White piece not found")
                    end
                    color = :light_cyan
                end
                if bgc == "w"
                    @printf(Crayon(bold=true, foreground=color,
                        background=:dark_gray), "%s ", c)
                    bgc = "b"
                else
                    @printf(Crayon(bold=true, foreground=color,
                        background=:default), "%s ", c)
                    bgc = "w"
                end

            end
        end

        if i > 1 && i < 8
            label = labels[i-1]
            piece = pieces[label]
            @printf(Crayon(reset=true), "| %s %s \n", piece, label)
        else
            @printf(Crayon(reset=true), "|\n")
        end
        if bgc == "w"
            bgc = "b"
        else
            bgc = "w"
        end
    end
    @printf(Crayon(reset=true), "  o-------------------------o\n")

    if player_color == "white"
        @printf("     a  b  c  d  e  f  g  h\n")
    else
        @printf("     h  g  f  e  d  c  b  a\n")
    end
end


text:
using Dates
using DataFrames
using Statistics
using Test
using TSx

# constants
DATA_SIZE = 400

# global variables
data_vector = randn(DATA_SIZE)
data_array = Array([data_vector data_vector])

index_range = 1:DATA_SIZE
index_integer = collect(index_range)
index_timetype = Date(2007, 1,1) + Day.(0:(DATA_SIZE - 1))

df_vector = DataFrame([data_vector], ["data"])
df_integer_index = DataFrame(Index = index_integer, data = data_vector)
df_timetype_index = DataFrame(Index = index_timetype, data = data_vector)

@testset "TS()" begin
    include("TS.jl")
end

@testset "getindex()" begin
    include("getindex.jl")
end

@testset "index()" begin
    include("index.jl")
end

ext: jl


text:
# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function simple_polish_fastq(simple_kmer_graph, fastq_file; min_depth=3)
#     solid_vertices = filter(v -> simple_kmer_graph.vprops[v][:weight] >= min_depth, Graphs.vertices(simple_kmer_graph))
#     filtered_simple_kmer_graph, vertex_map = Graphs.induced_subgraph(simple_kmer_graph, solid_vertices)
# #     display(simple_kmer_graph)
# #     display(filtered_simple_kmer_graph)
#     kmers = sort!(graph_to_kmers(filtered_simple_kmer_graph))

#     old_kmers = graph_to_kmers(simple_kmer_graph)

#     k = filtered_simple_kmer_graph.gprops[:k]

#     polished_fastq_file = replace(fastq_file, ".fastq" => ".k$k.d$(min_depth).fastq")

#     transition_probabilities = Mycelia.initialize_transition_probabilities(filtered_simple_kmer_graph)
#     state_likelihoods = [Float64(filtered_simple_kmer_graph.vprops[v][:weight]) for v in Graphs.vertices(filtered_simple_kmer_graph)]
#     state_likelihoods ./= sum(state_likelihoods)

# #     @info "counting the number of records to establish runtime estimate"
#     number_of_records = 0
#     for fastq_record in FASTX.FASTQ.Reader(open(fastq_file))
#         number_of_records += 1
#     end
#     progress_bar = ProgressMeter.Progress(number_of_records, 1)

#     fastq_reader = FASTX.FASTQ.Reader(open(fastq_file))
#     fastq_writer = FASTX.FASTQ.Writer(open(polished_fastq_file, "w"))

#     for fastx_record in fastq_reader
#         ProgressMeter.next!(progress_bar)
#         bubble_start = 0
#         updated_path = Vector{Pair{Int, Bool}}()
#     #     @show FASTX.sequence(fastx_record)
#         for (i, kmer) in enumerate(BioSequences.each(BioSequences.BigDNAMer{k}, FASTX.sequence(fastx_record)))
#             canonical_kmer = min(kmer.fw, kmer.bw)
#             orientation = canonical_kmer == kmer.fw
#             kmer_index_range = searchsorted(kmers, canonical_kmer)

#             old_kmer_index = searchsortedfirst(old_kmers, canonical_kmer)

# #             FASTX.identifier(fastx_record) == "4" && @show kmer_index_range
# #             FASTX.identifier(fastx_record) == "4" && @show kmers[kmer_index_range]
# #             FASTX.identifier(fastx_record) == "4" && @show old_kmers[old_kmer_index]
#             kmer_is_solid = !isempty(kmer_index_range)
# #             FASTX.identifier(fastx_record) == "4" && @show canonical_kmer
# #             FASTX.identifier(fastx_record) == "4" &&  @show kmer_is_solid

#             if kmer_is_solid
#                 kmer_index = first(kmer_index_range)
# #                 FASTX.identifier(fastx_record) == "4" && @show filtered_simple_kmer_graph.vprops[kmer_index], simple_kmer_graph.vprops[old_kmer_index]
#             else
#                 kmer_index = 0
#             end

#             in_bubble = bubble_start > 0

#             if !kmer_is_solid
#                 if !in_bubble
# #                     FASTX.identifier(fastx_record) == "4" && @show "starting a bubble"
#                     bubble_start = i
#                 else
# #                     FASTX.identifier(fastx_record) == "4" && @show "continuing in a bubble"
#                 end
#             else
#                 if !in_bubble
# #                     FASTX.identifier(fastx_record) == "4" && @show "pushing solid kmer to updated path"
#                     push!(updated_path, kmer_index => orientation)
#                 else
#                     if bubble_start == 1
# #                         FASTX.identifier(fastx_record) == "4" && @show "ending an opening bubble"
#                         # we're in a bubble that began at the beginning of the read
#                         # we'll do nothing and just remove this
#                         # equivalent to tip clipping
# #                         FASTX.identifier(fastx_record) == "4" && @show "pushing solid kmer to updated path"
#                         push!(updated_path, kmer_index => orientation)
#                         bubble_start = 0
#                     else
# #                         FASTX.identifier(fastx_record) == "4" && @show "found end of an anchored bubble -- correcting"
#                         source_vertex, source_orientation = last(updated_path)
#                         destination_vertex, destination_orientation = kmer_index, orientation

#                         shortest_paths = Graphs.yen_k_shortest_paths(
#                             filtered_simple_kmer_graph,
#                             source_vertex,
#                             destination_vertex,
#                             Graphs.weights(filtered_simple_kmer_graph),
#                             10).paths

#                         if isempty(shortest_paths)
#                             @show source_vertex, destination_vertex
#                             @warn "no valid alternate paths found for $(FASTX.identifier(fastx_record)), continuing"
#                             break
# #                             @show fastx_record
# #                             error("try increasing min_depth")
#                         end
#                         candidate_path_probabilities = ones(length(shortest_paths))
#                         oriented_candidate_paths = [
#                             [last(updated_path)] for i in 1:length(shortest_paths)
#                         ]

#                         for (i, candidate_path) in enumerate(shortest_paths)
#                             for dest_vertex in candidate_path[2:end]
#                                 source_vertex, source_orientation = last(oriented_candidate_paths[i])
#                                 candidate_path_probabilities[i] *= transition_probabilities[source_orientation][source_vertex, dest_vertex]
#                                 candidate_path_probabilities[i] *= state_likelihoods[dest_vertex]
#                                 if candidate_path_probabilities[i] > 0
#                                     edge = Graphs.Edge(source_vertex, dest_vertex)
#                                     destination_orientation =
#                                     first(
#                                         filter(o -> o.source_orientation == source_orientation,
#                                             filtered_simple_kmer_graph.eprops[edge][:orientations])).destination_orientation
#                                     push!(oriented_candidate_paths[i], (dest_vertex => destination_orientation))
#                                 else
#                                     break # this path is no good, evaluate the next
#                                 end
#                             end
#                         end
#                         non_zero_indices = findall(p -> p > 0, candidate_path_probabilities)
#                         if isempty(non_zero_indices)
# #                             @show candidate_path_probabilities
#                             @warn "no resolution found for read $(FASTX.identifier(fastx_record))"
#                             break
# #                             error("no valid alternate path probabilities")
# #                             error("try increasing min_depth?")
#                         end

#                         candidate_path_probabilities = candidate_path_probabilities[non_zero_indices]
#                         oriented_candidate_paths = oriented_candidate_paths[non_zero_indices]

#                         # offset is for debugging
#                         # make sure that anchors on both sides are the same
#                         offset = 0
#                         observed_sequence = FASTX.sequence(fastx_record)[bubble_start+k-1-offset:i-1+offset]
#                         for (i, oriented_candidate_path) in enumerate(oriented_candidate_paths)
#                             candidate_sequence = Mycelia.oriented_path_to_sequence(
#                                 filtered_simple_kmer_graph,
#                                 oriented_candidate_path)
#                             candidate_sequence = candidate_sequence[k+1-offset:end-k+offset]
#                             alignment_result = BioAlignments.pairalign(
#                                 BioAlignments.LevenshteinDistance(),
#                                 candidate_sequence,
#                                 observed_sequence)
# #                             @show alignment_result
#     #                         @show alignment_result
#                             average_error_rate = Statistics.mean(Mycelia.q_value_to_error_rate.(FASTX.quality(fastx_record)))
#                             for error in 1:alignment_result.value
#                                 candidate_path_probabilities[i] *= average_error_rate
#                             end
#                             for match in 1:BioAlignments.count_matches(alignment_result.aln)
#                                 candidate_path_probabilities[i] *= (1 - average_error_rate)
#                             end
#                         end

#                         chosen_replacement = StatsBase.sample(oriented_candidate_paths, StatsBase.weights(candidate_path_probabilities))

#                         for i in 2:length(chosen_replacement)
#                             oriented_state = chosen_replacement[i]
#                             push!(updated_path, oriented_state)
#                         end
#                         bubble_start = 0
#                     end
#                 end
#             end
#         end
#     #     @show updated_path
#         sequence = Mycelia.oriented_path_to_sequence(filtered_simple_kmer_graph, updated_path)
#         alignment_result = BioAlignments.pairalign(
#             BioAlignments.LevenshteinDistance(),
#             sequence,
#             FASTX.sequence(fastx_record))
# #         if alignment_result.value > 0
# # #             @show alignment_result
# #         end
#         quality = StatsBase.sample(FASTX.quality(fastx_record), length(sequence), replace=true, ordered=true)
#         description =  join(filter(!isempty, (FASTX.description(fastx_record), "k$k.d$(min_depth)")), '.')
#         identifier = FASTX.identifier(fastx_record)
#         new_record = FASTX.FASTQ.Record(identifier, description, sequence, quality)
#         write(fastq_writer, new_record)
#     end
#     close(fastq_reader)
#     close(fastq_writer)
#     return polished_fastq_file
# end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function assess_observations(graph::Mycelia.KmerGraph{KMER_TYPE}, observations, error_rate; verbose = isinteractive()) where {KMER_TYPE}
#     k = last(KMER_TYPE.parameters)
#     total_edits_accepted = 0
#     total_bases_evaluated = 0
#     reads_processed = 0
#     maximum_likelihood_observations = Vector{BioSequences.LongDNASeq}(undef, length(observations))
#     for (observation_index, observation) in enumerate(observations)
#         if length(observation) >= k
#             optimal_path, edit_distance, relative_likelihood = Mycelia.viterbi_maximum_likelihood_path(graph, observation, error_rate)
#             maximum_likelihood_observation = Mycelia.oriented_path_to_sequence(optimal_path, graph.kmers)
#             maximum_likelihood_observations[observation_index] = maximum_likelihood_observation
#             reads_processed += 1
#             total_bases_evaluated += length(observation)
#             total_edits_accepted += edit_distance
#         else
#             maximum_likelihood_observations[observation_index] = observation
#         end
#     end
#     inferred_error_rate = round(total_edits_accepted / total_bases_evaluated, digits = 3)
#     if verbose
#         display("reads_processed = $(reads_processed)")
#         display("total_edits_accepted = $(total_edits_accepted)")
#         display("inferred_error_rate = $(inferred_error_rate)")
#     end
#     if total_edits_accepted == 0
#         has_converged = true
#     else
#         has_converged = false
#     end
#     return maximum_likelihood_observations, has_converged
# end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function iterate_until_convergence(ks, observations, error_rate)
#     graph = nothing
#     for k in ks
#         graph = Mycelia.KmerGraph(BioSequences.DNAMer{k}, observations)
#         observations, has_converged = assess_observations(graph, observations, error_rate; verbose = verbose)
#     end
#     return graph, observations
# end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function clip_low_coverage_tips(graph, observations)
#     connected_components = Graphs.connected_components(graph.graph)
#     vertices_to_keep = Int[]
#     for connected_component in connected_components

#         component_coverage = graph.counts[connected_component]
#         median = Statistics.median(component_coverage)
#         standard_deviation = Statistics.std(component_coverage)

#         for vertex in connected_component
#             keep = true
#             if Graphs.degree(graph.graph, vertex) == 1
#                 this_coverage = graph.counts[vertex]
#                 is_low_coverage = (graph.counts[vertex] == 1) ||
#                                     (median-this_coverage) > (3*standard_deviation)
#                 if is_low_coverage
#                     keep = false
#                 end
#             end
#             if keep
#                 push!(vertices_to_keep, vertex)
#             end
#         end
#     end

#     KmerType = first(typeof(graph).parameters)
#     pruned_graph = Mycelia.KmerGraph(KmerType, observations, graph.kmers[vertices_to_keep], graph.counts[vertices_to_keep])

#     return pruned_graph
# end

max_stars_repo_name: snystrom/Mycelia

max_issues_repo_name: snystrom/Mycelia

max_forks_repo_name: cjprybol/Eisenia


text:
#----------------- Get the A matrix
function getDivGrad(n1,n2,n3)
    # the Divergence
    D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
    D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
    D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))
    # DIV from faces to cell-centers

    Div = [D1 D2 D3]

    return Div*Div'
end
#----------------- 1D finite difference on staggered grid

# generate 1D derivatives
ddx(n)=spdiags(ones(n)*[-1 1],[0,1],n,n+1)

#------------- Build a diagonal matrix
function spdiags(B,d,m,n)
    # spdiags(B,d,m,n)
    # creates a sparse matrix from its diagonals
    d = d[:]
    p = length(d)

    len = zeros(p+1,1)
    for k = 1:p
        len[k+1] = Int(len[k]+length(max(1,1-d[k]):min(m,n-d[k])))
    end
    a = zeros(round(Int, len[p+1]), 3)
    for k = 1:p
        # Append new d[k]-th diagonal to compact form
        i = max(1,1-d[k]):min(m,n-d[k])
        a[(round(Int, len[k])+1):round(Int, len[k+1]),:] = [i i+d[k] B[i+(m >= n)*d[k], k]]
    end

    sparse(round(Int, a[:,1]), round(Int, a[:,2]), a[:,3], m, n)
end


text:
using PyCall

@pyimport gym

function play_env(n::Network, env, pcfg::Dict, seed::Int64, trial::Int64,
                  train::Bool=true)
    env[:seed](seed)
    ob = env[:reset]()
    ehigh = min.(env[:observation_space][:high], 20)
    elow = max.(env[:observation_space][:low], -20)
    ma_reward = -Inf
    max_reward = 0.0
    reward = -Inf
    done = false
    step = 0
    bad_step = 0

    while ~done
        nob = min.(max.(((ob - elow) ./ (ehigh - elow)), 0.0), 1.0) .* pcfg["fr"]
        ins = rand(n.nin, pcfg["tinput"])
        inputs = falses(n.nin, pcfg["tinput"])
        for i in 1:n.nin
            inputs[i, :] = ins[i, :] .< nob[mod(i, length(nob))+1]
        end
        for t in 1:pcfg["tinput"]
            outs = step!(n, inputs[:, t], train)
        end
        ocount = zeros(pcfg["n_actions"])
        for t in 1:pcfg["toutput"]
            outs = step!(n, train)
            for o in 1:length(outs)
                if outs[o]
                    ocount[mod(o, length(ocount))+1] += 1.0
                end
            end
        end
        for t in 1:pcfg["trest"]
            step!(n, train)
        end
        maxcount = maximum(ocount)
        mincount = minimum(ocount)
        if maxcount == mincount
            bad_step += 1
        else
            bad_step = 0
        end
        ob, r, done, _ = env[:step](indmax(ocount)-1)
        reward = r
        if pcfg["id"] == "Acrobot-v1"
            reward = acos(ob[1]) # max angle
            max_reward = max(max_reward, reward)
        elseif pcfg["id"] == "CartPole-v0"
            reward = 1.0 - acos(cos(ob[3])) # 1 - angle
            max_reward += 1.0
        elseif pcfg["id"] == "MountainCar-v0"
            reward = (ob[1] - elow[1]) / (ehigh[1] - elow[1])  # mountain car
            max_reward = max(max_reward, reward)
        else
            max_reward += reward
        end
        if ma_reward == -Inf
            ma_reward = reward
        else
            ma_reward = (1.0 - pcfg["ma_rate"]) * ma_reward + pcfg["ma_rate"] * reward
        end
        dop = 0.0
        if ma_reward != 0.0
            dop = reward / ma_reward
        end
        if pcfg["onpolicy"]
            if pcfg["dmethod"] == 1
                if dop > 1.0
                    n.da[1] += (dop - 1.0)
                end
            elseif pcfg["dmethod"] == 2
                n.da[1] *= dop
            elseif pcfg["dmethod"] == 3
                n.da[1] = dop
            end
        end
        step += 1
        Logging.info(@sprintf("S: %s %d %d %d %d %0.6f %0.6f %0.6f %0.6f",
                              pcfg["id"], trial, step, maxcount, mincount,
                              max_reward, reward, ma_reward, dop))
        if bad_step > 20
            return -1e4
        end
    end

    if (pcfg["id"] == "Acrobot-v1") || (pcfg["id"] == "MountainCar-v0")
        max_reward = max_reward - step
    end

    max_reward
end

function repeat_trials(n::Network, env, pcfg::Dict)
    init_weights = copy(n.weights)
    bad_trial = 0
    ma_reward = play_env(n, env, pcfg, 0, 0, false)
    if ma_reward == -1e4
        return ma_reward
    end

    for i in 1:pcfg["n_trials"]
        weights = copy(n.weights)
        reward = play_env(n, env, pcfg, i, i, true)
        if reward == -1e4
            return reward
        end
        dop = 0.0
        if ma_reward != 0.0
            dop = reward / ma_reward
        end
        if dop > 1.0 && ~pcfg["onpolicy"]
            n.da[1] += (dop - 1.0)
        end
        ma_reward = (1.0 - pcfg["ma_rate"]) * ma_reward + pcfg["ma_rate"] * reward
        dweight = sum(abs.(n.weights - weights)) / length(n.weights)
        Logging.info(@sprintf("T: %s %d %0.6f %0.6f %0.6f %e %e",
                              pcfg["id"], i, reward, ma_reward, dop,
                              sum(abs.(n.weights - init_weights)) / length(n.weights),
                              dweight))
        if dweight < 1e-4
            bad_trial += 1
        else
            bad_trial = 0
        end
        if bad_trial > 10
            return -1e4
        end
    end

    mean(x->play_env(n, env, pcfg, pcfg["seed"]+x, pcfg["n_trials"]+x, false), 1:3)
end


text:
"""
    cudnnSetTensor!(x, s)

Set all elements of tensor `x` to scalar `s` and return `x`.
"""
function cudnnSetTensor!(
    x, s::Real;
    format::cudnnTensorFormat_t = CUDNN_TENSOR_NCHW,
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x; format)
)
    cudnnSetTensor(handle(), xDesc, x, Ref(eltype(x)(s)))
    return x
end

"""
    cudnnScaleTensor(x, s)
    cudnnScaleTensor!(y, x, s)

Scale all elements of tensor `x` with scale `s` and return the result. `cudnnScaleTensor`
allocates a new array for the answer, `cudnnScaleTensor!` overwrites `y`.
"""
cudnnScaleTensor, cudnnScaleTensor!

function cudnnScaleTensor!(
    y, x, s::Real;
    format::cudnnTensorFormat_t = CUDNN_TENSOR_NCHW,
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x; format)
)
    y === x || copyto!(y, x)
    cudnnScaleTensor(handle(), xDesc, y, Ref(eltype(y)(s)))
    return y
end

cudnnScaleTensor(x, s::Real; o...) = cudnnScaleTensor!(similar(x), x, s; o...)

# cudnnAddTensor does not support all broadcasting dimensions, use cudnnOpTensor instead.
# Compared to libknet8 x .+ b it is ~2x slower for (1,1,100,100), ~30% faster for (14,14,256,32)
# CUDA.jl x .+ b is 2x slower than both

"""
    cudnnAddTensor(x, b; alpha)
    cudnnAddTensor!(y, x, b; alpha, beta)

Broadcast-add tensor `b` to tensor `x`. `alpha=1, beta=1` are used for scaling, i.e. `y .=
alpha * b .+ beta * x`.  `cudnnAddTensor` allocates a new array for the answer,
`cudnnAddTensor!` overwrites `y`. Does not support all valid broadcasting dimensions.  For
more flexible broadcast operations see `cudnnOpTensor`.
"""
cudnnAddTensor, cudnnAddTensor!

function cudnnAddTensor!(
    y, x, b;
    alpha::Real=1,
    beta::Real=1,
    format::cudnnTensorFormat_t = CUDNN_TENSOR_NCHW,
    bDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(b; format),
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x; format),
)
    T = eltype(x)
    alpha, beta = scalingParameter(T, alpha), scalingParameter(T, beta)
    y === x || copyto!(y, x)
    cudnnAddTensor(handle(), alpha, bDesc, b, beta, xDesc, y)
    return y
end

cudnnAddTensor(x, b; o...) = cudnnAddTensor!(similar(x), x, b; o...)


text:
using CMAEvolutionStrategy

if isinstalled("scipy.optimize")
    include("scipy_optimize.jl")
    using .SciPyOptimize
end

function local_search!(
        objective::Function,
        ip::Vector{Int64},
        population::Matrix{Float64},
        n_population::Int64,
        n_gene::Int64;
        method::String,
        n_children::Int64,
        maxiter::Int64)::Matrix{Float64}
    if method == "mutation"
        idx::BitArray{1} = trues(n_population)
        idx[ip[1]] = false

        children::Matrix{Float64} = zeros(n_children, n_gene + 1)

        for i in 1:n_children
            ip[2:end] = sample(
                collect(1:n_population)[idx], n_gene + 1, replace=false
            )
            children[i, :] = mutation(objective, population[ip, :], n_gene)
        end

        family::Matrix{Float64} = zeros(n_children + 1, n_gene + 1)
        @inbounds for i in 1:n_gene + 1
            @simd for j in 1:n_children
                family[j, i] = children[j, i]
            end
            family[n_children + 1, i] = population[ip[1], i]
        end

        family = sortslices(family, dims=1, by=x -> x[end])

        for i in 1:n_gene + 1
            @inbounds population[ip[1], i] = family[1, i]  # Best
        end
    elseif method == "powell"
        population = fmin_powell(objective, n_gene, population, ip, maxiter)
    elseif method == "de"
        population = fmin_de(objective, n_gene, population, ip, maxiter)
    elseif method == "cmaes"
        x0 = [population[ip[1], i] for i in 1:n_gene]
        lb = vec(minimum(population[:, 1:n_gene], dims=1))
        ub = vec(maximum(population[:, 1:n_gene], dims=1))
        s0 = 0.25 * median(ub - lb)
        result = CMAEvolutionStrategy.minimize(
            objective, x0, s0;
            lower=lb,
            upper=ub,
            verbosity=0,
            multi_threading=true,
            maxiter=maxiter,
            maxfevals=1000 * n_gene,
        )
        x_best_fit = xbest(result)
        obj_val = objective(x_best_fit)
        if obj_val < objective(x0)
            for i in 1:n_gene
                @inbounds population[ip[1], i] = x_best_fit[i]
            end
            population[ip[1], end] = obj_val
        end
    end
    population = sortslices(population, dims=1, by=x -> x[end])

    return population
end

function mutation(
        objective::Function,
        parents::Matrix{Float64},
        n_gene::Int64)::Vector{Float64}
    # Normal Distribution Mutation
    child::Vector{Float64} = zeros(n_gene + 1)

    GAMMA::Float64 = 0.35 / sqrt(n_gene)

    t2::Vector{Float64} = zeros(n_gene)
    centroid::Vector{Float64} = reshape(
        mean(parents[2:end, 1:n_gene], dims=1), n_gene
    )

    for i in 1:n_gene + 1
        t2 += randn() * GAMMA * (parents[i + 1, 1:n_gene] - centroid)
    end

    @inbounds for i in 1:n_gene
        child[i] = clamp(parents[1, i] + t2[i], 0.0, 1.0)
    end

    child[end] = objective(child[1:n_gene])

    return child
end


text:
"""
	SeisRemove(filename)

Delete a seis file (deletes the text file, binary data file, and binary header file if there is one)

*Credits: AS, 2015*

"""
function SeisRemove(filename::AbstractString)

	filename_data = ParseDataName(filename)
	filename_headers = ParseHeaderName(filename)
	rm(filename);
	rm(filename_data);
	if filename_headers != "NULL"
		rm(filename_headers);
	end

end

function SeisRemove(filename::Array{AbstractString,1})

	for j = 1 : length(filename)
		SeisRemove(filename[j])
	end

end


text:
using LinearAlgebra
using LinearModel
using SparseArrays, SharedArrays
using HDF5
using Profile
using Distributed
using Distributions
using Plots
using Printf

macro print_array(q)
    return :(Base.print_array(stdout,round.( $(esc(q)) ,sigdigits=4));println())
end

include("seif_fct.jl")

function seif()
    filename1 = "aa3_lsr2.h5"
    filename2 = "aa3_dr.h5"
    filename3 = "z.h5"
    include_dir = "./data/"

    dt = 25e-3
    graph = spzeros(Bool,1,1)
    qnoise = diagm(0 => [5, 0.02])

    speed, steering, timeit = h5open(include_dir * filename2, "r") do file
          (read(file, "/speed"), read(file, "/steering"), read(file, "/time"))
    end

    timelsr = h5open(include_dir * filename1, "r") do file
          read(file, "/timeLsr")
    end
    llsr = size(timelsr,1)

    z = h5open(include_dir * filename3, "r") do file
          read(file, "/zz")
    end

    μ = zeros(3)
    ξ = zeros(3)
    Ω = 10e4*diagm(0 => [1,1,1])
    m₀ = zeros(Int,0)
    n = 20

    stindex = 1
    j = findfirst(timelsr .>= timeit[stindex])[1]

    gr(reuse=true)
    plt = plot([0;.1],[0;.2];axis=false,
    leg=false,
    aspect_ratio=1,
    grid=false)
    xlims!((-200,200))
    ylims!((-50,300))

    poses = zeros(3,5000)
    #ξ_hist = zeros(3,500)
    t₁ = time()

    for i ∈ stindex:size(timeit,1)
        t₂ = time()

        if i > stindex && i%5000==0
            plot!(poses[1,:],poses[2,:],line = (:black, 1, 0.2))
            scatter!(μ[4:2:end],μ[5:2:end],
            markersize = 1 )
            Plots.gui()
            h5open(include_dir * "iter_$i.h5", "w") do file
                write(file, "/poses", poses)
                write(file, "/μ", μ)
                write(file, "/ξ", ξ)
                write(file, "/Ω", Ω)
                write(file, "/m₀", m₀)
                #write(file, "/ξ_hist", ξ_hist)
            end
            poses = zeros(3,5000)
            #ξ_hist = zeros(3,500)
        end

        motion(speed[i],steering[i],dt,m₀,μ,ξ,Ω,graph)

        estimate(m₀,μ,ξ,Ω)

        m₃ = zeros(Int,0)
        while j ≤ llsr && timelsr[j] < timeit[i]+dt*1000
             z₁ = z[2*j-1:2*j,z[2*j,:] .≠ 0]
             if !isempty(z₁)
                 c = correspondence(z₁,m₀,μ,Ω,graph)
                 Ω,graph = measurement(z₁,c,μ,ξ,Ω,graph)
                 m₃ = sort(union(m₃,c))
             end
             j += 1
        end
        if !isempty(m₃)
                m₄ = sort(setdiff(m₀,m₃))
                q = [m₄;unique(m₃)]
                tmp = q[max(1,length(q)-n+1):end]
                m₁ = setdiff(m₀,tmp)
                m₁ = union(m₁,setdiff(m₃,tmp))
                m₀ = tmp
                if !isempty(m₁)
                    sparsification(m₀,m₁,μ,ξ,Ω,graph)
                end
        end

        t₃ = time()
        println("\33[2J")
        @printf "iter: %d\n" i
        @printf "iter time: %.5f\n" t₃-t₂
        @printf "avg: %.5f\n" ((t₃-t₁)/(i-stindex+1))

        poses[:,1+i%5000] = μ[1:3]
    end
end #function seif()


text:
@testset "Order" begin
    @test DNAKmer("AA") < DNAKmer("AC") < DNAKmer("AG") < DNAKmer("AT") < DNAKmer("CA")
    @test RNAKmer("AA") < RNAKmer("AC") < RNAKmer("AG") < RNAKmer("AU") < RNAKmer("CA")
end


text:
"""
    get_parametric_type(i)

Return the types `T1`,..,`Tn` of a parametric type `i::T{T1,..,Tn}`.

!!! note
    - If the parametric type has only one parameter, use `get_parametric_type(i)[1]`.
    - If no parametric type, returns an array containing the type itself.
"""
function get_parametric_type(i)
    length(collect(typeof(i).parameters)) > 0 ? collect(typeof(i).parameters) : [typeof(i)]
end


text:
module GModel

using DocStringExtensions

# TODO: Remove build (which currently prevents segfault):
#using Pkg; Pkg.build()

using Random
using Sundials # CVODE_BDF() solver for ODE
using Distributions
using LinearAlgebra
using DifferentialEquations

export run_G
export run_G_ensemble

# TODO: It would be nice to have run_G_ensemble take a pointer to the
# function G (called run_G below), which maps the parameters u to G(u),
# as an input argument. In the example below, run_G runs the model, but
# in general run_G will be supplied by the user and run the specific
# model the user wants to do UQ with.
# So, the interface would ideally be something like:
#      run_G_ensemble(params, G) where G is user-defined
"""
    GSettings{FT<:AbstractFloat, KT, D}

Structure to hold all information to run the forward model G

# Fields
$(DocStringExtensions.FIELDS)
"""
struct GSettings{FT<:AbstractFloat, KT, D}
    "a kernel tensor"
    kernel::KT
    "a model distribution"
    dist::D
    "the moments of dist that model should return"
    moments::Array{FT, 1}
    "time period over which to run the model, e.g., `(0, 1)`"
    tspan::Tuple{FT, FT}
end

"""
    run_G_ensemble(params::Array{FT, 2},
                   settings::GSettings{FT},
                   update_params,
                   moment,
                   get_src;
                   rng_seed=42) where {FT<:AbstractFloat}

Run the forward model G for an array of parameters by iteratively
calling run_G for each of the N_ensemble parameter values.
Return g_ens, an array of size N_ensemble x N_data, where
g_ens[j,:] = G(params[j,:])

 - `params` - array of size N_ensemble x N_parameters containing the
              parameters for which G will be run
 - `settings` - a GSetttings struct

"""
function run_G_ensemble(params::Array{FT, 2},
                        settings::GSettings{FT},
                        update_params,
                        moment,
                        get_src;
                        rng_seed=42) where {FT<:AbstractFloat}

    N_ens = size(params, 2) # params is N_ens x N_params
    n_moments = length(settings.moments)
    g_ens = zeros(n_moments,N_ens)

    Random.seed!(rng_seed)
    for i in 1:N_ens
        # run the model with the current parameters, i.e., map θ to G(θ)
        g_ens[:, i] = run_G(params[:, i], settings, update_params, moment, get_src)
    end

    return g_ens
end

"""
    run_G(u::Array{FT, 1},
          settings::GSettings{FT},
          update_params,
          moment,
          get_src) where {FT<:AbstractFloat}

Return g_u = G(u), a vector of length N_data.

 - `u` - parameter vector of length N_parameters
 - `settings` - a GSettings struct

"""
function run_G(u::Array{FT, 1},
               settings::GSettings{FT},
               update_params,
               moment,
               get_src) where {FT<:AbstractFloat}

    # generate the initial distribution
    dist = update_params(settings.dist, u)

    # Numerical parameters
    tol = FT(1e-7)

    # Make sure moments are up to date. mom0 is the initial condition for the
    # ODE problem
    moments = settings.moments
    moments_init = fill(NaN, length(moments))
    for (i, mom) in enumerate(moments)
        moments_init[i] = moment(dist, convert(FT, mom))
    end

    # Set up ODE problem: dM/dt = f(M,p,t)
    rhs(M, p, t) = get_src(M, dist, settings.kernel)
    prob = ODEProblem(rhs, moments_init, settings.tspan)
    # Solve the ODE
    sol = solve(prob, CVODE_BDF(), alg_hints=[:stiff], reltol=tol, abstol=tol)
    # Return moments at last time step
    moments_final = vcat(sol.u'...)[end, :]

    return moments_final
end

end # module


text:
# Convert PramaterTable to DataFrame

# Convert EnsembleParameterTable to DataFrame


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import .Base: unsafe_convert, lock, trylock, unlock, islocked, wait, notify, AbstractLock

export SpinLock

# Important Note: these low-level primitives defined here
#   are typically not for general usage

##########################################
# Atomic Locks
##########################################

"""
    SpinLock()

Create a non-reentrant, test-and-test-and-set spin lock.
Recursive use will result in a deadlock.
This kind of lock should only be used around code that takes little time
to execute and does not block (e.g. perform I/O).
In general, [`ReentrantLock`](@ref) should be used instead.

Each [`lock`](@ref) must be matched with an [`unlock`](@ref).

Test-and-test-and-set spin locks are quickest up to about 30ish
contending threads. If you have more contention than that, different
synchronization approaches should be considered.
"""
mutable struct SpinLock <: AbstractLock
    # we make this much larger than necessary to minimize false-sharing
    @atomic owned::Int
    SpinLock() = new(0)
end

# Note: this cannot assert that the lock is held by the correct thread, because we do not
# track which thread locked it. Users beware.
Base.assert_havelock(l::SpinLock) = islocked(l) ? nothing : Base.concurrency_violation()

function lock(l::SpinLock)
    while true
        if @inline trylock(l)
            return
        end
        ccall(:jl_cpu_pause, Cvoid, ())
        # Temporary solution before we have gc transition support in codegen.
        ccall(:jl_gc_safepoint, Cvoid, ())
    end
end

function trylock(l::SpinLock)
    if l.owned == 0
        GC.disable_finalizers()
        p = @atomicswap :acquire l.owned = 1
        if p == 0
            return true
        end
        GC.enable_finalizers()
    end
    return false
end

function unlock(l::SpinLock)
    if (@atomicswap :release l.owned = 0) == 0
        error("unlock count must match lock count")
    end
    GC.enable_finalizers()
    ccall(:jl_cpu_wake, Cvoid, ())
    return
end

function islocked(l::SpinLock)
    return l.owned != 0
end

max_stars_repo_name: wentasah/julia

max_issues_repo_name: wentasah/julia

max_forks_repo_name: wentasah/julia


text:
using DifferentialEquations
using Test
using MRIgeneralizedBloch
using MRIgeneralizedBloch: apply_hamiltonian_graham_superlorentzian!

##
max_error = 5e-2

## set parameters
ω1 = π / 500e-6
ω0 = 200.0
B1 = 1.0
m0s = 0.15
R1f = 0.3
R1s = 2.0
R2f = 1 / 65e-3
T2s = 10e-6
Rx = 30.0
TRF = 500e-6
alg = Tsit5()

## baseline IDE solution
m0 = [0.5 * (1 - m0s), 0.0, 0.5 * (1 - m0s), m0s, 1.0]
Graham_sol = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian!, m0, (0, TRF), (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s)), alg)

## Analytical gradients (using ApproxFun)
grad_list = [grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]

Graham_sol_grad = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian!, [m0; zeros(5 * (length(grad_list)))], (0.0, TRF), (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s, grad_list)), Tsit5())

## FD derivative wrt. m0s
dm0s = 1e-9
Graham_sol_dm0s = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, (m0s + dm0s), R1f, R2f, Rx, R1s, T2s),
    ),
    alg,
)

t = 0:1e-5:TRF
dxf = similar(t)
dyf = similar(t)
dzf = similar(t)
dzs = similar(t)
dxf_fd = similar(t)
dyf_fd = similar(t)
dzf_fd = similar(t)
dzs_fd = similar(t)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dm0s(t[i])[1] - Graham_sol(t[i])[1]) / dm0s
    dyf_fd[i] = (Graham_sol_dm0s(t[i])[2] - Graham_sol(t[i])[2]) / dm0s
    dzf_fd[i] = (Graham_sol_dm0s(t[i])[3] - Graham_sol(t[i])[3]) / dm0s
    dzs_fd[i] = (Graham_sol_dm0s(t[i])[4] - Graham_sol(t[i])[4]) / dm0s

    dxf[i] = Graham_sol_grad(t[i])[6]
    dyf[i] = Graham_sol_grad(t[i])[7]
    dzf[i] = Graham_sol_grad(t[i])[8]
    dzs[i] = Graham_sol_grad(t[i])[9]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. R1f
dR1f = 1e-9
Graham_sol_dR1f = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, m0s, (R1f + dR1f), R2f, Rx, R1s, T2s),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dR1f(t[i])[1] - Graham_sol(t[i])[1]) / dR1f
    dyf_fd[i] = (Graham_sol_dR1f(t[i])[2] - Graham_sol(t[i])[2]) / dR1f
    dzf_fd[i] = (Graham_sol_dR1f(t[i])[3] - Graham_sol(t[i])[3]) / dR1f
    dzs_fd[i] = (Graham_sol_dR1f(t[i])[4] - Graham_sol(t[i])[4]) / dR1f

    dxf[i] = Graham_sol_grad(t[i])[11]
    dyf[i] = Graham_sol_grad(t[i])[12]
    dzf[i] = Graham_sol_grad(t[i])[13]
    dzs[i] = Graham_sol_grad(t[i])[14]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error * 10

## test derivative wrt. R2f
dR2f = 1e-9
Graham_sol_dR2f = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, m0s, R1f, (R2f + dR2f), Rx, R1s, T2s),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dR2f(t[i])[1] - Graham_sol(t[i])[1]) / dR2f
    dyf_fd[i] = (Graham_sol_dR2f(t[i])[2] - Graham_sol(t[i])[2]) / dR2f
    dzf_fd[i] = (Graham_sol_dR2f(t[i])[3] - Graham_sol(t[i])[3]) / dR2f
    dzs_fd[i] = (Graham_sol_dR2f(t[i])[4] - Graham_sol(t[i])[4]) / dR2f

    dxf[i] = Graham_sol_grad(t[i])[16]
    dyf[i] = Graham_sol_grad(t[i])[17]
    dzf[i] = Graham_sol_grad(t[i])[18]
    dzs[i] = Graham_sol_grad(t[i])[19]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error

## test derivative wrt. Rx
dRx = 1e-6
Graham_sol_dRx = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, m0s, R1f, R2f, (Rx + dRx), R1s, T2s),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dRx(t[i])[1] - Graham_sol(t[i])[1]) / dRx
    dyf_fd[i] = (Graham_sol_dRx(t[i])[2] - Graham_sol(t[i])[2]) / dRx
    dzf_fd[i] = (Graham_sol_dRx(t[i])[3] - Graham_sol(t[i])[3]) / dRx
    dzs_fd[i] = (Graham_sol_dRx(t[i])[4] - Graham_sol(t[i])[4]) / dRx

    dxf[i] = Graham_sol_grad(t[i])[21]
    dyf[i] = Graham_sol_grad(t[i])[22]
    dzf[i] = Graham_sol_grad(t[i])[23]
    dzs[i] = Graham_sol_grad(t[i])[24]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. R1s
dR1s = 1e-5
Graham_sol_dR1s = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, (R1s + dR1s), T2s),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dR1s(t[i])[1] - Graham_sol(t[i])[1]) / dR1s
    dyf_fd[i] = (Graham_sol_dR1s(t[i])[2] - Graham_sol(t[i])[2]) / dR1s
    dzf_fd[i] = (Graham_sol_dR1s(t[i])[3] - Graham_sol(t[i])[3]) / dR1s
    dzs_fd[i] = (Graham_sol_dR1s(t[i])[4] - Graham_sol(t[i])[4]) / dR1s

    dxf[i] = Graham_sol_grad(t[i])[26]
    dyf[i] = Graham_sol_grad(t[i])[27]
    dzf[i] = Graham_sol_grad(t[i])[28]
    dzs[i] = Graham_sol_grad(t[i])[29]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error * 10

## test derivative wrt. T2s
dT2s = 1e-14
Graham_sol_dT2s = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, (T2s + dT2s)),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dT2s(t[i])[1] - Graham_sol(t[i])[1]) / dT2s
    dyf_fd[i] = (Graham_sol_dT2s(t[i])[2] - Graham_sol(t[i])[2]) / dT2s
    dzf_fd[i] = (Graham_sol_dT2s(t[i])[3] - Graham_sol(t[i])[3]) / dT2s
    dzs_fd[i] = (Graham_sol_dT2s(t[i])[4] - Graham_sol(t[i])[4]) / dT2s

    dxf[i] = Graham_sol_grad(t[i])[31]
    dyf[i] = Graham_sol_grad(t[i])[32]
    dzf[i] = Graham_sol_grad(t[i])[33]
    dzs[i] = Graham_sol_grad(t[i])[34]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. ω0
dω0 = 1
Graham_sol_dω0 = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, (ω0 + dω0), TRF, m0s, R1f, R2f, Rx, R1s, T2s),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dω0(t[i])[1] - Graham_sol(t[i])[1]) / dω0
    dyf_fd[i] = (Graham_sol_dω0(t[i])[2] - Graham_sol(t[i])[2]) / dω0
    dzf_fd[i] = (Graham_sol_dω0(t[i])[3] - Graham_sol(t[i])[3]) / dω0
    dzs_fd[i] = (Graham_sol_dω0(t[i])[4] - Graham_sol(t[i])[4]) / dω0

    dxf[i] = Graham_sol_grad(t[i])[36]
    dyf[i] = Graham_sol_grad(t[i])[37]
    dzf[i] = Graham_sol_grad(t[i])[38]
    dzs[i] = Graham_sol_grad(t[i])[39]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error

## test derivative wrt. B1
dB1 = 1e-9
Graham_sol_dB1 = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, (B1 + dB1), ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dB1(t[i])[1] - Graham_sol(t[i])[1]) / dB1
    dyf_fd[i] = (Graham_sol_dB1(t[i])[2] - Graham_sol(t[i])[2]) / dB1
    dzf_fd[i] = (Graham_sol_dB1(t[i])[3] - Graham_sol(t[i])[3]) / dB1
    dzs_fd[i] = (Graham_sol_dB1(t[i])[4] - Graham_sol(t[i])[4]) / dB1

    dxf[i] = Graham_sol_grad(t[i])[41]
    dyf[i] = Graham_sol_grad(t[i])[42]
    dzf[i] = Graham_sol_grad(t[i])[43]
    dzs[i] = Graham_sol_grad(t[i])[44]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## ##############################################################
# R1a test
#################################################################
R1a = 0.8
# baseline IDE solution
m0 = [0.5 * (1 - m0s), 0.0, 0.5 * (1 - m0s), m0s, 1.0]
Graham_sol = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian!, m0, (0, TRF), (ω1, B1, ω0, TRF, m0s, R1a, R2f, Rx, R1a, T2s)), alg)

## Analytical gradients (using ApproxFun)
grad_list = [grad_R1a()]

Graham_sol_grad = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian!, [m0; zeros(5 * (length(grad_list)))], (0.0, TRF), (ω1, B1, ω0, TRF, m0s, R1a, R2f, Rx, R1a, T2s, grad_list)), Tsit5())

## FD derivative wrt. R1a
dR1a = 1e-6
Graham_sol_dm0s = solve(
    ODEProblem(
        apply_hamiltonian_graham_superlorentzian!,
        m0,
        (0.0, TRF),
        (ω1, B1, ω0, TRF, m0s, R1a + dR1a, R2f, Rx, R1a + dR1a, T2s),
    ),
    alg,
)

t = 0:1e-5:TRF
dxf = similar(t)
dyf = similar(t)
dzf = similar(t)
dzs = similar(t)
dxf_fd = similar(t)
dyf_fd = similar(t)
dzf_fd = similar(t)
dzs_fd = similar(t)

for i = 1:length(t)
    dxf_fd[i] = (Graham_sol_dm0s(t[i])[1] - Graham_sol(t[i])[1]) / dR1a
    dyf_fd[i] = (Graham_sol_dm0s(t[i])[2] - Graham_sol(t[i])[2]) / dR1a
    dzf_fd[i] = (Graham_sol_dm0s(t[i])[3] - Graham_sol(t[i])[3]) / dR1a
    dzs_fd[i] = (Graham_sol_dm0s(t[i])[4] - Graham_sol(t[i])[4]) / dR1a

    dxf[i] = Graham_sol_grad(t[i])[6]
    dyf[i] = Graham_sol_grad(t[i])[7]
    dzf[i] = Graham_sol_grad(t[i])[8]
    dzs[i] = Graham_sol_grad(t[i])[9]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## plot the functions for debug purposes
# using Plots
# plotlyjs()
# theme(:lime);
# plot(t, dxf, ticks=:native, label="xf", legend=:bottomleft)
# plot!(t, dxf_fd, label="xf_fd")
# plot!(t, dyf, ticks=:native, label="yf")
# plot!(t, dyf_fd, label="yf_fd")
# plot!(t, dzf, ticks=:native, label="zf")
# plot!(t, dzf_fd, label="zf_fd")
# plot!(t, dzs, label="zs")
# plot!(t, dzs_fd, label="zs_fd")


text:
function with_temp_dir(f::Function)
    original_directory = pwd()
    _temp_dir = mktempdir()
    atexit(() -> rm(_temp_dir; force = true, recursive = true))
    cd(_temp_dir)
    result = f(_temp_dir)
    cd(original_directory)
    rm(_temp_dir; force = true, recursive = true)
    return result
end


text:
export  initialiseFromExcel,
        initialiseFromDB,
        runFromExcel

"""
"""
function initialiseFromExcel( mpSim::MPsim, filename::AbstractString;
    initDB::Bool=true, showInfo::Bool=true )::Bool

    # Add extension .xlsx if necessary.
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    # Check if the file exists.
    if !ispath( filename )
        @warn string( "The file '", filename,
            "' does not exist. The simulation is incorrectly configured ",
            "and should not be relied upon." )
        return false
    end  # if !ispath( filename )

    setSimulationShowInfo!( mpSim, showInfo )
    configSource = "Excel"
    isOkay = true

    XLSX.openxlsx( filename ) do xf
        # Check if the file has a tab "General".
        if !XLSX.hassheet( xf, "General" )
            @warn "Excel file has no tab 'General', cannot perform basic simulation setup. The simulation is incorrectly configured and should not be relied upon."
            isOkay = false
            return
        end  # if !XLSX.hassheet( xf, "General" )

        # Read database parameters.
        filePath = filename[1:(end-5)]
        configSource = readDBpars( mpSim, xf["General"], filePath )

        # Don't read Excel parameters if the source of the configuration is a
        #   database.
        if configSource != "Excel"
            return
        end  # if !isDBnew

        # Read general parameters.
        readGeneralPars( mpSim, xf["General"], filePath )

        XLSX.openxlsx( mpSim.catFileName ) do catXF
            readAttritionSchemes( mpSim, catXF )
            readAttributes( mpSim, xf, catXF )
            readBaseNodes( mpSim, xf, catXF )
            readCompoundNodes( mpSim, xf, catXF )
        end  # XLSX.openxlsx( mpSim.catFileName ) do catXF

        readTransitions( mpSim, xf )
        readRetirement( mpSim, xf )
        readInitPop( mpSim, xf, filePath )
    end  # XLSX.openxlsx( filename ) do xf

    return isOkay

end  # initialiseFromExcel( mpSim, filename, initDB, showInfo )

"""
```
initialiseFromExcel(
    filename::String;
    initDB::Bool=true,
    showOutput::Bool=true )
```
"""
function initialiseFromExcel( filename::String; initDB::Bool=true,
    showInfo::Bool=true )::Tuple{MPsim, Bool}

    mpSim = MPsim()
    isOkay = initialiseFromExcel( mpSim, filename, initDB=initDB,
        showInfo=showInfo )

    return (mpSim, isOkay)

end  # initialiseFromExcel( filename, initDB, showInfo )

"""
```
initialiseFromDB(
    mpSim::MPsim,
    filename::String;
    configName::String="config",
    showInfo::Bool=true )
```
"""
function initialiseFromDB( mpSim::MPsim, filename::String;
    configName::String="config", showInfo::Bool=true )::Bool

    # Add extension .sqlite if necessary.
    filename = string( filename, endswith( filename, ".sqlite" ) ? "" :
        ".sqlite" )

    # Check if the file exists.
    if !ispath( filename )
        @warn string( "The file '", filename,
            "' does not exist. The simulation is incorrectly configured ",
            "and should not be relied upon." )
        return false
    end  # if !ispath( filename )

    configDB = SQLite.DB( filename )

    # Check if the database has the required table.
    if configName ∉ SQLite.tables( configDB )[:name]
        @warn string( "The database '", filename, "' does not have a table '",
            configName, "'. The simulation is incorrectly configured ",
            "and should not be relied upon." )
        return false
    end  # if configName ∉ SQLite.tables( configDB )[:name]

    setSimulationShowInfo!( mpSim, showInfo )
    isOkay = true

    readGeneralPars( mpSim, configDB, configName )
    readAttritionSchemes( mpSim, configDB, configName )
    readAttributes( mpSim, configDB, configName )
    readBaseNodes( mpSim, configDB, configName )
    readCompoundNodes( mpSim, configDB, configName )
    readRecruitment( mpSim, configDB, configName )
    readTransitions( mpSim, configDB, configName )
    readRetirement( mpSim, configDB, configName )

    return isOkay

end  # initialiseFromDB( mpSim, filename, configName, showInfo )

"""
```
initialiseFromDB(
    filename::String;
    configName::String="config",
    showInfo::Bool=true )
```
"""
function initialiseFromDB( filename::String; configName::String="config",
    showInfo::Bool=true )::Tuple{MPsim, Bool}

    mpSim = MPsim()
    isOkay = initialiseFromDB( mpSim, filename, configName=configName,
        showInfo=showInfo )
    return (mpSim, isOkay)

end  # initialiseFromDB( filename, configName, showInfo )

function runFromExcel( mpSim::MPsim, filename::String;
    initDB::Bool=true, showInfo::Bool=true )::Bool

    isOkay = initialiseFromExcel( mpSim, filename, initDB=initDB,
        showInfo=showInfo )

    if !isOkay
        return false
    end  # if !isOkay

    # Add extension .xlsx if necessary.
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )
    return runSim( mpSim, filename, showInfo )

end  # runFromExcel( mpSim, filename, initDB, showInfo )

function runFromExcel( filename; initDB::Bool=true,
    showInfo::Bool=true )::Tuple{MPsim, Bool}

    mpSim = MPsim()
    isOkay = runFromExcel( mpSim, filename, initDB=initDB, showInfo=showInfo )
    return (mpSim, isOkay)

end  # runFromExcel( fileName, initDB, showInfo )


text:
## Initialize
import Pkg
Pkg.activate("/global/homes/a/aramreye/.julia/environments/v1.5/")
using Plots, BenchmarkTools, Profile
gr()
Pkg.activate("/global/homes/a/aramreye/RamirezReyes_Yang_2020_SpontaneousCyclogenesis/Project.toml")
using AvailablePotentialEnergyFramework
using Statistics:mean
using NCDatasets
import ImageSegmentation: seeded_region_growing, labels_map, region_adjacency_graph

## Read test file
test_file_with_cyclones = "/global/cscratch1/sd/aramreye/sam3d/subsetsForLin/f5e-4_2km_1000km_control_2d.nc"
test_file_with_cyclones_3d = "/global/cscratch1/sd/aramreye/sam3d/subsetsForLin/f5e-4_2km_1000km_control_3d.nc"
ds = Dataset(test_file_with_cyclones);
ds_3d = Dataset(test_file_with_cyclones_3d);
pres = ds["PSFC"][:,:,end- 10:end] :: Array{Float32,3};
U = ds_3d["U"][:,:,:,end - 10:end] :: Array{Float32,4};
V = ds_3d["V"][:,:,:, end - 10: end] :: Array{Float32,4};
x = ds_3d["x"][:] :: Array{Float32,1};
z = ds_3d["z"][:] :: Array{Float32,1};
close(ds)
close(ds_3d)
windspeed = hypot.(U,V);
pres_anomaly = pres[:,:,:] .- mean(pres[:,:,:],dims=(1,2));

function speedtest1()
  test_file_with_cyclones_3d = "/global/cscratch1/sd/aramreye/sam3d/subsetsForLin/f5e-4_2km_1000km_control_3d.nc"
  ds_3d = Dataset(test_file_with_cyclones_3d);
  U = ds_3d["U"][:,:,:,end - 10:end] :: Array{Float32,4};
  close(ds_3d)
end

function speedtest2()
  test_file_with_cyclones_3d = "/global/cscratch1/sd/aramreye/sam3d/subsetsForLin/f5e-4_2km_1000km_control_3d.nc"
  ds_3d = Dataset(test_file_with_cyclones_3d);
  U = variable(ds_3d,"U")[:,:,:,end - 10:end] :: Array{Float32,4};
  close(ds_3d)
end

@btime speedtest1()

@btime speedtest2()

## Test center detection
function test_centerdetection(pres_anomaly,windspeed)
    anim = @animate for t in 1:size(windspeed,4)
        centers = @views AvailablePotentialEnergyFramework.findcyclonecenters_aspressureminima(pres_anomaly[:,:,t],-5,2000);
        centers_x= [centers[i][1] for i in 1:length(centers)];
        centers_y = [centers[i][2] for i in 1:length(centers)];
        p1 = @views heatmap(permutedims(windspeed[:,:,1,t]),title="Wind speed (m/s)",c=:viridis)
        p2 = @views heatmap(permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa) and centers , timestep = $t")
        scatter!(p2,centers_x,centers_y,markersize=10,markershape=:cross,c=:red)
        plot(p1,p2,size=(800,380))
    end
    gif(anim, "testresults/center_detection.gif", fps = 1)
end

test_centerdetection(pres_anomaly,windspeed)

## Test cyclone detection and segmentation
function test_cyclonesegmentation(pressure,pres_anomaly)
    anim = @animate for t in 1:size(windspeed,4)
        centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,t],-5,2000)
        labels_map2 = copy(labels_map(cyclones));
        labels_map2[labels_map2 .== 1000] .= 0;
        p1 = @views heatmap(permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
        p2 = heatmap(permutedims(labels_map2),reuse=false,title="Segmented cyclones timestep = $t")
        plot(p1,p2,size=(800,380))
    end
    gif(anim, "testresults/cyclone_segmentation.gif", fps = 1)
end

test_cyclonesegmentation(pres,pres_anomaly)

#Test cyclone averaging
function test_cycloneaddition_withmask!(presaddition,windspeedaddition,buf2d,buf3d,pres,windspeed,pres_anomaly)
    totalcyclonecount = 0
    anim = @animate for t in 1:size(windspeed,4)
        centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,t],-5,2000)
        if !isnothing(centers_and_labels[1])
            count2d = AvailablePotentialEnergyFramework.add_allcyclones!(presaddition,buf2d,pres[:,:,t],cyclones,centers_and_labels)
            count3d = AvailablePotentialEnergyFramework.add_allcyclones!(windspeedaddition,buf3d,windspeed[:,:,:,t],cyclones,centers_and_labels)
            totalcyclonecount += count2d
            p1 = @views heatmap(1e-3x,1e-3x,permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
            p2 = @views heatmap(1e-3x,1e-3x,permutedims(presaddition[:,:])./totalcyclonecount,title="Averaged surface pressure (hPa)",clims=(950,1010))
            p3 = @views contourf(1e-3x,1e-3z,permutedims(windspeedaddition[:,256,:])./totalcyclonecount,title="Averaged wind speed (m/s)")
            plot(p1,p2,p3,size=(1200,390), layout = (1,3))
        end
    end
    gif(anim, "testresults/cycloneaveraging_withmask.gif", fps = 1)
end

size_pres = size(pres)[1:2]
size_windspeed = size(windspeed)[1:3]
presaddition = zeros(size_pres)
windspeedaddition = zeros(size_windspeed)
buf2d = similar(presaddition)
buf3d = similar(windspeedaddition)
test_cycloneaddition_withmask!(presaddition,windspeedaddition,buf2d,buf3d,pres,windspeed,pres_anomaly)

#Test cyclone averaging
function test_cycloneaddition_withoutmask!(presaddition,windspeedaddition,buf2d,buf3d,pres,windspeed,pres_anomaly)
    totalcyclonecount = 0
    anim = @animate for t in 1:size(windspeed,4)
        centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,t],-5,2000)
        if !isnothing(centers_and_labels[1])
            count2d = AvailablePotentialEnergyFramework.add_allcyclones!(presaddition,buf2d,pres[:,:,t],cyclones,centers_and_labels,false)
            count3d = AvailablePotentialEnergyFramework.add_allcyclones!(windspeedaddition,buf3d,windspeed[:,:,:,t],cyclones,centers_and_labels,false)
            totalcyclonecount += count2d
            p1 = @views heatmap(1e-3x,1e-3x,permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
            p2 = @views heatmap(1e-3x,1e-3x,permutedims(presaddition[:,:])./totalcyclonecount,title="Averaged surface pressure (hPa)",clims=(950,1010))
            p3 = @views contourf(1e-3x,1e-3z,permutedims(windspeedaddition[:,256,:])./totalcyclonecount,title="Averaged wind speed (m/s)")
            plot(p1,p2,p3,size=(1200,390), layout = (1,3))
        end
    end
    gif(anim, "testresults/cycloneaveraging_withoutmask.gif", fps = 1)
end

size_pres = size(pres)[1:2]
size_windspeed = size(windspeed)[1:3]
presaddition = zeros(size_pres)
windspeedaddition = zeros(size_windspeed)
buf2d = similar(presaddition)
buf3d = similar(windspeedaddition)
test_cycloneaddition_withoutmask!(presaddition,windspeedaddition,buf2d,buf3d,pres,windspeed,pres_anomaly)

#Test cyclone averaging
function test_radialaverage_withoutmask!(pressureradial,windspeedradial,presaddition,windspeedaddition,buf2d,buf3d,radiusbins,pres,windspeed,pres_anomaly)
    totalcyclonecount = 0
    anim = @animate for t in 1:size(windspeed,4)
        centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,t],-5,2000)
        if !isnothing(centers_and_labels[1])
            count2d = AvailablePotentialEnergyFramework.add_allcyclones!(presaddition,buf2d,pres[:,:,t],cyclones,centers_and_labels,false)
            count3d = AvailablePotentialEnergyFramework.add_allcyclones!(windspeedaddition,buf3d,windspeed[:,:,:,t],cyclones,centers_and_labels,false)
            totalcyclonecount += count2d
            for idx in 1:(length(radiusbins) - 1)
                windspeedradial[idx,:] .= AvailablePotentialEnergyFramework.averageallindistance( (radiusbins[idx],radiusbins[idx+1]),windspeedaddition./totalcyclonecount,(256,256), 2000 )
                pressureradial[idx] = AvailablePotentialEnergyFramework.averageallindistance( (radiusbins[idx],radiusbins[idx+1]),presaddition./totalcyclonecount,(256,256), 2000 )
            end
            p1 = @views heatmap(1e-3x,1e-3x,permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
            p2 = @views plot(1e-3(radiusbins[1:end-1]),pressureradial,title="Averaged surface pressure (hPa)",clims=(950,1010))
            p3 = @views contourf(1e-3(radiusbins[1:end-1]),1e-3z,permutedims(windspeedradial),title="Averaged wind speed (m/s)")
            plot(p1,p2,p3,size=(1200,390), layout = (1,3))
        end
    end
    gif(anim, "testresults/composite_radial_withoutmask.gif", fps = 1)
end

#Test azimuthal averaging
size_pres = size(pres)[1:2]
size_windspeed = size(windspeed)[1:3]
presaddition = zeros(size_pres)
windspeedaddition = zeros(size_windspeed)
buf2d = similar(presaddition)
buf3d = similar(windspeedaddition)
radiusbins = 1000:2000:512000;
presradial = zeros(length(radiusbins) - 1)
windspeedradial = zeros(length(radiusbins) - 1,size(windspeed,3))
test_radialaverage_withoutmask!(presradial,windspeedradial,presaddition,windspeedaddition,buf2d,buf3d,radiusbins,pres,windspeed,pres_anomaly)

#Test cyclone averaging
function test_radialaverage_withmask!(pressureradial,windspeedradial,presaddition,windspeedaddition,buf2d,buf3d,radiusbins,pres,windspeed,pres_anomaly)
    totalcyclonecount = 0
    anim = @animate for t in 1:size(windspeed,4)
    #for t in 1:size(windspeed,4)
        centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,t],-5,2000)
        if !isnothing(centers_and_labels[1])
            #@info centers_and_labels[1]
            count2d = AvailablePotentialEnergyFramework.add_allcyclones!(presaddition,buf2d,pres[:,:,t],cyclones,centers_and_labels,true)
            count3d = AvailablePotentialEnergyFramework.add_allcyclones!(windspeedaddition,buf3d,windspeed[:,:,:,t],cyclones,centers_and_labels,true)
            totalcyclonecount += count2d
            for idx in 1:(length(radiusbins) - 1)
                if !iszero(totalcyclonecount)
                    windspeedradial[idx,:] .= AvailablePotentialEnergyFramework.averageallindistance( (radiusbins[idx],radiusbins[idx+1]),windspeedaddition./totalcyclonecount,(256,256), 2000 )
                    pressureradial[idx] = AvailablePotentialEnergyFramework.averageallindistance( (radiusbins[idx],radiusbins[idx+1]),presaddition./totalcyclonecount,(256,256), 2000 )
                end
            end
            p1 = @views heatmap(1e-3x,1e-3x,permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
            p2 = @views plot(1e-3.*(radiusbins[1:end-1]),pressureradial,title="Averaged surface pressure (hPa)",clims=(950,1010))
            p3 = @views contourf(1e-3.*(radiusbins[1:end-1]),1e-3z,permutedims(windspeedradial),title="Averaged wind speed (m/s)")
            plot(p1,p2,p3,size=(1200,390), layout = (1,3))
        end
    end
    gif(anim, "testresults/composite_radial_withmask.gif", fps = 1)
end

#Test azimuthal averaging
size_pres = size(pres)[1:2]
size_windspeed = size(windspeed)[1:3]
presaddition = zeros(size_pres)
windspeedaddition = zeros(size_windspeed)
buf2d = similar(presaddition)
buf3d = similar(windspeedaddition)
radiusbins = 1000:2000:512000;
presradial = ones(length(radiusbins) - 1)
windspeedradial = ones(length(radiusbins) - 1,size(windspeed,3))
test_radialaverage_withmask!(presradial,windspeedradial,presaddition,windspeedaddition,buf2d,buf3d,radiusbins,pres,windspeed,pres_anomaly)


text:
@safetestset "isascii" begin
    using Validator

    @test Validator.isascii("foobar") === true
    @test Validator.isascii("0987654321") === true
    @test Validator.isascii("test@example.com") === true
    @test Validator.isascii("1234abcDEF") === true

    @test Validator.isascii("ｆｏｏbar") === false
    @test Validator.isascii("ｘｙｚ０９８") === false
    @test Validator.isascii("１２３456") === false
    @test Validator.isascii("ｶﾀｶﾅ") === false
end


text:
@defcomp impactforests begin
    regions = Index()

    forests = Variable(index=[time,regions])

    forbm = Parameter(index=[regions])
    gdp90 = Parameter(index=[regions])
    pop90 = Parameter(index=[regions])

    acco2 = Parameter(index=[time])

    income = Parameter(index=[time,regions])
    population = Parameter(index=[time,regions])
    temp = Parameter(index=[time,regions])

    forel = Parameter(default = -0.31)
    fornl = Parameter(default = 1)
    co2pre = Parameter()
    forco2 = Parameter(default = 2.2972)

    function run_timestep(p, v, d, t)

        if !is_first(t)
            for r in d.regions
                ypc = 1000.0 * p.income[t, r] / p.population[t, r]
                ypc90 = p.gdp90[r] / p.pop90[r] * 1000.0

                # TODO -oDavid Anthoff: RT uses -lP.forel for ppp case
                v.forests[t, r] = p.forbm[r] * p.income[t, r] * (ypc / ypc90)^p.forel * (0.5 * p.temp[t, r]^p.fornl + 0.5 * log(p.acco2[t - 1] / p.co2pre) * p.forco2)

                if v.forests[t, r] > 0.1 * p.income[t, r]
                    v.forests[t, r] = 0.1 * p.income[t, r]
                end
            end
        end
    end
end


text:
function Error(experimental_data_array,simulation_time_array,simulation_state_array,output_index,species_symbol,data_dictionary)

  # initialize -
  total_error_value = Inf

  # Interpole the simulation data to the experimental time scale -
  interpolated_simulated_trajectory = np.interp(experimental_data_array[:,1],simulation_time_array,simulation_state_array[:,output_index]);

  # Compute the error -
  # Error is the sum of a shape function, a scale function and a background_copy_number function -

  # Compute the shape term -
  scaled_data_array = alpha_scaling(experimental_data_array[:,2]);
  scaled_simulation_array = alpha_scaling(interpolated_simulated_trajectory);
  shape_term = sum((scaled_data_array - scaled_simulation_array).^2);

  # Compute the initial value term -
  background_copy_number_dictionary = data_dictionary["background_copy_number_dictionary"];
  background_copy_number = background_copy_number_dictionary[species_symbol]
  cell_volume = data_dictionary["cell_volume"];
  av_number = data_dictionary["av_number"];
  initial_concentration = background_copy_number*(1/av_number)*(1/cell_volume)*1e9; # nM
  initial_value_term = ((interpolated_simulated_trajectory[1] - initial_concentration)^2);

  # Compute the scale term -
  simulated_fold_change = (interpolated_simulated_trajectory)./(interpolated_simulated_trajectory[1]);
  scale_term = sum((simulated_fold_change - experimental_data_array[:,2]).^2);

  # @show (interpolated_simulated_trajectory[1],initial_concentration,interpolated_simulated_trajectory[end],species_symbol)
  # @show (shape_term,scale_term,initial_value_term)

  # compute the total error -
  total_error_value = shape_term+scale_term+initial_value_term;

  # return -
  return total_error_value
end

function Error_E2F(experimental_data_array,simulation_time_array,simulation_state_array,output_index,species_symbol,data_dictionary)

  # initialize -
  total_error_value = Inf

  # Interpole the simulation data to the experimental time scale -
  interpolated_simulated_trajectory = np.interp(experimental_data_array[:,1],simulation_time_array,simulation_state_array[:,output_index]);

  # Compute G1G0 -
  simulated_G1G0 = 1./interpolated_simulated_trajectory;

  # Compute the shape term -
  scaled_data_array = alpha_scaling(experimental_data_array[:,2]);
  scaled_simulation_array = alpha_scaling(simulated_G1G0);
  shape_term = sum((scaled_data_array - scaled_simulation_array).^2);

  # Compute the initial value term -
  background_copy_number_dictionary = data_dictionary["background_copy_number_dictionary"];
  background_copy_number = background_copy_number_dictionary[species_symbol]
  cell_volume = data_dictionary["cell_volume"];
  av_number = data_dictionary["av_number"];
  initial_concentration = background_copy_number*(1/av_number)*(1/cell_volume)*1e9; # nM
  initial_value_term = ((interpolated_simulated_trajectory[1] - initial_concentration)^2);

  # Compute the fold change term -
  scaling_factor = simulated_G1G0[1];
  simulated_fold_change = (simulated_G1G0)./(scaling_factor);
  scale_term = sum((simulated_fold_change - experimental_data_array[:,2]).^2);

  # compute the total error -
  total_error_value = shape_term+scale_term+initial_value_term;

  # return -
  return total_error_value;
end

function alpha_scaling(data_array)

  # get the dimension -
  number_of_points = length(data_array)

  # what is the min index?
  min_index = indmin(data_array)
  min_value = data_array[min_index]

  # What is the max index?
  max_index = indmax(data_array)
  max_value = data_array[max_index]

  # compute the scaling -
  scaled_array = Float64[]
  for index = 1:number_of_points

    value = data_array[index]
    scaled_value = (value - min_value)/(max_value - min_value)
    push!(scaled_array,scaled_value)
  end

  return scaled_array
end


text:
# Internal function to exchange to references in the repository

function _swap(repo::Repo{T},posa::Int,posb::Int) where T
    ida=repo.ids[posa]
    idb=repo.ids[posb]
    repo.pos[ida]=posb
    repo.pos[idb]=posa
    repo.ids[posa]=idb
    repo.ids[posb]=ida
end

function _print_class(repo::Repo{T},c::Int,max::Int) where T
    N=repo.counts[c]
    maxh=fld1(max,2)
    if N<=max || maxh<1
        for i=1:N
            id=repo_id(repo,c,i)
            println("    $i - $id: $(repo[id])")
        end
    else
        for i=1:maxh
            id=repo_id(repo,c,i)
            println("    $i - $id: $(repo[id])")
        end
        println("    ...")
        for i=(N-maxh+1):N
            id=repo_id(repo,c,i)
            println("    $i - $id: $(repo[id])")
        end
    end
end

function _print_class(io::IO,repo::Repo{T},c::Int,max::Int) where T
    N=repo.counts[c]
    maxh=fld1(max,2)
    if N<=max || maxh<1
        for i=1:N
            id=repo_id(repo,c,i)
            println(io,"    $i - $id: $(repo[id])")
        end
    else
        for i=1:maxh
            id=repo_id(repo,c,i)
            println(io,"    $i - $id: $(repo[id])")
        end
        println("    ...")
        for i=(N-maxh+1):N
            id=repo_id(repo,c,i)
            println(io,"    $i - $id: $(repo[id])")
        end
    end
end


text:
# Geometric.
geo(p::Float64) = trace(:flip, Bernoulli(p)) ? 0 : 1 + trace(:geo, geo, p)

@testset "Trace" begin
    ret, cl = simulate(geo, 0.5)
    @test haskey(cl.trace, :flip)
end

@testset "Generate" begin
    sel = target([(:flip, ) => true])
    ret, cl, w = generate(sel, geo, 0.5)
    @test haskey(cl.trace, :flip)
    @test cl[:flip] == true
end

@testset "Update" begin
    ret, cl = simulate(geo, 0.5)
    sel = target([(:flip, ) => true])
    ret, cl, w, retdiff, d = update(sel, cl)
    @test cl[:flip] == true
end


text:
function parse(::Type{Time}, s::String, ::Type{IGCTimeFormat})
    return Time(s, IGC_TIME_FMT)
end

function parse(::Type{Date}, s::String, ::Type{IGCDateFormat})
    m = match(DATE_PATTERN, s)
    day = parse(Int, m[1])
    month = parse(Int, m[2])
    year = parse(Int, m[3])
    year = 2000 + year
    return Date(year, month, day)
end

function parse(::Type{IGCLatitude}, s::String)
    d = parse(Int64, s[1:2])
    m = parse(Float64, s[3:7]) / 1000.
    ordinal = s[8]

    latitude = d + m / 60.

    if !(0. <= latitude <= 90.)
        throw(ArgumentError("Invalid latitude format $(repr(s)) - range error"))
    end

    if ordinal ∉ ('S', 'N')
        throw(ArgumentError("Invalid latitude format $(repr(s)) - ordinal error"))
    end

    if ordinal == 'S'
        latitude = -latitude
    end

    return IGCLatitude(latitude)
end

function parse(::Type{IGCLongitude}, s::String)
    d = parse(Float64, s[1:3])
    m = parse(Float64, s[4:8]) / 1000.
    ordinal = s[9]

    longitude = d + m / 60.

    if !(0. <= longitude <= 180.)
        throw(ArgumentError("Invalid longitude format $(repr(s)) - range error"))
    end

    if ordinal ∉ ('W', 'E')
        throw(ArgumentError("Invalid longitude format $(repr(s)) - ordinal error"))
    end

    if ordinal == 'W'
        longitude = -longitude
    end

    return IGCLongitude(longitude)
end

function parse(::Type{IGCPressureAltitude}, s::String)
    return IGCPressureAltitude(parse(Int, s))
end

function parse(::Type{IGCGpsAltitude}, s::String)
    return IGCGpsAltitude(parse(Int, s))
end

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

function parse(::Type{A_record}, line::String)
    @assert line[1] == 'A'
    id_addition = length(line) == 7 ? "" : strip(line[8:end])
    return A_record(line[2:4], line[5:7], id_addition)
end

function parse(::Type{B_record}, line::String)
    @assert line[1] == 'B'

    return B_record(
        parse(Time, line[2:7], IGCTimeFormat),
        parse(IGCLatitude, line[8:15]),
        parse(IGCLongitude, line[16:24]),
        FixValidity.parse(line[25]),
        parse(IGCPressureAltitude, line[26:30]),
        parse(IGCGpsAltitude, line[31:35]),
        36,
        strip(line[36:end])
    )
end

function parse(::Type{Abstract_C_record}, line::String)
    @assert line[1] == 'C'

    istaskinfo = isdigit(line[9])

    if istaskinfo
        return C_record_task_info(
            parse(Date, line[2:7], IGCDateFormat),
            parse(Time, line[8:13], IGCTimeFormat),
            parse(Date, line[14:19], IGCDateFormat),
            line[20:23],
            parse(Int64, line[24:25]),
            strip(line[26:end])
        )
    else
        return C_record_waypoint_info(
            parse(IGCLatitude, line[2:9]),
            parse(IGCLongitude, line[10:18]),
            strip(line[19:end])
        )
    end
end

function parse(::Type{D_record}, line::String)
    @assert line[1] == 'D'
    qualifier = GpsQualifier.parse(line[2])
    station_id = line[3:6]
    return D_record(
        qualifier,
        station_id
    )
end

function parse(::Type{E_record}, line::String)
    @assert line[1] == 'E'
    time = parse(Time, line[2:7], IGCTimeFormat)
    tlc = line[8:10]
    extension_string = strip(line[11:end])
    return E_record(
        time,
        tlc,
        extension_string
    )
end

function parse(::Type{F_record}, line::String)
    @assert line[1] == 'F'

    time = parse(Time, line[2:7], IGCTimeFormat)

    # each satellite ID should have two digits
    if (length(strip(line)) - 7) % 2 != 0
        throw(ArgumentError("F record formatting is incorrect"))
    end

    satellites = []
    nb_satellites = floor(Int, (length(strip(line)) - 7) / 2)

    starting_byte = 8
    for satelite_index in 1:nb_satellites
        push!(satellites, line[starting_byte:starting_byte + 1])
        starting_byte += 2
    end

    return F_record(
        time,  # time
        satellites  # satellites
    )
end

function parse(::Type{G_record}, line::String)
    @assert line[1] == 'G'

    security_code = strip(line[2:end])

    return G_record(
        security_code
    )
end

function parse(::Type{IGCExtension}, line::String)
    nb_extensions = parse(Int64, line[2:3])

    if nb_extensions * 7 + 3 != length(strip(line))
        throw(ArgumentError("I record contains incorrect number of digits"))
    end

    extensions = IGCExtension[]

    for extension_index in 0:nb_extensions - 1
        extension_str = line[extension_index * 7 + 4:(extension_index + 1) * 7 + 3]
        start_byte = parse(Int64, extension_str[1:2])
        end_byte = parse(Int64, extension_str[3:4])
        tlc = extension_str[5:7]

        extension = IGCExtension(start_byte:end_byte, tlc)
        push!(extensions, extension)
    end

    return extensions
end

function parse(::Type{I_record}, line::String)
    @assert line[1] == 'I'

    line = string(strip(line))

    extensions = parse(IGCExtension, line)

    return I_record(extensions)
end

function parse(::Type{J_record}, line::String)
    @assert line[1] == 'J'

    line = string(strip(line))

    extensions = parse(IGCExtension, line)

    return J_record(extensions)
end

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

function parse(::Type{H_record_FiXAccuracy}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.FXA(.*)"i
    m = match(pattern, line)
    fix_accuracy = m[1]
    if length(fix_accuracy) == 0
        fix_accuracy = missing
        return H_record_FiXAccuracy(source, fix_accuracy)
    else
        fix_accuracy = parse(Int, fix_accuracy)
        return H_record_FiXAccuracy(source, fix_accuracy)
    end
end

function parse(::Type{H_record_DaTE}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.DTE(\d\d\d\d\d\d)"i
    m = match(pattern, line)
    date_str = string(m[1])
    date = parse(Date, date_str, IGCDateFormat)
    return H_record_DaTE(source, date)
end

function parse(::Type{H_record_PiLoT}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.PLT.*: ?(.*)"i
    m = match(pattern, line)
    name = m[1]
    return H_record_PiLoT(source, name)
end

function parse(::Type{H_record_Copilot}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.CM2CREW2: ?(.*)"i
    m = match(pattern, line)
    name = m[1]
    return H_record_Copilot(source, name)
end

function parse(::Type{H_record_GliderType}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.GTYGLIDERTYPE: ?(.*)"i
    m = match(pattern, line)
    return H_record_GliderType(source, m[1])
end

function parse(::Type{H_record_GliderRegistration}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.GIDGLIDERID: ?(.*)"i
    m = match(pattern, line)
    return H_record_GliderRegistration(source, m[1])
end

function parse(::Type{H_record_GpsDatum}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.DTM100GPSDATUM: ?(.*)"i
    m = match(pattern, line)
    return H_record_GpsDatum(source, m[1])
end

function parse(::Type{H_record_FirmwareRevision}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.RFWFIRMWAREVERSION: ?(.*)"i
    m = match(pattern, line)
    return H_record_FirmwareRevision(source, m[1])
end

function parse(::Type{H_record_HardwareRevision}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.RHWHARDWAREVERSION: ?(.*)"i
    m = match(pattern, line)
    return H_record_HardwareRevision(source, m[1])
end

function parse(::Type{H_record_ManufacturerModel}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.FTYFRTYPE: ?(.*),(.*)"i
    m = match(pattern, line)
    manufacturer = string(strip(m[1]))
    model = strip(m[2])
    return H_record_ManufacturerModel(source, manufacturer, model)
end

function parse(::Type{H_record_GpsReceiver}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.GPS: ?(.*)"i
    m = match(pattern, line)
    s_gps_receiver = m[1]
    t = split(s_gps_receiver, ",")

    if length(t) == 1
        manufacturer = t[1]
        model = ""
        channels = ""
        max_alt = ""
    elseif length(t) == 2
        manufacturer, model = t
        channels = ""
        max_alt = ""
    elseif length(t) == 3
        manufacturer, channels, max_alt = t
        model = ""
    elseif length(t) == 4
        manufacturer, model, channels, max_alt = t
    else
        throw(ArgumentError("Invalid GpsReceiver $line"))
    end

    manufacturer = strip(manufacturer)
    model = strip(model)
    channels = strip(channels)
    max_alt = strip(max_alt)

    if length(channels) > 0
        pattern = r"(\d*)(ch)?"i
        m = match(pattern, channels)
        channels = parse(Int, m[1])
    else
        channels = missing
    end

    if length(max_alt) > 0
        pattern = r"(max)?(\d*)(\D*)"i
        m = match(pattern, max_alt)
        max_alt_value = parse(Int, strip(m[2]))

        max_alt_unit = m[3]
        if !ismissing(max_alt_value) && length(max_alt_unit) == 0
            max_alt_unit = "m"
        end
    else
        max_alt_value = missing
        max_alt_unit = ""
    end

    max_alt = MaxAlt(max_alt_value, max_alt_unit)

    return H_record_GpsReceiver(source,
        manufacturer,
        model,
        channels,
        max_alt
    )
end

function parse(::Type{H_record_PressureAltitudeSensor}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.PRSPRESSALTSENSOR: ?(.*)"i
    m = match(pattern, line)
    s = m[1]
    t = split(s, ",")

    if length(t) == 1
        manufacturer = t[1]
        model = ""
        max_alt = ""
    elseif length(t) == 2
        manufacturer, model = t
        max_alt = ""
    elseif length(t) == 3
        manufacturer, model, max_alt = t
    else
        throw(ArgumentError("Invalid PressureAltitudeSensor $line"))
    end

    manufacturer = string(strip(manufacturer))
    model = string(strip(model))
    max_alt = strip(max_alt)

    if length(max_alt) > 0
        pattern = r"(max)?(\d*)(\D*)"i
        m = match(pattern, max_alt)
        max_alt_value = parse(Int, strip(m[2]))

        max_alt_unit = m[3]
        if !ismissing(max_alt_value) && length(max_alt_unit) == 0
            max_alt_unit = "m"
        end
    else
        max_alt_value = missing
        max_alt_unit = ""
    end

    max_alt = MaxAlt(max_alt_value, max_alt_unit)

    return H_record_PressureAltitudeSensor(source,
        manufacturer,
        model,
        max_alt
    )
end

function parse(::Type{H_record_CompetitionId}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.CIDCOMPETITIONID: ?(.*)"i
    m = match(pattern, line)
    s = strip(m[1])
    return H_record_CompetitionId(source, s)
end

function parse(::Type{H_record_CompetitionClass}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.CCLCOMPETITIONCLASS: ?(.*)"i
    m = match(pattern, line)
    s = strip(m[1])
    return H_record_CompetitionClass(source, s)
end

function parse(::Type{H_record_TimeZoneOffset}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.TZNTIMEZONE: ?(.*)"i
    m = match(pattern, line)
    h = string(strip(m[1]))
    if occursin(".", h)
        h, m = split(h, ".")
        h = parse(Int64, h)
        m = rpad(m, 2, "0")
        m = sign(h) * parse(Int64, m)
    else
        h = parse(Int64, h)
        m = 0
    end
    minutes = round(Int, h * 60 + m)
    return H_record_TimeZoneOffset(minutes, source=source)
end

function parse(::Type{H_record_MeansOfPropulsion}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.MOPSENSOR: ?(.*)"i
    m = match(pattern, line)
    s = m[1]
    return H_record_MeansOfPropulsion(source, s)
end

function parse(::Type{H_record_Site}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.SITSite: ?(.*)"i
    m = match(pattern, line)
    s = m[1]
    return H_record_Site(source, s)
end

function parse(::Type{H_record_UnitsOfMeasure}, line::String)
    line = strip(line)
    source = line[2]
    pattern = r"H.UNTUnits: ?(.*)"i
    m = match(pattern, line)
    s = m[1]
    return H_record_UnitsOfMeasure(source, split(s, ","))
end

const D_TLC_H_RECORD = Dict(
    "DTE" => H_record_DaTE,
    "FXA" => H_record_FiXAccuracy,
    "PLT" => H_record_PiLoT,
    "CM2" => H_record_Copilot,
    "GTY" => H_record_GliderType,
    "GID" => H_record_GliderRegistration,
    "DTM" => H_record_GpsDatum,
    "RFW" => H_record_FirmwareRevision,
    "RHW" => H_record_HardwareRevision,
    "FTY" => H_record_ManufacturerModel,
    "GPS" => H_record_GpsReceiver,
    "PRS" => H_record_PressureAltitudeSensor,
    "CID" => H_record_CompetitionId,
    "CCL" => H_record_CompetitionClass,
    "TZN" => H_record_TimeZoneOffset,
    "MOP" => H_record_MeansOfPropulsion,
    "SIT" => H_record_Site,
    "TZO" => H_record_TimeZoneOffset,
    "UNT" => H_record_UnitsOfMeasure
)

function parse(::Type{Abstract_H_record}, line::String)
    @assert line[1] == 'H'

    # three letter code
    tlc = line[3:5]

    try
        T_H_record = D_TLC_H_RECORD[tlc]
        value = parse(T_H_record, line)
        return value
    catch e
        if isa(e, KeyError)
            throw(ArgumentError("Invalid h-record '$(strip(line))'"))
        else
            throw(e)
        end
    end
end

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

function parse(::Type{K_record}, line::String)
    @assert line[1] == 'K'
    line = string(strip(line))
    time = parse(Time, line[2:7], IGCTimeFormat)
    value_string = line[8:end]
    start_index = 7
    return K_record(
        time,
        value_string,
        start_index
    )
end

function parse(::Type{L_record}, line::String)
    @assert line[1] == 'L'
    line = strip(line)
    source = line[2:4]
    comment = line[5:end]
    return L_record(
        source,
        comment
    )
end

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

function parse(::Type{Vector{Abstract_IGC_record}}, s, sep=EOL_DEFAULT)
    return parse.(Abstract_IGC_record, string.(split(s, sep)))
end

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

function update(record::B_record, fix_record_extensions::AbstractExtensionsRecord)
    # process/update record with fix record extensions

    i = record.start_index_extensions
    ext = record.extensions_string

    for extension in fix_record_extensions.extensions
        start_byte, end_byte = extension.bytes
        start_byte = start_byte - i - 1
        end_byte = end_byte - i - 1

        update(record, extension)

    end

    new_record = B_record(
        record.time,
        record.latitude,
        record.longitude,
        record.validity,
        record.pressure_alt,
        record.gps_alt,
        -1,  # start_index_extensions
        "",  # extensions_string
        IGCExtension[]
    )

    return new_record
end

function update(record::B_record, extension::IGCExtension)
    # todo
    #record.update(
    #    {extension['extension_type']: int(ext[start_byte:end_byte + 1])}
    #)
end

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

function read(fname::AbstractString, ::Type{IGCDocument}; parsing_mode=ParsingMode.DEFAULT, store_all_records=DEFAULT_STORE_ALL_RECORDS)
    stream = open(fname)
    return parse(IGCDocument, stream, parsing_mode=parsing_mode, store_all_records=store_all_records)
end

function parse(::Type{IGCDocument}, s; parsing_mode=ParsingMode.DEFAULT, store_all_records=DEFAULT_STORE_ALL_RECORDS)
    stream = IOBuffer(s)
    return parse(IGCDocument, stream::IO, parsing_mode=parsing_mode, store_all_records=store_all_records)
end

function parse(::Type{Abstract_IGC_record}, line::String)
    line = string(line)
    record_char = uppercase(line[1])
    T_record = D_RECORD_CHAR[record_char]
    return parse(T_record, line)
end

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

function update!(igcdoc::IGCDocument, rec::A_record)
    fieldname = :logger_id
    if isempty(getfield(igcdoc, fieldname))
        setfield!(igcdoc, fieldname, rec)
        write(igcdoc.stream, rec)
    else
        throw(IGCNonUniqueRecordException("Non-unique $(typeof(rec)) prev: $(getfield(igcdoc, fieldname)) new: $rec"))
    end
end

#=
function update!(igcdoc::IGCDocument, rec::B_record)
    push!(igcdoc.fix_records, rec)
end

function update!(igcdoc::IGCDocument, rec::G_record)
    push!(igcdoc.security_records, rec)
end
=#

function update!(igcdoc::IGCDocument, rec::T) where {T <: Abstract_H_record}
    fieldname = D_HEADER_RECORD_FIELD[T]
    if igcdoc.parsing_mode == ParsingMode.STRICT && !ismissing(getfield(igcdoc.header, fieldname))
        msg = "Non-unique $(typeof(rec)) prev: $(getfield(igcdoc.header, fieldname)) new: $rec"
        exc = IGCParseException(msg)
        throw(exc)
    end
    setfield!(igcdoc.header, fieldname, rec)
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

function update!(igcdoc::IGCDocument, rec::H_record_DaTE)
    fieldname = :date
    if igcdoc.parsing_mode == ParsingMode.STRICT && !ismissing(getfield(igcdoc.header, fieldname))
        msg = "Non-unique $(typeof(rec)) prev: $(getfield(igcdoc.header, fieldname)) new: $rec"
        exc = IGCParseException(msg)
        throw(exc)
    end
    #dt = ZonedDateTime(DateTime(rec.value), TZ)
    #igcdoc.dt_first = dt
    #igcdoc.dt_last = dt
    setfield!(igcdoc.header, fieldname, rec)
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

function update!(igcdoc::IGCDocument, rec::C_record_task_info)
    task = Task(rec)
    # push!(igcdoc.tasks, task)
    igcdoc.task = task
    #dt = ZonedDateTime(task.info.declaration, TZ)
    #if dt > igcdoc.dt_last
    #    igcdoc.dt_last = dt
    #end
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

function update!(igcdoc::IGCDocument, rec::C_record_waypoint_info)
    #task = igcdoc.tasks[end]
    task = igcdoc.task
    push!(task.waypoints, rec)
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

function update!(igcdoc::IGCDocument, records::Vector{C_record_waypoint_info})
    for rec in records
        update!(igcdoc, rec)
    end
end

# timestamped record
function update!(igcdoc::IGCDocument, rec::T) where {T <: Union{B_record, E_record, F_record, K_record}}
    fieldname = D_RECORD_FIELD[T]
    #current_date = Date(DateTime(igcdoc.dt_last))
    #last_time = Time(DateTime(igcdoc.dt_last))
    #if rec.time < last_time
    #    dt = ZonedDateTime(current_date + Dates.day(1) + rec.time, TZ)
    #else
    #    dt = ZonedDateTime(current_date + rec.time, TZ)
    #end
    push!(getfield(igcdoc, fieldname), rec)
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
    #println("$dt $(igcdoc.dt_last)")
    #igcdoc.dt_last = dt
end

function update!(igcdoc::IGCDocument, rec::T) where {T <: Union{D_record, G_record, L_record}}
    fieldname = D_RECORD_FIELD[T]
    push!(getfield(igcdoc, fieldname), rec)
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

function update!(igcdoc::IGCDocument, rec::T) where {T <: Union{I_record, J_record}}
    fieldname = D_RECORD_FIELD[T]
    setfield!(igcdoc, fieldname, rec)
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

#=
function update!(igcdoc::IGCDocument, rec::I_record)
    igcdoc.fix_record_extensions = rec
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end

function update!(igcdoc::IGCDocument, rec::J_record)
    igcdoc.k_record_extensions = rec
    write(igcdoc.stream, igcdoc.eol)
    write(igcdoc.stream, rec)
end
=#

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

function update!(igcdoc::IGCDocument, records::Vector{T}) where {T <: Abstract_IGC_record}
    for rec in records
        update!(igcdoc, rec)
    end
end

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

function parse(::Type{IGCDocument}, stream::IO; parsing_mode=ParsingMode.STRICT, store_all_records=DEFAULT_STORE_ALL_RECORDS)
    igcdoc = IGCDocument(parsing_mode=parsing_mode)
    for (line_id, line) in enumerate(eachline(stream))
        try
            rec = parse(Abstract_IGC_record, line)
            if store_all_records
                push!(igcdoc.all_records, rec)
            end
            update!(igcdoc, rec)
        catch e
            exc = IGCParseException(line_id, line, e)
            if parsing_mode == ParsingMode.UNSTRICT
                push!(igcdoc.errors, exc)
            else
                throw(exc)
            end
        end
    end
    return igcdoc
end


text:
function parseInput(filename)
    values = readlines(filename*".txt")
    values = split.(values, "")
    return values
end

function rightNeighbor(i,m)
    if (i+1) % m == 0
        return i+1
    else return (i+1) % m
    end
end

function downNeighbor(i,n)
    if (i+1)%n == 0
        return i+1
    else
        return (i+1) % n
    end
end

function solve(filename)
    grid = parseInput(filename)
    no_changes = Inf
    no_steps = 0
    n = length(grid)
    m = length(grid[1])

    while no_changes > 0
        no_steps += 1
        no_changes = 0
        newGrid = Vector()
        for (i,line) in enumerate(grid)
            newLine = Vector{String}()
            for k in 1:m
                newLine = push!(newLine,line[k])
            end
            for (j,char) in enumerate(line)
                if char == ">" && line[rightNeighbor(j,m)] == "."
                    newLine[rightNeighbor(j,m)] = ">"
                    newLine[j] = "."
                    no_changes += 1
                end
            end
            newGrid = push!(newGrid, newLine)
        end
        grid = newGrid

        newGrid = Vector()
        for line in grid
            newLine = Vector{String}()
            for char in line
                newLine = push!(newLine, char)
            end
            newGrid = push!(newGrid, newLine)
        end

        for (i,line) in enumerate(grid)
            for (j,char) in enumerate(line)
                if char == "v" && grid[downNeighbor(i,n)][j] == "."
                    newGrid[downNeighbor(i,n)][j] = "v"
                    newGrid[i][j] = "."
                    no_changes += 1
                end
            end
        end
        grid = newGrid
    end
    println("Stopped after "*string(no_steps)*" steps.")
end

solve("input")


text:
module DateParsing

## RFC3339
include("tryparse.jl")
include("fastparse.jl")

end # module