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

import Base.Checked: add_with_overflow, mul_with_overflow

## string to integer functions ##

"""
    parse(type, str; base)

Parse a string as a number. For `Integer` types, a base can be specified
(the default is 10). For floating-point types, the string is parsed as a decimal
floating-point number.  `Complex` types are parsed from decimal strings
of the form `"R±Iim"` as a `Complex(R,I)` of the requested type; `"i"` or `"j"` can also be
used instead of `"im"`, and `"R"` or `"Iim"` are also permitted.
If the string does not contain a valid number, an error is raised.

```jldoctest
julia> parse(Int, "1234")
1234

julia> parse(Int, "1234", base = 5)
194

julia> parse(Int, "afc", base = 16)
2812

julia> parse(Float64, "1.2e-3")
0.0012

julia> parse(Complex{Float64}, "3.2e-1 + 4.5im")
0.32 + 4.5im
```
"""
parse(T::Type, str; base = Int)

function parse(::Type{T}, c::AbstractChar; base::Integer = 10) where T<:Integer
    a::Int = (base <= 36 ? 10 : 36)
    2 <= base <= 62 || throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
    d = '0' <= c <= '9' ? c-'0'    :
        'A' <= c <= 'Z' ? c-'A'+10 :
        'a' <= c <= 'z' ? c-'a'+a  : throw(ArgumentError("invalid digit: $(repr(c))"))
    d < base || throw(ArgumentError("invalid base $base digit $(repr(c))"))
    convert(T, d)
end

function parseint_next(s::AbstractString, startpos::Int, endpos::Int)
    (0 < startpos <= endpos) || (return Char(0), 0, 0)
    j = startpos
    c, startpos = next(s,startpos)
    c, startpos, j
end

function parseint_preamble(signed::Bool, base::Int, s::AbstractString, startpos::Int, endpos::Int)
    c, i, j = parseint_next(s, startpos, endpos)

    while isspace(c)
        c, i, j = parseint_next(s,i,endpos)
    end
    (j == 0) && (return 0, 0, 0)

    sgn = 1
    if signed
        if c == '-' || c == '+'
            (c == '-') && (sgn = -1)
            c, i, j = parseint_next(s,i,endpos)
        end
    end

    while isspace(c)
        c, i, j = parseint_next(s,i,endpos)
    end
    (j == 0) && (return 0, 0, 0)

    if base == 0
        if c == '0' && !done(s,i)
            c, i = next(s,i)
            base = c=='b' ? 2 : c=='o' ? 8 : c=='x' ? 16 : 10
            if base != 10
                c, i, j = parseint_next(s,i,endpos)
            end
        else
            base = 10
        end
    end
    return sgn, base, j
end

function tryparse_internal(::Type{T}, s::AbstractString, startpos::Int, endpos::Int, base_::Integer, raise::Bool) where T<:Integer
    sgn, base, i = parseint_preamble(T<:Signed, Int(base_), s, startpos, endpos)
    if sgn == 0 && base == 0 && i == 0
        raise && throw(ArgumentError("input string is empty or only contains whitespace"))
        return nothing
    end
    if !(2 <= base <= 62)
        raise && throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
        return nothing
    end
    if i == 0
        raise && throw(ArgumentError("premature end of integer: $(repr(SubString(s,startpos,endpos)))"))
        return nothing
    end
    c, i = parseint_next(s,i,endpos)
    if i == 0
        raise && throw(ArgumentError("premature end of integer: $(repr(SubString(s,startpos,endpos)))"))
        return nothing
    end

    base = convert(T,base)
    m::T = div(typemax(T)-base+1,base)
    n::T = 0
    a::Int = base <= 36 ? 10 : 36
    while n <= m
        d::T = '0' <= c <= '9' ? c-'0'    :
               'A' <= c <= 'Z' ? c-'A'+10 :
               'a' <= c <= 'z' ? c-'a'+a  : base
        if d >= base
            raise && throw(ArgumentError("invalid base $base digit $(repr(c)) in $(repr(SubString(s,startpos,endpos)))"))
            return nothing
        end
        n *= base
        n += d
        if i > endpos
            n *= sgn
            return n
        end
        c, i = next(s,i)
        isspace(c) && break
    end
    (T <: Signed) && (n *= sgn)
    while !isspace(c)
        d::T = '0' <= c <= '9' ? c-'0'    :
        'A' <= c <= 'Z' ? c-'A'+10 :
            'a' <= c <= 'z' ? c-'a'+a  : base
        if d >= base
            raise && throw(ArgumentError("invalid base $base digit $(repr(c)) in $(repr(SubString(s,startpos,endpos)))"))
            return nothing
        end
        (T <: Signed) && (d *= sgn)

        n, ov_mul = mul_with_overflow(n, base)
        n, ov_add = add_with_overflow(n, d)
        if ov_mul | ov_add
            raise && throw(OverflowError("overflow parsing $(repr(SubString(s,startpos,endpos)))"))
            return nothing
        end
        (i > endpos) && return n
        c, i = next(s,i)
    end
    while i <= endpos
        c, i = next(s,i)
        if !isspace(c)
            raise && throw(ArgumentError("extra characters after whitespace in $(repr(SubString(s,startpos,endpos)))"))
            return nothing
        end
    end
    return n
end

function tryparse_internal(::Type{Bool}, sbuff::Union{String,SubString{String}},
        startpos::Int, endpos::Int, base::Integer, raise::Bool)
    if isempty(sbuff)
        raise && throw(ArgumentError("input string is empty"))
        return nothing
    end

    orig_start = startpos
    orig_end   = endpos

    # Ignore leading and trailing whitespace
    while isspace(sbuff[startpos]) && startpos <= endpos
        startpos = nextind(sbuff, startpos)
    end
    while isspace(sbuff[endpos]) && endpos >= startpos
        endpos = prevind(sbuff, endpos)
    end

    len = endpos - startpos + 1
    p   = pointer(sbuff) + startpos - 1
    GC.@preserve sbuff begin
        (len == 4) && (0 == ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
                                  p, "true", 4)) && (return true)
        (len == 5) && (0 == ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
                                  p, "false", 5)) && (return false)
    end

    if raise
        substr = SubString(sbuff, orig_start, orig_end) # show input string in the error to avoid confusion
        if all(isspace, substr)
            throw(ArgumentError("input string only contains whitespace"))
        else
            throw(ArgumentError("invalid Bool representation: $(repr(substr))"))
        end
    end
    return nothing
end

@inline function check_valid_base(base)
    if 2 <= base <= 62
        return base
    end
    throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
end

"""
    tryparse(type, str; base)

Like [`parse`](@ref), but returns either a value of the requested type,
or [`nothing`](@ref) if the string does not contain a valid number.
"""
function tryparse(::Type{T}, s::AbstractString; base::Union{Nothing,Integer} = nothing) where {T<:Integer}
    # Zero base means, "figure it out"
    tryparse_internal(T, s, firstindex(s), lastindex(s), base===nothing ? 0 : check_valid_base(base), false)
end

function parse(::Type{T}, s::AbstractString; base::Union{Nothing,Integer} = nothing) where {T<:Integer}
    convert(T, tryparse_internal(T, s, firstindex(s), lastindex(s),
                                 base===nothing ? 0 : check_valid_base(base), true))
end

## string to float functions ##

function tryparse(::Type{Float64}, s::String)
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
                          (Ptr{UInt8},Csize_t,Csize_t), s, 0, sizeof(s))
    hasvalue ? val : nothing
end
function tryparse(::Type{Float64}, s::SubString{String})
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
                          (Ptr{UInt8},Csize_t,Csize_t), s.string, s.offset, s.ncodeunits)
    hasvalue ? val : nothing
end
function tryparse_internal(::Type{Float64}, s::String, startpos::Int, endpos::Int)
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
                          (Ptr{UInt8},Csize_t,Csize_t), s, startpos-1, endpos-startpos+1)
    hasvalue ? val : nothing
end
function tryparse_internal(::Type{Float64}, s::SubString{String}, startpos::Int, endpos::Int)
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
                          (Ptr{UInt8},Csize_t,Csize_t), s.string, s.offset+startpos-1, endpos-startpos+1)
    hasvalue ? val : nothing
end
function tryparse(::Type{Float32}, s::String)
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
                          (Ptr{UInt8},Csize_t,Csize_t), s, 0, sizeof(s))
    hasvalue ? val : nothing
end
function tryparse(::Type{Float32}, s::SubString{String})
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
                          (Ptr{UInt8},Csize_t,Csize_t), s.string, s.offset, s.ncodeunits)
    hasvalue ? val : nothing
end
function tryparse_internal(::Type{Float32}, s::String, startpos::Int, endpos::Int)
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
                          (Ptr{UInt8},Csize_t,Csize_t), s, startpos-1, endpos-startpos+1)
    hasvalue ? val : nothing
end
function tryparse_internal(::Type{Float32}, s::SubString{String}, startpos::Int, endpos::Int)
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
                          (Ptr{UInt8},Csize_t,Csize_t), s.string, s.offset+startpos-1, endpos-startpos+1)
    hasvalue ? val : nothing
end
tryparse(::Type{T}, s::AbstractString) where {T<:Union{Float32,Float64}} = tryparse(T, String(s))
tryparse(::Type{Float16}, s::AbstractString) =
    convert(Union{Float16, Nothing}, tryparse(Float32, s))
tryparse_internal(::Type{Float16}, s::AbstractString, startpos::Int, endpos::Int) =
    convert(Union{Float16, Nothing}, tryparse_internal(Float32, s, startpos, endpos))

## string to complex functions ##

function tryparse_internal(::Type{Complex{T}}, s::Union{String,SubString{String}}, i::Int, e::Int, raise::Bool) where {T<:Real}
    # skip initial whitespace
    while i ≤ e && isspace(s[i])
        i = nextind(s, i)
    end
    if i > e
        raise && throw(ArgumentError("input string is empty or only contains whitespace"))
        return nothing
    end

    # find index of ± separating real/imaginary parts (if any)
    i₊ = coalesce(findnext(in(('+','-')), s, i), 0)
    if i₊ == i # leading ± sign
        i₊ = coalesce(findnext(in(('+','-')), s, i₊+1), 0)
    end
    if i₊ != 0 && s[i₊-1] in ('e','E') # exponent sign
        i₊ = coalesce(findnext(in(('+','-')), s, i₊+1), 0)
    end

    # find trailing im/i/j
    iᵢ = coalesce(findprev(in(('m','i','j')), s, e), 0)
    if iᵢ > 0 && s[iᵢ] == 'm' # im
        iᵢ -= 1
        if s[iᵢ] != 'i'
            raise && throw(ArgumentError("expected trailing \"im\", found only \"m\""))
            return nothing
        end
    end

    if i₊ == 0 # purely real or imaginary value
        if iᵢ > 0 # purely imaginary
            x = tryparse_internal(T, s, i, iᵢ-1, raise)
            x === nothing && return nothing
            return Complex{T}(zero(x),x)
        else # purely real
            x = tryparse_internal(T, s, i, e, raise)
            x === nothing && return nothing
            return Complex{T}(x)
        end
    end

    if iᵢ < i₊
        raise && throw(ArgumentError("missing imaginary unit"))
        return nothing # no imaginary part
    end

    # parse real part
    re = tryparse_internal(T, s, i, i₊-1, raise)
    re === nothing && return nothing

    # parse imaginary part
    im = tryparse_internal(T, s, i₊+1, iᵢ-1, raise)
    im === nothing && return nothing

    return Complex{T}(re, s[i₊]=='-' ? -im : im)
end

# the ±1 indexing above for ascii chars is specific to String, so convert:
tryparse_internal(T::Type{Complex{S}}, s::AbstractString, i::Int, e::Int, raise::Bool) where S<:Real =
    tryparse_internal(T, String(s), i, e, raise)

# fallback methods for tryparse_internal
tryparse_internal(::Type{T}, s::AbstractString, startpos::Int, endpos::Int) where T<:Real =
    startpos == firstindex(s) && endpos == lastindex(s) ? tryparse(T, s) : tryparse(T, SubString(s, startpos, endpos))
function tryparse_internal(::Type{T}, s::AbstractString, startpos::Int, endpos::Int, raise::Bool) where T<:Real
    result = tryparse_internal(T, s, startpos, endpos)
    if raise && result === nothing
        throw(ArgumentError("cannot parse $(repr(s[startpos:endpos])) as $T"))
    end
    return result
end
tryparse_internal(::Type{T}, s::AbstractString, startpos::Int, endpos::Int, raise::Bool) where T<:Integer =
    tryparse_internal(T, s, startpos, endpos, 10, raise)

parse(::Type{T}, s::AbstractString) where T<:Union{Real,Complex} =
    convert(T, tryparse_internal(T, s, firstindex(s), lastindex(s), true))

tryparse(T::Type{Complex{S}}, s::AbstractString) where S<:Real =
    tryparse_internal(T, s, firstindex(s), lastindex(s), false)

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using KilliFeeder
using Test

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

ext: jl


text:
using ..ChemistryFeaturization.AbstractType: AbstractCodec
import ..ChemistryFeaturization.encode
import ..ChemistryFeaturization.decode

"""
    DirectCodec

Codec type whose encoding function is some constant * the identity function.
"""
struct DirectCodec <: AbstractCodec
    scale::Number
end

encode(dc::DirectCodec, val) = dc.scale .* val
decode(dc::DirectCodec, encoded) = encoded ./ dc.scale


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule cilkrts_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("cilkrts")
JLLWrappers.@generate_main_file("cilkrts", UUID("71772805-00bc-5a29-9044-a26d819b7806"))
end  # module cilkrts_jll


text:
using POMDPs
using DESPOT
using POMDPToolbox

include("rockSample.jl")
include("rockSampleBounds.jl")

function main(;
                grid_size::Int64            = 4,
                num_rocks::Int64            = 4,
                n_reps::Int64               = 1,
                n_particles::Int64          = 500, # for both solver and belief updater
                main_seed::Int64            = 42,
                discount::Float64           = 0.95,
                search_depth::Int64         = 90,
                time_per_move::Float64      = 1.0,
                pruning_constant::Float64   = 0.0,
                eta::Float64                = 0.95,
                sim_len::Int64              = -1,
                max_trials::Int64           = -1,
                approximate_ubound::Bool    = false,
                debug::Int64                = 0
                )

    total_sim_steps::Int64                  = 0
    total_discounted_return::Float64        = 0.
    total_undiscounted_return::Float64      = 0.
    total_run_time::Float64                 = 0.

    # Optional parameters can be adjusted, as shown below.
    # Performance tip: control use of computational resources either by
    # limiting time_per_move, by limiting the number of trials per move, or both.
    # Setting either parameter to 0 or a negative number disables that limit.

    search_depth::Int64 = 90 #default: 90
    time_per_move::Float64 = -1.0 # sec, default: 1, unlimited: -1
    pruning_constant::Float64 = 0.0
    eta::Float64 = 0.95 # default: 0.95
    sim_len::Int64 = -1 # default: -1
    max_trials::Int64 = 100 # default: -1
    approximate_ubound::Bool = false
    tiny::Float64 = 1e-6
    debug::Int64 = 0

    for i in 1:n_reps
        @printf("\n\n\n\n================= Run %d =================\n", i)
        sim_steps,
        discounted_return,
        undiscounted_return,
        run_time = execute(
                    grid_size = grid_size,
                    num_rocks = num_rocks,
                    n_particles = n_particles,
                    main_seed = main_seed,
                    discount = discount,
                    search_depth = search_depth,
                    time_per_move = time_per_move,
                    pruning_constant = pruning_constant,
                    eta = eta,
                    sim_len = sim_len,
                    max_trials = max_trials,
                    approximate_ubound = approximate_ubound,
                    debug = debug
                    )

        total_sim_steps               += sim_steps
        total_discounted_return       += discounted_return
        total_undiscounted_return     += undiscounted_return
        total_run_time                += run_time
    end

    avg_sim_steps               = total_sim_steps/n_reps
    avg_discounted_return       = total_discounted_return/n_reps
    avg_undiscounted_return     = total_undiscounted_return/n_reps
    avg_run_time                = total_run_time/n_reps

    if (n_reps > 1)
        @printf("\n================= Batch Averages =================\n")
        @printf("Number of steps = %d\n", avg_sim_steps)
        @printf("Discounted return = %.2f\n", avg_discounted_return)
        @printf("Undiscounted return = %.2f\n", avg_undiscounted_return)
        @printf("Runtime = %.2f sec\n", avg_run_time)
    end
    return avg_sim_steps, avg_discounted_return, avg_undiscounted_return, avg_run_time
end

function execute(;
                grid_size::Int64            = 4,
                num_rocks::Int64            = 4,
                n_particles::Int64          = 500, # for both solver and belief updater
                main_seed::Int64            = 42,
                discount::Float64           = 0.95,
                search_depth::Int64         = 90,
                time_per_move::Float64      = 1.0,
                pruning_constant::Float64   = 0.0,
                eta::Float64                = 0.95,
                sim_len::Int64              = -1,
                max_trials::Int64           = -1,
                approximate_ubound::Bool    = false,
                debug::Int64                = 0
                )

    rand_max::Int64 = 2^31-1 # 2147483647

    # generate unique random seeds (optional, if not supplied, default values will be used)
    seed  ::UInt32   = convert(UInt32, main_seed)  # the main random seed that's used to set the other seeds
    w_seed::UInt32   = seed ⊻  n_particles      # world seed, used in the overall simulation
    b_seed::UInt32   = seed ⊻ (n_particles + 1) # belief seed, used for belief particle sampling, among other things
    m_seed::UInt32   = seed ⊻ (n_particles + 2) # model seed, used to initialize the problem model

    pomdp   = RockSample(
                        grid_size,
                        num_rocks,
                        rand_max = rand_max,      # optional, default: 2^31-1
                        belief_seed = b_seed,     # optional, default: 479
                        model_seed  = m_seed,     # optional, default: 476
                        discount    = discount)   # optional, default: 0.95

    # construct a belief updater and specify some of the optional keyword parameters
    bu = DESPOTBeliefUpdater{RockSampleState,
                             RockSampleAction,
                             RockSampleObs}(
                             pomdp::POMDP,
                             seed = seed,
                             rand_max = rand_max,
                             n_particles = n_particles)

    # create initial belief and allocate an updated belief object
    initial_states = initial_state_distribution(pomdp)
    current_belief = initialize_belief(bu, initial_states)
    updated_belief = create_belief(bu)
    custom_bounds = RockSampleBounds(pomdp)

    solver::DESPOTSolver = DESPOTSolver{RockSampleState,
                                        RockSampleAction,
                                        RockSampleObs,
                                        RockSampleBounds,
                                        RandomStreams
                                       }(
                               # specify the optional keyword parameters
                               bounds = custom_bounds,
                               search_depth = search_depth,
                               main_seed = seed, # specify the main random seed
                               time_per_move = time_per_move,
                               n_particles = n_particles,
                               pruning_constant = pruning_constant,
                               eta = eta,
                               sim_len = sim_len,
                               approximate_ubound = approximate_ubound,
                               max_trials =  max_trials,
                               rand_max = rand_max,
                               debug = debug)

    state::RockSampleState       = start_state(pomdp)
    next_state::RockSampleState  = RockSampleState()
    obs::RockSampleObs   = RockSampleObs()
    rewards::Array{Float64}      = Array{Float64}(0)

    rng::DESPOTDefaultRNG = DESPOTDefaultRNG(w_seed, rand_max) # used to advance the state of the simulation (world)
    policy::DESPOTPolicy = POMDPs.solve(solver, pomdp)

    sim_steps::Int64 = 0
    r::Float64 = 0.0
    println("\nSTARTING STATE: $state")
    show_state(pomdp, state) #TODO: wrap RockSample in a module
    tic() # start the clock
    while !isterminal(pomdp, state) &&
        (solver.config.sim_len == -1 || sim_steps < solver.config.sim_len)
        println("\n*************** STEP $(sim_steps+1) ***************")
        action = POMDPs.action(policy, current_belief)
        transition_distribution = POMDPs.transition(pomdp, state, action)
        next_state = POMDPs.rand(rng, transition_distribution, next_state) # update state to next state
        observation_distribution = POMDPs.observation(pomdp, state, action, next_state)
        observation_distribution.debug = 1 #TODO: remove -debug
        obs = POMDPs.rand(rng, observation_distribution)
        r = POMDPs.reward(pomdp, state, action)
        push!(rewards, r)
        state = next_state
        updated_belief = POMDPs.update(bu, current_belief, action, obs)
        current_belief = deepcopy(updated_belief) #TODO: perhaps this could be done better
        println("Action = $action")
        println("State = $next_state"); show_state(pomdp, next_state) #TODO: change once abstract types are introduced
        print(  "Observation = "); show_obs(pomdp, obs) #TODO: change once abstract types are introduced
        println("Reward = $r")
        sim_steps += 1
    end
    run_time::Float64 = toq() # stop the clock

    # Compute discounted reward
    discounted_reward::Float64 = 0.0
    multiplier::Float64 = 1.0
    for r in rewards
        discounted_reward += multiplier * r
        multiplier *= pomdp.discount
    end

    println("\n********** EXECUTION SUMMARY **********")
    @printf("Number of steps = %d\n", sim_steps)
    @printf("Undiscounted return = %.2f\n", sum(rewards))
    @printf("Discounted return = %.2f\n", discounted_reward)
    @printf("Runtime = %.2f sec\n", run_time)

    return sim_steps, sum(rewards), discounted_reward, run_time
end


text:
__precompile__()
module DataProcessingHierarchyTools
using ProgressMeter
using Glob
using MAT
using StrTables, LaTeX_Entities
using StableHashes
import StableHashes.shash
import Base:filter,show, convert

include("types.jl")

const def = LaTeX_Entities.default

function git_annex()
    cmd = nothing
    try
        cmd = readchomp(`which git-annex`)
    catch ee
    end
    return cmd
end

const levels = ["subjects", "subject", "day", "session", "array", "channel", "cell"]
const level_patterns = [r"[0-9A-Za-z]*", r"[0-9]{8}", r"session[0-9]{2}", r"array[0-9]{2}", r"channel[0-9]{3}", r"cell[0-9]{2}"]
const level_patterns_s = ["*", "*", "[0-9]*", "session[0-9]*", "array[0-9]*", "channel[0-9]*", "cell[0-9]*"]

#TODO: Find a way to (automatically) track dependency betweeen types. We could also just do this manually, i.e. have a depends_on function that lists types that a given type depends (directly) on. E.g. depends_on(Raster) = PSTH. If a type does not depend on any other types, just return an empty list
depends_on(::Type{T}) where T <: DPHData = DataType[]

"""
Returns a list of files on which `args` depends.
"""
dependencies(args::T) where T <: DPHDataArgs = String[]

function check_args(X1::T, X2::T) where T <: DPHDataArgs
    matches = true
    for f in fieldnames(T)
        x1 = getfield(X1,f)
        x2 = getfield(X2,f)
        if typeof(x1) <: AbstractVector
            if length(x1) != length(x2)
                matches = false
                break
            end
            for (_x1,_x2) in zip(x1,x2)
                if _x1 != _x2
                    matches = false
                    break
                end
            end
            if !matches
                break
            end
        else
            if !(x1 ≈ x2)
                matches = false
                break
            end
        end
    end
    matches
end

function check_args(X::T, args...) where T <: DPHDataArgs
    matches = true
    for (a0,a1) in zip(fieldnames(T), args)
        x = getfield(X, a0)
        if typeof(x) <: AbstractVector
            if length(x) != length(a1)
                matches = false
                break
            end
            for (x1,x2) in zip(x,a1)
                if x1 != x2
                    matches = false
                    break
                end
            end
            if !matches
                break
            end
        else
            if !(x ≈ a1)
                matches = false
                break
            end
        end
    end
    matches
end

function get_numbers(ss::String)
    filter(isdigit,ss)
end

shortnames = Dict("subjects" => x->"",
                  "subject" => x->x[1:1],
                  "day" => x->x,
                  "session" => x->"s$(get_numbers(x))",
                  "array" => x->"a$(get_numbers(x))",
                  "channel" => x->"g$(get_numbers(x))",
                  "cell" => x->"c$(get_numbers(x))")

function get_shortname(ss::String)
    this_level = level(ss)
    this_idx = findfirst(l->this_level==l, levels)
    _r, _p = splitdir(ss)
    qs = String[]
    while !isempty(_p) && this_idx > 1  # Stop once we have reached the bottom level
        push!(qs, shortnames[level(_p)](_p))
        _r, _p = splitdir(_r)
        this_idx -= 1 # Keep track of where we are in the hierarchy
    end
    join(reverse(qs),"")
end

"""
Returns the full name of the current level
"""
function get_fullname(ss=pwd())
    this_level = level(ss)
    this_idx = findfirst(l->this_level==l, levels)
    pp = [get_level_name(this_level, ss)]
    while this_idx > 2
        this_idx -= 1
        this_level = levels[this_idx]
        push!(pp, get_level_name(this_level, ss))
    end
    joinpath(reverse(pp)...)
end

level() = level(pwd())
level(::Type{DPHData}) = error("Not implemented")
filename(::Type{DPHData}) = error("Not implemented")
datatype(::Type{DPHDataArgs}) = error("Not implemented")
datatype(X::T) where T <: DPHDataArgs = datatype(T)
version(X::DPHDataArgs) = "UNKNOWN"

function filename(args::T) where T <: DPHDataArgs
    fname = filename(datatype(T))
    h = string(shash(args),base=16)
    bn, ext = splitext(fname)
    fname = join([bn, "_", h, ext])
    fname
end

function filename(args::Vector{T}) where T <: DPHDataArgs
    fname = filename(datatype(T))
    h = string(shash(args),base=16)
    bn, ext = splitext(fname)
    fname = join([bn, "_", h, ext])
    fname
end

matname(::Type{DPHData}) = error("Not implemented")

function load(args::T) where T <: DPHDataArgs
    fname = filename(args)
    load(datatype(T), fname)
end

"""
Returns `true` if the data described by `args` has already been
computed
"""
function computed(args::T) where T <: DPHDataArgs
    fname = filename(args)
    return isfile(fname) || islink(fname)
end

function load(args::Vector{T}) where T <: DPHDataArgs
    fname = filename(args)
    if isfile(fname)
        return load(Vector{datatype(T)}, fname)
    end
    error("No data exist with the specified arguments")
end

function plot_data(::Type{T},fig, args::T2, plotargs::T3) where T <: DPHData where T2 <: DPHDataArgs where T3 <: DPHPlotArgs
    error("Not implemented")
end

export DPHData, level, filename, plot_data, datatype, BootstrappedDataArgs

"""
Get the level of the directory represented by `cwd`.
"""
function level(cwd::String)
    numbers = map(x->first(string(x)), 0:9)
    dd = last(splitdir(cwd))
    ss = string(rstrip(dd, numbers))
    if isempty(ss)
        # only numbers; assume this is a date
        ss = "day"
    elseif dd == ss
        #no numbers, this is the subject name
        ss = "subject"
    end
    return ss
end

"""
Get the path of `dir` up `target_level`
"""
function get_level_path(target_level::String, dir=pwd())
    parts = splitpath(dir)
    new_parts = String[]
    target_idx = findfirst(levels .== target_level)
    for p in parts
        this_idx = findfirst(levels .== level(p))
        if this_idx <= target_idx
            push!(new_parts, p)
        end
    end
    return joinpath(new_parts...)
end

"""
Get the name of the requested level
"""
function get_level_name(target_level::String, dir=pwd())
    this_level = level(dir)
    this_idx = findfirst(l->this_level==l, levels)
    target_idx = findfirst(l->target_level==l, levels)
    i = this_idx
    cwd = dir
    pp = ""
    while i >= target_idx
        cwd, pp = splitdir(cwd)
        i -= 1
    end
    pp
end

"""
Get all directories corresponding to `target_level` under the current hierarchy
"""
function get_level_dirs(target_level::String, dir=pwd())
    dirs = cd(dir) do
        this_level = level()
        this_idx = findfirst(l->this_level==l, levels)
        target_idx = findfirst(l->target_level==l, levels)
        if target_idx == this_idx
            dirs = ["."]
        elseif target_idx < this_idx
            rel_path = process_level(target_level, dir)
            dirs = glob(joinpath(rel_path, "..", level_patterns_s[target_idx]))
        else
            dirs = glob(joinpath(level_patterns_s[this_idx+1:target_idx]...))
        end
    end
    if dir != pwd()
        dirs = [joinpath(dir, d) for d in dirs]
    end
    dirs
end

"""
Get all unique level directories contained in `dirs`.
"""
function get_level_dirs(level::String, dirs::Vector{String})
    cwd = pwd()
    level_dirs = String[]
    for c in dirs
        cd(c) do
            cd(process_level(level))
            _dir = pwd()
            _dir = strip(replace(_dir, cwd => ""), '/')
            if !(_dir in level_dirs)
                push!(level_dirs, _dir)
            end
        end
    end
    level_dirs
end

"""
Returns the relative path to an object of type `T`, using `dir` as the starting point.
"""
function process_level(::Type{T}, dir=pwd();kvs...) where T <: DPHData
    target_level = level(T)
    process_level(target_level, dir;kvs...)
end

"""
Returns the path relative to `dir` of the level `target_level`.
"""
function process_level(target_level::String, dir=pwd();kvs...)
    # get the current level
    this_level = level(dir)
    this_idx = findfirst(l->this_level==l, levels)
    target_idx = findfirst(l->target_level==l, levels)
    for lidx in [this_idx, target_idx]
        if lidx == nothing || !(0 < lidx <= length(levels))
            throw(ArgumentError("Unknown level"))
        end
    end
    pl = ["."]
    append!(pl, [".." for i in 1:(this_idx - target_idx)])
    dirstring = joinpath(pl...)
end

"""
Returns the path of the current directory relative to `level_dir`

"""
function get_relative_path(level_dir::String,dir=pwd())
    this_level = level(dir)
    this_idx = findfirst(l->this_level==l, levels)
    target_idx = findfirst(l->level_dir==l, levels)
    if target_idx+1 <= this_idx
        parts = String[]
        for ii in target_idx+1:this_idx
            push!(parts, get_level_name(levels[ii],dir))
        end
    elseif target_idx == this_idx
        parts = ["."]
    else
        error("Target level must be below the current level")
    end
    joinpath(parts...)
end

"""
Process each directory in `dirs`, creating an object of type `T`, and returning a concatenation of those objects.
"""
function process_dirs(::Type{T}, dirs::Vector{String}, args...;kvs...) where T <: DPHData
    pp = cd(dirs[1]) do
        T(args...;kvs...)
    end
    @showprogress 1 "Processing dirs..." for d in dirs[2:end]
        _pp = cd(d) do
            T(args...;kvs...)
        end
        pp = hcat(pp, _pp)
    end
    return pp
end

"""
Visit each directory in `dirs`, instantiating type `T` with argumments `args` and keyword arguments `kvs`. Note that this function is similar to `process_dirs`, except unlike that function, `visit_dirs` does not return any results.
"""
function visit_dirs(::Type{T}, dirs::Vector{String}, args...;kvs...) where T <: DPHData
    skipped_dirs = String[]
    @showprogress 1 "Processing dirs..." for d in dirs
        cd(d) do
            try
                T(args...;kvs...)
            catch
                push!(skipped_dirs, d)
            end
        end
    end
    skipped_dirs
end

function visit_dirs(func::Function, dirs::Vector{String}, args...;kvs...)
    skipped_dirs = String[]
    @showprogress 1 "Processing dirs..." for d in dirs
        cd(d) do
         #   try
                func(args...;kvs...)
         #   catch
         #       push!(skipped_dirs, d)
         #   end
        end
    end
    skipped_dirs
end

"""
Process each directory in `dirs` by running the function `func`.
"""
function process_dirs(func::Function, dirs::Vector{String}, args...;kvs...)
    Q = Vector{Any}(undef, length(dirs))
    @showprogress 1 "Processing dirs..." for (i,d) in enumerate(dirs)
        Q[i] = cd(d) do
            func(args...;kvs...)
        end
    end
    Q
end

"""
Load an object of type `T` from the current directory, using additional constructor arguments `args`.
"""
function load(::Type{T}, args...;kvs...) where T <: DPHData
    dir = process_level(T)
    qq = cd(dir) do
        if isfile(filename(T))
            qq = T()
        else
            qq = zero(T)
        end
        qq
    end
    qq
end

"""
Convert some unicde symbols to their latex equivalent before saving
"""
function sanitise(ss::String)
    oo = String[]
    for (i,c) in enumerate(ss)
        sc = string(c)
        m = matches(def, sc)
        if !isempty(m)
            push!(oo, m[1])
        else
            push!(oo, sc)
        end
    end
    join(oo, "")
end

function desanitise(ss::String)
    oo = String[]
    ss_split = split(ss, "_")
    for _ss in ss_split
        _nss = filter(!isdigit, _ss)
        ll = lookupname(def, _nss)
        if isempty(ll)
            ll = _nss
        end
        lln = replace(_ss, _nss => ll)
        push!(oo, lln)
    end
    join(oo, "_")
end

function save(X::T, fname=filename(X.args);overwrite=false) where T <: DPHData
    if isfile(fname) && !overwrite
        error("File $fname already exists")
    end
    Q = convert(Dict{String,Any}, X)
    MAT.matwrite(fname,Q)
end

function save(X::Vector{T}) where T <: DPHData
    fname = filename([x.args for x in X])
    Q = Dict{String, Dict{String,Any}}()
    for (i,x) in enumerate(X)
        Q["idx$(i)"] = convert(Dict{String,Any}, x)
    end
    MAT.matwrite(fname,Q)
end

function Base.convert(::Type{Dict{String,Any}}, X::T) where T <: Union{DPHData, DPHDataArgs}
    Q = Dict{String,Any}()
    for f in fieldnames(T)
        v = getfield(X, f)
        fs = string(f)
        fs = sanitise(fs)
        if typeof(v) <: AbstractVector
            Q[fs] = collect(v)
        elseif typeof(v) <: DPHDataArgs
            Q[fs] = convert(Dict{String,Any}, v)
        elseif typeof(v) <: Symbol
            Q[fs] = string(v)
        else
            Q[fs] = v
        end
    end
    Q
end

function load(::Type{T}, fname=filename(T)) where T <: DPHData
    if git_annex() != nothing
        is_annex = false
        try
            run(pipeline(`$(git_annex()) status`, stdout=devnull, stderr=devnull))
            is_annex = true
        catch
        end
        if is_annex
            run(`$(git_annex()) get $fname`)
        end
    end
    if isfile(fname)
        Q = MAT.matread(fname)
    else
        error("No data exist with the specified arguments")
    end
    convert(T, Q)
end

function load(::Type{Vector{T}}, fname=filename{T}) where T <: DPHData
    Q = MAT.matread(fname)
    X = Vector{T}(undef, length(Q))
    for (k,v) in Q
        ii = parse(Int64, replace(k,"idx" => ""))
        X[ii] = convert(T,v)
    end
    X
end

function Base.convert(::Type{T}, Q::Dict{String, Any}) where T <: Union{DPHData, DPHDataArgs}
    a_args = Any[]
    for f in fieldnames(T)
        tt = fieldtype(T,f)
        fs = string(f)
        fs = sanitise(fs)
        if tt <: Symbol
            vv = Symbol(Q[fs])
        elseif tt <: DPHDataArgs
            #handle arguments to other types here
            vv = convert(tt, Q[fs])
        elseif tt <: AbstractVector && !(typeof(Q[fs]) <: AbstractVector)
            vv = eltype(tt)[Q[fs];]
        elseif tt <: AbstractMatrix && !(typeof(Q[fs]) <: AbstractMatrix)
            vv = fill(Q[fs], 1,1)
        else
            vv = Q[fs]
        end
        push!(a_args, vv)
    end
    T(a_args...)
end

"""
Return those directories among `dirs` where `func`,using arguments `args`,  returns true for an object whose arguments are compatible with  `typeargs`.
"""
function Base.filter(func::Function, typeargs::T2, dirs::Vector{String}, args...;verbose=0)  where T2 <: DPHDataArgs
    outdirs = String[]
    for d in dirs
        cd(d) do
            fname = filename(typeargs)
            if isfile(fname)
                aa = false
                try
                    X = load(typeargs)
                    aa = func(X,args...)
                catch
                    if verbose > 0
                        rethrow()
                    end
                end
                if aa
                    push!(outdirs, d)
                end
            end
        end
    end
    outdirs
end

function Base.show(io::IO, X::T) where T <: DPHDataArgs
    compact = get(io, :compact, false)
    print("$T with fields:\n")
    for f in fieldnames(T)
        v = getfield(X, f)
        if fieldtype(T, f) <: String
            print(io, "\t$f = \"$v\"\n")
        elseif fieldtype(T, f) <: Symbol
            print(io, "\t$f = :$v\n")
        else
            print(io, "\t$f = $v\n")
        end
    end
end

function findargs(::Type{T}, cwd=pwd();kvs...) where T <: DPHData
	fname = filename(T)
	fname = replace(fname, ".mat" => "*.mat")
	arg_type = fieldtype(T, :args)
	targs = filter(k->k[1] in fieldnames(arg_type), kvs)
	args = cd(cwd) do
		files = glob(fname)
		args = arg_type[]
		for (ii,ff) in enumerate(files)
			X = load(T, ff)
			found = true
			for (k,v) in targs
				if getfield(X.args,k) != v
					found = false
					break
				end
			end
			if found
				push!(args, X.args)
			end
		end
		args
	end
	args
end

function reset!(args::DPHDataArgs)
    reset!(filename(args))
end

"""
Unlocks the file pointed to be `fname` if it is under git annex control, so that it can be overwritten.
"""
function reset!(fname::String)
    if islink(fname)
        if git_annex != nothing
            run(`$(git_annex()) get $fname`)
            run(`$(git_annex()) unlock $fname`)
        end
    end
end
end # module


text:
using Documenter, NormalSplines

makedocs(
    sitename = "NormalSplines.jl",
	format = Documenter.HTML(),
	authors = "Igor Kohanovsky",
    pages = [
				"Home" => "index.md",
				"Example Usage" => "Usage.md",
				"Public API" => "Public-API.md",
			]
)

deploydocs(
	     repo = "github.com/IgorKohan/NormalSplines.jl.git",
	    )


text:
"""
    MeshUtilModule

Module for mesh utility functions used in other meshing modules.
"""
module MeshUtilModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import Base.BitSet
import LinearAlgebra: norm

"""
    HyperFaceContainer

Type of a hyper face data to store a newly generated node.
"""
mutable struct HyperFaceContainer
    o::BitSet # numbers of the other nodes on the hyperface
    n::Int # new node number generated on the hyperface
end

"""
    makecontainer()

Make hyper face container.

This is a dictionary of hyper faces, indexed with an anchor node. The anchor
node is the smallest node number within the connectivity of the hyper face.
"""
makecontainer() = Dict{FInt, Array{HyperFaceContainer}}();

"""
    addhyperface!(container,hyperface,newn)

Add a hyper face to the container.

The new node is stored in hyper face data in the container and can be retrieved
later.
"""
function addhyperface!(container,hyperface,newn)
    h=sort([i for i in hyperface])
    anchor=h[1]; other= BitSet(h[2:end]);
    C=get(container,anchor,HyperFaceContainer[]);
    fnd=false;
    for k=1:length(C)
        if C[k].o == other
            fnd=true; break;
        end
    end
    if fnd!=true
        push!(C,HyperFaceContainer(other,newn)); newn=newn+1;#new node added
    end
    container[anchor]= C; # set the new array
    return newn # return the number of the added node, possibly unchanged
end

"""
    findhyperface!(container,hyperface)

Find a hyper face in the container.

If the container holds the indicated hyper face, it is returned together with
the stored new node number.
"""
function findhyperface!(container,hyperface)
    h=sort([i for i in hyperface])
    anchor=h[1]; other= BitSet(h[2:end]);
    C=get(container,anchor,HyperFaceContainer[]);
    for k=1:length(C)
        if C[k].o == other
            return [anchor, C[k].o], C[k].n
        end
    end
    return [], 0
end

function ontosphere(xyz::FFltMat,radius::FFlt)
# % Project nodes onto a sphere of given radius.
# %
# % function fens= onto_sphere(fens,radius,list)
# %
# % fens = finite element node set,
# % radius = radius of the sphere,
# % list = optional argument, if not empty then only
# %	     the nodes in the list are to be moved; otherwise all nodes are moved.

    for j=1:size(xyz,1)
        xyz[j,:] =xyz[j,:]*radius/norm(xyz[j,:]);
    end
    return xyz;
end

"""
    linearspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate linear space.

Generate a linear sequence of `N` numbers between  `start` and `stop` (i. e.
sequence  of number with uniform intervals inbetween).

# Example
```
julia> linearspace(2.0, 3.0, 5)
2.0:0.25:3.0
```
"""
function linearspace(start::T, stop::T, N::Int)  where {T<:Number}
    return range(start, stop = stop, length = N)
end

"""
    gradedspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate quadratic space.

Generate a quadratic sequence of `N` numbers between `start` and `stop`. This
sequence corresponds to separation of adjacent numbers that increases linearly
from start to finish.

# Example
```
julia> gradedspace(2.0, 3.0, 5)
5-element Array{Float64,1}:
 2.0
 2.0625
 2.25
 2.5625
 3.0
```
"""
function gradedspace(start::T, stop::T, N::Int, strength=2)  where {T<:Number}
    x = range(0.0, stop = 1.0, length = N);
    x = x.^strength
    # for i = 1:strength
    #     x = cumsum(x);
    # end
    x = x/maximum(x);
    out = start .* (1.0 .- x) .+ stop .* x;
end

# function inbox(box::FFltVec,sdim::FInt,x::FFltVec)
#     for i=1:sdim
#         if (!inrange(box[2*i-1],box[2*i],x[i]))
#             return false
#         end
#     end
#     return true
# end

"""
    logspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate logarithmic space.

Generate a logarithmic sequence of `N` numbers between  `start` and `stop` (i. e.
sequence  of number with uniform intervals inbetween in the log space).

# Example
```
julia> logspace(2.0, 3.0, 5)
5-element Array{Float64,1}:
  100.0
  177.82794100389228
  316.2277660168379
  562.341325190349
 1000.0
```
"""
function logspace(start::T, stop::T, N::Int)  where {T<:Number}
    es = range(start, stop = stop, length = N)
    return [10^e for e in es]
end

end


text:
module GridapP4est
  using FillArrays
  using Gridap
  using PartitionedArrays
  using GridapDistributed
  using P4est_wrapper
  const PArrays = PartitionedArrays
  include("UniformlyRefinedForestOfOctreesDiscreteModels.jl")
  export UniformlyRefinedForestOfOctreesDiscreteModel
end


text:
module PointCloud
using jInv.Mesh
using jInv.Utils
using jInv.InverseSolve
using ShapeReconstructionPaLS.Utils
using ShapeReconstructionPaLS.ParamLevelSet
using MAT
using SparseArrays
using Distributed
import jInv.ForwardShare.getData
import jInv.ForwardShare.getSensTMatVec
import jInv.ForwardShare.getSensMatVec
import jInv.ForwardShare.ForwardProbType

## The rotation planes are:
# first rotation (theta) is for planes x1 and x2
# second rotation (phi) is for planes x1 and x3
# Sampling is done by summing planes of x1 and x2 for all values of x3.

export PointCloudParam, getPointCloudParam
mutable struct PointCloudParam <: ForwardProbType
	Mesh      				:: RegularMesh
	P						:: Array{Array{Float64,2}}
	Normals					:: Array{Array{Float64,2}}
	npcAll					:: Int64
	workerSubIdxs
    theta_phi_rad			:: Array{Float64,2}
	b						:: Array{Float64,2}
	method					:: String
	Jacobian				:: SparseMatrixCSC{Float32,Int32}
end

function getPointCloudParamInternal(Mesh::RegularMesh, P:: Array{Array{Float64,2}}, Normals::Array{Array{Float64,2}},npcAll::Int64,workerSubIdxs,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},method)
	return PointCloudParam(Mesh::RegularMesh,P::Array{Array{Float64,2}}, Normals::Array{Array{Float64,2}}, npcAll, workerSubIdxs, theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2}, method,spzeros(Float32,Int32,0,0));
end

function getPointCloudParam(Mesh::RegularMesh,P::Array{Array{Float64,2}}, Normals::Array{Array{Float64,2}},theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},numWorkers::Int64,method = MATBased)
	## This function does use the parallel mechanism of jInv (i.e., returns a RemoteChannel), even if numWorkers=1.
	if numWorkers > nworkers()
		numWorkers = nworkers();
	end
	SourcesSubInd = Array{Array{Int64,1}}(undef,numWorkers);
	ActualWorkers = workers();
	if numWorkers < nworkers()
		ActualWorkers = ActualWorkers[1:numWorkers];
	end
	pFor   = Array{RemoteChannel}(undef,numWorkers);
	i = 1; nextidx() = (idx=i; i+=1; idx)
	idx=i;
	npc  = 2
	# send out jobs
	@sync begin
		for w = ActualWorkers
			@async begin
				while true
					idx = nextidx();
					if idx > numWorkers
						break
					end
					I_p = getIndicesOfKthWorker(numWorkers,idx,npc);
					tmp = initRemoteChannel(getPointCloudParamInternal,w,Mesh::RegularMesh,P::Array{Array{Float64,2}}, Normals::Array{Array{Float64,2}},npc ,I_p,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},method);
					pFor[1]  = initRemoteChannel(getPointCloudParamInternal,w,Mesh::RegularMesh,P::Array{Array{Float64,2}}, Normals::Array{Array{Float64,2}},npc,I_p,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},method);
					wait(pFor[1]);

				end
			end
		end
	end
	return pFor # Array of Remote Refs
end

import jInv.Utils.clear!
function clear!(pFor::PointCloudParam)
	pFor.SampleMat = spzeros(0);
	clear!(pFor.Mesh);
	return pFor;
end

include("prepareSyntheticPointCloudData.jl");
include("getData.jl");
include("sensitivityFuncs.jl");
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
#################### Mathematical Operators ####################
function cummean(x::AbstractArray)
    return mapslices(cummean, x, dims = 1)
end

function cummean(x::AbstractVector)
    y = similar(x, Float64)
    xs = 0.0
    fill!(y, xs)

    c = 0
    for i in 1:length(x)
        if !ismissing(x[i])
            c += 1
            xs += x[i]
        end
        y[i] = xs / c
    end
    return y
end

## Csorgo S and Faraway JJ. The exact and asymptotic distributions of the
## Cramer-von Mises statistic. Journal of the Royal Statistical Society,
## Series B, 58: 221-234, 1996.
function pcramer(q::Real)
    p = 0.0
    for k in 0:3
        c1 = 4.0 * k + 1.0
        c2 = c1^2 / (16.0 * q)
        p += gamma(k + 0.5) / factorial(k) * sqrt(c1) * exp(-c2) * besselk(0.25, c2)
    end
    return p / (pi^1.5 * sqrt(q))
end


text:
module ReadICARTT

import DataStructures: OrderedDict;
using Dates;
using Unitful;

using ..ICARTTUnits;
using ..ICARTTExceptions: ICARTTParsingException, ICARTTNotImplementedException;

export read_icartt_file;

##############
# Exceptions #
##############

_no_default(name) = error("$name is a required argument");

"""
    _HeaderField(fields, readfxn)

A struct used internally to decribed how to read a field in the ICARTT header,
before the "bulk" part. `fields` must be a tuple of strings that are the key(s)
in the metadata dict to store the value(s). `readfxn` must be a function that,
given the line read in from the ICARTT file, parses it and returns the value(s)
to store. The number of values stored must equal the number of fields, i.e. in:

```
julia> _HeaderField(("Acquisition_date", "Processing_date"), _icartt_read_data_dates)
```

`_icartt_read_data_dates` must return two values; the first is stored in the metadata
dict under "Acquisition_date" and the second under "Processing_date". This is
how you handle when multiple pieces of information are defined on one line.
"""
struct _HeaderField
    fields::Tuple
    readfxn::Function
end

"""
    MergeDataField{T <: Unitful.Quantity}(name, unit, fill, scale, values)

This struct is used to represent the data in the ICARTT file. It has fields "name",
"unit", "fill", "scale", and "values".

    * "name" is the name of the variable, as a String
    * "unit" is the units of the variable, as a Unitful.Units instance
    * "fill" is the fill value of the variable, as a Float
    * "scale" is the scale factor of the variable, as a Float
    * "values" is a 1D array with elements of type T. It is required that T be
      a child of Unitful.Quantity; this enforces keeping the units of the measurements
      with them.

This struct has both positional and named constructors for clarity., e.g. the two
methods will create the same struct:
```
julia> MergeDataField("Altitude", u"m", -999999.0, 1.0, Array{1.0u"m", 1}(...))
julia> MergeDataField(name="Altitude", unit=u"m", fill=-999999.0, scale=1.0, values=Array{1.0u"m", 1}(...))
```
Note however that the two cannot be mixed and matched.
"""
struct MergeDataField{T <: Unitful.Quantity}
    name::String;
    unit::Unitful.Units;
    fill::Float64;
    scale::Float64;
    values::Array{T, 1};  # in theory, I think the type could be inferred from "unit"
    MergeDataField(name, unit, fill, scale, values::Array{T,1}) where T <: Unitful.Quantity = new{T}(name, unit, fill, scale, values);
    MergeDataField(;name=_no_default("name"), unit=_no_default("unit"), fill=_no_default("fill"),
                    scale=_no_default("scale"), values::Array{T,1}=_no_default("values")) where T <: Unitful.Quantity = new{T}(name, unit, fill, scale, values);
end

"""
    AirMerge(data, metadata)

A struct that represents a whole ICARTT file. `data` is the ordered dict of data
fields returned by `_read_icartt_data` and `metadata` the ordered dict of metadata
returned by `_read_icartt_metadata`.
"""
struct AirMerge
    data::OrderedDict{String, MergeDataField};
    metadata::OrderedDict{String, Any};
end

"""
    _icartt_read_simple(line)

Function to read in header line that needs no conversion except to be returned
in a tuple so that indexing the returned result gives the whole line.
"""
function _icartt_read_simple(line)
    return (line,)
end

"""
    _icartt_read_data_dates(line)

Specialized function to read the two data dates from the ICARTT header. Since both
the acquisition and processing dates are given on the same line, this separates them
and returns two separate `Date` structs.
"""
function _icartt_read_data_dates(line)
    acq_yr, acq_mn, acq_day, proc_yr, proc_mn, proc_dy = _split_icartt_line(line)
    acq_date = Date("$acq_yr-$acq_mn-$acq_day", "y-m-d")
    proc_date = Date("$proc_yr-$proc_mn-$proc_dy", "y-m-d")
    return acq_date, proc_date
end

#############
# Constants #
#############

# Define the line numbers on which each of these pieces of information should
# be. This is only for the one-line pieces of information; multiline pieces are
# handled by their own custom functions.
const _ICARTT_HDR_PI = 2
const _ICARTT_HDR_PI_AFFIL = _ICARTT_HDR_PI + 1
const _ICARTT_HDR_DESCRIPTION = _ICARTT_HDR_PI_AFFIL + 1
const _ICARTT_HDR_MISSION = _ICARTT_HDR_DESCRIPTION + 1
const _ICARTT_HDR_FILE_VOL_NUM = _ICARTT_HDR_MISSION + 1
const _ICARTT_HDR_DATA_UTC_DATES = _ICARTT_HDR_FILE_VOL_NUM + 1
const _ICARTT_HDR_DATA_INTERVAL_SEC = _ICARTT_HDR_DATA_UTC_DATES + 1
const _ICARTT_HDR_INDEPENDENT_VAR = _ICARTT_HDR_DATA_INTERVAL_SEC + 1
const _ICARTT_HDR_N_VAR = _ICARTT_HDR_INDEPENDENT_VAR + 1
const _ICARTT_HDR_SCALE_FACTORS = _ICARTT_HDR_N_VAR + 1
const _ICARTT_HDR_FILL_VALS = _ICARTT_HDR_SCALE_FACTORS + 1
const _ICARTT_HDR_START_HEADER_BULK = _ICARTT_HDR_FILL_VALS + 1

# This dictionary maps the line number to the piece of information it contains.
# See the HeaderField struct docstring for more info.
_ICARTT_HDR_SPECIALS = Dict(
    _ICARTT_HDR_PI => _HeaderField(("PI",), _icartt_read_simple),
    _ICARTT_HDR_PI_AFFIL => _HeaderField(("PI_affiliation",), _icartt_read_simple),
    _ICARTT_HDR_DESCRIPTION => _HeaderField(("Data_description",), _icartt_read_simple),
    _ICARTT_HDR_MISSION => _HeaderField(("Mission",), _icartt_read_simple),
    _ICARTT_HDR_DATA_UTC_DATES => _HeaderField(("Acquisition_date", "Processing_date"), _icartt_read_data_dates),
    _ICARTT_HDR_DATA_INTERVAL_SEC => _HeaderField(("Data_interval",), _icartt_read_simple),
    _ICARTT_HDR_INDEPENDENT_VAR => _HeaderField(("Independent_variable",), _icartt_read_simple),
)

# This tuple lists the categories of information expected in the normal comments
# section. These will be used in a regular expression, so regex wildcards can be
# used.
_ICARTT_NORMAL_COMMENTS = ("PI_CONTACT_INFO", "PLATFORM", "LOCATION", "ASSOCIATED_DATA",
                           "INSTRUMENT_INFO", "DATA_INFO", "UNCERTAINTY",
                           "ULOD_FLAG", "ULOD_VALUE", "LLOD_FLAG", "LLOD_VALUE",
                           "DM_CONTACT_INFO", "PROJECT_INFO", "STIPULATIONS_ON_USE",
                           "OTHER_COMMENTS", "REVISION", "R\\d+");

"""
    read_icartt_file(filename::String; extra_aliases=nothing, extra_alias_files=nothing, alias_dict_overwrite=false, no_default_aliases=false, verbose=0)

Reads data files formatted to the ICARTT standard defined at
https://www-air.larc.nasa.gov/missions/etc/IcarttDataFormat.htm
and
https://cdn.earthdata.nasa.gov/conduit/upload/6158/ESDS-RFC-029v2.pdf
Returns an `AirMerge` structure containing the ICARTT metadata and data.

`verbose` controls the level of printing to the console, for debugging purposes.
Set to higher numbers to give more information.

The data from the ICARTT file is stored as arrays of Unitful `Quantity`s that
combine the value and unit from the file. The helps ensure that calculations
don't miss conversion factors, but requires that the units defined in the ICARTT
header are written in a way that Unitful can understand. To aid that, the units
are reformatted by `ICARTTUnits.sanitize_raw_unit_strings()`. One of the steps
it takes is to replace undefined units with their Unitful equivalent. What units
it replaces are defined by the alias dictionary, which can be modified by the
following keyword arguments.

`no_default_aliases` will eliminate the default aliases defined in `ICARTTUnits`.

`extra_aliases`, `extra_alias_files`, and `alias_dict_overwrite` allow you to
specify additional unit aliases.

An example is when time is given in units of "hours"; Unitful can't understand
"hours" so we need to replace it with the abbreviation that it does understand,
"hr". This function has three ways to specify what unit strings should be
replaced:

1. ICARTT.ICARTTUnits has some built-in aliases. These are stored in the
   dictionary `ICARTT.ICARTTUnits.unit_aliases`, or can be printed by
   `ICARTT.ICARTTUnits.list_default_unit_aliases()`. You can modify this
   dictionary to add new units globally, though this is not the recommended
   method. The format of the dictionary is the same
2. Pass as `extra_aliases` a dictionary that defines the desired units as the keys
   as the strings that should be replaced in arrays as the values. That is the
   following dict would define replacing "hour" with "hr" and any of "deg", "degs",
   or "Degs" with "°":

```
   extra_aliases = Dict("hr" => ["hour"], "°" => ["deg", "degs", "Degs"]);
```

   Note that the behavior when one of the units defined in the `extra_aliases` is
   already defined depends on the value of `alias_dict_overwrite`. If it is `false`
   (default), the definitions in `extra_aliases` are added to any existing definitions.
   If `true`, then existing definitions are replaced. For example, the default
   aliases include one to replace "nanometers" with "nm". If

```
   extra_aliases = Dict("nm" => "nanometres")
```

   i.e. it uses the British spelling, then with `alias_dict_overwrite = false`,
   both "nanometers" and "nanometres" will be replaced by "nm". If `alias_dict_overwrite = true`,
   then *only* "nanometres" will be replaced.

3. Pass names of one or more text files that describe the same mapping as the
   dictionary as `extra_alias_files`. May pass a single file as a string, or more
   than one as an array of strings. The format of these files is:

```
   replace with: to replace[, to replace[, to replace]][: {replace,append}]
```

   Each line is the same as a key-value pair in the `extra_aliases` dictionary.
   The standard aliases are defined in such a file, the first few lines are:

   °: degs, deg, Degs
   hr: hour
   percent: %
   std_m: std m
   nm: nanometers
   : #

   This says to replace "degs", "deg", or "Degs" with "°", "hour" with "hr", and so
   on. Note that "std m" (including the space) gets replaced with "std_m", and "#"
   with nothing. For both the strings to replace and to replace them with, spaces
   inside of letters matter, so e.g. in "std m" even though the full string between
   the colon and comma is " std m", the leading (and if applicable, trailing) spaces
   are stripped.

   For your text files, the final "append" or "replace" after an optional second
   colon controls what happens if your file declares aliases for a unit that already
   had aliases. If your file had a line:

```
   °: degree, degrees
```

   then any of "degs", "deg", "Degs", "degree", or "degrees" would be replaced
   with "°" (the first three from the defaults, the last two from your file). If
   instead the line was:

```
   °: degree, degrees : replace
```

   then *only* "degree" and "degrees" would be replaced with "°". This allows you
   to override defaults selectively. Note that the order you pass the files in
   matters. If you had two files:

```
   # File A.txt:
   °: degree, degrees : replace

   # File B.txt:
   °: Deg
```

   then passing `extra_alias_files=["File A.txt", "File B.txt"]` would mean that
   "degree", "degrees", and "Deg" would all be replaced by "°", while
   `extra_alias_files=["File B.txt", "File A.txt"]` means *only* "degree" and
   "degrees" get replaced by "°" - the replace line in File A wipes out the
   corresponding line in File B as well as the defaults.
"""
function read_icartt_file(filename::String; extra_aliases=nothing, extra_alias_files=nothing,
        alias_dict_overwrite=false, no_default_aliases=false, verbose=0)
    metadata, data = open(filename, "r") do io
        metadata, variable_info = _read_icartt_metadata(io; extra_aliases=extra_aliases,
            extra_alias_files=extra_alias_files, alias_dict_overwrite=alias_dict_overwrite,
            no_default_aliases=no_default_aliases, verbose=verbose);
        data = _read_icartt_data(io, variable_info; verbose=verbose);
        return metadata, data
    end
    return AirMerge(data, metadata);
end

"""
    _read_icartt_metadata(io::IO; extra_aliases=nothing, extra_alias_files=nothing, alias_dict_overwrite=false, no_default_aliases=false, verbose=0)

Given a handle to an open ICARTT file (`io`), reads in the metadata in the header
of the file and returns it as a dictionary. Note that this assumes that the file
pointed to by `io` has been opened, but nothing read from it (or that the position
of the IO cursor has been returned to the start of the file).

`verbose` controls the level of printing to the console, for debugging purposes.
Set to higher numbers to give more information.

See `read_icartt_file` for a description of `extra_aliases`, `extra_alias_files`,
and `alias_dict_overwrite`.
"""
function _read_icartt_metadata(io::IO; extra_aliases=nothing, extra_alias_files=nothing, alias_dict_overwrite=false, no_default_aliases=false, verbose=0)
    if verbose > 0
        println("Reading metadata")
    end

    metadict = OrderedDict{String, Any}();
    aliases_dict = _setup_unit_aliases(extra_aliases, extra_alias_files; dict_overwrite=alias_dict_overwrite, no_default_aliases=no_default_aliases, verbose=verbose);

    # read the first line to get the number of lines in the header and the file format index
    n_header_lines, ffi = _split_icartt_line(readline(io));
    n_header_lines = parse(Int64, n_header_lines);
    metadict["file_format_index"] = ffi;

    # now read in the rest of the header lines
    n_line = 2;
    bulk_mode = "variable_names";
    # init variables to be set in the if statement due to scoping rules
    n_var = nothing;
    scale_factors = nothing;
    fill_values = nothing;
    variables = nothing;
    table_header = nothing;

    while n_line <= n_header_lines
        if n_line < _ICARTT_HDR_START_HEADER_BULK
            # While we're in the initial part of the header, where there is one
            # piece of information per line, we read in each line first, then
            # decide what to do with it.
            line = readline(io);
            if verbose > 1
                println("Read line: $line")
            end

            # since ICARTT isn't a fully self-describing format, we have to hard-code
            # some of the meanings of the header elements.
            if n_line in keys(_ICARTT_HDR_SPECIALS)
                # Deal with the one-line header information fields that we save into
                # the metadata.
                fields = _ICARTT_HDR_SPECIALS[n_line].fields;
                readfxn = _ICARTT_HDR_SPECIALS[n_line].readfxn;
                values = readfxn(line);
                for i_val = 1:length(fields)
                    metadict[fields[i_val]] = values[i_val];
                end
            elseif n_line == _ICARTT_HDR_N_VAR
                n_var = parse(Int64, line);
            elseif n_line == _ICARTT_HDR_SCALE_FACTORS
                scale_factors = _split_icartt_array(Int64, line);
            elseif n_line == _ICARTT_HDR_FILL_VALS
                fill_values = _split_icartt_array(Float64, line);
            else
                if verbose > 0
                    println("ICARTT parsing: ignoring line $n_line ($line)");
                end
            end

            n_line += 1;
        else
            # Once we get into what I call the "bulk" part of the header, where
            # information for one particular category is spread across multiple
            # lines, it makes more sense to let the subfunctions handle reading
            # the lines, since it's cleaner to let them determine how many
            # lines to read.
            if bulk_mode == "variable_names";
                if verbose > 0
                    println("Reading variables");
                end
                variables, n_line = _parse_variable_names(io, n_var, n_line; aliases_dict=aliases_dict, verbose=verbose);
                bulk_mode = "special_comments";

            elseif bulk_mode == "special_comments"
                if verbose > 0
                    println("Reading special comments");
                end
                metadict["Special_comments"], n_line = _parse_special_comments(io, n_line);
                bulk_mode = "normal_comments";

            elseif bulk_mode == "normal_comments"
                if verbose > 0
                    println("Reading normal comments");
                end
                comments, table_header, n_line = _parse_normal_comments(io, n_line);
                metadict = merge(metadict, comments);
                bulk_mode == "complete";

            elseif bulk_mode == "complete"
                throw(ICARTTParsingException("Still trying to parse the header,
                but have already finished parsing the last section (normal comments)"));

            else
                throw(ICARTTNotImplementedException("No behavior defined for
                bulk_mode = '$bulk_mode'"));
            end
        end
    end

    # Almost done. The last three things we want to do are 1) merge the lists of
    # scale factors and fill values into the array of descriptive tuples,
    # 2) insert the primary variable into the front of the list of variables and
    # their units, then 3) check that the list of variables matches the table header.
    variables = _merge_variable_info(variables, fill_values, scale_factors)

    temp_tup = _parse_variable(metadict["Independent_variable"]);
    primary_unit_tup = (field=temp_tup.field, unit=temp_tup.unit, fill=NaN, scale=1.0)
    splice!(variables, 1:0, [primary_unit_tup]);

    _check_variables_vs_header(variables, table_header);

    return metadict, variables
end

"""
    _read_icartt_data(io::IO, variable_info; verbose=0)

Reads the data table section of an ICARTT file, given a handle to an open ICARTT
file (`io`) and an array of named tuples (`variable_info`) describing the
name ("field", String), unit ("unit", Unitful.Units), fill value ("fill", Float),
and scale factor ("scale", Float) of each variable, in order.

`verbose` controls the level of printing to the console, for debugging purposes.
Set to higher numbers to give more information.

Returns a dictionary with String keys (variable names) and MergeDataField values.
"""
function _read_icartt_data(io::IO, variable_info; verbose=0)
    # go ahead and read in the rest of the data. we can use this to figure
    # out the array length to preallocate
    if verbose > 0
        println("Reading data.")
    end

    table_lines = [_split_icartt_array(Float64, l) for l in readlines(io)];
    n_meas = length(table_lines);
    # setup a dict of the final structs to receive the data
    data = OrderedDict{String, MergeDataField}();

    for (i, info) in enumerate(variable_info)
        data_array = Array{typeof(Unitful.Quantity(1.0,info.unit)), 1}(undef, n_meas);
        if verbose > 1
            println("Reading variable #$i: $(info.field) $(info.unit)")
        end
        for j = 1:n_meas
            data_array[j] = Unitful.Quantity(table_lines[j][i], info.unit);
        end

        data[info.field] = MergeDataField(name=info.field,
                                          unit=info.unit,
                                          fill=info.fill,
                                          scale=info.scale,
                                          values=data_array);
    end

    return data;
end

#####################
# Utility functions #
#####################

"""
    _split_icartt_line(line)

Split a line in an ICARTT file into an array of strings. Assumes that separate
values are separated by commas, and automatically strips whitespace from each
component.
"""
function _split_icartt_line(line)
    return [strip(s) for s in split(line, ",")];
end

"""
    _split_icartt_array(T, line)

Takes a line of an ICARTT file (`line`) and converts it to an array of type `T`
using _split_icartt_line to separate distinct values.
"""
function _split_icartt_array(T, line)
    return [parse(T, s) for s in _split_icartt_line(line)];
end

"""
    _setup_unit_aliases(extra_aliases::Nothing, extra_alias_files; dict_overwrite=false, no_default_aliases=false, verbose=0)
    _setup_unit_aliases(extra_aliases::AbstractDict, extra_alias_files::Nothing; dict_overwrite=false, no_default_aliases=false, verbose=0)
    _setup_unit_aliases(extra_aliases::AbstractDict, extra_alias_files::AbstractString; dict_overwrite=false, no_default_aliases=false, verbose=0)
    _setup_unit_aliases(extra_aliases::AbstractDict, extra_alias_files::AbstractArray; dict_overwrite=false, no_default_aliases=false, verbose=0)

Combine the standard unit aliases dictionary with additional configuration files
and/or an additional dictionary. The different methods handle different cases;
if no extra dictionary is desired, pass `nothing` for `extra_aliases`. If no
additional alias files desired, pass `nothing` for `extra_alias_files` as well.
Otherwise, `extra_alias_files` can be a string or array of strings.

`dict_overwrite` will cause entries in the `extra_aliases` to override those
defined in the defaults or the extra files. `no_default_aliases` will omit the
predefined aliases from the mapping.
"""
function _setup_unit_aliases(extra_aliases::Nothing, extra_alias_files; dict_overwrite=false, no_default_aliases=false, verbose=0)
    return _setup_unit_aliases(Dict{String, Array{AbstractString,1}}(), extra_alias_files; dict_overwrite=dict_overwrite, verbose=verbose);
end

# A note on types for future: there is a distinction between specifying a type
# as in extra_alias_files::AbstractString and the "internal" types of a dict or
# array. In the former, it seems that Julia's dispatch will match if extra_alias_files
# is any subtype of AbstractString and there is not a more specific method
# available. However, in the latter, for e.g. extra_aliases::AbstractDict, that will
# also match for any subtype of AbstractDict, but if we want to specify the internal
# types of the dict, then e.g. AbstractDict{AbstractString, AbstractString} will *not*
# match any subtype of AbstractDict with any subtype of AbstractString for the keys
# and values; in this case, we must *explicitly* indicate that subtypes are acceptable
# with the AbstractDict{<:AbstractString,<:AbstractString} notation.
#
# This particular case is even messier b/c we want a dict with any type of array
# containing any type of strings as the values. Consider the following:
#
# julia> isa(Dict{String,Array{String,1}}(), Dict{String,<:AbstractArray{String,1}})
# true
#
# julia> isa(Dict{String,Array{String,1}}(), Dict{String,<:AbstractArray{AbstractString,1}})
# false
#
# In the first, the concrete dict has arrays as values, and since arrays are subtypes of
# AbstractArray *and* have the *exact same type* for their values, the concrete dict
# is a subtype of the given datatype. In the second though, since the expected type
# expects some type of array with AbstractStrings, it doesn't match. We must do:
#
# julia> isa(Dict{String,Array{String,1}}(), Dict{String,<:AbstractArray{<:AbstractString,1}})
# true

function _setup_unit_aliases(extra_aliases::AbstractDict{<:AbstractString, <:AbstractArray{<:AbstractString,1}}, extra_alias_files::Nothing;
        dict_overwrite=false, no_default_aliases=false, verbose=0)
    return _setup_unit_aliases(extra_aliases, Array{String,1}(); verbose=verbose);
end

function _setup_unit_aliases(extra_aliases::AbstractDict{<:AbstractString, <:AbstractArray{<:AbstractString,1}}, extra_alias_files::AbstractString;
        dict_overwrite=false, no_default_aliases=false, verbose=0)
    return _setup_unit_aliases(extra_aliases, [extra_alias_files]; verbose=verbose);
end

function _setup_unit_aliases(extra_aliases::AbstractDict{<:AbstractString, <:AbstractArray{<:AbstractString,1}},
                             extra_alias_files::AbstractArray{<:AbstractString,1};
                             dict_overwrite=false, no_default_aliases=false, verbose=0)
    if no_default_aliases
        aliases = Dict{String, Array{AbstractString,1}}();
    else
        aliases = copy(ICARTTUnits.unit_aliases);
    end
    for f in extra_alias_files
        if verbose > 1
            println("Reading alias file $f...")
        end
        ICARTTUnits._read_alias_config!(aliases, f; verbose=verbose);
    end

    if dict_overwrite
        merge!(aliases, extra_aliases);
    else
        for (key, val) in extra_aliases
            if key in keys(aliases)
                append!(aliases[key], val);
            else
                aliases[key] = val;
            end
        end
    end
    return aliases;
end

"""
    _parse_variable_names(io, n_vars, n_lines; aliases_dict=nothing, verbose=0)

Parse the list of variables in the ICARTT header file pointed to by `io`. Assumes
that the next call to `readline(io)` will return the first line defining a variable.
`n_vars` must be an integer specifying how many variables there are to read.
`n_line` is the current line number of the ICARTT file; it will be advanced properly
and returned with the variables.

`aliases_dict` is the dictionary that defines what parts of the unit string to
replace so that Unitful understands the definition.

`verbose` controls the level of printing to the console, for debugging purposes.
Set to higher numbers to give more information.

Returns:
    1) an array of named tuples with fields "field" and "unit" describing the
name and unit of the variable. "unit" will be a Unitful.Units instance.
    2) the line number at the end of the variables section of the header.
"""
function _parse_variable_names(io::IO, n_vars, n_line; aliases_dict=nothing, verbose=0)
    # must retain order to be able to match them up with their fill values
    varnames = Array{NamedTuple{(:field, :unit), Tuple{String,Unitful.Units}}, 1}(undef, n_vars)
    for i = 1:n_vars
        n_line += 1
        # assuming the variable name lines are of the form "name, unit"
        line = readline(io);
        if verbose > 1
            println("$i: $line");
        end
        varnames[i] = _parse_variable(line, n_line; aliases_dict=aliases_dict);
    end
    return varnames, n_line;
end

"""
    _parse_variable(line, n_line=-1; aliases_dict=nothing)

Given a line read in from an ICARTT file containing a variable name and unit
separated by a comma, creates a named tuple containing the variable name as a
string and the unit as a Unitful.Units instance. The fields are named "field" and
"unit", respectively.

`n_line` is optional, if given, it should specify the line number of the ICARTT
file that the variable is being read from; it will be included in an error message
if there is a problem interpreting the variable.

`aliases_dict` is the dictionary that defines what parts of the unit string to
replace so that Unitful understands the definition.
"""
function _parse_variable(line, n_line=-1; aliases_dict=nothing)
    name, unit = _split_icartt_line(line);
    punit = nothing;
    try
        punit = ICARTTUnits.parse_unit_string(unit; aliases_dict=aliases_dict)
    catch err
        # Unitful doesn't use a specific error type for unit conversion
        # errors so the best course of action is to print a message explaining
        # what unit failed and rethrow the error
        printstyled("Problem while converting unit \"$unit\" on line $n_line:\n"; color=:red)
        rethrow(err)
    end
    return (field=name, unit=punit);
end

"""
    _parse_special_comments(io, n_line)

Parse the special comments section of the ICARTT file pointed to by IO handle
`io`. `n_line` is the current line number, it will be advanced as the special
comments are read in.

Returns the special comments as a single string and the line number at the end
of the special comments section.
"""
function _parse_special_comments(io, n_line)
    # The number of special comment lines is the next line to read.
    n_cmt_lines = parse(Int64, readline(io));
    n_line += 1;
    special_comments = "";
    for i = 1:n_cmt_lines
        n_line += 1;
        line = readline(io, keep=true);
        special_comments *= line;
    end
    return special_comments, n_line;
end

"""
    _parse_normal_comments(io, n_line)

Parse the normal comments section of the ICARTT file pointed to by IO handle
`io`. `n_line` is the current line number, it will be advanced as the normal
comments are read in.

Returns the comments as an OrderedDict containing the standard comment categories
as keys and their values as the string values, also returns the line number at
the end of the normal comments section.
"""
function _parse_normal_comments(io, n_line)
    n_cmt_lines = parse(Int64, readline(io));
    n_line += 1;
    comments = OrderedDict{String, Any}();

    # Unlike the special comments section, which is intended to be pretty open
    # for whatever comments are needed, the normal comments section has specific
    # piece of information that are expected to be included. This should always
    # start their line and be followed immedately by a colon. After the line
    # that indicates the number of normal comments to read, we start searching
    # for this pattern.
    #
    # Note that section 2.3.2.17 of https://cdn.earthdata.nasa.gov/conduit/upload/6158/ESDS-RFC-029v2.pdf
    # indicates that the normal comments section may include a free form text
    # section that is not currently implemented by this reader. This section may
    # include both unformatted text and custom keyword-value pairs.
    current_cmt = nothing;
    categories = join(_ICARTT_NORMAL_COMMENTS, "|");
    re = Regex("^($categories)(?=:)");
    for i = 1:(n_cmt_lines - 1)
        line = readline(io);
        n_line += 1;

        # See if this line starts a new category:
        m = match(re, line);
        if m != nothing
            # This line starts a new category. Set that as the category we'll
            # store the value in and initialize that value in the dictionary with
            # the rest of this line.
            current_cmt = m.match;
            comments[current_cmt] = strip(split(line, ":")[2]);
        else
            # This line does not start a new category.
            if current_cmt == nothing
                # If we don't have a current category, then there's something wrong
                # either with the format of the file or our parsing of it.
                throw(ICARTTParsingException("Trying to read a normal comment, but have not found a comment keyword yet!"))
            else
                # This must be a multiline comment, so append it to the current
                # value.
                comments[current_cmt] *= " " * line;
            end
        end
    end

    # The last line will be the table header
    line = readline(io);
    n_line += 1;
    header = _split_icartt_line(line);
    return comments, header, n_line
end

"""
    _merge_variable_info(variables, fill_vals, scale_factors)

Combine the fill values, scale factors, and units for each variable into a
single named tuple with fields "field", "unit", "fill", and "scale". Returns an
array of these tuples, one per variable.

`variable` must be the array of tuples containing the variable name and units.
`fill_vals` must be an array of the same length containing the fill values for
each unit. `scale_factors` is likewise the array of scale factors. Order matters,
the order must be the same for all three inputs.
"""
function _merge_variable_info(variables, fill_vals, scale_factors)
    new_vars = Array{NamedTuple{(:field, :unit, :fill, :scale), Tuple{String, Unitful.Units, Float64, Float64}}, 1}(undef, length(variables));
    for i in 1:length(new_vars)
        new_tup = (field=variables[i].field, unit=variables[i].unit, fill=fill_vals[i], scale=scale_factors[i]);
        new_vars[i] = new_tup;
    end
    return new_vars;
end

"""
    _check_variables_vs_header(vars, header)

Check the array of variable information against the table header. `vars` is the
array of named tuples created by `_merge_variable_info`, with the primary unit
added. `header` is the array of variable names resulting from reading the table
header which is the last line of the file header.

Throws an ICARTTParsingException if the two do not match (either are different
lengths or have a variable name that is different).
"""
function _check_variables_vs_header(vars, header)
    if length(vars) != length(header)
        throw(ICARTTParsingException("The number of variables defined ($(length(vars)))  different from the number in the table header ($(length(header)))"));
    end

    for i=1:length(vars)
        if vars[i].field != header[i]
            throw(ICARTTParsingException("Variable number $i has a different name ($(vars[i].field)) than in the table header ($(header[i]))"));
        end
    end
end

end

max_stars_repo_name: joshua-laughner/Julia-ICARTT

max_issues_repo_name: joshua-laughner/Julia-ICARTT

max_forks_repo_name: joshua-laughner/Julia-ICARTT


text:
function write_vtk(M, file_path::String)

	file_stream::IOStream = open(file_path, "w")
	@printf(file_stream, "%s", "# vtk DataFile Version 3.0\n")
	@printf(file_stream, "%s", "Generated from within Julia\n")
	@printf(file_stream, "%s", "ASCII\n")
	@printf(file_stream, "%s", "DATASET STRUCTURED_POINTS\n")
	@printf(file_stream, "%s", join(("DIMENSIONS ", string(size(M, 1)), " ", string(size(M, 2)), " ", string(size(M, 3)), "\n")))
	@printf(file_stream, "%s", "ORIGIN 0 0 0\n")
	@printf(file_stream, "%s", "SPACING 1 1 1\n")
	@printf(file_stream, "%s", join(("POINT_DATA ", string(prod(size(M))), "\n")))
	@printf(file_stream, "%s", "SCALARS values float 1\n")
	@printf(file_stream, "%s", "LOOKUP_TABLE default\n")
	for i = 1:length(M)
		@printf(file_stream, "%s\n", M[i])
	end
	close(file_stream)

	nothing
end

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
@testset "Minimal Example" begin
    @testset "One dimension" begin
        input, mom, force, dist = minimalsetup1d(10)
        @test isa(input, JuSwalbe.Inputconstants)
        @test isa(mom, JuSwalbe.Macroquant{Vector{Float64},Vector{Float64}})
        @test isa(force, JuSwalbe.Forces{Vector{Float64}})
        @test isa(dist, JuSwalbe.DistributionD1Q3)
        @test length(mom.height) == 10
        @test length(force.slip) == 10
        @test length(dist.f1) == 10
    end
    @testset "Two dimension" begin
        input, mom, force, dist = minimalsetup2d(10,5)
        @test isa(input, JuSwalbe.Inputconstants)
        @test isa(mom, JuSwalbe.Macroquant{Matrix{Float64},JuSwalbe.Twovector{Matrix{Float64}}})
        @test isa(force, JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{Float64}}})
        @test isa(dist, JuSwalbe.DistributionD2Q9)
        @test size(mom.height) == (10,5)
        @test size(force.slip.x) == (10,5)
        @test size(dist.f1) == (10,5)
    end
end

@testset "Simple constructions" begin
    @testset "Simple Moment one dimension" begin
        for type in [Float16 Float32 Float64]
            mom = simplemoment1d(5, T=type)
            @test isa(mom, JuSwalbe.Macroquant{Vector{type}, Vector{type}})
            @test length(mom.height) == 5
            for i in mom.height
                @test i == type(1.0)
            end
            @test isa(mom.velocity, Vector{type})
            @test length(mom.velocity) == 5
            for i in mom.velocity
                @test i == type(0.1)
            end
            @test length(mom.pressure) == 5
            @test length(mom.energy) == 5
        end
    end
    @testset "Simple Twovector" begin
        for type in [Float16 Float32 Float64]
            xy = simpleTwovector(5,5; T=type)
            @test isa(xy, JuSwalbe.Twovector{Matrix{type}})
            @test size(xy.x) == (5,5)
            @test size(xy.y) == (5,5)
            for i in xy.x
                @test i == type(0.1)
            end
            for i in xy.y
                @test i == type(-0.1)
            end
        end
    end
    @testset "Simple Moment two dimensions" begin
        for type in [Float16 Float32 Float64]
            mom = simplemoment2d(5,5, T=type)
            @test isa(mom, JuSwalbe.Macroquant{Matrix{type}, JuSwalbe.Twovector{Matrix{type}}})
            @test size(mom.height) == (5,5)
            for i in mom.height
                @test i == type(1.0)
            end
            @test isa(mom.velocity, JuSwalbe.Twovector{Matrix{type}})
            @test size(mom.velocity.x) == (5,5)
            for i in mom.velocity.x
                @test i == type(0.1)
            end
            @test size(mom.velocity.y) == (5,5)
            for i in mom.velocity.y
                @test i == type(-0.1)
            end
            @test size(mom.pressure) == (5,5)
            @test size(mom.energy) == (5,5)
        end
    end
    @testset "Simple D2Q9 Distribution" begin
        for type in [Float16 Float32 Float64]
            soldict = Dict(:f0 => type(1.0), :f1 => type(0.1), :f2 => type(-0.1), :f3 => type(0.01), :f4 => type(-0.01),
                           :f5 => type(0.2), :f6 => type(0.02), :f7 => type(0.02), :f8 => type(0.5))

                           xy = simpledistD2Q9(5,5; T=type)
            @test isa(xy, JuSwalbe.DistributionD2Q9{Matrix{type}})
            xydict = struct2dict(xy)
            for key in keys(xydict)
                @test size(xydict[key]) == (5,5)
                for element in xydict[key]
                    @test element == soldict[key]
                end
            end

        end
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
using CubicEoS
using CubicEoSDatabase
using Test

#= Utilities =#

"Non-recursively compare field values of `x` and `y`."
function compare_structs(x::T, y::T) where {T}
    for field in fieldnames(T)
        @test getfield(x, field) == getfield(y, field)
    end
end

@testset "CubicEoS.jl" begin
    include("dbload.jl")
    include("basic_thermo.jl")
    include("chempotential.jl")
    include("solvecubic.jl")
    include("newton.jl")
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


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

struct StepperOptions
    init_step_size::Float64
    min_step_size::Float64
    max_step_size::Float64
    min_step_growth::Float64
    max_step_growth::Float64
    step_safety_factor::Float64
end

function StepperOptions(;
    init_step_size = 0.0, # indicates recalculation
    min_step_size = 1e-12,
    max_step_size = 1.0,
    min_step_growth = 0.1,
    max_step_growth = 5.0,
    step_safety_factor = 0.9)

    return StepperOptions(
        init_step_size,
        min_step_size,
        max_step_size,
        min_step_growth,
        max_step_growth,
        step_safety_factor)
end

abstract type AbstractStepSizeControl
end

struct ClassicalStepSizeControl <: AbstractStepSizeControl
end

struct PIStepSizeControl <: AbstractStepSizeControl
    alpha_coeff::Float64
    beta_coeff::Float64
end
function PIStepSizeControl()
    alpha_coeff = 0.7
    beta_coeff  = 0.3
    return PIStepSizeControl(alpha_coeff, beta_coeff)
end

function step_size_control(
    step_size_controller::SC,
    ode_solver,
    stepper_temp,
    error_estimation::ErrorEstimation,
    rejection_break_flag,
    stepper_options::StepperOptions) where {SC <: AbstractStepSizeControl}
    error("Default implementation")
end

function step_size_control(
    step_size_controller::ClassicalStepSizeControl,
    ode_solver,
    stepper_temp,
    error_estimation::ErrorEstimation,
    rejection_break_flag,
    stepper_options::StepperOptions)

    solver_order    = order(ode_solver)
    step_size       = stepper_temp.dt
	err_est         = error_estimation.current + eps(error_estimation.current)

    s               = stepper_options.step_safety_factor
    e               = 1 / (1 + solver_order)
	step_growth     = s * (1 / err_est) ^ e

    step_size       = __bound_step_size(step_growth, step_size, stepper_options)
    stepper_temp.dt = step_size

    return step_size
end

function step_size_control(
    step_size_controller::PIStepSizeControl,
    ode_solver,
    stepper_temp,
    error_estimation::ErrorEstimation,
    rejection_break_flag,
    stepper_options::StepperOptions)

    solver_order    = order(ode_solver)
    step_size       = stepper_temp.dt
	err_est         = error_estimation.current + eps(error_estimation.current)
    err_est_prev    = error_estimation.last    + eps(error_estimation.last)
    alpha           = step_size_controller.alpha_coeff / solver_order
    beta            = step_size_controller.beta_coeff / solver_order
    s               = stepper_options.step_safety_factor
	step_growth     = s * (1 / err_est)^alpha * err_est_prev^beta
    step_size       = __bound_step_size(step_growth, step_size, stepper_options)
    stepper_temp.dt = step_size

    return step_size
end

function __bound_step_size(step_growth, step_size, stepper_options::StepperOptions)
    step_growth = max(step_growth, stepper_options.min_step_growth)
	step_growth = min(step_growth, stepper_options.max_step_growth)
	step_size   = step_size * step_growth
    step_size   = max(step_size, stepper_options.min_step_size)
	step_size   = min(step_size, stepper_options.max_step_size)
    return step_size
end

# =========================================================================== #
# StepperTemporaries
# =========================================================================== #

abstract type AbstractStepperTemporaries
end

mutable struct StepperTemporaries <: AbstractStepperTemporaries
    t_old::Float64
    t_new::Float64
    dt::Float64

    x_old::Vector{Float64}
    x_new::Vector{Float64}
    x_est::Vector{Float64}
    der_x::Matrix{Float64}
    x_tmp::Matrix{Float64}

    function StepperTemporaries(n_dim::Int, n_intermediates::Int, n_tmp::Int)
        t_old = 0.0
        t_new = 0.0
        dt    = 0.0
        x_old = Vector{Float64}(undef, n_dim)
        x_new = Vector{Float64}(undef, n_dim)
        x_est = Vector{Float64}(undef, n_dim)
        der_x = Matrix{Float64}(undef, n_dim, n_intermediates)
        x_tmp = Matrix{Float64}(undef, n_dim, n_tmp)

        return new(t_old, t_new, dt, x_old, x_new, x_est, der_x, x_tmp)
    end
end

function n_dim(stepper_temp::S) where {S <: AbstractStepperTemporaries}
    return length(stepper_temp.x_old)
end

function n_intermediates(stepper_temp::S) where {S <: AbstractStepperTemporaries}
    return size(stepper_temp.der_x, 2)
end

function n_temporaries(stepper_temp::S) where {S <: AbstractStepperTemporaries}
    return size(stepper_temp.x_tmp, 2)
end

function set_t_old!(stepper_temp::S, t_old::Float64) where {S <: AbstractStepperTemporaries}
    stepper_temp.t_old = t_old
    return stepper_temp
end

function set_t_new!(stepper_temp::S, t_new::Float64) where {S <: AbstractStepperTemporaries}
    stepper_temp.t_new = t_new
    return stepper_temp
end

function set_dt!(stepper_temp::S, dt::Float64) where {S <: AbstractStepperTemporaries}
    stepper_temp.dt = dt
    return stepper_temp
end

function set_x_old!(stepper_temp::S, x_old) where {S <: AbstractStepperTemporaries}
    copy!(stepper_temp.x_old, x_old)
    return stepper_temp
end

function set_x_new!(stepper_temp::S, x_new) where {S <: AbstractStepperTemporaries}
    stepper_temp.x_new = x_new
    return stepper_temp
end

function set_x_est!(stepper_temp::S, x_est) where {S <: AbstractStepperTemporaries}
    stepper_temp.x_est = x_est
    return stepper_temp
end

function swap_x!(stepper_temp::S) where {S <: AbstractStepperTemporaries}
    stepper_temp.x_old = stepper_temp.x_new
    return stepper_temp
end

function swap_fsal!(stepper_temp::S) where {S <: AbstractStepperTemporaries}
    # copy the last column into the first column
    L = n_intermediates(stepper_temp)
    for ii = 1:n_dim(stepper_temp)
        stepper_temp.der_x[ii, 1] = stepper_temp.der_x[ii, L]
    end
    return stepper_temp
end

mutable struct RosenbrockStepperTemporaries <: AbstractStepperTemporaries
    t_old::Float64
    t_new::Float64
    dt::Float64

    x_old::Vector{Float64}
    x_new::Vector{Float64}
    x_est::Vector{Float64}
    x_tmp::Vector{Float64}
    der_x::Matrix{Float64}
    jac::Matrix{Float64}
    dfdt::Vector{Float64}
    tmp::Matrix{Float64}

    function RosenbrockStepperTemporaries(n_dim::Int, n_intermediates::Int, n_tmp::Int)
        t_old = 0.0
        t_new = 0.0
        dt    = 0.0

        x_old = Vector{Float64}(undef, n_dim)
        x_new = Vector{Float64}(undef, n_dim)
        x_est = Vector{Float64}(undef, n_dim)
        x_tmp = Vector{Float64}(undef, n_dim)
        der_x = Matrix{Float64}(undef, n_dim, n_intermediates)
        jac   = Matrix{Float64}(undef, n_dim, n_dim)
        dfdt  = Vector{Float64}(undef, n_dim)
        tmp   = Matrix{Float64}(undef, n_dim, n_tmp)

        return new(t_old, t_new, dt, x_old, x_new, x_est, x_tmp, der_x, jac, dfdt, tmp)
    end
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#!/usr/bin/env julia

doc = """extract_mclade.jl: extracts maximum unrooted clades and the induced trees from mul-trees

Usage:
    extract_mclade.jl <gene_trees> [<out>]

Options:
    -h --help   Show this screen.
"""
include("phylo.jl")
include("external.jl")
include("phylo_extra.jl")

using DocOpt, Logging, ProgressMeter

args = docopt(doc)
inpath = args["<gene_trees>"]
outpath = isnothing(args["<out>"]) ? inpath : args["<out>"]

@info "Now parsing input $inpath"
gene_trees = parse_newicks(open(inpath))
@info "Finished parsing!"

maximal_trees = Tree[]
@showprogress "Computing maximal clades" for (i, gt) in enumerate(gene_trees)
    singly_trees = Tree[]
    for t=singly_cuts(gt)
        if !istrivial(t) && !issinglysubtree(t)
            push!(singly_trees, t)
        end
    end
    mt = sort(collect(maximal_elements(singly_trees));　by=it -> length(it._clades))
    append!(maximal_trees, mt)
end

open(outpath * ".maximal.tre", "w+") do f
    writetrees(f, maximal_trees)
end

open(outpath * ".maximal.clades", "w+") do f
    for t in maximal_trees
        println(f, t._clades)
    end
end

max_stars_repo_name: RuneBlaze/clade-phylo

max_issues_repo_name: RuneBlaze/clade-phylo

max_forks_repo_name: RuneBlaze/clade-phylo


text:
using OctreeBH
using StaticArrays
using PyPlot
#plt.style.use("dark_background")
plt.style.use("default")

function plot_tree_grid(node::Node{N,T}, ix, iy, cc, ax) where {N,T}
	@assert N>=2
	if isLeaf(node)
    	xmin = node.center[ix] - 0.5*node.length[ix]
    	xmax = node.center[ix] + 0.5*node.length[ix]
    	ymin = node.center[iy] - 0.5*node.length[iy]
    	ymax = node.center[iy] + 0.5*node.length[iy]
    	color=cc
    	ax.plot([xmin,xmin],[ymin,ymax], c=color)
    	ax.plot([xmin,xmax],[ymin,ymin], c=color)
    	ax.plot([xmax,xmax],[ymin,ymax], c=color)
    	ax.plot([xmin,xmax],[ymax,ymax], c=color)
	else
		@inbounds for i in 1:2^N
	        #println("open this node")
			plot_tree_grid(node.child[i], ix, iy, cc, ax)
	    end
	end
end

function draw_circle_nonperiodic(x, radius, ix, iy, cc, ax)
    circle = plt.Circle((x[ix]            , x[iy])            , radius, fill=false, color="tab:red")
    ax.add_artist(circle)
end

function draw_circle_periodic(x, radius, ix, iy, boxsizes, cc, ax)
	#draw a circle at point x and the four nearest periodic images
    circle = plt.Circle((x[ix]            , x[iy])            , radius, fill=false, color=cc)
    ax.add_artist(circle)
    circle = plt.Circle((x[ix]+boxsizes[1], x[iy])            , radius, fill=false, color=cc)
    ax.add_artist(circle)
    circle = plt.Circle((x[ix]-boxsizes[1], x[iy])            , radius, fill=false, color=cc)
    ax.add_artist(circle)
    circle = plt.Circle((x[ix]            , x[iy]+boxsizes[2]), radius, fill=false, color=cc)
    ax.add_artist(circle)
    circle = plt.Circle((x[ix]            , x[iy]-boxsizes[2]), radius, fill=false, color=cc)
    ax.add_artist(circle)
end

const BOXSIZE = 1.0

const T = Float64
const N = 2 #spatial dimension
const Npart = 100 #number of particles
const hsml0 = 0.15 #searching radius

using Random
Random.seed!(1114)

#randomly distributing particles
X = [@SVector rand(N) for _ in 1:Npart]

topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree
center = topnode_length .* 0.5

hsml = ones(Npart) .* hsml0
mass = ones(Npart) #particle mass
part = [Data{N,T}(SVector(X[i]), i, hsml[i], mass[i]) for i in eachindex(X)]

#build the tree
tree = buildtree(part, center, topnode_length);

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

clf()
fig, ax = subplots(1, 1, figsize=(6, 6))
ms = 6
ix, iy = 1, 2

ax.plot(getindex.(X,ix), getindex.(X,iy), ".", c="tab:blue", ms=ms)

@time plot_tree_grid(tree, ix, iy, "grey", ax)

X0 = @SVector[0.1, 0.2]
ax.plot(X0[ix], X0[iy], "*", c="tab:green", ms=ms)
idx_ngbs = get_gather_ngb_tree(X0, hsml0, tree, boxsizes)
ax.plot(getindex.(X[idx_ngbs],1), getindex.(X[idx_ngbs],2), ".", c="tab:red", ms=ms)
draw_circle_periodic(X0, hsml0, ix, iy, boxsizes, "tab:green", ax)

X1 = @SVector[0.7, 0.8]
ax.plot(X1[ix], X1[iy], "*", c="tab:green", ms=ms)
idx_ngbs = get_scatter_ngb_tree(X1, tree, boxsizes)
ax.plot(getindex.(X[idx_ngbs],1), getindex.(X[idx_ngbs],2), ".", c="tab:red", ms=ms)
draw_circle_periodic(X1, hsml0, ix, iy, boxsizes, "tab:green", ax)
for i in eachindex(idx_ngbs)
	draw_circle_periodic(X[idx_ngbs[i]], hsml[idx_ngbs[i]], ix, iy, boxsizes, "tab:red", ax)
end

ax.axis([0,boxsizes[ix],0,boxsizes[iy]])
ax.set_aspect(1)